`solver`

direct and iterative solver

# Description

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);

# Direct versus iterative

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.

# Computing the determinant

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.

# Automatic choice and customization

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.

# Preconditionners for iterative solvers

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.

# Implementation

This documentation has been generated from file linalg/lib/solver.h

The `solver` class is simply an alias to the `solver_basic` class

typedef solver_basic<Float> solver;
rheolef::solver
solver_basic< Float > solver
Definition: solver.h:292

The `solver_basic` class provides an interface to a data container:

template <class T, class M = rheo_default_memory_model>
class solver_basic: public smart_pointer_clone<solver_abstract_rep<T,M> > {
public:
// typedefs:
typedef solver_abstract_rep<T,M> rep;
typedef smart_pointer_clone<rep> base;
typedef typename rep::size_type size_type;
typedef typename rep::determinant_type determinant_type;
// allocators:
solver_basic ();
explicit solver_basic (const csr<T,M>& a, const solver_option& opt = solver_option());
void update_values (const csr<T,M>& a);
// accessors:
vec<T,M> trans_solve (const vec<T,M>& b) const;
vec<T,M> solve (const vec<T,M>& b) const;
determinant_type det() const;
const solver_option& option() const;
void set_preconditioner (const solver_basic<T,M>&);
bool initialized() const;
std::string name() const;
};
rheolef::solver_abstract_rep::size_type
csr< T, M >::size_type size_type
Definition: solver.h:193
rheolef::solver_basic::determinant_type
rep::determinant_type determinant_type
Definition: solver.h:271
rheolef::solver_basic::set_preconditioner
void set_preconditioner(const solver_basic< T, M > &)
Definition: solver.h:331
rheolef::solver_basic::name
std::string name() const
Definition: solver.h:366
rheolef::solver_basic::det
determinant_type det() const
Definition: solver.h:338
rheolef::solver_basic::update_values
void update_values(const csr< T, M > &a)
Definition: solver.h:313
rheolef::solver_basic::solve
vec< T, M > solve(const vec< T, M > &b) const
Definition: solver.h:345
rheolef::solver_basic::trans_solve
vec< T, M > trans_solve(const vec< T, M > &b) const
Definition: solver.h:352
rheolef::solver_basic::option
const solver_option & option() const
Definition: solver.h:324
rheolef::solver_basic::rep
solver_abstract_rep< T, M > rep
Definition: solver.h:268
rheolef::solver_basic::size_type
rep::size_type size_type
Definition: solver.h:270
rheolef::solver_basic::base
smart_pointer_clone< rep > base
Definition: solver.h:269
rheolef::solver_basic::solver_basic
solver_basic()
Definition: solver.h:300
rheolef::solver_basic::initialized
bool initialized() const
Definition: solver.h:359
mkgeo_ball.b
int b
Definition: mkgeo_ball.sh:152
a
Definition: diffusion_isotropic.h:25