/******************************************************************************
 *
 * AMDiS - Adaptive multidimensional simulations
 *
 * Copyright (C) 2013 Dresden University of Technology. All Rights Reserved.
 * Web: https://fusionforge.zih.tu-dresden.de/projects/amdis
 *
 * Authors: 
 * Simon Vey, Thomas Witkowski, Andreas Naumann, Simon Praetorius, et al.
 *
 * This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE
 * WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE.
 *
 *
 * This file is part of AMDiS
 *
 * See also license.opensource.txt in the distribution.
 * 
 ******************************************************************************/


/** \file PetscSolver.h */

#ifndef AMDIS_PETSC_SOLVER_SEQ_H
#define AMDIS_PETSC_SOLVER_SEQ_H

#ifdef HAVE_SEQ_PETSC

#include "solver/LinearSolver.h"
#include "solver/MatrixStreams.h"
#include "solver/PetscTypes.h"
#include "Timer.h"
#include <vector>
#include <iostream>
#include <boost/mpl/bool.hpp>

#include <petsc.h>
#include <petscksp.h>
#include <petscmat.h> 
#include <petscvec.h>
#include <petscsys.h>
#include <petscao.h>

namespace AMDiS {
  
  class PetscRunner : public OEMRunner
  {
  public:
    /// Constructor of standard PETSc runner. Reads ksp and pc parameters from initfile.
    PetscRunner(LinearSolver* oem_);

    typedef SolverMatrix<Matrix<DOFMatrix*> > BlockMatrix;
    
    /// initialize the solver \ref ksp and preconditioner \ref pc
    void init(const BlockMatrix& A, const Mat& fullMatrix);
      
    /// solve the linear equation \f$ A\cdot x = b \f$ by applying the PETSc solver \ref ksp
    int solve(const Mat& A, Vec& x, const Vec& b);
    
    /// destroy solver \ref ksp and preconditioner \ref pc
    virtual void exit()
    {
      KSPDestroy(&ksp);
    }

    /// get the PETSc solver \ref ksp
    KSP getSolver() 
    { 
      return ksp; 
    }

    /// get the PETSc preconditioner \ref pc
    PC getPc() 
    { 
      return pc; 
    }

    ~PetscRunner() { }
    
  protected:    
    static void createSubSolver(KSP &ksp_, Mat m, std::string kspPrefix_);
    
    static void setSolver(KSP ksp_, std::string kspPrefix_,
			  KSPType kspType, PCType pcType, 
			  PetscReal rtol = PETSC_DEFAULT,
			  PetscReal atol = PETSC_DEFAULT,
			  PetscInt maxIt = PETSC_DEFAULT);

    LinearSolver& oem;
    
  //private:    
    /// PETSc solver object
    KSP ksp;

    /// PETSc preconditioner object
    PC pc;

    /// KSP database prefix
    std::string kspPrefix;
    
    bool zeroStartVector;
  };

  
  /**
   * \ingroup Solver
   * 
   * \brief
   * Wrapper for the external PETSc solver:
   *   http://www.mcs.anl.gov/petsc/
   *
   * This is a suite of data structures and routines for the 
   * scalable (parallel) solution of scientific applications.
   */
  template< typename Runner = PetscRunner >
  class PetscSolver : public LinearSolver 
  {
  public:
    /// Creator class used in the LinearSolverMap.
    class Creator : public LinearSolverCreator
    {
    public:
      virtual ~Creator() {}
      
      /// Returns a new PetscSolver object.
      LinearSolver* create() 
      { 
	return new PetscSolver<Runner>(this->name); 
      }
    };
    
    PetscSolver(std::string n)
    : LinearSolver(n),
      runner(this)
    {}
    
    ~PetscSolver() {}
    
    virtual OEMRunner* getRunner()
    {
      return &runner;
    }
    
    void setNestedVectors(bool nested = true)
    {
      vecSol.isNested = nested;
      vecRhs.isNested = nested;
    }
 
    void setNestedMatrix(bool nested = true)
    {
      petscMat.isNested=nested;
    } 

    void setNested(bool n)
    {
      setNestedVectors(n); 
      setNestedMatrix(n); 
    }
    
  protected:    
    Runner runner;

    int solveLinearSystem(const SolverMatrix<Matrix<DOFMatrix*> >& A,
			  SystemVector& x, 
			  SystemVector& b,
			  bool createMatrixData,
			  bool storeMatrixData) 
    { FUNCNAME("PetscSolver::solveLinearSystem()");

      Timer t;
      // transfer matrix and rhs-vector to PETSc data-structures
      if (createMatrixData) {
	petscMat << A;
	runner.init(A, petscMat.matrix);
      }

      vecSol << x;
      vecRhs << b;
      INFO(info, 8)("fill PETSc matrix needed %.5f seconds\n", t.elapsed());

      // solve the linear system using PETSc solvers
      t.reset();
      error = runner.solve(petscMat.matrix, vecSol.vector, vecRhs.vector);

      // transfer solution from PETSc vector to SystemVector
      vecSol.vector >> x;

      vecSol.destroy();
      vecRhs.destroy();
      if (!storeMatrixData){
	petscMat.destroy();
	runner.exit();
      }
      return error;
    }
    
  private:
    /// PETSc System-Matrix
    PetscMatrix petscMat;

    /// Solution and RHS vectors.
    PetscVector vecSol, vecRhs;
  };
  
} // end namespace AMDiS

#endif // HAVE_SEQ_PETSC

#endif // AMDIS_PETSC_SOLVER_SEQ_H