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

       csr(2), csr(2)

COPYRIGHT

       Copyright (C) 2000-2018 Pierre Saramito <Pierre.Saramito@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.