Dune Core Modules (2.6.0)

solvers.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
4#ifndef DUNE_ISTL_SOLVERS_HH
5#define DUNE_ISTL_SOLVERS_HH
6
7#include <array>
8#include <cmath>
9#include <complex>
10#include <iostream>
11#include <memory>
12#include <type_traits>
13#include <vector>
14
15#include <dune/common/conditional.hh>
17#include <dune/common/hybridutilities.hh>
19#include <dune/common/std/type_traits.hh>
20#include <dune/common/timer.hh>
21
22#include <dune/istl/allocator.hh>
24#include <dune/istl/eigenvalue/arpackpp.hh>
25#include <dune/istl/istlexception.hh>
27#include <dune/istl/preconditioner.hh>
29#include <dune/istl/solver.hh>
30
31namespace Dune {
43 //=====================================================================
44 // Implementation of this interface
45 //=====================================================================
46
55 template<class X>
56 class LoopSolver : public IterativeSolver<X,X> {
57 public:
62
63 // copy base class constructors
65
66 // don't shadow four-argument version of apply defined in the base class
67 using IterativeSolver<X,X>::apply;
68
70 virtual void apply (X& x, X& b, InverseOperatorResult& res)
71 {
72 // clear solver statistics
73 res.clear();
74
75 // start a timer
76 Timer watch;
77
78 // prepare preconditioner
79 _prec->pre(x,b);
80
81 // overwrite b with defect
82 _op->applyscaleadd(-1,x,b);
83
84 // compute norm, \todo parallelization
85 real_type def0 = _sp->norm(b);
86
87 // printing
88 if (_verbose>0)
89 {
90 std::cout << "=== LoopSolver" << std::endl;
91 if (_verbose>1)
92 {
93 this->printHeader(std::cout);
94 this->printOutput(std::cout,0,def0);
95 }
96 }
97
98 // allocate correction vector
99 X v(x);
100
101 // iteration loop
102 int i=1; real_type def=def0;
103 for ( ; i<=_maxit; i++ )
104 {
105 v = 0; // clear correction
106 _prec->apply(v,b); // apply preconditioner
107 x += v; // update solution
108 _op->applyscaleadd(-1,v,b); // update defect
109 real_type defnew=_sp->norm(b); // comp defect norm
110 if (_verbose>1) // print
111 this->printOutput(std::cout,i,defnew,def);
112 //std::cout << i << " " << defnew << " " << defnew/def << std::endl;
113 def = defnew; // update norm
114 if (all_true(def<def0*_reduction) || max_value(def)<1E-30) // convergence check
115 {
116 res.converged = true;
117 break;
118 }
119 }
120
121 //correct i which is wrong if convergence was not achieved.
122 i=std::min(_maxit,i);
123
124 // print
125 if (_verbose==1)
126 this->printOutput(std::cout,i,def);
127
128 // postprocess preconditioner
129 _prec->post(x);
130
131 // fill statistics
132 res.iterations = i;
133 res.reduction = static_cast<double>(max_value(def/def0));
134 res.conv_rate = pow(res.reduction,1.0/i);
135 res.elapsed = watch.elapsed();
136
137 // final print
138 if (_verbose>0)
139 {
140 std::cout << "=== rate=" << res.conv_rate
141 << ", T=" << res.elapsed
142 << ", TIT=" << res.elapsed/i
143 << ", IT=" << i << std::endl;
144 }
145 }
146
147 protected:
148 using IterativeSolver<X,X>::_op;
149 using IterativeSolver<X,X>::_prec;
150 using IterativeSolver<X,X>::_sp;
151 using IterativeSolver<X,X>::_reduction;
152 using IterativeSolver<X,X>::_maxit;
153 using IterativeSolver<X,X>::_verbose;
154 };
155
156
157 // all these solvers are taken from the SUMO library
159 template<class X>
160 class GradientSolver : public IterativeSolver<X,X> {
161 public:
165 using typename IterativeSolver<X,X>::real_type;
166
167 // copy base class constructors
169
170 // don't shadow four-argument version of apply defined in the base class
171 using IterativeSolver<X,X>::apply;
172
178 virtual void apply (X& x, X& b, InverseOperatorResult& res)
179 {
180 res.clear(); // clear solver statistics
181 Timer watch; // start a timer
182 _prec->pre(x,b); // prepare preconditioner
183 _op->applyscaleadd(-1,x,b); // overwrite b with defect
184
185 X p(x); // create local vectors
186 X q(b);
187
188 real_type def0 = _sp->norm(b); // compute norm
189
190 if (_verbose>0) // printing
191 {
192 std::cout << "=== GradientSolver" << std::endl;
193 if (_verbose>1)
194 {
195 this->printHeader(std::cout);
196 this->printOutput(std::cout,0,def0);
197 }
198 }
199
200 int i=1; real_type def=def0; // loop variables
201 field_type lambda;
202 for ( ; i<=_maxit; i++ )
203 {
204 p = 0; // clear correction
205 _prec->apply(p,b); // apply preconditioner
206 _op->apply(p,q); // q=Ap
207 lambda = _sp->dot(p,b)/_sp->dot(q,p); // minimization
208 x.axpy(lambda,p); // update solution
209 b.axpy(-lambda,q); // update defect
210
211 real_type defnew=_sp->norm(b); // comp defect norm
212 if (_verbose>1) // print
213 this->printOutput(std::cout,i,defnew,def);
214
215 def = defnew; // update norm
216 if (all_true(def<def0*_reduction) || max_value(def)<1E-30) // convergence check
217 {
218 res.converged = true;
219 break;
220 }
221 }
222
223 //correct i which is wrong if convergence was not achieved.
224 i=std::min(_maxit,i);
225
226 if (_verbose==1) // printing for non verbose
227 this->printOutput(std::cout,i,def);
228
229 _prec->post(x); // postprocess preconditioner
230 res.iterations = i; // fill statistics
231 res.reduction = static_cast<double>(max_value(def/def0));
232 res.conv_rate = pow(res.reduction,1.0/i);
233 res.elapsed = watch.elapsed();
234 if (_verbose>0) // final print
235 std::cout << "=== rate=" << res.conv_rate
236 << ", T=" << res.elapsed
237 << ", TIT=" << res.elapsed/i
238 << ", IT=" << i << std::endl;
239 }
240
241 protected:
242 using IterativeSolver<X,X>::_op;
243 using IterativeSolver<X,X>::_prec;
244 using IterativeSolver<X,X>::_sp;
245 using IterativeSolver<X,X>::_reduction;
246 using IterativeSolver<X,X>::_maxit;
247 using IterativeSolver<X,X>::_verbose;
248 };
249
250
251
253 template<class X>
254 class CGSolver : public IterativeSolver<X,X> {
255 public:
259 using typename IterativeSolver<X,X>::real_type;
260
261 // copy base class constructors
263
264 private:
266
267 protected:
268
269 using enableConditionEstimate_t = Dune::Std::bool_constant<(std::is_same<field_type,float>::value || std::is_same<field_type,double>::value)>;
270
271 public:
272
273 // don't shadow four-argument version of apply defined in the base class
274 using IterativeSolver<X,X>::apply;
275
284 scalar_real_type reduction, int maxit, int verbose, bool condition_estimate) : IterativeSolver<X,X>(op, prec, reduction, maxit, verbose),
285 condition_estimate_(condition_estimate)
286 {
287 if (condition_estimate && !(enableConditionEstimate_t{})) {
288 condition_estimate_ = false;
289 std::cerr << "WARNING: Condition estimate was disabled. It is only available for double and float field types!" << std::endl;
290 }
291 }
292
301 scalar_real_type reduction, int maxit, int verbose, bool condition_estimate) : IterativeSolver<X,X>(op, sp, prec, reduction, maxit, verbose),
302 condition_estimate_(condition_estimate)
303 {
304 if (condition_estimate && !(std::is_same<field_type,float>::value || std::is_same<field_type,double>::value)) {
305 condition_estimate_ = false;
306 std::cerr << "WARNING: Condition estimate was disabled. It is only available for double and float field types!" << std::endl;
307 }
308 }
309
310
322 virtual void apply (X& x, X& b, InverseOperatorResult& res)
323 {
324 using std::isfinite;
325
326 res.clear(); // clear solver statistics
327 Timer watch; // start a timer
328 _prec->pre(x,b); // prepare preconditioner
329 _op->applyscaleadd(-1,x,b); // overwrite b with defect
330
331 X p(x); // the search direction
332 X q(x); // a temporary vector
333
334 real_type def0 = _sp->norm(b); // compute norm
335
336 if (!all_true(isfinite(def0))) // check for inf or NaN
337 {
338 if (_verbose>0)
339 std::cout << "=== CGSolver: abort due to infinite or NaN initial defect"
340 << std::endl;
341 DUNE_THROW(SolverAbort, "CGSolver: initial defect=" << def0
342 << " is infinite or NaN");
343 }
344
345 if (max_value(def0)<1E-30) // convergence check
346 {
347 res.converged = true;
348 res.iterations = 0; // fill statistics
349 res.reduction = 0;
350 res.conv_rate = 0;
351 res.elapsed=0;
352 if (_verbose>0) // final print
353 std::cout << "=== rate=" << res.conv_rate
354 << ", T=" << res.elapsed << ", TIT=" << res.elapsed
355 << ", IT=0" << std::endl;
356 return;
357 }
358
359 if (_verbose>0) // printing
360 {
361 std::cout << "=== CGSolver" << std::endl;
362 if (_verbose>1) {
363 this->printHeader(std::cout);
364 this->printOutput(std::cout,0,def0);
365 }
366 }
367
368 // Remember lambda and beta values for condition estimate
369 std::vector<real_type> lambdas(0);
370 std::vector<real_type> betas(0);
371
372 // some local variables
373 real_type def=def0; // loop variables
374 field_type rho,rholast,lambda,alpha,beta;
375
376 // determine initial search direction
377 p = 0; // clear correction
378 _prec->apply(p,b); // apply preconditioner
379 rholast = _sp->dot(p,b); // orthogonalization
380
381 // the loop
382 int i=1;
383 for ( ; i<=_maxit; i++ )
384 {
385 // minimize in given search direction p
386 _op->apply(p,q); // q=Ap
387 alpha = _sp->dot(p,q); // scalar product
388 lambda = rholast/alpha; // minimization
389 Hybrid::ifElse(enableConditionEstimate_t{}, [&](auto id) {
390 if (condition_estimate_)
391 lambdas.push_back(std::real(id(lambda)));
392 });
393 x.axpy(lambda,p); // update solution
394 b.axpy(-lambda,q); // update defect
395
396 // convergence test
397 real_type defnew=_sp->norm(b); // comp defect norm
398
399 if (_verbose>1) // print
400 this->printOutput(std::cout,i,defnew,def);
401
402 def = defnew; // update norm
403 if (!all_true(isfinite(def))) // check for inf or NaN
404 {
405 if (_verbose>0)
406 std::cout << "=== CGSolver: abort due to infinite or NaN defect"
407 << std::endl;
409 "CGSolver: defect=" << def << " is infinite or NaN");
410 }
411
412 if (all_true(def<def0*_reduction) || max_value(def)<1E-30) // convergence check
413 {
414 res.converged = true;
415 break;
416 }
417
418 // determine new search direction
419 q = 0; // clear correction
420 _prec->apply(q,b); // apply preconditioner
421 rho = _sp->dot(q,b); // orthogonalization
422 beta = rho/rholast; // scaling factor
423 Hybrid::ifElse(enableConditionEstimate_t{}, [&](auto id) {
424 if (condition_estimate_)
425 betas.push_back(std::real(id(beta)));
426 });
427 p *= beta; // scale old search direction
428 p += q; // orthogonalization with correction
429 rholast = rho; // remember rho for recurrence
430 }
431
432 //correct i which is wrong if convergence was not achieved.
433 i=std::min(_maxit,i);
434
435 if (_verbose==1) // printing for non verbose
436 this->printOutput(std::cout,i,def);
437
438 _prec->post(x); // postprocess preconditioner
439 res.iterations = i; // fill statistics
440 res.reduction = static_cast<double>(max_value(max_value(def/def0)));
441 res.conv_rate = pow(res.reduction,1.0/i);
442 res.elapsed = watch.elapsed();
443
444 if (_verbose>0) // final print
445 {
446 std::cout << "=== rate=" << res.conv_rate
447 << ", T=" << res.elapsed
448 << ", TIT=" << res.elapsed/i
449 << ", IT=" << i << std::endl;
450 }
451
452
453 if (condition_estimate_) {
454#if HAVE_ARPACKPP
455 Hybrid::ifElse(enableConditionEstimate_t{}, [&](auto id) {
456
457 // Build T matrix which has extreme eigenvalues approximating
458 // those of the original system
459 // (see Y. Saad, Iterative methods for sparse linear systems)
460
462
463 for (auto row = T.createbegin(); row != T.createend(); ++row) {
464 if (row.index() > 0)
465 row.insert(row.index()-1);
466 row.insert(row.index());
467 if (row.index() < T.N() - 1)
468 row.insert(row.index()+1);
469 }
470 for (int row = 0; row < i; ++row) {
471 if (row > 0) {
472 T[row][row-1] = std::sqrt(id(betas[row-1])) / lambdas[row-1];
473 }
474
475 T[row][row] = 1.0 / id(lambdas[row]);
476 if (row > 0) {
477 T[row][row] += betas[row-1] / lambdas[row-1];
478 }
479
480 if (row < i - 1) {
481 T[row][row+1] = std::sqrt(id(betas[row])) / lambdas[row];
482 }
483 }
484
485 // Compute largest and smallest eigenvalue of T matrix and return as estimate
487
488 real_type eps = 0.0;
489 COND_VEC eigv;
490 real_type min_eigv, max_eigv;
491 arpack.computeSymMinMagnitude (eps, eigv, min_eigv);
492 arpack.computeSymMaxMagnitude (eps, eigv, max_eigv);
493
494 res.condition_estimate = id(max_eigv / min_eigv);
495
496 if (this->_verbose > 0) {
497 std::cout << "Min eigv estimate: " << min_eigv << std::endl;
498 std::cout << "Max eigv estimate: " << max_eigv << std::endl;
499 std::cout << "Condition estimate: " << max_eigv / min_eigv << std::endl;
500 }
501 });
502#else
503 std::cerr << "WARNING: Condition estimate was requested. This requires ARPACK, but ARPACK was not found!" << std::endl;
504#endif
505 }
506 }
507
508 private:
509 bool condition_estimate_ = false;
510
511 // Matrix and vector types used for condition estimate
514
515 protected:
516 using IterativeSolver<X,X>::_op;
517 using IterativeSolver<X,X>::_prec;
518 using IterativeSolver<X,X>::_sp;
519 using IterativeSolver<X,X>::_reduction;
520 using IterativeSolver<X,X>::_maxit;
521 using IterativeSolver<X,X>::_verbose;
522 };
523
524
525 // Ronald Kriemanns BiCG-STAB implementation from Sumo
527 template<class X>
528 class BiCGSTABSolver : public IterativeSolver<X,X> {
529 public:
533 using typename IterativeSolver<X,X>::real_type;
534
535 // copy base class constructors
537
538 // don't shadow four-argument version of apply defined in the base class
539 using IterativeSolver<X,X>::apply;
540
548 virtual void apply (X& x, X& b, InverseOperatorResult& res)
549 {
550 using std::abs;
551 const real_type EPSILON=1e-80;
552 using std::abs;
553 double it;
554 field_type rho, rho_new, alpha, beta, h, omega;
555 real_type norm, norm_old, norm_0;
556
557 //
558 // get vectors and matrix
559 //
560 X& r=b;
561 X p(x);
562 X v(x);
563 X t(x);
564 X y(x);
565 X rt(x);
566
567 //
568 // begin iteration
569 //
570
571 // r = r - Ax; rt = r
572 res.clear(); // clear solver statistics
573 Timer watch; // start a timer
574 _prec->pre(x,r); // prepare preconditioner
575 _op->applyscaleadd(-1,x,r); // overwrite b with defect
576
577 rt=r;
578
579 norm = norm_old = norm_0 = _sp->norm(r);
580
581 p=0;
582 v=0;
583
584 rho = 1;
585 alpha = 1;
586 omega = 1;
587
588 if (_verbose>0) // printing
589 {
590 std::cout << "=== BiCGSTABSolver" << std::endl;
591 if (_verbose>1)
592 {
593 this->printHeader(std::cout);
594 this->printOutput(std::cout,0,norm_0);
595 //std::cout << " Iter Defect Rate" << std::endl;
596 //std::cout << " 0" << std::setw(14) << norm_0 << std::endl;
597 }
598 }
599
600 if ( all_true(norm < (_reduction * norm_0)) || max_value(norm)<1E-30)
601 {
602 res.converged = 1;
603 _prec->post(x); // postprocess preconditioner
604 res.iterations = 0; // fill statistics
605 res.reduction = 0;
606 res.conv_rate = 0;
607 res.elapsed = watch.elapsed();
608 return;
609 }
610
611 //
612 // iteration
613 //
614
615 for (it = 0.5; it < _maxit; it+=.5)
616 {
617 //
618 // preprocess, set vecsizes etc.
619 //
620
621 // rho_new = < rt , r >
622 rho_new = _sp->dot(rt,r);
623
624 // look if breakdown occurred
625 if (all_true(abs(rho) <= EPSILON))
626 DUNE_THROW(SolverAbort,"breakdown in BiCGSTAB - rho "
627 << rho << " <= EPSILON " << max_value(EPSILON)
628 << " after " << it << " iterations");
629 if (all_true(abs(omega) <= EPSILON))
630 DUNE_THROW(SolverAbort,"breakdown in BiCGSTAB - omega "
631 << omega << " <= EPSILON " << max_value(EPSILON)
632 << " after " << it << " iterations");
633
634
635 if (it<1)
636 p = r;
637 else
638 {
639 beta = ( rho_new / rho ) * ( alpha / omega );
640 p.axpy(-omega,v); // p = r + beta (p - omega*v)
641 p *= beta;
642 p += r;
643 }
644
645 // y = W^-1 * p
646 y = 0;
647 _prec->apply(y,p); // apply preconditioner
648
649 // v = A * y
650 _op->apply(y,v);
651
652 // alpha = rho_new / < rt, v >
653 h = _sp->dot(rt,v);
654
655 if ( all_true(abs(h) < EPSILON) )
656 DUNE_THROW(SolverAbort,"abs(h) < EPSILON in BiCGSTAB - abs(h) "
657 << abs(h) << " < EPSILON " << max_value(EPSILON)
658 << " after " << it << " iterations");
659
660 alpha = rho_new / h;
661
662 // apply first correction to x
663 // x <- x + alpha y
664 x.axpy(alpha,y);
665
666 // r = r - alpha*v
667 r.axpy(-alpha,v);
668
669 //
670 // test stop criteria
671 //
672
673 norm = _sp->norm(r);
674
675 if (_verbose>1) // print
676 {
677 this->printOutput(std::cout,it,norm,norm_old);
678 }
679
680 if ( all_true(norm < (_reduction * norm_0)) )
681 {
682 res.converged = 1;
683 break;
684 }
685 it+=.5;
686
687 norm_old = norm;
688
689 // y = W^-1 * r
690 y = 0;
691 _prec->apply(y,r);
692
693 // t = A * y
694 _op->apply(y,t);
695
696 // omega = < t, r > / < t, t >
697 omega = _sp->dot(t,r)/_sp->dot(t,t);
698
699 // apply second correction to x
700 // x <- x + omega y
701 x.axpy(omega,y);
702
703 // r = s - omega*t (remember : r = s)
704 r.axpy(-omega,t);
705
706 rho = rho_new;
707
708 //
709 // test stop criteria
710 //
711
712 norm = _sp->norm(r);
713
714 if (_verbose > 1) // print
715 {
716 this->printOutput(std::cout,it,norm,norm_old);
717 }
718
719 if ( all_true(norm < (_reduction * norm_0)) || max_value(norm)<1E-30)
720 {
721 res.converged = 1;
722 break;
723 }
724
725 norm_old = norm;
726 } // end for
727
728 //correct i which is wrong if convergence was not achieved.
729 it=std::min(static_cast<double>(_maxit),it);
730
731 if (_verbose==1) // printing for non verbose
732 this->printOutput(std::cout,it,norm);
733
734 _prec->post(x); // postprocess preconditioner
735 res.iterations = static_cast<int>(std::ceil(it)); // fill statistics
736 res.reduction = static_cast<double>(max_value(norm/norm_0));
737 res.conv_rate = pow(res.reduction,1.0/it);
738 res.elapsed = watch.elapsed();
739 if (_verbose>0) // final print
740 std::cout << "=== rate=" << res.conv_rate
741 << ", T=" << res.elapsed
742 << ", TIT=" << res.elapsed/it
743 << ", IT=" << it << std::endl;
744 }
745
746 protected:
747 using IterativeSolver<X,X>::_op;
748 using IterativeSolver<X,X>::_prec;
749 using IterativeSolver<X,X>::_sp;
750 using IterativeSolver<X,X>::_reduction;
751 using IterativeSolver<X,X>::_maxit;
752 using IterativeSolver<X,X>::_verbose;
753 };
754
761 template<class X>
762 class MINRESSolver : public IterativeSolver<X,X> {
763 public:
767 using typename IterativeSolver<X,X>::real_type;
768
769 // copy base class constructors
771
772 // don't shadow four-argument version of apply defined in the base class
773 using IterativeSolver<X,X>::apply;
774
780 virtual void apply (X& x, X& b, InverseOperatorResult& res)
781 {
782 using std::sqrt;
783 using std::abs;
784 // clear solver statistics
785 res.clear();
786 // start a timer
787 Dune::Timer watch;
788 watch.reset();
789 // prepare preconditioner
790 _prec->pre(x,b);
791 // overwrite rhs with defect
792 _op->applyscaleadd(-1,x,b);
793
794 // compute residual norm
795 real_type def0 = _sp->norm(b);
796
797 // printing
798 if(_verbose > 0) {
799 std::cout << "=== MINRESSolver" << std::endl;
800 if(_verbose > 1) {
801 this->printHeader(std::cout);
802 this->printOutput(std::cout,0,def0);
803 }
804 }
805
806 // check for convergence
807 if( max_value(def0) < 1e-30 ) {
808 res.converged = true;
809 res.iterations = 0;
810 res.reduction = 0;
811 res.conv_rate = 0;
812 res.elapsed = 0.0;
813 // final print
814 if(_verbose > 0)
815 std::cout << "=== rate=" << res.conv_rate
816 << ", T=" << res.elapsed
817 << ", TIT=" << res.elapsed
818 << ", IT=" << res.iterations
819 << std::endl;
820 return;
821 }
822
823 // the defect norm
824 real_type def = def0;
825 // recurrence coefficients as computed in Lanczos algorithm
826 field_type alpha, beta;
827 // diagonal entries of givens rotation
828 std::array<real_type,2> c{{0.0,0.0}};
829 // off-diagonal entries of givens rotation
830 std::array<field_type,2> s{{0.0,0.0}};
831
832 // recurrence coefficients (column k of tridiag matrix T_k)
833 std::array<field_type,3> T{{0.0,0.0,0.0}};
834
835 // the rhs vector of the min problem
836 std::array<field_type,2> xi{{1.0,0.0}};
837
838 // some temporary vectors
839 X z(b), dummy(b);
840
841 // initialize and clear correction
842 z = 0.0;
843 _prec->apply(z,b);
844
845 // beta is real and positive in exact arithmetic
846 // since it is the norm of the basis vectors (in unpreconditioned case)
847 beta = sqrt(_sp->dot(b,z));
848 field_type beta0 = beta;
849
850 // the search directions
851 std::array<X,3> p{{b,b,b}};
852 p[0] = 0.0;
853 p[1] = 0.0;
854 p[2] = 0.0;
855
856 // orthonormal basis vectors (in unpreconditioned case)
857 std::array<X,3> q{{b,b,b}};
858 q[0] = 0.0;
859 q[1] *= real_type(1.0)/beta;
860 q[2] = 0.0;
861
862 z *= real_type(1.0)/beta;
863
864 // the loop
865 int i = 1;
866 for( ; i<=_maxit; i++) {
867
868 dummy = z;
869 int i1 = i%3,
870 i0 = (i1+2)%3,
871 i2 = (i1+1)%3;
872
873 // symmetrically preconditioned Lanczos algorithm (see Greenbaum p.121)
874 _op->apply(z,q[i2]); // q[i2] = Az
875 q[i2].axpy(-beta,q[i0]);
876 // alpha is real since it is the diagonal entry of the hermitian tridiagonal matrix
877 // from the Lanczos Algorithm
878 // so the order in the scalar product doesn't matter even for the complex case
879 alpha = _sp->dot(z,q[i2]);
880 q[i2].axpy(-alpha,q[i1]);
881
882 z = 0.0;
883 _prec->apply(z,q[i2]);
884
885 // beta is real and positive in exact arithmetic
886 // since it is the norm of the basis vectors (in unpreconditioned case)
887 beta = sqrt(_sp->dot(q[i2],z));
888
889 q[i2] *= real_type(1.0)/beta;
890 z *= real_type(1.0)/beta;
891
892 // QR Factorization of recurrence coefficient matrix
893 // apply previous givens rotations to last column of T
894 T[1] = T[2];
895 if(i>2) {
896 T[0] = s[i%2]*T[1];
897 T[1] = c[i%2]*T[1];
898 }
899 if(i>1) {
900 T[2] = c[(i+1)%2]*alpha - s[(i+1)%2]*T[1];
901 T[1] = c[(i+1)%2]*T[1] + s[(i+1)%2]*alpha;
902 }
903 else
904 T[2] = alpha;
905
906 // update QR factorization
907 generateGivensRotation(T[2],beta,c[i%2],s[i%2]);
908 // to last column of T_k
909 T[2] = c[i%2]*T[2] + s[i%2]*beta;
910 // and to the rhs xi of the min problem
911 xi[i%2] = -s[i%2]*xi[(i+1)%2];
912 xi[(i+1)%2] *= c[i%2];
913
914 // compute correction direction
915 p[i2] = dummy;
916 p[i2].axpy(-T[1],p[i1]);
917 p[i2].axpy(-T[0],p[i0]);
918 p[i2] *= real_type(1.0)/T[2];
919
920 // apply correction/update solution
921 x.axpy(beta0*xi[(i+1)%2],p[i2]);
922
923 // remember beta_old
924 T[2] = beta;
925
926 // check for convergence
927 // the last entry in the rhs of the min-problem is the residual
928 real_type defnew = abs(beta0*xi[i%2]);
929
930 if(_verbose > 1)
931 this->printOutput(std::cout,i,defnew,def);
932
933 def = defnew;
934 if(all_true(def < def0*_reduction)
935 || max_value(def) < 1e-30 || i == _maxit ) {
936 res.converged = true;
937 break;
938 }
939 } // end for
940
941 if(_verbose == 1)
942 this->printOutput(std::cout,i,def);
943
944 // postprocess preconditioner
945 _prec->post(x);
946 // fill statistics
947 res.iterations = i;
948 res.reduction = static_cast<double>(max_value(def/def0));
949 res.conv_rate = pow(res.reduction,1.0/i);
950 res.elapsed = watch.elapsed();
951
952 // final print
953 if(_verbose > 0) {
954 std::cout << "=== rate=" << res.conv_rate
955 << ", T=" << res.elapsed
956 << ", TIT=" << res.elapsed/i
957 << ", IT=" << i << std::endl;
958 }
959 }
960
961 private:
962
963 // helper function to extract the real value of a real or complex number
964 inline
965 real_type to_real(const real_type & v)
966 {
967 return v;
968 }
969
970 inline
971 real_type to_real(const std::complex<real_type> & v)
972 {
973 return v.real();
974 }
975
976 void generateGivensRotation(field_type &dx, field_type &dy, real_type &cs, field_type &sn)
977 {
978 using std::sqrt;
979 using std::abs;
980 using std::max;
981 using std::min;
982 const real_type eps = 1e-15;
983 real_type norm_dx = abs(dx);
984 real_type norm_dy = abs(dy);
985 real_type norm_max = max(norm_dx, norm_dy);
986 real_type norm_min = min(norm_dx, norm_dy);
987 real_type temp = norm_min/norm_max;
988 // we rewrite the code in a vectorizable fashion
989 cs = cond(norm_dy < eps,
990 real_type(1.0),
991 cond(norm_dx < eps,
992 real_type(0.0),
993 cond(norm_dy > norm_dx,
994 real_type(1.0)/sqrt(real_type(1.0) + temp*temp)*temp,
995 real_type(1.0)/sqrt(real_type(1.0) + temp*temp)
996 )));
997 sn = cond(norm_dy < eps,
998 field_type(0.0),
999 cond(norm_dx < eps,
1000 field_type(1.0),
1001 cond(norm_dy > norm_dx,
1002 // dy and dx are real in exact arithmetic
1003 // thus dx*dy is real so we can explicitly enforce it
1004 field_type(1.0)/sqrt(real_type(1.0) + temp*temp)*dx*dy/norm_dx/norm_dy,
1005 // dy and dx is real in exact arithmetic
1006 // so we don't have to conjugate both of them
1007 field_type(1.0)/sqrt(real_type(1.0) + temp*temp)*dy/dx
1008 )));
1009 }
1010
1011 protected:
1012 using IterativeSolver<X,X>::_op;
1013 using IterativeSolver<X,X>::_prec;
1014 using IterativeSolver<X,X>::_sp;
1015 using IterativeSolver<X,X>::_reduction;
1016 using IterativeSolver<X,X>::_maxit;
1017 using IterativeSolver<X,X>::_verbose;
1018 };
1019
1033 template<class X, class Y=X, class F = Y>
1035 {
1036 public:
1037 using typename IterativeSolver<X,Y>::domain_type;
1038 using typename IterativeSolver<X,Y>::range_type;
1039 using typename IterativeSolver<X,Y>::field_type;
1040 using typename IterativeSolver<X,Y>::real_type;
1041
1042 private:
1044
1046 using fAlloc = ReboundAllocatorType<X,field_type>;
1048 using rAlloc = ReboundAllocatorType<X,real_type>;
1049
1050 public:
1051
1058 RestartedGMResSolver (LinearOperator<X,Y>& op, Preconditioner<X,Y>& prec, scalar_real_type reduction, int restart, int maxit, int verbose) :
1059 IterativeSolver<X,Y>::IterativeSolver(op,prec,reduction,maxit,verbose),
1060 _restart(restart)
1061 {}
1062
1069 RestartedGMResSolver (LinearOperator<X,Y>& op, ScalarProduct<X>& sp, Preconditioner<X,Y>& prec, scalar_real_type reduction, int restart, int maxit, int verbose) :
1070 IterativeSolver<X,Y>::IterativeSolver(op,sp,prec,reduction,maxit,verbose),
1071 _restart(restart)
1072 {}
1073
1082 virtual void apply (X& x, Y& b, InverseOperatorResult& res)
1083 {
1084 apply(x,b,max_value(_reduction),res);
1085 }
1086
1095 virtual void apply (X& x, Y& b, double reduction, InverseOperatorResult& res)
1096 {
1097 using std::abs;
1098 const real_type EPSILON = 1e-80;
1099 const int m = _restart;
1100 real_type norm, norm_old = 0.0, norm_0;
1101 int j = 1;
1102 std::vector<field_type,fAlloc> s(m+1), sn(m);
1103 std::vector<real_type,rAlloc> cs(m);
1104 // need copy of rhs if GMRes has to be restarted
1105 Y b2(b);
1106 // helper vector
1107 Y w(b);
1108 std::vector< std::vector<field_type,fAlloc> > H(m+1,s);
1109 std::vector<F> v(m+1,b);
1110
1111 // start timer
1112 Dune::Timer watch;
1113 watch.reset();
1114
1115 // clear solver statistics and set res.converged to false
1116 res.clear();
1117 _prec->pre(x,b);
1118
1119 // calculate defect and overwrite rhs with it
1120 _op->applyscaleadd(-1.0,x,b); // b -= Ax
1121 // calculate preconditioned defect
1122 v[0] = 0.0; _prec->apply(v[0],b); // r = W^-1 b
1123 norm_0 = _sp->norm(v[0]);
1124 norm = norm_0;
1125 norm_old = norm;
1126
1127 // print header
1128 if(_verbose > 0)
1129 {
1130 std::cout << "=== RestartedGMResSolver" << std::endl;
1131 if(_verbose > 1) {
1132 this->printHeader(std::cout);
1133 this->printOutput(std::cout,0,norm_0);
1134 }
1135 }
1136
1137 if(all_true(norm_0 < EPSILON)) {
1138 _prec->post(x);
1139 res.converged = true;
1140 if(_verbose > 0) // final print
1141 print_result(res);
1142 }
1143
1144 while(j <= _maxit && res.converged != true) {
1145
1146 int i = 0;
1147 v[0] *= real_type(1.0)/norm;
1148 s[0] = norm;
1149 for(i=1; i<m+1; i++)
1150 s[i] = 0.0;
1151
1152 for(i=0; i < m && j <= _maxit && res.converged != true; i++, j++) {
1153 w = 0.0;
1154 // use v[i+1] as temporary vector
1155 v[i+1] = 0.0;
1156 // do Arnoldi algorithm
1157 _op->apply(v[i],v[i+1]);
1158 _prec->apply(w,v[i+1]);
1159 for(int k=0; k<i+1; k++) {
1160 // notice that _sp->dot(v[k],w) = v[k]\adjoint w
1161 // so one has to pay attention to the order
1162 // in the scalar product for the complex case
1163 // doing the modified Gram-Schmidt algorithm
1164 H[k][i] = _sp->dot(v[k],w);
1165 // w -= H[k][i] * v[k]
1166 w.axpy(-H[k][i],v[k]);
1167 }
1168 H[i+1][i] = _sp->norm(w);
1169 if(all_true(abs(H[i+1][i]) < EPSILON))
1171 "breakdown in GMRes - |w| == 0.0 after " << j << " iterations");
1172
1173 // normalize new vector
1174 v[i+1] = w; v[i+1] *= real_type(1.0)/H[i+1][i];
1175
1176 // update QR factorization
1177 for(int k=0; k<i; k++)
1178 applyPlaneRotation(H[k][i],H[k+1][i],cs[k],sn[k]);
1179
1180 // compute new givens rotation
1181 generatePlaneRotation(H[i][i],H[i+1][i],cs[i],sn[i]);
1182 // finish updating QR factorization
1183 applyPlaneRotation(H[i][i],H[i+1][i],cs[i],sn[i]);
1184 applyPlaneRotation(s[i],s[i+1],cs[i],sn[i]);
1185
1186 // norm of the defect is the last component the vector s
1187 norm = abs(s[i+1]);
1188
1189 // print current iteration statistics
1190 if(_verbose > 1) {
1191 this->printOutput(std::cout,j,norm,norm_old);
1192 }
1193
1194 norm_old = norm;
1195
1196 // check convergence
1197 if(all_true(norm < real_type(reduction) * norm_0))
1198 res.converged = true;
1199
1200 } // end for
1201
1202 // calculate update vector
1203 w = 0.0;
1204 update(w,i,H,s,v);
1205 // and current iterate
1206 x += w;
1207
1208 // restart GMRes if convergence was not achieved,
1209 // i.e. linear defect has not reached desired reduction
1210 // and if j < _maxit
1211 if( res.converged != true && j <= _maxit ) {
1212
1213 if(_verbose > 0)
1214 std::cout << "=== GMRes::restart" << std::endl;
1215 // get saved rhs
1216 b = b2;
1217 // calculate new defect
1218 _op->applyscaleadd(-1.0,x,b); // b -= Ax;
1219 // calculate preconditioned defect
1220 v[0] = 0.0;
1221 _prec->apply(v[0],b);
1222 norm = _sp->norm(v[0]);
1223 norm_old = norm;
1224 }
1225
1226 } //end while
1227
1228 // postprocess preconditioner
1229 _prec->post(x);
1230
1231 // save solver statistics
1232 res.iterations = j-1; // it has to be j-1!!!
1233 res.reduction = static_cast<double>(max_value(norm/norm_0));
1234 res.conv_rate = pow(res.reduction,1.0/(j-1));
1235 res.elapsed = watch.elapsed();
1236
1237 if(_verbose>0)
1238 print_result(res);
1239
1240 }
1241
1242 private :
1243
1244 void print_result(const InverseOperatorResult& res) const {
1245 int k = res.iterations>0 ? res.iterations : 1;
1246 std::cout << "=== rate=" << res.conv_rate
1247 << ", T=" << res.elapsed
1248 << ", TIT=" << res.elapsed/k
1249 << ", IT=" << res.iterations
1250 << std::endl;
1251 }
1252
1253 void update(X& w, int i,
1254 const std::vector<std::vector<field_type,fAlloc> >& H,
1255 const std::vector<field_type,fAlloc>& s,
1256 const std::vector<X>& v) {
1257 // solution vector of the upper triangular system
1258 std::vector<field_type,fAlloc> y(s);
1259
1260 // backsolve
1261 for(int a=i-1; a>=0; a--) {
1262 field_type rhs(s[a]);
1263 for(int b=a+1; b<i; b++)
1264 rhs -= H[a][b]*y[b];
1265 y[a] = rhs/H[a][a];
1266
1267 // compute update on the fly
1268 // w += y[a]*v[a]
1269 w.axpy(y[a],v[a]);
1270 }
1271 }
1272
1273 template<typename T>
1274 typename std::enable_if<std::is_same<field_type,real_type>::value,T>::type conjugate(const T& t) {
1275 return t;
1276 }
1277
1278 template<typename T>
1279 typename std::enable_if<!std::is_same<field_type,real_type>::value,T>::type conjugate(const T& t) {
1280 using std::conj;
1281 return conj(t);
1282 }
1283
1284 // helper function to extract the real value of a real or complex number
1285 inline
1286 real_type to_real(const real_type & v)
1287 {
1288 return v;
1289 }
1290
1291 inline
1292 real_type to_real(const std::complex<real_type> & v)
1293 {
1294 return v.real();
1295 }
1296
1297 void
1298 generatePlaneRotation(field_type &dx, field_type &dy, real_type &cs, field_type &sn)
1299 {
1300 using std::sqrt;
1301 using std::abs;
1302 using std::max;
1303 using std::min;
1304 const real_type eps = 1e-15;
1305 real_type norm_dx = abs(dx);
1306 real_type norm_dy = abs(dy);
1307 real_type norm_max = max(norm_dx, norm_dy);
1308 real_type norm_min = min(norm_dx, norm_dy);
1309 real_type temp = norm_min/norm_max;
1310 // we rewrite the code in a vectorizable fashion
1311 cs = cond(norm_dy < eps,
1312 real_type(1.0),
1313 cond(norm_dx < eps,
1314 real_type(0.0),
1315 cond(norm_dy > norm_dx,
1316 real_type(1.0)/sqrt(real_type(1.0) + temp*temp)*temp,
1317 real_type(1.0)/sqrt(real_type(1.0) + temp*temp)
1318 )));
1319 sn = cond(norm_dy < eps,
1320 field_type(0.0),
1321 cond(norm_dx < eps,
1322 field_type(1.0),
1323 cond(norm_dy > norm_dx,
1324 field_type(1.0)/sqrt(real_type(1.0) + temp*temp)*dx*conjugate(dy)/norm_dx/norm_dy,
1325 field_type(1.0)/sqrt(real_type(1.0) + temp*temp)*conjugate(dy/dx)
1326 )));
1327 }
1328
1329
1330 void
1331 applyPlaneRotation(field_type &dx, field_type &dy, real_type &cs, field_type &sn)
1332 {
1333 field_type temp = cs * dx + sn * dy;
1334 dy = -conjugate(sn) * dx + cs * dy;
1335 dx = temp;
1336 }
1337
1338 using IterativeSolver<X,Y>::_op;
1339 using IterativeSolver<X,Y>::_prec;
1340 using IterativeSolver<X,Y>::_sp;
1341 using IterativeSolver<X,Y>::_reduction;
1342 using IterativeSolver<X,Y>::_maxit;
1343 using IterativeSolver<X,Y>::_verbose;
1344 int _restart;
1345 };
1346
1347
1361 template<class X>
1363 {
1364 public:
1365 using typename IterativeSolver<X,X>::domain_type;
1366 using typename IterativeSolver<X,X>::range_type;
1367 using typename IterativeSolver<X,X>::field_type;
1368 using typename IterativeSolver<X,X>::real_type;
1369
1370 private:
1372
1374 using fAlloc = ReboundAllocatorType<X,field_type>;
1375
1376 public:
1377
1378 // don't shadow four-argument version of apply defined in the base class
1379 using IterativeSolver<X,X>::apply;
1380
1387 GeneralizedPCGSolver (LinearOperator<X,X>& op, Preconditioner<X,X>& prec, scalar_real_type reduction, int maxit, int verbose, int restart = 10) :
1388 IterativeSolver<X,X>::IterativeSolver(op,prec,reduction,maxit,verbose),
1389 _restart(restart)
1390 {}
1391
1399 GeneralizedPCGSolver (LinearOperator<X,X>& op, ScalarProduct<X>& sp, Preconditioner<X,X>& prec, scalar_real_type reduction, int maxit, int verbose, int restart = 10) :
1400 IterativeSolver<X,X>::IterativeSolver(op,sp,prec,reduction,maxit,verbose),
1401 _restart(restart)
1402 {}
1403
1409 virtual void apply (X& x, X& b, InverseOperatorResult& res)
1410 {
1411 res.clear(); // clear solver statistics
1412 Timer watch; // start a timer
1413 _prec->pre(x,b); // prepare preconditioner
1414 _op->applyscaleadd(-1,x,b); // overwrite b with defect
1415
1416 std::vector<std::shared_ptr<X> > p(_restart);
1417 std::vector<field_type,fAlloc> pp(_restart);
1418 X q(x); // a temporary vector
1419 X prec_res(x); // a temporary vector for preconditioner output
1420
1421 p[0].reset(new X(x));
1422
1423 real_type def0 = _sp->norm(b); // compute norm
1424 if ( max_value(def0) < 1E-30 ) // convergence check
1425 {
1426 res.converged = true;
1427 res.iterations = 0; // fill statistics
1428 res.reduction = 0;
1429 res.conv_rate = 0;
1430 res.elapsed=0;
1431 if (_verbose>0) // final print
1432 std::cout << "=== rate=" << res.conv_rate
1433 << ", T=" << res.elapsed << ", TIT=" << res.elapsed
1434 << ", IT=0" << std::endl;
1435 return;
1436 }
1437
1438 if (_verbose>0) // printing
1439 {
1440 std::cout << "=== GeneralizedPCGSolver" << std::endl;
1441 if (_verbose>1) {
1442 this->printHeader(std::cout);
1443 this->printOutput(std::cout,0,def0);
1444 }
1445 }
1446 // some local variables
1447 real_type def=def0; // loop variables
1448 field_type rho, lambda;
1449
1450 int i=0;
1451 int ii=0;
1452 // determine initial search direction
1453 *(p[0]) = 0; // clear correction
1454 _prec->apply(*(p[0]),b); // apply preconditioner
1455 rho = _sp->dot(*(p[0]),b); // orthogonalization
1456 _op->apply(*(p[0]),q); // q=Ap
1457 pp[0] = _sp->dot(*(p[0]),q); // scalar product
1458 lambda = rho/pp[0]; // minimization
1459 x.axpy(lambda,*(p[0])); // update solution
1460 b.axpy(-lambda,q); // update defect
1461
1462 // convergence test
1463 real_type defnew=_sp->norm(b); // comp defect norm
1464 ++i;
1465 if (_verbose>1) // print
1466 this->printOutput(std::cout,i,defnew,def);
1467 def = defnew; // update norm
1468 if (all_true(def<def0*_reduction) || max_value(def)<1E-30) // convergence check
1469 {
1470 res.converged = true;
1471 if (_verbose>0) // final print
1472 {
1473 std::cout << "=== rate=" << res.conv_rate
1474 << ", T=" << res.elapsed
1475 << ", TIT=" << res.elapsed
1476 << ", IT=" << 1 << std::endl;
1477 }
1478 return;
1479 }
1480
1481 while(i<_maxit) {
1482 // the loop
1483 int end=std::min(_restart, _maxit-i+1);
1484 for (ii=1; ii<end; ++ii )
1485 {
1486 //std::cout<<" ii="<<ii<<" i="<<i<<std::endl;
1487 // compute next conjugate direction
1488 prec_res = 0; // clear correction
1489 _prec->apply(prec_res,b); // apply preconditioner
1490
1491 p[ii].reset(new X(prec_res));
1492 _op->apply(prec_res, q);
1493
1494 for(int j=0; j<ii; ++j) {
1495 rho =_sp->dot(q,*(p[j]))/pp[j];
1496 p[ii]->axpy(-rho, *(p[j]));
1497 }
1498
1499 // minimize in given search direction
1500 _op->apply(*(p[ii]),q); // q=Ap
1501 pp[ii] = _sp->dot(*(p[ii]),q); // scalar product
1502 rho = _sp->dot(*(p[ii]),b); // orthogonalization
1503 lambda = rho/pp[ii]; // minimization
1504 x.axpy(lambda,*(p[ii])); // update solution
1505 b.axpy(-lambda,q); // update defect
1506
1507 // convergence test
1508 defnew=_sp->norm(b); // comp defect norm
1509
1510 ++i;
1511 if (_verbose>1) // print
1512 this->printOutput(std::cout,i,defnew,def);
1513
1514 def = defnew; // update norm
1515 if (all_true(def<def0*_reduction) || max_value(def)<1E-30) // convergence check
1516 {
1517 res.converged = true;
1518 break;
1519 }
1520 }
1521 if(res.converged)
1522 break;
1523 if(end==_restart) {
1524 *(p[0])=*(p[_restart-1]);
1525 pp[0]=pp[_restart-1];
1526 }
1527 }
1528
1529 // postprocess preconditioner
1530 _prec->post(x);
1531
1532 // fill statistics
1533 res.iterations = i;
1534 res.reduction = static_cast<double>(max_value(def/def0));
1535 res.conv_rate = pow(res.reduction,1.0/i);
1536 res.elapsed = watch.elapsed();
1537
1538 if (_verbose>0) // final print
1539 {
1540 std::cout << "=== rate=" << res.conv_rate
1541 << ", T=" << res.elapsed
1542 << ", TIT=" << res.elapsed/i
1543 << ", IT=" << i+1 << std::endl;
1544 }
1545 }
1546
1547 private:
1548 using IterativeSolver<X,X>::_op;
1549 using IterativeSolver<X,X>::_prec;
1550 using IterativeSolver<X,X>::_sp;
1551 using IterativeSolver<X,X>::_reduction;
1552 using IterativeSolver<X,X>::_maxit;
1553 using IterativeSolver<X,X>::_verbose;
1554 int _restart;
1555 };
1556
1559} // end namespace
1560
1561#endif
Implementation of the BCRSMatrix class.
Wrapper to use a range of ARPACK++ eigenvalue solvers.
Definition: arpackpp.hh:242
void computeSymMaxMagnitude(const Real &epsilon, BlockVector &x, Real &lambda) const
Assume the matrix to be square, symmetric and perform IRLM to compute an approximation lambda of its ...
Definition: arpackpp.hh:286
void computeSymMinMagnitude(const Real &epsilon, BlockVector &x, Real &lambda) const
Assume the matrix to be square, symmetric and perform IRLM to compute an approximation lambda of its ...
Definition: arpackpp.hh:388
A sparse block matrix with compressed row storage.
Definition: bcrsmatrix.hh:423
@ row_wise
Build in a row-wise manner.
Definition: bcrsmatrix.hh:480
CreateIterator createend()
get create iterator pointing to one after the last block
Definition: bcrsmatrix.hh:1062
CreateIterator createbegin()
get initial create iterator
Definition: bcrsmatrix.hh:1056
size_type N() const
number of rows (counted in blocks)
Definition: bcrsmatrix.hh:1894
Bi-conjugate Gradient Stabilized (BiCG-STAB)
Definition: solvers.hh:528
virtual void apply(X &x, X &b, InverseOperatorResult &res)
Apply inverse operator.
Definition: solvers.hh:548
A vector of blocks with memory management.
Definition: bvector.hh:317
conjugate gradient method
Definition: solvers.hh:254
CGSolver(LinearOperator< X, X > &op, Preconditioner< X, X > &prec, scalar_real_type reduction, int maxit, int verbose, bool condition_estimate)
Constructor to initialize a CG solver.
Definition: solvers.hh:283
CGSolver(LinearOperator< X, X > &op, ScalarProduct< X > &sp, Preconditioner< X, X > &prec, scalar_real_type reduction, int maxit, int verbose, bool condition_estimate)
Constructor to initialize a CG solver.
Definition: solvers.hh:300
virtual void apply(X &x, X &b, InverseOperatorResult &res)
Apply inverse operator.
Definition: solvers.hh:322
Generalized preconditioned conjugate gradient solver.
Definition: solvers.hh:1363
GeneralizedPCGSolver(LinearOperator< X, X > &op, Preconditioner< X, X > &prec, scalar_real_type reduction, int maxit, int verbose, int restart=10)
Set up nonlinear preconditioned conjugate gradient solver.
Definition: solvers.hh:1387
virtual void apply(X &x, X &b, InverseOperatorResult &res)
Apply inverse operator.
Definition: solvers.hh:1409
GeneralizedPCGSolver(LinearOperator< X, X > &op, ScalarProduct< X > &sp, Preconditioner< X, X > &prec, scalar_real_type reduction, int maxit, int verbose, int restart=10)
Set up nonlinear preconditioned conjugate gradient solver.
Definition: solvers.hh:1399
gradient method
Definition: solvers.hh:160
virtual void apply(X &x, X &b, InverseOperatorResult &res)
Apply inverse operator.
Definition: solvers.hh:178
void printHeader(std::ostream &s) const
helper function for printing header of solver output
Definition: solver.hh:155
void printOutput(std::ostream &s, const CountType &iter, const DataType &norm, const DataType &norm_old) const
helper function for printing solver output
Definition: solver.hh:164
Y range_type
Type of the range of the operator to be inverted.
Definition: solver.hh:97
SimdScalar< real_type > scalar_real_type
scalar type underlying the field_type
Definition: solver.hh:106
X domain_type
Type of the domain of the operator to be inverted.
Definition: solver.hh:94
X::field_type field_type
The field type of the operator.
Definition: solver.hh:100
FieldTraits< field_type >::real_type real_type
The real type of the field type (is the same if using real numbers, but differs for std::complex)
Definition: solver.hh:103
Base class for all implementations of iterative solvers.
Definition: solver.hh:195
IterativeSolver(LinearOperator< X, X > &op, Preconditioner< X, X > &prec, scalar_real_type reduction, int maxit, int verbose)
General constructor to initialize an iterative solver.
Definition: solver.hh:222
Preconditioned loop solver.
Definition: solvers.hh:56
virtual void apply(X &x, X &b, InverseOperatorResult &res)
Apply inverse operator,.
Definition: solvers.hh:70
Minimal Residual Method (MINRES)
Definition: solvers.hh:762
virtual void apply(X &x, X &b, InverseOperatorResult &res)
Apply inverse operator.
Definition: solvers.hh:780
implements the Generalized Minimal Residual (GMRes) method
Definition: solvers.hh:1035
RestartedGMResSolver(LinearOperator< X, Y > &op, Preconditioner< X, Y > &prec, scalar_real_type reduction, int restart, int maxit, int verbose)
Set up RestartedGMResSolver solver.
Definition: solvers.hh:1058
RestartedGMResSolver(LinearOperator< X, Y > &op, ScalarProduct< X > &sp, Preconditioner< X, Y > &prec, scalar_real_type reduction, int restart, int maxit, int verbose)
Set up RestartedGMResSolver solver.
Definition: solvers.hh:1069
virtual void apply(X &x, Y &b, double reduction, InverseOperatorResult &res)
Apply inverse operator.
Definition: solvers.hh:1095
virtual void apply(X &x, Y &b, InverseOperatorResult &res)
Apply inverse operator.
Definition: solvers.hh:1082
Base class for scalar product and norm computation.
Definition: scalarproducts.hh:46
Thrown when a solver aborts due to some problem.
Definition: istlexception.hh:43
A simple stop watch.
Definition: timer.hh:41
void reset() noexcept
Reset timer while keeping the running/stopped state.
Definition: timer.hh:55
double elapsed() const noexcept
Get elapsed user-time from last reset until now/last stop in seconds.
Definition: timer.hh:75
A few common exception classes.
std::integral_constant< bool, value > bool_constant
A template alias for std::integral_constant<bool, value>
Definition: type_traits.hh:118
#define DUNE_THROW(E, m)
Definition: exceptions.hh:216
decltype(auto) ifElse(const Condition &condition, IfFunc &&ifFunc, ElseFunc &&elseFunc)
A conditional expression.
Definition: hybridutilities.hh:389
Dune namespace.
Definition: alignedallocator.hh:10
const T1 cond(bool b, const T1 &v1, const T2 &v2)
conditional evaluate
Definition: conditional.hh:26
Define general, extensible interface for operators. The available implementation wraps a matrix.
Utilities for reduction like operations on ranges.
Define base class for scalar product and norm.
Define general, extensible interface for inverse operators.
Statistics about the application of an inverse operator.
Definition: solver.hh:41
double condition_estimate
Estimate of condition number.
Definition: solver.hh:71
double elapsed
Elapsed time in seconds.
Definition: solver.hh:74
int iterations
Number of iterations.
Definition: solver.hh:59
double reduction
Reduction achieved: .
Definition: solver.hh:62
void clear()
Resets all data.
Definition: solver.hh:49
double conv_rate
Convergence rate (average reduction per step)
Definition: solver.hh:68
bool converged
True if convergence criterion has been met.
Definition: solver.hh:65
A simple timing class.
Creative Commons License   |  Legal Statements / Impressum  |  Hosted by TU Dresden  |  generated with Hugo v0.111.3 (Nov 24, 23:30, 2024)