gmreslinearsolver.h

Go to the documentation of this file.

#ifndef GMRES_LINEAR_SOLVER_H_
#define GMRES_LINEAR_SOLVER_H_

#include <deal.II/lac/vector.h>
#include <deal.II/lac/block_sparsity_pattern.h>
#include <deal.II/lac/block_sparse_matrix.h>
#include <deal.II/lac/compressed_simple_sparsity_pattern.h>
#include <deal.II/lac/solver_gmres.h>
#include <deal.II/lac/precondition.h>
#include <deal.II/lac/full_matrix.h>

#include <deal.II/dofs/dof_tools.h>

#include <deal.II/numerics/vector_tools.h>

#include <vector>

namespace DOpE
{

  template <typename PRECONDITIONER, typename SPARSITYPATTERN, typename MATRIX, typename VECTOR>
    class GMRESLinearSolverWithMatrix
  {
  public:
    GMRESLinearSolverWithMatrix( ParameterReader &param_reader);
    ~GMRESLinearSolverWithMatrix();

    static void declare_params(ParameterReader &param_reader);
    
    template<typename PROBLEM>
      void ReInit(PROBLEM& pde);

    template<typename PROBLEM, typename INTEGRATOR>
      void Solve(PROBLEM& pde,INTEGRATOR& integr, VECTOR &rhs, VECTOR &solution, bool force_matrix_build=false);
     
  protected:

  private:
    SPARSITYPATTERN sparsity_pattern_;
    MATRIX matrix_;

    double linear_global_tol_, linear_tol_;
    int  linear_maxiter_, no_tmp_vectors_;
  };

/*********************************Implementation************************************************/

template <typename PRECONDITIONER,typename SPARSITYPATTERN, typename MATRIX, typename VECTOR>
 void GMRESLinearSolverWithMatrix<PRECONDITIONER,SPARSITYPATTERN,MATRIX,VECTOR>::declare_params(ParameterReader &param_reader)
  {
    param_reader.SetSubsection("gmres_withmatrix parameters");
    param_reader.declare_entry("linear_global_tol", "1.e-10",Patterns::Double(0),"global tolerance for the gmres iteration");
    param_reader.declare_entry("linear_maxiter", "1000",Patterns::Integer(0),"maximal number of gmres steps");
    param_reader.declare_entry("no_tmp_vectors", "100",Patterns::Integer(0),"Number of temporary vectors");
  }
/******************************************************/

  template <typename PRECONDITIONER,typename SPARSITYPATTERN, typename MATRIX, typename VECTOR>
    GMRESLinearSolverWithMatrix<PRECONDITIONER,SPARSITYPATTERN,MATRIX,VECTOR>
    ::GMRESLinearSolverWithMatrix(ParameterReader &param_reader) 
{
  param_reader.SetSubsection("gmres_withmatrix parameters");
  linear_global_tol_ = param_reader.get_double ("linear_global_tol");
  linear_maxiter_    = param_reader.get_integer ("linear_maxiter"); 
  no_tmp_vectors_    = param_reader.get_integer ("no_tmp_vectors"); 

}

/******************************************************/

template <typename PRECONDITIONER,typename SPARSITYPATTERN, typename MATRIX, typename VECTOR>
  GMRESLinearSolverWithMatrix<PRECONDITIONER,SPARSITYPATTERN,MATRIX,VECTOR>::~GMRESLinearSolverWithMatrix()
{
}

/******************************************************/

template <typename PRECONDITIONER,typename SPARSITYPATTERN, typename MATRIX, typename VECTOR>
  template<typename PROBLEM>
  void  GMRESLinearSolverWithMatrix<PRECONDITIONER,SPARSITYPATTERN,MATRIX,VECTOR>::ReInit(PROBLEM& pde)
{
  matrix_.clear();
  pde.ComputeSparsityPattern(sparsity_pattern_);
  matrix_.reinit(sparsity_pattern_);
}

/******************************************************/

template <typename PRECONDITIONER,typename SPARSITYPATTERN, typename MATRIX, typename VECTOR>
  template<typename PROBLEM, typename INTEGRATOR>
  void GMRESLinearSolverWithMatrix<PRECONDITIONER,SPARSITYPATTERN,MATRIX,VECTOR>::Solve(PROBLEM& pde, 
                                                                                        INTEGRATOR& integr,
                                                                                        VECTOR &rhs, 
                                                                                        VECTOR &solution, 
                                                                                        bool force_matrix_build)
{
  if(force_matrix_build)
  { 
    integr.ComputeMatrix (pde,matrix_);
  }
 
  integr.ApplyNewtonBoundaryValues(pde,matrix_,rhs,solution);

  dealii::SolverControl solver_control (linear_maxiter_, linear_global_tol_,false,false);

  // This is gmres specific
  dealii::GrowingVectorMemory<VECTOR> vector_memory;
  typename dealii::SolverGMRES<VECTOR>::AdditionalData gmres_data;
  gmres_data.max_n_tmp_vectors = no_tmp_vectors_;


  dealii::SolverGMRES<VECTOR> gmres (solver_control, vector_memory, gmres_data);
  PRECONDITIONER precondition;
  precondition.initialize(matrix_);
  gmres.solve (matrix_, solution, rhs,
    precondition);

  pde.GetDoFConstraints().distribute(solution);
}


}
#endif