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. Otherwhise, for 1D or 2D problems, the direct method is prefered.

These default choices can be supersetted by using explicit options:

        solver_option_type 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_type 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 usefull 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_type 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_type& opt = solver_option_type());
  const solver_option_type& 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_type& opt = solver_option_type());
template <class T, class M>
solver_basic<T,M> lu  (const csr<T,M>& a, const solver_option_type& opt = solver_option_type());
// 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_type& opt = solver_option_type());
template <class T, class M>
solver_basic<T,M> ilu0(const csr<T,M>& a, const solver_option_type& opt = solver_option_type());
template <class T, class M>
solver_basic<T,M> ldlt_seq(const csr<T,M>& a, const solver_option_type& opt = solver_option_type());
template <class T, class M>
solver_basic<T,M> lu_seq  (const csr<T,M>& a, const solver_option_type& opt = solver_option_type());
typedef solver_basic<Float> solver;