Conjugate Gradient Method

First write a serial vesion of CG that works correctly. Because we know a priori that a parallel version will be needed, here's how I would break down the code into functions/subroutines. This uses Matlab syntax, which allows multiple return arguments from function. If you use C, you'll have to modify things to use addresses.

•  [A, b, startrow, nrows, ncols, nonzeros, x0] = 
              getAb(fname_A, fname_b, fname_xinit, fname_solution)
      

  reads the files with the given names. startrow is the starting row number for the block row that a process is assigned. In general that would be rank*n/p + 1. For a serial version, startrow = 1. Notice that the row numbers are indexed starting with 1, not 0. If you have trouble with returning memory allocated in a subprogram, add a task input to getAB:

  1. task = 1: only read the first line: nrows, ncols, nnz, startrow

  2. [allocate the arrays in the main program here]

  3. task = 2: read the full data

  but you should not need to do that.

•  [maxits, tol] = setCGparameters() 

•  [d, r0, rho] = InitializeCG(A, b, x0, nrow, ncols, nnz, startrow)

  This carries out iteration 0, before the while loop that does most of the iteration.

•  [x, pltret] = cg(A, b, x0, maxits, tol, nrow, ncols, nnz, startrow) 

  implements the main loop.

•  OK = checksolution(A, b, x) 

•  OK = writesolutionfile(x) 

The Matlab code is not broken down that way (it was never intended to be parallel). The files are

• drivecg.m: main driver program that calls the other functions and then plots the resulting residual norm history.
• cg.m: the combined initialization and the iteration phases
• setA.m: creates A and b. We'll need to change that to read in A and b from a file via getAb()
• writearray.m:
• writearray.m: writes out an ASCII file containing the data for a given array.
• testwrite.m: short test function for writearray, which you don't need.

• cg.pdf: the notes from Monday, 06 Apr 2009. You don't really need these for writing a CG code, just to explain any seemingly bizarre results. It does include a slide with the statement of the algorithm, if you can't read Matlab scripts.
• cleanup: shell script that deletes the files created by drivecg.m
• As a start for the larger problem sizes, there is a version that only creates the systems and writes the files without invoking cg.m. It uses writearray.m and setA.m from above, and replaces drivecg.m with SetupSystem.m. SetupSystem.m calls Matlab's spy() plot, so comment that part out if you don't want to see the sparsity pattern of the coefficient matrix. All of this is bundled in a
• gzipped tarfile

New Matlab file additions

1. 08 Apr 2009: have put an additional tarball with setup methods This creates a directory parsetup, with one new file: SetupParallel.m. This does not run the CG algorithm, all it does is create symmetric positive definite systems and write them out to files with the naming convention coeffmatrix001, coeffmatrix002, ..., coeffmatrix***. Those files you can then ship over to odin via scp, and put one in /tmp on each node you have an allocation for. That is a temporary file location and gets cleaned up periodically, but is guaranteed to be on each node's local hard drive, not on a network mounted filesystem. Then your MPI code can each read in the data it needs from /tmp/coeffmatrix*** in parallel. The rhs vector is also partitioned, but the initial x is not - still, every process needs a copy of the full x vector.

2. 08 Apr 2009: Update to writearray.m: SetupParallel.m calls a version of writearray.m that is not the same as the original one listed above. This one takes a third (optional) argument to allow specifying where the row numbering will start when arrays are being written out. You can use this in place of the older one, since this rowoffset argument is an optional one. But the older one will not work if called from SetupParallel.m.

3. Added, 09 Apr 2009: A version of the Matlab linear system generator that uses sparse matrices. This is one that you can use to generate large systems - it took me less than 40 seconds to create the files for a 60001 x 60001 matrix partitioned into 64 block rows. Storage for that problem size is about 7.4 Mbytes. The systems are stored in COO format: a header line with nrow, ncols, nnz, and startrow number, followed by all the nonzero entries, one per line. Each line has a row number, a column number, and the value. All of this is 1-based indexing.

   Beware one thing: some of the same names have been set up for the dense and sparse versions, because the Matlab code is nearly identical in either case (and cg.m required zero changes, since in Matlab w = A*d works correctly for dense or sparse matrices A). Notice that the Matlab codes are only for creating test cases for a parallel code to work on - Matlab is not parallel itself. The utilities in some of the files below may help you in debugging, but you do not have to implement anything but solver, and a driver that reads in the partitioned systems and writes out the solution and residual history.

   1. The entire tarball is in parsetup2.tar which creates a subdirectory called parsetup2.
   2. SetupParallel.m has only minor differences from the previous one, but does have some efficiency considerations in addition to using sparse matrices.
   3. setA.m similarly barely differs from the previous version. Do a diff (or better yet, an xxdiff) on the files to see what changed.
   4. writearraysparse.m required only four line changes from the original writearray.m. To keep both available, the new parsetup2 directory has the dense matrix version in writearraydense.m, with a link to the actual one to use.

4. Added, 09 Apr 2009: to go along with the sparse matrix generator and partitioner, a set of Matlab codes that read in, assemble the matrix A from the partitioned files, and then runs CG is tar'd up in assemble.tar. It contains
   1. runcg.m which calls
   2. readsystem.m that in turn uses
   3. getarray.m and
   4. getparameters.m

   You do not need to use the codes in the "assemble" directory - they are just there to help with any debugging you need to do.

More Example Files

• Added, 21 Apr 2009: if it helps you with BLAS linking here is a serial CG gzipped tar file with a small 4900 x4900 test problem. This is mainly for those wanting a working makefile, and does not follow the specifications for the CG version in this project. If you download or use this, be sure to read the file "info" in the created directory. Otherwise you can get seriously fouled up and confused. And yes, I intentionally put it together from Fortran codes so you'd have yet another example for reading that language.