24 #ifndef QMRS_LINEAR_SOLVER_H_
25 #define QMRS_LINEAR_SOLVER_H_
27 #include <deal.II/lac/vector.h>
28 #include <deal.II/lac/block_sparsity_pattern.h>
29 #include <deal.II/lac/block_sparse_matrix.h>
30 #include <deal.II/lac/compressed_simple_sparsity_pattern.h>
31 #include <deal.II/lac/solver_qmrs.h>
32 #include <deal.II/lac/precondition.h>
33 #include <deal.II/lac/full_matrix.h>
34 #include <deal.II/dofs/dof_tools.h>
35 #include <deal.II/numerics/vector_tools.h>
54 template <
typename PRECONDITIONER,
typename SPARSITYPATTERN,
typename MATRIX,
typename VECTOR>
67 template<
typename PROBLEM>
88 template<
typename PROBLEM,
typename INTEGRATOR>
89 void Solve(PROBLEM& pde, INTEGRATOR& integr, VECTOR &rhs, VECTOR &solution,
bool force_matrix_build=
false);
94 SPARSITYPATTERN sparsity_pattern_;
97 double linear_global_tol_, linear_tol_;
103 template <
typename PRECONDITIONER,
typename SPARSITYPATTERN,
typename MATRIX,
typename VECTOR>
106 param_reader.
SetSubsection(
"qmrslinearsolver_withmatrix parameters");
107 param_reader.
declare_entry(
"linear_global_tol",
"1.e-16",Patterns::Double(0),
"global tolerance for the cg iteration");
108 param_reader.
declare_entry(
"linear_tol",
"1.e-12",Patterns::Double(0),
"relative tolerance for the cg iteration");
109 param_reader.
declare_entry(
"linear_maxiter",
"1000",Patterns::Integer(0),
"maximal number of cg steps");
113 template <
typename PRECONDITIONER,
typename SPARSITYPATTERN,
typename MATRIX,
typename VECTOR>
117 param_reader.
SetSubsection(
"qmrslinearsolver_withmatrix parameters");
118 linear_global_tol_ = param_reader.
get_double (
"linear_global_tol");
119 linear_tol_ = param_reader.
get_double (
"linear_tol");
120 linear_maxiter_ = param_reader.
get_integer (
"linear_maxiter");
126 template <
typename PRECONDITIONER,
typename SPARSITYPATTERN,
typename MATRIX,
typename VECTOR>
133 template <
typename PRECONDITIONER,
typename SPARSITYPATTERN,
typename MATRIX,
typename VECTOR>
134 template<
typename PROBLEM>
138 pde.ComputeSparsityPattern(sparsity_pattern_);
139 matrix_.reinit(sparsity_pattern_);
143 template <
typename PRECONDITIONER,
typename SPARSITYPATTERN,
typename MATRIX,
typename VECTOR>
144 template<
typename PROBLEM,
typename INTEGRATOR>
149 bool force_matrix_build)
151 if(force_matrix_build)
153 integr.ComputeMatrix (pde,matrix_);
156 integr.ApplyNewtonBoundaryValues(pde,matrix_,rhs,solution);
158 dealii::SolverControl solver_control (linear_maxiter_, linear_global_tol_,
false,
true);
159 dealii::SolverQMRS<VECTOR> qmres (solver_control);
160 PRECONDITIONER precondition;
161 precondition.initialize(matrix_);
162 qmres.solve (matrix_, solution, rhs,
164 VECTOR tmp(solution.size());
165 matrix_.vmult(tmp,solution);
167 std::cout<<
"XXX"<<tmp.linfty_norm()<<
" ---- "<<tmp.l2_norm()<<std::endl;
168 pde.GetDoFConstraints().distribute(solution);
double get_double(const std::string &entry_name)
Definition: parameterreader.h:115
void Solve(PROBLEM &pde, INTEGRATOR &integr, VECTOR &rhs, VECTOR &solution, bool force_matrix_build=false)
Definition: qmrslinearsolver.h:145
Definition: parameterreader.h:36
Definition: qmrslinearsolver.h:55
void declare_entry(const std::string &entry, const std::string &default_value, const Patterns::PatternBase &pattern=Patterns::Anything(), const std::string &documentation=std::string())
Definition: parameterreader.h:98
static void declare_params(ParameterReader ¶m_reader)
Definition: qmrslinearsolver.h:104
void ReInit(PROBLEM &pde)
Definition: qmrslinearsolver.h:135
int get_integer(const std::string &entry_name)
Definition: parameterreader.h:126
void SetSubsection(const std::string subsection)
Definition: parameterreader.h:93
~QMRSLinearSolverWithMatrix()
Definition: qmrslinearsolver.h:127
QMRSLinearSolverWithMatrix(ParameterReader ¶m_reader)
Definition: qmrslinearsolver.h:115