librheolef-dev 6.5-1build1 / usr / include / rheolef / solver.h

#ifndef _RHEOLEF_SOLVER_H
#define _RHEOLEF_SOLVER_H
///
/// This file is part of Rheolef.
///
/// Copyright (C) 2000-2009 Pierre Saramito <Pierre.Saramito@imag.fr>
///
/// Rheolef is free software; you can redistribute it and/or modify
/// it under the terms of the GNU General Public License as published by
/// the Free Software Foundation; either version 2 of the License, or
/// (at your option) any later version.
///
/// Rheolef is distributed in the hope that it will be useful,
/// but WITHOUT ANY WARRANTY; without even the implied warranty of
/// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
/// GNU General Public License for more details.
///
/// You should have received a copy of the GNU General Public License
/// along with Rheolef; if not, write to the Free Software
/// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
///
/// =========================================================================
// direct solver interface
//
#include "rheolef/csr.h"

namespace rheolef {

// =======================================================================
// options for the direct solver
// =======================================================================

class solver_option_type {
public:
  typedef std::size_t size_type;

  static const long int  decide = -1;         // Indicate the solver will try to choose the best method

// data:

  mutable long int   iterative;          // Indicate if we want to use iterative solver
  size_type          level_of_fill;      // Level of fill [1:5] for incomplete factorisation
  size_type          amalgamation;       // Level of amalgamation [10:70] for Kass                         
  size_type          ooc;                // out-of-core limit (Mo/percent depending on compilation options) 
  double	     tol;		 // tolerance when using iterative methode
  size_type	     max_iter;		 // maximim number of iteration when using iterative method
  size_type          verbose_level;      // 0, 1, 2, 3
  bool               do_check;

// allocator with default values:

  solver_option_type ()
   : iterative(decide),
     level_of_fill(1),
     amalgamation(10),
     ooc(20000),
     tol(1e-10),
     max_iter(1000),
     verbose_level(0),
     do_check(false)
  {}
};
// =======================================================================
// solver_abstract_rep: an abstract interface for {pastix,trilinos} libs
// =======================================================================
template <class T, class M>
class solver_abstract_rep {
public:
  solver_abstract_rep (const solver_option_type& opt) : _opt(opt) {}
  virtual ~solver_abstract_rep () {}
  virtual void update_values (const csr<T,M>& a) = 0;
  virtual vec<T,M> trans_solve (const vec<T,M>& b) const = 0;
  virtual vec<T,M> solve       (const vec<T,M>& b) const = 0;
  solver_option_type& option() const { return _opt; }
// data:
  mutable solver_option_type _opt;
};
// =======================================================================
// solver_rep: a pointer to the abstract representation
// =======================================================================
template <class T, class M>
class solver_rep {
public:
  explicit solver_rep   ();
  explicit solver_rep   (const csr<T,M>& a, const solver_option_type& opt = solver_option_type());
  void build_lu  (const csr<T,M>& a, const solver_option_type& opt);
  void build_ilu (const csr<T,M>& a, const solver_option_type& opt);
  ~solver_rep ();
  void update_values (const csr<T,M>& a);
  vec<T,M> trans_solve (const vec<T,M>& b) const;
  vec<T,M> solve       (const vec<T,M>& b) const;
protected:
  typedef solver_abstract_rep<T,M>  rep;
  rep* _ptr;
}; 
// =======================================================================
// the direct solver interface
// =======================================================================
/*Class:solver
NAME:   @code{solver} - direct or interative solver interface
@clindex solver
@clindex csr
@clindex vec
@clindex permutation
@toindex pastix, multifrontal solver library
@cindex Choleski factorization
@cindex direct solver
@cindex supernodal factorization
DESCRIPTION:
 @noindent
 The class implements a matrix factorization:
 LU factorization for an unsymmetric matrix and
 Choleski fatorisation for a symmetric one.

 Let @var{a} be a square invertible matrix in 
 @code{csr} format (@pxref{csr class}).
 @example
     	csr<Float> a;
 @end example
 @noindent
 We get the factorization by:
 @example
     	solver<Float> sa (a);
 @end example
 @noindent
 Each call to the direct solver for a*x = b writes either:
 @example
     	vec<Float> x = sa.solve(b);
 @end example
 When the matrix is modified in a computation loop but
 conserves its sparsity pattern, an efficient re-factorization
 writes:
 @example
     	sa.update_values (new_a);
     	x = sa.solve(b);
 @end example
 This approach skip the long step of the symbolic factization step.

ITERATIVE SOLVER:
 The factorization can also be incomplete, i.e. a pseudo-inverse,
 suitable for preconditionning iterative methods.
 In that case, the sa.solve(b) call runs a conjugate gradient
 when the matrix is symmetric, or a generalized minimum residual
 algorithm when the matrix is unsymmetric.

AUTOMATIC CHOICE AND CUSTOMIZATION:
 The symmetry of the matrix is tested via the a.is_symmetric() property
 (@pxref{csr class}) while the choice between direct or iterative solver
 is switched 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:
 @example
	solver_option_type opt;
	opt.iterative = true;
     	solver<Float> sa (a, opt);
 @end example
 See the solver.h header for the complete list of available options.

IMPLEMENTATION NOTE:       
 The implementation bases on the pastix library.

AUTHOR: Pierre.Saramito@imag.fr
DATE:   4 march 2011
End:
*/
//<verbatim:
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;

// allocator:

  solver_basic ();
  explicit solver_basic (const csr<T,M>& a, const solver_option_type& opt = solver_option_type());
  void update_values (const csr<T,M>& a);

// accessors:

  vec<T,M> trans_solve (const vec<T,M>& b) const;
  vec<T,M> solve       (const vec<T,M>& b) 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());
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());

typedef solver_basic<Float> solver;
//>verbatim:

// =======================================================================
// solver_basic: inlined
// =======================================================================
template <class T, class M>
inline
solver_basic<T,M>::solver_basic ()
 : base (new_macro(rep))
{
}
template <class T, class M>
inline
solver_basic<T,M>::solver_basic (const csr<T,M>& a, const solver_option_type& opt)
 : base (new_macro(rep(a,opt)))
{
}
template <class T, class M>
inline
solver_basic<T,M>
lu (const csr<T,M>& a, const solver_option_type& opt)
{
  solver_basic<T,M> sa;
  sa.data().build_lu (a, opt);
  return sa;
}
template <class T, class M>
inline
solver_basic<T,M>
ldlt (const csr<T,M>& a, const solver_option_type& opt)
{
  return lu(a,opt); // symmetry flag is indicated in the matrix a
}
template <class T, class M>
inline
solver_basic<T,M>
ilu0 (const csr<T,M>& a, const solver_option_type& opt)
{
  solver_basic<T,M> sa;
  sa.data().build_ilu (a, opt);
  return sa;
}
template <class T, class M>
inline
solver_basic<T,M>
ic0 (const csr<T,M>& a, const solver_option_type& opt)
{
  return ilu0(a,opt); // symmetry flag is indicated in the matrix a
}
template <class T, class M>
inline
void
solver_basic<T,M>::update_values (const csr<T,M>& a)
{
  base::data().update_values (a);
}
template <class T, class M>
inline
vec<T,M>
solver_basic<T,M>::solve (const vec<T,M>& b) const
{
  return base::data().solve (b);
}
template <class T, class M>
inline
vec<T,M>
solver_basic<T,M>::trans_solve (const vec<T,M>& b) const
{
  return base::data().trans_solve (b);
}

} // namespace rheolef
#endif // _RHEOLEF_SOLVER_H