Rheolef
7.2
an efficient C++ finite element environment
|
The class implements the numerical resolution of a linear system. Let a
be a square and invertible matrix in the csr
sparse format. The construction of a solver writes:
solver sa (a);
and the resolution of a*x = b
expresses simply:
vec<Float> x = sa.solve(b);
When the matrix is modified in a computation loop, the solver could be re-initialized by:
sa.update_values (new_a); vec<Float> x = sa.solve(b);
The choice between a direct or an iterative method for solving the linear system is by default performed automatically: it depends upon the sparsity pattern of the matrix, in order to achieve the best performances. The solver_option
class allows one to change this default behavior.
solver_option sopt; sopt.iterative = true; solver sa (a, sopt); vec<Float> x = sa.solve(b);
The direct approach bases on the Choleski factorization for a symmetric definite positive matrix, and on the LU one otherwise. Conversely, the iterative approach bases on the cg
conjugate gradient algorithm for a symmetric definite positive matrix, and on the gmres
algorithm otherwise.
This feature is useful e.g. when tracking a change of sign in the determinant of a matrix, e.g. during the continuation
algorithm. When using a direct method, the determinant of the matrix can be computed as:
solver_option sopt; sopt.iterative = false; solver sa (a, sopt); cout << sa.det().mantissa << "*" << sa.det().base << "^" << sa.det().exponant << endl;
The sa.det()
method returns an object of type solver::determinant_type
that contains a mantissa and an exponent in a given base (generally 2 or 10). For some rare direct solvers, the computation of the determinant is not yet fully supported: it is the case for the Cholesky factorization from the eigen
library. In you find such a problem, please switch to another solver library, see the solver_option
class.
When the matrix is obtained from the finite element discretization of 3D partial differential problems, the iterative solver is the default choice. Otherwise, the direct solver is selected.
More precisely, the choice between direct or iterative solver depends upon the a.pattern_dimension()
property of the csr
sparse matrix. When this pattern dimension is 3, an iterative method is faster and less memory consuming than a direct one. See User's guide for a discussion on this subject.
The symmetry-positive-definiteness of the matrix is tested via the a.is_symmetric()
and a.is_definite_positive()
properties of the csr
sparse matrix. These properties determine the choices between Cholesky/LU methods for the direct case, and between cg/gmres
algorithms for the iterative one. Most of the time, these properties are automatically well initialized by the finite element assembly procedure, via the integrate
function.
Nevertheless, in some special cases, e.g. a linear combination of matrix, or when the matrix has been read from a file, it could be necessary to force either the symmetry or the positive-definiteness property by the appropriate csr
member function before to send the matrix to the solver.
When using an iterative method, the default is to perform no preconditionning. Several preconditionners are available: the mic
modified incomplete Cholesky for symmetric matrix and the ilut
incomplete LU one for unsymmetric matrix and the do-nothing eye
identity preconditionner. A preconditionner can be supplied via:
solver_option sopt; sopt.iterative = true; solver sa (a, sopt); sa.set_preconditionner (ilut(a)); vec<Float> x = sa.solve(b);
Note also the eye
that returns the solver for the identity matrix: it could be be used for specifying that we do not use a preconditionner. This is the default behavior. The set_preconditionner
member function should be called before the first call to the solve
method: if no preconditioner has been defined at the first call to solve
, the default eye
preconditionner is selected.
This documentation has been generated from file linalg/lib/solver.h
The solver
class is simply an alias to the solver_basic
class
The solver_basic
class provides an interface to a data container: