**TOPS' Sparse Direct Solver for Ill-Conditioned Systems**

**PIs:** X.-S. Li^1,2^, J. W. Demmel^2^, A. Pinar^1^, E. G. Ng^1^

^1^Lawrence Berkeley National Lab, ^2^U. California-Berkeley

Summary

***Sparse systems of linear equations arise at the heart of many
large-scale simulations. The Terascale Optimal PDE Simulations (TOPS)
project is providing high-performance sparse direct solvers, which have
had significant impacts on an accelerator physics application and a
fusion project.***

The kernels of many large-scale simulations often require the solution
of large sparse systems of linear equations. Some of these linear
systems can be very ill-conditioned and will require direct, as opposed
to iterative, solution.

SuperLU is a leading scalable solver for sparse linear systems using
direct methods. It is especially targeted for nonsymmetric, indefinite
problems. The library routines perform an LU factorization with
numerical pivoting and triangular solutions through forward and back
substitutions. The factorization can be applied to non-square matrices.
In addition, the driver routines contain the functionalities of
equilibrating the matrix, reordering the rows and columns of the matrix
for stability and sparsity, iterative refinement, estimating the
condition number, and computing the forward and backward error bounds.
The solver supports both real and complex data types.

The collection of routines has been carefully designed for optimal
performance on modern architectures. The sequential solver has achieved
up to 40% of the peak megaflop rate on many hierarchical memory
machines. The parallel solver has also achieved high performance, up to
a hundred-fold speedup on large matrices.

SuperLU is widely adopted in research and commercial use, with thousands
of users. It has been integrated into other TOPS solvers, such as PETSc
and Hypre. DOE application codes that use SuperLU include NIMROD,
Omega3P, DSpice, NIKE, Trilinos, and several computational plasma
physics projects at PPPL. SuperLU is used commercially in Mathematica,
FEMLAB, Python, and the HP Mathematical Library, among others.

Several significant scientific achievements have been made with the help
of SuperLU.

One achievement is the solution of a long-standing problem of scattering
in a quantum system of three charged particles. Parallel SuperLU played
a critical role in building the preconditioners for solving the complex,
nonsymmetric, and very ill-conditioned linear systems, up to the order
of 8 million, enabling a breakthrough result that was reported in a
cover article of *Science* (24 December 1999).

The second achievement is in the design of accelerators, which requires
large-scale eigenanalysis in electromagnetic modal analysis calculations
using the Omega3P code. SuperLU is paired with PARPACK, a parallel
implementation of the Implicitly Restarted Arnoldi Method, to obtain a
parallel implementation of the exact shift-and-invert algorithm (ESIL).
This ESIL solver accurately computed the selected interior eigenvalues
and vectors of the generalized eigenvalue problems and is usually faster
than the algorithm previously used in Omega3P. The largest system solved
is of order 7.5 million, with 300 million nonzeros.

A recent noteworthy achievement was made possible by the integration of
SuperLU into the NIMROD production code. NIMROD is used to simulate
fusion reactor plasmas by solving the continuum equations of
magnetohydrodynamics, using discretizations of higher than second order
on an unstructured mesh. The algorithm requires the solution of several
large sparse systems in parallel at every time step in a simulation. The
stiffness inherent in the physical system leads to the desire to solve
implicit systems from which the fastest waves are filtered. The implicit
systems suffer from elliptic ill conditioning and they generally lack
"M-matrix" structure due to the higher order discretization. These
properties represent challenges for traditional incomplete
factorization-preconditioned Krylov methods. For two-dimensional linear
calculations of MHD instabilities, NIMROD runs two orders of magnitude
faster with SuperLU than it does with the previously employed
preconditioned conjugate gradient solver. For cutting-edge
three-dimensional, nonlinear tokamak simulations (Figure 1), which
require supercomputing resources, NIMROD with SuperLU runs four to five
times faster than before. This improved efficiency translates directly
into four to five times more physics results for the same amount of time
on the computer --- a major improvement in scientific productivity. It
also makes it possible to contemplate scientifically necessary increases
in grid resolution.

![](media/image1.wmf){width="3.0631944444444446in"
height="2.1465277777777776in"}

*Figure 1. Full 3D numerical simulation of plasma particle drift orbits
in a tokamak.*

Future work on SuperLU includes improving the performance of the
triangular solution phase, parallelizing the symbolic factorization
algorithm, and leveraging the memory hierarchy tuning techniques being
developed in TOPS.

The performance of SuperLU depends in part on the sparsity of the
triangular factors. Thus, an important piece of ongoing work is the
development of ordering algorithms, which find column and/or row
permutations that will keep the triangular factors sparse. Both greedy
(such as minimum-degree and minimum-deficiency) algorithms and nested
dissection algorithms are under investigation. Recent work includes the
development of a nested dissection algorithm specifically for
nonsymmetric matrices.

The TOPS project webpage may be found at
[http://www.tops-scidac.org](http://www.tops-scidac.org/).

**For further information on this subject contact:**

Prof. David E. Keyes, Project Lead

Columbia University

Phone: 212-854-1120

david.keyes@columbia.edu