Problem Set #5 (Part II)¶

Block Partitioning + Threads¶

Recall the approach to parallelization outlined in the Design Patterns for Parallel Programming. The first step is to find concurrency – how can we decompose the problem into independent (ish) pieces and then (second step) build an algorithm around that?

The first approach we want to look at with norm is to simply separate the problem into blocks. Computing the sum of squares of the first half (say) of a vector can be done independently of the second half. We can combine the results to get our final result.

The program pnorm.cpp is a driver program that calls the function partitioned_two_norm_a in the file pnorm.hpp. The program as provided uses C++ threads to launch workers, each of which executes a specified block of the vector – from begin to end. Unfortunately, the program is not correct – more specifically, there is a race that causes it to give the wrong answer.

Fix the implementation of partitioned_two_norm_a so that it no longer has a race and so that it achieves some parallel speedup. (Eliminating the race by replacing the parallel solution with a fully sequential solution is not the answer I am looking for.)

Answer the following questions in Questions.rst.

• What was the data race?

• What did you do to fix the data race? Explain why the race is actually eliminated (rather than, say, just made less likely).

• How much parallel speedup do you see for 1, 2, 4, and 8 threads?

Block Partitioning + Tasks¶

One of the C++ core guidelines that we presented in lecture related to concurrent programming was to prefer using tasks to using threads – to think in terms of the work to be done rather than orchestrating the worker to do the work.

The program fnorm.cpp is the driver program for this problem. It calls the functions partitioned_two_norm_a and partitioned_two_norm_b, of which skeletons are provided.

Complete the function partitioned_two_norm_a, using std::async to create futures. You should save the futures in a vector and then gather the returned results (using get) to get your final answer. Use std::launch::async as your launch policy. The skeletons in fnorm.hpp suggest using helper functions as the tasks, but you can use lambda or function objects if you like.

Once you are satisfied with partitioned_two_norm_a, copy it verbatim to partitioned_two_norm_b and change the launch policy from std::launch::async to std::launch::deferred.

Answer the following questions in Questions.rst.

• How much parallel speedup do you see for 1, 2, 4, and 8 threads for partitioned_two_norm_a?

• How much parallel speedup do you see for 1, 2, 4, and 8 threads for partitioned_two_norm_b?

• Explain the differences you see between partitioned_two_norm_a and partitioned_two_norm_b.

Cyclic Partitioning + Your Choice¶

Another option for decomposing the two norm computation is to use a “cyclic” or “strided” decomposition. That is, rather than each task operating on contiguous chunks of the problem, each task operates on interleaved sections of the vector.

Complete the body of cyclic_two_norm in cnorm.hpp. The implementation mechanism you use is up to you (threads or tasks). You can use separate helper functions, lambdas, or function objects. The requirements are that the function parallelizes according to the specified number of threads, that it uses a cyclic partitioning (strides), and that it computes the correct answer (within the bounds of floating point accuracy).

Answer the following questions in Questions.rst.

• How much parallel speedup do you see for 1, 2, 4, and 8 threads?

• How does the performance of cyclic partitioning compare to blocked? Explain any significant differences, referring to, say, performance models or CPU architectural models.

Divide and Conquer + Tasks (AMATH 583 ONLY)¶

Another approach to decomposing a problem – particularly if you can’t (or don’t want to) try to divide it into equally sized chunks all at once – is to “divide and conquer.” In this approach, our function breaks the problem to be solved into two pieces (“divide”), solves the two pieces, and then combines the results. The twist is that to solve each of the two pieces (“conquer”), the function calls itself.

To see what we mean, suppose we want to take the sum of squares of the elements of a vector (as a precursor for computing the Euclidean norm, say) using a function sum_of_squares. Well, we can abstract that away by pretending we can compute the sum of squares of the first half of the vector and the sum of squares of the second half of the vector. And, we abstract that away by saying we can compute those halves by calling sum_of_sqares.

double sum_of_squares(const Vector& x, size_t begin, size_t end) {
  return sum_of_squares(x, begin, begin+(end-begin)/2)
       + sum_of_squares(x, begin+(end-begin)/2, end) ;
}

To use an over-used phrase – see what I did there? We’ve defined the function as the result of that same function on two smaller pieces. This is completely correct mathematically, by the way, the sum of squares of a vector is the sum of the sum of squares of two halves of the vector.

Almost. It isn’t correct for the case where end - begin is equal to unity – we can no longer partition the vector. So, in this code example, as with the mathematical example in the sidebar, we need some non-self-referential value (a base case).

With this computation there are several options for a base case. We could stop when the difference between begin and end is unity. We could stop when the difference is sufficiently small. Or, we could stop when we have achieved a certain amount of concurrency (have descended enough levels).

In the latter case:

double sum_of_squares(const Vector& x, size_t begin, size_t end, size_t level) {
  if (level == 0) {
    double sum = 0;
    // code that computes the sum of squares of x from begin to end into sum
    return sum;
  } else {
    return sum_of_squares(x, begin, begin+(end-begin)/2, level-1)
         + sum_of_squares(x, begin+(end-begin)/2, end , level-1) ;
  }
}

Now, we can parallelize this divide and conquer approach by launching asynchronous tasks for the recursive calls and getting the values from the futures.

The driver program for this problem is in rnorm.cpp.

Complete the implementation for recursive_two_norm_a in rnorm.hpp. Your implementation should use the approach outlined above.

Answer these questions in Questions.rst.

• How much parallel speedup do you see for 1, 2, 4, and 8 threads?

• What will happen if you use std:::launch::deferred instead of std:::launch::async when launching tasks? When will the computations happen? Will you see any speedup? For your convenience, the driver program will also call recursive_two_norm_b – which you can implement as a copy of recursive_two_norm_a but with the launch policy changed.

General Questions¶

Answer these questions in Questions.rst.

• For the different approaches to parallelization, were there any major differences in how much parallel speedup that you saw?

• You may have seen the speedup slowing down as the problem sizes got larger – if you didn’t keep trying larger problem sizes. What is limiting parallel speedup for two_norm (regardless of approach)? What would determine the problem sizes where you should see ideal speedup? (Hint: Roofline model.)

Partitioning Matrix Vector Product¶

But, before diving right into that, let’s apply our parallelization design patterns. One way to split up the problem would be to divide it up like we did with two_norm:

(For illustration purposes we represent the matrix as a dense matrix.) Conceptually this seems nice, but mathematically it is incorrect. Since we are partitioning the matrix so that each task has all of the columns of A we also need to be able to access all of the columns of x – so we can’t really partition x though we can partition y and we can partition A – row wise.

If we did want to split up x, we could partition A column wise but then each task has all of the rows of A so we can’t partition y.

(Answer not required in Questions.rst but you should think about this.) What does this situation look like if we are computing the transpose product? If we partition A by rows and want to take the transpose product against x can we partition x or y? Similarly, if we partition A by columns can we partition x or y?

Partitioning Sparse Matrix Vector Product¶

Now let’s return to thinking about partitioning sparse matrix-vector product. We have three sparse matrix representations: COO, CSR, and CSC. We can really only partition CSR and CSC and we can really only partition CSR by rows and CSC by columns. (The looping over the elements in these formats is over rows and columns for CSR and CSC, respectively in the outer loop – which can be partitioned. The loop orders can’t be changed, and the inner loop can’t really be partitioned based on index range.)

So that leaves four algorithms with partitionings – CSR partitioned by row for matvec and transpose matvec and CSC partitioned by column for matvec and transposed matvec. For each of the matrix formats, one of the operations allowed us to partition x and the other to partition y.

(Answer not required in Questions.rst but you should think about this and discuss on Piazza.) Are there any potential problems with not being able to partition x? Are there any potential problems with not being able to partition y? (In parallel computing, “potential problem” usually means race condition.)

Implement the safe-to-implement parallel algorithms as the appropriate overloaded member functions in CSCMatrix.hpp and CSRMatrix.hpp. The prototype should be the same as the existing member functions, but with an additional argument specifying the number of partitions. You will need to add the appropriate #include statements for any additional C++ library functionality you use. You will also need to modify the appropriate dispatch function in amath583sparse.cpp (these are already overloaded with number of partitions but only call the sequential member functions) – search for “Fix me”.

Answer the following in Questions.rst.

• Which methods did you implement?

• How much parallel speedup do you see for the methods that you implemented for 1, 2, 4, and 8 threads?

Graph Data Files¶

I have uploaded a set of data files that are somewhat larger – enough to be able to see meaningful speedup with parallelization but not so large as to overwhelm your disk space or network connection (hopefully). The Canvas announcement about this assignment has a link to the Google Drive folder containing the files.

Several of the graphs are available in ASCII format (“.mtx”) or binary format (“.mmio”). I suggest using the binary format as that will load much more quickly when you are doing experiments.

If you want to find out more specific information about any of these graphs, they were downloaded from https://sparse.tamu.edu where you can search on the file name to get all of the details.

Profiling¶

Run pagerank.exe as it is (unaccelerated) on some of the large-ish graphs. You should be seeing compute times on the order of a few seconds to maybe ten or so seconds.

I checked the pagerank program with a profiler and found that it spends a majority of its time in mult, more specifically in CSR::t_matvec. This is the function we want to speed up.

(And by speeding up – I mean, maybe, use – since we already sped it up above. If you add the number of partitions as an argument to mult, the program should “work”.)