Loading...
Searching...
No Matches
poweriteration.hh
Go to the documentation of this file.
// 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
blocklevel.hh
Helper functions for determining the vector/matrix block level.
solvers.hh
Implementations of the inverse operator interface.
io.hh
Some generic functions for pretty printing vectors and matrices.
istlexception.hh
operators.hh
Define general, extensible interface for operators. The available implementation wraps a matrix.
solvertype.hh
Templates characterizing the type of a solver.
solvercategory.hh
Dune::printvector
void printvector(std::ostream &s, const V &v, std::string title, std::string rowtext, int columns=1, int width=10, int precision=2)
Print an ISTL vector.
Definition io.hh:89
Dune
Definition allocator.hh:11
Dune::blockLevel
constexpr std::size_t blockLevel()
Determine the block level of a possibly nested vector/matrix type.
Definition blocklevel.hh:176
Dune::BCRSMatrix
A sparse block matrix with compressed row storage.
Definition bcrsmatrix.hh:467
Dune::BCRSMatrix::size_type
A::size_type size_type
The type for the index access and the size.
Definition bcrsmatrix.hh:498
Dune::BlockVector
A vector of blocks with memory management.
Definition bvector.hh:392
Dune::BlockVector::field_type
typename Imp::BlockTraits< B >::field_type field_type
export the type representing the field
Definition bvector.hh:398
Dune::PowerIteration_Algorithms::itMatrix_
std::unique_ptr< BCRSMatrix > itMatrix_
Definition poweriteration.hh:1034
Dune::PowerIteration_Algorithms::PowerIteration_Algorithms
PowerIteration_Algorithms(const PowerIteration_Algorithms &)=delete
Dune::PowerIteration_Algorithms::ScalingOperator
Impl::ScalingLinearOperator< BlockVector > ScalingOperator
Definition poweriteration.hh:181
Dune::PowerIteration_Algorithms::blank_
const std::string blank_
Definition poweriteration.hh:1043
Dune::PowerIteration_Algorithms::OperatorSum
Impl::LinearOperatorSum< MatrixOperator, ScalingOperator > OperatorSum
Definition poweriteration.hh:182
Dune::PowerIteration_Algorithms::MatrixOperator
Dune::MatrixAdapter< BCRSMatrix, BlockVector, BlockVector > MatrixOperator
Definition poweriteration.hh:180
Dune::PowerIteration_Algorithms::applyInverseIteration
void applyInverseIteration(const Real &epsilon, ISTLLinearSolver &solver, BlockVector &x, Real &lambda) const
Perform the inverse iteration algorithm to compute an approximation lambda of the least dominant (i....
Definition poweriteration.hh:355
Dune::PowerIteration_Algorithms::applyTLIMEIteration
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
Perform the "two-level iterative method for eigenvalue calculations (TLIME)" iteration algorith...
Definition poweriteration.hh:691
Dune::PowerIteration_Algorithms::getIterationOperator
IterationOperator & getIterationOperator()
Return the iteration operator (m_ - mu_*I).
Definition poweriteration.hh:925
Dune::PowerIteration_Algorithms::itOperator_
OperatorSum itOperator_
Definition poweriteration.hh:1030
Dune::PowerIteration_Algorithms::m_
const BCRSMatrix & m_
Definition poweriteration.hh:1018
Dune::PowerIteration_Algorithms::PowerIteration_Algorithms
PowerIteration_Algorithms(const BCRSMatrix &m, const unsigned int nIterationsMax=1000, const unsigned int verbosity_level=0)
Construct from required parameters.
Definition poweriteration.hh:206
Dune::PowerIteration_Algorithms::nIterationsMax_
const unsigned int nIterationsMax_
Definition poweriteration.hh:1019
Dune::PowerIteration_Algorithms::applyPowerIteration
void applyPowerIteration(const Real &epsilon, BlockVector &x, Real &lambda) const
Perform the power iteration algorithm to compute an approximation lambda of the dominant (i....
Definition poweriteration.hh:260
Dune::PowerIteration_Algorithms::IterationOperator
OperatorSum IterationOperator
Type of iteration operator (m_ - mu_*I).
Definition poweriteration.hh:189
Dune::PowerIteration_Algorithms::applyRayleighQuotientIteration
void applyRayleighQuotientIteration(const Real &epsilon, ISTLLinearSolver &solver, BlockVector &x, Real &lambda) const
Perform the Rayleigh quotient iteration algorithm to compute an approximation lambda of an eigenvalue...
Definition poweriteration.hh:533
Dune::PowerIteration_Algorithms::applyInverseIteration
void applyInverseIteration(const Real &gamma, const Real &epsilon, ISTLLinearSolver &solver, BlockVector &x, Real &lambda) const
Perform the inverse iteration with shift algorithm to compute an approximation lambda of the eigenval...
Definition poweriteration.hh:394
Dune::PowerIteration_Algorithms::scalingOperator_
const ScalingOperator scalingOperator_
Definition poweriteration.hh:1029
Dune::PowerIteration_Algorithms::updateShiftMu
void updateShiftMu(const Real &mu, ISTLLinearSolver &solver) const
Update shift mu_, i.e. update iteration operator/matrix (m_ - mu_*I).
Definition poweriteration.hh:983
Dune::PowerIteration_Algorithms::operator=
PowerIteration_Algorithms & operator=(const PowerIteration_Algorithms &)=delete
Dune::PowerIteration_Algorithms::getIterationCount
unsigned int getIterationCount() const
Return the number of iterations in last application of an algorithm.
Definition poweriteration.hh:959
Dune::PowerIteration_Algorithms::matrixOperator_
const MatrixOperator matrixOperator_
Definition poweriteration.hh:1028
Dune::PowerIteration_Algorithms::getIterationMatrix
const BCRSMatrix & getIterationMatrix() const
Return the iteration matrix (m_ - mu_*I), provided on demand when needed (e.g. for direct solvers or ...
Definition poweriteration.hh:945
Dune::PowerIteration_Algorithms::nIterations_
unsigned int nIterations_
Definition poweriteration.hh:1039
Dune::PowerIteration_Algorithms::verbosity_level_
const unsigned int verbosity_level_
Definition poweriteration.hh:1022
Dune::PowerIteration_Algorithms::title_
const std::string title_
Definition poweriteration.hh:1042
Dune::PowerIteration_Algorithms::Real
BlockVector::field_type Real
Type of underlying field.
Definition poweriteration.hh:186
Dune::PowerIteration_Algorithms::mu_
Real mu_
Definition poweriteration.hh:1025
Dune::ISTLError
derive error class from the base class in common
Definition istlexception.hh:19
Dune::LinearOperator< X, X >::category
virtual SolverCategory::Category category() const=0
Dune::MatrixAdapter
Adapter to turn a matrix into a linear operator.
Definition operators.hh:135
Dune::InverseOperatorResult
Statistics about the application of an inverse operator.
Definition solver.hh:50
Dune::SolverHelper::setMatrix
static void setMatrix(ISTLLinearSolver &solver, const BCRSMatrix &matrix)
Definition solver.hh:526
Dune::SolverCategory::Category
Category
Definition solvercategory.hh:23
Dune::SolverCategory::sequential
@ sequential
Category for sequential solvers.
Definition solvercategory.hh:25