Iterative Methods for Solving Systems of Equations

1. CS 584

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 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. Consider • Are systems of equations and finite element methods related?

7. 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

8. 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

9. 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

10. 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; }

11. 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.

12. Inner Product • Suppose data is distributed row-wise • Inner product is simply dot product • IP = Sum(x[j] * x[j]) • This only requires a global sum collapse • O(log n)

13. 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(n log n)

14. 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)

15. 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

16. æ ö - - 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

17. 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.

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

19. 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.

20. Partitioning P0 P1 P2 P3

21. P0 P1 P2 P3 Communication

22. 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

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

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

25. 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