/***************************************************************************
 *   Copyright (C) 2005 by Ricardo Ortiz                                   *
 *   rortizro@math.uiowa.edu                                               *
 *                                                                         *
 *   This program 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.                                   *
 *                                                                         *
 *   This program 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 this program; if not, write to the                         *
 *   Free Software Foundation, Inc.,                                       *
 *   59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.             *
 ***************************************************************************/
#ifndef MY_MG_H
#define MY_MG_H

#ifdef DEBUG
#define checkpt() (cout << "Checkpoint: file\"" << __FILE__ << "\", line " \
                   << __LINE__ << endl)
#define DBG(stmt) (stmt)
#else
#define checkpt()
#define DBG(stmt)
#endif // DEBUG

#include <iostream>

using std::cout;
using std::endl;

// Note: level = 0 corresponds to the **finest** grid;
//       increasing level corresponds to coarsening the grid.

// Vector operations assumed: Vector x, y, out;
//   Vector(int_type dim)                                -- constructor
//   Vector()                                            -- constructor
//   x.size()                                            -- size (=dimension) of x
//   x.resize(int_type new_dim)                          -- resize to new_dim
//   copy(const Vector &x, Vector &out)                  -- out <- x
//   add(const Vector &x, const Vector &y, Vector &out)  -- out <- x+y
//   scale(Vector &x, value_type scalar)                 -- x <- scalar*x

namespace MY_MG
{


///////////////////////////////////////////////////////////////////////////
///
/// @struct TypeBinder
/// @brief  This class contain all types used in Multigrid
///
///////////////////////////////////////////////////////////////////////////

  template< typename T, class VectorType, 
            class MatrixType, class PermType >
  class TypeBinder
  {
   public:
    typedef TypeBinder< T, VectorType, MatrixType, PermType > Types;
  
    typedef T            value_type;  // Scalar Type
    typedef unsigned int int_type;    // Integer Type
    typedef MatrixType   Matrix;      // Matrix (dense)
    typedef VectorType   Vector;      // Vector (dense)
    typedef PermType     Perm;        // Permutation or pivot
    typedef typename Vector::size_type size_type;
  
    // Grid class contains vectors for internal use
    struct Grid
    {
      Vector x, b, r;
    
      // Constructor of dimension dim
      Grid(size_type dim)
      {  Vector x(dim), b(dim), r(dim);  }
      // Empty constructor
      Grid() { }
      // resize -- sets new dimension
      void resize(size_type new_dim)
      {
        x.resize(new_dim);
        b.resize(new_dim);
        r.resize(new_dim);
      }
    }; // class Grid
    
  }; //TypeBinder
  
  
  
  ///////////////////////////////////////////////////////////////////////////
  ///
  /// @brief This class implements a multigrid V-cycle scheme
  ///
  ///////////////////////////////////////////////////////////////////////////
  template < typename Types >
  class multigrid
  {
   protected:
    // typedef's from Types for better readability
    typedef typename Types::value_type value_type; // Scalar type
    typedef typename Types::int_type   int_type;   // Integer type
    typedef typename Types::size_type  size_type;  // size type
    typedef typename Types::Vector     Vector;     // Vector type
    typedef typename Types::Grid       Grid;       // Grid
    
   private:
    // Private data items
    Grid     *grid;           // Array of grid structures
    int_type *dims;           // Array of dimension of each grid
    int_type  coarsest_level; // Coarsest grid
    int_type  nu[2];          // Relaxation parameters
    
   public:
    // Constructor: builds a multigrid method with n_levels, 
    // and level k has dimension _dims[k].
    // multigrid(int_type n_levels, const int_type _dims[], 
    //           const int_type _nu[])
    multigrid(int_type n_levels, const int_type _dims[], const int_type _nu[])
    {
      dims = new int_type[n_levels];
      grid = new Grid[n_levels];
      coarsest_level = n_levels - 1;
      for ( int i = 0; i < n_levels; i++ ) {
        dims[i] = _dims[i];
        grid[i].resize( dims[i] );
      }
      nu[0] = _nu[0];    nu[1] = _nu[1];
    };
    
  
    // Constructor: need to use resize() later
    // multigrid(int_type n_levels)
    multigrid(int_type n_levels)
    {
      dims = new int_type[n_levels];
      grid = new Grid[n_levels];
      coarsest_level = n_levels - 1;
      nu[0] = 1;    nu[1] = 1;
    }
    
  
    // resize -- Note: dims must have length >= coarsest_level + 1
    // void resize(const int_type _dims[])
    void resize(const int_type _dims[])
    {
      for ( int_type i = 0; i <= coarsest_level; ++i )
        {
          dims[i] = _dims[i];
          grid[i].resize( dims[i] );
        }
    }
    
  
    // set_nu -- update nu[] array; _nu must have length >= 2
    // void set_nu(const int_type _nu[])
    void set_nu(const int_type _nu[])
    {
      nu[0] = _nu[0];    nu[1] = _nu[1];
    }
    
  
    // Destructor: free up grid and dims
    // ~multigrid()
    ~multigrid()
    {
      delete [] dims;
      delete [] grid;
    };
    
  
    int_type const get_max_level()         const 
    { return coarsest_level; }
    int_type const get_dim(int_type level) const 
    { assert(level >= 0 && level <= coarsest_level );
      return dims[level]; };
    int_type const get_nu(int_type i)      const 
    { assert(i >= 0 && i < 2); return nu[i]; };
  
    // residual -- out <- b-A[level]*x
    // Vector &residual(const Vector &x, const Vector &b, 
    //                      int_type level, Vector &out)
    Vector &residual(const Vector &x, const Vector &b, 
                     int_type level, Vector &out)
    {
      assert(level <= coarsest_level && level >= 0);
      assert(x.size() == dims[level] && b.size() == dims[level]
             && out.size() == dims[level]);
      apply_operator(x,level,out);
      scale(out,value_type(-1));
      add(out,b,out);
      return out;
    }
    
  
    // The "ingredients" for multigrid (apply_operator, pre_relax,
    // post_relax, interpolate, restriction, and solve)
    // These are the main "ingredients" for constructing 
    // the multigrid method.  These should return the output vector, 
    // which is the Vector not marked with "const"
    virtual Vector &pre_relax (Vector &v, const Vector &b, int_type level) = 0;
    virtual Vector &post_relax(Vector &v, const Vector &b, int_type level) = 0;
    virtual Vector &solve     (Vector &v, const Vector &b, int_type level) = 0;
    
    virtual Vector &interpolate(const Vector &v, int_type level, 
                                Vector &out) = 0;
    virtual Vector &restriction(const Vector &v, int_type level, 
                                Vector &out) = 0;
    virtual Vector &apply_operator(const Vector &v, int_type level, 
                                   Vector &out) = 0;
    
  
    // MGSolve -- multigrid solver for solving A[0]*x = rhs
    // Uses levels up to "top_level"
    // Vector &MGSolve( int_type top_level, Vector &x, const Vector& rhs )
    Vector &MGSolve( int_type top_level, Vector &x, const Vector& rhs )
    {
      assert(top_level <= coarsest_level && top_level >= 0);
      Vector interpolant(dims[0]);
    
      copy(rhs, grid[0].b);
      copy(x,   grid[0].x);
    
      for ( int_type j = 0;j < top_level; ++j ) {
        assert(grid[j].x.size() == dims[j]);
        pre_relax( grid[j].x, grid[j].b, j );
        residual( grid[j].x, grid[j].b, j, grid[j].r);
        restriction( grid[j].r, j+1, grid[j+1].b );
        // We also need to set grid[j+1].x to zero as well...
        for ( int_type i = 0; i < grid[j+1].x.size(); ++i )
          grid[j+1].x[i] = 0;
        assert(grid[j].r.size() == dims[j]);
        assert(grid[j+1].b.size() == dims[j+1]);
      }
    
      solve( grid[top_level].x, grid[top_level].b, top_level);
    
      for ( int_type j = top_level; j > 0; --j ) {
        interpolant.resize( dims[j-1] );
        interpolate( grid[j].x, j-1, interpolant );
        add( interpolant , grid[j-1].x );
        post_relax( grid[j-1].x, grid[j-1].b, j-1);
      }
      copy(grid[0].x, x);
      return x;
    }
    
  }; // multigrid
  
  



  ///////////////////////////////////////////////////////////////////////////
  ///
  /// @brief This class implements a Conjugate Gradient algorithm(the solver)
  ///
  ///////////////////////////////////////////////////////////////////////////
  /**********************************************************
  template< typename VectorType >
  class Conjugate_Gradient : public Operators<VectorType>
  {
  public:
    //--- Convenience stuff for better readability
    typedef VectorType Vector;
    typedef typename Vector::size_type   int_type;  ///< Integer type
    typedef typename Vector::value_type value_type; ///< Scalar type



    Vector &solve(Vector &v, const Vector &b, int_type level)
    {
      if(v.size() == 1) {
        v[0]=0.5*b[0];
        return v;
      }
      int_type n = v.size();
      Vector r(n),p(n),Ap(n);
      value_type beta,alpha,dot_r;
      apply_operator(v,level,r);
      scale(r,-1);
      add(b,r,r);
      //            MG::copy(b,r);
      copy(r,p);
      for (int_type i = 0; i < n; ++i) {
        dot_r = dot(r,r);
        if( dot_r < 0.00001) return v;
        apply_operator(p,level,Ap);
        alpha = dot_r/dot(Ap,p);
        scale(p,alpha);
        add(p,v,v);
        scale(Ap,-alpha);
        add(Ap,r,r);

        beta = dot(r,r)/dot_r;
        scale(p,beta/alpha);
        add(p,r,p);
      }

      return v;
    };
    void testing(){};
  }
  ; // Conjugate_Gradient


    /// My default types
  template <typename Vector>
  struct default_algorithms
  {
    typedef Gauss_Seidel<Vector>       Relaxation;
    typedef Conjugate_Gradient<Vector> Solver;
    typedef Operators<Vector>          Operators;
  };

  /// This class defines the types to be overiden
  template< typename Vector, typename Algorithms = default_algorithms<Vector> >
  class DefaultAlgorithms
  {
  public:
    typedef typename Algorithms::Relaxation relaxation; ///< Default Iterative Solver (Gauss_Seidel)
    typedef typename Algorithms::Solver     solver;     ///< Default solver (Conjugate_Gradient)
    typedef typename Algorithms::Operators  operators;  ///< Default operators (See Operators class)
  };

  /// class to define a use of the DefaultAlgorithms types
  /// avoids ambiguities if we derive from DefaultAlgorithm more than once
  template <typename Vector>
  class DefaultAlgorithmsArgs : virtual public DefaultAlgorithms< Vector >
  {}
  ;

  /// @brief Helper Classes
  /// Usage: Relaxation< Jacobi<Vector> >
  template < typename  Algorithm >
  class Relaxation: virtual public DefaultAlgorithms<typename Algorithm::Vector>
  {
  public:
    typedef Algorithm relaxation; ///< overriding typedef in DefaultAlgorithms
  };

  /// @brief Helper Classes
  /// Usage: Solver< Conjugate_Gradient<Vector> >
  template < typename Algorithm >
  class Solver: virtual public DefaultAlgorithms<typename Algorithm::Vector>
  {
  public:
    typedef Algorithm solver;     ///< overriding typedef in DefaultAlgorithms
  };

  /// @brief Helper Classes
  /// Usage: Operations< Operators<Vector> >
  template < typename Algorithm >
  class Operations: virtual public DefaultAlgorithms<typename Algorithm::Vector>
  {
  public:
    typedef Algorithm operators;  ///< overriding typedef in DefaultAlgorithms
  };

  /// @brief This class set the default base types
  /// It is needed to allow the various Setter types to be identical.
  /// This class is used by AlgorithmSelector.
  /// See http://www.informit.com/articles/article.asp?p=31473&rl=1.
  template <typename Base, int t>
  class  Algorithm_Discriminator : public Base
  {}
  ;

  /// @brief merge the different template arguments into a single type that overrides default typedef members with whichever non-defaults were specified
  template <typename Setter1, typename Setter2, typename Setter3>
  class AlgorithmSelector : public Algorithm_Discriminator<Setter1,1>,
                            public Algorithm_Discriminator<Setter2,2>,
                            public Algorithm_Discriminator<Setter3,3>
  {}
  ;


  ///////////////////////////////////////////////////////////////////////////
  ///
  /// @struct TypeBinder
  /// @brief  This class contain all types used in Multigrid
  ///
  ///////////////////////////////////////////////////////////////////////////
  template<
    typename T,
    class VectorType,
    class MatrixType,
    template< typename > class MultiGridType,
    typename AlgorithmSetter1 = DefaultAlgorithmsArgs<VectorType>,
    typename AlgorithmSetter2 = DefaultAlgorithmsArgs<VectorType>,
    typename AlgorithmSetter3 = DefaultAlgorithmsArgs<VectorType>
  >
  class TypeBinder
  {
  public:
    typedef TypeBinder< T, VectorType, MatrixType, MultiGridType, AlgorithmSetter1, AlgorithmSetter2, AlgorithmSetter3 > Types;
    typedef MultiGridType<Types>   Multigrid;

    typedef AlgorithmSelector< AlgorithmSetter1, AlgorithmSetter2, AlgorithmSetter3 > Algorithms;

    typedef typename Algorithms::relaxation  Relaxation;
    typedef typename Algorithms::solver      Solver;
    typedef typename Algorithms::operators   Operators;
    typedef T                      value_type;    ///< Scalar Type
    typedef unsigned int           int_type;      ///< Integer Type
    typedef MatrixType             Matrix;        ///< nxn Matrix Container
    typedef VectorType             Vector;        ///< Vector Container
    typedef typename Matrix::size_type size_type; ///< Matrix::size_type

    / **
     *  @struct Grid
     *
     *   @brief Contain intergrid vectors and operators
     * /
    struct Grid
    {
      Vector x;   ///< Iteration vector
      Vector b;   ///< Load vector
      Vector r;

      / **
       * @brief Constructor
       * @param dim 
       * @return 
       * /
      Grid(size_type dim)
      {
        Vector x(dim);
        Vector b(dim);
        Vector r(dim);
      }
      / **
       * @brief Empty Constructor
       * @return 
       * /

      Grid() { }

      void resize(size_type new_dim)
      {
        x.resize(new_dim);
        b.resize(new_dim);
        r.resize(new_dim);
      }
    }
    ;

  }
  ; //TypeBinder
  ****************************************************************/ 

} //namespace MY_MG
#endif // MY_MG_H