// ============================================================================
// ==                                                                        ==
// == AMDiS - Adaptive multidimensional simulations                          ==
// ==                                                                        ==
// ============================================================================
// ==                                                                        ==
// ==  TU Dresden                                                            ==
// ==                                                                        ==
// ==  Institut für Wissenschaftliches Rechnen                               ==
// ==  Zellescher Weg 12-14                                                  ==
// ==  01069 Dresden                                                         ==
// ==  germany                                                               ==
// ==                                                                        ==
// ============================================================================
// ==                                                                        ==
// ==  https://gforge.zih.tu-dresden.de/projects/amdis/                      ==
// ==                                                                        ==
// ============================================================================

/** \file ITL_Preconditioner.h */

#ifndef AMDIS_ITL_PRECONDITIONER_H
#define AMDIS_ITL_PRECONDITIONER_H

#include "DOFMatrix.h"
#include "CreatorInterface.h"

#include <boost/numeric/itl/itl.hpp>
#include <boost/numeric/itl/pc/ilu_0.hpp>
#include <boost/numeric/itl/pc/ic_0.hpp>

namespace AMDiS {

  // ============================================================================
  // ===== class ITL_Preconditioner ==============================================
  // ============================================================================

 
  /**
   * \ingroup Solver
   * 
   * \brief Common base class for wrappers to use ITL preconditioners in AMDiS.
   */
  class ITL_BasePreconditioner 
  {
  public:
      typedef DOFMatrix::value_type value_type;

      virtual ~ITL_BasePreconditioner() {}

      virtual mtl::dense_vector<value_type>
      member_solve(const mtl::dense_vector<value_type>& vin) const = 0;

      virtual mtl::dense_vector<value_type>
      member_adjoint_solve(const mtl::dense_vector<value_type>& vin) const = 0;
  };

    inline mtl::dense_vector<DOFMatrix::value_type>
    solve(const ITL_BasePreconditioner& P, const mtl::dense_vector<DOFMatrix::value_type>& vin)
    {
	return P.member_solve(vin);
    }

    inline mtl::dense_vector<DOFMatrix::value_type>
    adjoint_solve(const ITL_BasePreconditioner& P, const mtl::dense_vector<DOFMatrix::value_type>& vin)
    {
	return P.member_adjoint_solve(vin);
    }

  /**
   * \ingroup Solver
   * 
   * \brief Wrapper for using ITL preconditioners in AMDiS.
   */
  template <typename Preconditioner>
  class ITL_Preconditioner : public ITL_BasePreconditioner
  {
    typedef DOFMatrix::value_type value_type;
  public:
    ITL_Preconditioner(const DOFMatrix::base_matrix_type& A) : precon(A) {}
   
    mtl::dense_vector<value_type>
    member_solve(const mtl::dense_vector<value_type>& vin) const
    {
      return solve(precon, vin);
    }

    mtl::dense_vector<value_type>
    member_adjoint_solve(const mtl::dense_vector<value_type>& vin) const
    {
      return adjoint_solve(precon, vin);
    }

    /// Creator class used in the OEMSolverMap.
    class Creator : public CreatorInterface<ITL_BasePreconditioner>
    {
    public:
      virtual ~Creator() {}
      
      /** \brief
       * Creates an ITL preconditioner
       */
      ITL_BasePreconditioner *create(const DOFMatrix::base_matrix_type& A) { 
	return new ITL_Preconditioner<Preconditioner>(A);
      }

      ITL_BasePreconditioner *create() {
	ERROR_EXIT("Must not be called! Only defined to avoid abstract function."); return 0;
      }
    };

  private:
      Preconditioner precon;
  };


  // ============================================================================
  // ===== class DiagonalPreconditioner =========================================
  // ============================================================================

  /**
   * \ingroup Solver
   * 
   * \brief
   * Diagonal preconditioner. 
   */
  class DiagonalPreconditioner 
      : public ITL_Preconditioner<itl::pc::diagonal<DOFMatrix::base_matrix_type> >
  {
      typedef ITL_Preconditioner<itl::pc::diagonal<DOFMatrix::base_matrix_type> > base;

      DiagonalPreconditioner(const DOFMatrix::base_matrix_type& A) : base(A) {}
  };

  /**
   * \ingroup Solver
   * 
   * \brief
   * Identity preconditioner. Behaves like no preconditioning.
   */
  class IdentityPreconditioner 
      : public ITL_Preconditioner<itl::pc::identity<DOFMatrix::base_matrix_type> >
  {
      typedef ITL_Preconditioner<itl::pc::identity<DOFMatrix::base_matrix_type> > base;
  public:
      IdentityPreconditioner(const DOFMatrix::base_matrix_type& A) : base(A) {}
  };


  // ============================================================================
  // ===== class ILUPreconditioner ==============================================
  // ============================================================================

  /**
   * \ingroup Solver
   * 
   * \brief
   * ILU (Incomplete LU factorization) preconditioner.
   *
   * The preconditioner is used from ITL. It corresponds for instance to "Iterative Methods for
   * Sparce Linear Systems", second edition, Yousef Saad. The preconditioner is
   * described in chapter 10.3 (algorithm 10.4).
   */
  class ILUPreconditioner 
      : public ITL_Preconditioner< itl::pc::ilu_0<DOFMatrix::base_matrix_type> >
  {};

  class ICPreconditioner 
      : public ITL_Preconditioner< itl::pc::ic_0<DOFMatrix::base_matrix_type> >
  {};


} // namespace AMDiS

#endif // AMDIS_ITL_PRECONDITIONER_H