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>
49 const field_type& mutable_scaling)
54 virtual void apply (
const X& x, Y& y)
const 81 template <
class OP1,
class OP2>
84 typename OP1::range_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);
109 const domain_type& x, range_type& y)
const 113 op2_.applyscaleadd(1.0,x,temp);
130 template <
typename ISTLLinearSolver,
typename BCRSMatrix>
137 static const bool is_direct_solver
148 template <
bool is_direct_solver,
typename Dummy =
void>
160 template <
typename Dummy>
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;
395 applyInverseIteration(gamma,epsilon,solver,x,lambda);
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;
449 updateShiftMu(gamma,solver);
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 = " 602 updateShiftMu(lambda,solver);
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 <<
").");
782 updateShiftMu(mu,solver);
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_);
1084 #endif // DUNE_ISTL_EIGENVALUE_POWERITERATION_HH BlockVector::field_type Real
Type of underlying field.
Definition: poweriteration.hh:221
const OP2 & op2_
Definition: poweriteration.hh:119
LinearOperatorSum(const OP1 &op1, const OP2 &op2)
Definition: poweriteration.hh:93
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
A linear operator.
Definition: operators.hh:62
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
ScalingLinearOperator< BlockVector > ScalingOperator
Definition: poweriteration.hh:216
A::size_type size_type
The type for the index access and the size.
Definition: bcrsmatrix.hh:457
A vector of blocks with memory management.
Definition: bvector.hh:312
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
OP1::range_type range_type
Definition: poweriteration.hh:88
virtual void applyscaleadd(field_type alpha, const X &x, Y &y) const
apply operator to x, scale and add:
Definition: poweriteration.hh:60
Implementations of the inverse operator interface.
const field_type immutable_scaling_
Definition: poweriteration.hh:68
std::unique_ptr< BCRSMatrix > itMatrix_
Definition: poweriteration.hh:1069
block_vector_unmanaged & axpy(const field_type &a, const block_vector_unmanaged &y)
vector space axpy operation
Definition: bvector.hh:124
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 itOperator_
Definition: poweriteration.hh:1065
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
Templates characterizing the type of a solver.
A linear operator scaling vectors by a scalar value. The scalar value can be changed as it is given i...
Definition: poweriteration.hh:39
Y range_type
Definition: poweriteration.hh:43
X::field_type field_type
Definition: poweriteration.hh:44
Category for sequential solvers.
Definition: solvercategory.hh:21
Some generic functions for pretty printing vectors and matrices.
static void setMatrix(ISTLLinearSolver &, const BCRSMatrix &)
Definition: poweriteration.hh:151
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
OperatorSum IterationOperator
Type of iteration operator (m_ - mu_*I)
Definition: poweriteration.hh:224
const BCRSMatrix & m_
Definition: poweriteration.hh:1053
ScalingLinearOperator(field_type immutable_scaling, const field_type &mutable_scaling)
Definition: poweriteration.hh:48
FieldTraits< field_type >::real_type two_norm() const
two norm sqrt(sum over squared values of entries)
Definition: bvector.hh:193
const std::string blank_
Definition: poweriteration.hh:1078
A sparse block matrix with compressed row storage.
Definition: bcrsmatrix.hh:422
const MatrixOperator matrixOperator_
Definition: poweriteration.hh:1063
A linear operator representing the sum of two linear operators.
Definition: poweriteration.hh:82
domain_type::field_type field_type
Definition: poweriteration.hh:89
Definition: basearray.hh:19
const field_type & mutable_scaling_
Definition: poweriteration.hh:69
void updateShiftMu(const Real &mu, ISTLLinearSolver &solver) const
Update shift mu_, i.e. update iteration operator/matrix (m_ - mu_*I).
Definition: poweriteration.hh:1018
const unsigned int nIterationsMax_
Definition: poweriteration.hh:1054
Adapter to turn a matrix into a linear operator.
Definition: operators.hh:121
Definition: poweriteration.hh:46
OP1::domain_type domain_type
Definition: poweriteration.hh:87
Definition: solvertype.hh:13
The number of blocklevels the matrix contains.
Definition: bcrsmatrix.hh:465
virtual void applyscaleadd(field_type alpha, const domain_type &x, range_type &y) const
Definition: poweriteration.hh:108
const unsigned int verbosity_level_
Definition: poweriteration.hh:1057
Real mu_
Definition: poweriteration.hh:1060
unsigned int nIterations_
Definition: poweriteration.hh:1074
const OP1 & op1_
Definition: poweriteration.hh:118
B::field_type field_type
export the type representing the field
Definition: bvector.hh:319
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
Implementation that works together with iterative ISTL solvers, e.g. Dune::CGSolver or Dune::BiCGSTAB...
Definition: poweriteration.hh:149
X domain_type
Definition: poweriteration.hh:42
Define general, extensible interface for operators. The available implementation wraps a matrix...
LinearOperatorSum< MatrixOperator, ScalingOperator > OperatorSum
Definition: poweriteration.hh:217
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
derive error class from the base class in common
Definition: istlexception.hh:16
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
static void setMatrix(ISTLLinearSolver &solver, const BCRSMatrix &matrix)
Definition: poweriteration.hh:134
const ScalingOperator scalingOperator_
Definition: poweriteration.hh:1064
A class template for performing some iterative eigenvalue algorithms based on power iteration...
Definition: poweriteration.hh:210
static void setMatrix(ISTLLinearSolver &solver, const BCRSMatrix &matrix)
Definition: poweriteration.hh:163
Helper class for notifying a DUNE-ISTL linear solver about a change of the iteration matrix object in...
Definition: poweriteration.hh:131
Statistics about the application of an inverse operator.
Definition: solver.hh:31
const std::string title_
Definition: poweriteration.hh:1077
unsigned int getIterationCount() const
Return the number of iterations in last application of an algorithm.
Definition: poweriteration.hh:994
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