03/23/07 CS267 Lecture 20 1CS 267: MultigridKathy Yelickwww.cs.berkeley.edu/~yelick/cs267_s0703/23/07 CS267 Lecture 16 2Algorithms for 2D (3D) Poisson Equation (N = n2(n3) vars)Algorithm Serial PRAM Memory #Procs° Dense LU N3NN2N2° Band LU N2 (N7/3)N N3/2 (N5/3)N (N4/3)° Jacobi N2(N5/3) N (N2/3) N N° Explicit Inv. N2log N N2N2° Conj.Gradients N3/2 (N4/3) N1/2(1/3)*log N N N° Red/Black SOR N3/2 (N4/3) N1/2 (N1/3) N N° Sparse LU N3/2 (N2)N1/2 N*log N (N4/3) N° FFT N*log N log N N N° Multigrid N log2NN N° Lower bound N log N NPRAM is an idealized parallel model with zero cost communicationReference: James Demmel, Applied Numerical Linear Algebra, SIAM, 1997.03/23/07 CS267 Lecture 16 3Algorithms for 2D (3D) Poisson Equation (N = n2(n3) vars)Algorithm Serial PRAM Memory #Procs° Dense LU N3NN2N2° Band LU N2 (N7/3)N N3/2 (N5/3)N (N4/3)° Jacobi N2(N5/3) N (N2/3) N N° Explicit Inv. N2log N N2N2° Conj.Gradients N3/2 (N4/3) N1/2(1/3)*log N N N° Red/Black SOR N3/2 (N4/3) N1/2 (N1/3) N N° Sparse LU N3/2 (N2)N1/2 N*log N (N4/3) N° FFT N*log N log N N N° Multigrid N log2NN N° Lower bound N log N NPRAM is an idealized parallel model with zero cost communicationReference: James Demmel, Applied Numerical Linear Algebra, SIAM, 1997.03/23/07 CS267 Lecture 16 4Review of Previous Lectures° Review Poisson equation° Overview of Methods for Poisson Equation° Jacobi’s method° Red-Black SOR method° Conjugate Gradients° FFT° Multigrid° Comparison of methodsReduce to sparse-matrix-vector multiplyNeed them to understand Multigrid03/23/07 CS267 Lecture 16 5Poisson’s equation in 1D: ∂2u/∂x2= f(x)2 -1 -1 2 -1-1 2 -1-1 2 -1-1 2T =2-1 -1Graph and “stencil”03/23/07 CS267 Lecture 16 62D Poisson’s equation° Similar to the 1D case, but the matrix T is now° 3D is analogous4 -1 -1-1 4 -1 -1-1 4 -1-1 4 -1 -1-1 -1 4 -1 -1 -1 -1 4 -1-1 4 -1-1 -1 4 -1-1 -1 4T =4-1-1-1-1Graph and “stencil”03/23/07 CS267 Lecture 16 7Multigrid Motivation° Recall that Jacobi, SOR, CG, or any other sparse-matrix-vector-multiply-based algorithm can only move information one grid cell at a time• Take at least n steps to move information across n x n grid° Can show that decreasing error by fixed factor c<1 takes Ω(log n) steps• Convergence to fixed error < 1 takes Ω(log n ) steps° Therefore, converging in O(1) steps requires moving information across grid faster than to one neighboring grid cell per step03/23/07 CS267 Lecture 16 8Multigrid Methods• We studied several iterative methods• Jacobi, SOR, Guass-Seidel, Red-Black variations, Conjugate Gradients (CG)• All use sparse matrix-vector multiply (nearest neighbor communication on grid)• Key problem with iterative methods is that:• detail (short wavelength) is correct • convergence controlled by coarse (long wavelength) structure • In simple methods one needs of order N2iterations to get good results• Ironically, one goes to large N (fine mesh) to get detail • If all you wanted was coarse structure, a smaller mesh would be fine• Basic idea in multigrid is key in many areas of science• Solve a problem at multiple scales• We get coarse structure from small N and fine detail from large N• Good qualitative idea but how do we implement? Slide source: Geoffrey Fox and(indirectly) Ulrich Ruede03/23/07 CS267 Lecture 16 9Gauss Seidel is Slow I° Take Laplace’s equation in the Unit Square with initial guess as φ = 0 and boundary conditions that are zero except on one side° For N=31 x 31 Grid it takes around 1000 (N2) iterations to get a reasonable answerBoundary Conditions Exact SolutionSlide source: Geoffrey Fox and(indirectly) Ulrich Ruede03/23/07 CS267 Lecture 16 10Gauss Seidel is Slow II1 Iteration 10 Iterations100 Iterations1000 IterationsSlide source: Geoffrey Fox and (indirectly) Ulrich Ruede03/23/07 CS267 Lecture 16 11Multigrid Overview° Basic Algorithm:• Replace problem on fine grid by an approximation on a coarser grid• Solve the coarse grid problem approximately, and use the solution as a starting guess for the fine-grid problem, which is then iteratively updated• Solve the coarse grid problem recursively, i.e. by using a still coarser grid approximation, etc.° Success depends on coarse grid solution being a good approximation to the fine grid03/23/07 CS267 Lecture 16 12Same Big Idea used elsewhere° Replace fine problem by coarse approximation, recursively° Multilevel Graph Partitioning (METIS): • Replace graph to be partitioned by a coarser graph, obtained via Maximal Independent Set• Given partitioning of coarse grid, refine using Kernighan-Lin° Barnes-Hut (and Fast Multipole Method) for computing gravitational forces on n particles in O(n log n) time:• Approximate particles in box by total mass, center of gravity• Good enough for distant particles; for close ones, divide box recursively° All examples depend on coarse approximation being accurate enough (at least if we are far enough away)03/23/07 CS267 Lecture 16 13Multigrid uses Divide-and-Conquer in 2 Ways° First way:• Solve problem on a given grid by calling Multigrid on a coarse approximation to get a good guess to refine° Second way:• Think of error as a sum of sine curves of different frequencies• Same idea as FFT solution, but not explicit in algorithm• Each call to Multigrid responsible for suppressing coefficients of sine curves of the lower half of the frequencies in the error (pictures later)03/23/07 CS267 Lecture 16 14Multigrid Sketch in 1D• Consider a 2m+1 grid in 1D for simplicity• Let P(i)be the problem of solving the discrete Poisson equation on a 2i+1 grid in 1D• Write linear system as T(i) * x(i) = b(i)• P(m), P(m-1), … , P(1)is sequence of problems from finest to coarsest03/23/07 CS267 Lecture 16 15Multigrid Sketch in 2D• Consider a 2m+1 by 2m+1 grid in 2D• Let P(i)be the problem of solving the discrete Poisson equation on a 2i+1 by 2i+1 grid in 2D• Write linear system as T(i) * x(i) = b(i)• P(m), P(m-1), … , P(1)is sequence of problems from finest to coarsest03/23/07 CS267 Lecture 16 16Multigrid HierarchyRelaxInterpolateRestrictRelaxRelaxRelaxRelaxSlide source: Geoffrey Fox03/23/07 CS267 Lecture 16 17Basic Multigrid Ideas• In picture, relax
View Full Document
Unlocking...