Solver

BYOS

The numerical linear algebra solver that you use will ultimately be the bottleneck of your large scale inversion. To be the most flexible, SimPEG provides wrappers rather than a comprehensive set of solvers (i.e. BYOS).

The interface is as follows:

A # Where A is a sparse matrix (or linear operator)
Ainv = Solver(A, **solverOpts) # Create a solver object with key word arguments
x = Ainv * b # Where b is a numpy array of shape (n,) or (n,*)
Ainv.clean() # This cleans the memory footprint(if any)

Note

This is somewhat an abuse of notation for solvers as we never actually create A inverse. Instead we are creating an object that acts like A inverse, whether that be a Krylov subspace solver or an LU decomposition.

To wrap up solvers in scipy.sparse.linalg it takes one line of code:

Solver   = SolverWrapD(sp.linalg.spsolve, factorize=False)
SolverLU = SolverWrapD(sp.linalg.splu, factorize=True)
SolverCG = SolverWrapI(sp.linalg.cg)

Note

The above solvers are loaded into the base name space of SimPEG.

See also

• https://bitbucket.org/petsc/petsc4py

• https://github.com/bfroehle/pymumps

• https://github.com/rowanc1/pymatsolver

The API

SimPEG.utils.solver_utils.SolverWrapD(fun, factorize=True, checkAccuracy=True, accuracyTol=1e-06, name=None)

Wraps a direct Solver.

import scipy.sparse as sp
Solver   = solver_utils.SolverWrapD(sp.linalg.spsolve, factorize=False)
SolverLU = solver_utils.SolverWrapD(sp.linalg.splu, factorize=True)

SimPEG.utils.solver_utils.SolverWrapI(fun, checkAccuracy=True, accuracyTol=1e-05, name=None)

Wraps an iterative Solver.

import scipy.sparse as sp
SolverCG = solver_utils.SolverWrapI(sp.linalg.cg)