poweriteration.hh

// SPDX-FileCopyrightText: Copyright © DUNE Project contributors, see file LICENSE.md in module root
// SPDX-License-Identifier: LicenseRef-GPL-2.0-only-with-DUNE-exception
// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
// vi: set et ts=4 sw=2 sts=2:
#ifndef DUNE_ISTL_EIGENVALUE_POWERITERATION_HH
#define DUNE_ISTL_EIGENVALUE_POWERITERATION_HH
 
#include <cstddef>  // provides std::size_t
#include <cmath>    // provides std::sqrt, std::abs
 
#include <type_traits>  // provides std::is_same
#include <iostream>     // provides std::cout, std::endl
#include <limits>       // provides std::numeric_limits
#include <ios>          // provides std::left, std::ios::left
#include <iomanip>      // provides std::setw, std::resetiosflags
#include <memory>       // provides std::unique_ptr
#include <string>       // provides std::string
 
#include <dune/common/exceptions.hh>  // provides DUNE_THROW(...)
 
#include <dune/istl/blocklevel.hh>      // provides Dune::blockLevel
#include <dune/istl/operators.hh>       // provides Dune::LinearOperator
#include <dune/istl/solvercategory.hh>  // provides Dune::SolverCategory::sequential
#include <dune/istl/solvertype.hh>      // provides Dune::IsDirectSolver
#include <dune/istl/operators.hh>       // provides Dune::MatrixAdapter
#include <dune/istl/istlexception.hh>   // provides Dune::ISTLError
#include <dune/istl/io.hh>              // provides Dune::printvector(...)
#include <dune/istl/solvers.hh>         // provides Dune::InverseOperatorResult
 
namespace Dune
{
 
  namespace Impl {
    template <class X, class Y = X>
    class ScalingLinearOperator : public Dune::LinearOperator<X,Y>
    {
    public:
      typedef X domain_type;
      typedef Y range_type;
      typedef typename X::field_type field_type;
 
      ScalingLinearOperator (field_type immutable_scaling,
        const field_type& mutable_scaling)
        : immutable_scaling_(immutable_scaling),
          mutable_scaling_(mutable_scaling)
      {}
 
      void apply (const X& x, Y& y) const override
      {
        y = x;
        y *= immutable_scaling_*mutable_scaling_;
      }
 
      void applyscaleadd (field_type alpha, const X& x, Y& y) const override
      {
        X temp(x);
        temp *= immutable_scaling_*mutable_scaling_;
        y.axpy(alpha,temp);
      }
 
      SolverCategory::Category category() const override
      {
        return SolverCategory::sequential;
      }
 
    protected:
      const field_type immutable_scaling_;
      const field_type& mutable_scaling_;
    };
 
 
    template <class OP1, class OP2>
    class LinearOperatorSum
      : public Dune::LinearOperator<typename OP1::domain_type,
                                    typename OP1::range_type>
    {
    public:
      typedef typename OP1::domain_type domain_type;
      typedef typename OP1::range_type range_type;
      typedef typename domain_type::field_type field_type;
 
      LinearOperatorSum (const OP1& op1, const OP2& op2)
        : op1_(op1), op2_(op2)
      {
        static_assert(std::is_same<typename OP2::domain_type,domain_type>::value,
          "Domain type of both operators doesn't match!");
        static_assert(std::is_same<typename OP2::range_type,range_type>::value,
          "Range type of both operators doesn't match!");
      }
 
      void apply (const domain_type& x, range_type& y) const override
      {
        op1_.apply(x,y);
        op2_.applyscaleadd(1.0,x,y);
      }
 
      void applyscaleadd (field_type alpha,
        const domain_type& x, range_type& y) const override
      {
        range_type temp(y);
        op1_.apply(x,temp);
        op2_.applyscaleadd(1.0,x,temp);
        y.axpy(alpha,temp);
      }
 
      SolverCategory::Category category() const override
      {
        return SolverCategory::sequential;
      }
 
    protected:
      const OP1& op1_;
      const OP2& op2_;
    };
  } // end namespace Impl
 
  template <typename BCRSMatrix, typename BlockVector>
  class PowerIteration_Algorithms
  {
  protected:
    // Type definitions for type of iteration operator (m_ - mu_*I)
    typedef typename Dune::MatrixAdapter<BCRSMatrix,BlockVector,BlockVector>
      MatrixOperator;
    typedef Impl::ScalingLinearOperator<BlockVector> ScalingOperator;
    typedef Impl::LinearOperatorSum<MatrixOperator,ScalingOperator> OperatorSum;
 
  public:
    typedef typename BlockVector::field_type Real;
 
    typedef OperatorSum IterationOperator;
 
  public:
    PowerIteration_Algorithms (const BCRSMatrix& m,
                               const unsigned int nIterationsMax = 1000,
                               const unsigned int verbosity_level = 0)
      : m_(m), nIterationsMax_(nIterationsMax),
        verbosity_level_(verbosity_level),
        mu_(0.0),
        matrixOperator_(m_),
        scalingOperator_(-1.0,mu_),
        itOperator_(matrixOperator_,scalingOperator_),
        nIterations_(0),
        title_("    PowerIteration_Algorithms: "),
        blank_(title_.length(),' ')
    {
      // assert that BCRSMatrix type has blocklevel 2
      static_assert
        (blockLevel<BCRSMatrix>() == 2,
         "Only BCRSMatrices with blocklevel 2 are supported.");
 
      // assert that BCRSMatrix type has square blocks
      static_assert
        (BCRSMatrix::block_type::rows == BCRSMatrix::block_type::cols,
         "Only BCRSMatrices with square blocks are supported.");
 
      // assert that m_ is square
      const int nrows = m_.M() * BCRSMatrix::block_type::rows;
      const int ncols = m_.N() * BCRSMatrix::block_type::cols;
      if (nrows != ncols)
        DUNE_THROW(Dune::ISTLError,"Matrix is not square ("
                   << nrows << "x" << ncols << ").");
    }
 
    PowerIteration_Algorithms (const PowerIteration_Algorithms&) = delete;
 
    PowerIteration_Algorithms&
      operator= (const PowerIteration_Algorithms&) = delete;
 
    inline void applyPowerIteration (const Real& epsilon,
                                     BlockVector& x, Real& lambda) const
    {
      // print verbosity information
      if (verbosity_level_ > 0)
        std::cout << title_
                  << "Performing power iteration approximating "
                  << "the dominant eigenvalue." << std::endl;
 
      // allocate memory for auxiliary variables
      BlockVector y(x);
      BlockVector temp(x);
 
      // perform power iteration
      x *= (1.0 / x.two_norm());
      m_.mv(x,y);
      Real r_norm = std::numeric_limits<Real>::max();
      nIterations_ = 0;
      while (r_norm > epsilon)
      {
        // update and check number of iterations
        if (++nIterations_ > nIterationsMax_)
          DUNE_THROW(Dune::ISTLError,"Power iteration did not converge "
                     << "in " << nIterationsMax_ << " iterations "
                     << "(║residual║_2 = " << r_norm << ", epsilon = "
                     << epsilon << ").");
 
        // do one iteration of the power iteration algorithm
        // (use that y = m_ * x)
        x = y;
        x *= (1.0 / y.two_norm());
 
        // get approximated eigenvalue lambda via the Rayleigh quotient
        m_.mv(x,y);
        lambda = x * y;
 
        // get norm of residual (use that y = m_ * x)
        temp = y;
        temp.axpy(-lambda,x);
        r_norm = temp.two_norm();
 
        // print verbosity information
        if (verbosity_level_ > 1)
          std::cout << blank_ << std::left
                    << "iteration " << std::setw(3) << nIterations_
                    << " (║residual║_2 = " << std::setw(11) << r_norm
                    << "): λ = " << lambda << std::endl
                    << std::resetiosflags(std::ios::left);
      }
 
      // print verbosity information
      if (verbosity_level_ > 0)
      {
        std::cout << blank_ << "Result ("
                  << "#iterations = " << nIterations_ << ", "
                  << "║residual║_2 = " << r_norm << "): "
                  << "λ = " << lambda << std::endl;
        if (verbosity_level_ > 2)
        {
          // print approximated eigenvector via DUNE-ISTL I/O methods
          Dune::printvector(std::cout,x,blank_+"x",blank_+"row");
        }
      }
    }
 
    template <typename ISTLLinearSolver,
              bool avoidLinSolverCrime = false>
    inline void applyInverseIteration (const Real& epsilon,
                                       ISTLLinearSolver& solver,
                                       BlockVector& x, Real& lambda) const
    {
      constexpr Real gamma = 0.0;
      applyInverseIteration(gamma,epsilon,solver,x,lambda);
    }
 
    template <typename ISTLLinearSolver,
              bool avoidLinSolverCrime = false>
    inline void applyInverseIteration (const Real& gamma,
                                       const Real& epsilon,
                                       ISTLLinearSolver& solver,
                                       BlockVector& x, Real& lambda) const
    {
      // print verbosity information
      if (verbosity_level_ > 0)
      {
        std::cout << title_;
        if (gamma == 0.0)
          std::cout << "Performing inverse iteration approximating "
                    << "the least dominant eigenvalue." << std::endl;
        else
          std::cout << "Performing inverse iteration with shift "
                    << "gamma = " << gamma << " approximating the "
                    << "eigenvalue closest to gamma." << std::endl;
      }
 
      // initialize iteration operator,
      // initialize iteration matrix when needed
      updateShiftMu(gamma,solver);
 
      // allocate memory for linear solver statistics
      Dune::InverseOperatorResult solver_statistics;
 
      // allocate memory for auxiliary variables
      BlockVector y(x);
      Real y_norm;
      BlockVector temp(x);
 
      // perform inverse iteration with shift
      x *= (1.0 / x.two_norm());
      Real r_norm = std::numeric_limits<Real>::max();
      nIterations_ = 0;
      while (r_norm > epsilon)
      {
        // update and check number of iterations
        if (++nIterations_ > nIterationsMax_)
          DUNE_THROW(Dune::ISTLError,"Inverse iteration "
                     << (gamma != 0.0 ? "with shift " : "") << "did not "
                     << "converge in " << nIterationsMax_ << " iterations "
                     << "(║residual║_2 = " << r_norm << ", epsilon = "
                     << epsilon << ").");
 
        // do one iteration of the inverse iteration with shift algorithm,
        // part 1: solve (m_ - gamma*I) * y = x for y
        // (protect x from being changed)
        temp = x;
        solver.apply(y,temp,solver_statistics);
 
        // get norm of y
        y_norm = y.two_norm();
 
        // compile time switch between accuracy and efficiency
        if (avoidLinSolverCrime)
        {
          // get approximated eigenvalue lambda via the Rayleigh quotient
          // (use that x_new = y / y_norm)
          m_.mv(y,temp);
          lambda = (y * temp) / (y_norm * y_norm);
 
          // get norm of residual
          // (use that x_new = y / y_norm, additionally use that temp = m_ * y)
          temp.axpy(-lambda,y);
          r_norm = temp.two_norm() / y_norm;
        }
        else
        {
          // get approximated eigenvalue lambda via the Rayleigh quotient
          // (use that x_new = y / y_norm and use that (m_ - gamma*I) * y = x)
          lambda = gamma + (y * x) / (y_norm * y_norm);
 
          // get norm of residual
          // (use that x_new = y / y_norm and use that (m_ - gamma*I) * y = x)
          temp = x; temp.axpy(gamma-lambda,y);
          r_norm = temp.two_norm() / y_norm;
        }
 
        // do one iteration of the inverse iteration with shift algorithm,
        // part 2: update x
        x = y;
        x *= (1.0 / y_norm);
 
        // print verbosity information
        if (verbosity_level_ > 1)
          std::cout << blank_ << std::left
                    << "iteration " << std::setw(3) << nIterations_
                    << " (║residual║_2 = " << std::setw(11) << r_norm
                    << "): λ = " << lambda << std::endl
                    << std::resetiosflags(std::ios::left);
      }
 
      // print verbosity information
      if (verbosity_level_ > 0)
      {
        std::cout << blank_ << "Result ("
                  << "#iterations = " << nIterations_ << ", "
                  << "║residual║_2 = " << r_norm << "): "
                  << "λ = " << lambda << std::endl;
        if (verbosity_level_ > 2)
        {
          // print approximated eigenvector via DUNE-ISTL I/O methods
          Dune::printvector(std::cout,x,blank_+"x",blank_+"row");
        }
      }
    }
 
    template <typename ISTLLinearSolver,
              bool avoidLinSolverCrime = false>
    inline void applyRayleighQuotientIteration (const Real& epsilon,
                                                ISTLLinearSolver& solver,
                                                BlockVector& x, Real& lambda) const
    {
      // print verbosity information
      if (verbosity_level_ > 0)
        std::cout << title_
                  << "Performing Rayleigh quotient iteration for "
                  << "estimated eigenvalue " << lambda << "." << std::endl;
 
      // allocate memory for linear solver statistics
      Dune::InverseOperatorResult solver_statistics;
 
      // allocate memory for auxiliary variables
      BlockVector y(x);
      Real y_norm;
      Real lambda_update;
      BlockVector temp(x);
 
      // perform Rayleigh quotient iteration
      x *= (1.0 / x.two_norm());
      Real r_norm = std::numeric_limits<Real>::max();
      nIterations_ = 0;
      while (r_norm > epsilon)
      {
        // update and check number of iterations
        if (++nIterations_ > nIterationsMax_)
          DUNE_THROW(Dune::ISTLError,"Rayleigh quotient iteration did not "
                     << "converge in " << nIterationsMax_ << " iterations "
                     << "(║residual║_2 = " << r_norm << ", epsilon = "
                     << epsilon << ").");
 
        // update iteration operator,
        // update iteration matrix when needed
        updateShiftMu(lambda,solver);
 
        // do one iteration of the Rayleigh quotient iteration algorithm,
        // part 1: solve (m_ - lambda*I) * y = x for y
        // (protect x from being changed)
        temp = x;
        solver.apply(y,temp,solver_statistics);
 
        // get norm of y
        y_norm = y.two_norm();
 
        // compile time switch between accuracy and efficiency
        if (avoidLinSolverCrime)
        {
          // get approximated eigenvalue lambda via the Rayleigh quotient
          // (use that x_new = y / y_norm)
          m_.mv(y,temp);
          lambda = (y * temp) / (y_norm * y_norm);
 
          // get norm of residual
          // (use that x_new = y / y_norm, additionally use that temp = m_ * y)
          temp.axpy(-lambda,y);
          r_norm = temp.two_norm() / y_norm;
        }
        else
        {
          // get approximated eigenvalue lambda via the Rayleigh quotient
          // (use that x_new = y / y_norm and use that (m_ - lambda_old*I) * y = x)
          lambda_update = (y * x) / (y_norm * y_norm);
          lambda += lambda_update;
 
          // get norm of residual
          // (use that x_new = y / y_norm and use that (m_ - lambda_old*I) * y = x)
          temp = x; temp.axpy(-lambda_update,y);
          r_norm = temp.two_norm() / y_norm;
        }
 
        // do one iteration of the Rayleigh quotient iteration algorithm,
        // part 2: update x
        x = y;
        x *= (1.0 / y_norm);
 
        // print verbosity information
        if (verbosity_level_ > 1)
          std::cout << blank_ << std::left
                    << "iteration " << std::setw(3) << nIterations_
                    << " (║residual║_2 = " << std::setw(11) << r_norm
                    << "): λ = " << lambda << std::endl
                    << std::resetiosflags(std::ios::left);
      }
 
      // print verbosity information
      if (verbosity_level_ > 0)
      {
        std::cout << blank_ << "Result ("
                  << "#iterations = " << nIterations_ << ", "
                  << "║residual║_2 = " << r_norm << "): "
                  << "λ = " << lambda << std::endl;
        if (verbosity_level_ > 2)
        {
          // print approximated eigenvector via DUNE-ISTL I/O methods
          Dune::printvector(std::cout,x,blank_+"x",blank_+"row");
        }
      }
    }
 
    template <typename ISTLLinearSolver,
              bool avoidLinSolverCrime = false>
    inline void applyTLIMEIteration (const Real& gamma, const Real& eta,
                                     const Real& epsilon,
                                     ISTLLinearSolver& solver,
                                     const Real& delta, const std::size_t& m,
                                     bool& extrnl,
                                     BlockVector& x, Real& lambda) const
    {
      // use same variable names as in [Szyld, 1988]
      BlockVector& x_s = x;
      Real& mu_s = lambda;
 
      // print verbosity information
      if (verbosity_level_ > 0)
        std::cout << title_
                  << "Performing TLIME iteration for "
                  << "estimated eigenvalue in the "
                  << "interval (" << gamma - eta << ","
                  << gamma + eta << ")." << std::endl;
 
      // allocate memory for linear solver statistics
      Dune::InverseOperatorResult solver_statistics;
 
      // allocate memory for auxiliary variables
      bool doRQI;
      Real mu;
      BlockVector y(x_s);
      Real omega;
      Real mu_s_old;
      Real mu_s_update;
      BlockVector temp(x_s);
      Real q_norm, r_norm;
 
      // perform TLIME iteration
      x_s *= (1.0 / x_s.two_norm());
      extrnl = true;
      doRQI = false;
      r_norm = std::numeric_limits<Real>::max();
      nIterations_ = 0;
      while (r_norm > epsilon)
      {
        // update and check number of iterations
        if (++nIterations_ > nIterationsMax_)
          DUNE_THROW(Dune::ISTLError,"TLIME iteration did not "
                     << "converge in " << nIterationsMax_
                     << " iterations (║residual║_2 = " << r_norm
                     << ", epsilon = " << epsilon << ").");
 
        // set shift for next iteration according to inverse iteration
        // with shift (II) resp. Rayleigh quotient iteration (RQI)
        if (doRQI)
          mu = mu_s;
        else
          mu = gamma;
 
        // update II/RQI iteration operator,
        // update II/RQI iteration matrix when needed
        updateShiftMu(mu,solver);
 
        // do one iteration of the II/RQI algorithm,
        // part 1: solve (m_ - mu*I) * y = x for y
        temp = x_s;
        solver.apply(y,temp,solver_statistics);
 
        // do one iteration of the II/RQI algorithm,
        // part 2: compute omega
        omega = (1.0 / y.two_norm());
 
        // backup the old Rayleigh quotient
        mu_s_old = mu_s;
 
        // compile time switch between accuracy and efficiency
        if (avoidLinSolverCrime)
        {
          // update the Rayleigh quotient mu_s, i.e. the approximated eigenvalue
          // (use that x_new = y * omega)
          m_.mv(y,temp);
          mu_s = (y * temp) * (omega * omega);
 
          // get norm of "the residual with respect to the shift used by II",
          // use normal representation of q
          // (use that x_new = y * omega, use that temp = m_ * y)
          temp.axpy(-gamma,y);
          q_norm = temp.two_norm() * omega;
 
          // get norm of "the residual with respect to the Rayleigh quotient"
          r_norm = q_norm*q_norm - (gamma-mu_s)*(gamma-mu_s);
          // prevent that truncation errors invalidate the norm
          // (we don't want to calculate sqrt of a negative number)
          if (r_norm >= 0)
          {
            // use relation between the norms of r and q for efficiency
            r_norm = std::sqrt(r_norm);
          }
          else
          {
            // use relation between r and q
            // (use that x_new = y * omega, use that temp = (m_ - gamma*I) * y = q / omega)
            temp.axpy(gamma-mu_s,y);
            r_norm = temp.two_norm() * omega;
          }
        }
        else
        {
          // update the Rayleigh quotient mu_s, i.e. the approximated eigenvalue
          if (!doRQI)
          {
            // (use that x_new = y * omega, additionally use that (m_ - gamma*I) * y = x_s)
            mu_s = gamma + (y * x_s) * (omega * omega);
          }
          else
          {
            // (use that x_new = y * omega, additionally use that (m_ - mu_s_old*I) * y = x_s)
            mu_s_update = (y * x_s) * (omega * omega);
            mu_s += mu_s_update;
          }
 
          // get norm of "the residual with respect to the shift used by II"
          if (!doRQI)
          {
            // use special representation of q in the II case
            // (use that x_new = y * omega, additionally use that (m_ - gamma*I) * y = x_s)
            q_norm = omega;
          }
          else
          {
            // use special representation of q in the RQI case
            // (use that x_new = y * omega, additionally use that (m_ - mu_s_old*I) * y = x_s)
            temp = x_s; temp.axpy(mu_s-gamma,y);
            q_norm = temp.two_norm() * omega;
          }
 
          // get norm of "the residual with respect to the Rayleigh quotient"
          // don't use efficient relation between the norms of r and q, as
          // this relation seems to yield a less accurate r_norm in the case
          // where linear solver crime is admitted
          if (!doRQI)
          {
            // (use that x_new = y * omega and use that (m_ - gamma*I) * y = x_s)
            temp = x_s; temp.axpy(gamma-lambda,y);
            r_norm = temp.two_norm() * omega;
          }
          else
          {
            // (use that x_new = y * omega and use that (m_ - mu_s_old*I) * y = x_s)
            temp = x_s; temp.axpy(-mu_s_update,y);
            r_norm = temp.two_norm() * omega;
          }
        }
 
        // do one iteration of the II/RQI algorithm,
        // part 3: update x
        x_s = y; x_s *= omega;
 
        // // for relative residual norm mode, scale with mu_s^{-1}
        // r_norm /= std::abs(mu_s);
 
        // print verbosity information
        if (verbosity_level_ > 1)
          std::cout << blank_ << "iteration "
                    << std::left << std::setw(3) << nIterations_
                    << " (" << (doRQI ? "RQI," : "II, ")
                    << " " << (doRQI ? "—>" : "  ") << " "
                    << "║r║_2 = " << std::setw(11) << r_norm
                    << ", " << (doRQI ? "  " : "—>") << " "
                    << "║q║_2 = " << std::setw(11) << q_norm
                    << "): λ = " << lambda << std::endl
                    << std::resetiosflags(std::ios::left);
 
        // check if the eigenvalue closest to gamma lies in J
        if (!doRQI && q_norm < eta)
        {
          // J is not free of eigenvalues
          extrnl = false;
 
          // by theory we know now that mu_s also lies in J
          assert(std::abs(mu_s-gamma) < eta);
 
          // switch to RQI
          doRQI = true;
        }
 
        // revert to II if J is not free of eigenvalues but
        // at some point mu_s falls back again outside J
        if (!extrnl && doRQI && std::abs(mu_s-gamma) >= eta)
          doRQI = false;
 
        // if eigenvalue closest to gamma does not lie in J use RQI
        // solely to accelerate the convergence to this eigenvalue
        // when II has become stationary
        if (extrnl && !doRQI)
        {
          // switch to RQI if the relative change of the Rayleigh
          // quotient indicates that II has become stationary
          if (nIterations_ >= m &&
              std::abs(mu_s - mu_s_old) / std::abs(mu_s) < delta)
            doRQI = true;
        }
      }
 
      // // compute final residual and lambda again (paranoia....)
      // m_.mv(x_s,temp);
      // mu_s = x_s * temp;
      // temp.axpy(-mu_s,x_s);
      // r_norm = temp.two_norm();
      // // r_norm /= std::abs(mu_s);
 
      // print verbosity information
      if (verbosity_level_ > 0)
      {
        if (extrnl)
          std::cout << blank_ << "Interval "
                    << "(" << gamma - eta << "," << gamma + eta
                    << ") is free of eigenvalues, approximating "
                    << "the closest eigenvalue." << std::endl;
        std::cout << blank_ << "Result ("
                  << "#iterations = " << nIterations_ << ", "
                  << "║residual║_2 = " << r_norm << "): "
                  << "λ = " << lambda << std::endl;
        if (verbosity_level_ > 2)
        {
          // print approximated eigenvector via DUNE-ISTL I/O methods
          Dune::printvector(std::cout,x,blank_+"x",blank_+"row");
        }
      }
    }
 
    inline IterationOperator& getIterationOperator ()
    {
      // return iteration operator
      return itOperator_;
    }
 
    inline const BCRSMatrix& getIterationMatrix () const
    {
      // create iteration matrix on demand
      if (!itMatrix_)
        itMatrix_ = std::make_unique<BCRSMatrix>(m_);
 
      // return iteration matrix
      return *itMatrix_;
    }
 
    inline unsigned int getIterationCount () const
    {
      if (nIterations_ == 0)
        DUNE_THROW(Dune::ISTLError,"No algorithm applied, yet.");
 
      return nIterations_;
    }
 
  protected:
    template <typename ISTLLinearSolver>
    inline void updateShiftMu (const Real& mu,
                               ISTLLinearSolver& solver) const
    {
      // do nothing if new shift equals the old one
      if (mu == mu_) return;
 
      // update shift mu_, i.e. update iteration operator
      mu_ = mu;
 
      // update iteration matrix when needed
      if (itMatrix_)
      {
        // iterate over entries in iteration matrix diagonal
        constexpr int rowBlockSize = BCRSMatrix::block_type::rows;
        constexpr int colBlockSize = BCRSMatrix::block_type::cols;
        for (typename BCRSMatrix::size_type i = 0;
             i < itMatrix_->M()*rowBlockSize; ++i)
        {
          // access m_[i,i] where i is the flat index of a row/column
          const Real& m_entry = m_
            [i/rowBlockSize][i/colBlockSize][i%rowBlockSize][i%colBlockSize];
          // access *itMatrix[i,i] where i is the flat index of a row/column
          Real& entry = (*itMatrix_)
            [i/rowBlockSize][i/colBlockSize][i%rowBlockSize][i%colBlockSize];
          // change current entry in iteration matrix diagonal
          entry = m_entry - mu_;
        }
        // notify linear solver about change of the iteration matrix object
        SolverHelper<ISTLLinearSolver,BCRSMatrix>::setMatrix
          (solver,*itMatrix_);
      }
    }
 
  protected:
    // parameters related to iterative eigenvalue algorithms
    const BCRSMatrix& m_;
    const unsigned int nIterationsMax_;
 
    // verbosity setting
    const unsigned int verbosity_level_;
 
    // shift mu_ used by iteration operator/matrix (m_ - mu_*I)
    mutable Real mu_;
 
    // iteration operator (m_ - mu_*I), passing shift mu_ by reference
    const MatrixOperator matrixOperator_;
    const ScalingOperator scalingOperator_;
    OperatorSum itOperator_;
 
    // iteration matrix (m_ - mu_*I), provided on demand when needed
    // (e.g. for preconditioning)
    mutable std::unique_ptr<BCRSMatrix> itMatrix_;
 
    // memory for storing temporary variables (mutable as they shall
    // just be effectless auxiliary variables of the const apply*(...)
    // methods)
    mutable unsigned int nIterations_;
 
    // constants for printing verbosity information
    const std::string title_;
    const std::string blank_;
  };
 
}  // namespace Dune
 
#endif  // DUNE_ISTL_EIGENVALUE_POWERITERATION_HH