qmrslinearsolver.h

Go to the documentation of this file.

#ifndef _QMRS_LINEAR_SOLVER_H_
#define _QMRS_LINEAR_SOLVER_H_

#include <lac/vector.h>
#include <lac/block_sparsity_pattern.h>
#include <lac/block_sparse_matrix.h>
#include <lac/compressed_simple_sparsity_pattern.h>
#include <lac/solver_qmrs.h>
#include <lac/precondition.h>
#include <lac/full_matrix.h>

#include <dofs/dof_tools.h>

#include <numerics/vector_tools.h>

#include <vector>

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

    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;
  };

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

template <typename PRECONDITIONER,typename SPARSITYPATTERN, typename MATRIX, typename VECTOR>
 void QMRSLinearSolverWithMatrix<PRECONDITIONER,SPARSITYPATTERN,MATRIX,VECTOR>::declare_params(ParameterReader &param_reader)
  {
    param_reader.SetSubsection("qmrslinearsolver_withmatrix parameters");
    param_reader.declare_entry("linear_global_tol", "1.e-16",Patterns::Double(0),"global tolerance for the cg iteration");
    param_reader.declare_entry("linear_tol", "1.e-12",Patterns::Double(0),"relative tolerance for the cg iteration");
    param_reader.declare_entry("linear_maxiter", "1000",Patterns::Integer(0),"maximal number of cg steps");
  }
/******************************************************/

  template <typename PRECONDITIONER,typename SPARSITYPATTERN, typename MATRIX, typename VECTOR>
    QMRSLinearSolverWithMatrix<PRECONDITIONER,SPARSITYPATTERN,MATRIX,VECTOR>
    ::QMRSLinearSolverWithMatrix(ParameterReader &param_reader) 
{
  param_reader.SetSubsection("qmrslinearsolver_withmatrix parameters");
  _linear_global_tol = param_reader.get_double ("linear_global_tol");
  _linear_tol        = param_reader.get_double ("linear_tol"); 
  _linear_maxiter    = param_reader.get_integer ("linear_maxiter"); 

}

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

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

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

template <typename PRECONDITIONER,typename SPARSITYPATTERN, typename MATRIX, typename VECTOR>
  template<typename PROBLEM>
  void  QMRSLinearSolverWithMatrix<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 QMRSLinearSolverWithMatrix<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,true);//letzte Arg = false!
  dealii::SolverQMRS<VECTOR> qmres (solver_control);
  PRECONDITIONER precondition;
  precondition.initialize(_matrix);
  qmres.solve (_matrix, solution, rhs,
            precondition);
  VECTOR tmp(solution.size());
  _matrix.vmult(tmp,solution);
  tmp-= rhs;
  std::cout<<"XXX"<<tmp.linfty_norm()<<" ---- "<<tmp.l2_norm()<<std::endl;
  pde.GetDoFConstraints().distribute(solution);
}
/******************************************************/



}
#endif