diff --git a/CHANGELOG.md b/CHANGELOG.md
index 01150cf1e238f4fc3e89f4ecaea26267c1374133..186604ce2f493b0087d6d7897978e6f2d256e56d 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -1,5 +1,9 @@
 # Master (will become release 2.9)
 
+- Add `const` qualifier to `LinearOperator` and `ScalarProduct` in
+  `IterativeSolver`. In particular, the constructors of iterative solvers have
+  changed.
+
 - Solvers are more robust if used with multiple right-hand sides and one lane starts with the exact solution.
 
 - Added a function to write nested matrices as SVG objects: `writeSVGMatrix(...)`
diff --git a/dune/istl/solver.hh b/dune/istl/solver.hh
index c4b6ad7c06dc5ee9c79c66266d678381071c25dd..f4d97cb145e25a92a9a966cde14f8706d3ffe661 100644
--- a/dune/istl/solver.hh
+++ b/dune/istl/solver.hh
@@ -225,7 +225,7 @@ namespace Dune
        <li> 2 : print line for each iteration </li>
        </ul>
      */
-    IterativeSolver (LinearOperator<X,Y>& op, Preconditioner<X,Y>& prec, scalar_real_type reduction, int maxit, int verbose) :
+    IterativeSolver (const LinearOperator<X,Y>& op, Preconditioner<X,Y>& prec, scalar_real_type reduction, int maxit, int verbose) :
       _op(stackobject_to_shared_ptr(op)),
       _prec(stackobject_to_shared_ptr(prec)),
       _sp(new SeqScalarProduct<X>),
@@ -257,7 +257,7 @@ namespace Dune
         <li> 2 : print line for each iteration </li>
         </ul>
      */
-    IterativeSolver (LinearOperator<X,Y>& op, ScalarProduct<X>& sp, Preconditioner<X,Y>& prec,
+    IterativeSolver (const LinearOperator<X,Y>& op, const ScalarProduct<X>& sp, Preconditioner<X,Y>& prec,
       scalar_real_type reduction, int maxit, int verbose) :
       _op(stackobject_to_shared_ptr(op)),
       _prec(stackobject_to_shared_ptr(prec)),
@@ -285,7 +285,7 @@ namespace Dune
 
        See \ref ISTL_Factory for the ParameterTree layout and examples.
      */
-    IterativeSolver (std::shared_ptr<LinearOperator<X,Y> > op, std::shared_ptr<Preconditioner<X,X> > prec, const ParameterTree& configuration) :
+    IterativeSolver (std::shared_ptr<const LinearOperator<X,Y> > op, std::shared_ptr<Preconditioner<X,X> > prec, const ParameterTree& configuration) :
       IterativeSolver(op,std::make_shared<SeqScalarProduct<X>>(),prec,
         configuration.get<real_type>("reduction"),
         configuration.get<int>("maxit"),
@@ -308,7 +308,7 @@ namespace Dune
 
        See \ref ISTL_Factory for the ParameterTree layout and examples.
      */
-    IterativeSolver (std::shared_ptr<LinearOperator<X,Y> > op, std::shared_ptr<ScalarProduct<X> > sp, std::shared_ptr<Preconditioner<X,X> > prec, const ParameterTree& configuration) :
+    IterativeSolver (std::shared_ptr<const LinearOperator<X,Y> > op, std::shared_ptr<const ScalarProduct<X> > sp, std::shared_ptr<Preconditioner<X,X> > prec, const ParameterTree& configuration) :
       IterativeSolver(op,sp,prec,
         configuration.get<scalar_real_type>("reduction"),
         configuration.get<int>("maxit"),
@@ -335,8 +335,8 @@ namespace Dune
         <li> 2 : print line for each iteration </li>
         </ul>
      */
-    IterativeSolver (std::shared_ptr<LinearOperator<X,Y>> op,
-                     std::shared_ptr<ScalarProduct<X>> sp,
+    IterativeSolver (std::shared_ptr<const LinearOperator<X,Y>> op,
+                     std::shared_ptr<const ScalarProduct<X>> sp,
                      std::shared_ptr<Preconditioner<X,Y>> prec,
                      scalar_real_type reduction, int maxit, int verbose) :
       _op(op),
@@ -499,9 +499,9 @@ namespace Dune
     };
 
   protected:
-    std::shared_ptr<LinearOperator<X,Y>> _op;
+    std::shared_ptr<const LinearOperator<X,Y>> _op;
     std::shared_ptr<Preconditioner<X,Y>> _prec;
-    std::shared_ptr<ScalarProduct<X>> _sp;
+    std::shared_ptr<const ScalarProduct<X>> _sp;
     scalar_real_type _reduction;
     int _maxit;
     int _verbose;
diff --git a/dune/istl/solvers.hh b/dune/istl/solvers.hh
index 011e3eaa41fb1ad4b630707941209208e161b15d..24ae04feb9f992daf4798546c3e658be2de72dd4 100644
--- a/dune/istl/solvers.hh
+++ b/dune/istl/solvers.hh
@@ -212,12 +212,12 @@ namespace Dune {
 
     /*!
       \brief Constructor to initialize a CG solver.
-      \copydetails IterativeSolver::IterativeSolver(LinearOperator<X,Y>&, Preconditioner<X,Y>&, real_type, int, int)
+      \copydetails IterativeSolver::IterativeSolver(const LinearOperator<X,Y>&, Preconditioner<X,Y>&, real_type, int, int)
       \param condition_estimate Whether to calculate an estimate of the condition number.
                                 The estimate is given in the InverseOperatorResult returned by apply().
                                 This is only supported for float and double field types.
     */
-    CGSolver (LinearOperator<X,X>& op, Preconditioner<X,X>& prec,
+    CGSolver (const LinearOperator<X,X>& op, Preconditioner<X,X>& prec,
       scalar_real_type reduction, int maxit, int verbose, bool condition_estimate) : IterativeSolver<X,X>(op, prec, reduction, maxit, verbose),
       condition_estimate_(condition_estimate)
     {
@@ -229,12 +229,12 @@ namespace Dune {
 
     /*!
       \brief Constructor to initialize a CG solver.
-      \copydetails IterativeSolver::IterativeSolver(LinearOperator<X,Y>&, ScalarProduct<X>&, Preconditioner<X,Y>&, real_type, int, int)
+      \copydetails IterativeSolver::IterativeSolver(const LinearOperator<X,Y>&, const ScalarProduct<X>&, Preconditioner<X,Y>&, real_type, int, int)
       \param condition_estimate Whether to calculate an estimate of the condition number.
                                 The estimate is given in the InverseOperatorResult returned by apply().
                                 This is only supported for float and double field types.
     */
-    CGSolver (LinearOperator<X,X>& op, ScalarProduct<X>& sp, Preconditioner<X,X>& prec,
+    CGSolver (const LinearOperator<X,X>& op, const ScalarProduct<X>& sp, Preconditioner<X,X>& prec,
       scalar_real_type reduction, int maxit, int verbose, bool condition_estimate) : IterativeSolver<X,X>(op, sp, prec, reduction, maxit, verbose),
       condition_estimate_(condition_estimate)
     {
@@ -246,12 +246,12 @@ namespace Dune {
 
     /*!
       \brief Constructor to initialize a CG solver.
-      \copydetails IterativeSolver::IterativeSolver(std::shared_ptr<LinearOperator<X,Y>>, std::shared_ptr<ScalarProduct<X>>, std::shared_ptr<Preconditioner<X,Y>>, real_type, int, int)
+      \copydetails IterativeSolver::IterativeSolver(std::shared_ptr<const LinearOperator<X,Y>>, std::shared_ptr<ScalarProduct<X>>, std::shared_ptr<Preconditioner<X,Y>>, real_type, int, int)
       \param condition_estimate Whether to calculate an estimate of the condition number.
                                 The estimate is given in the InverseOperatorResult returned by apply().
                                 This is only supported for float and double field types.
     */
-    CGSolver (std::shared_ptr<LinearOperator<X,X>> op, std::shared_ptr<ScalarProduct<X>> sp,
+    CGSolver (std::shared_ptr<const LinearOperator<X,X>> op, std::shared_ptr<ScalarProduct<X>> sp,
               std::shared_ptr<Preconditioner<X,X>> prec,
               scalar_real_type reduction, int maxit, int verbose, bool condition_estimate)
       : IterativeSolver<X,X>(op, sp, prec, reduction, maxit, verbose),
@@ -842,10 +842,10 @@ namespace Dune {
     /*!
        \brief Set up RestartedGMResSolver solver.
 
-       \copydoc LoopSolver::LoopSolver(L&,P&,double,int,int)
+       \copydoc LoopSolver::LoopSolver(const L&,P&,double,int,int)
        \param restart number of GMRes cycles before restart
      */
-    RestartedGMResSolver (LinearOperator<X,Y>& op, Preconditioner<X,Y>& prec, scalar_real_type reduction, int restart, int maxit, int verbose) :
+    RestartedGMResSolver (const LinearOperator<X,Y>& op, Preconditioner<X,Y>& prec, scalar_real_type reduction, int restart, int maxit, int verbose) :
       IterativeSolver<X,Y>::IterativeSolver(op,prec,reduction,maxit,verbose),
       _restart(restart)
     {}
@@ -853,10 +853,10 @@ namespace Dune {
     /*!
        \brief Set up RestartedGMResSolver solver.
 
-       \copydoc LoopSolver::LoopSolver(L&,S&,P&,double,int,int)
+       \copydoc LoopSolver::LoopSolver(const L&, const S&,P&,double,int,int)
        \param restart number of GMRes cycles before restart
      */
-    RestartedGMResSolver (LinearOperator<X,Y>& op, ScalarProduct<X>& sp, Preconditioner<X,Y>& prec, scalar_real_type reduction, int restart, int maxit, int verbose) :
+    RestartedGMResSolver (const LinearOperator<X,Y>& op, const ScalarProduct<X>& sp, Preconditioner<X,Y>& prec, scalar_real_type reduction, int restart, int maxit, int verbose) :
       IterativeSolver<X,Y>::IterativeSolver(op,sp,prec,reduction,maxit,verbose),
       _restart(restart)
     {}
@@ -864,7 +864,7 @@ namespace Dune {
     /*!
        \brief Constructor.
 
-       \copydoc IterativeSolver::IterativeSolver(L&,S&,P&,const ParameterTree&)
+       \copydoc IterativeSolver::IterativeSolver(const L&, const S&,P&,const ParameterTree&)
 
        Additional parameter:
        ParameterTree Key | Meaning
@@ -873,12 +873,12 @@ namespace Dune {
 
        See \ref ISTL_Factory for the ParameterTree layout and examples.
      */
-    RestartedGMResSolver (std::shared_ptr<LinearOperator<X,Y> > op, std::shared_ptr<Preconditioner<X,X> > prec, const ParameterTree& configuration) :
+    RestartedGMResSolver (std::shared_ptr<const LinearOperator<X,Y> > op, std::shared_ptr<Preconditioner<X,X> > prec, const ParameterTree& configuration) :
       IterativeSolver<X,Y>::IterativeSolver(op,prec,configuration),
       _restart(configuration.get<int>("restart"))
     {}
 
-    RestartedGMResSolver (std::shared_ptr<LinearOperator<X,Y> > op, std::shared_ptr<ScalarProduct<X> > sp, std::shared_ptr<Preconditioner<X,X> > prec, const ParameterTree& configuration) :
+    RestartedGMResSolver (std::shared_ptr<const LinearOperator<X,Y> > op, std::shared_ptr<const ScalarProduct<X> > sp, std::shared_ptr<Preconditioner<X,X> > prec, const ParameterTree& configuration) :
       IterativeSolver<X,Y>::IterativeSolver(op,sp,prec,configuration),
       _restart(configuration.get<int>("restart"))
     {}
@@ -886,11 +886,11 @@ namespace Dune {
     /*!
       \brief Set up RestartedGMResSolver solver.
 
-      \copydoc LoopSolver::LoopSolver(std::shared_ptr<L>,std::shared_ptr<S>,std::shared_ptr<P>,double,int,int)
+      \copydoc LoopSolver::LoopSolver(std::shared_ptr<const L>,std::shared_ptr<const S>,std::shared_ptr<P>,double,int,int)
        \param restart number of GMRes cycles before restart
      */
-    RestartedGMResSolver (std::shared_ptr<LinearOperator<X,Y>> op,
-                          std::shared_ptr<ScalarProduct<X>> sp,
+    RestartedGMResSolver (std::shared_ptr<const LinearOperator<X,Y>> op,
+                          std::shared_ptr<const ScalarProduct<X>> sp,
                           std::shared_ptr<Preconditioner<X,Y>> prec,
                           scalar_real_type reduction, int restart, int maxit, int verbose) :
       IterativeSolver<X,Y>::IterativeSolver(op,sp,prec,reduction,maxit,verbose),
@@ -1323,10 +1323,10 @@ private:
     /*!
        \brief Set up nonlinear preconditioned conjugate gradient solver.
 
-       \copydoc LoopSolver::LoopSolver(L&,P&,double,int,int)
+       \copydoc LoopSolver::LoopSolver(const L&,P&,double,int,int)
        \param restart number of GMRes cycles before restart
      */
-    GeneralizedPCGSolver (LinearOperator<X,X>& op, Preconditioner<X,X>& prec, scalar_real_type reduction, int maxit, int verbose, int restart = 10) :
+    GeneralizedPCGSolver (const LinearOperator<X,X>& op, Preconditioner<X,X>& prec, scalar_real_type reduction, int maxit, int verbose, int restart = 10) :
       IterativeSolver<X,X>::IterativeSolver(op,prec,reduction,maxit,verbose),
       _restart(restart)
     {}
@@ -1334,11 +1334,11 @@ private:
     /*!
        \brief Set up nonlinear preconditioned conjugate gradient solver.
 
-       \copydoc LoopSolver::LoopSolver(L&,S&,P&,double,int,int)
+       \copydoc LoopSolver::LoopSolver(const L&, const S&,P&,double,int,int)
        \param restart When to restart the construction of
        the Krylov search space.
      */
-    GeneralizedPCGSolver (LinearOperator<X,X>& op, ScalarProduct<X>& sp, Preconditioner<X,X>& prec, scalar_real_type reduction, int maxit, int verbose, int restart = 10) :
+    GeneralizedPCGSolver (const LinearOperator<X,X>& op, const ScalarProduct<X>& sp, Preconditioner<X,X>& prec, scalar_real_type reduction, int maxit, int verbose, int restart = 10) :
       IterativeSolver<X,X>::IterativeSolver(op,sp,prec,reduction,maxit,verbose),
       _restart(restart)
     {}
@@ -1347,7 +1347,7 @@ private:
      /*!
        \brief Constructor.
 
-       \copydoc IterativeSolver::IterativeSolver(L&,S&,P&,const ParameterTree&)
+       \copydoc IterativeSolver::IterativeSolver(const L&,const S&,P&,const ParameterTree&)
 
        Additional parameter:
        ParameterTree Key | Meaning
@@ -1356,24 +1356,24 @@ private:
 
        See \ref ISTL_Factory for the ParameterTree layout and examples.
      */
-    GeneralizedPCGSolver (std::shared_ptr<LinearOperator<X,X> > op, std::shared_ptr<Preconditioner<X,X> > prec, const ParameterTree& configuration) :
+    GeneralizedPCGSolver (std::shared_ptr<const LinearOperator<X,X> > op, std::shared_ptr<Preconditioner<X,X> > prec, const ParameterTree& configuration) :
       IterativeSolver<X,X>::IterativeSolver(op,prec,configuration),
       _restart(configuration.get<int>("restart"))
     {}
 
-    GeneralizedPCGSolver (std::shared_ptr<LinearOperator<X,X> > op, std::shared_ptr<ScalarProduct<X> > sp, std::shared_ptr<Preconditioner<X,X> > prec, const ParameterTree& configuration) :
+    GeneralizedPCGSolver (std::shared_ptr<const LinearOperator<X,X> > op, std::shared_ptr<const ScalarProduct<X> > sp, std::shared_ptr<Preconditioner<X,X> > prec, const ParameterTree& configuration) :
       IterativeSolver<X,X>::IterativeSolver(op,sp,prec,configuration),
       _restart(configuration.get<int>("restart"))
     {}
     /*!
       \brief Set up nonlinear preconditioned conjugate gradient solver.
 
-      \copydoc LoopSolver::LoopSolver(std::shared_ptr<L>,std::shared_ptr<S>,std::shared_ptr<P>,double,int,int)
+      \copydoc LoopSolver::LoopSolver(std::shared_ptr<const L>,std::shared_ptr<const S>,std::shared_ptr<P>,double,int,int)
       \param restart When to restart the construction of
       the Krylov search space.
     */
-    GeneralizedPCGSolver (std::shared_ptr<LinearOperator<X,X>> op,
-                          std::shared_ptr<ScalarProduct<X>> sp,
+    GeneralizedPCGSolver (std::shared_ptr<const LinearOperator<X,X>> op,
+                          std::shared_ptr<const ScalarProduct<X>> sp,
                           std::shared_ptr<Preconditioner<X,X>> prec,
                           scalar_real_type reduction, int maxit, int verbose,
                           int restart = 10) :
@@ -1511,31 +1511,31 @@ private:
     using IterativeSolver<X,X>::apply;
     /*!
       \brief Constructor to initialize a RestartedFCG solver.
-      \copydetails IterativeSolver::IterativeSolver(LinearOperator<X,Y>&, Preconditioner<X,Y>&, real_type, int, int, int)
+      \copydetails IterativeSolver::IterativeSolver(const LinearOperator<X,Y>&, Preconditioner<X,Y>&, real_type, int, int, int)
       \param mmax is the maximal number of previous vectors which are orthogonalized against the new search direction.
     */
-    RestartedFCGSolver (LinearOperator<X,X>& op, Preconditioner<X,X>& prec,
+    RestartedFCGSolver (const LinearOperator<X,X>& op, Preconditioner<X,X>& prec,
                         scalar_real_type reduction, int maxit, int verbose, int mmax = 10) : IterativeSolver<X,X>(op, prec, reduction, maxit, verbose), _mmax(mmax)
     {
     }
 
     /*!
       \brief Constructor to initialize a RestartedFCG solver.
-      \copydetails IterativeSolver::IterativeSolver(LinearOperator<X,Y>&, ScalarProduct<X>&, Preconditioner<X,Y>&, real_type, int, int,int)
+      \copydetails IterativeSolver::IterativeSolver(const LinearOperator<X,Y>&, const ScalarProduct<X>&, Preconditioner<X,Y>&, real_type, int, int,int)
       \param mmax is the maximal number of previous vectors which are orthogonalized against the new search direction.
     */
-    RestartedFCGSolver (LinearOperator<X,X>& op, ScalarProduct<X>& sp, Preconditioner<X,X>& prec,
+    RestartedFCGSolver (const LinearOperator<X,X>& op, const ScalarProduct<X>& sp, Preconditioner<X,X>& prec,
                         scalar_real_type reduction, int maxit, int verbose, int mmax = 10) : IterativeSolver<X,X>(op, sp, prec, reduction, maxit, verbose), _mmax(mmax)
     {
     }
 
     /*!
       \brief Constructor to initialize a RestartedFCG solver.
-      \copydetails IterativeSolver::IterativeSolver(std::shared_ptr<LinearOperator<X,Y>>, std::shared_ptr<ScalarProduct<X>>, std::shared_ptr<Preconditioner<X,Y>>, real_type, int, int,int)
+      \copydetails IterativeSolver::IterativeSolver(std::shared_ptr<const LinearOperator<X,Y>>, std::shared_ptr<const ScalarProduct<X>>, std::shared_ptr<Preconditioner<X,Y>>, real_type, int, int,int)
       \param mmax is the maximal number of previous vectors which are orthogonalized against the new search direction.
     */
-    RestartedFCGSolver (std::shared_ptr<LinearOperator<X,X>> op,
-                        std::shared_ptr<ScalarProduct<X>> sp,
+    RestartedFCGSolver (std::shared_ptr<const LinearOperator<X,X>> op,
+                        std::shared_ptr<const ScalarProduct<X>> sp,
                         std::shared_ptr<Preconditioner<X,X>> prec,
                         scalar_real_type reduction, int maxit, int verbose,
                         int mmax = 10)
@@ -1545,7 +1545,7 @@ private:
     /*!
        \brief Constructor.
 
-       \copydoc IterativeSolver::IterativeSolver(L&,S&,P&,const ParameterTree&)
+       \copydoc IterativeSolver::IterativeSolver(const L&, const S&,P&,const ParameterTree&)
 
        Additional parameter:
        ParameterTree Key | Meaning
@@ -1554,14 +1554,14 @@ private:
 
        See \ref ISTL_Factory for the ParameterTree layout and examples.
      */
-    RestartedFCGSolver (std::shared_ptr<LinearOperator<X,X>> op,
+    RestartedFCGSolver (std::shared_ptr<const LinearOperator<X,X>> op,
                         std::shared_ptr<Preconditioner<X,X>> prec,
                         const ParameterTree& config)
       : IterativeSolver<X,X>(op, prec, config), _mmax(config.get("mmax", 10))
     {}
 
-    RestartedFCGSolver (std::shared_ptr<LinearOperator<X,X>> op,
-                        std::shared_ptr<ScalarProduct<X>> sp,
+    RestartedFCGSolver (std::shared_ptr<const LinearOperator<X,X>> op,
+                        std::shared_ptr<const ScalarProduct<X>> sp,
                         std::shared_ptr<Preconditioner<X,X>> prec,
                         const ParameterTree& config)
       : IterativeSolver<X,X>(op, sp, prec, config), _mmax(config.get("mmax", 10))
diff --git a/dune/istl/test/solvertest.cc b/dune/istl/test/solvertest.cc
index 54a2c9d859262187bda4c7208736603b41e08f4d..69cfb5dd1daebfff7329b88fd8bbdd4e54bcea8a 100644
--- a/dune/istl/test/solvertest.cc
+++ b/dune/istl/test/solvertest.cc
@@ -68,7 +68,7 @@ int main(int argc, char** argv)
   typedef Dune::MatrixAdapter<BCRSMat,BVector,BVector> Operator;
 
   BCRSMat mat;
-  Operator fop(mat);
+  const Operator fop(mat);
   BVector b(N*N), x(N*N);
 
   setupLaplacian(mat,N);