dune-istl  2.8.0
poweriteration.hh
Go to the documentation of this file.
1 // -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
2 // vi: set et ts=4 sw=2 sts=2:
3 #ifndef DUNE_ISTL_EIGENVALUE_POWERITERATION_HH
4 #define DUNE_ISTL_EIGENVALUE_POWERITERATION_HH
5 
6 #include <cstddef> // provides std::size_t
7 #include <cmath> // provides std::sqrt, std::abs
8 
9 #include <type_traits> // provides std::is_same
10 #include <iostream> // provides std::cout, std::endl
11 #include <limits> // provides std::numeric_limits
12 #include <ios> // provides std::left, std::ios::left
13 #include <iomanip> // provides std::setw, std::resetiosflags
14 #include <memory> // provides std::unique_ptr
15 #include <string> // provides std::string
16 
17 #include <dune/common/exceptions.hh> // provides DUNE_THROW(...)
18 
19 #include <dune/istl/blocklevel.hh> // provides Dune::blockLevel
20 #include <dune/istl/operators.hh> // provides Dune::LinearOperator
21 #include <dune/istl/solvercategory.hh> // provides Dune::SolverCategory::sequential
22 #include <dune/istl/solvertype.hh> // provides Dune::IsDirectSolver
23 #include <dune/istl/operators.hh> // provides Dune::MatrixAdapter
24 #include <dune/istl/istlexception.hh> // provides Dune::ISTLError
25 #include <dune/istl/io.hh> // provides Dune::printvector(...)
26 #include <dune/istl/solvers.hh> // provides Dune::InverseOperatorResult
27 
28 namespace Dune
29 {
30 
35  namespace Impl {
43  template <class X, class Y = X>
44  class ScalingLinearOperator : public Dune::LinearOperator<X,Y>
45  {
46  public:
47  typedef X domain_type;
48  typedef Y range_type;
49  typedef typename X::field_type field_type;
50 
51  ScalingLinearOperator (field_type immutable_scaling,
52  const field_type& mutable_scaling)
53  : immutable_scaling_(immutable_scaling),
54  mutable_scaling_(mutable_scaling)
55  {}
56 
57  virtual void apply (const X& x, Y& y) const
58  {
59  y = x;
60  y *= immutable_scaling_*mutable_scaling_;
61  }
62 
63  virtual void applyscaleadd (field_type alpha, const X& x, Y& y) const
64  {
65  X temp(x);
66  temp *= immutable_scaling_*mutable_scaling_;
67  y.axpy(alpha,temp);
68  }
69 
71  virtual SolverCategory::Category category() const
72  {
74  }
75 
76  protected:
77  const field_type immutable_scaling_;
78  const field_type& mutable_scaling_;
79  };
80 
81 
90  template <class OP1, class OP2>
91  class LinearOperatorSum
92  : public Dune::LinearOperator<typename OP1::domain_type,
93  typename OP1::range_type>
94  {
95  public:
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;
99 
100  LinearOperatorSum (const OP1& op1, const OP2& op2)
101  : op1_(op1), op2_(op2)
102  {
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!");
107  }
108 
109  virtual void apply (const domain_type& x, range_type& y) const
110  {
111  op1_.apply(x,y);
112  op2_.applyscaleadd(1.0,x,y);
113  }
114 
115  virtual void applyscaleadd (field_type alpha,
116  const domain_type& x, range_type& y) const
117  {
118  range_type temp(y);
119  op1_.apply(x,temp);
120  op2_.applyscaleadd(1.0,x,temp);
121  y.axpy(alpha,temp);
122  }
123 
125  virtual SolverCategory::Category category() const
126  {
128  }
129 
130  protected:
131  const OP1& op1_;
132  const OP2& op2_;
133  };
134  } // end namespace Impl
135 
172  template <typename BCRSMatrix, typename BlockVector>
174  {
175  protected:
176  // Type definitions for type of iteration operator (m_ - mu_*I)
179  typedef Impl::ScalingLinearOperator<BlockVector> ScalingOperator;
180  typedef Impl::LinearOperatorSum<MatrixOperator,ScalingOperator> OperatorSum;
181 
182  public:
184  typedef typename BlockVector::field_type Real;
185 
188 
189  public:
205  const unsigned int nIterationsMax = 1000,
206  const unsigned int verbosity_level = 0)
207  : m_(m), nIterationsMax_(nIterationsMax),
208  verbosity_level_(verbosity_level),
209  mu_(0.0),
211  scalingOperator_(-1.0,mu_),
213  nIterations_(0),
214  title_(" PowerIteration_Algorithms: "),
215  blank_(title_.length(),' ')
216  {
217  // assert that BCRSMatrix type has blocklevel 2
218  static_assert
219  (blockLevel<BCRSMatrix>() == 2,
220  "Only BCRSMatrices with blocklevel 2 are supported.");
221 
222  // assert that BCRSMatrix type has square blocks
223  static_assert
224  (BCRSMatrix::block_type::rows == BCRSMatrix::block_type::cols,
225  "Only BCRSMatrices with square blocks are supported.");
226 
227  // assert that m_ is square
228  const int nrows = m_.M() * BCRSMatrix::block_type::rows;
229  const int ncols = m_.N() * BCRSMatrix::block_type::cols;
230  if (nrows != ncols)
231  DUNE_THROW(Dune::ISTLError,"Matrix is not square ("
232  << nrows << "x" << ncols << ").");
233  }
234 
239 
245 
258  inline void applyPowerIteration (const Real& epsilon,
259  BlockVector& x, Real& lambda) const
260  {
261  // print verbosity information
262  if (verbosity_level_ > 0)
263  std::cout << title_
264  << "Performing power iteration approximating "
265  << "the dominant eigenvalue." << std::endl;
266 
267  // allocate memory for auxiliary variables
268  BlockVector y(x);
269  BlockVector temp(x);
270 
271  // perform power iteration
272  x *= (1.0 / x.two_norm());
273  m_.mv(x,y);
274  Real r_norm = std::numeric_limits<Real>::max();
275  nIterations_ = 0;
276  while (r_norm > epsilon)
277  {
278  // update and check number of iterations
280  DUNE_THROW(Dune::ISTLError,"Power iteration did not converge "
281  << "in " << nIterationsMax_ << " iterations "
282  << "(║residual║_2 = " << r_norm << ", epsilon = "
283  << epsilon << ").");
284 
285  // do one iteration of the power iteration algorithm
286  // (use that y = m_ * x)
287  x = y;
288  x *= (1.0 / y.two_norm());
289 
290  // get approximated eigenvalue lambda via the Rayleigh quotient
291  m_.mv(x,y);
292  lambda = x * y;
293 
294  // get norm of residual (use that y = m_ * x)
295  temp = y;
296  temp.axpy(-lambda,x);
297  r_norm = temp.two_norm();
298 
299  // print verbosity information
300  if (verbosity_level_ > 1)
301  std::cout << blank_ << std::left
302  << "iteration " << std::setw(3) << nIterations_
303  << " (║residual║_2 = " << std::setw(11) << r_norm
304  << "): λ = " << lambda << std::endl
305  << std::resetiosflags(std::ios::left);
306  }
307 
308  // print verbosity information
309  if (verbosity_level_ > 0)
310  {
311  std::cout << blank_ << "Result ("
312  << "#iterations = " << nIterations_ << ", "
313  << "║residual║_2 = " << r_norm << "): "
314  << "λ = " << lambda << std::endl;
315  if (verbosity_level_ > 2)
316  {
317  // print approximated eigenvector via DUNE-ISTL I/O methods
318  Dune::printvector(std::cout,x,blank_+"x",blank_+"row");
319  }
320  }
321  }
322 
351  template <typename ISTLLinearSolver,
352  bool avoidLinSolverCrime = false>
353  inline void applyInverseIteration (const Real& epsilon,
354  ISTLLinearSolver& solver,
355  BlockVector& x, Real& lambda) const
356  {
357  constexpr Real gamma = 0.0;
358  applyInverseIteration(gamma,epsilon,solver,x,lambda);
359  }
360 
390  template <typename ISTLLinearSolver,
391  bool avoidLinSolverCrime = false>
392  inline void applyInverseIteration (const Real& gamma,
393  const Real& epsilon,
394  ISTLLinearSolver& solver,
395  BlockVector& x, Real& lambda) const
396  {
397  // print verbosity information
398  if (verbosity_level_ > 0)
399  {
400  std::cout << title_;
401  if (gamma == 0.0)
402  std::cout << "Performing inverse iteration approximating "
403  << "the least dominant eigenvalue." << std::endl;
404  else
405  std::cout << "Performing inverse iteration with shift "
406  << "gamma = " << gamma << " approximating the "
407  << "eigenvalue closest to gamma." << std::endl;
408  }
409 
410  // initialize iteration operator,
411  // initialize iteration matrix when needed
412  updateShiftMu(gamma,solver);
413 
414  // allocate memory for linear solver statistics
415  Dune::InverseOperatorResult solver_statistics;
416 
417  // allocate memory for auxiliary variables
418  BlockVector y(x);
419  Real y_norm;
420  BlockVector temp(x);
421 
422  // perform inverse iteration with shift
423  x *= (1.0 / x.two_norm());
424  Real r_norm = std::numeric_limits<Real>::max();
425  nIterations_ = 0;
426  while (r_norm > epsilon)
427  {
428  // update and check number of iterations
430  DUNE_THROW(Dune::ISTLError,"Inverse iteration "
431  << (gamma != 0.0 ? "with shift " : "") << "did not "
432  << "converge in " << nIterationsMax_ << " iterations "
433  << "(║residual║_2 = " << r_norm << ", epsilon = "
434  << epsilon << ").");
435 
436  // do one iteration of the inverse iteration with shift algorithm,
437  // part 1: solve (m_ - gamma*I) * y = x for y
438  // (protect x from being changed)
439  temp = x;
440  solver.apply(y,temp,solver_statistics);
441 
442  // get norm of y
443  y_norm = y.two_norm();
444 
445  // compile time switch between accuracy and efficiency
446  if (avoidLinSolverCrime)
447  {
448  // get approximated eigenvalue lambda via the Rayleigh quotient
449  // (use that x_new = y / y_norm)
450  m_.mv(y,temp);
451  lambda = (y * temp) / (y_norm * y_norm);
452 
453  // get norm of residual
454  // (use that x_new = y / y_norm, additionally use that temp = m_ * y)
455  temp.axpy(-lambda,y);
456  r_norm = temp.two_norm() / y_norm;
457  }
458  else
459  {
460  // get approximated eigenvalue lambda via the Rayleigh quotient
461  // (use that x_new = y / y_norm and use that (m_ - gamma*I) * y = x)
462  lambda = gamma + (y * x) / (y_norm * y_norm);
463 
464  // get norm of residual
465  // (use that x_new = y / y_norm and use that (m_ - gamma*I) * y = x)
466  temp = x; temp.axpy(gamma-lambda,y);
467  r_norm = temp.two_norm() / y_norm;
468  }
469 
470  // do one iteration of the inverse iteration with shift algorithm,
471  // part 2: update x
472  x = y;
473  x *= (1.0 / y_norm);
474 
475  // print verbosity information
476  if (verbosity_level_ > 1)
477  std::cout << blank_ << std::left
478  << "iteration " << std::setw(3) << nIterations_
479  << " (║residual║_2 = " << std::setw(11) << r_norm
480  << "): λ = " << lambda << std::endl
481  << std::resetiosflags(std::ios::left);
482  }
483 
484  // print verbosity information
485  if (verbosity_level_ > 0)
486  {
487  std::cout << blank_ << "Result ("
488  << "#iterations = " << nIterations_ << ", "
489  << "║residual║_2 = " << r_norm << "): "
490  << "λ = " << lambda << std::endl;
491  if (verbosity_level_ > 2)
492  {
493  // print approximated eigenvector via DUNE-ISTL I/O methods
494  Dune::printvector(std::cout,x,blank_+"x",blank_+"row");
495  }
496  }
497  }
498 
529  template <typename ISTLLinearSolver,
530  bool avoidLinSolverCrime = false>
531  inline void applyRayleighQuotientIteration (const Real& epsilon,
532  ISTLLinearSolver& solver,
533  BlockVector& x, Real& lambda) const
534  {
535  // print verbosity information
536  if (verbosity_level_ > 0)
537  std::cout << title_
538  << "Performing Rayleigh quotient iteration for "
539  << "estimated eigenvalue " << lambda << "." << std::endl;
540 
541  // allocate memory for linear solver statistics
542  Dune::InverseOperatorResult solver_statistics;
543 
544  // allocate memory for auxiliary variables
545  BlockVector y(x);
546  Real y_norm;
547  Real lambda_update;
548  BlockVector temp(x);
549 
550  // perform Rayleigh quotient iteration
551  x *= (1.0 / x.two_norm());
552  Real r_norm = std::numeric_limits<Real>::max();
553  nIterations_ = 0;
554  while (r_norm > epsilon)
555  {
556  // update and check number of iterations
558  DUNE_THROW(Dune::ISTLError,"Rayleigh quotient iteration did not "
559  << "converge in " << nIterationsMax_ << " iterations "
560  << "(║residual║_2 = " << r_norm << ", epsilon = "
561  << epsilon << ").");
562 
563  // update iteration operator,
564  // update iteration matrix when needed
565  updateShiftMu(lambda,solver);
566 
567  // do one iteration of the Rayleigh quotient iteration algorithm,
568  // part 1: solve (m_ - lambda*I) * y = x for y
569  // (protect x from being changed)
570  temp = x;
571  solver.apply(y,temp,solver_statistics);
572 
573  // get norm of y
574  y_norm = y.two_norm();
575 
576  // compile time switch between accuracy and efficiency
577  if (avoidLinSolverCrime)
578  {
579  // get approximated eigenvalue lambda via the Rayleigh quotient
580  // (use that x_new = y / y_norm)
581  m_.mv(y,temp);
582  lambda = (y * temp) / (y_norm * y_norm);
583 
584  // get norm of residual
585  // (use that x_new = y / y_norm, additionally use that temp = m_ * y)
586  temp.axpy(-lambda,y);
587  r_norm = temp.two_norm() / y_norm;
588  }
589  else
590  {
591  // get approximated eigenvalue lambda via the Rayleigh quotient
592  // (use that x_new = y / y_norm and use that (m_ - lambda_old*I) * y = x)
593  lambda_update = (y * x) / (y_norm * y_norm);
594  lambda += lambda_update;
595 
596  // get norm of residual
597  // (use that x_new = y / y_norm and use that (m_ - lambda_old*I) * y = x)
598  temp = x; temp.axpy(-lambda_update,y);
599  r_norm = temp.two_norm() / y_norm;
600  }
601 
602  // do one iteration of the Rayleigh quotient iteration algorithm,
603  // part 2: update x
604  x = y;
605  x *= (1.0 / y_norm);
606 
607  // print verbosity information
608  if (verbosity_level_ > 1)
609  std::cout << blank_ << std::left
610  << "iteration " << std::setw(3) << nIterations_
611  << " (║residual║_2 = " << std::setw(11) << r_norm
612  << "): λ = " << lambda << std::endl
613  << std::resetiosflags(std::ios::left);
614  }
615 
616  // print verbosity information
617  if (verbosity_level_ > 0)
618  {
619  std::cout << blank_ << "Result ("
620  << "#iterations = " << nIterations_ << ", "
621  << "║residual║_2 = " << r_norm << "): "
622  << "λ = " << lambda << std::endl;
623  if (verbosity_level_ > 2)
624  {
625  // print approximated eigenvector via DUNE-ISTL I/O methods
626  Dune::printvector(std::cout,x,blank_+"x",blank_+"row");
627  }
628  }
629  }
630 
687  template <typename ISTLLinearSolver,
688  bool avoidLinSolverCrime = false>
689  inline void applyTLIMEIteration (const Real& gamma, const Real& eta,
690  const Real& epsilon,
691  ISTLLinearSolver& solver,
692  const Real& delta, const std::size_t& m,
693  bool& extrnl,
694  BlockVector& x, Real& lambda) const
695  {
696  // use same variable names as in [Szyld, 1988]
697  BlockVector& x_s = x;
698  Real& mu_s = lambda;
699 
700  // print verbosity information
701  if (verbosity_level_ > 0)
702  std::cout << title_
703  << "Performing TLIME iteration for "
704  << "estimated eigenvalue in the "
705  << "interval (" << gamma - eta << ","
706  << gamma + eta << ")." << std::endl;
707 
708  // allocate memory for linear solver statistics
709  Dune::InverseOperatorResult solver_statistics;
710 
711  // allocate memory for auxiliary variables
712  bool doRQI;
713  Real mu;
714  BlockVector y(x_s);
715  Real omega;
716  Real mu_s_old;
717  Real mu_s_update;
718  BlockVector temp(x_s);
719  Real q_norm, r_norm;
720 
721  // perform TLIME iteration
722  x_s *= (1.0 / x_s.two_norm());
723  extrnl = true;
724  doRQI = false;
725  r_norm = std::numeric_limits<Real>::max();
726  nIterations_ = 0;
727  while (r_norm > epsilon)
728  {
729  // update and check number of iterations
731  DUNE_THROW(Dune::ISTLError,"TLIME iteration did not "
732  << "converge in " << nIterationsMax_
733  << " iterations (║residual║_2 = " << r_norm
734  << ", epsilon = " << epsilon << ").");
735 
736  // set shift for next iteration according to inverse iteration
737  // with shift (II) resp. Rayleigh quotient iteration (RQI)
738  if (doRQI)
739  mu = mu_s;
740  else
741  mu = gamma;
742 
743  // update II/RQI iteration operator,
744  // update II/RQI iteration matrix when needed
745  updateShiftMu(mu,solver);
746 
747  // do one iteration of the II/RQI algorithm,
748  // part 1: solve (m_ - mu*I) * y = x for y
749  temp = x_s;
750  solver.apply(y,temp,solver_statistics);
751 
752  // do one iteration of the II/RQI algorithm,
753  // part 2: compute omega
754  omega = (1.0 / y.two_norm());
755 
756  // backup the old Rayleigh quotient
757  mu_s_old = mu_s;
758 
759  // compile time switch between accuracy and efficiency
760  if (avoidLinSolverCrime)
761  {
762  // update the Rayleigh quotient mu_s, i.e. the approximated eigenvalue
763  // (use that x_new = y * omega)
764  m_.mv(y,temp);
765  mu_s = (y * temp) * (omega * omega);
766 
767  // get norm of "the residual with respect to the shift used by II",
768  // use normal representation of q
769  // (use that x_new = y * omega, use that temp = m_ * y)
770  temp.axpy(-gamma,y);
771  q_norm = temp.two_norm() * omega;
772 
773  // get norm of "the residual with respect to the Rayleigh quotient"
774  r_norm = q_norm*q_norm - (gamma-mu_s)*(gamma-mu_s);
775  // prevent that truncation errors invalidate the norm
776  // (we don't want to calculate sqrt of a negative number)
777  if (r_norm >= 0)
778  {
779  // use relation between the norms of r and q for efficiency
780  r_norm = std::sqrt(r_norm);
781  }
782  else
783  {
784  // use relation between r and q
785  // (use that x_new = y * omega, use that temp = (m_ - gamma*I) * y = q / omega)
786  temp.axpy(gamma-mu_s,y);
787  r_norm = temp.two_norm() * omega;
788  }
789  }
790  else
791  {
792  // update the Rayleigh quotient mu_s, i.e. the approximated eigenvalue
793  if (!doRQI)
794  {
795  // (use that x_new = y * omega, additionally use that (m_ - gamma*I) * y = x_s)
796  mu_s = gamma + (y * x_s) * (omega * omega);
797  }
798  else
799  {
800  // (use that x_new = y * omega, additionally use that (m_ - mu_s_old*I) * y = x_s)
801  mu_s_update = (y * x_s) * (omega * omega);
802  mu_s += mu_s_update;
803  }
804 
805  // get norm of "the residual with respect to the shift used by II"
806  if (!doRQI)
807  {
808  // use special representation of q in the II case
809  // (use that x_new = y * omega, additionally use that (m_ - gamma*I) * y = x_s)
810  q_norm = omega;
811  }
812  else
813  {
814  // use special representation of q in the RQI case
815  // (use that x_new = y * omega, additionally use that (m_ - mu_s_old*I) * y = x_s)
816  temp = x_s; temp.axpy(mu_s-gamma,y);
817  q_norm = temp.two_norm() * omega;
818  }
819 
820  // get norm of "the residual with respect to the Rayleigh quotient"
821  // don't use efficient relation between the norms of r and q, as
822  // this relation seems to yield a less accurate r_norm in the case
823  // where linear solver crime is admitted
824  if (!doRQI)
825  {
826  // (use that x_new = y * omega and use that (m_ - gamma*I) * y = x_s)
827  temp = x_s; temp.axpy(gamma-lambda,y);
828  r_norm = temp.two_norm() * omega;
829  }
830  else
831  {
832  // (use that x_new = y * omega and use that (m_ - mu_s_old*I) * y = x_s)
833  temp = x_s; temp.axpy(-mu_s_update,y);
834  r_norm = temp.two_norm() * omega;
835  }
836  }
837 
838  // do one iteration of the II/RQI algorithm,
839  // part 3: update x
840  x_s = y; x_s *= omega;
841 
842  // // for relative residual norm mode, scale with mu_s^{-1}
843  // r_norm /= std::abs(mu_s);
844 
845  // print verbosity information
846  if (verbosity_level_ > 1)
847  std::cout << blank_ << "iteration "
848  << std::left << std::setw(3) << nIterations_
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);
856 
857  // check if the eigenvalue closest to gamma lies in J
858  if (!doRQI && q_norm < eta)
859  {
860  // J is not free of eigenvalues
861  extrnl = false;
862 
863  // by theory we know now that mu_s also lies in J
864  assert(std::abs(mu_s-gamma) < eta);
865 
866  // switch to RQI
867  doRQI = true;
868  }
869 
870  // revert to II if J is not free of eigenvalues but
871  // at some point mu_s falls back again outside J
872  if (!extrnl && doRQI && std::abs(mu_s-gamma) >= eta)
873  doRQI = false;
874 
875  // if eigenvalue closest to gamma does not lie in J use RQI
876  // solely to accelerate the convergence to this eigenvalue
877  // when II has become stationary
878  if (extrnl && !doRQI)
879  {
880  // switch to RQI if the relative change of the Rayleigh
881  // quotient indicates that II has become stationary
882  if (nIterations_ >= m &&
883  std::abs(mu_s - mu_s_old) / std::abs(mu_s) < delta)
884  doRQI = true;
885  }
886  }
887 
888  // // compute final residual and lambda again (paranoia....)
889  // m_.mv(x_s,temp);
890  // mu_s = x_s * temp;
891  // temp.axpy(-mu_s,x_s);
892  // r_norm = temp.two_norm();
893  // // r_norm /= std::abs(mu_s);
894 
895  // print verbosity information
896  if (verbosity_level_ > 0)
897  {
898  if (extrnl)
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 ("
904  << "#iterations = " << nIterations_ << ", "
905  << "║residual║_2 = " << r_norm << "): "
906  << "λ = " << lambda << std::endl;
907  if (verbosity_level_ > 2)
908  {
909  // print approximated eigenvector via DUNE-ISTL I/O methods
910  Dune::printvector(std::cout,x,blank_+"x",blank_+"row");
911  }
912  }
913  }
914 
924  {
925  // return iteration operator
926  return itOperator_;
927  }
928 
943  inline const BCRSMatrix& getIterationMatrix () const
944  {
945  // create iteration matrix on demand
946  if (!itMatrix_)
947  itMatrix_ = std::make_unique<BCRSMatrix>(m_);
948 
949  // return iteration matrix
950  return *itMatrix_;
951  }
952 
957  inline unsigned int getIterationCount () const
958  {
959  if (nIterations_ == 0)
960  DUNE_THROW(Dune::ISTLError,"No algorithm applied, yet.");
961 
962  return nIterations_;
963  }
964 
965  protected:
980  template <typename ISTLLinearSolver>
981  inline void updateShiftMu (const Real& mu,
982  ISTLLinearSolver& solver) const
983  {
984  // do nothing if new shift equals the old one
985  if (mu == mu_) return;
986 
987  // update shift mu_, i.e. update iteration operator
988  mu_ = mu;
989 
990  // update iteration matrix when needed
991  if (itMatrix_)
992  {
993  // iterate over entries in iteration matrix diagonal
994  constexpr int rowBlockSize = BCRSMatrix::block_type::rows;
995  constexpr int colBlockSize = BCRSMatrix::block_type::cols;
996  for (typename BCRSMatrix::size_type i = 0;
997  i < itMatrix_->M()*rowBlockSize; ++i)
998  {
999  // access m_[i,i] where i is the flat index of a row/column
1000  const Real& m_entry = m_
1001  [i/rowBlockSize][i/colBlockSize][i%rowBlockSize][i%colBlockSize];
1002  // access *itMatrix[i,i] where i is the flat index of a row/column
1003  Real& entry = (*itMatrix_)
1004  [i/rowBlockSize][i/colBlockSize][i%rowBlockSize][i%colBlockSize];
1005  // change current entry in iteration matrix diagonal
1006  entry = m_entry - mu_;
1007  }
1008  // notify linear solver about change of the iteration matrix object
1010  (solver,*itMatrix_);
1011  }
1012  }
1013 
1014  protected:
1015  // parameters related to iterative eigenvalue algorithms
1016  const BCRSMatrix& m_;
1017  const unsigned int nIterationsMax_;
1018 
1019  // verbosity setting
1020  const unsigned int verbosity_level_;
1021 
1022  // shift mu_ used by iteration operator/matrix (m_ - mu_*I)
1023  mutable Real mu_;
1024 
1025  // iteration operator (m_ - mu_*I), passing shift mu_ by reference
1029 
1030  // iteration matrix (m_ - mu_*I), provided on demand when needed
1031  // (e.g. for preconditioning)
1032  mutable std::unique_ptr<BCRSMatrix> itMatrix_;
1033 
1034  // memory for storing temporary variables (mutable as they shall
1035  // just be effectless auxiliary variables of the const apply*(...)
1036  // methods)
1037  mutable unsigned int nIterations_;
1038 
1039  // constants for printing verbosity information
1040  const std::string title_;
1041  const std::string blank_;
1042  };
1043 
1046 } // namespace Dune
1047 
1048 #endif // DUNE_ISTL_EIGENVALUE_POWERITERATION_HH
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