Linear System Solvers¶

Solving a system of linear equations is a key operation in many scientific applications. More specifically, for a matrix \(\mathbf{A} \in \mathbb{C}^{m \times n}\) and a vector \(\vec{b} \in \mathbb{C}^{m}\), we want to find a solution vector \(\vec{x} \in \mathbb{C}^{n}\) that satisfies

\[\mathbf{A} \vec{x} = \vec{b}\]

If \(m < n\) or \(m > n\), then the system do not have an exact solution, and thus, we are typically looking for the vector \(\vec{x}\) that minimize the norm \(\|\vec{b} - \mathbf{A} \vec{x} \|\) (i.e., the least-square solution). If the matrix \(\mathbf{A}\) is dense and relatively small, then we can use routines from the LAPACK to solve the system with a high degree of accuracy. However, if the matrix is sparse and/or very large, then we need to look for alternative numerical methods.

Sparse Direct Solvers¶

template<Number T> class Pardiso¶

PARDISO solver from Intel MKL. It finds the solution \( \mathbf{X} \) for the sparse linear system

\[ \mathbf{AX} = \mathbf{B}, \]

where \( \mathbf{A} \) is a complex or real square matrix, using the LU (or similar) factorization, i.e., \( \mathbf{PA} = \mathbf{LU}\) with \( \mathbf{P}\) as the permutation matrix to reduce the fill-in, \( \mathbf{L}\) as the lower triangular factor and \( \mathbf{U}\) as the upper triangular factor.

Usage:

CSRMatrix<double> A;
Vector<double> b;
Vector<double> x;
Pardiso<double> pardiso;
// ====== Initialize A, b ====== //
pardiso.factorize(A);                   // Compute the factorization PA = LU
x = pardiso.solve(b);                   // Solve the linear system using the LU factors

Warning

This solver is only available when using the Intel MKL backend.

Template Parameters:: T – Datatype used by the solver.

Public Functions

Pardiso()¶: Creates a new Pardiso solver.

Pardiso(const CSRMap<T, Host> &A)¶

Create a new Pardiso solver and then factorizes the matrix A.

Parameters:: A – [in] input matrix

virtual ~Pardiso()¶: Default destructor.

std::array<index_t, 64> &params()¶

Returns:: a reference to the parameters of the solver.

void factorize(const CSRMap<T, Host> &A)¶

Factorizes the matrix A as PA = LU.

Parameters:: A – [in] input matrix

PardisoVecResult solve(VectorMap<T, Host> b)¶

Solves the linear system \( \mathbf{A} \vec{x} = \vec{b} \) using the previous calculated P, L, U factors.

Parameters:: b – [in] right hand side

PardisoMatResult solve(MatrixMap<T, Host> B)¶

Solves the linear system \( \mathbf{AX} = \mathbf{B} \) using the previous calculated P, L, U factors.

Parameters:: B – [in] right hand side

PardisoVecResult adjoint_solve(VectorMap<T, Host> b)¶

Solves the linear system \( \mathbf{A}^* \vec{x} = \vec{b} \) using the previous calculated P, L, U factors.

Parameters:: b – [in] right hand side

PardisoMatResult adjoint_solve(MatrixMap<T, Host> B)¶

Solves the linear system \( \mathbf{A}^* \mathbf{X} = \mathbf{B} \) using the previous calculated P, L, U factors.

Parameters:: B – [in] right hand side

AMGCL¶

The iterative solvers are provided by AMGCL. AMGCL is a header-only C++ library for solving large sparse linear systems with algebraic multigrid (AMG) method. AMG is one of the most effective iterative methods for solution of equation systems arising, for example, from discretizing PDEs on unstructured grids. The method can be used as a black-box solver for various computational problems, since it does not require any information about the underlying geometry. AMG is often used not as a standalone solver but as a preconditioner within an iterative solver (e.g. Conjugate Gradients, BiCGStab, or GMRES). See the documentation of AMGCL here.

Usage:

using namespace HighFM;
CSRMatrix<double> A;
Vector<double> x, b;
//========== Initialize A, b ==============//

// Create a preconditioner based on the AMG, Smoothed Aggregation and SPAI0 methods
AMGSA<double, SPAI0> P(A);

// Create a solver using the GMRES algorithm
GMRES<double> gmres(A.rows());

// Solve the system
auto [steps, error] = solve_linsys(P, A, x, b, gmres);

API:

template<typename Backend> using ILU0 = amgcl::relaxation::ilu0<Backend>¶: Incomplete LU relaxation with zero fill-in. Used as parameter for the AMGClassic, AMGSA or RelaxPreconditioner.

template<typename Backend> using ILUT = amgcl::relaxation::ilut<Backend>¶: Incomplete LU relaxation that drops the matrix entries based on its magnitude. The nonzero structure is determined dynamically. Used as parameter for the AMGClassic, AMGSA or RelaxPreconditioner.

template<typename Backend> using SPAI0 = amgcl::relaxation::spai0<Backend>¶: Sparse approximate inverse (SPAI) with a diagonal matrix M. Used as parameter for the AMGClassic, AMGSA or RelaxPreconditioner.

template<typename Backend> using SPAI1 = amgcl::relaxation::spai1<Backend>¶: Sparse approximate inverse (SPAI) with a matrix M that has the same sparsity pattern as A. Used as parameter for the AMGClassic, AMGSA or RelaxPreconditioner.

template<Number T, template<typename> typename Relax, DataLocation L = Host> using AMGClassic = amgcl::amg<amgcl::backend::HighFMAdaptor<T, L>, amgcl::coarsening::ruge_stuben, Relax>¶

Algebraic Multigrid Preconditioner using direct Ruge-Stuben coarsening and direct interpolation.

Template Parameters:

T – Data type used by the solver.
Relax – Relaxation method.
L – Location of the data.

template<Number T, template<typename> typename Relax, DataLocation L = Host> using AMGSA = amgcl::amg<amgcl::backend::HighFMAdaptor<T, L>, amgcl::coarsening::smoothed_aggregation, Relax>¶

Algebraic Multigrid using a Smoothed Aggregation coarsening.

Template Parameters:

T – Data type used by the solver.
Relax – Relaxation method.
L – Location of the data.

template<Number T, DataLocation L = Host> using GMRESOptions = typename amgcl::solver::gmres<amgcl::backend::HighFMAdaptor<T, L>>::params¶: Parameters for the GMRES solver.

template<Number T, DataLocation L = Host> using GMRES = amgcl::solver::gmres<amgcl::backend::HighFMAdaptor<T, L>>¶

Generalized Minimal Residual Method (GMRES) is a popular iterative method for solving nonsymmetric linear systems. This solver uses restarting to limit computational and memory cost of the Arnoldi procedure.

Template Parameters:

T – Data type used by the solver.
L – Location of the data.

template<Number T, DataLocation L = Host> using ConjugateGradientOptions = typename amgcl::solver::cg<amgcl::backend::HighFMAdaptor<T, L>>::params¶: Parameters for the ConjugateGradient solver.

template<Number T, DataLocation L = Host> using ConjugateGradient = amgcl::solver::cg<amgcl::backend::HighFMAdaptor<T, L>>¶

Conjugate Gradient is an iterative method for solving linear whose coefficient matrix is positive-semidefinite (i.e., the matrix is symmetric/hermitian and all its eigenvalues are nonnegatives).

Template Parameters:

T – Data type used by the solver.
L – Location of the data.

template<Number T, DataLocation L, template<typename, typename> typename SA, template<typename, typename> typename SX, template<typename, typename> typename SB, typename Solver, typename Preconditioner> std::pair<index_t, double> solve_linsys(Preconditioner &&P, const CSRMatrix<T, L, SA> &A, Vector<T, L, SX> &x, const Vector<T, L, SB> &b, Solver &solver)¶

template<Number T, DataLocation L, template<typename, typename> typename SA, template<typename, typename> typename SX, template<typename, typename> typename SB, typename Solver, typename Preconditioner> std::pair<index_t, double> solve_linsys(Preconditioner &P, const CSRMatrix<T, L, SA> &A, Vector<T, L, SX> &x, Vector<T, L, SB> &b, Solver &solver)¶

template<Number T, DataLocation L, template<typename, typename> typename SA, template<typename, typename> typename SX, template<typename, typename> typename SB, typename Solver, typename Preconditioner> std::pair<index_t, double> solve_linsys(Preconditioner &&P, const CSRMatrix<T, L, SA> &A, Vector<T, L, SX> &x, Vector<T, L, SB> &b, Solver &&solver)¶

template<Number T, DataLocation L, template<typename, typename> typename SA, template<typename, typename> typename SX, template<typename, typename> typename SB, typename Solver, typename Preconditioner> std::pair<index_t, double> solve_linsys(Preconditioner &P, const CSRMatrix<T, L, SA> &A, Vector<T, L, SX> &x, Vector<T, L, SB> &b, Solver &&solver)¶

Solves the linear system Ax = b using the iterative solver and the preconditioner P from the amgcl package. On entry, the vector x must contain the initial conditions, and it will be overwritten by the solution of the linear system.

Parameters:

P – [inout] Preconditioner
A – [in] Coefficient matrix
x – [inout] The solution vector
b – [in] Right-hand side
solver – [inout] The solver used for finding the solution of the linear system

Returns:

a tuple containing the number of steps and an estimate of the error.

Least-Square¶

template<Number T, DataLocation L = Host> class LSMR¶

Iterative solver for least-square problems using the LSMR algorithm [1]. It solves the system of linear equations \( \mathbf{A}\vec{x} = \vec{b} \). If the system is inconsistent, it solves the least-square problem

\[ \min_{\vec{x}} \| \mathbf{A}\vec{x} - \vec{b} \|_2 \]

where A is a real rectangular matrix of dimension m-by-n and b is a real vector of length m. The matrix A may be dense or sparse (usually sparse).

To improve the convergence rate of the solver, you can specify a right preconditioner \( \mathbf{P}^{-1} \), such that now LSMR solves the following least square problem

\[ \min_{\vec{x}} \| (\mathbf{AP})^\dagger \; \mathbf{P} \vec{x} - \vec{b} \|_2 \]

More specifically, the Preconditioner object is a function that computes \( \vec{y} = \mathbf{P}^{-1} \mathbf{P}^{-T} \vec{x}\) for an input vector x and output vector y.

The code was adapted from the SciPy package (https://github.com/scipy/scipy.git).

References:

[1] D. C.-L. Fong and M. Saunders, ‘LSMR: An Iterative Algorithm for Sparse Least-Squares Problems’, SIAM J. Sci. Comput., vol. 33, no. 5, pp. 2950–2971, Jan. 2011, doi: 10.1137/10079687X.

Template Parameters:

T – Data type used by the solver. Must be either float or double.
L – Location of the data.

Public Functions

LSMR(LSMROptions options = LSMROptions())¶

Creates a new LSMR solver.

Parameters:: options – [in] solver parameters. See LSMROptions for more information.

virtual ~LSMR() = default¶: Default destructor.

inline LSMRMatResult operator()(MatrixMap<T, L> A, VectorMap<T, L> b, std::optional<VectorMap<T, L>> x0 = std::nullopt, std::optional<Preconditioner> precond = std::nullopt)¶

Solves the least square problem Ax = b.

Parameters:

A – [in] coefficient matrix
b – [in] right-hand side of the equation
x0 – [in] initial guess for the solution (optional)
precond – [in] a Preconditioner to the system (optional). Suppose that \( \mathbf{L}^{-1}\) is the preconditioner of system, then precond must apply \( \mathbf{L}^{-1} \mathbf{L}^{-T}\) over a vector x.

Returns:

An abstract object representing the solution of the linear system.

LSMRReport report()¶

Returns:: Some relevant information about the last run. See LSMRReport for more information.

void set_params(LSMROptions options)¶

Set the parameters of the solver.

Parameters:: options – [in] solver parameters. See LSMROptions for more information.

struct LSMROptions¶

Parameters of the LSMR solver.

Public Members

double damping_factor = 0¶

Damping factor \( \lambda \) for regularized least-squares problems. If \( \lambda \) is not 0, LSMR solves the regularized least-squares problem::

\[\begin{split} \min_x \| \begin{bmatrix} \vec{b} \\ 0 \end{bmatrix} - \begin{bmatrix} \mathbf{A}\vec{x} \\ \lambda \mathbf{I} \end{bmatrix} \|_2 \end{split}\]

double atol = 1E-6¶: Stopping criteria. If the system Ax = b seems to be consistent, LSMR terminates when \( \| \vec{r}\|_2 \leq atol * \| \mathbf{A}\vec{x} \|_2 + btol * \| \vec{b}\|\). Otherwise, it terminates when \( \left \| \mathbf{A}^H \vec{r}\|_2 \leq atol * \| \mathbf{A}\vec{r} \right \|_2 \).

double cond_limit = 1E8¶: LSMR terminates when the estimate for cond(A) exceeds cond_limit.

index_t max_iter = kAuto ¶: LSMR terminates when the number of iterations reaches max_iter. The default is max_iter = min(m, n).

int verbose = 0¶

Set the verbosity level. 0 - No report. 1 - Report only the residual

\[ \| \mathbf{A}\vec{x} - \vec{b} \|_2 \]

2 - Full iteration log.

struct LSMRReport¶

Report of the last run of the LSMR solver.

Public Members

std::string stop_condition¶: Stop condition.

index_t num_iter¶: Number of iterations.

double cond_A¶: Estimated condition number of the matrix A

double norm_residual¶: Norm of the residual \( \| \vec{b} - \mathbf{A}\vec{x}\|_2 \).