3 #ifndef DUNE_ISTL_EIGENVALUE_POWERITERATION_HH 4 #define DUNE_ISTL_EIGENVALUE_POWERITERATION_HH 17 #include <dune/common/exceptions.hh> 42 template <
class X,
class Y = X>
50 ScalingLinearOperator (field_type immutable_scaling,
51 const field_type& mutable_scaling)
52 : immutable_scaling_(immutable_scaling),
53 mutable_scaling_(mutable_scaling)
56 virtual void apply (
const X& x, Y& y)
const 59 y *= immutable_scaling_*mutable_scaling_;
62 virtual void applyscaleadd (field_type alpha,
const X& x, Y& y)
const 65 temp *= immutable_scaling_*mutable_scaling_;
76 const field_type immutable_scaling_;
77 const field_type& mutable_scaling_;
89 template <
class OP1,
class OP2>
90 class LinearOperatorSum
92 typename OP1::range_type>
95 typedef typename OP1::domain_type domain_type;
96 typedef typename OP1::range_type range_type;
97 typedef typename domain_type::field_type field_type;
99 LinearOperatorSum (
const OP1& op1,
const OP2& op2)
100 : op1_(op1), op2_(op2)
102 static_assert(std::is_same<typename OP2::domain_type,domain_type>::value,
103 "Domain type of both operators doesn't match!");
104 static_assert(std::is_same<typename OP2::range_type,range_type>::value,
105 "Range type of both operators doesn't match!");
108 virtual void apply (
const domain_type& x, range_type& y)
const 111 op2_.applyscaleadd(1.0,x,y);
114 virtual void applyscaleadd (field_type alpha,
115 const domain_type& x, range_type& y)
const 119 op2_.applyscaleadd(1.0,x,temp);
171 template <
typename BCRSMatrix,
typename BlockVector>
179 typedef Impl::LinearOperatorSum<MatrixOperator,ScalingOperator>
OperatorSum;
204 const unsigned int nIterationsMax = 1000,
205 const unsigned int verbosity_level = 0)
206 : m_(m), nIterationsMax_(nIterationsMax),
207 verbosity_level_(verbosity_level),
210 scalingOperator_(-1.0,mu_),
211 itOperator_(matrixOperator_,scalingOperator_),
213 title_(
" PowerIteration_Algorithms: "),
214 blank_(title_.length(),
' ')
219 "Only BCRSMatrices with blocklevel 2 are supported.");
223 (BCRSMatrix::block_type::rows == BCRSMatrix::block_type::cols,
224 "Only BCRSMatrices with square blocks are supported.");
227 const int nrows = m_.M() * BCRSMatrix::block_type::rows;
228 const int ncols = m_.N() * BCRSMatrix::block_type::cols;
231 << nrows <<
"x" << ncols <<
").");
261 if (verbosity_level_ > 0)
263 <<
"Performing power iteration approximating " 264 <<
"the dominant eigenvalue." << std::endl;
271 x *= (1.0 / x.two_norm());
273 Real r_norm = std::numeric_limits<Real>::max();
275 while (r_norm > epsilon)
278 if (++nIterations_ > nIterationsMax_)
280 <<
"in " << nIterationsMax_ <<
" iterations " 281 <<
"(║residual║_2 = " << r_norm <<
", epsilon = " 287 x *= (1.0 / y.two_norm());
295 temp.axpy(-lambda,x);
296 r_norm = temp.two_norm();
299 if (verbosity_level_ > 1)
300 std::cout << blank_ << std::left
301 <<
"iteration " << std::setw(3) << nIterations_
302 <<
" (║residual║_2 = " << std::setw(11) << r_norm
303 <<
"): λ = " << lambda << std::endl
304 << std::resetiosflags(std::ios::left);
308 if (verbosity_level_ > 0)
310 std::cout << blank_ <<
"Result (" 311 <<
"#iterations = " << nIterations_ <<
", " 312 <<
"║residual║_2 = " << r_norm <<
"): " 313 <<
"λ = " << lambda << std::endl;
314 if (verbosity_level_ > 2)
350 template <
typename ISTLLinearSolver,
351 bool avoidLinSolverCrime =
false>
353 ISTLLinearSolver& solver,
356 constexpr Real gamma = 0.0;
357 applyInverseIteration(gamma,epsilon,solver,x,lambda);
389 template <
typename ISTLLinearSolver,
390 bool avoidLinSolverCrime =
false>
393 ISTLLinearSolver& solver,
397 if (verbosity_level_ > 0)
401 std::cout <<
"Performing inverse iteration approximating " 402 <<
"the least dominant eigenvalue." << std::endl;
404 std::cout <<
"Performing inverse iteration with shift " 405 <<
"gamma = " << gamma <<
" approximating the " 406 <<
"eigenvalue closest to gamma." << std::endl;
411 updateShiftMu(gamma,solver);
422 x *= (1.0 / x.two_norm());
423 Real r_norm = std::numeric_limits<Real>::max();
425 while (r_norm > epsilon)
428 if (++nIterations_ > nIterationsMax_)
430 << (gamma != 0.0 ?
"with shift " :
"") <<
"did not " 431 <<
"converge in " << nIterationsMax_ <<
" iterations " 432 <<
"(║residual║_2 = " << r_norm <<
", epsilon = " 439 solver.apply(y,temp,solver_statistics);
442 y_norm = y.two_norm();
445 if (avoidLinSolverCrime)
450 lambda = (y * temp) / (y_norm * y_norm);
454 temp.axpy(-lambda,y);
455 r_norm = temp.two_norm() / y_norm;
461 lambda = gamma + (y * x) / (y_norm * y_norm);
465 temp = x; temp.axpy(gamma-lambda,y);
466 r_norm = temp.two_norm() / y_norm;
475 if (verbosity_level_ > 1)
476 std::cout << blank_ << std::left
477 <<
"iteration " << std::setw(3) << nIterations_
478 <<
" (║residual║_2 = " << std::setw(11) << r_norm
479 <<
"): λ = " << lambda << std::endl
480 << std::resetiosflags(std::ios::left);
484 if (verbosity_level_ > 0)
486 std::cout << blank_ <<
"Result (" 487 <<
"#iterations = " << nIterations_ <<
", " 488 <<
"║residual║_2 = " << r_norm <<
"): " 489 <<
"λ = " << lambda << std::endl;
490 if (verbosity_level_ > 2)
528 template <
typename ISTLLinearSolver,
529 bool avoidLinSolverCrime =
false>
531 ISTLLinearSolver& solver,
535 if (verbosity_level_ > 0)
537 <<
"Performing Rayleigh quotient iteration for " 538 <<
"estimated eigenvalue " << lambda <<
"." << std::endl;
550 x *= (1.0 / x.two_norm());
551 Real r_norm = std::numeric_limits<Real>::max();
553 while (r_norm > epsilon)
556 if (++nIterations_ > nIterationsMax_)
558 <<
"converge in " << nIterationsMax_ <<
" iterations " 559 <<
"(║residual║_2 = " << r_norm <<
", epsilon = " 564 updateShiftMu(lambda,solver);
570 solver.apply(y,temp,solver_statistics);
573 y_norm = y.two_norm();
576 if (avoidLinSolverCrime)
581 lambda = (y * temp) / (y_norm * y_norm);
585 temp.axpy(-lambda,y);
586 r_norm = temp.two_norm() / y_norm;
592 lambda_update = (y * x) / (y_norm * y_norm);
593 lambda += lambda_update;
597 temp = x; temp.axpy(-lambda_update,y);
598 r_norm = temp.two_norm() / y_norm;
607 if (verbosity_level_ > 1)
608 std::cout << blank_ << std::left
609 <<
"iteration " << std::setw(3) << nIterations_
610 <<
" (║residual║_2 = " << std::setw(11) << r_norm
611 <<
"): λ = " << lambda << std::endl
612 << std::resetiosflags(std::ios::left);
616 if (verbosity_level_ > 0)
618 std::cout << blank_ <<
"Result (" 619 <<
"#iterations = " << nIterations_ <<
", " 620 <<
"║residual║_2 = " << r_norm <<
"): " 621 <<
"λ = " << lambda << std::endl;
622 if (verbosity_level_ > 2)
686 template <
typename ISTLLinearSolver,
687 bool avoidLinSolverCrime =
false>
690 ISTLLinearSolver& solver,
691 const Real& delta,
const std::size_t& m,
700 if (verbosity_level_ > 0)
702 <<
"Performing TLIME iteration for " 703 <<
"estimated eigenvalue in the " 704 <<
"interval (" << gamma - eta <<
"," 705 << gamma + eta <<
")." << std::endl;
721 x_s *= (1.0 / x_s.two_norm());
724 r_norm = std::numeric_limits<Real>::max();
726 while (r_norm > epsilon)
729 if (++nIterations_ > nIterationsMax_)
731 <<
"converge in " << nIterationsMax_
732 <<
" iterations (║residual║_2 = " << r_norm
733 <<
", epsilon = " << epsilon <<
").");
744 updateShiftMu(mu,solver);
749 solver.apply(y,temp,solver_statistics);
753 omega = (1.0 / y.two_norm());
759 if (avoidLinSolverCrime)
764 mu_s = (y * temp) * (omega * omega);
770 q_norm = temp.two_norm() * omega;
773 r_norm = q_norm*q_norm - (gamma-mu_s)*(gamma-mu_s);
779 r_norm = std::sqrt(r_norm);
785 temp.axpy(gamma-mu_s,y);
786 r_norm = temp.two_norm() * omega;
795 mu_s = gamma + (y * x_s) * (omega * omega);
800 mu_s_update = (y * x_s) * (omega * omega);
815 temp = x_s; temp.axpy(mu_s-gamma,y);
816 q_norm = temp.two_norm() * omega;
826 temp = x_s; temp.axpy(gamma-lambda,y);
827 r_norm = temp.two_norm() * omega;
832 temp = x_s; temp.axpy(-mu_s_update,y);
833 r_norm = temp.two_norm() * omega;
839 x_s = y; x_s *= omega;
845 if (verbosity_level_ > 1)
846 std::cout << blank_ <<
"iteration " 847 << std::left << std::setw(3) << nIterations_
848 <<
" (" << (doRQI ?
"RQI," :
"II, ")
849 <<
" " << (doRQI ?
"—>" :
" ") <<
" " 850 <<
"║r║_2 = " << std::setw(11) << r_norm
851 <<
", " << (doRQI ?
" " :
"—>") <<
" " 852 <<
"║q║_2 = " << std::setw(11) << q_norm
853 <<
"): λ = " << lambda << std::endl
854 << std::resetiosflags(std::ios::left);
857 if (!doRQI && q_norm < eta)
863 assert(std::abs(mu_s-gamma) < eta);
871 if (!extrnl && doRQI && std::abs(mu_s-gamma) >= eta)
877 if (extrnl && !doRQI)
881 if (nIterations_ >= m &&
882 std::abs(mu_s - mu_s_old) / std::abs(mu_s) < delta)
895 if (verbosity_level_ > 0)
898 std::cout << blank_ <<
"Interval " 899 <<
"(" << gamma - eta <<
"," << gamma + eta
900 <<
") is free of eigenvalues, approximating " 901 <<
"the closest eigenvalue." << std::endl;
902 std::cout << blank_ <<
"Result (" 903 <<
"#iterations = " << nIterations_ <<
", " 904 <<
"║residual║_2 = " << r_norm <<
"): " 905 <<
"λ = " << lambda << std::endl;
906 if (verbosity_level_ > 2)
946 itMatrix_ = std::unique_ptr<BCRSMatrix>(
new BCRSMatrix(m_));
958 if (nIterations_ == 0)
979 template <
typename ISTLLinearSolver>
981 ISTLLinearSolver& solver)
const 984 if (mu == mu_)
return;
993 constexpr
int rowBlockSize = BCRSMatrix::block_type::rows;
994 constexpr
int colBlockSize = BCRSMatrix::block_type::cols;
996 i < itMatrix_->M()*rowBlockSize; ++i)
999 const Real& m_entry = m_
1000 [i/rowBlockSize][i/colBlockSize][i%rowBlockSize][i%colBlockSize];
1002 Real& entry = (*itMatrix_)
1003 [i/rowBlockSize][i/colBlockSize][i%rowBlockSize][i%colBlockSize];
1005 entry = m_entry - mu_;
1009 (solver,*itMatrix_);
1047 #endif // DUNE_ISTL_EIGENVALUE_POWERITERATION_HH OperatorSum IterationOperator
Type of iteration operator (m_ - mu_*I)
Definition: poweriteration.hh:186
const std::string blank_
Definition: poweriteration.hh:1040
const MatrixOperator matrixOperator_
Definition: poweriteration.hh:1025
Definition: allocator.hh:7
void updateShiftMu(const Real &mu, ISTLLinearSolver &solver) const
Update shift mu_, i.e. update iteration operator/matrix (m_ - mu_*I).
Definition: poweriteration.hh:980
Impl::ScalingLinearOperator< BlockVector > ScalingOperator
Definition: poweriteration.hh:178
const unsigned int nIterationsMax_
Definition: poweriteration.hh:1016
A sparse block matrix with compressed row storage.
Definition: bcrsmatrix.hh:422
unsigned int nIterations_
Definition: poweriteration.hh:1036
Some generic functions for pretty printing vectors and matrices.
B::field_type field_type
export the type representing the field
Definition: bvector.hh:323
const unsigned int verbosity_level_
Definition: poweriteration.hh:1019
Real mu_
Definition: poweriteration.hh:1022
A linear operator.
Definition: operators.hh:64
Implementations of the inverse operator interface.
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:688
IterationOperator & getIterationOperator()
Return the iteration operator (m_ - mu_*I).
Definition: poweriteration.hh:922
PowerIteration_Algorithms(const BCRSMatrix &m, const unsigned int nIterationsMax=1000, const unsigned int verbosity_level=0)
Construct from required parameters.
Definition: poweriteration.hh:203
const ScalingOperator scalingOperator_
Definition: poweriteration.hh:1026
virtual void apply(const X &x, Y &y) const =0
apply operator to x: The input vector is consistent and the output must also be consistent on the in...
unsigned int getIterationCount() const
Return the number of iterations in last application of an algorithm.
Definition: poweriteration.hh:956
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:391
virtual void applyscaleadd(field_type alpha, const X &x, Y &y) const =0
apply operator to x, scale and add:
Define general, extensible interface for operators. The available implementation wraps a matrix...
const std::string title_
Definition: poweriteration.hh:1039
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:352
Adapter to turn a matrix into a linear operator.
Definition: operators.hh:133
BlockVector::field_type Real
Type of underlying field.
Definition: poweriteration.hh:183
Iterative eigenvalue algorithms based on power iteration.
Definition: poweriteration.hh:172
static void setMatrix(ISTLLinearSolver &solver, const BCRSMatrix &matrix)
Definition: solver.hh:320
derive error class from the base class in common
Definition: istlexception.hh:16
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:530
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
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:942
Dune::MatrixAdapter< BCRSMatrix, BlockVector, BlockVector > MatrixOperator
Definition: poweriteration.hh:177
std::unique_ptr< BCRSMatrix > itMatrix_
Definition: poweriteration.hh:1031
A::size_type size_type
The type for the index access and the size.
Definition: bcrsmatrix.hh:457
The number of blocklevels the matrix contains.
Definition: bcrsmatrix.hh:465
Category
Definition: solvercategory.hh:21
X domain_type
The type of the domain of the operator.
Definition: operators.hh:67
A vector of blocks with memory management.
Definition: bvector.hh:316
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:257
OperatorSum itOperator_
Definition: poweriteration.hh:1027
Y range_type
The type of the range of the operator.
Definition: operators.hh:69
Statistics about the application of an inverse operator.
Definition: solver.hh:40
Category for sequential solvers.
Definition: solvercategory.hh:23
const BCRSMatrix & m_
Definition: poweriteration.hh:1015
virtual SolverCategory::Category category() const =0
Category of the linear operator (see SolverCategory::Category)
X::field_type field_type
The field type of the operator.
Definition: operators.hh:71
Templates characterizing the type of a solver.
Impl::LinearOperatorSum< MatrixOperator, ScalingOperator > OperatorSum
Definition: poweriteration.hh:179