5#ifndef DUNE_ISTL_EIGENVALUE_POWERITERATION_HH
6#define DUNE_ISTL_EIGENVALUE_POWERITERATION_HH
23#include <dune/istl/solvercategory.hh>
26#include <dune/istl/istlexception.hh>
45 template <
class X,
class Y = X>
53 ScalingLinearOperator (field_type immutable_scaling,
54 const field_type& mutable_scaling)
55 : immutable_scaling_(immutable_scaling),
56 mutable_scaling_(mutable_scaling)
59 void apply (
const X& x, Y& y)
const override
62 y *= immutable_scaling_*mutable_scaling_;
65 void applyscaleadd (field_type alpha,
const X& x, Y& y)
const override
68 temp *= immutable_scaling_*mutable_scaling_;
79 const field_type immutable_scaling_;
80 const field_type& mutable_scaling_;
92 template <
class OP1,
class OP2>
93 class LinearOperatorSum
95 typename OP1::range_type>
98 typedef typename OP1::domain_type domain_type;
99 typedef typename OP1::range_type range_type;
100 typedef typename domain_type::field_type field_type;
102 LinearOperatorSum (
const OP1& op1,
const OP2& op2)
103 : op1_(op1), op2_(op2)
105 static_assert(std::is_same<typename OP2::domain_type,domain_type>::value,
106 "Domain type of both operators doesn't match!");
107 static_assert(std::is_same<typename OP2::range_type,range_type>::value,
108 "Range type of both operators doesn't match!");
111 void apply (
const domain_type& x, range_type& y)
const override
114 op2_.applyscaleadd(1.0,x,y);
118 const domain_type& x, range_type& y)
const override
122 op2_.applyscaleadd(1.0,x,temp);
174 template <
typename BCRSMatrix,
typename BlockVector>
181 typedef Impl::ScalingLinearOperator<BlockVector> ScalingOperator;
182 typedef Impl::LinearOperatorSum<MatrixOperator,ScalingOperator> OperatorSum;
207 const unsigned int nIterationsMax = 1000,
208 const unsigned int verbosity_level = 0)
209 : m_(m), nIterationsMax_(nIterationsMax),
210 verbosity_level_(verbosity_level),
213 scalingOperator_(-1.0,mu_),
214 itOperator_(matrixOperator_,scalingOperator_),
216 title_(
" PowerIteration_Algorithms: "),
217 blank_(title_.length(),
' ')
221 (blockLevel<BCRSMatrix>() == 2,
222 "Only BCRSMatrices with blocklevel 2 are supported.");
226 (BCRSMatrix::block_type::rows == BCRSMatrix::block_type::cols,
227 "Only BCRSMatrices with square blocks are supported.");
230 const int nrows = m_.
M() * BCRSMatrix::block_type::rows;
231 const int ncols = m_.
N() * BCRSMatrix::block_type::cols;
234 << nrows <<
"x" << ncols <<
").");
264 if (verbosity_level_ > 0)
266 <<
"Performing power iteration approximating "
267 <<
"the dominant eigenvalue." << std::endl;
274 x *= (1.0 / x.two_norm());
278 while (r_norm > epsilon)
281 if (++nIterations_ > nIterationsMax_)
283 <<
"in " << nIterationsMax_ <<
" iterations "
284 <<
"(║residual║_2 = " << r_norm <<
", epsilon = "
290 x *= (1.0 / y.two_norm());
298 temp.axpy(-lambda,x);
299 r_norm = temp.two_norm();
302 if (verbosity_level_ > 1)
303 std::cout << blank_ << std::left
304 <<
"iteration " << std::setw(3) << nIterations_
305 <<
" (║residual║_2 = " << std::setw(11) << r_norm
306 <<
"): λ = " << lambda << std::endl
307 << std::resetiosflags(std::ios::left);
311 if (verbosity_level_ > 0)
313 std::cout << blank_ <<
"Result ("
314 <<
"#iterations = " << nIterations_ <<
", "
315 <<
"║residual║_2 = " << r_norm <<
"): "
316 <<
"λ = " << lambda << std::endl;
317 if (verbosity_level_ > 2)
353 template <
typename ISTLLinearSolver,
354 bool avoidLinSolverCrime =
false>
356 ISTLLinearSolver& solver,
359 constexpr Real gamma = 0.0;
392 template <
typename ISTLLinearSolver,
393 bool avoidLinSolverCrime =
false>
396 ISTLLinearSolver& solver,
400 if (verbosity_level_ > 0)
404 std::cout <<
"Performing inverse iteration approximating "
405 <<
"the least dominant eigenvalue." << std::endl;
407 std::cout <<
"Performing inverse iteration with shift "
408 <<
"gamma = " << gamma <<
" approximating the "
409 <<
"eigenvalue closest to gamma." << std::endl;
425 x *= (1.0 / x.two_norm());
428 while (r_norm > epsilon)
431 if (++nIterations_ > nIterationsMax_)
433 << (gamma != 0.0 ?
"with shift " :
"") <<
"did not "
434 <<
"converge in " << nIterationsMax_ <<
" iterations "
435 <<
"(║residual║_2 = " << r_norm <<
", epsilon = "
442 solver.apply(y,temp,solver_statistics);
445 y_norm = y.two_norm();
448 if (avoidLinSolverCrime)
453 lambda = (y * temp) / (y_norm * y_norm);
457 temp.axpy(-lambda,y);
458 r_norm = temp.two_norm() / y_norm;
464 lambda = gamma + (y * x) / (y_norm * y_norm);
468 temp = x; temp.axpy(gamma-lambda,y);
469 r_norm = temp.two_norm() / y_norm;
478 if (verbosity_level_ > 1)
479 std::cout << blank_ << std::left
480 <<
"iteration " << std::setw(3) << nIterations_
481 <<
" (║residual║_2 = " << std::setw(11) << r_norm
482 <<
"): λ = " << lambda << std::endl
483 << std::resetiosflags(std::ios::left);
487 if (verbosity_level_ > 0)
489 std::cout << blank_ <<
"Result ("
490 <<
"#iterations = " << nIterations_ <<
", "
491 <<
"║residual║_2 = " << r_norm <<
"): "
492 <<
"λ = " << lambda << std::endl;
493 if (verbosity_level_ > 2)
531 template <
typename ISTLLinearSolver,
532 bool avoidLinSolverCrime =
false>
534 ISTLLinearSolver& solver,
538 if (verbosity_level_ > 0)
540 <<
"Performing Rayleigh quotient iteration for "
541 <<
"estimated eigenvalue " << lambda <<
"." << std::endl;
553 x *= (1.0 / x.two_norm());
556 while (r_norm > epsilon)
559 if (++nIterations_ > nIterationsMax_)
561 <<
"converge in " << nIterationsMax_ <<
" iterations "
562 <<
"(║residual║_2 = " << r_norm <<
", epsilon = "
573 solver.apply(y,temp,solver_statistics);
576 y_norm = y.two_norm();
579 if (avoidLinSolverCrime)
584 lambda = (y * temp) / (y_norm * y_norm);
588 temp.axpy(-lambda,y);
589 r_norm = temp.two_norm() / y_norm;
595 lambda_update = (y * x) / (y_norm * y_norm);
596 lambda += lambda_update;
600 temp = x; temp.axpy(-lambda_update,y);
601 r_norm = temp.two_norm() / y_norm;
610 if (verbosity_level_ > 1)
611 std::cout << blank_ << std::left
612 <<
"iteration " << std::setw(3) << nIterations_
613 <<
" (║residual║_2 = " << std::setw(11) << r_norm
614 <<
"): λ = " << lambda << std::endl
615 << std::resetiosflags(std::ios::left);
619 if (verbosity_level_ > 0)
621 std::cout << blank_ <<
"Result ("
622 <<
"#iterations = " << nIterations_ <<
", "
623 <<
"║residual║_2 = " << r_norm <<
"): "
624 <<
"λ = " << lambda << std::endl;
625 if (verbosity_level_ > 2)
689 template <
typename ISTLLinearSolver,
690 bool avoidLinSolverCrime =
false>
693 ISTLLinearSolver& solver,
694 const Real& delta,
const std::size_t& m,
703 if (verbosity_level_ > 0)
705 <<
"Performing TLIME iteration for "
706 <<
"estimated eigenvalue in the "
707 <<
"interval (" << gamma - eta <<
","
708 << gamma + eta <<
")." << std::endl;
724 x_s *= (1.0 / x_s.two_norm());
729 while (r_norm > epsilon)
732 if (++nIterations_ > nIterationsMax_)
734 <<
"converge in " << nIterationsMax_
735 <<
" iterations (║residual║_2 = " << r_norm
736 <<
", epsilon = " << epsilon <<
").");
752 solver.apply(y,temp,solver_statistics);
756 omega = (1.0 / y.two_norm());
762 if (avoidLinSolverCrime)
767 mu_s = (y * temp) * (omega * omega);
773 q_norm = temp.two_norm() * omega;
776 r_norm = q_norm*q_norm - (gamma-mu_s)*(gamma-mu_s);
782 r_norm = std::sqrt(r_norm);
788 temp.axpy(gamma-mu_s,y);
789 r_norm = temp.two_norm() * omega;
798 mu_s = gamma + (y * x_s) * (omega * omega);
803 mu_s_update = (y * x_s) * (omega * omega);
818 temp = x_s; temp.axpy(mu_s-gamma,y);
819 q_norm = temp.two_norm() * omega;
829 temp = x_s; temp.axpy(gamma-lambda,y);
830 r_norm = temp.two_norm() * omega;
835 temp = x_s; temp.axpy(-mu_s_update,y);
836 r_norm = temp.two_norm() * omega;
842 x_s = y; x_s *= omega;
848 if (verbosity_level_ > 1)
849 std::cout << blank_ <<
"iteration "
850 << std::left << std::setw(3) << nIterations_
851 <<
" (" << (doRQI ?
"RQI," :
"II, ")
852 <<
" " << (doRQI ?
"—>" :
" ") <<
" "
853 <<
"║r║_2 = " << std::setw(11) << r_norm
854 <<
", " << (doRQI ?
" " :
"—>") <<
" "
855 <<
"║q║_2 = " << std::setw(11) << q_norm
856 <<
"): λ = " << lambda << std::endl
857 << std::resetiosflags(std::ios::left);
860 if (!doRQI && q_norm < eta)
866 assert(std::abs(mu_s-gamma) < eta);
874 if (!extrnl && doRQI && std::abs(mu_s-gamma) >= eta)
880 if (extrnl && !doRQI)
884 if (nIterations_ >= m &&
885 std::abs(mu_s - mu_s_old) / std::abs(mu_s) < delta)
898 if (verbosity_level_ > 0)
901 std::cout << blank_ <<
"Interval "
902 <<
"(" << gamma - eta <<
"," << gamma + eta
903 <<
") is free of eigenvalues, approximating "
904 <<
"the closest eigenvalue." << std::endl;
905 std::cout << blank_ <<
"Result ("
906 <<
"#iterations = " << nIterations_ <<
", "
907 <<
"║residual║_2 = " << r_norm <<
"): "
908 <<
"λ = " << lambda << std::endl;
909 if (verbosity_level_ > 2)
949 itMatrix_ = std::make_unique<BCRSMatrix>(m_);
961 if (nIterations_ == 0)
982 template <
typename ISTLLinearSolver>
984 ISTLLinearSolver& solver)
const
987 if (mu == mu_)
return;
996 constexpr int rowBlockSize = BCRSMatrix::block_type::rows;
997 constexpr int colBlockSize = BCRSMatrix::block_type::cols;
999 i < itMatrix_->M()*rowBlockSize; ++i)
1002 const Real& m_entry = m_
1003 [i/rowBlockSize][i/colBlockSize][i%rowBlockSize][i%colBlockSize];
1005 Real& entry = (*itMatrix_)
1006 [i/rowBlockSize][i/colBlockSize][i%rowBlockSize][i%colBlockSize];
1008 entry = m_entry - mu_;
1012 (solver,*itMatrix_);
1019 const unsigned int nIterationsMax_;
1022 const unsigned int verbosity_level_;
1028 const MatrixOperator matrixOperator_;
1029 const ScalingOperator scalingOperator_;
1030 OperatorSum itOperator_;
1034 mutable std::unique_ptr<BCRSMatrix> itMatrix_;
1039 mutable unsigned int nIterations_;
1042 const std::string title_;
1043 const std::string blank_;
Helper functions for determining the vector/matrix block level.
A sparse block matrix with compressed row storage.
Definition: bcrsmatrix.hh:466
A::size_type size_type
The type for the index access and the size.
Definition: bcrsmatrix.hh:497
size_type M() const
number of columns (counted in blocks)
Definition: bcrsmatrix.hh:2007
void mv(const X &x, Y &y) const
y = A x
Definition: bcrsmatrix.hh:1641
size_type N() const
number of rows (counted in blocks)
Definition: bcrsmatrix.hh:2001
A vector of blocks with memory management.
Definition: bvector.hh:392
typename Imp::BlockTraits< B >::field_type field_type
export the type representing the field
Definition: bvector.hh:398
derive error class from the base class in common
Definition: istlexception.hh:19
A linear operator.
Definition: operators.hh:68
X::field_type field_type
The field type of the operator.
Definition: operators.hh:75
virtual void applyscaleadd(field_type alpha, const X &x, X &y) const=0
apply operator to x, scale and add:
virtual SolverCategory::Category category() const=0
Category of the linear operator (see SolverCategory::Category)
X range_type
The type of the range of the operator.
Definition: operators.hh:73
virtual void apply(const X &x, X &y) const=0
apply operator to x: The input vector is consistent and the output must also be consistent on the in...
X domain_type
The type of the domain of the operator.
Definition: operators.hh:71
Adapter to turn a matrix into a linear operator.
Definition: operators.hh:135
Iterative eigenvalue algorithms based on power iteration.
Definition: poweriteration.hh:176
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:355
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:691
IterationOperator & getIterationOperator()
Return the iteration operator (m_ - mu_*I).
Definition: poweriteration.hh:925
PowerIteration_Algorithms(const BCRSMatrix &m, const unsigned int nIterationsMax=1000, const unsigned int verbosity_level=0)
Construct from required parameters.
Definition: poweriteration.hh:206
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
OperatorSum IterationOperator
Type of iteration operator (m_ - mu_*I)
Definition: poweriteration.hh:189
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
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
void updateShiftMu(const Real &mu, ISTLLinearSolver &solver) const
Update shift mu_, i.e. update iteration operator/matrix (m_ - mu_*I).
Definition: poweriteration.hh:983
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:959
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
BlockVector::field_type Real
Type of underlying field.
Definition: poweriteration.hh:186
Helper class for notifying a DUNE-ISTL linear solver about a change of the iteration matrix object in...
Definition: solver.hh:524
A few common exception classes.
Some generic functions for pretty printing vectors and matrices.
#define DUNE_THROW(E, m)
Definition: exceptions.hh:218
constexpr auto max
Function object that returns the greater of the given values.
Definition: hybridutilities.hh:484
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 namespace.
Definition: alignedallocator.hh:13
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:50
Category
Definition: solvercategory.hh:23
@ sequential
Category for sequential solvers.
Definition: solvercategory.hh:25