arpackpp.hh

// SPDX-FileCopyrightText: Copyright (C) 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_ARPACKPP_HH
#define DUNE_ISTL_EIGENVALUE_ARPACKPP_HH
 
#if HAVE_ARPACKPP || defined DOXYGEN
 
#include <cmath>  // provides std::abs, std::pow, std::sqrt
 
#include <iostream>  // provides std::cout, std::endl
#include <string>    // provides std::string
 
#include <dune/common/fvector.hh>     // provides Dune::FieldVector
#include <dune/common/exceptions.hh>  // provides DUNE_THROW(...)
 
#include <dune/istl/blocklevel.hh>     // provides Dune::blockLevel
#include <dune/istl/bvector.hh>        // provides Dune::BlockVector
#include <dune/istl/istlexception.hh>  // provides Dune::ISTLError
#include <dune/istl/io.hh>             // provides Dune::printvector(...)
 
#ifdef Status
#undef Status        // prevent preprocessor from damaging the ARPACK++
                     // code when "X11/Xlib.h" is included (the latter
                     // defines Status as "#define Status int" and
                     // ARPACK++ provides a class with a method called
                     // Status)
#endif
#include "arssym.h"  // provides ARSymStdEig
 
namespace Dune
{
 
  namespace Impl {
    template <class BCRSMatrix>
    class ArPackPlusPlus_BCRSMatrixWrapper
    {
    public:
      typedef typename BCRSMatrix::field_type Real;
 
    public:
      ArPackPlusPlus_BCRSMatrixWrapper (const BCRSMatrix& A)
        : A_(A),
          m_(A_.M() * mBlock), n_(A_.N() * nBlock)
      {
        // assert that BCRSMatrix type has blocklevel 2
        static_assert
          (blockLevel<BCRSMatrix>() == 2,
            "Only BCRSMatrices with blocklevel 2 are supported.");
 
        // allocate memory for auxiliary block vector objects
        // which are compatible to matrix rows / columns
        domainBlockVector.resize(A_.N());
        rangeBlockVector.resize(A_.M());
      }
 
      inline void multMv (Real* v, Real* w)
      {
        // get vector v as an object of appropriate type
        arrayToDomainBlockVector(v,domainBlockVector);
 
        // perform matrix-vector product
        A_.mv(domainBlockVector,rangeBlockVector);
 
        // get vector w from object of appropriate type
        rangeBlockVectorToArray(rangeBlockVector,w);
      };
 
      inline void multMtMv (Real* v, Real* w)
      {
        // get vector v as an object of appropriate type
        arrayToDomainBlockVector(v,domainBlockVector);
 
        // perform matrix-vector product
        A_.mv(domainBlockVector,rangeBlockVector);
        A_.mtv(rangeBlockVector,domainBlockVector);
 
        // get vector w from object of appropriate type
        domainBlockVectorToArray(domainBlockVector,w);
      };
 
      inline void multMMtv (Real* v, Real* w)
      {
        // get vector v as an object of appropriate type
        arrayToRangeBlockVector(v,rangeBlockVector);
 
        // perform matrix-vector product
        A_.mtv(rangeBlockVector,domainBlockVector);
        A_.mv(domainBlockVector,rangeBlockVector);
 
        // get vector w from object of appropriate type
        rangeBlockVectorToArray(rangeBlockVector,w);
      };
 
      inline int nrows () const { return m_; }
 
      inline int ncols () const { return n_; }
 
    protected:
      // Number of rows and columns in each block of the matrix
      constexpr static int mBlock = BCRSMatrix::block_type::rows;
      constexpr static int nBlock = BCRSMatrix::block_type::cols;
 
      // Type of vectors in the domain of the linear map associated with
      // the matrix, i.e. block vectors compatible to matrix rows
      constexpr static int dbvBlockSize = nBlock;
      typedef Dune::FieldVector<Real,dbvBlockSize> DomainBlockVectorBlock;
      typedef Dune::BlockVector<DomainBlockVectorBlock> DomainBlockVector;
 
      // Type of vectors in the range of the linear map associated with
      // the matrix, i.e. block vectors compatible to matrix columns
      constexpr static int rbvBlockSize = mBlock;
      typedef Dune::FieldVector<Real,rbvBlockSize> RangeBlockVectorBlock;
      typedef Dune::BlockVector<RangeBlockVectorBlock> RangeBlockVector;
 
      // Types for vector index access
      typedef typename DomainBlockVector::size_type dbv_size_type;
      typedef typename RangeBlockVector::size_type rbv_size_type;
      typedef typename DomainBlockVectorBlock::size_type dbvb_size_type;
      typedef typename RangeBlockVectorBlock::size_type rbvb_size_type;
 
      // Get vector v from a block vector object which is compatible to
      // matrix rows
      static inline void
      domainBlockVectorToArray (const DomainBlockVector& dbv, Real* v)
      {
        for (dbv_size_type block = 0; block < dbv.N(); ++block)
          for (dbvb_size_type iBlock = 0; iBlock < dbvBlockSize; ++iBlock)
            v[block*dbvBlockSize + iBlock] = dbv[block][iBlock];
      }
 
      // Get vector v from a block vector object which is compatible to
      // matrix columns
      static inline void
      rangeBlockVectorToArray (const RangeBlockVector& rbv, Real* v)
      {
        for (rbv_size_type block = 0; block < rbv.N(); ++block)
          for (rbvb_size_type iBlock = 0; iBlock < rbvBlockSize; ++iBlock)
            v[block*rbvBlockSize + iBlock] = rbv[block][iBlock];
      }
 
    public:
      static inline void arrayToDomainBlockVector (const Real* v,
        DomainBlockVector& dbv)
      {
        for (dbv_size_type block = 0; block < dbv.N(); ++block)
          for (dbvb_size_type iBlock = 0; iBlock < dbvBlockSize; ++iBlock)
            dbv[block][iBlock] = v[block*dbvBlockSize + iBlock];
      }
 
      static inline void arrayToRangeBlockVector (const Real* v,
        RangeBlockVector& rbv)
      {
        for (rbv_size_type block = 0; block < rbv.N(); ++block)
          for (rbvb_size_type iBlock = 0; iBlock < rbvBlockSize; ++iBlock)
            rbv[block][iBlock] = v[block*rbvBlockSize + iBlock];
      }
 
    protected:
      // The DUNE-ISTL BCRSMatrix
      const BCRSMatrix& A_;
 
      // Number of rows and columns in the matrix
      const int m_, n_;
 
      // Auxiliary block vector objects which are
      // compatible to matrix rows / columns
      mutable DomainBlockVector domainBlockVector;
      mutable RangeBlockVector rangeBlockVector;
    };
  } // end namespace Impl
 
  template <typename BCRSMatrix, typename BlockVector>
  class ArPackPlusPlus_Algorithms
  {
  public:
    typedef typename BlockVector::field_type Real;
 
  public:
    ArPackPlusPlus_Algorithms (const BCRSMatrix& m,
                               const unsigned int nIterationsMax = 100000,
                               const unsigned int verbosity_level = 0)
      : m_(m), nIterationsMax_(nIterationsMax),
        verbosity_level_(verbosity_level),
        nIterations_(0),
        title_("    ArPackPlusPlus_Algorithms: "),
        blank_(title_.length(),' ')
    {}
 
    inline void computeSymMaxMagnitude (const Real& epsilon,
                                        BlockVector& x, Real& lambda) const
    {
      // print verbosity information
      if (verbosity_level_ > 0)
        std::cout << title_ << "Computing an approximation of "
                  << "the dominant eigenvalue of a matrix which "
                  << "is assumed to be symmetric." << std::endl;
 
      // use type ArPackPlusPlus_BCRSMatrixWrapper to store matrix information
      // and to perform the product A*v (LU decomposition is not used)
      typedef Impl::ArPackPlusPlus_BCRSMatrixWrapper<BCRSMatrix> WrappedMatrix;
      WrappedMatrix A(m_);
 
      // get number of rows and columns in A
      const int nrows = A.nrows();
      const int ncols = A.ncols();
 
      // assert that A is square
      if (nrows != ncols)
        DUNE_THROW(Dune::ISTLError,"Matrix is not square ("
                   << nrows << "x" << ncols << ").");
 
      // allocate memory for variables, set parameters
      const int nev = 1;                     // Number of eigenvalues to compute
      int ncv = std::min(20, nrows);         // Number of Arnoldi vectors generated at each iteration (0 == auto)
      const Real tol = epsilon;              // Stopping tolerance (relative accuracy of Ritz values) (0 == machine precision)
      const int maxit = nIterationsMax_*nev; // Maximum number of Arnoldi update iterations allowed   (0 == 100*nev)
      Real* ev = new Real[nev];              // Computed eigenvalues of A
      const bool ivec = true;                // Flag deciding if eigenvectors shall be determined
      int nconv;                             // Number of converged eigenvalues
 
      // define what we need: eigenvalues with largest magnitude
      char which[] = "LM";
      ARSymStdEig<Real,WrappedMatrix>
        dprob(nrows, nev, &A, &WrappedMatrix::multMv, which, ncv, tol, maxit);
 
      // set ARPACK verbosity mode if requested
      if (verbosity_level_ > 3) dprob.Trace();
 
      // find eigenvalues and eigenvectors of A, obtain the eigenvalues
      nconv = dprob.Eigenvalues(ev,ivec);
 
      // obtain approximated dominant eigenvalue of A
      lambda = ev[nev-1];
 
      // obtain associated approximated eigenvector of A
      Real* x_raw = dprob.RawEigenvector(nev-1);
      WrappedMatrix::arrayToDomainBlockVector(x_raw,x);
 
      // obtain number of Arnoldi update iterations actually taken
      nIterations_ = dprob.GetIter();
 
      // compute residual norm
      BlockVector r(x);
      Real* Ax_raw = new Real[nrows];
      A.multMv(x_raw,Ax_raw);
      WrappedMatrix::arrayToDomainBlockVector(Ax_raw,r);
      r.axpy(-lambda,x);
      const Real r_norm = r.two_norm();
 
      // print verbosity information
      if (verbosity_level_ > 0)
      {
        if (verbosity_level_ > 1)
        {
          // print some information about the problem
          std::cout << blank_ << "Obtained eigenvalues of A by solving "
                    << "A*x = λ*x using the ARPACK++ class ARSym"
                    << "StdEig:" << std::endl;
          std::cout << blank_ << "       converged eigenvalues of A: "
                    << nconv << " / " << nev << std::endl;
          std::cout << blank_ << "         dominant eigenvalue of A: "
                    << lambda << 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");
        }
      }
 
      // free dynamically allocated memory
      delete[] Ax_raw;
      delete[] ev;
    }
 
    inline void computeSymMinMagnitude (const Real& epsilon,
                                        BlockVector& x, Real& lambda) const
    {
      // print verbosity information
      if (verbosity_level_ > 0)
        std::cout << title_ << "Computing an approximation of the "
                  << "least dominant eigenvalue of a matrix which "
                  << "is assumed to be symmetric." << std::endl;
 
      // use type ArPackPlusPlus_BCRSMatrixWrapper to store matrix information
      // and to perform the product A*v (LU decomposition is not used)
      typedef Impl::ArPackPlusPlus_BCRSMatrixWrapper<BCRSMatrix> WrappedMatrix;
      WrappedMatrix A(m_);
 
      // get number of rows and columns in A
      const int nrows = A.nrows();
      const int ncols = A.ncols();
 
      // assert that A is square
      if (nrows != ncols)
        DUNE_THROW(Dune::ISTLError,"Matrix is not square ("
                   << nrows << "x" << ncols << ").");
 
      // allocate memory for variables, set parameters
      const int nev = 1;                     // Number of eigenvalues to compute
      int ncv = std::min(20, nrows);         // Number of Arnoldi vectors generated at each iteration (0 == auto)
      const Real tol = epsilon;              // Stopping tolerance (relative accuracy of Ritz values) (0 == machine precision)
      const int maxit = nIterationsMax_*nev; // Maximum number of Arnoldi update iterations allowed   (0 == 100*nev)
      Real* ev = new Real[nev];              // Computed eigenvalues of A
      const bool ivec = true;                // Flag deciding if eigenvectors shall be determined
      int nconv;                             // Number of converged eigenvalues
 
      // define what we need: eigenvalues with smallest magnitude
      char which[] = "SM";
      ARSymStdEig<Real,WrappedMatrix>
        dprob(nrows, nev, &A, &WrappedMatrix::multMv, which, ncv, tol, maxit);
 
      // set ARPACK verbosity mode if requested
      if (verbosity_level_ > 3) dprob.Trace();
 
      // find eigenvalues and eigenvectors of A, obtain the eigenvalues
      nconv = dprob.Eigenvalues(ev,ivec);
 
      // obtain approximated least dominant eigenvalue of A
      lambda = ev[nev-1];
 
      // obtain associated approximated eigenvector of A
      Real* x_raw = dprob.RawEigenvector(nev-1);
      WrappedMatrix::arrayToDomainBlockVector(x_raw,x);
 
      // obtain number of Arnoldi update iterations actually taken
      nIterations_ = dprob.GetIter();
 
      // compute residual norm
      BlockVector r(x);
      Real* Ax_raw = new Real[nrows];
      A.multMv(x_raw,Ax_raw);
      WrappedMatrix::arrayToDomainBlockVector(Ax_raw,r);
      r.axpy(-lambda,x);
      const Real r_norm = r.two_norm();
 
      // print verbosity information
      if (verbosity_level_ > 0)
      {
        if (verbosity_level_ > 1)
        {
          // print some information about the problem
          std::cout << blank_ << "Obtained eigenvalues of A by solving "
                    << "A*x = λ*x using the ARPACK++ class ARSym"
                    << "StdEig:" << std::endl;
          std::cout << blank_ << "       converged eigenvalues of A: "
                    << nconv << " / " << nev << std::endl;
          std::cout << blank_ << "   least dominant eigenvalue of A: "
                    << lambda << 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");
        }
      }
 
      // free dynamically allocated memory
      delete[] Ax_raw;
      delete[] ev;
    }
 
    inline void computeSymCond2 (const Real& epsilon, Real& cond_2) const
    {
      // print verbosity information
      if (verbosity_level_ > 0)
        std::cout << title_ << "Computing an approximation of the "
                  << "spectral condition number of a matrix which "
                  << "is assumed to be symmetric." << std::endl;
 
      // use type ArPackPlusPlus_BCRSMatrixWrapper to store matrix information
      // and to perform the product A*v (LU decomposition is not used)
      typedef Impl::ArPackPlusPlus_BCRSMatrixWrapper<BCRSMatrix> WrappedMatrix;
      WrappedMatrix A(m_);
 
      // get number of rows and columns in A
      const int nrows = A.nrows();
      const int ncols = A.ncols();
 
      // assert that A is square
      if (nrows != ncols)
        DUNE_THROW(Dune::ISTLError,"Matrix is not square ("
                   << nrows << "x" << ncols << ").");
 
      // allocate memory for variables, set parameters
      const int nev = 2;                     // Number of eigenvalues to compute
      int ncv = std::min(20, nrows);         // Number of Arnoldi vectors generated at each iteration (0 == auto)
      const Real tol = epsilon;              // Stopping tolerance (relative accuracy of Ritz values) (0 == machine precision)
      const int maxit = nIterationsMax_*nev; // Maximum number of Arnoldi update iterations allowed   (0 == 100*nev)
      Real* ev = new Real[nev];              // Computed eigenvalues of A
      const bool ivec = true;                // Flag deciding if eigenvectors shall be determined
      int nconv;                             // Number of converged eigenvalues
 
      // define what we need: eigenvalues from both ends of the spectrum
      char which[] = "BE";
      ARSymStdEig<Real,WrappedMatrix>
        dprob(nrows, nev, &A, &WrappedMatrix::multMv, which, ncv, tol, maxit);
 
      // set ARPACK verbosity mode if requested
      if (verbosity_level_ > 3) dprob.Trace();
 
      // find eigenvalues and eigenvectors of A, obtain the eigenvalues
      nconv = dprob.Eigenvalues(ev,ivec);
 
      // obtain approximated dominant and least dominant eigenvalue of A
      const Real& lambda_max = ev[nev-1];
      const Real& lambda_min = ev[0];
 
      // obtain associated approximated eigenvectors of A
      Real* x_max_raw = dprob.RawEigenvector(nev-1);
      Real* x_min_raw = dprob.RawEigenvector(0);
 
      // obtain approximated spectral condition number of A
      cond_2 = std::abs(lambda_max / lambda_min);
 
      // obtain number of Arnoldi update iterations actually taken
      nIterations_ = dprob.GetIter();
 
      // compute each residual norm
      Real* Ax_max_raw = new Real[nrows];
      Real* Ax_min_raw = new Real[nrows];
      A.multMv(x_max_raw,Ax_max_raw);
      A.multMv(x_min_raw,Ax_min_raw);
      Real r_max_norm = 0.0;
      Real r_min_norm = 0.0;
      for (int i = 0; i < nrows; ++i)
      {
        r_max_norm += std::pow(Ax_max_raw[i] - lambda_max * x_max_raw[i],2);
        r_min_norm += std::pow(Ax_min_raw[i] - lambda_min * x_min_raw[i],2);
      }
      r_max_norm = std::sqrt(r_max_norm);
      r_min_norm = std::sqrt(r_min_norm);
 
      // print verbosity information
      if (verbosity_level_ > 0)
      {
        if (verbosity_level_ > 1)
        {
          // print some information about the problem
          std::cout << blank_ << "Obtained eigenvalues of A by solving "
                    << "A*x = λ*x using the ARPACK++ class ARSym"
                    << "StdEig:" << std::endl;
          std::cout << blank_ << "       converged eigenvalues of A: "
                    << nconv << " / " << nev << std::endl;
          std::cout << blank_ << "         dominant eigenvalue of A: "
                    << lambda_max << std::endl;
          std::cout << blank_ << "   least dominant eigenvalue of A: "
                    << lambda_min << std::endl;
          std::cout << blank_ << "   spectral condition number of A: "
                    << cond_2 << std::endl;
        }
        std::cout << blank_ << "Result ("
                  << "#iterations = " << nIterations_ << ", "
                  << "║residual║_2 = {" << r_max_norm << ","
                  << r_min_norm << "}, " << "λ = {"
                  << lambda_max << "," << lambda_min
                  << "}): cond_2 = " << cond_2 << std::endl;
      }
 
      // free dynamically allocated memory
      delete[] Ax_min_raw;
      delete[] Ax_max_raw;
      delete[] ev;
    }
 
    inline void computeNonSymMax (const Real& epsilon,
                                  BlockVector& x, Real& sigma) const
    {
      // print verbosity information
      if (verbosity_level_ > 0)
        std::cout << title_ << "Computing an approximation of the "
                  << "largest singular value of a matrix which "
                  << "is assumed to be nonsymmetric." << std::endl;
 
      // use type ArPackPlusPlus_BCRSMatrixWrapper to store matrix information
      // and to perform the product A^T*A*v (LU decomposition is not used)
      typedef Impl::ArPackPlusPlus_BCRSMatrixWrapper<BCRSMatrix> WrappedMatrix;
      WrappedMatrix A(m_);
 
      // get number of rows and columns in A
      const int nrows = A.nrows();
      const int ncols = A.ncols();
 
      // assert that A has more rows than columns (extend code later to the opposite case!)
      if (nrows < ncols)
        DUNE_THROW(Dune::ISTLError,"Matrix has less rows than "
                   << "columns (" << nrows << "x" << ncols << ")."
                   << " This case is not implemented, yet.");
 
      // allocate memory for variables, set parameters
      const int nev = 1;                     // Number of eigenvalues to compute
      int ncv = std::min(20, nrows);         // Number of Arnoldi vectors generated at each iteration (0 == auto)
      const Real tol = epsilon;              // Stopping tolerance (relative accuracy of Ritz values) (0 == machine precision)
      const int maxit = nIterationsMax_*nev; // Maximum number of Arnoldi update iterations allowed   (0 == 100*nev)
      Real* ev = new Real[nev];              // Computed eigenvalues of A^T*A
      const bool ivec = true;                // Flag deciding if eigenvectors shall be determined
      int nconv;                             // Number of converged eigenvalues
 
      // define what we need: eigenvalues with largest algebraic value
      char which[] = "LA";
      ARSymStdEig<Real,WrappedMatrix>
        dprob(ncols, nev, &A, &WrappedMatrix::multMtMv, which, ncv, tol, maxit);
 
      // set ARPACK verbosity mode if requested
      if (verbosity_level_ > 3) dprob.Trace();
 
      // find eigenvalues and eigenvectors of A^T*A, obtain the eigenvalues
      nconv = dprob.Eigenvalues(ev,ivec);
 
      // obtain approximated largest eigenvalue of A^T*A
      const Real& lambda = ev[nev-1];
 
      // obtain associated approximated eigenvector of A^T*A
      Real* x_raw = dprob.RawEigenvector(nev-1);
      WrappedMatrix::arrayToDomainBlockVector(x_raw,x);
 
      // obtain number of Arnoldi update iterations actually taken
      nIterations_ = dprob.GetIter();
 
      // compute residual norm
      BlockVector r(x);
      Real* AtAx_raw = new Real[ncols];
      A.multMtMv(x_raw,AtAx_raw);
      WrappedMatrix::arrayToDomainBlockVector(AtAx_raw,r);
      r.axpy(-lambda,x);
      const Real r_norm = r.two_norm();
 
      // calculate largest singular value of A (note that
      // x is right-singular / left-singular vector of A)
      sigma = std::sqrt(lambda);
 
      // print verbosity information
      if (verbosity_level_ > 0)
      {
        if (verbosity_level_ > 1)
        {
          // print some information about the problem
          std::cout << blank_ << "Obtained singular values of A by sol"
                    << "ving (A^T*A)*x = σ²*x using the ARPACK++ "
                    << "class ARSymStdEig:" << std::endl;
          std::cout << blank_ << "   converged eigenvalues of A^T*A: "
                    << nconv << " / " << nev << std::endl;
          std::cout << blank_ << "      largest eigenvalue of A^T*A: "
                    << lambda << std::endl;
          std::cout << blank_ << "   => largest singular value of A: "
                    << sigma << std::endl;
        }
        std::cout << blank_ << "Result ("
                  << "#iterations = " << nIterations_ << ", "
                  << "║residual║_2 = " << r_norm << "): "
                  << "σ = " << sigma << std::endl;
        if (verbosity_level_ > 2)
        {
          // print approximated right-singular / left-singular vector
          // via DUNE-ISTL I/O methods
          Dune::printvector(std::cout,x,blank_+"x",blank_+"row");
        }
      }
 
      // free dynamically allocated memory
      delete[] AtAx_raw;
      delete[] ev;
    }
 
    inline void computeNonSymMin (const Real& epsilon,
                                  BlockVector& x, Real& sigma) const
    {
      // print verbosity information
      if (verbosity_level_ > 0)
        std::cout << title_ << "Computing an approximation of the "
                  << "smallest singular value of a matrix which "
                  << "is assumed to be nonsymmetric." << std::endl;
 
      // use type ArPackPlusPlus_BCRSMatrixWrapper to store matrix information
      // and to perform the product A^T*A*v (LU decomposition is not used)
      typedef Impl::ArPackPlusPlus_BCRSMatrixWrapper<BCRSMatrix> WrappedMatrix;
      WrappedMatrix A(m_);
 
      // get number of rows and columns in A
      const int nrows = A.nrows();
      const int ncols = A.ncols();
 
      // assert that A has more rows than columns (extend code later to the opposite case!)
      if (nrows < ncols)
        DUNE_THROW(Dune::ISTLError,"Matrix has less rows than "
                   << "columns (" << nrows << "x" << ncols << ")."
                   << " This case is not implemented, yet.");
 
      // allocate memory for variables, set parameters
      const int nev = 1;                     // Number of eigenvalues to compute
      int ncv = std::min(20, nrows);         // Number of Arnoldi vectors generated at each iteration (0 == auto)
      const Real tol = epsilon;              // Stopping tolerance (relative accuracy of Ritz values) (0 == machine precision)
      const int maxit = nIterationsMax_*nev; // Maximum number of Arnoldi update iterations allowed   (0 == 100*nev)
      Real* ev = new Real[nev];              // Computed eigenvalues of A^T*A
      const bool ivec = true;                // Flag deciding if eigenvectors shall be determined
      int nconv;                             // Number of converged eigenvalues
 
      // define what we need: eigenvalues with smallest algebraic value
      char which[] = "SA";
      ARSymStdEig<Real,WrappedMatrix>
        dprob(ncols, nev, &A, &WrappedMatrix::multMtMv, which, ncv, tol, maxit);
 
      // set ARPACK verbosity mode if requested
      if (verbosity_level_ > 3) dprob.Trace();
 
      // find eigenvalues and eigenvectors of A^T*A, obtain the eigenvalues
      nconv = dprob.Eigenvalues(ev,ivec);
 
      // obtain approximated smallest eigenvalue of A^T*A
      const Real& lambda = ev[nev-1];
 
      // obtain associated approximated eigenvector of A^T*A
      Real* x_raw = dprob.RawEigenvector(nev-1);
      WrappedMatrix::arrayToDomainBlockVector(x_raw,x);
 
      // obtain number of Arnoldi update iterations actually taken
      nIterations_ = dprob.GetIter();
 
      // compute residual norm
      BlockVector r(x);
      Real* AtAx_raw = new Real[ncols];
      A.multMtMv(x_raw,AtAx_raw);
      WrappedMatrix::arrayToDomainBlockVector(AtAx_raw,r);
      r.axpy(-lambda,x);
      const Real r_norm = r.two_norm();
 
      // calculate smallest singular value of A (note that
      // x is right-singular / left-singular vector of A)
      sigma = std::sqrt(lambda);
 
      // print verbosity information
      if (verbosity_level_ > 0)
      {
        if (verbosity_level_ > 1)
        {
          // print some information about the problem
          std::cout << blank_ << "Obtained singular values of A by sol"
                    << "ving (A^T*A)*x = σ²*x using the ARPACK++ "
                    << "class ARSymStdEig:" << std::endl;
          std::cout << blank_ << "   converged eigenvalues of A^T*A: "
                    << nconv << " / " << nev << std::endl;
          std::cout << blank_ << "     smallest eigenvalue of A^T*A: "
                    << lambda << std::endl;
          std::cout << blank_ << "  => smallest singular value of A: "
                    << sigma << std::endl;
        }
        std::cout << blank_ << "Result ("
                  << "#iterations = " << nIterations_ << ", "
                  << "║residual║_2 = " << r_norm << "): "
                  << "σ = " << sigma << std::endl;
        if (verbosity_level_ > 2)
        {
          // print approximated right-singular / left-singular vector
          // via DUNE-ISTL I/O methods
          Dune::printvector(std::cout,x,blank_+"x",blank_+"row");
        }
      }
 
      // free dynamically allocated memory
      delete[] AtAx_raw;
      delete[] ev;
    }
 
    inline void computeNonSymCond2 (const Real& epsilon, Real& cond_2) const
    {
      // print verbosity information
      if (verbosity_level_ > 0)
        std::cout << title_ << "Computing an approximation of the "
                  << "spectral condition number of a matrix which "
                  << "is assumed to be nonsymmetric." << std::endl;
 
      // use type ArPackPlusPlus_BCRSMatrixWrapper to store matrix information
      // and to perform the product A^T*A*v (LU decomposition is not used)
      typedef Impl::ArPackPlusPlus_BCRSMatrixWrapper<BCRSMatrix> WrappedMatrix;
      WrappedMatrix A(m_);
 
      // get number of rows and columns in A
      const int nrows = A.nrows();
      const int ncols = A.ncols();
 
      // assert that A has more rows than columns (extend code later to the opposite case!)
      if (nrows < ncols)
        DUNE_THROW(Dune::ISTLError,"Matrix has less rows than "
                   << "columns (" << nrows << "x" << ncols << ")."
                   << " This case is not implemented, yet.");
 
      // allocate memory for variables, set parameters
      const int nev = 2;                     // Number of eigenvalues to compute
      int ncv = std::min(20, nrows);         // Number of Arnoldi vectors generated at each iteration (0 == auto)
      const Real tol = epsilon;              // Stopping tolerance (relative accuracy of Ritz values) (0 == machine precision)
      const int maxit = nIterationsMax_*nev; // Maximum number of Arnoldi update iterations allowed   (0 == 100*nev)
      Real* ev = new Real[nev];              // Computed eigenvalues of A^T*A
      const bool ivec = true;                // Flag deciding if eigenvectors shall be determined
      int nconv;                             // Number of converged eigenvalues
 
      // define what we need: eigenvalues from both ends of the spectrum
      char which[] = "BE";
      ARSymStdEig<Real,WrappedMatrix>
        dprob(ncols, nev, &A, &WrappedMatrix::multMtMv, which, ncv, tol, maxit);
 
      // set ARPACK verbosity mode if requested
      if (verbosity_level_ > 3) dprob.Trace();
 
      // find eigenvalues and eigenvectors of A^T*A, obtain the eigenvalues
      nconv = dprob.Eigenvalues(ev,ivec);
 
      // obtain approximated largest and smallest eigenvalue of A^T*A
      const Real& lambda_max = ev[nev-1];
      const Real& lambda_min = ev[0];
 
      // obtain associated approximated eigenvectors of A^T*A
      Real* x_max_raw = dprob.RawEigenvector(nev-1);
      Real* x_min_raw = dprob.RawEigenvector(0);
 
      // obtain number of Arnoldi update iterations actually taken
      nIterations_ = dprob.GetIter();
 
      // compute each residual norm
      Real* AtAx_max_raw = new Real[ncols];
      Real* AtAx_min_raw = new Real[ncols];
      A.multMtMv(x_max_raw,AtAx_max_raw);
      A.multMtMv(x_min_raw,AtAx_min_raw);
      Real r_max_norm = 0.0;
      Real r_min_norm = 0.0;
      for (int i = 0; i < ncols; ++i)
      {
        r_max_norm += std::pow(AtAx_max_raw[i] - lambda_max * x_max_raw[i],2);
        r_min_norm += std::pow(AtAx_min_raw[i] - lambda_min * x_min_raw[i],2);
      }
      r_max_norm = std::sqrt(r_max_norm);
      r_min_norm = std::sqrt(r_min_norm);
 
      // calculate largest and smallest singular value of A
      const Real sigma_max = std::sqrt(lambda_max);
      const Real sigma_min = std::sqrt(lambda_min);
 
      // obtain approximated spectral condition number of A
      cond_2 = sigma_max / sigma_min;
 
      // print verbosity information
      if (verbosity_level_ > 0)
      {
        if (verbosity_level_ > 1)
        {
          // print some information about the problem
          std::cout << blank_ << "Obtained singular values of A by sol"
                    << "ving (A^T*A)*x = σ²*x using the ARPACK++ "
                    << "class ARSymStdEig:" << std::endl;
          std::cout << blank_ << "   converged eigenvalues of A^T*A: "
                    << nconv << " / " << nev << std::endl;
          std::cout << blank_ << "      largest eigenvalue of A^T*A: "
                    << lambda_max << std::endl;
          std::cout << blank_ << "     smallest eigenvalue of A^T*A: "
                    << lambda_min << std::endl;
          std::cout << blank_ << "  =>  largest singular value of A: "
                    << sigma_max << std::endl;
          std::cout << blank_ << "  => smallest singular value of A: "
                    << sigma_min << std::endl;
        }
        std::cout << blank_ << "Result ("
                  << "#iterations = " << nIterations_ << ", "
                  << "║residual║_2 = {" << r_max_norm << ","
                  << r_min_norm << "}, " << "σ = {"
                  << sigma_max << "," << sigma_min
                  << "}): cond_2 = " << cond_2 << std::endl;
      }
 
      // free dynamically allocated memory
      delete[] AtAx_min_raw;
      delete[] AtAx_max_raw;
      delete[] ev;
    }
 
    inline unsigned int getIterationCount () const
    {
      if (nIterations_ == 0)
        DUNE_THROW(Dune::ISTLError,"No algorithm applied, yet.");
 
      return nIterations_;
    }
 
  protected:
    // parameters related to iterative eigenvalue algorithms
    const BCRSMatrix& m_;
    const unsigned int nIterationsMax_;
 
    // verbosity setting
    const unsigned int verbosity_level_;
 
    // 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  // HAVE_ARPACKPP
 
#endif  // DUNE_ISTL_EIGENVALUE_ARPACKPP_HH