3 #ifndef DUNE_ISTL_EIGENVALUE_POWERITERATION_HH
4 #define DUNE_ISTL_EIGENVALUE_POWERITERATION_HH
17 #include <dune/common/exceptions.hh>
43 template <
class X,
class Y = X>
51 ScalingLinearOperator (field_type immutable_scaling,
52 const field_type& mutable_scaling)
53 : immutable_scaling_(immutable_scaling),
54 mutable_scaling_(mutable_scaling)
57 virtual void apply (
const X& x, Y& y)
const
60 y *= immutable_scaling_*mutable_scaling_;
63 virtual void applyscaleadd (field_type alpha,
const X& x, Y& y)
const
66 temp *= immutable_scaling_*mutable_scaling_;
77 const field_type immutable_scaling_;
78 const field_type& mutable_scaling_;
90 template <
class OP1,
class OP2>
91 class LinearOperatorSum
93 typename OP1::range_type>
96 typedef typename OP1::domain_type domain_type;
97 typedef typename OP1::range_type range_type;
98 typedef typename domain_type::field_type field_type;
100 LinearOperatorSum (
const OP1& op1,
const OP2& op2)
101 : op1_(op1), op2_(op2)
103 static_assert(std::is_same<typename OP2::domain_type,domain_type>::value,
104 "Domain type of both operators doesn't match!");
105 static_assert(std::is_same<typename OP2::range_type,range_type>::value,
106 "Range type of both operators doesn't match!");
109 virtual void apply (
const domain_type& x, range_type& y)
const
112 op2_.applyscaleadd(1.0,x,y);
116 const domain_type& x, range_type& y)
const
120 op2_.applyscaleadd(1.0,x,temp);
172 template <
typename BCRSMatrix,
typename BlockVector>
180 typedef Impl::LinearOperatorSum<MatrixOperator,ScalingOperator>
OperatorSum;
205 const unsigned int nIterationsMax = 1000,
206 const unsigned int verbosity_level = 0)
214 title_(
" PowerIteration_Algorithms: "),
219 (blockLevel<BCRSMatrix>() == 2,
220 "Only BCRSMatrices with blocklevel 2 are supported.");
224 (BCRSMatrix::block_type::rows == BCRSMatrix::block_type::cols,
225 "Only BCRSMatrices with square blocks are supported.");
228 const int nrows =
m_.
M() * BCRSMatrix::block_type::rows;
229 const int ncols =
m_.
N() * BCRSMatrix::block_type::cols;
232 << nrows <<
"x" << ncols <<
").");
264 <<
"Performing power iteration approximating "
265 <<
"the dominant eigenvalue." << std::endl;
272 x *= (1.0 / x.two_norm());
274 Real r_norm = std::numeric_limits<Real>::max();
276 while (r_norm > epsilon)
282 <<
"(║residual║_2 = " << r_norm <<
", epsilon = "
288 x *= (1.0 / y.two_norm());
296 temp.axpy(-lambda,x);
297 r_norm = temp.two_norm();
301 std::cout <<
blank_ << std::left
303 <<
" (║residual║_2 = " << std::setw(11) << r_norm
304 <<
"): λ = " << lambda << std::endl
305 << std::resetiosflags(std::ios::left);
311 std::cout <<
blank_ <<
"Result ("
313 <<
"║residual║_2 = " << r_norm <<
"): "
314 <<
"λ = " << lambda << std::endl;
351 template <
typename ISTLLinearSolver,
352 bool avoidLinSolverCrime =
false>
354 ISTLLinearSolver& solver,
357 constexpr
Real gamma = 0.0;
390 template <
typename ISTLLinearSolver,
391 bool avoidLinSolverCrime =
false>
394 ISTLLinearSolver& solver,
402 std::cout <<
"Performing inverse iteration approximating "
403 <<
"the least dominant eigenvalue." << std::endl;
405 std::cout <<
"Performing inverse iteration with shift "
406 <<
"gamma = " << gamma <<
" approximating the "
407 <<
"eigenvalue closest to gamma." << std::endl;
423 x *= (1.0 / x.two_norm());
424 Real r_norm = std::numeric_limits<Real>::max();
426 while (r_norm > epsilon)
431 << (gamma != 0.0 ?
"with shift " :
"") <<
"did not "
433 <<
"(║residual║_2 = " << r_norm <<
", epsilon = "
440 solver.apply(y,temp,solver_statistics);
443 y_norm = y.two_norm();
446 if (avoidLinSolverCrime)
451 lambda = (y * temp) / (y_norm * y_norm);
455 temp.axpy(-lambda,y);
456 r_norm = temp.two_norm() / y_norm;
462 lambda = gamma + (y * x) / (y_norm * y_norm);
466 temp = x; temp.axpy(gamma-lambda,y);
467 r_norm = temp.two_norm() / y_norm;
477 std::cout <<
blank_ << std::left
479 <<
" (║residual║_2 = " << std::setw(11) << r_norm
480 <<
"): λ = " << lambda << std::endl
481 << std::resetiosflags(std::ios::left);
487 std::cout <<
blank_ <<
"Result ("
489 <<
"║residual║_2 = " << r_norm <<
"): "
490 <<
"λ = " << lambda << std::endl;
529 template <
typename ISTLLinearSolver,
530 bool avoidLinSolverCrime =
false>
532 ISTLLinearSolver& solver,
538 <<
"Performing Rayleigh quotient iteration for "
539 <<
"estimated eigenvalue " << lambda <<
"." << std::endl;
551 x *= (1.0 / x.two_norm());
552 Real r_norm = std::numeric_limits<Real>::max();
554 while (r_norm > epsilon)
560 <<
"(║residual║_2 = " << r_norm <<
", epsilon = "
571 solver.apply(y,temp,solver_statistics);
574 y_norm = y.two_norm();
577 if (avoidLinSolverCrime)
582 lambda = (y * temp) / (y_norm * y_norm);
586 temp.axpy(-lambda,y);
587 r_norm = temp.two_norm() / y_norm;
593 lambda_update = (y * x) / (y_norm * y_norm);
594 lambda += lambda_update;
598 temp = x; temp.axpy(-lambda_update,y);
599 r_norm = temp.two_norm() / y_norm;
609 std::cout <<
blank_ << std::left
611 <<
" (║residual║_2 = " << std::setw(11) << r_norm
612 <<
"): λ = " << lambda << std::endl
613 << std::resetiosflags(std::ios::left);
619 std::cout <<
blank_ <<
"Result ("
621 <<
"║residual║_2 = " << r_norm <<
"): "
622 <<
"λ = " << lambda << std::endl;
687 template <
typename ISTLLinearSolver,
688 bool avoidLinSolverCrime =
false>
691 ISTLLinearSolver& solver,
692 const Real& delta,
const std::size_t& m,
703 <<
"Performing TLIME iteration for "
704 <<
"estimated eigenvalue in the "
705 <<
"interval (" << gamma - eta <<
","
706 << gamma + eta <<
")." << std::endl;
722 x_s *= (1.0 / x_s.two_norm());
725 r_norm = std::numeric_limits<Real>::max();
727 while (r_norm > epsilon)
733 <<
" iterations (║residual║_2 = " << r_norm
734 <<
", epsilon = " << epsilon <<
").");
750 solver.apply(y,temp,solver_statistics);
754 omega = (1.0 / y.two_norm());
760 if (avoidLinSolverCrime)
765 mu_s = (y * temp) * (omega * omega);
771 q_norm = temp.two_norm() * omega;
774 r_norm = q_norm*q_norm - (gamma-mu_s)*(gamma-mu_s);
780 r_norm = std::sqrt(r_norm);
786 temp.axpy(gamma-mu_s,y);
787 r_norm = temp.two_norm() * omega;
796 mu_s = gamma + (y * x_s) * (omega * omega);
801 mu_s_update = (y * x_s) * (omega * omega);
816 temp = x_s; temp.axpy(mu_s-gamma,y);
817 q_norm = temp.two_norm() * omega;
827 temp = x_s; temp.axpy(gamma-lambda,y);
828 r_norm = temp.two_norm() * omega;
833 temp = x_s; temp.axpy(-mu_s_update,y);
834 r_norm = temp.two_norm() * omega;
840 x_s = y; x_s *= omega;
847 std::cout <<
blank_ <<
"iteration "
849 <<
" (" << (doRQI ?
"RQI," :
"II, ")
850 <<
" " << (doRQI ?
"—>" :
" ") <<
" "
851 <<
"║r║_2 = " << std::setw(11) << r_norm
852 <<
", " << (doRQI ?
" " :
"—>") <<
" "
853 <<
"║q║_2 = " << std::setw(11) << q_norm
854 <<
"): λ = " << lambda << std::endl
855 << std::resetiosflags(std::ios::left);
858 if (!doRQI && q_norm < eta)
864 assert(std::abs(mu_s-gamma) < eta);
872 if (!extrnl && doRQI && std::abs(mu_s-gamma) >= eta)
878 if (extrnl && !doRQI)
883 std::abs(mu_s - mu_s_old) / std::abs(mu_s) < delta)
899 std::cout <<
blank_ <<
"Interval "
900 <<
"(" << gamma - eta <<
"," << gamma + eta
901 <<
") is free of eigenvalues, approximating "
902 <<
"the closest eigenvalue." << std::endl;
903 std::cout <<
blank_ <<
"Result ("
905 <<
"║residual║_2 = " << r_norm <<
"): "
906 <<
"λ = " << lambda << std::endl;
980 template <
typename ISTLLinearSolver>
982 ISTLLinearSolver& solver)
const
985 if (mu ==
mu_)
return;
994 constexpr
int rowBlockSize = BCRSMatrix::block_type::rows;
995 constexpr
int colBlockSize = BCRSMatrix::block_type::cols;
1001 [i/rowBlockSize][i/colBlockSize][i%rowBlockSize][i%colBlockSize];
1003 Real& entry = (*itMatrix_)
1004 [i/rowBlockSize][i/colBlockSize][i%rowBlockSize][i%colBlockSize];
1006 entry = m_entry -
mu_;
Helper functions for determining the vector/matrix block level.
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:86
Definition: allocator.hh:9
A sparse block matrix with compressed row storage.
Definition: bcrsmatrix.hh:464
A::size_type size_type
The type for the index access and the size.
Definition: bcrsmatrix.hh:498
size_type M() const
number of columns (counted in blocks)
Definition: bcrsmatrix.hh:1976
void mv(const X &x, Y &y) const
y = A x
Definition: bcrsmatrix.hh:1610
size_type N() const
number of rows (counted in blocks)
Definition: bcrsmatrix.hh:1970
A vector of blocks with memory management.
Definition: bvector.hh:393
typename Imp::BlockTraits< B >::field_type field_type
export the type representing the field
Definition: bvector.hh:399
Iterative eigenvalue algorithms based on power iteration.
Definition: poweriteration.hh:174
std::unique_ptr< BCRSMatrix > itMatrix_
Definition: poweriteration.hh:1032
PowerIteration_Algorithms(const PowerIteration_Algorithms &)=delete
Impl::ScalingLinearOperator< BlockVector > ScalingOperator
Definition: poweriteration.hh:179
const std::string blank_
Definition: poweriteration.hh:1041
Impl::LinearOperatorSum< MatrixOperator, ScalingOperator > OperatorSum
Definition: poweriteration.hh:180
Dune::MatrixAdapter< BCRSMatrix, BlockVector, BlockVector > MatrixOperator
Definition: poweriteration.hh:178
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:353
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:689
OperatorSum itOperator_
Definition: poweriteration.hh:1028
const BCRSMatrix & m_
Definition: poweriteration.hh:1016
PowerIteration_Algorithms(const BCRSMatrix &m, const unsigned int nIterationsMax=1000, const unsigned int verbosity_level=0)
Construct from required parameters.
Definition: poweriteration.hh:204
const unsigned int nIterationsMax_
Definition: poweriteration.hh:1017
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:258
OperatorSum IterationOperator
Type of iteration operator (m_ - mu_*I)
Definition: poweriteration.hh:187
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:531
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:392
const ScalingOperator scalingOperator_
Definition: poweriteration.hh:1027
void updateShiftMu(const Real &mu, ISTLLinearSolver &solver) const
Update shift mu_, i.e. update iteration operator/matrix (m_ - mu_*I).
Definition: poweriteration.hh:981
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:943
unsigned int getIterationCount() const
Return the number of iterations in last application of an algorithm.
Definition: poweriteration.hh:957
const MatrixOperator matrixOperator_
Definition: poweriteration.hh:1026
unsigned int nIterations_
Definition: poweriteration.hh:1037
const unsigned int verbosity_level_
Definition: poweriteration.hh:1020
IterationOperator & getIterationOperator()
Return the iteration operator (m_ - mu_*I).
Definition: poweriteration.hh:923
const std::string title_
Definition: poweriteration.hh:1040
PowerIteration_Algorithms & operator=(const PowerIteration_Algorithms &)=delete
BlockVector::field_type Real
Type of underlying field.
Definition: poweriteration.hh:184
Real mu_
Definition: poweriteration.hh:1023
derive error class from the base class in common
Definition: istlexception.hh:17
A linear operator.
Definition: operators.hh:65
X::field_type field_type
The field type of the operator.
Definition: operators.hh:72
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:70
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:68
Adapter to turn a matrix into a linear operator.
Definition: operators.hh:135
Statistics about the application of an inverse operator.
Definition: solver.hh:46
static void setMatrix(ISTLLinearSolver &solver, const BCRSMatrix &matrix)
Definition: solver.hh:522
Category
Definition: solvercategory.hh:21
@ sequential
Category for sequential solvers.
Definition: solvercategory.hh:23