4/12/2004 CS267 Lectures 19 1CS 267 Applications of Parallel ComputersSolving Linear Systems Arising from PDEsUsing MultigridKatherine Yelickwww.cs.berkeley.edu/~yelick/cs2674/12/2004 CS267 Lectures 19 2Where Was I Last Week?! Whirlwind tour of Japanese Supercomputing! Computer facilities! The Earth Simulator (40 Tflops)! A Supercluster of PCs (11 Tflops)! Three government agencies! MEXT, CSTP, METI! Universities and Research Labs! University of Tokyo! University to Tsukuba! RIST, AIST! Vendors! NEC, Hitachi, Fujitsu, Sony4/12/2004 CS267 Lectures 19 3Review of Earlier Lecture on Solving PDEs• Review Poisson equation• Overview of Methods for Poisson Equation• Jacobi’s method• Red-Black SOR method• Conjugate Gradients• FFTToday:• MultigridReduce to sparse-matrix-vector multiplyNeed them to understand Multigrid4/12/2004 CS267 Lectures 19 42D Poisson’s equation: ∂∂∂∂2u/∂∂∂∂x2+ ∂∂∂∂2u/∂∂∂∂y2 = f(x,y)! The the matrix T for a regular mesh is! 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”4/12/2004 CS267 Lectures 19 5Algorithms for 2D Poisson with N UnknownsAlgorithm Serial PRAM Memory #Procs! Dense LU N3NN2N2! Band LU N2NN3/2N! Jacobi N2NNN! Explicit Inv. N2log N N2N2! Conj.Grad. N 3/2N 1/2*log N N N! RB SOR N 3/2N 1/2NN! Sparse LU N 3/2N 1/2N*log N 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 communication4/12/2004 CS267 Lectures 19 6Multigrid Motivation! The FFT is O(log N), but has an O(N) transpose that prevents scaling! Recall that Jacobi, SOR, CG, or any other sparse-matrix-vector-multiply-based algorithm can only move information one grid call at a time! Takes N1/2steps to get information across grid! Can show that decreasing error by fixed factor c<1 takes ΩΩΩΩ(log n) steps ! Convergence to fixed error < 1 takes ΩΩΩΩ(log n ) steps! SOR and CG take N1/2steps! Therefore, converging in O(1) steps requires moving information across grid faster than to one neighboring grid cell per step4/12/2004 CS267 Lectures 19 7Multigrid Motivation4/12/2004 CS267 Lectures 19 8Multigrid 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 be 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 grid4/12/2004 CS267 Lectures 19 9Same Big Idea Used in Other Problems! 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):! Approximate leaf node of QuadTree containing many particles by Center-of-Mass and Total-Mass (or multipole expansion)! Approximate internal nodes of QuadTree by combining expansions of children! Force in any particular node approximated by combining expansions of a small set of other nodes! All examples depend on coarse approximation being accurate enough (at least if we are far enough away)4/12/2004 CS267 Lectures 19 10Multigrid 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 errordifferent grid levels will damp different frequency errors4/12/2004 CS267 Lectures 19 11Multigrid Sketch on a Regular 1D Mesh! 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 a sequence of problems from finest to coarsest4/12/2004 CS267 Lectures 19 12Multigrid Sketch on a Regular 2D Mesh! Consider a 2m+1 by 2m+1 grid! 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 a sequence of problems from finest to coarsest4/12/2004 CS267 Lectures 19 13Multigrid Operators! For problem P(i):! b(i) is the RHS and ! x(i) is the current estimated solution ! (T(i) is implicit in the operators below.)! Multigrid will be defined by a repeated sequence of operators! The multigrid operators average values on neighboring grid points! Neighboring grid points on coarse problems are far away in fine problems, so information moves quickly on coarse problems! Levels will be constructed explicitly! The next slide gives an overview of the operators; details will come laterboth live on grids of size 2i-14/12/2004 CS267 Lectures 19 14Multigrid Operators! For problem P(i):! b(i) is the RHS and ! x(i) is the current estimated solution ! The restriction operator R(i) maps P(i)to P(i-1)! Restricts problem on fine grid P(i)to coarse grid P(i-1)! Uses sampling or averaging! Right-hand sided is also restricted: b(i-1)= R(i) (b(i))! The interpolation operator In(i-1) maps solution x(i-1) to x(i)! Maps one approximate solution to another! Interpolates solution on coarse grid P(i-1)to fine grid P(i)! x(i) = In(i-1)(x(i-1))! The solution operator S(i) takes P(i)and improves solution x(i) ! Uses “weighted” Jacobi or SOR on a single level of the grid! x improved(i) = S(i) (b(i), x(i))both live on grids of size 2i-14/12/2004 CS267 Lectures 19 15Multigrid V-Cycle AlgorithmFunction MGV ( b(i), x(i) )… Solve T(i)*x(i) = b(i) given b(i) and an initial guess for x(i)… return an improved x(i)if (i = 1) compute exact solution x(1) of P(1)only 1 unknownreturn x(1)else x(i) = S(i) (b(i), x(i)) improve solution by damping high frequency errorr(i) = T(i)*x(i) - b(i) compute residuald(i) = In(i-1) ( MGV( R(i) ( r(i) ), 0 ) ) solve T(i)*d(i) = r(i) recursively x(i) = x(i) - d(i) correct fine grid solutionx(i) = S(i) ( b(i), x(i) ) improve
View Full Document
Unlocking...