Iterative Methods in Linear Algebra part 1 Stan Tomov Innovative Computing Laboratory Computer Science Department The University of Tennessee Wednesday April 8 2009 CS 594 04 08 2009 CS 594 04 08 2009 Outline Part I Discussion Motivation for iterative methods Part II Stationary Iterative Methods Part III Nonstationary Iterative Methods CS 594 04 08 2009 Part I Discussion CS 594 04 08 2009 Discussion Part 1 CS 594 04 08 2009 Discussion Part 2 CS 594 04 08 2009 Discussion Original matrix load matrix output S spconvert matrix spy S Reordered matrix load A data S spconvert A spy S CS 594 04 08 2009 Discussion Original sub matrix load matrix output S spconvert matrix spy S 1 8660 1 8660 Reordered sub matrix load A data S spconvert A spy S 8602 8700 8602 8700 CS 594 04 08 2009 Discussion Note Columns have to be sorted CS 594 04 08 2009 Discussion Results on torc0 Pentium III 933 MHz 16 KB L1 and 256 KB L2 cache Results on woodstock Intel Xeon 5160 Woodcrest 3 0GHz L1 32 KB L2 4 MB Original matrix Mflop s 42 4568 Original matrix Mflop s 386 8 Reordered matrix Mflop s 45 3954 Reordered matrix Mflop s 387 6 BCRS Mflop s 72 17374 BCRS Mflop s 894 2 CS 594 04 08 2009 Homework 9 CS 594 04 08 2009 Questions Homework 9 PART I PETSc Many iterative solvers We will see how is projection used in iterative methods PART II hybrid CPU GPUs computations CS 594 04 08 2009 Hybrid CPU GPU computations CS 594 04 08 2009 Iterative Methods CS 594 04 08 2009 Motivation So far we have discussed and have seen Many engineering and physics simulations lead to sparse matrices e g PDE based simulations and their discretizations based on FEM FVM finite differences Petrov Galerkin type of conditions etc How to optimize performance on sparse matrix computations Some software how to solve sparse matrix problems e g PETSc The question is Can we solve sparse matrix problems faster than using direct sparse methods In certain cases Yes using iterative methods CS 594 04 08 2009 Sparse direct methods Sparse direct methods These are based on Gaussian elimination LU Performed in sparse format Are there limitations of the approach Yes they have fill ins which lead to more memory requirements more flops being performed Fill ins can become prohibitively high CS 594 04 08 2009 Sparse direct methods Consider LU for the matrix below a nonzero is represented by a 1st step of LU factorization will introduce fill ins marked by an F CS 594 04 08 2009 Sparse direct methods Fill ins can be improved by reordering Remember we talked about it in slightly different context for speed and parallelization Consider the following reordering These were extreme cases but still the problem exists CS 594 04 08 2009 Sparse methods Sparse direct vs dense methods Dense takes O n2 storage O n3 flops runs within peak performance Sparse direct can take O nonz storage and O nonz flops but these can also grow due to fill ins and performance is bad with nonz n2 and proper ordering it pays off to do sparse direct Software table from Julien Langou http www netlib org utk people JackDongarra la sw html Package DSCPACK HSL MFACT MUMPS PSPASES SPARSE SPOOLES SuperLU TAUCS UMFPACK Y12M Support yes yes yes yes yes yes yes yes Real X X X X X X X X X X X Type Complex f77 X X X X X X X X X X X X Language c c X X X X X X X X X Mode Seq Dist X M X X X M M X X M X M X X X SPD X X X X X X X X Gen X X X X X X X X CS 594 04 08 2009 Sparse methods What about Iterative Methods Think for example xi 1 xi P b Axi Pluses Storage is only O nonz for the matrix and a few working vectors Can take only a few iterations to converge e g n for P A 1 it takes 1 iteration check In general easy to parallelize CS 594 04 08 2009 Sparse methods What about Iterative Methods Think for example xi 1 xi P b Axi Minuses Performance can be bad as we saw in Lecture 11 and today s discussion Convergence can be slow or even stagnate But can be improved with preconditioning Think of P as a preconditioner an operator matrix P A 1 Optimal preconditioners e g multigrid can be lead to convergence in O 1 iterations CS 594 04 08 2009 Part II Stationary Iterative Methods CS 594 04 08 2009 Stationary Iterative Methods Can be expressed in the form xi 1 Bxi c where B and c do not depend on i older simpler easy to implement but usually not as effective as nonstationary examples Richardson Jacobi Gauss Seidel SOR etc next CS 594 04 08 2009 Richardson Method Richardson iteration xi 1 xi b Axi I A xi b 1 i e B from the definition above is B I A Denote ei x xi and rewrite 1 x xi 1 x xi Ax Axi ei 1 ei Aei I A ei ei 1 I A ei I A ei I A ei 1 2 i I A e0 i e for convergence ei 0 we need I A 1 for some norm e g when A 1 CS 594 04 08 2009 Jacobi Method Jacobi Method xi 1 xi D 1 b Axi I D 1 A xi D 1 b 2 where D is the diagonal of A assumed nonzero B I D 1 A from the definition Denote ei x xi and rewrite 2 1 x xi 1 x xi D ei 1 ei D 1 I D 1 ei 1 I D I D Ax Axi Aei A ei 1 A ei I D 1 1 A ei I D 1 2 A ei 1 i A e0 i e for convergence ei 0 we need I D 1 A 1 for some norm e g when D 1 A 1 CS 594 04 08 2009 Gauss Seidel Method Gauss Seidel Method xi 1 D L 1 Uxi b 3 where L is the lower and U the upper triangular part of A and D is the diagonal Equivalently 8 i 1 x1 i 1 x2 i 1 x3 i 1 xn i i i a1n xn a11 i a2n xn a22 i a3n xn a22 i 1 an n 1 xn 1 ann b1 a12 x2 a13 x3 a14 x4 i 1 i i b2 a21 x1 a23 x3 a24 x4 i 1 i 1 i b2 a31 x1 a32 x2 a34 x4 i 1 bn an1 x1 i 1 an2 x2 an3 x3 i i 1 CS 594 04 08 2009 A comparison from Julien Langou slides Gauss Seidel method 10 B 0 B A B 1 0 1 0 1 12 2 3 2 2 1 9 1 0 0 3 1 10 …
View Full Document
Unlocking...