NAME
solver - direct or interative solver interface
DESCRIPTION
The class implements a matrix factorization: LU factorization for an unsymmetric matrix and Choleski fatorisation for a symmetric one.
Let a be a square invertible matrix in csr format (see csr(2)).
csr<Float> a;
We get the factorization by:
solver sa (a);
Each call to the direct solver for a*x = b writes either:
vec<Float> x = sa.solve(b);
When the matrix is modified in a computation loop but conserves its sparsity pattern, an efficient re-factorization writes:
sa.update_values
(new_a);
x = sa.solve(b);
AUTOMATIC CHOICE AND CUSTOMIZATION
The symmetry of the matrix is tested via the a.is_symmetric() property (see csr(2)) while the choice between direct or iterative solver is switched by default from the a.pattern_dimension() value. When the pattern is 3D, an iterative method is faster and less memory consuming. Otherwise, for 1D or 2D problems, the direct method is preferred.
These default choices can be supersetted by using explicit options:
solver_option
opt;
opt.iterative = true;
solver sa (a, opt);
When using an iterative method, the sa.solve(b) call the conjugate gradient when the matrix is symmetric, or the generalized minimum residual algorithm when the matrix is unsymmetric.
COMPUTATION OF THE DETERMINANT
When using a direct method, the determinant of the a matrix can be computed as:
solver_option
opt;
opt.iterative = false;
solver sa (a, opt);
cout << sa.det().mantissa << "*2^"
<< sa.det().exponant << endl;
The sa.det() method returns an object of type solver::determinant_type that contains a mantissa and an exponent in base 2. This feature is useful e.g. when tracking a change of sign in the determinant of a matrix.
PRECONDITIONNERS FOR ITERATIVE SOLVER
When using an iterative method, the default is to do no preconditionning. A suitable preconditionner can be supplied via:
solver_option
opt;
opt.iterative = true;
solver sa (a, opt);
sa.set_preconditionner (pa);
x = sa.solve(b);
The supplied pa variable is also of type solver. A library of commonly used preconditionners is still in development.
IMPLEMENTATION
template
<class T, class M = rheo_default_memory_model>
class solver_basic : public
smart_pointer<solver_rep<T,M> > {
public:
// typedefs:
typedef
solver_rep<T,M> rep;
typedef smart_pointer<rep> base;
typedef typename rep::size_type size_type;
typedef typename rep::determinant_type determinant_type;
// allocator:
solver_basic
();
explicit solver_basic (const csr<T,M>& a, const
solver_option& opt = solver_option());
const solver_option& option() const;
void update_values (const csr<T,M>& a);
void set_preconditionner (const
solver_basic<T,M>&);
// 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;
};
// factorizations:
template <class T, class M>
solver_basic<T,M> ldlt(const csr<T,M>& a,
const solver_option& opt = solver_option());
template <class T, class M>
solver_basic<T,M> lu (const csr<T,M>& a,
const solver_option& opt = solver_option());
//
preconditionners:
template <class T, class M =
rheo_default_memory_model>
solver_basic<T,M> eye_basic();
inline solver_basic<Float> eye() { return
eye_basic<Float>(); }
template <class T, class M>
solver_basic<T,M> ic0 (const csr<T,M>& a,
const solver_option& opt = solver_option());
template <class T, class M>
solver_basic<T,M> ilu0(const csr<T,M>& a,
const solver_option& opt = solver_option());
template <class T, class M>
solver_basic<T,M> ldlt_seq(const csr<T,M>&
a, const solver_option& opt = solver_option());
template <class T, class M>
solver_basic<T,M> lu_seq (const csr<T,M>& a,
const solver_option& opt = solver_option());
typedef solver_basic<Float> solver;
SEE ALSO
COPYRIGHT
Copyright (C) 2000-2018 Pierre Saramito <Pierre.Saramito [AT] imag.fr> GPLv3+: GNU GPL version 3 or later <http://gnu.org/licenses/gpl.html>;. This is free software: you are free to change and redistribute it. There is NO WARRANTY, to the extent permitted by law.