# Presentation: 233806

## Use of Parallel SuperLU in Large-Scale Scienticfic Calculations

- X. Sherry Li
- Lawrence Berkeley National Lab
- SIAM Annual Meeting, July 12-16, 2004

**Notes:**
25 minutes, 4:30-4:55, Tuesday

## Outline

- Algorithms and implementation in SuperLU_DIST
- Compare with sequential SuperLU
- In quantum mechanics (linear systems)
- With M. Baertschy, C. W. McCurdy, T. N. Rescigno, W. A. Isaacs
- In accelerator design (eigenvalue problems)
- With P. Husbands, C. Yang, L. Lee
- New developments

- Saint Girons

- 06/11/2003

## Introduction

- Direct method (Gaussian elimination) very useful for ill-conditioned systems
- Uses of direct solver in large-scale applications:
- Solve the system completely
- How accurate is the solution?
- Pivot growth, condition number, error bounds
- Preconditioner for iterative solvers
- linear system: M-1A x = M-1b
- shift-and-invert eigensolver: Kx = x  (K -  I)-1 x = x where  = 1 / ( - )

- Saint Girons

- 06/11/2003

- SuperLU (and _MT)
- [Demmel/Gilbert/Li]
- Partial pivoting
- With diagonal preference
- Sparsity order: A’*A-based
- Fix column order
- “BLAS–2.5”
- Left-looking (fan-in)
- More “read” than “write”

- SuperLU_DIST
- [Li/Demmel]
- “Static” pivoting
- Based on A
- Sparsity order: A’+A-based
- Fix row & column order
- BLAS–3
- Right-looking (fan-out)
- More parallelism
- 2D block-cyclic layout

- Saint Girons

- 06/11/2003

**Notes:**
“write” is more expensive for cache-coherent SMP
SuperLU_DIST was almost start from scratch

## How to distribute the matrices?

- Matrices involved:
- A, B (turned into X) – input, users manipulate them
- L, U – output, users do not need to see them
- A (sparse) and B (dense) are distributed by block rows
- 						Local A stored in
- 						***Compressed Row***
- Natural for users, and consistent with other popular packages: PETSc, Aztec, etc.

- A

- B

**x     x      x    **

**x     x      x**

**x      x          **

**x      x           **

**P0**

**P1**

**P2**

- Saint Girons

- 06/11/2003

## Distribution of L & U

- 2D (non-uniform) block cyclic distribution
- Better for GE scalability, load balance
- Provide re-distribution routine to change A from 1D block layout into 2D block cyclic layout

**3**

**1**

- Saint Girons

- 06/11/2003

## Quantum Mechanics

- Scattering in a quantum system of  three charged particles
- Simplest example is ionization of a hydrogen atom by collision with an electron:  e- + H  H+ + 2e-
- Seek the particles’ wave functions represented by the time-independent Schrodinger equation
- First solution to this long-standing unsolved problem [Recigno, et. al. *Science*, 24 Dec 1999]

- Saint Girons

- 06/11/2003

## Quantum Mechanics, cont.

- New *Exterior Complex Scaling* formulism to represent scattering states which simplifies boundary conditions
- Becomes computationally tractable
- Finite difference leads to complex, unsymmetric systems
- Diagonal blocks have the structure of 2D finite difference Laplacian matrices
- 	Very sparse: nonzeros per row <= 13
- Off-diagonal block is a diagonal matrix
- Between 6 to 24 blocks, each of order
- 	between 200K and 2M
- Total dimension up to 8.4 M
- Too much fill if use direct method . . .

- Saint Girons

- 06/11/2003

**Notes:**
Long range connectivity due to coupled diagonals  cannot find smaller separator  any ordering method fail

## Matrix Problem

- Example: N = 1024, NNZ = 12544, Cond. Num. 4e+03
- No iterative solvers without preconditioners or with simple preconditioners converge

- Saint Girons

- 06/11/2003

## SuperLU_DIST as Preconditioner

- Block-diagonal preconditioner for CGS iteration:
- 	M-1 A x = M-1 b
- 	M = diag(A11, A22, A33, ...)
- Use SuperLU_DIST for each
- 	diagonal block
- Need create 3 process grids,
- 	same logical ranks (0:3),
- 	different physical ranks
- Each grid has its own MPI communicator
- No pivoting, nor iterative refinement

- Saint Girons

- 06/11/2003

**Notes:**
Ax is fast; direct solver is the bottleneck.
96-280 iterations

## Timings on IBM SP

- Complex, unsymmetric
- N = 2 M, NNZ = 26 M
- Fill-ins using Metis: 1.3 G
- 	(50x fill)
- Factorization speed
- 10x speedup (4 to 128 P)
- Up to 30 Gflops

- Saint Girons

- 06/11/2003

## Accelerator Cavity Design

- Calculate cavity mode frequencies and field vectors
- Solve Maxwell equation in electromagnetic field
- Omega3P simulation code (SLAC)

- Omega3P model of a 47-cell section of the 206-cell Next Linear Collider accelerator structure

- Individual cells used in accelerator structure

- Saint Girons

- 06/11/2003

## Accelerator, cont.

- Finite element methods lead to large sparse generalized eigensystem  Kx = Mx
- Real symmetric for lossless cavities
- 	Complex symmetric when lossy in cavities
- Seek interior eigenvalues (tightly clustered) that are relatively small in magnitude

- Saint Girons

- 06/11/2003

## Accelerator, cont.

- Speed up Lanczos convergence by shift-invert
- 	 Seek largest eigenvalues, well separated, of the transformed system
- 		M (K -  M)-1 x =  M x
- 	 	 = 1 / ( - )
- The Filtering algorithm [Y. Sun]
- Inexact shift-invert Lanczos + JOCC
- Exact shift-invert Lanczos (ESIL)
- PARPACK for Lanczos
- SuperLU_DIST for shifted linear system
- No pivoting, nor iterative refinement

- Saint Girons

- 06/11/2003

## DDS47, Linear Elements

- Total eigensolver time: N = 1.3 M, NNZ = 20 M

- Saint Girons

- 06/11/2003

## Largest Problem

- DDS47, quadratic elements
- N = 7.5 M, NNZ = 304 M
- 6 G fill-ins using Metis
- 24 processors (8x3)
- Factor: 3,347 s
- 1 Solve: 61 s
- Eigensolver: 9,259 s  (~2.5 hrs)
- 1 shift, 10 eigenvalues, 55 solves

- Saint Girons

- 06/11/2003

**Notes:**
No result from Filtering algo
What’s the limit of ESIL?

## Parallel Symbolic Factorization[L. Grigori]

- Parallel ordering with ParMETIS on G(A’+A)
- Separator tree (binary) to guide computation
- Each step: one row of U, one column of L
- Within each separator: 1D block cyclic distribution
- Send necessary contribution to parent processor
- Preliminary results:
- Reasonable speedup: up to 6x
- Need improve memory balance

- Saint Girons

- 06/11/2003

**Notes:**
Disadvantage: tied with ParMETIS ordering
Integrate into SuperLU_DIST in half a years

## Improve Triangular Solve

- Important in preconditioning . . .
- “Switch-to-dense” towards the end
- Dense part can be of order
- Remove overhead of sparse data structure
-    and control structure
- On uniprocessors, achieved 1.8x speedup
-    together with register blocking [Vuduc et al.]
- Better speedup expected in parallel code

- Saint Girons

- 06/11/2003

## Improve Triangular Solve

- Partitioned inverse [Alvarado,Pothen,Schreiber]
- Preprocessing: find a partitioned representation of L
- Solve: m sparse matrix-vector multiply
- Requirements:
- Each     is invertible in place (no more fills generated)
- Find m as small as possible by proper reordering
- Trivial partitions:
- m = n,  i.e., each     is an elementary matrix from GE
- m = # of supernodes,  i.e., supernode partition

- Saint Girons

- 06/11/2003

## Extras

## GE with Static Pivoting

- Before factorization, scale and permute A to maximize diagonal: A  Pr Dr A Dc
- MC64 [Duff/Koster ‘99]
- Parallel auction algorithm on the way [Riedy]
- Find a permutation Pc to preserve sparsity:
- 	A  Pc A Pc’
- While factorizing A = LU, replace tiny pivots by
- 	||A||, without changing structures of L & U
- Most matrices insensitive to this perturbation
- If needed, use iterative refinement or an iterative solver to improve first solution

- Saint Girons

- 06/11/2003

## Numerical Results

- 68 unsymmetric matrices
- If no pivoting at all
- 26 contain zeros on diagonal to begin with
- 2 create new zeros on diagonal during elimination
- bbmat, orsreg_1
- Most of the others get large errors due to pivot growth
- Working precision iterative refinement (64 bits)
- Detailed table at www.nersc.gov/~xiaoye/SuperLU/GESP

- Saint Girons

- 06/11/2003

## Sparse Backward Error

- Saint Girons

- 06/11/2003

**Notes:**
Termination criterion

## Accuracy: GESP vs GEPP

- Saint Girons

- 06/11/2003

**Notes:**
Stopping criterion

## Refinement Steps

- Saint Girons

- 06/11/2003

## Times of Various Phases

- Saint Girons

- 06/11/2003

## Better Sparsity Ordering

- SuperLU_DIST: fill-ins in L+U (10^6)
- Diagonal Markowitz (based on A): 10-15% better
- [Amestoy/Li/Ng ’03]

- Saint Girons

- 06/11/2003

## BLAS 3 Kernel

- Copying is almost free
- 20-40% better than BLAS 2.5 on IBM SP

- Saint Girons

- 06/11/2003

## Default Setting May Not Be Good

- Diagonal perturbation is bad for 4 matrices
- Ex11, fidap011, inaccura, raefsky4
- Work well by either one of the 2 methods
- Turn off diagonal perturbation, or
- Use LU as preconditioner for GMRES [SPARSKIT]

- Saint Girons

- 06/11/2003

## SuperLU_DIST as Preconditioner

- SuperLU_DIST as block-diagonal preconditioner for CGS iteration
- 	M-1A x = M-1b
- 	M = diag(A11, A22, A33, ...)
- Run multiple SuperLU_DIST simultaneously
- No pivoting, nor iterative refinement
- 96 to 280 iterations @ 1 ~ 2 minute/iteration using 64 IBM SP processors
- Total time ~ 1 to 8 hours

- Saint Girons

- 06/11/2003

## SuperLU_DIST Enhancements

- 64-bit addressing
- Some integers need to be 64-bit for the largest problem (otherwise integer overflow)
- Will large memory SMP nodes help us?
- IBM SP nodes
- 	16-way 375MHz POWER3GB/node: 16, 32 (64 nodes), 64 (4 nodes) “Colony” switch

- Saint Girons

- 06/11/2003

**Notes:**
A description of how we improved SuperLU
Use of 64-bit
Use of Shared Memory
Change “Issues”

## Use of Shared Memory

- Symbolic factorization (7 GB for largest problem) still serial and replicated on all processors
- Save communication by not having to broadcast output
- Better to perform 1 symbfact() per SMP node and share the results [Husbands]
- Also used in High Performance Linpack

- Saint Girons

- 06/11/2003

## Some Open Problems

- Parameter setting is configurable by users
- Need an automatic strategy to choose a proper configuration
- Run the default configuration first (mostly works)
- Compute forward & backward error bounds
- If large errors, invoke an iterative solver
- How does extra-precision iterative refinement help?
- Investigating dense LU using XBLAS
- Guaranteed stability if variable precision is used during GE [Demmel]

- Saint Girons

- 06/11/2003

**Notes:**
Manual vs automatic
Software integration: PETSc, or other iterative package

## Future Work

- Parallel symbolic factorization to enhance memory scalability
- Better triangular solve (latency-bound)
- Include ILU capability
- Optimizations for SMP clusters
- Hybrid approach with threads in each MPI task

- Saint Girons

- 06/11/2003