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