Prof David Bindel

Please click the play button below.

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

```
% 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

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
- Call LAPACK/ScaLAPACK!
- For \(n\) large: wait a while

```
% 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

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!

- You
*Order your unknowns*for sparsity- Again, good to use someone else’s software!

- For \(n\) large, 3D: get lots of memory and wait

```
% 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

- Dense: LAPACK, ScaLAPACK, PLAPACK, PLASMA, MAGMA
- Sparse direct: UMFPACK, TAUCS, SuperLU, MUMPS, Pardiso, SPOOLES, ...
- Sparse iterative: too many!

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

(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

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

Reading: Templates book