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)