bionic (2) solver.2rheolef.gz

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

SEE ALSO

       csr(2), csr(2)