1 / 30

Iterative Methods for Solving Systems of Equations

Iterative Methods, Gaussian Elimination, Indirect Methods, Convergence, Accuracy, Sparse Matrix, Stable, Jacobi Iteration, Gauss-Seidel Iteration, Parallelization

mgrisham
Download Presentation

Iterative Methods for Solving Systems of Equations

An Image/Link below is provided (as is) to download presentation Download Policy: Content on the Website is provided to you AS IS for your information and personal use and may not be sold / licensed / shared on other websites without getting consent from its author. Content is provided to you AS IS for your information and personal use only. Download presentation by click this link. While downloading, if for some reason you are not able to download a presentation, the publisher may have deleted the file from their server. During download, if you can't get a presentation, the file might be deleted by the publisher.

E N D

Presentation Transcript


  1. CS 484

  2. Iterative Methods • Gaussian elimination is considered to be a direct method to solve a system. • An indirect method produces a sequence of values that converges to the solution of the system. • Computation is halted in an indirect algorithm when a specified accuracy is reached.

  3. Why Iterative Methods? • Sometimes we don't need to be exact. • Input has inaccuracy, etc. • Only requires a few iterations • If the system is sparse, the matrix can be stored in a different format. • Iterative methods are usually stable. • Errors are dampened

  4. 7 -6 x1 3 -8 9 x2 -4 = Iterative Methods • Consider the following system. • Now write out the equations • solve for the ith unknown in the ith equation x1 = 6/7 x2 + 3/7 x2 = 8/9 x1 - 4/9

  5. Iterative Methods • Come up with an initial guess for each xi • Hopefully, the equations will produce better values and which can be used to calculate better values, etc. and will converge to the answer.

  6. Iterative Methods • Jacobi iteration • Use all old values to compute new values k x1 x2 0 0.00000 0.00000 10 0.14865 -0.19820 20 0.18682 -0.24908 30 0.19662 -0.26215 40 0.19913 -0.26551 50 0.19977 -0.26637

  7. Jacobi Iteration • The ith equation has the form - 1 n å = A [ i , j ] x [ j ] b [ i ] = 0 j • Which can be rewritten as: æ ö 1 å ç ÷ = - x [ i ] b [ i ] A [ i , j ] x [ j ] ç ÷ A [ i , i ] è ø ¹ j i

  8. r [ i ] = + x [ i ] x [ i ] A [ i , i ] Jacobi Iteration • The vector (b - Ax) is zero when and if we get the exact answer. • Define this vector to be the residual r • Now rewrite the solution equation

  9. void Jacobi(float A[][], float b[], float x[], float epsilon) { int k = 0; float x1[]; float r[]; float r0norm; // Randomly select an initial x vector r = b - Ax; // This involves matrix-vector mult etc. r0norm = ||r||2; // This is basically the magnitude while (||r||2 > epsilon * r0norm) { for (j = 0; j < n; j++) x1[j] = r[j] / A[j,j] + x[j]; r = b - Ax; } x = x1; }

  10. Parallelization of Jacobi • 3 main computations per iteration • Inner product (2 norm) • Loop calculating x[j]s • Matrix-vector mult. to calculate r • If A[j,j], b[j], & r[j] are on the same proc. • Loop requires no communication • Inner product and Matrix-vector mult require communication.

  11. Inner Product (2 norm of residual) • Suppose data is distributed row-wise • Inner product (2 norm) is simply dot product • IP = Sum(x[j] * x[j]) • This only requires a global sum collapse • O(log p) A x = b P0 P1 P2 P3

  12. Matrix-Vector Multiplication • Again data is distributed row-wise • Each proc. requires all of the elements in the vector to calculate their part of the resulting answer. • This results in all to all gather • O(p log p)

  13. Jacobi Iteration • Resulting cost for float (4 bytes) • Tcomm = #iterations * (TIP + TMVM) • TIP = log p * (ts + tw * 4) • TMVM = p log p * (ts + tw * nrows/p * 4)

  14. Iterative Methods • Gauss-Seidel • Use the new values as soon as available k x1 x2 0 0.00000 0.00000 10 0.21977 -0.24909 20 0.20130 -0.26531 30 0.20008 -0.26659 40 0.20000 -0.26666

  15. æ ö - - 1 1 i n 1 ç ÷ å å = - - x [ i ] b [ i ] x [ j ] A [ i , j ] x [ j ] A [ i , j ] ç ÷ - 1 k k k A [ i , i ] è ø = = + 0 1 j j i Gauss-Seidel Iteration • The basic Gauss-Seidel iteration is

  16. Gauss-Seidel Iteration • Rule: Always use the newest values of the previously computed variables. • Problem: Sequential? • Gauss-Seidel is indeed sequential if the matrix is dense. • Parallelism is a function of the sparsity and ordering of the equation matrix.

  17. Gauss-Seidel Iteration • We can increase possible parallelism by changing the numbering of a system.

  18. Parallelizing Red-Black GSI • Partitioning? Block checkerboard. • Communication? 2 phases per iteration. 1- compute red cells using values from black cells 2- compute black cells using values from red cells Communication is required for each phase.

  19. Partitioning P0 P1 P2 P3

  20. P0 P1 P2 P3 Communication

  21. Procedure Gauss-SeidelRedBlack while ( error > limit ) send black values to neighbors recv black values from neighbors compute red values send red values to neighbors recv red values from neighbors compute black values compute error /* only do every so often */ endwhile

  22. Extending Red-Black Coloring • Goal: Produce a graph coloring scheme such that no node has a neighbor of the same color. • Simple finite element and finite difference methods produce graphs with only 4 neighbors. • Two colors suffice • What about more complex graphs?

  23. More complex graphs • Use graph coloring heuristics. • Number the nodes one color at a time.

  24. Successive Overrelaxation • Devised to speed up the convergence of Gauss-Seidel • apply extrapolation using a weighted average between current & previous iterates • Choose weighting that will accelerate the rate of convergence

  25. æ æ ö ö - - - - 1 1 1 1 i i n n 1 1 ç ç ÷ ÷ å å å å = = - - - - ∂ x [ [ i i ] ] b b [ [ i i ] ] x x [ [ j j ] ] A A [ [ i i , , j j ] ] x x [ [ j j ] ] A A [ [ i i , , j j ] ] ç ç ÷ ÷ - - 1 1 k k k k k A A [ [ i i , , i i ] ] è è ø ø = = = = + + 0 0 1 1 j j j j i i SOR • Gauss-Seidel iteration • Don’t compute directly into x • Compute weighted average

  26. SOR • SOR Equation • Choose w in the range [0, 2] • technically w < 1 is underrelaxation • if you choose w > 2 it won’t converge • if you choose w incorrectly it won’t converge

  27. SOR Algorithm Choose an initial guess for x[i] for k = 1,2….. for i = 1,2, … n ∂ = 0 for j = 1,2,…., i-1 ∂ = ∂ + a[i,j] * xk[j] end for j = i+1, … , n ∂ = ∂ + a[i,j] * xk-1[j] end ∂ = (b[i] - ∂) / a[i,i] xk[i] = xk-1[i] + w * (∂ - xk-1[i]) end check for convergence end

  28. Parallelizing SOR • Just like Gauss-Seidel • Create the dependency graph • Color • Number by color • Phases • Communicate nodes of previous phase • Compute nodes of this phase • Move to next phase

  29. Conclusion • Iterative methods are used when an exact answer is not computable or needed. • Gauss-Seidel converges faster than Jacobi but parallelism is trickier. • Finite element codes are simply systems of equations • Solve with either Jacobi or Gauss-Seidel

  30. Consider • Are systems of equations and finite element methods related?

More Related