Linear Solver

Algorithms to solve linear systems of equations. More...

Namespaces
namespace	LinearSolver
	Functionality for solving linear systems of equations.
namespace	LinearSolver::Tags
	The DataBox tags associated with the linear solver.
namespace	LinearSolver::InnerProductImpls
	Implementations of LinearSolver::inner_product.

Classes
struct	LinearSolver::cg::ConjugateGradient< Metavariables, FieldsTag, OptionsGroup, SourceTag >
	A conjugate gradient solver for linear systems of equations \(Ax=b\) where the operator \(A\) is symmetric. More...
struct	LinearSolver::gmres::Gmres< Metavariables, FieldsTag, OptionsGroup, Preconditioned, SourceTag, ArraySectionIdTag >
	A GMRES solver for nonsymmetric linear systems of equations \(Ax=b\). More...
struct	LinearSolver::Richardson::Richardson< FieldsTag, OptionsGroup, SourceTag, ArraySectionIdTag >
	A simple Richardson scheme for solving a system of linear equations \(Ax=b\). More...
struct	LinearSolver::Schwarz::Schwarz< FieldsTag, OptionsGroup, SubdomainOperator, SubdomainPreconditioners, SourceTag, ArraySectionIdTag >
	An additive Schwarz subdomain solver for linear systems of equations \(Ax=b\). More...

Functions
template<typename Lhs, typename Rhs>
auto	LinearSolver::inner_product (const Lhs &lhs, const Rhs &rhs)
	The local part of the Euclidean inner product on the vector space w.r.t. which the addition and scalar multiplication of both Lhs and Rhs is defined.
template<typename T>
double	LinearSolver::magnitude_square (const T &vector)
	The local part of the Euclidean inner product of a vector with itself.
int	lapack::general_matrix_linear_solve (gsl::not_null< DataVector * > rhs_in_solution_out, gsl::not_null< Matrix * > matrix_operator, int number_of_rhs=0)
	Wrapper for LAPACK dgesv, which solves the general linear equation \(A x = b\) by LUP (Lower-triangular, upper-triangular, and permutation) decomposition.
int	lapack::general_matrix_linear_solve (gsl::not_null< DataVector * > rhs_in_solution_out, gsl::not_null< std::vector< int > * > pivots, gsl::not_null< Matrix * > matrix_operator, int number_of_rhs=0)
	Wrapper for LAPACK dgesv, which solves the general linear equation \(A x = b\) by LUP (Lower-triangular, upper-triangular, and permutation) decomposition.
int	lapack::general_matrix_linear_solve (gsl::not_null< DataVector * > solution, gsl::not_null< Matrix * > matrix_operator, const DataVector &rhs, int number_of_rhs=0)
	Wrapper for LAPACK dgesv, which solves the general linear equation \(A x = b\) by LUP (Lower-triangular, upper-triangular, and permutation) decomposition.
int	lapack::general_matrix_linear_solve (gsl::not_null< DataVector * > solution, gsl::not_null< std::vector< int > * > pivots, gsl::not_null< Matrix * > matrix_operator, const DataVector &rhs, int number_of_rhs=0)
	Wrapper for LAPACK dgesv, which solves the general linear equation \(A x = b\) by LUP (Lower-triangular, upper-triangular, and permutation) decomposition.
int	lapack::general_matrix_linear_solve (gsl::not_null< DataVector * > rhs_in_solution_out, const Matrix &matrix_operator, int number_of_rhs=0)
	Wrapper for LAPACK dgesv, which solves the general linear equation \(A x = b\) by LUP (Lower-triangular, upper-triangular, and permutation) decomposition.
int	lapack::general_matrix_linear_solve (gsl::not_null< DataVector * > solution, const Matrix &matrix_operator, const DataVector &rhs, int number_of_rhs=0)
	Wrapper for LAPACK dgesv, which solves the general linear equation \(A x = b\) by LUP (Lower-triangular, upper-triangular, and permutation) decomposition.

Detailed Description

Algorithms to solve linear systems of equations.

Details

In a way, the linear solver is for elliptic systems what time stepping is for the evolution code. This is because the DG scheme for an elliptic system reduces to a linear system of equations of the type \(Ax=b\), where \(A\) is a global matrix representing the DG discretization of the problem. Since this is one equation for each node in the computational domain it becomes unfeasible to numerically invert the global matrix \(A\). Instead, we solve the problem iteratively so that we never need to construct \(A\) globally but only need \(Ax\) that can be evaluated locally by virtue of the DG formulation. This action of the operator is what we have to supply in each step of the iterative algorithms implemented here. It is where most of the computational cost goes and usually involves computing a volume contribution for each element and communicating fluxes with neighboring elements. Since the iterative algorithms typically scale badly with increasing grid size, a preconditioner \(P\) is needed in order to make \(P^{-1}A\) easier to invert.

Note

The smallest possible residual magnitude the linear solver can reach is the product between the machine epsilon and the condition number of the linear operator that is being inverted. Smaller residuals are numerical artifacts. Requiring an absolute or relative residual below this limit will likely make the linear solver run until it reaches its maximum number of iterations.

Remember that when the linear operator \(A\) corresponds to a PDE discretization, decreasing the linear solver residual below the discretization error will not improve the numerical solution any further. I.e. the error \(e_k=x_k-x_\mathrm{analytic}\) to an analytic solution will be dominated by the linear solver residual at first, but even if the discretization \(Ax_k=b\) was exactly solved after some iteration \(k\), the discretization residual \(Ae_k=b-Ax_\mathrm{analytic}=r_\mathrm{discretization}\) would still remain. Therefore, ideally choose the absolute or relative residual criteria based on an estimate of the discretization residual.

In the iterative algorithms we usually don't work with the physical field \(x\) directly. Instead we need to apply the operator to an internal variable defined by the respective algorithm. This variable is exposed as the LinearSolver::Tags::Operand prefix, and the algorithm expects that the computed operator action is written into db::add_tag_prefix<LinearSolver::Tags::OperatorAppliedTo, LinearSolver::Tags::Operand<...>> in each step.

Function Documentation

general_matrix_linear_solve() [1/6]

int lapack::general_matrix_linear_solve	(	gsl::not_null< DataVector * >	rhs_in_solution_out,
		const Matrix &	matrix_operator,
		int	number_of_rhs = 0 )

Wrapper for LAPACK dgesv, which solves the general linear equation \(A x = b\) by LUP (Lower-triangular, upper-triangular, and permutation) decomposition.

Details

Several interfaces are provided with combinations of input parameters due to the versatility of the LAPACK utility.

• rhs_in_solution_out or solution : If the rhs DataVector is not separately supplied, this single pass-by-pointer DataVector acts as the input \( b\), and is overwritten with the solution \(x\); in this case rhs_in_solution_out must be sufficiently large to contain the solution. If rhs is separately supplied, solution will contain \(x\) at the end of the algorithm, and must be sufficiently large to contain the solution.
• rhs (optional) : An input DataVector representing \(b\) that is not overwritten.
• matrix_operator: The Matrix \(A\). If passed by pointer, this matrix will be 'destroyed' and overwritten with LU decomposition information. Optionally, a pivots vector may be passed by pointer so that the full LUP decomposition information can be recovered from the LAPACK output. If passed by const reference, the wrapper will make a copy so that LAPACK does not overwrite the supplied matrix.
• pivots (optional) : a std::vector<int>, passed by pointer, which will contain the permutation information necessary to reassemble the full matrix information if LAPACK is permitted to modify the input matrix in-place.
• number_of_rhs : The number of columns of \(x\) and \(b\). This is thought of as the 'number of right-hand-sides' to be solved for. In most cases, that additional dimension can be inferred from the dimensionality of the matrix and the input vector(s), which occurs if number_of_rhs is set to 0 (default). If the provided DataVector is an inappropriate size for that inference (i.e. is not a multiple of the largest dimension of the matrix), or if you only want to solve for the first number_of_rhs columns, the parameter number_of_rhs must be supplied.

The function return int is the value provided by the INFO field of the LAPACK call. It is 0 for a successful linear solve, and nonzero values code for types of failures of the algorithm. See LAPACK documentation for further details about dgesv: http://www.netlib.org/lapack/

general_matrix_linear_solve() [2/6]

int lapack::general_matrix_linear_solve	(	gsl::not_null< DataVector * >	rhs_in_solution_out,
		gsl::not_null< Matrix * >	matrix_operator,
		int	number_of_rhs = 0 )

Wrapper for LAPACK dgesv, which solves the general linear equation \(A x = b\) by LUP (Lower-triangular, upper-triangular, and permutation) decomposition.

Details

Several interfaces are provided with combinations of input parameters due to the versatility of the LAPACK utility.

• rhs_in_solution_out or solution : If the rhs DataVector is not separately supplied, this single pass-by-pointer DataVector acts as the input \( b\), and is overwritten with the solution \(x\); in this case rhs_in_solution_out must be sufficiently large to contain the solution. If rhs is separately supplied, solution will contain \(x\) at the end of the algorithm, and must be sufficiently large to contain the solution.
• rhs (optional) : An input DataVector representing \(b\) that is not overwritten.
• matrix_operator: The Matrix \(A\). If passed by pointer, this matrix will be 'destroyed' and overwritten with LU decomposition information. Optionally, a pivots vector may be passed by pointer so that the full LUP decomposition information can be recovered from the LAPACK output. If passed by const reference, the wrapper will make a copy so that LAPACK does not overwrite the supplied matrix.
• pivots (optional) : a std::vector<int>, passed by pointer, which will contain the permutation information necessary to reassemble the full matrix information if LAPACK is permitted to modify the input matrix in-place.
• number_of_rhs : The number of columns of \(x\) and \(b\). This is thought of as the 'number of right-hand-sides' to be solved for. In most cases, that additional dimension can be inferred from the dimensionality of the matrix and the input vector(s), which occurs if number_of_rhs is set to 0 (default). If the provided DataVector is an inappropriate size for that inference (i.e. is not a multiple of the largest dimension of the matrix), or if you only want to solve for the first number_of_rhs columns, the parameter number_of_rhs must be supplied.

The function return int is the value provided by the INFO field of the LAPACK call. It is 0 for a successful linear solve, and nonzero values code for types of failures of the algorithm. See LAPACK documentation for further details about dgesv: http://www.netlib.org/lapack/

general_matrix_linear_solve() [3/6]

int lapack::general_matrix_linear_solve	(	gsl::not_null< DataVector * >	rhs_in_solution_out,
		gsl::not_null< std::vector< int > * >	pivots,
		gsl::not_null< Matrix * >	matrix_operator,
		int	number_of_rhs = 0 )

Wrapper for LAPACK dgesv, which solves the general linear equation \(A x = b\) by LUP (Lower-triangular, upper-triangular, and permutation) decomposition.

Details

Several interfaces are provided with combinations of input parameters due to the versatility of the LAPACK utility.

• rhs_in_solution_out or solution : If the rhs DataVector is not separately supplied, this single pass-by-pointer DataVector acts as the input \( b\), and is overwritten with the solution \(x\); in this case rhs_in_solution_out must be sufficiently large to contain the solution. If rhs is separately supplied, solution will contain \(x\) at the end of the algorithm, and must be sufficiently large to contain the solution.
• rhs (optional) : An input DataVector representing \(b\) that is not overwritten.
• matrix_operator: The Matrix \(A\). If passed by pointer, this matrix will be 'destroyed' and overwritten with LU decomposition information. Optionally, a pivots vector may be passed by pointer so that the full LUP decomposition information can be recovered from the LAPACK output. If passed by const reference, the wrapper will make a copy so that LAPACK does not overwrite the supplied matrix.
• pivots (optional) : a std::vector<int>, passed by pointer, which will contain the permutation information necessary to reassemble the full matrix information if LAPACK is permitted to modify the input matrix in-place.
• number_of_rhs : The number of columns of \(x\) and \(b\). This is thought of as the 'number of right-hand-sides' to be solved for. In most cases, that additional dimension can be inferred from the dimensionality of the matrix and the input vector(s), which occurs if number_of_rhs is set to 0 (default). If the provided DataVector is an inappropriate size for that inference (i.e. is not a multiple of the largest dimension of the matrix), or if you only want to solve for the first number_of_rhs columns, the parameter number_of_rhs must be supplied.

The function return int is the value provided by the INFO field of the LAPACK call. It is 0 for a successful linear solve, and nonzero values code for types of failures of the algorithm. See LAPACK documentation for further details about dgesv: http://www.netlib.org/lapack/

general_matrix_linear_solve() [4/6]

int lapack::general_matrix_linear_solve	(	gsl::not_null< DataVector * >	solution,
		const Matrix &	matrix_operator,
		const DataVector &	rhs,
		int	number_of_rhs = 0 )

Wrapper for LAPACK dgesv, which solves the general linear equation \(A x = b\) by LUP (Lower-triangular, upper-triangular, and permutation) decomposition.

Details

Several interfaces are provided with combinations of input parameters due to the versatility of the LAPACK utility.

• rhs_in_solution_out or solution : If the rhs DataVector is not separately supplied, this single pass-by-pointer DataVector acts as the input \( b\), and is overwritten with the solution \(x\); in this case rhs_in_solution_out must be sufficiently large to contain the solution. If rhs is separately supplied, solution will contain \(x\) at the end of the algorithm, and must be sufficiently large to contain the solution.
• rhs (optional) : An input DataVector representing \(b\) that is not overwritten.
• matrix_operator: The Matrix \(A\). If passed by pointer, this matrix will be 'destroyed' and overwritten with LU decomposition information. Optionally, a pivots vector may be passed by pointer so that the full LUP decomposition information can be recovered from the LAPACK output. If passed by const reference, the wrapper will make a copy so that LAPACK does not overwrite the supplied matrix.
• pivots (optional) : a std::vector<int>, passed by pointer, which will contain the permutation information necessary to reassemble the full matrix information if LAPACK is permitted to modify the input matrix in-place.
• number_of_rhs : The number of columns of \(x\) and \(b\). This is thought of as the 'number of right-hand-sides' to be solved for. In most cases, that additional dimension can be inferred from the dimensionality of the matrix and the input vector(s), which occurs if number_of_rhs is set to 0 (default). If the provided DataVector is an inappropriate size for that inference (i.e. is not a multiple of the largest dimension of the matrix), or if you only want to solve for the first number_of_rhs columns, the parameter number_of_rhs must be supplied.

The function return int is the value provided by the INFO field of the LAPACK call. It is 0 for a successful linear solve, and nonzero values code for types of failures of the algorithm. See LAPACK documentation for further details about dgesv: http://www.netlib.org/lapack/

general_matrix_linear_solve() [5/6]

int lapack::general_matrix_linear_solve	(	gsl::not_null< DataVector * >	solution,
		gsl::not_null< Matrix * >	matrix_operator,
		const DataVector &	rhs,
		int	number_of_rhs = 0 )

Wrapper for LAPACK dgesv, which solves the general linear equation \(A x = b\) by LUP (Lower-triangular, upper-triangular, and permutation) decomposition.

Details

Several interfaces are provided with combinations of input parameters due to the versatility of the LAPACK utility.

• rhs_in_solution_out or solution : If the rhs DataVector is not separately supplied, this single pass-by-pointer DataVector acts as the input \( b\), and is overwritten with the solution \(x\); in this case rhs_in_solution_out must be sufficiently large to contain the solution. If rhs is separately supplied, solution will contain \(x\) at the end of the algorithm, and must be sufficiently large to contain the solution.
• rhs (optional) : An input DataVector representing \(b\) that is not overwritten.
• matrix_operator: The Matrix \(A\). If passed by pointer, this matrix will be 'destroyed' and overwritten with LU decomposition information. Optionally, a pivots vector may be passed by pointer so that the full LUP decomposition information can be recovered from the LAPACK output. If passed by const reference, the wrapper will make a copy so that LAPACK does not overwrite the supplied matrix.
• pivots (optional) : a std::vector<int>, passed by pointer, which will contain the permutation information necessary to reassemble the full matrix information if LAPACK is permitted to modify the input matrix in-place.
• number_of_rhs : The number of columns of \(x\) and \(b\). This is thought of as the 'number of right-hand-sides' to be solved for. In most cases, that additional dimension can be inferred from the dimensionality of the matrix and the input vector(s), which occurs if number_of_rhs is set to 0 (default). If the provided DataVector is an inappropriate size for that inference (i.e. is not a multiple of the largest dimension of the matrix), or if you only want to solve for the first number_of_rhs columns, the parameter number_of_rhs must be supplied.

The function return int is the value provided by the INFO field of the LAPACK call. It is 0 for a successful linear solve, and nonzero values code for types of failures of the algorithm. See LAPACK documentation for further details about dgesv: http://www.netlib.org/lapack/

general_matrix_linear_solve() [6/6]

int lapack::general_matrix_linear_solve	(	gsl::not_null< DataVector * >	solution,
		gsl::not_null< std::vector< int > * >	pivots,
		gsl::not_null< Matrix * >	matrix_operator,
		const DataVector &	rhs,
		int	number_of_rhs = 0 )

Wrapper for LAPACK dgesv, which solves the general linear equation \(A x = b\) by LUP (Lower-triangular, upper-triangular, and permutation) decomposition.

Details

Several interfaces are provided with combinations of input parameters due to the versatility of the LAPACK utility.

• rhs_in_solution_out or solution : If the rhs DataVector is not separately supplied, this single pass-by-pointer DataVector acts as the input \( b\), and is overwritten with the solution \(x\); in this case rhs_in_solution_out must be sufficiently large to contain the solution. If rhs is separately supplied, solution will contain \(x\) at the end of the algorithm, and must be sufficiently large to contain the solution.
• rhs (optional) : An input DataVector representing \(b\) that is not overwritten.
• matrix_operator: The Matrix \(A\). If passed by pointer, this matrix will be 'destroyed' and overwritten with LU decomposition information. Optionally, a pivots vector may be passed by pointer so that the full LUP decomposition information can be recovered from the LAPACK output. If passed by const reference, the wrapper will make a copy so that LAPACK does not overwrite the supplied matrix.
• pivots (optional) : a std::vector<int>, passed by pointer, which will contain the permutation information necessary to reassemble the full matrix information if LAPACK is permitted to modify the input matrix in-place.
• number_of_rhs : The number of columns of \(x\) and \(b\). This is thought of as the 'number of right-hand-sides' to be solved for. In most cases, that additional dimension can be inferred from the dimensionality of the matrix and the input vector(s), which occurs if number_of_rhs is set to 0 (default). If the provided DataVector is an inappropriate size for that inference (i.e. is not a multiple of the largest dimension of the matrix), or if you only want to solve for the first number_of_rhs columns, the parameter number_of_rhs must be supplied.

The function return int is the value provided by the INFO field of the LAPACK call. It is 0 for a successful linear solve, and nonzero values code for types of failures of the algorithm. See LAPACK documentation for further details about dgesv: http://www.netlib.org/lapack/

inner_product()

template<typename Lhs, typename Rhs>

auto LinearSolver::inner_product	(	const Lhs &	lhs,
		const Rhs &	rhs )

The local part of the Euclidean inner product on the vector space w.r.t. which the addition and scalar multiplication of both Lhs and Rhs is defined.

Details

The linear solver works under the following assumptions:

• The data represented by lhs and rhs can each be interpreted as the local chunk of a vector of the same vector space \(V\). Local means there are vectors \(q, p\in V\) such that lhs and rhs represent the components of these vectors w.r.t. a subset \(B_i\) of a basis \(B\subset V\).
• The * and + operators of Lhs and Rhs implement the scalar multiplication and addition in the vector space locally, i.e. restricted to \(B_i\) in the above sense.
• The inner product is the local part \(\langle p,q\rangle|_{B_i}\) of the standard Euclidean dot product in the vector space so that globally it is \(\langle p,q\rangle=\sum_{i}\langle p,q\rangle|_{B_i}\) for \(B=\mathop{\dot{\bigcup}}_i B_i\).

In practice this means that the full vectors \(p\) and \(q\) can be distributed on many elements, where each only holds local chunks lhs and rhs of the components. Scalar multiplication and addition can be performed locally as expected, but computing the full inner product requires a global reduction over all elements that sums their local inner_products.

magnitude_square()

template<typename T>

double LinearSolver::magnitude_square

(

const T &

vector

)

The local part of the Euclidean inner product of a vector with itself.

See also

LinearSolver::inner_product for details on the inner product. This function just returns the inner product of a vector with itself, and ensures that the result is real.