The object of this project is to write a parallel subroutine to multiply a sparse matrix by a dense vector (yielding another dense vector), and then to use that "matvec" routine in a routine to solve a symmetric positive definite system of linear equations by the conjugate gradient method, also known as "CG".
CG is an algorithm for solving a linear equation system of the form Ax = b that works when the matrix A is symmetric and has positive eigenvalues. The theory behind CG, while mathematically very pretty, is too far afield from this course for us to go into. (If you want to go into the theory here is a good description.)
All you need to know for this homework is what the algorithm does, which is pretty simple. CG is an iterative method, which means that it starts out with a guess at the value of x (usually you start with the vector of all zeros), and then it does "iterations", one after another, improving the guess at each one, until the guess is close enough. The main work of an iteration involves the multiplication of one (dense) vector by the (sparse) matrix A, and two computations of the dot product of two vectors. The algorithm is outlined on page 7 of these slides.
The matvec routine is partly an exercise in data structures for irregular, sparse objects. You'll use a compressed row data structure that can store any sparse matrix in space proportional to the number of nonzero elements in the matrix.
The data generator we supply will give you the matrix data (only the nonzero entries) distributed over the processors by rows, with n/p rows of an n-by-n matrix on each of p processors. You may rearrange the matrix data if you wish in order to make the sequence of matvecs faster.
There is a sequential Matlab code on the course web site. The Matlab code includes a harness with a data generator and validator, which you can use to check the correctness of your code. The data generator can produce an identity matrix, or the matrix whose graph is a regular k-by-k-by-k grid corresponding to the discretized Poisson equation in three dimensions. (We will also supply you with a parallel data generator that runs on DataStar to generate these matrices in parallel.) The Matlab routine "cgsolve.m" is the conjugate gradient solver. The sparse matvec is built into Matlab; it is just the line "Ad = A*d;"
The data distribution and communication for this problem are up to you. You will probably want to think about not communicating more of the vector than necessary in the matvec. You may want to experiment with rearranging the rows and columns of the matrix to reduce the amount of communication; for fairly small matrices this is unlikely to help but it might be useful for the biggest matrices you can run. You should be able to do pretty large problems on a parallel machine; on my laptop, the sequential Matlab code takes about 2 minutes for a 1,000,000-by-1,000,000 matrix.
We'd like you to turn in the same sort of scaled-speedup experiment and analysis as you did for Homework 2, using the Poisson equation matrix. In this case, you won't expect to see constant time as you scale up the problem, because the number of iterations to converge will grow as the problem size grows. The necessary number of iterations should theoretically grow proportionally to k, the linear dimension of the mesh, which is the cube root of n. Since the time for the matvec is O(n), the overall work grows as O(n^(4/3)).
Grading on this homework will be split between correctness and performance as it was on Homework 2. You should use the largest problems you can (within a run time of a few minutes).
The FAQ on the course home page will be updated as necessary over the next two weeks.