.\"
.de Id
..
.de Sp
.if n .sp
.if t .sp 0.4
..
.TH solver 2rheolef "rheolef-6.7" "rheolef-6.7" "rheolef-6.7"
.\" label: /*Class:solver
.SH NAME
\fBsolver\fP - direct or interative solver interface
.\" skip: @clindex solver
.\" skip: @clindex csr
.\" skip: @clindex vec
.SH DESCRIPTION

The class implements a matrix factorization:
LU factorization for an unsymmetric matrix and
Choleski fatorisation for a symmetric one.
.PP
Let \fIa\fP be a square invertible matrix in
\fBcsr\fP format (see csr(2)).
.\" begin_example
.Sp
.nf
        csr<Float> a;
.Sp
.fi
.\" end_example
We get the factorization by:
.\" begin_example
.Sp
.nf
        solver sa (a);
.Sp
.fi
.\" end_example
Each call to the direct solver for a*x = b writes either:
.\" begin_example
.Sp
.nf
        vec<Float> x = sa.solve(b);
.Sp
.fi
.\" end_example
When the matrix is modified in a computation loop but
conserves its sparsity pattern, an efficient re-factorization
writes:
.\" begin_example
.Sp
.nf
        sa.update_values (new_a);
        x = sa.solve(b);
.Sp
.fi
.\" end_example
.PP
.SH AUTOMATIC CHOICE AND CUSTOMIZATION

.\" skip: @cindex direct solver
.\" skip: @cindex iterative solver
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.
.PP
These default choices can be supersetted by using explicit options:
.\" begin_example
.Sp
.nf
        solver_option_type opt;
        opt.iterative = true;
        solver sa (a, opt);
.Sp
.fi
.\" end_example
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.
.PP
.SH COMPUTATION OF THE DETERMINANT

.\" skip: @cindex determinant
When using a direct method, the determinant of the \fBa\fP matrix
can be computed as:
.\" begin_example
.Sp
.nf
        solver_option_type opt;
        opt.iterative = false;
        solver sa (a, opt);
        cout << sa.det().mantissa << "*2^" << sa.det().exponant << endl;
.Sp
.fi
.\" end_example
The \fBsa.det()\fP method returns an object of type \fBsolver::determinant_type\fP
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.
.PP
.SH PRECONDITIONNERS FOR ITERATIVE SOLVER

.\" skip: @cindex preconditionner
When using an iterative method, the default is to do no preconditionning.
A suitable preconditionner can be supplied via:
.\" begin_example
.Sp
.nf
        solver_option_type opt;
        opt.iterative = true;
        solver sa (a, opt);
        sa.set_preconditionner (pa);
        x = sa.solve(b);
.Sp
.fi
.\" end_example
The supplied \fBpa\fP variable is also of type \fBsolver\fP.
A library of commonly used preconditionners is still in development.
.PP
.\" skip start:AUTHOR:
.\" skip start:DATE:
.\" END
.SH IMPLEMENTATION
.\" begin_example
.Sp
.nf
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;
.Sp
.fi
.\" end_example

.\" LENGTH = 2
.SH SEE ALSO
csr(2), csr(2)