| Date | Description |
| 18-Nov-10 | Original posting |
| 24-Nov-10 | New positive define matrix with improved condition number |
In this assignment you'll implement the Conjugate gradient algorithm. Groups of 2 will have a choice of implementing the algorithm in MPI (on Triton) or in CUDA (on Lincoln). Groups of 3 will implement the algorithm both in MPI and in CUDA. On Triton, you'll conduct scaling studies on with hundreds of cores. Try and make the code go as fast as you can.
The Conjugate Gradient Method (CG) is an iterative method for solving systems of linear equations of the form
Ax= b
where x is unknown (a vector), b is a given vector (the "right hand side") and A is a given matrix with a special form. In particular, A is square, symmetric and positive-definite. A matrix is symmetric if for all i and j, element (i,j) is equal to element (j,i). A matrix is positive definite if all of its eigenvalues are positive; equivalently, xTAx > 0 for all non-zero vectors x, or if the following two properties hold:
If you are interested in digging deeper into the theory about CG, read An Introduction to the Conjugate Gradient Method Without the Agonizing Pain, by J. Shewchuk.
CG is an iterative algorithm, as shown in the following pseudo-code, where b, xk, rk, pk and dk are n-element vectors, A is an n×n element matrix.
x0 =0, r0=b, d0=r0
for k = 1, 2, ....
αk = (rTk-1rk-1)/ (dTk-1Adk-1) // Step length
xk = xk-1 + αk dk-1 // Approximate solution
rk = rk-1 - αk Adk-1 // Residual
βk = (rTkrk)/ (rTk-1rk-1) // Improvement
dk = rk + βk dk-1 // Search Direction
end for
During each iteration, the algorithm performs one matrix-vector multiplication (Adk-1) and two dot products (product (rTk-1rk-1) will have been saved in the previous iteration). Minimize the number of temporaries needed to store quantities from previous iterations. In theory, i.e. with exact arithmetic, the algorithm will converge in n steps. In practice, i.e. with finite arithmetic, roundoff errors will slow the convergence rate down according to the condition number, which is the ratio of the maximum and minimum eigenvalues of the matrix A. (You can dig more deeply into this in the Shewchuck report cited above.)
You will be provided with code to generate positive definite matrices. Run for a fixed number of iterations; if you run to convegence, the running times will grow as you increase n, i.e. when performing a strong scaling study.
If you are working in a group of two, implement cg in either MPI or CUDA. If working in a group of, implement both variants.
Serial C code for CG has been provided on Lincoln, Triton and Lilliput in $(PUB)/HW/A5/CG/C. There is also a matlab version. The C code outputs the residual, which you can use to check correctness:
residual = |Ax - b|2
where |.|2 is the L2 norm computed as the vector length of Ax-b. We have provided a routine to generate positive definite matrices, that correspond to the system of equations obtained by solving the discrete Poisson equation on a 2D mesh (See Solving the Discrete Poisson Equation, by Jim Demmel).
Note that the The generated matrices are of size n2 × n2. When doing week scaling studies with P increasing in increments of 4, you'll have to pick values of n such that the resultant values of n2 are close to integer powers of 2, e.g. for n=32, we have n2=1024. But in order to scale to 2048 × 2048 matrices we will pick n=45, since 452 = 2025 is closest to 2048.
MPI coders: generate the A matrix on a distinguished root process and disribute to the other processes. Do not include the time taken to generate and distribute the A matrix.
The code accepts various command line arguments, (see cmdLine.cpp). You may or may not need all the options, since some make sense for the MPI or CUDA, but not both.
If you would like to add other options that you think may be generally useful to your classmates, post to the A5 moodle forum (there are two forums for this assignment: A5-MPI and A5-CUDA) and we will respond.
The provided code can run for a fixed number of iterations (using the -i option) or until convergence. See the discussion in cg_solve.cpp
CUDA coders: generate the matrix A on the host using the provided code. When taking timings, do not include the time to initialize the A matrix on the host or to transmit to the device, nor the time to return the solution vector to the host.
MPI Coders. Generate the A matrix as described above, and do not include the time to generate and distribute it to the other cores.
The program provides various command line options which are parsed in the source file cmdLine.C. Be sure to look over this file.
When optimizing your code in MPI, try out the Intel Math Kernel Library (MKL) if you have time. Among other things, MKL provides heavily optimized vector (BLAS level 1) and vector-matrix (BLAS level 2) operations used in the CG code (The letter 'd' indicates double precision arithmetic). The Makefile for the provided serial C code include the appropriate setting to load MKL.
The arch file "arch.intel" has been set up to build an executable that uses MKL (see further instructions for compiling jobs on Triton. There is no need to use the MKL Link line Advisor). The interface to the BLAS is found in this section of the Intel Math Kernel Library Reference Manual. Keep in mind that the interface is in Fortran, and by convention, arrays are stored in column major order.
If implementing your code under CUDA, you may try out the CUBLAS, and note any improvment. However, to receive full credit, you must implement a version that does NOT use the CUBLAS.
Pick a problem size such that your runs take a about 5 seconds to complete. A good starting point is somewhere around N2=4Kx4K per core so experiment to determine a reasonable value that suits you.
If you are in a group of 3, compare your MPI and CUDA variants. In particular, determine the number of Triton cores required to achieve performance comparable to your optimized CUDA implementation.
Collect timings in batch mode to avoid processor contention. As described above, do not include the time taken to initialize the A matrix.
Recall that the running time of the application is that of the last core to finish. Don't time individual iterations, but rather, time an entire run, but excluding initialization, calls to MPI_Init() and MPI_Finalize() and any other activities that are not representative of a run with many iterations. communication time.
In the spirit of friendly competition, post your performance results to the appropriate A5 Moodle forum: A5-MPI and A5-CUDA.
For the purposes of grading, we will evaluate your MPI code on Triton and GPU code on Lincoln.
Your Makefile should build 2 targets. (You may name the source modules however you like.)
As in Assignment #3, specify the thread geometry on the make command line, e.g. for 16 × 16 thread blocks
make bx=16 by=16
To obtain different geometries you will need to change the values after the equal sign =. Be sure to run make clean before rebuilding the code with different thread block parameters (and similarly if you are changing the floating point precision).
When you specify bx= and by= make command line arguments, your code will be compiled with correspondingly defined C pre-processor command line arguments
DIM_X = -DBLOCKDIM_X=$(bx) DIM_Y = -DBLOCKDIM_Y=$(by)
where $(bx) and $(by) are the values you specified on the make line. When you specify your thread geometry in your CUDA source you'll refer to these C pre-processor macros as follows
dim3 threads(BLOCKDIM_X,BLOCKDIM_Y);
The Makefile you turn in should build the application with a 16x16 geometry. It's important that you stick to the above convention in order that your application build correctly for other geometries we will test. Assume that all values are integer powers of two.
It's important that you stick to the conventions described in this section so that your code build correctly when we run our auto-grader scripts.
Document your work in a well-written report of about 4-6 pages (DOUBLE SIDED), at the longer end of the range for groups of 3. Discuss your findings carefully. Be sure to explain your results. In particular, analyze, and offer insight into, what you observed. In the case of unexpected behavior or bottlenecks, pose an hypothesis to explain the phenomena and discuss how you would test the hypothesis if you had the time.
Your report should present a clear evaluation of the design of your code, including bottlenecks of the implementation, and describe any special coding or tuned parameters. Provide pseudo code listings that illustrates how you parallelized the code, but do not include full code listings as these will be in your turnin.
Cite any written works or software as appropriate in the text of the report. Citations should include any software you used that was written by someone else, or that was written by you for purposes other than this class.
Your turnin will consist of two parts.
To hand in your assignment electronically, follow instructions to be discussed in the Moodle A5 forum.
Your turnin must include a self-evaluation discussing the division of labor and other aspects of how you worked together as a team (A copy of the form is available HERE).
At the top level of your directory place an ASCII file called README.txt, that lists all the parts of your electronic submission. Or, you may create an html file called index.html, and from there link to a README file as well as the src and report directories.
Copyright © 2010 Scott B. Baden Fri [Thu Nov 18 10:04:34 PST 2010]