Dune Core Modules (2.5.0)

poweriteration.hh
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
28namespace Dune
29{
30
38 template <class X, class Y = X>
40 {
41 public:
42 typedef X domain_type;
43 typedef Y range_type;
44 typedef typename X::field_type field_type;
45
46 enum {category = Dune::SolverCategory::sequential};
47
48 ScalingLinearOperator (field_type immutable_scaling,
49 const field_type& mutable_scaling)
50 : immutable_scaling_(immutable_scaling),
51 mutable_scaling_(mutable_scaling)
52 {}
53
54 virtual void apply (const X& x, Y& y) const
55 {
56 y = x;
57 y *= immutable_scaling_*mutable_scaling_;
58 }
59
60 virtual void applyscaleadd (field_type alpha, const X& x, Y& y) const
61 {
62 X temp(x);
63 temp *= immutable_scaling_*mutable_scaling_;
64 y.axpy(alpha,temp);
65 }
66
67 protected:
68 const field_type immutable_scaling_;
69 const field_type& mutable_scaling_;
70 };
71
72
81 template <class OP1, class OP2>
83 : public Dune::LinearOperator<typename OP1::domain_type,
84 typename OP1::range_type>
85 {
86 public:
87 typedef typename OP1::domain_type domain_type;
88 typedef typename OP1::range_type range_type;
89 typedef typename domain_type::field_type field_type;
90
91 enum {category = Dune::SolverCategory::sequential};
92
93 LinearOperatorSum (const OP1& op1, const OP2& op2)
94 : op1_(op1), op2_(op2)
95 {
96 static_assert(std::is_same<typename OP2::domain_type,domain_type>::value,
97 "Domain type of both operators doesn't match!");
98 static_assert(std::is_same<typename OP2::range_type,range_type>::value,
99 "Range type of both operators doesn't match!");
100 }
101
102 virtual void apply (const domain_type& x, range_type& y) const
103 {
104 op1_.apply(x,y);
105 op2_.applyscaleadd(1.0,x,y);
106 }
107
108 virtual void applyscaleadd (field_type alpha,
109 const domain_type& x, range_type& y) const
110 {
111 range_type temp(y);
112 op1_.apply(x,temp);
113 op2_.applyscaleadd(1.0,x,temp);
114 y.axpy(alpha,temp);
115 }
116
117 protected:
118 const OP1& op1_;
119 const OP2& op2_;
120 };
121
122
130 template <typename ISTLLinearSolver, typename BCRSMatrix>
132 {
133 public:
134 static void setMatrix (ISTLLinearSolver& solver,
135 const BCRSMatrix& matrix)
136 {
137 static const bool is_direct_solver
138 = Dune::IsDirectSolver<ISTLLinearSolver>::value;
141 }
142
143 protected:
148 template <bool is_direct_solver, typename Dummy = void>
150 {
151 static void setMatrix (ISTLLinearSolver&,
152 const BCRSMatrix&)
153 {}
154 };
155
160 template <typename Dummy>
161 struct Implementation<true,Dummy>
162 {
163 static void setMatrix (ISTLLinearSolver& solver,
164 const BCRSMatrix& matrix)
165 {
166 solver.setMatrix(matrix);
167 }
168 };
169 };
170
171
209 template <typename BCRSMatrix, typename BlockVector>
211 {
212 protected:
213 // Type definitions for type of iteration operator (m_ - mu_*I)
218
219 public:
222
225
226 public:
242 const unsigned int nIterationsMax = 1000,
243 const unsigned int verbosity_level = 0)
244 : m_(m), nIterationsMax_(nIterationsMax),
245 verbosity_level_(verbosity_level),
246 mu_(0.0),
247 matrixOperator_(m_),
248 scalingOperator_(-1.0,mu_),
249 itOperator_(matrixOperator_,scalingOperator_),
250 nIterations_(0),
251 title_(" PowerIteration_Algorithms: "),
252 blank_(title_.length(),' ')
253 {
254 // assert that BCRSMatrix type has blocklevel 2
255 static_assert
257 "Only BCRSMatrices with blocklevel 2 are supported.");
258
259 // assert that BCRSMatrix type has square blocks
260 static_assert
261 (BCRSMatrix::block_type::rows == BCRSMatrix::block_type::cols,
262 "Only BCRSMatrices with square blocks are supported.");
263
264 // assert that m_ is square
265 const int nrows = m_.M() * BCRSMatrix::block_type::rows;
266 const int ncols = m_.N() * BCRSMatrix::block_type::cols;
267 if (nrows != ncols)
268 DUNE_THROW(Dune::ISTLError,"Matrix is not square ("
269 << nrows << "x" << ncols << ").");
270 }
271
276
282
295 inline void applyPowerIteration (const Real& epsilon,
296 BlockVector& x, Real& lambda) const
297 {
298 // print verbosity information
299 if (verbosity_level_ > 0)
300 std::cout << title_
301 << "Performing power iteration approximating "
302 << "the dominant eigenvalue." << std::endl;
303
304 // allocate memory for auxiliary variables
305 BlockVector y(x);
306 BlockVector temp(x);
307
308 // perform power iteration
309 x *= (1.0 / x.two_norm());
310 m_.mv(x,y);
311 Real r_norm = std::numeric_limits<Real>::max();
312 nIterations_ = 0;
313 while (r_norm > epsilon)
314 {
315 // update and check number of iterations
316 if (++nIterations_ > nIterationsMax_)
317 DUNE_THROW(Dune::ISTLError,"Power iteration did not converge "
318 << "in " << nIterationsMax_ << " iterations "
319 << "(║residual║_2 = " << r_norm << ", epsilon = "
320 << epsilon << ").");
321
322 // do one iteration of the power iteration algorithm
323 // (use that y = m_ * x)
324 x = y;
325 x *= (1.0 / y.two_norm());
326
327 // get approximated eigenvalue lambda via the Rayleigh quotient
328 m_.mv(x,y);
329 lambda = x * y;
330
331 // get norm of residual (use that y = m_ * x)
332 temp = y;
333 temp.axpy(-lambda,x);
334 r_norm = temp.two_norm();
335
336 // print verbosity information
337 if (verbosity_level_ > 1)
338 std::cout << blank_ << std::left
339 << "iteration " << std::setw(3) << nIterations_
340 << " (║residual║_2 = " << std::setw(11) << r_norm
341 << "): λ = " << lambda << std::endl
342 << std::resetiosflags(std::ios::left);
343 }
344
345 // print verbosity information
346 if (verbosity_level_ > 0)
347 {
348 std::cout << blank_ << "Result ("
349 << "#iterations = " << nIterations_ << ", "
350 << "║residual║_2 = " << r_norm << "): "
351 << "λ = " << lambda << std::endl;
352 if (verbosity_level_ > 2)
353 {
354 // print approximated eigenvector via DUNE-ISTL I/O methods
355 Dune::printvector(std::cout,x,blank_+"x",blank_+"row");
356 }
357 }
358 }
359
388 template <typename ISTLLinearSolver,
389 bool avoidLinSolverCrime = false>
390 inline void applyInverseIteration (const Real& epsilon,
391 ISTLLinearSolver& solver,
392 BlockVector& x, Real& lambda) const
393 {
394 constexpr Real gamma = 0.0;
395 applyInverseIteration(gamma,epsilon,solver,x,lambda);
396 }
397
427 template <typename ISTLLinearSolver,
428 bool avoidLinSolverCrime = false>
429 inline void applyInverseIteration (const Real& gamma,
430 const Real& epsilon,
431 ISTLLinearSolver& solver,
432 BlockVector& x, Real& lambda) const
433 {
434 // print verbosity information
435 if (verbosity_level_ > 0)
436 {
437 std::cout << title_;
438 if (gamma == 0.0)
439 std::cout << "Performing inverse iteration approximating "
440 << "the least dominant eigenvalue." << std::endl;
441 else
442 std::cout << "Performing inverse iteration with shift "
443 << "gamma = " << gamma << " approximating the "
444 << "eigenvalue closest to gamma." << std::endl;
445 }
446
447 // initialize iteration operator,
448 // initialize iteration matrix when needed
449 updateShiftMu(gamma,solver);
450
451 // allocate memory for linear solver statistics
452 Dune::InverseOperatorResult solver_statistics;
453
454 // allocate memory for auxiliary variables
455 BlockVector y(x);
456 Real y_norm;
457 BlockVector temp(x);
458
459 // perform inverse iteration with shift
460 x *= (1.0 / x.two_norm());
461 Real r_norm = std::numeric_limits<Real>::max();
462 nIterations_ = 0;
463 while (r_norm > epsilon)
464 {
465 // update and check number of iterations
466 if (++nIterations_ > nIterationsMax_)
467 DUNE_THROW(Dune::ISTLError,"Inverse iteration "
468 << (gamma != 0.0 ? "with shift " : "") << "did not "
469 << "converge in " << nIterationsMax_ << " iterations "
470 << "(║residual║_2 = " << r_norm << ", epsilon = "
471 << epsilon << ").");
472
473 // do one iteration of the inverse iteration with shift algorithm,
474 // part 1: solve (m_ - gamma*I) * y = x for y
475 // (protect x from being changed)
476 temp = x;
477 solver.apply(y,temp,solver_statistics);
478
479 // get norm of y
480 y_norm = y.two_norm();
481
482 // compile time switch between accuracy and efficiency
483 if (avoidLinSolverCrime)
484 {
485 // get approximated eigenvalue lambda via the Rayleigh quotient
486 // (use that x_new = y / y_norm)
487 m_.mv(y,temp);
488 lambda = (y * temp) / (y_norm * y_norm);
489
490 // get norm of residual
491 // (use that x_new = y / y_norm, additionally use that temp = m_ * y)
492 temp.axpy(-lambda,y);
493 r_norm = temp.two_norm() / y_norm;
494 }
495 else
496 {
497 // get approximated eigenvalue lambda via the Rayleigh quotient
498 // (use that x_new = y / y_norm and use that (m_ - gamma*I) * y = x)
499 lambda = gamma + (y * x) / (y_norm * y_norm);
500
501 // get norm of residual
502 // (use that x_new = y / y_norm and use that (m_ - gamma*I) * y = x)
503 temp = x; temp.axpy(gamma-lambda,y);
504 r_norm = temp.two_norm() / y_norm;
505 }
506
507 // do one iteration of the inverse iteration with shift algorithm,
508 // part 2: update x
509 x = y;
510 x *= (1.0 / y_norm);
511
512 // print verbosity information
513 if (verbosity_level_ > 1)
514 std::cout << blank_ << std::left
515 << "iteration " << std::setw(3) << nIterations_
516 << " (║residual║_2 = " << std::setw(11) << r_norm
517 << "): λ = " << lambda << std::endl
518 << std::resetiosflags(std::ios::left);
519 }
520
521 // print verbosity information
522 if (verbosity_level_ > 0)
523 {
524 std::cout << blank_ << "Result ("
525 << "#iterations = " << nIterations_ << ", "
526 << "║residual║_2 = " << r_norm << "): "
527 << "λ = " << lambda << std::endl;
528 if (verbosity_level_ > 2)
529 {
530 // print approximated eigenvector via DUNE-ISTL I/O methods
531 Dune::printvector(std::cout,x,blank_+"x",blank_+"row");
532 }
533 }
534 }
535
566 template <typename ISTLLinearSolver,
567 bool avoidLinSolverCrime = false>
568 inline void applyRayleighQuotientIteration (const Real& epsilon,
569 ISTLLinearSolver& solver,
570 BlockVector& x, Real& lambda) const
571 {
572 // print verbosity information
573 if (verbosity_level_ > 0)
574 std::cout << title_
575 << "Performing Rayleigh quotient iteration for "
576 << "estimated eigenvalue " << lambda << "." << std::endl;
577
578 // allocate memory for linear solver statistics
579 Dune::InverseOperatorResult solver_statistics;
580
581 // allocate memory for auxiliary variables
582 BlockVector y(x);
583 Real y_norm;
584 Real lambda_update;
585 BlockVector temp(x);
586
587 // perform Rayleigh quotient iteration
588 x *= (1.0 / x.two_norm());
589 Real r_norm = std::numeric_limits<Real>::max();
590 nIterations_ = 0;
591 while (r_norm > epsilon)
592 {
593 // update and check number of iterations
594 if (++nIterations_ > nIterationsMax_)
595 DUNE_THROW(Dune::ISTLError,"Rayleigh quotient iteration did not "
596 << "converge in " << nIterationsMax_ << " iterations "
597 << "(║residual║_2 = " << r_norm << ", epsilon = "
598 << epsilon << ").");
599
600 // update iteration operator,
601 // update iteration matrix when needed
602 updateShiftMu(lambda,solver);
603
604 // do one iteration of the Rayleigh quotient iteration algorithm,
605 // part 1: solve (m_ - lambda*I) * y = x for y
606 // (protect x from being changed)
607 temp = x;
608 solver.apply(y,temp,solver_statistics);
609
610 // get norm of y
611 y_norm = y.two_norm();
612
613 // compile time switch between accuracy and efficiency
614 if (avoidLinSolverCrime)
615 {
616 // get approximated eigenvalue lambda via the Rayleigh quotient
617 // (use that x_new = y / y_norm)
618 m_.mv(y,temp);
619 lambda = (y * temp) / (y_norm * y_norm);
620
621 // get norm of residual
622 // (use that x_new = y / y_norm, additionally use that temp = m_ * y)
623 temp.axpy(-lambda,y);
624 r_norm = temp.two_norm() / y_norm;
625 }
626 else
627 {
628 // get approximated eigenvalue lambda via the Rayleigh quotient
629 // (use that x_new = y / y_norm and use that (m_ - lambda_old*I) * y = x)
630 lambda_update = (y * x) / (y_norm * y_norm);
631 lambda += lambda_update;
632
633 // get norm of residual
634 // (use that x_new = y / y_norm and use that (m_ - lambda_old*I) * y = x)
635 temp = x; temp.axpy(-lambda_update,y);
636 r_norm = temp.two_norm() / y_norm;
637 }
638
639 // do one iteration of the Rayleigh quotient iteration algorithm,
640 // part 2: update x
641 x = y;
642 x *= (1.0 / y_norm);
643
644 // print verbosity information
645 if (verbosity_level_ > 1)
646 std::cout << blank_ << std::left
647 << "iteration " << std::setw(3) << nIterations_
648 << " (║residual║_2 = " << std::setw(11) << r_norm
649 << "): λ = " << lambda << std::endl
650 << std::resetiosflags(std::ios::left);
651 }
652
653 // print verbosity information
654 if (verbosity_level_ > 0)
655 {
656 std::cout << blank_ << "Result ("
657 << "#iterations = " << nIterations_ << ", "
658 << "║residual║_2 = " << r_norm << "): "
659 << "λ = " << lambda << std::endl;
660 if (verbosity_level_ > 2)
661 {
662 // print approximated eigenvector via DUNE-ISTL I/O methods
663 Dune::printvector(std::cout,x,blank_+"x",blank_+"row");
664 }
665 }
666 }
667
724 template <typename ISTLLinearSolver,
725 bool avoidLinSolverCrime = false>
726 inline void applyTLIMEIteration (const Real& gamma, const Real& eta,
727 const Real& epsilon,
728 ISTLLinearSolver& solver,
729 const Real& delta, const std::size_t& m,
730 bool& extrnl,
731 BlockVector& x, Real& lambda) const
732 {
733 // use same variable names as in [Szyld, 1988]
734 BlockVector& x_s = x;
735 Real& mu_s = lambda;
736
737 // print verbosity information
738 if (verbosity_level_ > 0)
739 std::cout << title_
740 << "Performing TLIME iteration for "
741 << "estimated eigenvalue in the "
742 << "interval (" << gamma - eta << ","
743 << gamma + eta << ")." << std::endl;
744
745 // allocate memory for linear solver statistics
746 Dune::InverseOperatorResult solver_statistics;
747
748 // allocate memory for auxiliary variables
749 bool doRQI;
750 Real mu;
751 BlockVector y(x_s);
752 Real omega;
753 Real mu_s_old;
754 Real mu_s_update;
755 BlockVector temp(x_s);
756 Real q_norm, r_norm;
757
758 // perform TLIME iteration
759 x_s *= (1.0 / x_s.two_norm());
760 extrnl = true;
761 doRQI = false;
762 r_norm = std::numeric_limits<Real>::max();
763 nIterations_ = 0;
764 while (r_norm > epsilon)
765 {
766 // update and check number of iterations
767 if (++nIterations_ > nIterationsMax_)
768 DUNE_THROW(Dune::ISTLError,"TLIME iteration did not "
769 << "converge in " << nIterationsMax_
770 << " iterations (║residual║_2 = " << r_norm
771 << ", epsilon = " << epsilon << ").");
772
773 // set shift for next iteration according to inverse iteration
774 // with shift (II) resp. Rayleigh quotient iteration (RQI)
775 if (doRQI)
776 mu = mu_s;
777 else
778 mu = gamma;
779
780 // update II/RQI iteration operator,
781 // update II/RQI iteration matrix when needed
782 updateShiftMu(mu,solver);
783
784 // do one iteration of the II/RQI algorithm,
785 // part 1: solve (m_ - mu*I) * y = x for y
786 temp = x_s;
787 solver.apply(y,temp,solver_statistics);
788
789 // do one iteration of the II/RQI algorithm,
790 // part 2: compute omega
791 omega = (1.0 / y.two_norm());
792
793 // backup the old Rayleigh quotient
794 mu_s_old = mu_s;
795
796 // compile time switch between accuracy and efficiency
797 if (avoidLinSolverCrime)
798 {
799 // update the Rayleigh quotient mu_s, i.e. the approximated eigenvalue
800 // (use that x_new = y * omega)
801 m_.mv(y,temp);
802 mu_s = (y * temp) * (omega * omega);
803
804 // get norm of "the residual with respect to the shift used by II",
805 // use normal representation of q
806 // (use that x_new = y * omega, use that temp = m_ * y)
807 temp.axpy(-gamma,y);
808 q_norm = temp.two_norm() * omega;
809
810 // get norm of "the residual with respect to the Rayleigh quotient"
811 r_norm = q_norm*q_norm - (gamma-mu_s)*(gamma-mu_s);
812 // prevent that truncation errors invalidate the norm
813 // (we don't want to calculate sqrt of a negative number)
814 if (r_norm >= 0)
815 {
816 // use relation between the norms of r and q for efficiency
817 r_norm = std::sqrt(r_norm);
818 }
819 else
820 {
821 // use relation between r and q
822 // (use that x_new = y * omega, use that temp = (m_ - gamma*I) * y = q / omega)
823 temp.axpy(gamma-mu_s,y);
824 r_norm = temp.two_norm() * omega;
825 }
826 }
827 else
828 {
829 // update the Rayleigh quotient mu_s, i.e. the approximated eigenvalue
830 if (!doRQI)
831 {
832 // (use that x_new = y * omega, additionally use that (m_ - gamma*I) * y = x_s)
833 mu_s = gamma + (y * x_s) * (omega * omega);
834 }
835 else
836 {
837 // (use that x_new = y * omega, additionally use that (m_ - mu_s_old*I) * y = x_s)
838 mu_s_update = (y * x_s) * (omega * omega);
839 mu_s += mu_s_update;
840 }
841
842 // get norm of "the residual with respect to the shift used by II"
843 if (!doRQI)
844 {
845 // use special representation of q in the II case
846 // (use that x_new = y * omega, additionally use that (m_ - gamma*I) * y = x_s)
847 q_norm = omega;
848 }
849 else
850 {
851 // use special representation of q in the RQI case
852 // (use that x_new = y * omega, additionally use that (m_ - mu_s_old*I) * y = x_s)
853 temp = x_s; temp.axpy(mu_s-gamma,y);
854 q_norm = temp.two_norm() * omega;
855 }
856
857 // get norm of "the residual with respect to the Rayleigh quotient"
858 // don't use efficient relation between the norms of r and q, as
859 // this relation seems to yield a less accurate r_norm in the case
860 // where linear solver crime is admitted
861 if (!doRQI)
862 {
863 // (use that x_new = y * omega and use that (m_ - gamma*I) * y = x_s)
864 temp = x_s; temp.axpy(gamma-lambda,y);
865 r_norm = temp.two_norm() * omega;
866 }
867 else
868 {
869 // (use that x_new = y * omega and use that (m_ - mu_s_old*I) * y = x_s)
870 temp = x_s; temp.axpy(-mu_s_update,y);
871 r_norm = temp.two_norm() * omega;
872 }
873 }
874
875 // do one iteration of the II/RQI algorithm,
876 // part 3: update x
877 x_s = y; x_s *= omega;
878
879 // // for relative residual norm mode, scale with mu_s^{-1}
880 // r_norm /= std::abs(mu_s);
881
882 // print verbosity information
883 if (verbosity_level_ > 1)
884 std::cout << blank_ << "iteration "
885 << std::left << std::setw(3) << nIterations_
886 << " (" << (doRQI ? "RQI," : "II, ")
887 << " " << (doRQI ? "—>" : " ") << " "
888 << "║r║_2 = " << std::setw(11) << r_norm
889 << ", " << (doRQI ? " " : "—>") << " "
890 << "║q║_2 = " << std::setw(11) << q_norm
891 << "): λ = " << lambda << std::endl
892 << std::resetiosflags(std::ios::left);
893
894 // check if the eigenvalue closest to gamma lies in J
895 if (!doRQI && q_norm < eta)
896 {
897 // J is not free of eigenvalues
898 extrnl = false;
899
900 // by theory we know now that mu_s also lies in J
901 assert(std::abs(mu_s-gamma) < eta);
902
903 // switch to RQI
904 doRQI = true;
905 }
906
907 // revert to II if J is not free of eigenvalues but
908 // at some point mu_s falls back again outside J
909 if (!extrnl && doRQI && std::abs(mu_s-gamma) >= eta)
910 doRQI = false;
911
912 // if eigenvalue closest to gamma does not lie in J use RQI
913 // solely to accelerate the convergence to this eigenvalue
914 // when II has become stationary
915 if (extrnl && !doRQI)
916 {
917 // switch to RQI if the relative change of the Rayleigh
918 // quotient indicates that II has become stationary
919 if (nIterations_ >= m &&
920 std::abs(mu_s - mu_s_old) / std::abs(mu_s) < delta)
921 doRQI = true;
922 }
923 }
924
925 // // compute final residual and lambda again (paranoia....)
926 // m_.mv(x_s,temp);
927 // mu_s = x_s * temp;
928 // temp.axpy(-mu_s,x_s);
929 // r_norm = temp.two_norm();
930 // // r_norm /= std::abs(mu_s);
931
932 // print verbosity information
933 if (verbosity_level_ > 0)
934 {
935 if (extrnl)
936 std::cout << blank_ << "Interval "
937 << "(" << gamma - eta << "," << gamma + eta
938 << ") is free of eigenvalues, approximating "
939 << "the closest eigenvalue." << std::endl;
940 std::cout << blank_ << "Result ("
941 << "#iterations = " << nIterations_ << ", "
942 << "║residual║_2 = " << r_norm << "): "
943 << "λ = " << lambda << std::endl;
944 if (verbosity_level_ > 2)
945 {
946 // print approximated eigenvector via DUNE-ISTL I/O methods
947 Dune::printvector(std::cout,x,blank_+"x",blank_+"row");
948 }
949 }
950 }
951
961 {
962 // return iteration operator
963 return itOperator_;
964 }
965
980 inline const BCRSMatrix& getIterationMatrix () const
981 {
982 // create iteration matrix on demand
983 if (!itMatrix_)
984 itMatrix_ = std::unique_ptr<BCRSMatrix>(new BCRSMatrix(m_));
985
986 // return iteration matrix
987 return *itMatrix_;
988 }
989
994 inline unsigned int getIterationCount () const
995 {
996 if (nIterations_ == 0)
997 DUNE_THROW(Dune::ISTLError,"No algorithm applied, yet.");
998
999 return nIterations_;
1000 }
1001
1002 protected:
1017 template <typename ISTLLinearSolver>
1018 inline void updateShiftMu (const Real& mu,
1019 ISTLLinearSolver& solver) const
1020 {
1021 // do nothing if new shift equals the old one
1022 if (mu == mu_) return;
1023
1024 // update shift mu_, i.e. update iteration operator
1025 mu_ = mu;
1026
1027 // update iteration matrix when needed
1028 if (itMatrix_)
1029 {
1030 // iterate over entries in iteration matrix diagonal
1031 constexpr int rowBlockSize = BCRSMatrix::block_type::rows;
1032 constexpr int colBlockSize = BCRSMatrix::block_type::cols;
1033 for (typename BCRSMatrix::size_type i = 0;
1034 i < itMatrix_->M()*rowBlockSize; ++i)
1035 {
1036 // access m_[i,i] where i is the flat index of a row/column
1037 const Real& m_entry = m_
1038 [i/rowBlockSize][i/colBlockSize][i%rowBlockSize][i%colBlockSize];
1039 // access *itMatrix[i,i] where i is the flat index of a row/column
1040 Real& entry = (*itMatrix_)
1041 [i/rowBlockSize][i/colBlockSize][i%rowBlockSize][i%colBlockSize];
1042 // change current entry in iteration matrix diagonal
1043 entry = m_entry - mu_;
1044 }
1045 // notify linear solver about change of the iteration matrix object
1047 (solver,*itMatrix_);
1048 }
1049 }
1050
1051 protected:
1052 // parameters related to iterative eigenvalue algorithms
1053 const BCRSMatrix& m_;
1054 const unsigned int nIterationsMax_;
1055
1056 // verbosity setting
1057 const unsigned int verbosity_level_;
1058
1059 // shift mu_ used by iteration operator/matrix (m_ - mu_*I)
1060 mutable Real mu_;
1061
1062 // iteration operator (m_ - mu_*I), passing shift mu_ by reference
1063 const MatrixOperator matrixOperator_;
1064 const ScalingOperator scalingOperator_;
1065 OperatorSum itOperator_;
1066
1067 // iteration matrix (m_ - mu_*I), provided on demand when needed
1068 // (e.g. for preconditioning)
1069 mutable std::unique_ptr<BCRSMatrix> itMatrix_;
1070
1071 // memory for storing temporary variables (mutable as they shall
1072 // just be effectless auxiliary variables of the const apply*(...)
1073 // methods)
1074 mutable unsigned int nIterations_;
1075
1076 // constants for printing verbosity information
1077 const std::string title_;
1078 const std::string blank_;
1079 };
1080
1081} // namespace Dune
1082
1083
1084#endif // DUNE_ISTL_EIGENVALUE_POWERITERATION_HH
A sparse block matrix with compressed row storage.
Definition: bcrsmatrix.hh:423
@ blocklevel
The number of blocklevels the matrix contains.
Definition: bcrsmatrix.hh:465
A::size_type size_type
The type for the index access and the size.
Definition: bcrsmatrix.hh:457
size_type M() const
number of columns (counted in blocks)
Definition: bcrsmatrix.hh:1909
void mv(const X &x, Y &y) const
y = A x
Definition: bcrsmatrix.hh:1580
size_type N() const
number of rows (counted in blocks)
Definition: bcrsmatrix.hh:1903
A vector of blocks with memory management.
Definition: bvector.hh:313
B::field_type field_type
export the type representing the field
Definition: bvector.hh:319
derive error class from the base class in common
Definition: istlexception.hh:16
A linear operator representing the sum of two linear operators.
Definition: poweriteration.hh:85
virtual void apply(const domain_type &x, range_type &y) const
apply operator to x: The input vector is consistent and the output must also be consistent on the in...
Definition: poweriteration.hh:102
A linear operator.
Definition: operators.hh:62
Adapter to turn a matrix into a linear operator.
Definition: operators.hh:122
A class template for performing some iterative eigenvalue algorithms based on power iteration.
Definition: poweriteration.hh:211
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:390
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:726
IterationOperator & getIterationOperator()
Return the iteration operator (m_ - mu_*I).
Definition: poweriteration.hh:960
PowerIteration_Algorithms(const BCRSMatrix &m, const unsigned int nIterationsMax=1000, const unsigned int verbosity_level=0)
Construct from required parameters.
Definition: poweriteration.hh:241
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:295
OperatorSum IterationOperator
Type of iteration operator (m_ - mu_*I)
Definition: poweriteration.hh:224
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:568
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:429
void updateShiftMu(const Real &mu, ISTLLinearSolver &solver) const
Update shift mu_, i.e. update iteration operator/matrix (m_ - mu_*I).
Definition: poweriteration.hh:1018
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:994
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:980
BlockVector::field_type Real
Type of underlying field.
Definition: poweriteration.hh:221
A linear operator scaling vectors by a scalar value. The scalar value can be changed as it is given i...
Definition: poweriteration.hh:40
Helper class for notifying a DUNE-ISTL linear solver about a change of the iteration matrix object in...
Definition: poweriteration.hh:132
FieldTraits< field_type >::real_type two_norm() const
two norm sqrt(sum over squared values of entries)
Definition: bvector.hh:193
block_vector_unmanaged & axpy(const field_type &a, const block_vector_unmanaged &y)
vector space axpy operation
Definition: bvector.hh:124
A few common exception classes.
#define DUNE_THROW(E, m)
Definition: exceptions.hh:216
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
Some generic functions for pretty printing vectors and matrices.
Dune namespace.
Definition: alignment.hh:11
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:32
@ sequential
Category for sequential solvers.
Definition: solvercategory.hh:21
Implementation that works together with iterative ISTL solvers, e.g. Dune::CGSolver or Dune::BiCGSTAB...
Definition: poweriteration.hh:150
Creative Commons License   |  Legal Statements / Impressum  |  Hosted by TU Dresden  |  generated with Hugo v0.111.3 (Nov 23, 23:29, 2024)