#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
//
// =========================================================================
// AUTHOR: Pierre.Saramito@imag.fr
// DATE:   4 march 2011

namespace rheolef {
/**
@linalgclassfile solver direct and iterative solver
@addindex preconditioner
@addindex determinant

Description
===========
The class implements the numerical resolution of a linear system.
Let `a` be a square and invertible matrix in 
the @ref csr_4 sparse format.
The construction of a solver writes:

        solver sa (a);

and the resolution of `a*x = b` expresses simply:

        vec<Float> x = sa.solve(b);

When the matrix is modified in a computation loop,
the solver could be re-initialized by:

        sa.update_values (new_a);
        vec<Float> x = sa.solve(b);

Direct versus iterative
=======================
The choice between a direct or an iterative method for 
solving the linear system 
is by default performed automatically:
it depends upon the sparsity pattern of the matrix,
in order to achieve the best performances.
The @ref solver_option_4 class allows one to change this default behavior.

        solver_option sopt;
        sopt.iterative = true;
        solver sa (a, sopt);
        vec<Float> x = sa.solve(b);

The direct approach bases on the Choleski factorization
for a symmetric definite positive matrix,
and on the LU one otherwise.
Conversely, the iterative approach bases on the
@ref cg_5 conjugate gradient algorithm
for a symmetric definite positive matrix,
and on the @ref gmres_5 algorithm otherwise.

Computing the determinant
=========================
This feature is useful e.g. when tracking a change of sign
in the determinant of a matrix, e.g. during 
the @ref continuation_3 algorithm.
When using a direct method, the determinant of the matrix
can be computed as:

        solver_option sopt;
        sopt.iterative = false;
        solver sa (a, sopt);
        cout << sa.det().mantissa << "*" 
             << sa.det().base << "^"
             << sa.det().exponant << endl;

The `sa.det()` method returns an object of type
`solver::determinant_type`
that contains a mantissa and an exponent in a given base (generally 2 or 10).
For some rare direct solvers, the computation of the determinant
is not yet fully supported:
it is the case for the Cholesky factorization from the `eigen` library. 
In you find such a problem, please switch to another solver library,
see the @ref solver_option_4 class.

Automatic choice and customization
==================================
When the matrix is obtained from the finite element
discretization of 3D partial differential problems,
the iterative solver is the default choice.
Otherwise, the direct solver is selected.

More precisely, the choice between direct or iterative solver
depends upon the `a.pattern_dimension()` property
of the @ref csr_4 sparse matrix.
When this pattern dimension is 3,
an iterative method is faster and less memory consuming
than a direct one.
See @ref usersguide_page for a discussion on this subject.

The symmetry-positive-definiteness of the matrix is tested via the
`a.is_symmetric()`
and
`a.is_definite_positive()`
properties of the @ref csr_4 sparse matrix.
These properties determine the choices
between Cholesky/LU methods for the direct case,
and between `cg/gmres` algorithms for the iterative one.
Most of the time, these properties are automatically
well initialized by the finite element assembly procedure,
via the @ref integrate_3 function.

Nevertheless, in some special cases, 
e.g. a linear combination of matrix,
or when the matrix has been read from a file,
it could be necessary to force either the symmetry
or the positive-definiteness property
by the appropriate @ref csr_4 member function
before to send the matrix to the solver.

Preconditionners for iterative solvers
======================================
When using an iterative method,
the default is to perform no preconditionning.
Several preconditionners are available:
the @ref mic_5 modified incomplete Cholesky
for symmetric matrix
and the @ref ilut_5 incomplete LU one
for unsymmetric matrix
and the do-nothing @ref eye_5 identity preconditionner.
A preconditionner can be supplied via:

        solver_option sopt;
        sopt.iterative = true;
        solver sa (a, sopt);
        sa.set_preconditionner (ilut(a));
        vec<Float> x = sa.solve(b);

Note also the @ref eye_5 that returns the solver for the identity matrix:
it could be be used for specifying that we do not use a preconditionner. 
This is the default behavior.
The `set_preconditionner` member function should be called
before the first call to the `solve` method:
if no preconditioner has been defined at the first call
to `solve`, the default @ref eye_5 preconditionner is selected.

Implementation
==============
@showfromfile
The `solver` class is simply an alias to the `solver_basic` class

@snippet solver.h verbatim_solver
@par

The `solver_basic` class provides an interface
to a data container:

@snippet solver.h verbatim_solver_basic
*/
} // namespace rheolef

#include "rheolef/csr.h"
#include "rheolef/solver_option.h"

namespace rheolef {

// =======================================================================
// solver_abstract_rep: an abstract interface for solvers
// =======================================================================
// forward declaration:
template <class T, class M> class solver_basic;

template <class T, class M>
class solver_abstract_rep {
public:
  typedef typename csr<T,M>::size_type  size_type;
  struct determinant_type;
  explicit solver_abstract_rep (const solver_option& opt) : _opt(opt) {}
  solver_abstract_rep (const solver_abstract_rep& x) : _opt(x._opt) {}
  const solver_option& option() const { return _opt; }
  virtual  solver_abstract_rep<T,M>* clone() const;
  virtual bool initialized() const { return false; }
  virtual ~solver_abstract_rep () {}
  virtual void update_values (const csr<T,M>& a) {}
  virtual vec<T,M> solve       (const vec<T,M>& b) const;
  virtual vec<T,M> trans_solve (const vec<T,M>& b) const;
  virtual void set_preconditioner (const solver_basic<T,M>&);
  virtual determinant_type det() const;
  virtual std::string name() const;
// internals:
  static solver_abstract_rep<T,M>* make_solver_ptr (const csr<T,M>& a, const solver_option& opt);
// data:
protected:
  mutable solver_option _opt;
};
// det = mantissa*base^exponent
template <class T, class M>
struct solver_abstract_rep<T,M>::determinant_type {
  T       mantissa, exponant, base;
  determinant_type(): mantissa(0), exponant(0), base(10) {} 
};
// -----------------------------------------------------------------------------
// solver_abstract_rep inlined
// -----------------------------------------------------------------------------
template <class T, class M>
inline
solver_abstract_rep<T,M>*
solver_abstract_rep<T,M>::clone() const
{
  typedef solver_abstract_rep<T,M> rep;
  return new_macro (rep(*this));
}
template <class T, class M>
inline
vec<T,M>
solver_abstract_rep<T,M>::solve (const vec<T,M>& b) const
{
  error_macro (name() << ": undefined solve() member function"); return vec<T,M>();
}
template <class T, class M>
inline
vec<T,M>
solver_abstract_rep<T,M>::trans_solve (const vec<T,M>& b) const
{
  error_macro (name() << ": undefined trans_solve() member function"); return vec<T,M>();
}
template <class T, class M>
inline
void
solver_abstract_rep<T,M>::set_preconditioner (const solver_basic<T,M>&)
{
  error_macro (name() << ": undefined set_preconditioner() member function");
}
template <class T, class M>
inline
typename solver_abstract_rep<T,M>::determinant_type
solver_abstract_rep<T,M>::det() const
{
  error_macro (name() << ": undefined det() member function");
  return determinant_type();
}
// =======================================================================
// the direct & iterative solver interface
// =======================================================================
// [verbatim_solver_basic]
template <class T, class M = rheo_default_memory_model>
class solver_basic: public smart_pointer_clone<solver_abstract_rep<T,M> > {
public:
// typedefs:

  typedef solver_abstract_rep<T,M>        rep;
  typedef smart_pointer_clone<rep>        base;
  typedef typename rep::size_type         size_type;
  typedef typename rep::determinant_type  determinant_type;

// allocators:

  solver_basic ();
  explicit solver_basic (const csr<T,M>& a, const solver_option& opt = solver_option());
  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;
  determinant_type det() const;
  const solver_option& option() const;
  void set_preconditioner (const solver_basic<T,M>&);
  bool initialized() const;
  std::string name() const;
};
// [verbatim_solver_basic]

// [verbatim_solver]
typedef solver_basic<Float> solver;
// [verbatim_solver]

// -----------------------------------------------------------------------
// solver_basic: inlined
// -----------------------------------------------------------------------
template <class T, class M>
inline
solver_basic<T,M>::solver_basic ()
 : base (new_macro(rep(solver_option())))
{
}
template <class T, class M>
inline
solver_basic<T,M>::solver_basic (const csr<T,M>& a, const solver_option& opt)
 : base (solver_abstract_rep<T,M>::make_solver_ptr(a,opt))
{
}
template <class T, class M>
inline
void
solver_basic<T,M>::update_values (const csr<T,M>& a)
{
  if (base::data().initialized()) {
    base::data().update_values (a);
  } else {
    solver_basic<T,M>::base::operator= (solver_abstract_rep<T,M>::make_solver_ptr(a,solver_option()));
  }
}
template <class T, class M>
inline
const solver_option&
solver_basic<T,M>::option() const
{
  return base::data().option();
}
template <class T, class M>
inline
void
solver_basic<T,M>::set_preconditioner (const solver_basic<T,M>& sa)
{
  base::data().set_preconditioner (sa);
}
template <class T, class M>
inline
typename solver_basic<T,M>::determinant_type
solver_basic<T,M>::det() const
{
  return base::data().det();
}
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);
}
template <class T, class M>
inline
bool
solver_basic<T,M>::initialized() const
{
  return base::data().initialized();
}
template <class T, class M>
inline
std::string
solver_basic<T,M>::name() const
{
  return base::data().name();
}

} // namespace rheolef
#endif // _RHEOLEF_SOLVER_H