CS 5220

Applications of Parallel Computers

Intro to Sparse LA

Prof David Bindel

Please click the play button below.

World of Linear Algebra

  • Dense methods (last week)
  • Sparse direct methods (Thurs)
  • Iterative methods (today and Thurs)

Dense methods

% Dense (LAPACK)
[L,U] = lu(A);
x = U\(L\b);
  • Direct representation of matrices with simple data structures (no need for indexing data structure)
  • Mostly \(O(n^3)\) factorization algorithms

Software Strategies: Dense Case

Assuming you want to use (vs develop) dense LA code:

  • Learn enough to identify right algorithm
    (e.g. is it symmetric? definite? banded? etc)
  • Learn high-level organizational ideas
  • Make sure you have a good BLAS
  • For \(n\) large: wait a while

Sparse direct methods

% Sparse direct (UMFPACK + COLAMD)
[L,U,P,Q] = lu(A);
x = Q*(U\(L\(P*b)));
  • Direct representation, keep only the nonzeros
  • Factorization costs depend on problem structure (1D cheap; 2D reasonable; 3D gets expensive; not easy to give a general rule, and NP hard to order for optimal sparsity)
  • Robust, but hard to scale to large 3D problems

Sparse Direct Strategies

Assuming you want to use (vs develop) sparse LA code

  • Identify right algorithm (mainly Cholesky vs LU)
  • Get a good solver (often from list)
    • You don’t want to roll your own!
  • Order your unknowns for sparsity
    • Again, good to use someone else’s software!
  • For \(n\) large, 3D: get lots of memory and wait

Sparse iterations

% Sparse iterative (PCG + incomplete Cholesky)
tol = 1e-6;
maxit = 500;
R = cholinc(A,'0');
x = pcg(A,b,tol,maxit,R',R);
  • Only need \(y = Ax\) (maybe \(y = A^T x\))
  • Produce successively better (?) approximations
  • Good convergence depends on preconditioning
  • Best preconditioners are often hard to parallelize

LA Software: the Wider World

  • Sparse direct: UMFPACK, TAUCS, SuperLU, MUMPS, Pardiso, SPOOLES, ...
  • Sparse iterative: too many!

Sparse Iterative Software

Assuming you want to use (vs develop) sparse LA software...

  • Identify a good algorithm (GMRES? CG?)
  • Pick a good preconditioner
    • Often helps to know the application
    • ... and to know how the solvers work!
  • Play with parameters, preconditioner variants, etc...
  • Swear until you get acceptable convergence?
  • Repeat for the next variation on the problem

Stacking Solvers

(Typical) example from a bone modeling package:

  • Outer load stepping loop
  • Newton method corrector for each load step
  • Preconditioned CG for linear system
  • Multigrid preconditioner
  • Sparse direct solver for coarse-grid solve (UMFPACK)
  • LAPACK/BLAS under that

Overview of the week

  • Up next: Stationary iterations
  • Thurs: Krylov and sparse direct

Reading: Templates book