diff --git a/.gitignore b/.gitignore
index d896e231f34910cb510864a1956a5904a0a628ba..d4c8ca60f31914b314371580fa2e0d4928656d94 100644
--- a/.gitignore
+++ b/.gitignore
@@ -26,6 +26,7 @@ ltmain.sh
am
.libs/
.deps/
+*~
*.o
test-driver
build-cmake/
diff --git a/dune/istl/CMakeLists.txt b/dune/istl/CMakeLists.txt
index 9209396cdbad0480d4c5e981485a79fb76e6b7e1..1f288b77a66664d1f13702b3358ab98664b02d02 100644
--- a/dune/istl/CMakeLists.txt
+++ b/dune/istl/CMakeLists.txt
@@ -28,12 +28,14 @@ install(FILES
overlappingschwarz.hh
owneroverlapcopy.hh
pardiso.hh
+ preconditoner.hh
preconditioners.hh
repartition.hh
scalarproducts.hh
scaledidmatrix.hh
schwarz.hh
solvercategory.hh
+ solver.hh
solvers.hh
solvertype.hh
superlu.hh
diff --git a/dune/istl/Makefile.am b/dune/istl/Makefile.am
index 8c2474f89520a2081897db6ccf5e9b719b9abe41..f36255c2653eb53d16d23bf2212efb7fbaa413c3 100644
--- a/dune/istl/Makefile.am
+++ b/dune/istl/Makefile.am
@@ -27,12 +27,14 @@ istl_HEADERS = basearray.hh \
overlappingschwarz.hh \
owneroverlapcopy.hh \
pardiso.hh \
+ preconditioner.hh \
preconditioners.hh \
repartition.hh \
scalarproducts.hh \
scaledidmatrix.hh \
schwarz.hh \
solvercategory.hh \
+ solver.hh \
solvers.hh \
solvertype.hh \
superlu.hh \
diff --git a/dune/istl/paamg/aggregates.hh b/dune/istl/paamg/aggregates.hh
index 769a76023302dd692add7ea16b37a52a46743777..e5005386f6e59aea6e435a4f1d19dbf2a9e344bd 100644
--- a/dune/istl/paamg/aggregates.hh
+++ b/dune/istl/paamg/aggregates.hh
@@ -1504,8 +1504,9 @@ namespace Dune
template<class G,class S>
inline void Aggregate<G,S>::add(std::vector<Vertex>& vertices)
{
+#ifndef NDEBUG
std::size_t oldsize = vertices_.size();
-
+#endif
typedef typename std::vector<Vertex>::iterator Iterator;
typedef typename VertexSet::iterator SIterator;
diff --git a/dune/istl/preconditioner.hh b/dune/istl/preconditioner.hh
new file mode 100644
index 0000000000000000000000000000000000000000..79452c3e01e287d75c47c87cb940b9cb763f6064
--- /dev/null
+++ b/dune/istl/preconditioner.hh
@@ -0,0 +1,82 @@
+// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
+// vi: set et ts=4 sw=2 sts=2:
+#ifndef DUNE_PRECONDITIONER_HH
+#define DUNE_PRECONDITIONER_HH
+namespace Dune {
+/**
+ * @addtogroup ISTL_Prec
+ * @{
+ */
+ //=====================================================================
+ /*! \brief Base class for matrix free definition of preconditioners.
+
+ Note that the operator, which is the basis for the preconditioning,
+ is supplied to the preconditioner from the outside in the
+ constructor or some other method.
+
+ This interface allows the encapsulation of all parallelization
+ aspects into the preconditioners.
+
+ \tparam X Type of the update
+ \tparam Y Type of the defect
+
+ */
+ //=====================================================================
+ template<class X, class Y>
+ class Preconditioner {
+ public:
+ //! \brief The domain type of the preconditioner.
+ typedef X domain_type;
+ //! \brief The range type of the preconditioner.
+ typedef Y range_type;
+ //! \brief The field type of the preconditioner.
+ typedef typename X::field_type field_type;
+
+ /*! \brief Prepare the preconditioner.
+
+ A solver solves a linear operator equation A(x)=b by applying
+ one or several steps of the preconditioner. The method pre()
+ is called before the first apply operation.
+ b and x are right hand side and solution vector of the linear
+ system respectively. It may. e.g., scale the system, allocate memory or
+ compute a (I)LU decomposition.
+ Note: The ILU decomposition could also be computed in the constructor
+ or with a separate method of the derived method if several
+ linear systems with the same matrix are to be solved.
+
+ \param x The left hand side of the equation.
+ \param b The right hand side of the equation.
+ */
+ virtual void pre (X& x, Y& b) = 0;
+
+ /*! \brief Apply one step of the preconditioner to the system A(v)=d.
+
+ On entry v=0 and d=b-A(x) (although this might not be
+ computed in that way. On exit v contains the update, i.e
+ one step computes \f$ v = M^{-1} d \f$ where \f$ M \f$ is the
+ approximate inverse of the operator \f$ A \f$ characterizing
+ the preconditioner.
+ \param[out] v The update to be computed
+ \param d The current defect.
+ */
+ virtual void apply (X& v, const Y& d) = 0;
+
+ /*! \brief Clean up.
+
+ This method is called after the last apply call for the
+ linear system to be solved. Memory may be deallocated safely
+ here. x is the solution of the linear equation.
+
+ \param x The right hand side of the equation.
+ */
+ virtual void post (X& x) = 0;
+
+ // every abstract base class has a virtual destructor
+ virtual ~Preconditioner () {}
+ };
+
+/**
+ * @}
+ */
+}
+#endif
diff --git a/dune/istl/preconditioners.hh b/dune/istl/preconditioners.hh
index 5cb894beeeae3c30d85b6ef7fc795e467bbe8657..7db4c775b31281a20405f9c317a0f11539760086 100644
--- a/dune/istl/preconditioners.hh
+++ b/dune/istl/preconditioners.hh
@@ -9,6 +9,8 @@
#include <iomanip>
#include <string>
+#include "preconditioner.hh"
+#include "solver.hh"
#include "solvercategory.hh"
#include "istlexception.hh"
#include "matrixutils.hh"
@@ -54,75 +56,55 @@ namespace Dune {
*/
- //=====================================================================
- /*! \brief Base class for matrix free definition of preconditioners.
-
- Note that the operator, which is the basis for the preconditioning,
- is supplied to the preconditioner from the outside in the
- constructor or some other method.
-
- This interface allows the encapsulation of all parallelization
- aspects into the preconditioners.
-
- \tparam X Type of the update
- \tparam Y Type of the defect
+ /**
+ * @brief Turns an InverseOperator into a Preconditioner.
+ * @tparam O The type of the inverse operator to wrap.
*/
- //=====================================================================
- template<class X, class Y>
- class Preconditioner {
+ template<class O, int c>
+ class InverseOperator2Preconditioner :
+ public Preconditioner<typename O::domain_type, typename O::range_type>
+ {
public:
//! \brief The domain type of the preconditioner.
- typedef X domain_type;
+ typedef typename O::domain_type domain_type;
//! \brief The range type of the preconditioner.
- typedef Y range_type;
+ typedef typename O::range_type range_type;
//! \brief The field type of the preconditioner.
- typedef typename X::field_type field_type;
-
- /*! \brief Prepare the preconditioner.
+ typedef typename range_type::field_type field_type;
+ typedef O InverseOperator;
- A solver solves a linear operator equation A(x)=b by applying
- one or several steps of the preconditioner. The method pre()
- is called before the first apply operation.
- b and x are right hand side and solution vector of the linear
- system respectively. It may. e.g., scale the system, allocate memory or
- compute a (I)LU decomposition.
- Note: The ILU decomposition could also be computed in the constructor
- or with a separate method of the derived method if several
- linear systems with the same matrix are to be solved.
-
- \param x The left hand side of the equation.
- \param b The right hand side of the equation.
- */
- virtual void pre (X& x, Y& b) = 0;
-
- /*! \brief Apply one step of the preconditioner to the system A(v)=d.
+ // define the category
+ enum {
+ //! \brief The category the preconditioner is part of.
+ category=c
+ };
- On entry v=0 and d=b-A(x) (although this might not be
- computed in that way. On exit v contains the update, i.e
- one step computes \f$ v = M^{-1} d \f$ where \f$ M \f$ is the
- approximate inverse of the operator \f$ A \f$ characterizing
- the preconditioner.
- \param[out] v The update to be computed
- \param d The current defect.
+ /**
+ * @brief Construct the preconditioner from the solver
+ * @param inverse_operator The inverse operator to wrap.
*/
- virtual void apply (X& v, const Y& d) = 0;
+ InverseOperator2Preconditioner(InverseOperator& inverse_operator)
+ : inverse_operator_(inverse_operator)
+ {}
- /*! \brief Clean up.
+ void pre(domain_type&,range_type&)
+ {}
- This method is called after the last apply call for the
- linear system to be solved. Memory may be deallocated safely
- here. x is the solution of the linear equation.
+ void apply(domain_type& v, const range_type& d)
+ {
+ InverseOperatorResult res;
+ range_type copy(d);
+ inverse_operator_.apply(v, copy, res);
+ }
- \param x The right hand side of the equation.
- */
- virtual void post (X& x) = 0;
+ void post(domain_type&)
+ {}
- // every abstract base class has a virtual destructor
- virtual ~Preconditioner () {}
+ private:
+ InverseOperator& inverse_operator_;
};
-
//=====================================================================
// Implementation of this interface for sequential ISTL-preconditioners
//=====================================================================
diff --git a/dune/istl/solver.hh b/dune/istl/solver.hh
new file mode 100644
index 0000000000000000000000000000000000000000..8ace9496b486558c489567bdfbbf8af0e844c7f0
--- /dev/null
+++ b/dune/istl/solver.hh
@@ -0,0 +1,157 @@
+// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
+// vi: set et ts=4 sw=2 sts=2:
+
+#ifndef DUNE_SOLVER_HH
+#define DUNE_SOLVER_HH
+
+#include <iomanip>
+#include <ostream>
+
+namespace Dune
+{
+/**
+ * @addtogroup ISTL_Solvers
+ * @{
+ */
+/** \file
+
+ \brief Define general, extensible interface for inverse operators.
+
+ Implementation here covers only inversion of linear operators,
+ but the implementation might be used for nonlinear operators
+ as well.
+ */
+ /**
+ \brief Statistics about the application of an inverse operator
+
+ The return value of an application of the inverse
+ operator delivers some important information about
+ the iteration.
+ */
+ struct InverseOperatorResult
+ {
+ /** \brief Default constructor */
+ InverseOperatorResult ()
+ {
+ clear();
+ }
+
+ /** \brief Resets all data */
+ void clear ()
+ {
+ iterations = 0;
+ reduction = 0;
+ converged = false;
+ conv_rate = 1;
+ elapsed = 0;
+ }
+
+ /** \brief Number of iterations */
+ int iterations;
+
+ /** \brief Reduction achieved: \f$ \|b-A(x^n)\|/\|b-A(x^0)\|\f$ */
+ double reduction;
+
+ /** \brief True if convergence criterion has been met */
+ bool converged;
+
+ /** \brief Convergence rate (average reduction per step) */
+ double conv_rate;
+
+ /** \brief Elapsed time in seconds */
+ double elapsed;
+ };
+
+
+ //=====================================================================
+ /*!
+ \brief Abstract base class for all solvers.
+
+ An InverseOperator computes the solution of \f$ A(x)=b\f$ where
+ \f$ A : X \to Y \f$ is an operator.
+ Note that the solver "knows" which operator
+ to invert and which preconditioner to apply (if any). The
+ user is only interested in inverting the operator.
+ InverseOperator might be a Newton scheme, a Krylov subspace method,
+ or a direct solver or just anything.
+ */
+ template<class X, class Y>
+ class InverseOperator {
+ public:
+ //! \brief Type of the domain of the operator to be inverted.
+ typedef X domain_type;
+
+ //! \brief Type of the range of the operator to be inverted.
+ typedef Y range_type;
+
+ /** \brief The field type of the operator. */
+ typedef typename X::field_type field_type;
+
+ /**
+ \brief Apply inverse operator,
+
+ \warning Note: right hand side b may be overwritten!
+
+ \param x The left hand side to store the result in.
+ \param b The right hand side
+ \param res Object to store the statistics about applying the operator.
+ */
+ virtual void apply (X& x, Y& b, InverseOperatorResult& res) = 0;
+
+ /*!
+ \brief apply inverse operator, with given convergence criteria.
+
+ \warning Right hand side b may be overwritten!
+
+ \param x The left hand side to store the result in.
+ \param b The right hand side
+ \param reduction The minimum defect reduction to achieve.
+ \param res Object to store the statistics about applying the operator.
+ */
+ virtual void apply (X& x, Y& b, double reduction, InverseOperatorResult& res) = 0;
+
+ //! \brief Destructor
+ virtual ~InverseOperator () {}
+
+ protected:
+ // spacing values
+ enum { iterationSpacing = 5 , normSpacing = 16 };
+
+ //! helper function for printing header of solver output
+ void printHeader(std::ostream& s) const
+ {
+ s << std::setw(iterationSpacing) << " Iter";
+ s << std::setw(normSpacing) << "Defect";
+ s << std::setw(normSpacing) << "Rate" << std::endl;
+ }
+
+ //! helper function for printing solver output
+ template <class DataType>
+ void printOutput(std::ostream& s,
+ const double iter,
+ const DataType& norm,
+ const DataType& norm_old) const
+ {
+ const DataType rate = norm/norm_old;
+ s << std::setw(iterationSpacing) << iter << " ";
+ s << std::setw(normSpacing) << norm << " ";
+ s << std::setw(normSpacing) << rate << std::endl;
+ }
+
+ //! helper function for printing solver output
+ template <class DataType>
+ void printOutput(std::ostream& s,
+ const double iter,
+ const DataType& norm) const
+ {
+ s << std::setw(iterationSpacing) << iter << " ";
+ s << std::setw(normSpacing) << norm << std::endl;
+ }
+ };
+
+/**
+ * @}
+ */
+}
+
+#endif
diff --git a/dune/istl/solvers.hh b/dune/istl/solvers.hh
index 64be4dd674a1959bee6a2804937cfccf3462abe6..27f191824ae9bd1a87ea5341fa65415e937d470f 100644
--- a/dune/istl/solvers.hh
+++ b/dune/istl/solvers.hh
@@ -9,11 +9,13 @@
#include <iostream>
#include <iomanip>
#include <string>
+#include <vector>
#include "istlexception.hh"
#include "operators.hh"
-#include "preconditioners.hh"
#include "scalarproducts.hh"
+#include "solver.hh"
+#include "preconditioner.hh"
#include <dune/common/timer.hh>
#include <dune/common/ftraits.hh>
#include <dune/common/shared_ptr.hh>
@@ -30,143 +32,14 @@ namespace Dune {
/** \file
- \brief Define general, extensible interface for inverse operators.
+ \brief Implementations of the inverse operator interface.
- Implementation here covers only inversion of linear operators,
- but the implementation might be used for nonlinear operators
- as well.
+ This file provides various preconditioned Krylov methods.
*/
- /**
- \brief Statistics about the application of an inverse operator
-
- The return value of an application of the inverse
- operator delivers some important information about
- the iteration.
- */
- struct InverseOperatorResult
- {
- /** \brief Default constructor */
- InverseOperatorResult ()
- {
- clear();
- }
-
- /** \brief Resets all data */
- void clear ()
- {
- iterations = 0;
- reduction = 0;
- converged = false;
- conv_rate = 1;
- elapsed = 0;
- }
-
- /** \brief Number of iterations */
- int iterations;
-
- /** \brief Reduction achieved: \f$ \|b-A(x^n)\|/\|b-A(x^0)\|\f$ */
- double reduction;
-
- /** \brief True if convergence criterion has been met */
- bool converged;
-
- /** \brief Convergence rate (average reduction per step) */
- double conv_rate;
-
- /** \brief Elapsed time in seconds */
- double elapsed;
- };
-
-
- //=====================================================================
- /*!
- \brief Abstract base class for all solvers.
-
- An InverseOperator computes the solution of \f$ A(x)=b\f$ where
- \f$ A : X \to Y \f$ is an operator.
- Note that the solver "knows" which operator
- to invert and which preconditioner to apply (if any). The
- user is only interested in inverting the operator.
- InverseOperator might be a Newton scheme, a Krylov subspace method,
- or a direct solver or just anything.
- */
- template<class X, class Y>
- class InverseOperator {
- public:
- //! \brief Type of the domain of the operator to be inverted.
- typedef X domain_type;
-
- //! \brief Type of the range of the operator to be inverted.
- typedef Y range_type;
-
- /** \brief The field type of the operator. */
- typedef typename X::field_type field_type;
-
- /**
- \brief Apply inverse operator,
- \warning Note: right hand side b may be overwritten!
- \param x The left hand side to store the result in.
- \param b The right hand side
- \param res Object to store the statistics about applying the operator.
- */
- virtual void apply (X& x, Y& b, InverseOperatorResult& res) = 0;
-
- /*!
- \brief apply inverse operator, with given convergence criteria.
-
- \warning Right hand side b may be overwritten!
-
- \param x The left hand side to store the result in.
- \param b The right hand side
- \param reduction The minimum defect reduction to achieve.
- \param res Object to store the statistics about applying the operator.
- */
- virtual void apply (X& x, Y& b, double reduction, InverseOperatorResult& res) = 0;
-
- //! \brief Destructor
- virtual ~InverseOperator () {}
-
- protected:
- // spacing values
- enum { iterationSpacing = 5 , normSpacing = 16 };
-
- //! helper function for printing header of solver output
- void printHeader(std::ostream& s) const
- {
- s << std::setw(iterationSpacing) << " Iter";
- s << std::setw(normSpacing) << "Defect";
- s << std::setw(normSpacing) << "Rate" << std::endl;
- }
-
- //! helper function for printing solver output
- template <class DataType>
- void printOutput(std::ostream& s,
- const double iter,
- const DataType& norm,
- const DataType& norm_old) const
- {
- const DataType rate = norm/norm_old;
- s << std::setw(iterationSpacing) << iter << " ";
- s << std::setw(normSpacing) << norm << " ";
- s << std::setw(normSpacing) << rate << std::endl;
- }
-
- //! helper function for printing solver output
- template <class DataType>
- void printOutput(std::ostream& s,
- const double iter,
- const DataType& norm) const
- {
- s << std::setw(iterationSpacing) << iter << " ";
- s << std::setw(normSpacing) << norm << std::endl;
- }
- };
-
-
- //=====================================================================
+ //=====================================================================
// Implementation of this interface
//=====================================================================
@@ -302,6 +175,9 @@ namespace Dune {
}
}
+ //correct i which is wrong if convergence was not achieved.
+ i=std::min(_maxit,i);
+
// print
if (_verbose==1)
this->printOutput(std::cout,i,def);
@@ -437,6 +313,9 @@ namespace Dune {
}
}
+ //correct i which is wrong if convergence was not achieved.
+ i=std::min(_maxit,i);
+
if (_verbose==1) // printing for non verbose
this->printOutput(std::cout,i,def);
@@ -598,6 +477,9 @@ namespace Dune {
rholast = rho; // remember rho for recurrence
}
+ //correct i which is wrong if convergence was not achieved.
+ i=std::min(_maxit,i);
+
if (_verbose==1) // printing for non verbose
this->printOutput(std::cout,i,def);
@@ -864,6 +746,9 @@ namespace Dune {
norm_old = norm;
} // end for
+ //correct i which is wrong if convergence was not achieved.
+ it=std::min((double)_maxit,it);
+
if (_verbose==1) // printing for non verbose
this->printOutput(std::cout,it,norm);
@@ -1110,6 +995,9 @@ namespace Dune {
}
}
+ //correct i which is wrong if convergence was not achieved.
+ i=std::min(_maxit,i);
+
if (_verbose==1) // printing for non verbose
this->printOutput(std::cout,i,def);
@@ -1369,6 +1257,9 @@ namespace Dune {
res.converged = false;
}
+ //correct i which is wrong if convergence was not achieved.
+ j=std::min(_maxit,j);
+
if (_verbose > 1) // print
{
this->printOutput(std::cout,j,norm,norm_old);
diff --git a/dune/istl/test/CMakeLists.txt b/dune/istl/test/CMakeLists.txt
index 9ea18d131f35caf6627856c9ca6fe91b1618ad86..60b3dfed0e2f741925ed79f5d91883128ed061e3 100644
--- a/dune/istl/test/CMakeLists.txt
+++ b/dune/istl/test/CMakeLists.txt
@@ -4,6 +4,7 @@ set(NORMALTEST
bcrsbuildtest
dotproducttest
iotest
+ inverseoperator2prectest
matrixiteratortest
matrixtest
matrixutilstest
@@ -50,6 +51,7 @@ add_executable(matrixiteratortest "matrixiteratortest.cc")
add_executable(mmtest mmtest.cc)
add_executable(mv "mv.cc")
add_executable(iotest "iotest.cc")
+add_executable(inverseoperator2prectest "inverseoperator2prectest.cc")
add_executable(scaledidmatrixtest "scaledidmatrixtest.cc")
add_executable(seqmatrixmarkettest "matrixmarkettest.cc")
#set_target_properties(seqmatrixmarkettest PROPERTIES COMPILE_FLAGS
diff --git a/dune/istl/test/Makefile.am b/dune/istl/test/Makefile.am
index f6d05b0759018b6c76029f1cea5d699b6b976cba..c0de4632f4cd9f65b52592e50e483bed7f994cd0 100644
--- a/dune/istl/test/Makefile.am
+++ b/dune/istl/test/Makefile.am
@@ -23,6 +23,7 @@ NORMALTESTS = basearraytest \
complexrhstest \
dotproducttest \
iotest \
+ inverseoperator2prectest \
matrixiteratortest \
matrixtest \
matrixutilstest \
@@ -103,6 +104,8 @@ mv_SOURCES = mv.cc
iotest_SOURCES = iotest.cc
+inverseoperator2prectest_SOURCES = inverseoperator2prectest.cc
+
scaledidmatrixtest_SOURCES = scaledidmatrixtest.cc
if MPI
diff --git a/dune/istl/test/complexrhstest.cc b/dune/istl/test/complexrhstest.cc
index 0ab63b801d4d9edce19d30c757587ddd3b2edff2..d1a9ee073e530c9afed72a7842a7985123b934cd 100644
--- a/dune/istl/test/complexrhstest.cc
+++ b/dune/istl/test/complexrhstest.cc
@@ -29,6 +29,7 @@
#include <dune/istl/bvector.hh>
#include <dune/istl/operators.hh>
#include <dune/istl/solvers.hh>
+#include <dune/istl/preconditioners.hh>
#include "laplacian.hh"
#if HAVE_SUPERLU
diff --git a/dune/istl/test/inverseoperator2prectest.cc b/dune/istl/test/inverseoperator2prectest.cc
new file mode 100644
index 0000000000000000000000000000000000000000..595b137432cce8d933c273cb915622b569c46821
--- /dev/null
+++ b/dune/istl/test/inverseoperator2prectest.cc
@@ -0,0 +1,61 @@
+// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
+// vi: set et ts=4 sw=2 sts=2:
+#include "config.h"
+#include<dune/istl/bvector.hh>
+#include<dune/istl/superlu.hh>
+#include<dune/istl/preconditioners.hh>
+#include<dune/istl/solvers.hh>
+#include<dune/istl/operators.hh>
+#include <dune/common/fmatrix.hh>
+#include <dune/common/fvector.hh>
+#include <laplacian.hh>
+
+int main(int argc, char** argv)
+{
+ const int BS=1;
+ int N=100;
+
+ if(argc>1)
+ N = atoi(argv[1]);
+ std::cout<<"testing for N="<<N<<" BS="<<1<<std::endl;
+
+ typedef Dune::FieldMatrix<double,BS,BS> MatrixBlock;
+ typedef Dune::BCRSMatrix<MatrixBlock> BCRSMat;
+ typedef Dune::FieldVector<double,BS> VectorBlock;
+ typedef Dune::BlockVector<VectorBlock> BVector;
+ typedef Dune::MatrixAdapter<BCRSMat,BVector,BVector> Operator;
+
+ BCRSMat mat;
+ Operator fop(mat);
+ BVector b(N*N), x(N*N);
+
+ setupLaplacian(mat,N);
+ b=0;
+ x=100;
+
+ Dune::InverseOperatorResult res, res1;
+ x=1;
+ mat.mv(x, b);
+ x=0;
+ Dune::SeqJac<BCRSMat,BVector,BVector> prec0(mat, 1,1.0);
+ const int category = Dune::SeqJac<BCRSMat,BVector,BVector>::category;
+ Dune::LoopSolver<BVector> solver0(fop, prec0, 1e-3,10,0);
+ Dune::InverseOperator2Preconditioner<Dune::LoopSolver<BVector>,category >
+ prec(solver0);
+ Dune::LoopSolver<BVector> solver(fop, prec, 1e-8,10,2);
+ solver.apply(x,b,res);
+
+ x=1;
+ mat.mv(x, b);
+ x=0;
+ std::cout<<"solver1"<<std::endl;
+ Dune::LoopSolver<BVector> solver1(fop, prec0, 1e-8,100,2);
+ solver1.apply(x,b,res1);
+
+ if(res1.iterations!=res.iterations*10)
+ {
+ std::cerr<<"Convergence rates do not match!"<<std::endl;
+ return 1;
+ }
+ return 0;
+}