# Presentation: 153128

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

- X. Sherry Li
- Lawrence Berkeley National Lab
- June 11, 2003

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

- Saint Girons

- 06/11/2003

- SuperLU (_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

## 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 ’01]

- 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

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

## 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 350K
- Total dimension up to 8.4 M
- Too much fill if use direct method

- Saint Girons

- 06/11/2003

**Notes:**
Long range connectivity

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

- 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

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

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

- Individual cells used in accelerating 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
- We added exact shift-invert Lanczos (ESIL)
- PARPACK for Lanczos
- SuperLU_DIST for shifted linear system
- No pivoting, nor iterative refinement

- 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

## 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)
- 10 eigenvalues, 1 shift, 55 solves
- Future: N=22.3 M, NNZ = 475 M, too big?

- Saint Girons

- 06/11/2003

**Notes:**
No result from Filtering algo

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