// ============================================================================
// ==                                                                        ==
// == AMDiS - Adaptive multidimensional simulations                          ==
// ==                                                                        ==
// ==  http://www.amdis-fem.org                                              ==
// ==                                                                        ==
// ============================================================================
//
// Software License for AMDiS
//
// Copyright (c) 2010 Dresden University of Technology 
// All rights reserved.
// Authors: Simon Vey, Thomas Witkowski et al.
//
// This file is part of AMDiS
//
// See also license.opensource.txt in the distribution.



/** \file ITL_Preconditioner.h */

#ifndef AMDIS_ITL_SOLVER_H
#define AMDIS_ITL_SOLVER_H

#include "ITL_OEMSolver.h"
#include <boost/numeric/itl/krylov/bicg.hpp>
#include <boost/numeric/itl/krylov/bicgstab_2.hpp>
#include <boost/numeric/itl/krylov/bicgstab_ell.hpp>
#include <boost/numeric/itl/krylov/bicgstab.hpp>
#include <boost/numeric/itl/krylov/cg.hpp>
#include <boost/numeric/itl/krylov/cgs.hpp>
#include <boost/numeric/itl/krylov/gmres.hpp>
#include <boost/numeric/itl/krylov/idr_s.hpp>
#include <boost/numeric/itl/krylov/qmr.hpp>
#include <boost/numeric/itl/krylov/tfqmr.hpp>
#include "itl/minres.hpp"


namespace AMDiS {

  // ============================================================================
  // ===== class CGSolver =======================================================
  // ============================================================================
  
  /// Helper class to establish a type from a function
  class cg_solver_type
  {
  public:
    template < class LinOp, class X, class B, class L, class R, class I >
    int operator()(const LinOp& A, X& x, const B& b, const L& l, const R& r, I& iter)
    { return itl::cg(A, x, b, l, r, iter); }
  };


  /**
   * \ingroup Solver
   * 
   * \brief
   * Solves a linear system by the conjugate gradient method (CG) and can be used for
   * symmetric positive definite system matrices.
   * Right preconditioner is ignored.
   */
  class CGSolver : public ITL_OEMSolver<cg_solver_type>
  {
  public:
    /// The constructor reads required parameters and sets solvers \ref name.
    CGSolver(std::string name) : ITL_OEMSolver<cg_solver_type>(name) {}
  };


  // ============================================================================
  // ===== class CGSSolver =======================================================
  // ============================================================================
  
  /// Helper class to establish a type from a function
  class cgs_solver_type
  {
  public:
    template < class LinOp, class X, class B, class L, class R, class I >
    int operator()(const LinOp& A, X& x, const B& b, const L& l, const R&, I& iter)
    { return itl::cgs(A, x, b, l, iter); }
  };


  /**
   * \ingroup Solver
   * 
   * \brief
   * Solves a linear system by the squared conjugate gradient method (CGS).
   * Right preconditioner is ignored.
   */
  class CGSSolver : public ITL_OEMSolver<cgs_solver_type>
  {
  public:
    /// The constructor reads required parameters and sets solvers \ref name.
    CGSSolver(std::string name) : ITL_OEMSolver<cgs_solver_type>(name) {}
  };


  // ============================================================================
  // ===== class BiCGSolver =====================================================
  // ============================================================================

  /// Helper class to establish a type from a function
  class bicg_solver_type
  {
  public:
    template < class LinOp, class X, class B, class L, class R, class I >
    int operator()(const LinOp& A, X& x, const B& b, const L& l, const R&, I& iter)
    { return itl::bicg(A, x, b, l, iter); }
  };

  /**  
   * \ingroup Solver
   *
   * \brief
   * Solves a linear system by a stabilized BiCG method and can be used for 
   * system matrices.
   */
  class BiCGSolver : public ITL_OEMSolver<bicg_solver_type>
  {
  public:
    /// The constructor reads required parameters and sets solvers \ref name.
    BiCGSolver(std::string name) : ITL_OEMSolver<bicg_solver_type>(name) {}
  };


  // ============================================================================
  // ===== class BiCGStab =======================================================
  // ============================================================================

  /// Helper class to establish a type from a function
  class bicgstab_type
  {
  public:
    template < class LinOp, class X, class B, class L, class R, class I >
    int operator()(const LinOp& A, X& x, const B& b, const L& l, const R&, I& iter)
    { return itl::bicgstab(A, x, b, l, iter); }
  };


  /**  
   * \ingroup Solver
   *
   * \brief
   * Solves a linear system by a stabilized BiCG method and can be used for 
   * system matrices.
   */
  class BiCGStabSolver : public ITL_OEMSolver<bicgstab_type>
  {
  public:
    /// The constructor reads required parameters and sets solvers \ref name.
    BiCGStabSolver(std::string name) : ITL_OEMSolver<bicgstab_type>(name) {}
  };


  // ============================================================================
  // ===== class BiCGStab2 ======================================================
  // ============================================================================

  /// Helper class to establish a type from a function
  class bicgstab2_type
  {
  public:
    template < class LinOp, class X, class B, class L, class R, class I >
    int operator()(const LinOp& A, X& x, const B& b, const L& l, const R&, I& iter)
    { return itl::bicgstab_2(A, x, b, l, iter); }
  };


  /**  
   * \ingroup Solver
   *
   * \brief
   * Solves a linear system by a stabilized BiCG method and can be used for 
   * system matrices.
   */
  class BiCGStab2Solver : public ITL_OEMSolver<bicgstab2_type>
  {
  public:
    /// The constructor reads required parameters and sets solvers \ref name.
    BiCGStab2Solver(std::string name) : ITL_OEMSolver<bicgstab2_type>(name) {}
  };


  // ============================================================================
  // ===== class QMRSolver =======================================================
  // ============================================================================
  
  /// Helper class to establish a type from a function
  class qmr_solver_type
  {
  public:
    template < class LinOp, class X, class B, class L, class R, class I >
    int operator()(const LinOp& A, X& x, const B& b, const L& l, const R& r, I& iter)
    { return itl::qmr(A, x, b, l, r, iter); }
  };


  /**
   * \ingroup Solver
   * 
   * \brief
   * Solves a linear system by the Quasi-Minimal Residual method (QMR).
   */
  class QMRSolver : public ITL_OEMSolver<qmr_solver_type>
  {
  public:
    /// The constructor reads required parameters and sets solvers \ref name.
    QMRSolver(std::string name) : ITL_OEMSolver<qmr_solver_type>(name) {}
  };


  // ============================================================================
  // ===== class TFQMRSolver =======================================================
  // ============================================================================
  
  /// Helper class to establish a type from a function
  class tfqmr_solver_type
  {
  public:
    template < class LinOp, class X, class B, class L, class R, class I >
    int operator()(const LinOp& A, X& x, const B& b, const L& l, const R& r, I& iter)
    { 
      return itl::tfqmr(A, x, b, l, r, iter); 
    }
  };


  /**
   * \ingroup Solver
   * 
   * \brief
   * Solves a linear system by the Transposed-Free Quasi-Minimal Residual method (TFQMR).
   * Does not use preconditioning currently.
   */
  class TFQMRSolver : public ITL_OEMSolver<tfqmr_solver_type>
  {
  public:
    /// The constructor reads required parameters and sets solvers \ref name.
    TFQMRSolver(std::string name) : ITL_OEMSolver<tfqmr_solver_type>(name) {}
  };



  // ============================================================================
  // ===== class BiCGStabEll ======================================================
  // ============================================================================

  /// Helper class to establish a type from a function
  class bicgstab_ell_type
  {
  public:
    template < class LinOp, class X, class B, class L, class R, class I >
    int operator()(const LinOp& A, X& x, const B& b, const L& l, const R& r, I& iter, int ell)
    { return itl::bicgstab_ell(A, x, b, l, r, iter, ell); }
  };

  /**  
   * \ingroup Solver
   *
   * \brief
   * Solves a linear system by a stabilized BiCG(l) method and can be used for 
   * system matrices.
   */
  class BiCGStabEllSolver : public ITL_OEMSolver_para<bicgstab_ell_type>
  {
  public:
    /// The constructor reads required parameters and sets solvers \ref name.
    BiCGStabEllSolver(std::string name) : ITL_OEMSolver_para<bicgstab_ell_type>(name) {}
  };

  // ============================================================================
  // ===== class GMRES ======================================================
  // ============================================================================

  /// Helper class to establish a type from a function
  class gmres_type
  {
  public:
    template < class LinOp, class X, class B, class L, class R, class I >
    int operator()(const LinOp& A, X& x, const B& b, const L& l, const R& r, I& iter, int ell)
    { return itl::gmres(A, x, b, l, r, iter, ell); }
  };

  /**  
   * \ingroup Solver
   *
   * \brief
   * Solves a linear system by the GMRES method.
   * The parameter ell is the maximal number of orthogonalized vectors.
   * The method is not preconditioned 
   */
  class GMResSolver : public ITL_OEMSolver_para<gmres_type>
  {
  public:
    /// The constructor reads required parameters and sets solvers \ref name.
    GMResSolver(std::string name) : ITL_OEMSolver_para<gmres_type>(name) {}
  };


  // ============================================================================
  // ===== class IDRs ======================================================
  // ============================================================================

  /// Helper class to establish a type from a function
  class idr_s_type
  {
  public:
    template < class LinOp, class X, class B, class L, class R, class I >
    int operator()(const LinOp& A, X& x, const B& b, const L& l, const R& r, I& iter, int ell)
    { return itl::idr_s(A, x, b, l, r, iter, ell); }
  };

  /**  
   * \ingroup Solver
   *
   * \brief
   * Solves a linear system by an Induced Dimension Reduction method and can be used for 
   * system matrices.
   * The parameter s can be specified as ell.
   * Peter Sonneveld and Martin B. van Gijzen, IDR(s): a family of simple and fast algorithms for solving large nonsymmetric linear systems. 
   * SIAM J. Sci. Comput. Vol. 31, No. 2, pp. 1035-1062 (2008). (copyright SIAM)
   * 
   */
  class IDRsSolver : public ITL_OEMSolver_para<idr_s_type>
  {
  public:
    /// The constructor reads required parameters and sets solvers \ref name.
    IDRsSolver(std::string name) : ITL_OEMSolver_para<idr_s_type>(name) {}
  };


  // ============================================================================
  // ===== class Minres  ========================================================
  // ============================================================================
  
  /// Helper class to establish a type from a function
  class minres_solver_type
  {
  public:
    template < class LinOp, class X, class B, class L, class R, class I >
    int operator()(const LinOp& A, X& x, const B& b, const L& l, const R& r, I& iter)
    { 
      return itl::minres(A, x, b, l, r, iter); 
    }
  };


  /**
   * \ingroup Solver
   * 
   * \brief
   * Solves a linear system by the Minres method. Can be used for symmetric 
   * indefinite systems.
   */
  class MinResSolver : public ITL_OEMSolver<minres_solver_type>
  {
  public:
    /// The constructor reads required parameters and sets solvers \ref name.
    MinResSolver(std::string name) : ITL_OEMSolver<minres_solver_type>(name) {}
  };



} // namespace AMDiS

#endif // AMDIS_ITL_SOLVER