3#ifndef DUNE_ISTL_EIGENVALUE_POWERITERATION_HH
4#define DUNE_ISTL_EIGENVALUE_POWERITERATION_HH
17#include <dune/common/exceptions.hh>
38 template <
class X,
class Y = X>
54 virtual void apply (
const X& x,
Y&
y)
const
81 template <
class OP1,
class OP2>
84 typename OP1::range_type>
96 static_assert(std::is_same<typename OP2::domain_type,domain_type>::value,
97 "Domain type of both operators doesn't match!");
98 static_assert(std::is_same<typename OP2::range_type,range_type>::value,
99 "Range type of both operators doesn't match!");
105 op2_.applyscaleadd(1.0,x,
y);
130 template <
typename ISTLLinearSolver,
typename BCRSMatrix>
148 template <
bool is_direct_solver,
typename Dummy =
void>
160 template <
typename Dummy>
209 template <
typename BCRSMatrix,
typename BlockVector>
251 title_(
" PowerIteration_Algorithms: "),
257 "Only BCRSMatrices with blocklevel 2 are supported.");
261 (BCRSMatrix::block_type::rows == BCRSMatrix::block_type::cols,
262 "Only BCRSMatrices with square blocks are supported.");
265 const int nrows =
m_.M() * BCRSMatrix::block_type::rows;
266 const int ncols =
m_.N() * BCRSMatrix::block_type::cols;
301 <<
"Performing power iteration approximating "
302 <<
"the dominant eigenvalue." << std::endl;
311 Real r_norm = std::numeric_limits<Real>::max();
319 <<
"(║residual║_2 = " <<
r_norm <<
", epsilon = "
325 x *= (1.0 /
y.two_norm());
338 std::cout <<
blank_ << std::left
340 <<
" (║residual║_2 = " << std::setw(11) <<
r_norm
341 <<
"): λ = " <<
lambda << std::endl
342 << std::resetiosflags(std::ios::left);
348 std::cout <<
blank_ <<
"Result ("
350 <<
"║residual║_2 = " <<
r_norm <<
"): "
351 <<
"λ = " <<
lambda << std::endl;
439 std::cout <<
"Performing inverse iteration approximating "
440 <<
"the least dominant eigenvalue." << std::endl;
442 std::cout <<
"Performing inverse iteration with shift "
443 <<
"gamma = " <<
gamma <<
" approximating the "
444 <<
"eigenvalue closest to gamma." << std::endl;
461 Real r_norm = std::numeric_limits<Real>::max();
468 << (
gamma != 0.0 ?
"with shift " :
"") <<
"did not "
470 <<
"(║residual║_2 = " <<
r_norm <<
", epsilon = "
514 std::cout <<
blank_ << std::left
516 <<
" (║residual║_2 = " << std::setw(11) <<
r_norm
517 <<
"): λ = " <<
lambda << std::endl
518 << std::resetiosflags(std::ios::left);
524 std::cout <<
blank_ <<
"Result ("
526 <<
"║residual║_2 = " <<
r_norm <<
"): "
527 <<
"λ = " <<
lambda << std::endl;
575 <<
"Performing Rayleigh quotient iteration for "
576 <<
"estimated eigenvalue " <<
lambda <<
"." << std::endl;
589 Real r_norm = std::numeric_limits<Real>::max();
597 <<
"(║residual║_2 = " <<
r_norm <<
", epsilon = "
646 std::cout <<
blank_ << std::left
648 <<
" (║residual║_2 = " << std::setw(11) <<
r_norm
649 <<
"): λ = " <<
lambda << std::endl
650 << std::resetiosflags(std::ios::left);
656 std::cout <<
blank_ <<
"Result ("
658 <<
"║residual║_2 = " <<
r_norm <<
"): "
659 <<
"λ = " <<
lambda << std::endl;
740 <<
"Performing TLIME iteration for "
741 <<
"estimated eigenvalue in the "
742 <<
"interval (" <<
gamma -
eta <<
","
743 <<
gamma +
eta <<
")." << std::endl;
759 x_s *= (1.0 /
x_s.two_norm());
762 r_norm = std::numeric_limits<Real>::max();
770 <<
" iterations (║residual║_2 = " <<
r_norm
771 <<
", epsilon = " <<
epsilon <<
").");
791 omega = (1.0 /
y.two_norm());
884 std::cout <<
blank_ <<
"iteration "
886 <<
" (" << (
doRQI ?
"RQI," :
"II, ")
887 <<
" " << (
doRQI ?
"—>" :
" ") <<
" "
888 <<
"║r║_2 = " << std::setw(11) <<
r_norm
889 <<
", " << (
doRQI ?
" " :
"—>") <<
" "
890 <<
"║q║_2 = " << std::setw(11) <<
q_norm
891 <<
"): λ = " <<
lambda << std::endl
892 << std::resetiosflags(std::ios::left);
936 std::cout <<
blank_ <<
"Interval "
938 <<
") is free of eigenvalues, approximating "
939 <<
"the closest eigenvalue." << std::endl;
940 std::cout <<
blank_ <<
"Result ("
942 <<
"║residual║_2 = " <<
r_norm <<
"): "
943 <<
"λ = " <<
lambda << std::endl;
1017 template <
typename ISTLLinearSolver>
1022 if (
mu ==
mu_)
return;
1031 constexpr int rowBlockSize = BCRSMatrix::block_type::rows;
1032 constexpr int colBlockSize = BCRSMatrix::block_type::cols;
1040 Real& entry = (*itMatrix_)
Some generic functions for pretty printing vectors and matrices.
Define general, extensible interface for operators. The available implementation wraps a matrix.
Implementations of the inverse operator interface.
Templates characterizing the type of a solver.
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:102
Definition basearray.hh:19
Statistics about compression achieved in implicit mode.
Definition bcrsmatrix.hh:81
A sparse block matrix with compressed row storage.
Definition bcrsmatrix.hh:412
A::size_type size_type
The type for the index access and the size.
Definition bcrsmatrix.hh:446
@ blocklevel
The number of blocklevels the matrix contains.
Definition bcrsmatrix.hh:454
FieldTraits< field_type >::real_type two_norm() const
two norm sqrt(sum over squared values of entries)
Definition bvector.hh:189
A vector of blocks with memory management.
Definition bvector.hh:309
B::field_type field_type
export the type representing the field
Definition bvector.hh:315
A linear operator scaling vectors by a scalar value. The scalar value can be changed as it is given i...
Definition poweriteration.hh:40
Y range_type
Definition poweriteration.hh:43
@ category
Definition poweriteration.hh:46
X::field_type field_type
Definition poweriteration.hh:44
virtual void apply(const X &x, Y &y) const
apply operator to x: The input vector is consistent and the output must also be consistent on the in...
Definition poweriteration.hh:54
const field_type & mutable_scaling_
Definition poweriteration.hh:69
const field_type immutable_scaling_
Definition poweriteration.hh:68
X domain_type
Definition poweriteration.hh:42
virtual void applyscaleadd(field_type alpha, const X &x, Y &y) const
apply operator to x, scale and add:
Definition poweriteration.hh:60
ScalingLinearOperator(field_type immutable_scaling, const field_type &mutable_scaling)
Definition poweriteration.hh:48
A linear operator representing the sum of two linear operators.
Definition poweriteration.hh:85
OP1::domain_type domain_type
Definition poweriteration.hh:87
LinearOperatorSum(const OP1 &op1, const OP2 &op2)
Definition poweriteration.hh:93
const OP1 & op1_
Definition poweriteration.hh:118
@ category
Definition poweriteration.hh:91
virtual void apply(const domain_type &x, range_type &y) const
apply operator to x: The input vector is consistent and the output must also be consistent on the in...
Definition poweriteration.hh:102
const OP2 & op2_
Definition poweriteration.hh:119
domain_type::field_type field_type
Definition poweriteration.hh:89
virtual void applyscaleadd(field_type alpha, const domain_type &x, range_type &y) const
Definition poweriteration.hh:108
OP1::range_type range_type
Definition poweriteration.hh:88
Helper class for notifying a DUNE-ISTL linear solver about a change of the iteration matrix object in...
Definition poweriteration.hh:132
static void setMatrix(ISTLLinearSolver &solver, const BCRSMatrix &matrix)
Definition poweriteration.hh:134
Implementation that works together with iterative ISTL solvers, e.g. Dune::CGSolver or Dune::BiCGSTAB...
Definition poweriteration.hh:150
static void setMatrix(ISTLLinearSolver &, const BCRSMatrix &)
Definition poweriteration.hh:151
static void setMatrix(ISTLLinearSolver &solver, const BCRSMatrix &matrix)
Definition poweriteration.hh:163
A class template for performing some iterative eigenvalue algorithms based on power iteration.
Definition poweriteration.hh:211
std::unique_ptr< BCRSMatrix > itMatrix_
Definition poweriteration.hh:1069
PowerIteration_Algorithms(const PowerIteration_Algorithms &)=delete
const std::string blank_
Definition poweriteration.hh:1078
Dune::MatrixAdapter< BCRSMatrix, BlockVector, BlockVector > MatrixOperator
Definition poweriteration.hh:215
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:390
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 algorit...
Definition poweriteration.hh:726
IterationOperator & getIterationOperator()
Return the iteration operator (m_ - mu_*I).
Definition poweriteration.hh:960
OperatorSum itOperator_
Definition poweriteration.hh:1065
LinearOperatorSum< MatrixOperator, ScalingOperator > OperatorSum
Definition poweriteration.hh:217
const BCRSMatrix & m_
Definition poweriteration.hh:1053
PowerIteration_Algorithms(const BCRSMatrix &m, const unsigned int nIterationsMax=1000, const unsigned int verbosity_level=0)
Construct from required parameters.
Definition poweriteration.hh:241
const unsigned int nIterationsMax_
Definition poweriteration.hh:1054
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:295
OperatorSum IterationOperator
Type of iteration operator (m_ - mu_*I)
Definition poweriteration.hh:224
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:568
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:429
const ScalingOperator scalingOperator_
Definition poweriteration.hh:1064
void updateShiftMu(const Real &mu, ISTLLinearSolver &solver) const
Update shift mu_, i.e. update iteration operator/matrix (m_ - mu_*I).
Definition poweriteration.hh:1018
PowerIteration_Algorithms & operator=(const PowerIteration_Algorithms &)=delete
unsigned int getIterationCount() const
Return the number of iterations in last application of an algorithm.
Definition poweriteration.hh:994
const MatrixOperator matrixOperator_
Definition poweriteration.hh:1063
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:980
unsigned int nIterations_
Definition poweriteration.hh:1074
const unsigned int verbosity_level_
Definition poweriteration.hh:1057
const std::string title_
Definition poweriteration.hh:1077
ScalingLinearOperator< BlockVector > ScalingOperator
Definition poweriteration.hh:216
BlockVector::field_type Real
Type of underlying field.
Definition poweriteration.hh:221
Real mu_
Definition poweriteration.hh:1060
derive error class from the base class in common
Definition istlexception.hh:16
A linear operator.
Definition operators.hh:62
Statistics about the application of an inverse operator.
Definition solver.hh:32
@ sequential
Category for sequential solvers.
Definition solvercategory.hh:21
Definition solvertype.hh:14