15-859E Programming Assignment 2:
Multigrid

Revision 2: 14 Oct. 98.

• Implement the multigrid algorithm for a simulation problem in two or more dimensions. By multigrid, I mean in particular the V-cycle or FMV-cycle or similar method. You could implement this in C++, C, FORTRAN, Matlab, Java, Mathematica, Maple, or other languages. (I used C++).
• You must implement the relaxation method, restriction, and interpolation yourself down to the grid element level, but you can use other people's code (from the web, from books) for other parts of the program.
• Your program need not be interactive.
• Get your program running on a large grid (at least 60000 unknowns). The larger, the better. In the following discussion, for concreteness I'll assume that your largest (finest) grid is 2-D and 256x256.
• Demonstrate that your program's solutions are accurate. If we define rel_residual(v,f) = norm(f-Av)/norm(f), where norm(x) is the L2 norm (square root of the sum of the squares), then shoot for rel_residual < 1e-6.
• Demonstrate that it is fast. Implement a simpler method for reference (this could be your relaxation method on the finest grid, for example) and compare the speed of your two methods at a range of sizes, say 8x8, 16x16, 32x32, ... 256x256. You probably won't have the patience to run your simple method on the finer grids (it would take hours), but try to make your speed comparison table as complete as possible by the deadline.
• When preparing this table, make sure that each of the runs is accurate within the above tolerance. It's not very meaningful to compare the speeds of methods that produce solutions of wildly different accuracies.
• Make a graph of number of unknowns (probably N^2) versus time for these two algorithms, labeled clearly with units, showing (hopefully) that your simple method is very slow, but multigrid is fast. You could use gnuplot for this.
• Fit a function to the empirical time cost curves. Are they approximately c*N^k for some k, 1<=k<=6? Can you determine the constants in the empirical cost formula?
• Make a picture of the solution function. You might try Maple, Mathematica, Matlab, or OpenGL for this. See the class software page for info on availability of the first three packages. The picture could be grayscale or color, 2-D (contour plot, say), or 3-D.
• Finish your code by the date listed at top. To submit your program, copy your source code, executable, and test files (if any) to students/yourlastname/multigrid. Include a README file explaining: what you did, whether you wrote all the code, and if not, what other code you used (give author, URL, program name, book and page... as appropriate), what system your executable runs on (e.g. 900 MHz Pentium 4). Also include the comparative table mentioned earlier, and discussion of the functions that fit it. If you find a text README file too restrictive, other formats are OK (Latex, Word, Framemaker, Mathematica...).
• During class on 22 Oct., turn in a printout of your README, a printout of your comparative time graph, and a picture of your solution. We'll discuss the outcome in class, but we probably won't do demos, since most of these programs will be non-interactive.

Tips

• If you choose to simulate Poisson's equation, see Briggs' equation 1.4. For variety, you might want to try various right hand sides (f).
• If you're more ambitious, you might try multigrid on another type of problem, say from computer vision, elasticity, heat diffusion, or fluid mechanics. In general, multigrid usually works well on elliptical partial differential equations with symmetric positive definite matrices. But don't get too carried away - I urge you to try multigrid on a square, power-of-two grid hierarchy before attempting new applications, irregular domains, or unstructured meshes. If you want to do new research using multigrid, you might want to save that and make it your project at the end of the semester.
• One way to get ambitious without as much risk is to do the basic assignment first, get it running fast enough (using optimization tips in Makefile plus code tuning and run it on a smallsh grid) that you can get a full solution in less than a second, then implement interactive control of the right hand side (f).
• I found that restriction with full weighting works converges far faster than injection, but you may want to do some experimentation. For a relaxation method, I recommend weighted Gauss-Seidel, but you may want to experiment.
• Note that it's never necessary to allocate space for the A matrix! In my implementation I had 2-D arrays for v, f, and r at power-of-two sizes, and those were the only large memory allocations. (Briggs doesn't talk about keeping an r grid, but I found that this simplified the restriction code).
• I'll assume in most of the following tips that you're programming in C++.
• Use double precision (64 bit) floating point, not single precision.
• The recommended starter code is in pubdir/src/multigrid, where pubdir = /afs/cs/project/classes-ph/859E/pub . Take a look at it here. simp.cc is a very short program showing how to create a grid (2-D array) of arbitrary size with SVL's "Mat" class, and how to time a program in software.
• If you choose Poisson's equation, it helps to make your grids (2^k+1)x(2^k+1) on a side, with the perimeter values of the v grid permanently set to zero, and only the center (2^k-1)x(2^k-1) sub-grid unknown. Again, see Briggs equation 1.4. The principal advantage of this is that the 3x3 neighborhood of any of the unknowns is within bounds, so for relaxation and restriction you won't need special case code to handle neighbors of an edge or corner grid node. The setup for other simulation problems may differ a bit.
• The small program levi.cc shows how to use a nifty argument parser I wrote. I find this very handy for adding command line options to a program. You could use it, for example, to set the grid size, number of iterations here and there (Briggs' parameters nu0, nu1, nu2), select which relaxation algorithm to use, set omega parameter for weighted Gauss-Seidel (if you use that), or turn on debugging. Source and documentation for the argument parser are here.
• I recommend you debug by working with small grids (9x9 or less, say) and printing them after each operation. To facilitate this, use the formatting features of printf (or ostream) to get things to line up, e.g.

       0.00   1.00   2.00              0 1 2
      90.33 140.67 100.05      not     90.33 140.67 100.05
       1.00   2.00   5.00              1 2 5
      

  The former was generated with printf(" %5.2f", m[i][j]), the latter with printf(" %g", m[i][j]).
• On an SGI, if you use SVL, debug with compiler options "-g -DVL_CHECKING" (the latter causes subscriptrange checking in SVL) and linker option "-lsvl.dbg" but compile for runs on big grids with "-O2" and without "-DVL_CHECKING" and link with "-lsvl". For me, this change made my program 5x faster. On Suns, the options are the same except for "-O" instead of "-O2". These issues are mentioned in Makefile.sgi and Makefile.sun.
• I found it helpful to have an option that runs the multigrid algorithm silently except for printing the relative residual (see above) of the final solution at the end. That way I could run it with different parameters (varying the number of relaxation steps, or switching relaxation methods, say) to find the combination that most quickly found an accurate solution. You could also make the whole thing interactive, if you wish.

15-859, Hierarchical Methods for Simulation
Paul Heckbert