// ============================================================================
// ==                                                                        ==
// == 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_PRECONDITIONER_H
#define AMDIS_ITL_PRECONDITIONER_H

#include "MTL4Types.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>
#include <boost/numeric/mtl/vector/assigner.hpp>

namespace AMDiS {

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

 
  /**
   * \ingroup Solver
   * 
   * \brief Common base class for wrappers to use ITL preconditioners in AMDiS.
   */
  template< class Vec >
  class ITL_BasePreconditioner 
  {
  public:
      virtual ~ITL_BasePreconditioner() {}

      // Old style AMDiS
      // virtual Vec member_solve(const Vec& vin) const =0 ;
      // virtual Vec member_adjoint_solve(const Vec& vin) const =0 ;

     // New style without copying the result vector
     virtual void solve(const Vec& x, Vec& y) const= 0;
     virtual void adjoint_solve(const Vec& x, Vec& y) const= 0;
  };



  template< typename Vec >
  inline //Vec
  itl::pc::solver<ITL_BasePreconditioner< Vec >, Vec, false>
  solve(const ITL_BasePreconditioner< Vec >& P, const Vec& vin)
  {
    return itl::pc::solver<ITL_BasePreconditioner< Vec >, Vec, false>(P, vin);
    //return P.member_solve(vin);
  }

  template< typename Vec >
  inline //Vec
  itl::pc::solver<ITL_BasePreconditioner< Vec >, Vec, true>
  adjoint_solve(const ITL_BasePreconditioner< Vec >& P, const Vec& vin)
  {
    return itl::pc::solver<ITL_BasePreconditioner< Vec >, Vec, true>(P, vin);
    //return P.member_adjoint_solve(vin);
  }

  /**
   * \ingroup Solver
   * 
   * \brief Wrapper for using ITL preconditioners in AMDiS.
   */
  template <typename Preconditioner, typename MTLVector, typename MTLMatrix >
  class ITL_Preconditioner : public ITL_BasePreconditioner< MTLVector >
  {
  public:
    ITL_Preconditioner(const MTLMatrix& A) : precon(A) {}
#if 0   
    MTLVector
    member_solve(const MTLVector& vin) const
    {
      return solve(precon, vin);
    }

    MTLVector
    member_adjoint_solve(const MTLVector& vin) const
    {
      return adjoint_solve(precon, vin);
    }
#else
      
    void solve(const MTLVector& vin, MTLVector& vout) const
    {
      precon.solve(vin, vout);
   }

    void adjoint_solve(const MTLVector& vin, MTLVector& vout) const
    {
      precon.adjoint_solve(vin, vout);
   }
#endif

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

  private:
      Preconditioner precon;
  };


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

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

  public:
    DiagonalPreconditioner(const MTLTypes::MTLMatrix& A) : base(A) {}
  };

  /**
   * \ingroup Solver
   * 
   * \brief
   * Identity preconditioner. Behaves like no preconditioning.
   */
  class IdentityPreconditioner 
    : public ITL_Preconditioner<itl::pc::identity<DOFMatrix::base_matrix_type>, MTLTypes::MTLVector, MTLTypes::MTLMatrix >
  {
    typedef ITL_Preconditioner<itl::pc::identity<DOFMatrix::base_matrix_type>, MTLTypes::MTLVector, MTLTypes::MTLMatrix > base;

  public:
      IdentityPreconditioner(const MTLTypes::MTLMatrix& 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>, MTLTypes::MTLVector, MTLTypes::MTLMatrix >
  {};

  class ICPreconditioner 
      : public ITL_Preconditioner< itl::pc::ic_0<DOFMatrix::base_matrix_type>, MTLTypes::MTLVector, MTLTypes::MTLMatrix >
  {};


} // namespace AMDiS

#endif // AMDIS_ITL_PRECONDITIONER_H