3#ifndef DUNE_ISTL_EIGENVALUE_POWERITERATION_HH
4#define DUNE_ISTL_EIGENVALUE_POWERITERATION_HH
20#include <dune/istl/solvercategory.hh>
23#include <dune/istl/istlexception.hh>
38 template <
class X,
class Y = X>
42 typedef X domain_type;
44 typedef typename X::field_type field_type;
49 const field_type& mutable_scaling)
50 : immutable_scaling_(immutable_scaling),
51 mutable_scaling_(mutable_scaling)
54 virtual void apply (
const X& x, Y& y)
const
57 y *= immutable_scaling_*mutable_scaling_;
60 virtual void applyscaleadd (field_type alpha,
const X& x, Y& y)
const
63 temp *= immutable_scaling_*mutable_scaling_;
68 const field_type immutable_scaling_;
69 const field_type& mutable_scaling_;
81 template <
class OP1,
class OP2>
84 typename OP1::range_type>
87 typedef typename OP1::domain_type domain_type;
88 typedef typename OP1::range_type range_type;
89 typedef typename domain_type::field_type field_type;
94 : op1_(op1), op2_(op2)
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!");
102 virtual void apply (
const domain_type& x, range_type& y)
const
105 op2_.applyscaleadd(1.0,x,y);
108 virtual void applyscaleadd (field_type alpha,
109 const domain_type& x, range_type& y)
const
113 op2_.applyscaleadd(1.0,x,temp);
130 template <
typename ISTLLinearSolver,
typename BCRSMatrix>
134 static void setMatrix (ISTLLinearSolver& solver,
137 static const bool is_direct_solver
138 = Dune::IsDirectSolver<ISTLLinearSolver>::value;
148 template <
bool is_direct_solver,
typename Dummy =
void>
151 static void setMatrix (ISTLLinearSolver&,
160 template <
typename Dummy>
163 static void setMatrix (ISTLLinearSolver& solver,
166 solver.setMatrix(matrix);
209 template <
typename BCRSMatrix,
typename BlockVector>
242 const unsigned int nIterationsMax = 1000,
243 const unsigned int verbosity_level = 0)
244 : m_(m), nIterationsMax_(nIterationsMax),
245 verbosity_level_(verbosity_level),
248 scalingOperator_(-1.0,mu_),
249 itOperator_(matrixOperator_,scalingOperator_),
251 title_(
" PowerIteration_Algorithms: "),
252 blank_(title_.length(),
' ')
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;
269 << nrows <<
"x" << ncols <<
").");
299 if (verbosity_level_ > 0)
301 <<
"Performing power iteration approximating "
302 <<
"the dominant eigenvalue." << std::endl;
311 Real r_norm = std::numeric_limits<Real>::max();
313 while (r_norm > epsilon)
316 if (++nIterations_ > nIterationsMax_)
318 <<
"in " << nIterationsMax_ <<
" iterations "
319 <<
"(║residual║_2 = " << r_norm <<
", epsilon = "
333 temp.
axpy(-lambda,x);
337 if (verbosity_level_ > 1)
338 std::cout << blank_ << std::left
339 <<
"iteration " << std::setw(3) << nIterations_
340 <<
" (║residual║_2 = " << std::setw(11) << r_norm
341 <<
"): λ = " << lambda << std::endl
342 << std::resetiosflags(std::ios::left);
346 if (verbosity_level_ > 0)
348 std::cout << blank_ <<
"Result ("
349 <<
"#iterations = " << nIterations_ <<
", "
350 <<
"║residual║_2 = " << r_norm <<
"): "
351 <<
"λ = " << lambda << std::endl;
352 if (verbosity_level_ > 2)
388 template <
typename ISTLLinearSolver,
389 bool avoidLinSolverCrime =
false>
391 ISTLLinearSolver& solver,
394 constexpr Real gamma = 0.0;
427 template <
typename ISTLLinearSolver,
428 bool avoidLinSolverCrime =
false>
431 ISTLLinearSolver& solver,
435 if (verbosity_level_ > 0)
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();
463 while (r_norm > epsilon)
466 if (++nIterations_ > nIterationsMax_)
468 << (gamma != 0.0 ?
"with shift " :
"") <<
"did not "
469 <<
"converge in " << nIterationsMax_ <<
" iterations "
470 <<
"(║residual║_2 = " << r_norm <<
", epsilon = "
477 solver.apply(y,temp,solver_statistics);
483 if (avoidLinSolverCrime)
488 lambda = (y * temp) / (y_norm * y_norm);
492 temp.
axpy(-lambda,y);
499 lambda = gamma + (y * x) / (y_norm * y_norm);
503 temp = x; temp.
axpy(gamma-lambda,y);
513 if (verbosity_level_ > 1)
514 std::cout << blank_ << std::left
515 <<
"iteration " << std::setw(3) << nIterations_
516 <<
" (║residual║_2 = " << std::setw(11) << r_norm
517 <<
"): λ = " << lambda << std::endl
518 << std::resetiosflags(std::ios::left);
522 if (verbosity_level_ > 0)
524 std::cout << blank_ <<
"Result ("
525 <<
"#iterations = " << nIterations_ <<
", "
526 <<
"║residual║_2 = " << r_norm <<
"): "
527 <<
"λ = " << lambda << std::endl;
528 if (verbosity_level_ > 2)
566 template <
typename ISTLLinearSolver,
567 bool avoidLinSolverCrime =
false>
569 ISTLLinearSolver& solver,
573 if (verbosity_level_ > 0)
575 <<
"Performing Rayleigh quotient iteration for "
576 <<
"estimated eigenvalue " << lambda <<
"." << std::endl;
589 Real r_norm = std::numeric_limits<Real>::max();
591 while (r_norm > epsilon)
594 if (++nIterations_ > nIterationsMax_)
596 <<
"converge in " << nIterationsMax_ <<
" iterations "
597 <<
"(║residual║_2 = " << r_norm <<
", epsilon = "
608 solver.apply(y,temp,solver_statistics);
614 if (avoidLinSolverCrime)
619 lambda = (y * temp) / (y_norm * y_norm);
623 temp.
axpy(-lambda,y);
630 lambda_update = (y * x) / (y_norm * y_norm);
631 lambda += lambda_update;
635 temp = x; temp.
axpy(-lambda_update,y);
645 if (verbosity_level_ > 1)
646 std::cout << blank_ << std::left
647 <<
"iteration " << std::setw(3) << nIterations_
648 <<
" (║residual║_2 = " << std::setw(11) << r_norm
649 <<
"): λ = " << lambda << std::endl
650 << std::resetiosflags(std::ios::left);
654 if (verbosity_level_ > 0)
656 std::cout << blank_ <<
"Result ("
657 <<
"#iterations = " << nIterations_ <<
", "
658 <<
"║residual║_2 = " << r_norm <<
"): "
659 <<
"λ = " << lambda << std::endl;
660 if (verbosity_level_ > 2)
724 template <
typename ISTLLinearSolver,
725 bool avoidLinSolverCrime =
false>
728 ISTLLinearSolver& solver,
729 const Real& delta,
const std::size_t& m,
738 if (verbosity_level_ > 0)
740 <<
"Performing TLIME iteration for "
741 <<
"estimated eigenvalue in the "
742 <<
"interval (" << gamma - eta <<
","
743 << gamma + eta <<
")." << std::endl;
762 r_norm = std::numeric_limits<Real>::max();
764 while (r_norm > epsilon)
767 if (++nIterations_ > nIterationsMax_)
769 <<
"converge in " << nIterationsMax_
770 <<
" iterations (║residual║_2 = " << r_norm
771 <<
", epsilon = " << epsilon <<
").");
787 solver.apply(y,temp,solver_statistics);
797 if (avoidLinSolverCrime)
802 mu_s = (y * temp) * (omega * omega);
811 r_norm = q_norm*q_norm - (gamma-mu_s)*(gamma-mu_s);
817 r_norm = std::sqrt(r_norm);
823 temp.
axpy(gamma-mu_s,y);
833 mu_s = gamma + (y * x_s) * (omega * omega);
838 mu_s_update = (y * x_s) * (omega * omega);
853 temp = x_s; temp.
axpy(mu_s-gamma,y);
864 temp = x_s; temp.
axpy(gamma-lambda,y);
870 temp = x_s; temp.
axpy(-mu_s_update,y);
877 x_s = y; x_s *= omega;
883 if (verbosity_level_ > 1)
884 std::cout << blank_ <<
"iteration "
885 << std::left << std::setw(3) << nIterations_
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);
895 if (!doRQI && q_norm < eta)
901 assert(std::abs(mu_s-gamma) < eta);
909 if (!extrnl && doRQI && std::abs(mu_s-gamma) >= eta)
915 if (extrnl && !doRQI)
919 if (nIterations_ >= m &&
920 std::abs(mu_s - mu_s_old) / std::abs(mu_s) < delta)
933 if (verbosity_level_ > 0)
936 std::cout << blank_ <<
"Interval "
937 <<
"(" << gamma - eta <<
"," << gamma + eta
938 <<
") is free of eigenvalues, approximating "
939 <<
"the closest eigenvalue." << std::endl;
940 std::cout << blank_ <<
"Result ("
941 <<
"#iterations = " << nIterations_ <<
", "
942 <<
"║residual║_2 = " << r_norm <<
"): "
943 <<
"λ = " << lambda << std::endl;
944 if (verbosity_level_ > 2)
984 itMatrix_ = std::unique_ptr<BCRSMatrix>(
new BCRSMatrix(m_));
996 if (nIterations_ == 0)
1017 template <
typename ISTLLinearSolver>
1019 ISTLLinearSolver& solver)
const
1022 if (mu == mu_)
return;
1031 constexpr int rowBlockSize = BCRSMatrix::block_type::rows;
1032 constexpr int colBlockSize = BCRSMatrix::block_type::cols;
1034 i < itMatrix_->M()*rowBlockSize; ++i)
1037 const Real& m_entry = m_
1038 [i/rowBlockSize][i/colBlockSize][i%rowBlockSize][i%colBlockSize];
1040 Real& entry = (*itMatrix_)
1041 [i/rowBlockSize][i/colBlockSize][i%rowBlockSize][i%colBlockSize];
1043 entry = m_entry - mu_;
1047 (solver,*itMatrix_);
1054 const unsigned int nIterationsMax_;
1057 const unsigned int verbosity_level_;
1063 const MatrixOperator matrixOperator_;
1064 const ScalingOperator scalingOperator_;
1065 OperatorSum itOperator_;
1069 mutable std::unique_ptr<BCRSMatrix> itMatrix_;
1074 mutable unsigned int nIterations_;
1077 const std::string title_;
1078 const std::string blank_;
A sparse block matrix with compressed row storage.
Definition: bcrsmatrix.hh:423
@ blocklevel
The number of blocklevels the matrix contains.
Definition: bcrsmatrix.hh:465
A::size_type size_type
The type for the index access and the size.
Definition: bcrsmatrix.hh:457
size_type M() const
number of columns (counted in blocks)
Definition: bcrsmatrix.hh:1909
void mv(const X &x, Y &y) const
y = A x
Definition: bcrsmatrix.hh:1580
size_type N() const
number of rows (counted in blocks)
Definition: bcrsmatrix.hh:1903
A vector of blocks with memory management.
Definition: bvector.hh:313
B::field_type field_type
export the type representing the field
Definition: bvector.hh:319
derive error class from the base class in common
Definition: istlexception.hh:16
A linear operator representing the sum of two linear operators.
Definition: poweriteration.hh:85
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
A linear operator.
Definition: operators.hh:62
Adapter to turn a matrix into a linear operator.
Definition: operators.hh:122
A class template for performing some iterative eigenvalue algorithms based on power iteration.
Definition: poweriteration.hh:211
PowerIteration_Algorithms(const PowerIteration_Algorithms &)=delete
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
PowerIteration_Algorithms(const BCRSMatrix &m, const unsigned int nIterationsMax=1000, const unsigned int verbosity_level=0)
Construct from required parameters.
Definition: poweriteration.hh:241
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
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 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
BlockVector::field_type Real
Type of underlying field.
Definition: poweriteration.hh:221
A linear operator scaling vectors by a scalar value. The scalar value can be changed as it is given i...
Definition: poweriteration.hh:40
Helper class for notifying a DUNE-ISTL linear solver about a change of the iteration matrix object in...
Definition: poweriteration.hh:132
FieldTraits< field_type >::real_type two_norm() const
two norm sqrt(sum over squared values of entries)
Definition: bvector.hh:193
block_vector_unmanaged & axpy(const field_type &a, const block_vector_unmanaged &y)
vector space axpy operation
Definition: bvector.hh:124
A few common exception classes.
#define DUNE_THROW(E, m)
Definition: exceptions.hh:216
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
Some generic functions for pretty printing vectors and matrices.
Dune namespace.
Definition: alignment.hh:11
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.
Statistics about the application of an inverse operator.
Definition: solver.hh:32
@ sequential
Category for sequential solvers.
Definition: solvercategory.hh:21
Implementation that works together with iterative ISTL solvers, e.g. Dune::CGSolver or Dune::BiCGSTAB...
Definition: poweriteration.hh:150