Numerical Linear Algebra and Fast Algorithms

Solutions of most practical numerical problems reduce to large matrix computations. They usually include the solutions of linear systems or eigen-value problems. Geophysics applications pose some unique challenges to modern direct and iterative solvers. Such challenges include high dimensions, ill conditioning, a large number of frequencies and sources (right-hand sides), modest accuracies, etc.

GMIG is known for its pioneering work on fast structured direct solvers, which produces approximate solutions with controllable accuracies and are suitable for multiple right-hand sides. Such solvers are applicable to Helmholtz equations in 2D and 3D, as well as more general sparse problems on irregular meshes. We show the existence of rank structures in the intermediate matrices in direct factorization or inversion, and then solve the problems with hierarchical semiseparable (HSS) methods integrated into a multifrontal framework. Our structured solvers have about O(n) and O(n) ~ O(n^4/3)complexity for 2D and 3D problems, respectively, with nearly O(n) memory requirement and nearly O(n) solution complexity. In comparison, standard factorizations cost O(n^3/2) in 2D and O(n²) in 3D, with O(n^4/3) memory in 3D.

Our methods involve multiple layers of hierarchical tree structures and advanced graph and sparse matrix techniques, and are highly scalable. They have been tested for problems with billions of variables on large clusters including NERSC's Cray XE6 through collaborations with LBL.

The solvers are used in the various seismic computations, especially 3D modeling of seismic wave propagation. Examples include time-harmonic elastic waves in anisotropic media and time-harmonic q^P-polarized waves in TTI media. Some of the work is conducted together with our industrial sponsors.

Our recent developments based on randomized sampling and selected inversion bring new insights into the factorization techniques, as well as pre-condition Gauss-Newton iterations and the computation of normal modes. New theories and methods are also developed for high dimensions and high frequencies.