13/16/2004 CS267 Lecture 15 1CS 267Applications of Parallel ComputersSolving Linear Systems Arising from PDEsKathy Yelickwww.cs.berkeley.edu/~yelick/cs2673/16/2004 CS267 Lecture 15 2Outline• Review Poisson equation• Overview of Methods for Poisson Equation• Jacobi’s method• Red-Black SOR method• Conjugate Gradient• FFT• Multigrid (later lecture this semester)• Reduce to sparse-matrix-vector multiply• Need them to understand Multigrid3/16/2004 CS267 Lecture 15 3Recap of “Sources of Parallelism” Lecture• Discrete event systems:• Examples: “Game of Life,” logic level circuit simulation. • Particle systems:• Examples: billiard balls, semiconductor device simulation, galaxies.• Lumped variables depending on continuous parameters:• ODEs, e.g., circuit simulation (Spice), structural mechanics, chemical kinetics.• Continuous variables depending on continuous parameters:• PDEs, e.g., heat, elasticity, electrostatics.• A given phenomenon can be modeled at multiple levels.• Many simulations combine more than one of these techniques.3/16/2004 CS267 Lecture 15 4Continuation of RecapExamples of PDE systems include• Hyperbolic problems (waves):• Quantum mechanics: Wave-function(position,time)• Elliptic (steady state) problems:• Electrostatic or Gravitational Potential: Potential(position)• Parabolic (time-dependent) problems:• Heat flow: Temperature(position, time)• Diffusion: Concentration(position, time)Many problems combine features of above• Fluid flow: Velocity,Pressure,Density(position,time)• Elasticity: Stress,Strain(position,time)3/16/2004 CS267 Lecture 15 5Terminology• Term hyperbolic, parabolic, elliptic, come from special cases of the general form of a second order linear PDEa * ∂∂∂∂2u/∂∂∂∂x2+ b * ∂∂∂∂2u/∂∂∂∂x∂∂∂∂y + c * ∂∂∂∂2u/∂∂∂∂y2 + d * ∂∂∂∂u/∂∂∂∂x + e * ∂∂∂∂u/∂∂∂∂y + f = 0(y denotes time for hyperbolic and parabolic equations)• Analog to solutions of general quadratic equationa * x2+ b * x*y + c * y2 + d * x + e * y + f = 0Ellipse: 4*a*c – b^2 > 0Hyperbola: 4*a*c – b^2 < 0Parabola: 4*a*c – b^2 = 03/16/2004 CS267 Lecture 15 6Solving PDEsSolution strategies:• Hyperbolic problems (waves):• Use explicit time-stepping • Solution at each point depends on neighbors at previous time• Elliptic (steady state) problems:• Everything depends on everything else• This means locality is harder to find than in hyperbolic problems• Parabolic (time-dependent) problems:• Involves an elliptic solve at each time-step• Focus on elliptic problems• Canonical example is the Poisson equation∂∂∂∂2u/∂∂∂∂x2+ ∂∂∂∂2u/∂∂∂∂y2 + ∂∂∂∂2u/∂∂∂∂z2 = f(x,y,z)23/16/2004 CS267 Lecture 15 7Poisson’s equation arises in many models• Heat flow: Temperature(position, time)• Diffusion: Concentration(position, time)• Electrostatic or Gravitational Potential: Potential(position)• Fluid flow: Velocity,Pressure,Density(position,time)• Quantum mechanics: Wave-function(position,time)• Elasticity: Stress,Strain(position,time)• Variations of pure Poisson involve other coefficients3D: ∂∂∂∂2u/∂∂∂∂x2+ ∂∂∂∂2u/∂∂∂∂y2 + ∂∂∂∂2u/∂∂∂∂z2 = f(x,y,z)2D: ∂∂∂∂2u/∂∂∂∂x2+ ∂∂∂∂2u/∂∂∂∂y2 = f(x,y)1D: ∂∂∂∂2u/∂∂∂∂x2= f(x)f represents the sources and boundary conditions3/16/2004 CS267 Lecture 15 8Relation of Poisson’s to Gravity, Electrostatics• Force on particle at (x,y,z) due to particle at 0 is-(x,y,z)/r^3, where r = sqrt(x +y +z )• Force is also gradient of potential V(x,y,z) = -1/rForce = -(d/dx V, d/dy V, d/dz V) = -grad V• V(x,y,z) satisfies Poisson’s equation (try it!)22 23/16/2004 CS267 Lecture 15 9Poisson’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”Discretize∂∂∂∂2u/∂∂∂∂x2= f(x) on regular mesh ui = u(i*h) to get[ u i+1–2*ui+ u i-1] / h2= f(x)Write as solving Tu = -h2* f for u where3/16/2004 CS267 Lecture 15 102D Poisson’s equation: ∂∂∂∂2u/∂∂∂∂x2+ ∂∂∂∂2u/∂∂∂∂y2 = f(x,y)• Similar to the 1D case, but the matrix T is now• Grid points numbered left to right, top row to bottom row• 3D is analogous• Similar “adjacency matrix” for arbitrary graph4 -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”3/16/2004 CS267 Lecture 15 11Algorithms 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 communication3/16/2004 CS267 Lecture 15 12Short explanations of Algorithms• Sorted in two orders (roughly):• from slowest to fastest on sequential machines• from most general (any matrix) to most specialized (matrices “like” Poisson)• Dense LU: Gaussian elimination; works on any N-by-N matrix• Band LU: exploit fact that T is nonzero only on sqrt(N) diagonals nearest main diagonal, so faster• Jacobi: essentially does matrix-vector multiply by T in inner loop of iterative algorithm. (May be used in Multigrid)• Explicit Inverse: assume we want to solve many systems with T, we can pre-compute and store inv(T) “amortized”, and just multiply by it. (Still expensive.)• Conjugate Gradients: uses matrix-vector multiplication, like Jacobi, but exploits mathematical properies of T that Jacobi does not• Red-Black SOR (Successive Overrelaxation): Variation of Jacobi that exploits yet different mathematical properties of T (May be used in Multigrid.)• Sparse LU: Gaussian elimination exploiting particular zero structure of T• FFT (Fast Fourier Transform): works only on matrices very like T• Multigrid: also works on matrices like T, that come from elliptic PDEs• Lower Bound: serial (time to print answer); parallel (time to combine N inputs)• Details in class notes and www.cs.berkeley.edu/~demmel/ma22133/16/2004 CS267 Lecture 15 13Comments on
View Full Document
Unlocking...