Dune Core Modules (2.4.1)

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 <cmath>
8#include <complex>
9#include <iostream>
10#include <iomanip>
11#include <memory>
12#include <string>
13#include <vector>
14
15#include "istlexception.hh"
16#include "operators.hh"
17#include "scalarproducts.hh"
18#include "solver.hh"
19#include "preconditioner.hh"
20#include <dune/common/array.hh>
22#include <dune/common/timer.hh>
25
26namespace Dune {
44 //=====================================================================
45 // Implementation of this interface
46 //=====================================================================
47
56 template<class X>
57 class LoopSolver : public InverseOperator<X,X> {
58 public:
60 typedef X domain_type;
62 typedef X range_type;
64 typedef typename X::field_type field_type;
66 typedef typename FieldTraits<field_type>::real_type real_type;
67
87 template<class L, class P>
88 LoopSolver (L& op, P& prec,
89 real_type reduction, int maxit, int verbose) :
90 ssp(), _op(op), _prec(prec), _sp(ssp), _reduction(reduction), _maxit(maxit), _verbose(verbose)
91 {
92 static_assert(static_cast<int>(L::category) == static_cast<int>(P::category),
93 "L and P have to have the same category!");
94 static_assert(static_cast<int>(L::category) == static_cast<int>(SolverCategory::sequential),
95 "L has to be sequential!");
96 }
97
118 template<class L, class S, class P>
119 LoopSolver (L& op, S& sp, P& prec,
120 real_type reduction, int maxit, int verbose) :
121 _op(op), _prec(prec), _sp(sp), _reduction(reduction), _maxit(maxit), _verbose(verbose)
122 {
123 static_assert(static_cast<int>(L::category) == static_cast<int>(P::category),
124 "L and P must have the same category!");
125 static_assert(static_cast<int>(L::category) == static_cast<int>(S::category),
126 "L and S must have the same category!");
127 }
128
129
131 virtual void apply (X& x, X& b, InverseOperatorResult& res)
132 {
133 // clear solver statistics
134 res.clear();
135
136 // start a timer
137 Timer watch;
138
139 // prepare preconditioner
140 _prec.pre(x,b);
141
142 // overwrite b with defect
143 _op.applyscaleadd(-1,x,b);
144
145 // compute norm, \todo parallelization
146 real_type def0 = _sp.norm(b);
147
148 // printing
149 if (_verbose>0)
150 {
151 std::cout << "=== LoopSolver" << std::endl;
152 if (_verbose>1)
153 {
154 this->printHeader(std::cout);
155 this->printOutput(std::cout,0,def0);
156 }
157 }
158
159 // allocate correction vector
160 X v(x);
161
162 // iteration loop
163 int i=1; real_type def=def0;
164 for ( ; i<=_maxit; i++ )
165 {
166 v = 0; // clear correction
167 _prec.apply(v,b); // apply preconditioner
168 x += v; // update solution
169 _op.applyscaleadd(-1,v,b); // update defect
170 real_type defnew=_sp.norm(b); // comp defect norm
171 if (_verbose>1) // print
172 this->printOutput(std::cout,i,defnew,def);
173 //std::cout << i << " " << defnew << " " << defnew/def << std::endl;
174 def = defnew; // update norm
175 if (def<def0*_reduction || def<1E-30) // convergence check
176 {
177 res.converged = true;
178 break;
179 }
180 }
181
182 //correct i which is wrong if convergence was not achieved.
183 i=std::min(_maxit,i);
184
185 // print
186 if (_verbose==1)
187 this->printOutput(std::cout,i,def);
188
189 // postprocess preconditioner
190 _prec.post(x);
191
192 // fill statistics
193 res.iterations = i;
194 res.reduction = def/def0;
195 res.conv_rate = pow(res.reduction,1.0/i);
196 res.elapsed = watch.elapsed();
197
198 // final print
199 if (_verbose>0)
200 {
201 std::cout << "=== rate=" << res.conv_rate
202 << ", T=" << res.elapsed
203 << ", TIT=" << res.elapsed/i
204 << ", IT=" << i << std::endl;
205 }
206 }
207
209 virtual void apply (X& x, X& b, double reduction, InverseOperatorResult& res)
210 {
211 real_type saved_reduction = _reduction;
212 _reduction = reduction;
213 (*this).apply(x,b,res);
214 _reduction = saved_reduction;
215 }
216
217 private:
220 Preconditioner<X,X>& _prec;
221 ScalarProduct<X>& _sp;
222 real_type _reduction;
223 int _maxit;
224 int _verbose;
225 };
226
227
228 // all these solvers are taken from the SUMO library
230 template<class X>
231 class GradientSolver : public InverseOperator<X,X> {
232 public:
234 typedef X domain_type;
236 typedef X range_type;
238 typedef typename X::field_type field_type;
240 typedef typename FieldTraits<field_type>::real_type real_type;
241
242
248 template<class L, class P>
249 GradientSolver (L& op, P& prec,
250 real_type reduction, int maxit, int verbose) :
251 ssp(), _op(op), _prec(prec), _sp(ssp), _reduction(reduction), _maxit(maxit), _verbose(verbose)
252 {
253 static_assert(static_cast<int>(L::category) == static_cast<int>(P::category),
254 "L and P have to have the same category!");
255 static_assert(static_cast<int>(L::category) == static_cast<int>(SolverCategory::sequential),
256 "L has to be sequential!");
257 }
263 template<class L, class S, class P>
264 GradientSolver (L& op, S& sp, P& prec,
265 real_type reduction, int maxit, int verbose) :
266 _op(op), _prec(prec), _sp(sp), _reduction(reduction), _maxit(maxit), _verbose(verbose)
267 {
268 static_assert(static_cast<int>(L::category) == static_cast<int>(P::category),
269 "L and P have to have the same category!");
270 static_assert(static_cast<int>(L::category) == static_cast<int>(S::category),
271 "L and S have to have the same category!");
272 }
273
279 virtual void apply (X& x, X& b, InverseOperatorResult& res)
280 {
281 res.clear(); // clear solver statistics
282 Timer watch; // start a timer
283 _prec.pre(x,b); // prepare preconditioner
284 _op.applyscaleadd(-1,x,b); // overwrite b with defect
285
286 X p(x); // create local vectors
287 X q(b);
288
289 real_type def0 = _sp.norm(b); // compute norm
290
291 if (_verbose>0) // printing
292 {
293 std::cout << "=== GradientSolver" << std::endl;
294 if (_verbose>1)
295 {
296 this->printHeader(std::cout);
297 this->printOutput(std::cout,0,def0);
298 }
299 }
300
301 int i=1; real_type def=def0; // loop variables
302 field_type lambda;
303 for ( ; i<=_maxit; i++ )
304 {
305 p = 0; // clear correction
306 _prec.apply(p,b); // apply preconditioner
307 _op.apply(p,q); // q=Ap
308 lambda = _sp.dot(p,b)/_sp.dot(q,p); // minimization
309 x.axpy(lambda,p); // update solution
310 b.axpy(-lambda,q); // update defect
311
312 real_type defnew=_sp.norm(b); // comp defect norm
313 if (_verbose>1) // print
314 this->printOutput(std::cout,i,defnew,def);
315
316 def = defnew; // update norm
317 if (def<def0*_reduction || def<1E-30) // convergence check
318 {
319 res.converged = true;
320 break;
321 }
322 }
323
324 //correct i which is wrong if convergence was not achieved.
325 i=std::min(_maxit,i);
326
327 if (_verbose==1) // printing for non verbose
328 this->printOutput(std::cout,i,def);
329
330 _prec.post(x); // postprocess preconditioner
331 res.iterations = i; // fill statistics
332 res.reduction = static_cast<double>(def/def0);
333 res.conv_rate = static_cast<double>(pow(res.reduction,1.0/i));
334 res.elapsed = watch.elapsed();
335 if (_verbose>0) // final print
336 std::cout << "=== rate=" << res.conv_rate
337 << ", T=" << res.elapsed
338 << ", TIT=" << res.elapsed/i
339 << ", IT=" << i << std::endl;
340 }
341
347 virtual void apply (X& x, X& b, double reduction, InverseOperatorResult& res)
348 {
349 real_type saved_reduction = _reduction;
350 _reduction = reduction;
351 (*this).apply(x,b,res);
352 _reduction = saved_reduction;
353 }
354
355 private:
358 Preconditioner<X,X>& _prec;
359 ScalarProduct<X>& _sp;
360 real_type _reduction;
361 int _maxit;
362 int _verbose;
363 };
364
365
366
368 template<class X>
369 class CGSolver : public InverseOperator<X,X> {
370 public:
372 typedef X domain_type;
374 typedef X range_type;
376 typedef typename X::field_type field_type;
378 typedef typename FieldTraits<field_type>::real_type real_type;
379
385 template<class L, class P>
386 CGSolver (L& op, P& prec, real_type reduction, int maxit, int verbose) :
387 ssp(), _op(op), _prec(prec), _sp(ssp), _reduction(reduction), _maxit(maxit), _verbose(verbose)
388 {
389 static_assert(static_cast<int>(L::category) == static_cast<int>(P::category),
390 "L and P must have the same category!");
391 static_assert(static_cast<int>(L::category) == static_cast<int>(SolverCategory::sequential),
392 "L must be sequential!");
393 }
399 template<class L, class S, class P>
400 CGSolver (L& op, S& sp, P& prec, real_type reduction, int maxit, int verbose) :
401 _op(op), _prec(prec), _sp(sp), _reduction(reduction), _maxit(maxit), _verbose(verbose)
402 {
403 static_assert(static_cast<int>(L::category) == static_cast<int>(P::category),
404 "L and P must have the same category!");
405 static_assert(static_cast<int>(L::category) == static_cast<int>(S::category),
406 "L and S must have the same category!");
407 }
408
414 virtual void apply (X& x, X& b, InverseOperatorResult& res)
415 {
416 res.clear(); // clear solver statistics
417 Timer watch; // start a timer
418 _prec.pre(x,b); // prepare preconditioner
419 _op.applyscaleadd(-1,x,b); // overwrite b with defect
420
421 X p(x); // the search direction
422 X q(x); // a temporary vector
423
424 real_type def0 = _sp.norm(b); // compute norm
425 if (def0<1E-30) // convergence check
426 {
427 res.converged = true;
428 res.iterations = 0; // fill statistics
429 res.reduction = 0;
430 res.conv_rate = 0;
431 res.elapsed=0;
432 if (_verbose>0) // final print
433 std::cout << "=== rate=" << res.conv_rate
434 << ", T=" << res.elapsed << ", TIT=" << res.elapsed
435 << ", IT=0" << std::endl;
436 return;
437 }
438
439 if (_verbose>0) // printing
440 {
441 std::cout << "=== CGSolver" << std::endl;
442 if (_verbose>1) {
443 this->printHeader(std::cout);
444 this->printOutput(std::cout,real_type(0),def0);
445 }
446 }
447
448 // some local variables
449 real_type def=def0; // loop variables
450 field_type rho,rholast,lambda,alpha,beta;
451
452 // determine initial search direction
453 p = 0; // clear correction
454 _prec.apply(p,b); // apply preconditioner
455 rholast = _sp.dot(p,b); // orthogonalization
456
457 // the loop
458 int i=1;
459 for ( ; i<=_maxit; i++ )
460 {
461 // minimize in given search direction p
462 _op.apply(p,q); // q=Ap
463 alpha = _sp.dot(p,q); // scalar product
464 lambda = rholast/alpha; // minimization
465 x.axpy(lambda,p); // update solution
466 b.axpy(-lambda,q); // update defect
467
468 // convergence test
469 real_type defnew=_sp.norm(b); // comp defect norm
470
471 if (_verbose>1) // print
472 this->printOutput(std::cout,real_type(i),defnew,def);
473
474 def = defnew; // update norm
475 if (def<def0*_reduction || def<1E-30) // convergence check
476 {
477 res.converged = true;
478 break;
479 }
480
481 // determine new search direction
482 q = 0; // clear correction
483 _prec.apply(q,b); // apply preconditioner
484 rho = _sp.dot(q,b); // orthogonalization
485 beta = rho/rholast; // scaling factor
486 p *= beta; // scale old search direction
487 p += q; // orthogonalization with correction
488 rholast = rho; // remember rho for recurrence
489 }
490
491 //correct i which is wrong if convergence was not achieved.
492 i=std::min(_maxit,i);
493
494 if (_verbose==1) // printing for non verbose
495 this->printOutput(std::cout,real_type(i),def);
496
497 _prec.post(x); // postprocess preconditioner
498 res.iterations = i; // fill statistics
499 res.reduction = static_cast<double>(def/def0);
500 res.conv_rate = static_cast<double>(pow(res.reduction,1.0/i));
501 res.elapsed = watch.elapsed();
502
503 if (_verbose>0) // final print
504 {
505 std::cout << "=== rate=" << res.conv_rate
506 << ", T=" << res.elapsed
507 << ", TIT=" << res.elapsed/i
508 << ", IT=" << i << std::endl;
509 }
510 }
511
517 virtual void apply (X& x, X& b, double reduction,
519 {
520 real_type saved_reduction = _reduction;
521 _reduction = reduction;
522 (*this).apply(x,b,res);
523 _reduction = saved_reduction;
524 }
525
526 private:
529 Preconditioner<X,X>& _prec;
530 ScalarProduct<X>& _sp;
531 real_type _reduction;
532 int _maxit;
533 int _verbose;
534 };
535
536
537 // Ronald Kriemanns BiCG-STAB implementation from Sumo
539 template<class X>
540 class BiCGSTABSolver : public InverseOperator<X,X> {
541 public:
543 typedef X domain_type;
545 typedef X range_type;
547 typedef typename X::field_type field_type;
549 typedef typename FieldTraits<field_type>::real_type real_type;
550
556 template<class L, class P>
557 BiCGSTABSolver (L& op, P& prec,
558 real_type reduction, int maxit, int verbose) :
559 ssp(), _op(op), _prec(prec), _sp(ssp), _reduction(reduction), _maxit(maxit), _verbose(verbose)
560 {
561 static_assert(static_cast<int>(L::category) == static_cast<int>(P::category),
562 "L and P must be of the same category!");
563 static_assert(static_cast<int>(L::category) == static_cast<int>(SolverCategory::sequential),
564 "L must be sequential!");
565 }
571 template<class L, class S, class P>
572 BiCGSTABSolver (L& op, S& sp, P& prec,
573 real_type reduction, int maxit, int verbose) :
574 _op(op), _prec(prec), _sp(sp), _reduction(reduction), _maxit(maxit), _verbose(verbose)
575 {
576 static_assert(static_cast<int>(L::category) == static_cast<int>(P::category),
577 "L and P must have the same category!");
578 static_assert(static_cast<int>(L::category) == static_cast<int>(S::category),
579 "L and S must have the same category!");
580 }
581
587 virtual void apply (X& x, X& b, InverseOperatorResult& res)
588 {
589 using std::abs;
590 const real_type EPSILON=1e-80;
591 double it;
592 field_type rho, rho_new, alpha, beta, h, omega;
593 real_type norm, norm_old, norm_0;
594
595 //
596 // get vectors and matrix
597 //
598 X& r=b;
599 X p(x);
600 X v(x);
601 X t(x);
602 X y(x);
603 X rt(x);
604
605 //
606 // begin iteration
607 //
608
609 // r = r - Ax; rt = r
610 res.clear(); // clear solver statistics
611 Timer watch; // start a timer
612 _prec.pre(x,r); // prepare preconditioner
613 _op.applyscaleadd(-1,x,r); // overwrite b with defect
614
615 rt=r;
616
617 norm = norm_old = norm_0 = _sp.norm(r);
618
619 p=0;
620 v=0;
621
622 rho = 1;
623 alpha = 1;
624 omega = 1;
625
626 if (_verbose>0) // printing
627 {
628 std::cout << "=== BiCGSTABSolver" << std::endl;
629 if (_verbose>1)
630 {
631 this->printHeader(std::cout);
632 this->printOutput(std::cout,0,norm_0);
633 //std::cout << " Iter Defect Rate" << std::endl;
634 //std::cout << " 0" << std::setw(14) << norm_0 << std::endl;
635 }
636 }
637
638 if ( norm < (_reduction * norm_0) || norm<1E-30)
639 {
640 res.converged = 1;
641 _prec.post(x); // postprocess preconditioner
642 res.iterations = 0; // fill statistics
643 res.reduction = 0;
644 res.conv_rate = 0;
645 res.elapsed = watch.elapsed();
646 return;
647 }
648
649 //
650 // iteration
651 //
652
653 for (it = 0.5; it < _maxit; it+=.5)
654 {
655 //
656 // preprocess, set vecsizes etc.
657 //
658
659 // rho_new = < rt , r >
660 rho_new = _sp.dot(rt,r);
661
662 // look if breakdown occured
663 if (abs(rho) <= EPSILON)
664 DUNE_THROW(ISTLError,"breakdown in BiCGSTAB - rho "
665 << rho << " <= EPSILON " << EPSILON
666 << " after " << it << " iterations");
667 if (abs(omega) <= EPSILON)
668 DUNE_THROW(ISTLError,"breakdown in BiCGSTAB - omega "
669 << omega << " <= EPSILON " << EPSILON
670 << " after " << it << " iterations");
671
672
673 if (it<1)
674 p = r;
675 else
676 {
677 beta = ( rho_new / rho ) * ( alpha / omega );
678 p.axpy(-omega,v); // p = r + beta (p - omega*v)
679 p *= beta;
680 p += r;
681 }
682
683 // y = W^-1 * p
684 y = 0;
685 _prec.apply(y,p); // apply preconditioner
686
687 // v = A * y
688 _op.apply(y,v);
689
690 // alpha = rho_new / < rt, v >
691 h = _sp.dot(rt,v);
692
693 if (abs(h) < EPSILON)
694 DUNE_THROW(ISTLError,"abs(h) < EPSILON in BiCGSTAB - abs(h) "
695 << abs(h) << " < EPSILON " << EPSILON
696 << " after " << it << " iterations");
697
698 alpha = rho_new / h;
699
700 // apply first correction to x
701 // x <- x + alpha y
702 x.axpy(alpha,y);
703
704 // r = r - alpha*v
705 r.axpy(-alpha,v);
706
707 //
708 // test stop criteria
709 //
710
711 norm = _sp.norm(r);
712
713 if (_verbose>1) // print
714 {
715 this->printOutput(std::cout,it,norm,norm_old);
716 }
717
718 if ( norm < (_reduction * norm_0) )
719 {
720 res.converged = 1;
721 break;
722 }
723 it+=.5;
724
725 norm_old = norm;
726
727 // y = W^-1 * r
728 y = 0;
729 _prec.apply(y,r);
730
731 // t = A * y
732 _op.apply(y,t);
733
734 // omega = < t, r > / < t, t >
735 omega = _sp.dot(t,r)/_sp.dot(t,t);
736
737 // apply second correction to x
738 // x <- x + omega y
739 x.axpy(omega,y);
740
741 // r = s - omega*t (remember : r = s)
742 r.axpy(-omega,t);
743
744 rho = rho_new;
745
746 //
747 // test stop criteria
748 //
749
750 norm = _sp.norm(r);
751
752 if (_verbose > 1) // print
753 {
754 this->printOutput(std::cout,it,norm,norm_old);
755 }
756
757 if ( norm < (_reduction * norm_0) || norm<1E-30)
758 {
759 res.converged = 1;
760 break;
761 }
762
763 norm_old = norm;
764 } // end for
765
766 //correct i which is wrong if convergence was not achieved.
767 it=std::min(static_cast<double>(_maxit),it);
768
769 if (_verbose==1) // printing for non verbose
770 this->printOutput(std::cout,it,norm);
771
772 _prec.post(x); // postprocess preconditioner
773 res.iterations = static_cast<int>(std::ceil(it)); // fill statistics
774 res.reduction = static_cast<double>(norm/norm_0);
775 res.conv_rate = static_cast<double>(pow(res.reduction,1.0/it));
776 res.elapsed = watch.elapsed();
777 if (_verbose>0) // final print
778 std::cout << "=== rate=" << res.conv_rate
779 << ", T=" << res.elapsed
780 << ", TIT=" << res.elapsed/it
781 << ", IT=" << it << std::endl;
782 }
783
789 virtual void apply (X& x, X& b, double reduction, InverseOperatorResult& res)
790 {
791 real_type saved_reduction = _reduction;
792 _reduction = reduction;
793 (*this).apply(x,b,res);
794 _reduction = saved_reduction;
795 }
796
797 private:
800 Preconditioner<X,X>& _prec;
801 ScalarProduct<X>& _sp;
802 real_type _reduction;
803 int _maxit;
804 int _verbose;
805 };
806
813 template<class X>
814 class MINRESSolver : public InverseOperator<X,X> {
815 public:
817 typedef X domain_type;
819 typedef X range_type;
821 typedef typename X::field_type field_type;
823 typedef typename FieldTraits<field_type>::real_type real_type;
824
830 template<class L, class P>
831 MINRESSolver (L& op, P& prec, real_type reduction, int maxit, int verbose) :
832 ssp(), _op(op), _prec(prec), _sp(ssp), _reduction(reduction), _maxit(maxit), _verbose(verbose)
833 {
834 static_assert(static_cast<int>(L::category) == static_cast<int>(P::category),
835 "L and P must have the same category!");
836 static_assert(static_cast<int>(L::category) == static_cast<int>(SolverCategory::sequential),
837 "L must be sequential!");
838 }
844 template<class L, class S, class P>
845 MINRESSolver (L& op, S& sp, P& prec, real_type reduction, int maxit, int verbose) :
846 _op(op), _prec(prec), _sp(sp), _reduction(reduction), _maxit(maxit), _verbose(verbose)
847 {
848 static_assert(static_cast<int>(L::category) == static_cast<int>(P::category),
849 "L and P must have the same category!");
850 static_assert(static_cast<int>(L::category) == static_cast<int>(S::category),
851 "L and S must have the same category!");
852 }
853
859 virtual void apply (X& x, X& b, InverseOperatorResult& res)
860 {
861 using std::sqrt;
862 using std::abs;
863 // clear solver statistics
864 res.clear();
865 // start a timer
866 Dune::Timer watch;
867 watch.reset();
868 // prepare preconditioner
869 _prec.pre(x,b);
870 // overwrite rhs with defect
871 _op.applyscaleadd(-1,x,b);
872
873 // compute residual norm
874 real_type def0 = _sp.norm(b);
875
876 // printing
877 if(_verbose > 0) {
878 std::cout << "=== MINRESSolver" << std::endl;
879 if(_verbose > 1) {
880 this->printHeader(std::cout);
881 this->printOutput(std::cout,0,def0);
882 }
883 }
884
885 // check for convergence
886 if(def0 < 1e-30 ) {
887 res.converged = true;
888 res.iterations = 0;
889 res.reduction = 0;
890 res.conv_rate = 0;
891 res.elapsed = 0.0;
892 // final print
893 if(_verbose > 0)
894 std::cout << "=== rate=" << res.conv_rate
895 << ", T=" << res.elapsed
896 << ", TIT=" << res.elapsed
897 << ", IT=" << res.iterations
898 << std::endl;
899 return;
900 }
901
902 // the defect norm
903 real_type def = def0;
904 // recurrence coefficients as computed in Lanczos algorithm
905 field_type alpha, beta;
906 // diagonal entries of givens rotation
907 Dune::array<real_type,2> c{{0.0,0.0}};
908 // off-diagonal entries of givens rotation
909 Dune::array<field_type,2> s{{0.0,0.0}};
910
911 // recurrence coefficients (column k of tridiag matrix T_k)
912 Dune::array<field_type,3> T{{0.0,0.0,0.0}};
913
914 // the rhs vector of the min problem
915 Dune::array<field_type,2> xi{{1.0,0.0}};
916
917 // some temporary vectors
918 X z(b), dummy(b);
919
920 // initialize and clear correction
921 z = 0.0;
922 _prec.apply(z,b);
923
924 // beta is real and positive in exact arithmetic
925 // since it is the norm of the basis vectors (in unpreconditioned case)
926 beta = sqrt(_sp.dot(b,z));
927 field_type beta0 = beta;
928
929 // the search directions
930 Dune::array<X,3> p{{b,b,b}};
931 p[0] = 0.0;
932 p[1] = 0.0;
933 p[2] = 0.0;
934
935 // orthonormal basis vectors (in unpreconditioned case)
936 Dune::array<X,3> q{{b,b,b}};
937 q[0] = 0.0;
938 q[1] *= 1.0/beta;
939 q[2] = 0.0;
940
941 z *= 1.0/beta;
942
943 // the loop
944 int i = 1;
945 for( ; i<=_maxit; i++) {
946
947 dummy = z;
948 int i1 = i%3,
949 i0 = (i1+2)%3,
950 i2 = (i1+1)%3;
951
952 // symmetrically preconditioned Lanczos algorithm (see Greenbaum p.121)
953 _op.apply(z,q[i2]); // q[i2] = Az
954 q[i2].axpy(-beta,q[i0]);
955 // alpha is real since it is the diagonal entry of the hermitian tridiagonal matrix
956 // from the Lanczos Algorithm
957 // so the order in the scalar product doesn't matter even for the complex case
958 alpha = _sp.dot(z,q[i2]);
959 q[i2].axpy(-alpha,q[i1]);
960
961 z = 0.0;
962 _prec.apply(z,q[i2]);
963
964 // beta is real and positive in exact arithmetic
965 // since it is the norm of the basis vectors (in unpreconditioned case)
966 beta = sqrt(_sp.dot(q[i2],z));
967
968 q[i2] *= 1.0/beta;
969 z *= 1.0/beta;
970
971 // QR Factorization of recurrence coefficient matrix
972 // apply previous givens rotations to last column of T
973 T[1] = T[2];
974 if(i>2) {
975 T[0] = s[i%2]*T[1];
976 T[1] = c[i%2]*T[1];
977 }
978 if(i>1) {
979 T[2] = c[(i+1)%2]*alpha - s[(i+1)%2]*T[1];
980 T[1] = c[(i+1)%2]*T[1] + s[(i+1)%2]*alpha;
981 }
982 else
983 T[2] = alpha;
984
985 // update QR factorization
986 generateGivensRotation(T[2],beta,c[i%2],s[i%2]);
987 // to last column of T_k
988 T[2] = c[i%2]*T[2] + s[i%2]*beta;
989 // and to the rhs xi of the min problem
990 xi[i%2] = -s[i%2]*xi[(i+1)%2];
991 xi[(i+1)%2] *= c[i%2];
992
993 // compute correction direction
994 p[i2] = dummy;
995 p[i2].axpy(-T[1],p[i1]);
996 p[i2].axpy(-T[0],p[i0]);
997 p[i2] *= 1.0/T[2];
998
999 // apply correction/update solution
1000 x.axpy(beta0*xi[(i+1)%2],p[i2]);
1001
1002 // remember beta_old
1003 T[2] = beta;
1004
1005 // check for convergence
1006 // the last entry in the rhs of the min-problem is the residual
1007 real_type defnew = abs(beta0*xi[i%2]);
1008
1009 if(_verbose > 1)
1010 this->printOutput(std::cout,i,defnew,def);
1011
1012 def = defnew;
1013 if(def < def0*_reduction || def < 1e-30 || i == _maxit ) {
1014 res.converged = true;
1015 break;
1016 }
1017 } // end for
1018
1019 if(_verbose == 1)
1020 this->printOutput(std::cout,i,def);
1021
1022 // postprocess preconditioner
1023 _prec.post(x);
1024 // fill statistics
1025 res.iterations = i;
1026 res.reduction = static_cast<double>(def/def0);
1027 res.conv_rate = static_cast<double>(pow(res.reduction,1.0/i));
1028 res.elapsed = watch.elapsed();
1029
1030 // final print
1031 if(_verbose > 0) {
1032 std::cout << "=== rate=" << res.conv_rate
1033 << ", T=" << res.elapsed
1034 << ", TIT=" << res.elapsed/i
1035 << ", IT=" << i << std::endl;
1036 }
1037 }
1038
1044 virtual void apply (X& x, X& b, double reduction, InverseOperatorResult& res)
1045 {
1046 real_type saved_reduction = _reduction;
1047 _reduction = reduction;
1048 (*this).apply(x,b,res);
1049 _reduction = saved_reduction;
1050 }
1051
1052 private:
1053
1054 void generateGivensRotation(field_type &dx, field_type &dy, real_type &cs, field_type &sn)
1055 {
1056 using std::sqrt;
1057 using std::abs;
1058 real_type norm_dx = abs(dx);
1059 real_type norm_dy = abs(dy);
1060 if(norm_dy < 1e-15) {
1061 cs = 1.0;
1062 sn = 0.0;
1063 } else if(norm_dx < 1e-15) {
1064 cs = 0.0;
1065 sn = 1.0;
1066 } else if(norm_dy > norm_dx) {
1067 real_type temp = norm_dx/norm_dy;
1068 cs = 1.0/sqrt(1.0 + temp*temp);
1069 sn = cs;
1070 cs *= temp;
1071 sn *= dx/norm_dx;
1072 // dy is real in exact arithmetic
1073 // so we don't need to conjugate here
1074 sn *= dy/norm_dy;
1075 } else {
1076 real_type temp = norm_dy/norm_dx;
1077 cs = 1.0/sqrt(1.0 + temp*temp);
1078 sn = cs;
1079 sn *= dy/dx;
1080 // dy and dx is real in exact arithmetic
1081 // so we don't have to conjugate both of them
1082 }
1083 }
1084
1085 SeqScalarProduct<X> ssp;
1086 LinearOperator<X,X>& _op;
1087 Preconditioner<X,X>& _prec;
1088 ScalarProduct<X>& _sp;
1089 real_type _reduction;
1090 int _maxit;
1091 int _verbose;
1092 };
1093
1107 template<class X, class Y=X, class F = Y>
1109 {
1110 public:
1112 typedef X domain_type;
1114 typedef Y range_type;
1116 typedef typename X::field_type field_type;
1118 typedef typename FieldTraits<field_type>::real_type real_type;
1120 typedef F basis_type;
1121
1122 template<class L, class P>
1123 DUNE_DEPRECATED_MSG("recalc_defect is a unused parameter! Use RestartedGMResSolver(L& op, P& prec, real_type reduction, int restart, int maxit, int verbose) instead")
1124 RestartedGMResSolver (L& op, P& prec, real_type reduction, int restart, int maxit, int verbose, bool recalc_defect)
1125 : _A(op)
1126 , _W(prec)
1127 , ssp()
1128 , _sp(ssp)
1129 , _restart(restart)
1130 , _reduction(reduction)
1131 , _maxit(maxit)
1132 , _verbose(verbose)
1133 {
1134 static_assert(static_cast<int>(P::category) == static_cast<int>(L::category),
1135 "P and L must be the same category!");
1136 static_assert(static_cast<int>(L::category) == static_cast<int>(SolverCategory::sequential),
1137 "L must be sequential!");
1138 }
1139
1140
1147 template<class L, class P>
1148 RestartedGMResSolver (L& op, P& prec, real_type reduction, int restart, int maxit, int verbose) :
1149 _A(op), _W(prec),
1150 ssp(), _sp(ssp), _restart(restart),
1151 _reduction(reduction), _maxit(maxit), _verbose(verbose)
1152 {
1153 static_assert(static_cast<int>(P::category) == static_cast<int>(L::category),
1154 "P and L must be the same category!");
1155 static_assert(static_cast<int>(L::category) == static_cast<int>(SolverCategory::sequential),
1156 "L must be sequential!");
1157 }
1158
1159 template<class L, class S, class P>
1160 DUNE_DEPRECATED_MSG("recalc_defect is a unused parameter! Use RestartedGMResSolver(L& op, S& sp, P& prec, real_type reduction, int restart, int maxit, int verbose) instead")
1161 RestartedGMResSolver(L& op, S& sp, P& prec, real_type reduction, int restart, int maxit, int verbose, bool recalc_defect)
1162 : _A(op)
1163 , _W(prec)
1164 , _sp(sp)
1165 , _restart(restart)
1166 , _reduction(reduction)
1167 , _maxit(maxit)
1168 , _verbose(verbose)
1169 {
1170 static_assert(static_cast<int>(P::category) == static_cast<int>(L::category),
1171 " P and L must have the same category!");
1172 static_assert(static_cast<int>(P::category) == static_cast<int>(S::category),
1173 "P and S must have the same category!");
1174 }
1175
1182 template<class L, class S, class P>
1183 RestartedGMResSolver (L& op, S& sp, P& prec, real_type reduction, int restart, int maxit, int verbose) :
1184 _A(op), _W(prec),
1185 _sp(sp), _restart(restart),
1186 _reduction(reduction), _maxit(maxit), _verbose(verbose)
1187 {
1188 static_assert(static_cast<int>(P::category) == static_cast<int>(L::category),
1189 "P and L must have the same category!");
1190 static_assert(static_cast<int>(P::category) == static_cast<int>(S::category),
1191 "P and S must have the same category!");
1192 }
1193
1195 virtual void apply (X& x, Y& b, InverseOperatorResult& res)
1196 {
1197 apply(x,b,_reduction,res);
1198 }
1199
1205 virtual void apply (X& x, Y& b, double reduction, InverseOperatorResult& res)
1206 {
1207 using std::abs;
1208 const real_type EPSILON = 1e-80;
1209 const int m = _restart;
1210 real_type norm, norm_old = 0.0, norm_0;
1211 int j = 1;
1212 std::vector<field_type> s(m+1), sn(m);
1213 std::vector<real_type> cs(m);
1214 // need copy of rhs if GMRes has to be restarted
1215 Y b2(b);
1216 // helper vector
1217 Y w(b);
1218 std::vector< std::vector<field_type> > H(m+1,s);
1219 std::vector<F> v(m+1,b);
1220
1221 // start timer
1222 Dune::Timer watch;
1223 watch.reset();
1224
1225 // clear solver statistics and set res.converged to false
1226 res.clear();
1227 _W.pre(x,b);
1228
1229 // calculate defect and overwrite rhs with it
1230 _A.applyscaleadd(-1.0,x,b); // b -= Ax
1231 // calculate preconditioned defect
1232 v[0] = 0.0; _W.apply(v[0],b); // r = W^-1 b
1233 norm_0 = _sp.norm(v[0]);
1234 norm = norm_0;
1235 norm_old = norm;
1236
1237 // print header
1238 if(_verbose > 0)
1239 {
1240 std::cout << "=== RestartedGMResSolver" << std::endl;
1241 if(_verbose > 1) {
1242 this->printHeader(std::cout);
1243 this->printOutput(std::cout,0,norm_0);
1244 }
1245 }
1246
1247 if(norm_0 < EPSILON) {
1248 _W.post(x);
1249 res.converged = true;
1250 if(_verbose > 0) // final print
1251 print_result(res);
1252 }
1253
1254 while(j <= _maxit && res.converged != true) {
1255
1256 int i = 0;
1257 v[0] *= 1.0/norm;
1258 s[0] = norm;
1259 for(i=1; i<m+1; i++)
1260 s[i] = 0.0;
1261
1262 for(i=0; i < m && j <= _maxit && res.converged != true; i++, j++) {
1263 w = 0.0;
1264 // use v[i+1] as temporary vector
1265 v[i+1] = 0.0;
1266 // do Arnoldi algorithm
1267 _A.apply(v[i],v[i+1]);
1268 _W.apply(w,v[i+1]);
1269 for(int k=0; k<i+1; k++) {
1270 // notice that _sp.dot(v[k],w) = v[k]\adjoint w
1271 // so one has to pay attention to the order
1272 // in the scalar product for the complex case
1273 // doing the modified Gram-Schmidt algorithm
1274 H[k][i] = _sp.dot(v[k],w);
1275 // w -= H[k][i] * v[k]
1276 w.axpy(-H[k][i],v[k]);
1277 }
1278 H[i+1][i] = _sp.norm(w);
1279 if(abs(H[i+1][i]) < EPSILON)
1281 "breakdown in GMRes - |w| == 0.0 after " << j << " iterations");
1282
1283 // normalize new vector
1284 v[i+1] = w; v[i+1] *= 1.0/H[i+1][i];
1285
1286 // update QR factorization
1287 for(int k=0; k<i; k++)
1288 applyPlaneRotation(H[k][i],H[k+1][i],cs[k],sn[k]);
1289
1290 // compute new givens rotation
1291 generatePlaneRotation(H[i][i],H[i+1][i],cs[i],sn[i]);
1292 // finish updating QR factorization
1293 applyPlaneRotation(H[i][i],H[i+1][i],cs[i],sn[i]);
1294 applyPlaneRotation(s[i],s[i+1],cs[i],sn[i]);
1295
1296 // norm of the defect is the last component the vector s
1297 norm = abs(s[i+1]);
1298
1299 // print current iteration statistics
1300 if(_verbose > 1) {
1301 this->printOutput(std::cout,j,norm,norm_old);
1302 }
1303
1304 norm_old = norm;
1305
1306 // check convergence
1307 if(norm < reduction * norm_0)
1308 res.converged = true;
1309
1310 } // end for
1311
1312 // calculate update vector
1313 w = 0.0;
1314 update(w,i,H,s,v);
1315 // and current iterate
1316 x += w;
1317
1318 // restart GMRes if convergence was not achieved,
1319 // i.e. linear defect has not reached desired reduction
1320 // and if j < _maxit
1321 if( res.converged != true && j <= _maxit ) {
1322
1323 if(_verbose > 0)
1324 std::cout << "=== GMRes::restart" << std::endl;
1325 // get saved rhs
1326 b = b2;
1327 // calculate new defect
1328 _A.applyscaleadd(-1.0,x,b); // b -= Ax;
1329 // calculate preconditioned defect
1330 v[0] = 0.0;
1331 _W.apply(v[0],b);
1332 norm = _sp.norm(v[0]);
1333 norm_old = norm;
1334 }
1335
1336 } //end while
1337
1338 // postprocess preconditioner
1339 _W.post(x);
1340
1341 // save solver statistics
1342 res.iterations = j-1; // it has to be j-1!!!
1343 res.reduction = static_cast<double>(norm/norm_0);
1344 res.conv_rate = static_cast<double>(pow(res.reduction,1.0/(j-1)));
1345 res.elapsed = watch.elapsed();
1346
1347 if(_verbose>0)
1348 print_result(res);
1349
1350 }
1351
1352 private :
1353
1354 void print_result(const InverseOperatorResult& res) const {
1355 int k = res.iterations>0 ? res.iterations : 1;
1356 std::cout << "=== rate=" << res.conv_rate
1357 << ", T=" << res.elapsed
1358 << ", TIT=" << res.elapsed/k
1359 << ", IT=" << res.iterations
1360 << std::endl;
1361 }
1362
1363 void update(X& w, int i,
1364 const std::vector<std::vector<field_type> >& H,
1365 const std::vector<field_type>& s,
1366 const std::vector<X>& v) {
1367 // solution vector of the upper triangular system
1368 std::vector<field_type> y(s);
1369
1370 // backsolve
1371 for(int a=i-1; a>=0; a--) {
1372 field_type rhs(s[a]);
1373 for(int b=a+1; b<i; b++)
1374 rhs -= H[a][b]*y[b];
1375 y[a] = rhs/H[a][a];
1376
1377 // compute update on the fly
1378 // w += y[a]*v[a]
1379 w.axpy(y[a],v[a]);
1380 }
1381 }
1382
1383 template<typename T>
1384 typename enable_if<is_same<field_type,real_type>::value,T>::type conjugate(const T& t) {
1385 return t;
1386 }
1387
1388 template<typename T>
1389 typename enable_if<!is_same<field_type,real_type>::value,T>::type conjugate(const T& t) {
1390 return conj(t);
1391 }
1392
1393 void
1394 generatePlaneRotation(field_type &dx, field_type &dy, real_type &cs, field_type &sn)
1395 {
1396 using std::sqrt;
1397 using std::abs;
1398 real_type norm_dx = abs(dx);
1399 real_type norm_dy = abs(dy);
1400 if(norm_dy < 1e-15) {
1401 cs = 1.0;
1402 sn = 0.0;
1403 } else if(norm_dx < 1e-15) {
1404 cs = 0.0;
1405 sn = 1.0;
1406 } else if(norm_dy > norm_dx) {
1407 real_type temp = norm_dx/norm_dy;
1408 cs = 1.0/sqrt(1.0 + temp*temp);
1409 sn = cs;
1410 cs *= temp;
1411 sn *= dx/norm_dx;
1412 sn *= conjugate(dy)/norm_dy;
1413 } else {
1414 real_type temp = norm_dy/norm_dx;
1415 cs = 1.0/sqrt(1.0 + temp*temp);
1416 sn = cs;
1417 sn *= conjugate(dy/dx);
1418 }
1419 }
1420
1421
1422 void
1423 applyPlaneRotation(field_type &dx, field_type &dy, real_type &cs, field_type &sn)
1424 {
1425 field_type temp = cs * dx + sn * dy;
1426 dy = -conjugate(sn) * dx + cs * dy;
1427 dx = temp;
1428 }
1429
1430 LinearOperator<X,Y>& _A;
1431 Preconditioner<X,Y>& _W;
1432 SeqScalarProduct<X> ssp;
1433 ScalarProduct<X>& _sp;
1434 int _restart;
1435 real_type _reduction;
1436 int _maxit;
1437 int _verbose;
1438 };
1439
1440
1454 template<class X>
1456 {
1457 public:
1459 typedef X domain_type;
1461 typedef X range_type;
1463 typedef typename X::field_type field_type;
1465 typedef typename FieldTraits<field_type>::real_type real_type;
1466
1474 template<class L, class P>
1475 GeneralizedPCGSolver (L& op, P& prec, real_type reduction, int maxit, int verbose,
1476 int restart=10) :
1477 ssp(), _op(op), _prec(prec), _sp(ssp), _reduction(reduction), _maxit(maxit),
1478 _verbose(verbose), _restart(std::min(maxit,restart))
1479 {
1480 static_assert(static_cast<int>(L::category) == static_cast<int>(P::category),
1481 "L and P have to have the same category!");
1482 static_assert(static_cast<int>(L::category) ==
1483 static_cast<int>(SolverCategory::sequential),
1484 "L has to be sequential!");
1485 }
1493 template<class L, class P, class S>
1494 GeneralizedPCGSolver (L& op, S& sp, P& prec,
1495 real_type reduction, int maxit, int verbose, int restart=10) :
1496 _op(op), _prec(prec), _sp(sp), _reduction(reduction), _maxit(maxit), _verbose(verbose),
1497 _restart(std::min(maxit,restart))
1498 {
1499 static_assert(static_cast<int>(L::category) == static_cast<int>(P::category),
1500 "L and P must have the same category!");
1501 static_assert(static_cast<int>(L::category) == static_cast<int>(S::category),
1502 "L and S must have the same category!");
1503 }
1509 virtual void apply (X& x, X& b, InverseOperatorResult& res)
1510 {
1511 res.clear(); // clear solver statistics
1512 Timer watch; // start a timer
1513 _prec.pre(x,b); // prepare preconditioner
1514 _op.applyscaleadd(-1,x,b); // overwrite b with defect
1515
1516 std::vector<std::shared_ptr<X> > p(_restart);
1517 std::vector<typename X::field_type> pp(_restart);
1518 X q(x); // a temporary vector
1519 X prec_res(x); // a temporary vector for preconditioner output
1520
1521 p[0].reset(new X(x));
1522
1523 real_type def0 = _sp.norm(b); // compute norm
1524 if (def0<1E-30) // convergence check
1525 {
1526 res.converged = true;
1527 res.iterations = 0; // fill statistics
1528 res.reduction = 0;
1529 res.conv_rate = 0;
1530 res.elapsed=0;
1531 if (_verbose>0) // final print
1532 std::cout << "=== rate=" << res.conv_rate
1533 << ", T=" << res.elapsed << ", TIT=" << res.elapsed
1534 << ", IT=0" << std::endl;
1535 return;
1536 }
1537
1538 if (_verbose>0) // printing
1539 {
1540 std::cout << "=== GeneralizedPCGSolver" << std::endl;
1541 if (_verbose>1) {
1542 this->printHeader(std::cout);
1543 this->printOutput(std::cout,0,def0);
1544 }
1545 }
1546 // some local variables
1547 real_type def=def0; // loop variables
1548 field_type rho, lambda;
1549
1550 int i=0;
1551 int ii=0;
1552 // determine initial search direction
1553 *(p[0]) = 0; // clear correction
1554 _prec.apply(*(p[0]),b); // apply preconditioner
1555 rho = _sp.dot(*(p[0]),b); // orthogonalization
1556 _op.apply(*(p[0]),q); // q=Ap
1557 pp[0] = _sp.dot(*(p[0]),q); // scalar product
1558 lambda = rho/pp[0]; // minimization
1559 x.axpy(lambda,*(p[0])); // update solution
1560 b.axpy(-lambda,q); // update defect
1561
1562 // convergence test
1563 real_type defnew=_sp.norm(b); // comp defect norm
1564 if (_verbose>1) // print
1565 this->printOutput(std::cout,++i,defnew,def);
1566 def = defnew; // update norm
1567 if (def<def0*_reduction || def<1E-30) // convergence check
1568 {
1569 res.converged = true;
1570 if (_verbose>0) // final print
1571 {
1572 std::cout << "=== rate=" << res.conv_rate
1573 << ", T=" << res.elapsed
1574 << ", TIT=" << res.elapsed
1575 << ", IT=" << 1 << std::endl;
1576 }
1577 return;
1578 }
1579
1580 while(i<_maxit) {
1581 // the loop
1582 int end=std::min(_restart, _maxit-i+1);
1583 for (ii=1; ii<end; ++ii )
1584 {
1585 //std::cout<<" ii="<<ii<<" i="<<i<<std::endl;
1586 // compute next conjugate direction
1587 prec_res = 0; // clear correction
1588 _prec.apply(prec_res,b); // apply preconditioner
1589
1590 p[ii].reset(new X(prec_res));
1591 _op.apply(prec_res, q);
1592
1593 for(int j=0; j<ii; ++j) {
1594 rho =_sp.dot(q,*(p[j]))/pp[j];
1595 p[ii]->axpy(-rho, *(p[j]));
1596 }
1597
1598 // minimize in given search direction
1599 _op.apply(*(p[ii]),q); // q=Ap
1600 pp[ii] = _sp.dot(*(p[ii]),q); // scalar product
1601 rho = _sp.dot(*(p[ii]),b); // orthogonalization
1602 lambda = rho/pp[ii]; // minimization
1603 x.axpy(lambda,*(p[ii])); // update solution
1604 b.axpy(-lambda,q); // update defect
1605
1606 // convergence test
1607 real_type defNew=_sp.norm(b); // comp defect norm
1608
1609 if (_verbose>1) // print
1610 this->printOutput(std::cout,++i,defNew,def);
1611
1612 def = defNew; // update norm
1613 if (def<def0*_reduction || def<1E-30) // convergence check
1614 {
1615 res.converged = true;
1616 break;
1617 }
1618 }
1619 if(res.converged)
1620 break;
1621 if(end==_restart) {
1622 *(p[0])=*(p[_restart-1]);
1623 pp[0]=pp[_restart-1];
1624 }
1625 }
1626
1627 // postprocess preconditioner
1628 _prec.post(x);
1629
1630 // fill statistics
1631 res.iterations = i;
1632 res.reduction = def/def0;
1633 res.conv_rate = pow(res.reduction,1.0/i);
1634 res.elapsed = watch.elapsed();
1635
1636 if (_verbose>0) // final print
1637 {
1638 std::cout << "=== rate=" << res.conv_rate
1639 << ", T=" << res.elapsed
1640 << ", TIT=" << res.elapsed/i
1641 << ", IT=" << i+1 << std::endl;
1642 }
1643 }
1644
1650 virtual void apply (X& x, X& b, double reduction,
1652 {
1653 real_type saved_reduction = _reduction;
1654 _reduction = reduction;
1655 (*this).apply(x,b,res);
1656 _reduction = saved_reduction;
1657 }
1658 private:
1661 Preconditioner<X,X>& _prec;
1662 ScalarProduct<X>& _sp;
1663 real_type _reduction;
1664 int _maxit;
1665 int _verbose;
1666 int _restart;
1667 };
1668
1671} // end namespace
1672
1673#endif
Fallback implementation of the std::array class (a static array)
Bi-conjugate Gradient Stabilized (BiCG-STAB)
Definition: solvers.hh:540
BiCGSTABSolver(L &op, S &sp, P &prec, real_type reduction, int maxit, int verbose)
Set up solver.
Definition: solvers.hh:572
X domain_type
The domain type of the operator to be inverted.
Definition: solvers.hh:543
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: solvers.hh:549
X range_type
The range type of the operator to be inverted.
Definition: solvers.hh:545
virtual void apply(X &x, X &b, InverseOperatorResult &res)
Apply inverse operator.
Definition: solvers.hh:587
X::field_type field_type
The field type of the operator to be inverted.
Definition: solvers.hh:547
virtual void apply(X &x, X &b, double reduction, InverseOperatorResult &res)
Apply inverse operator with given reduction factor.
Definition: solvers.hh:789
BiCGSTABSolver(L &op, P &prec, real_type reduction, int maxit, int verbose)
Set up solver.
Definition: solvers.hh:557
conjugate gradient method
Definition: solvers.hh:369
X range_type
The range type of the operator to be inverted.
Definition: solvers.hh:374
X domain_type
The domain type of the operator to be inverted.
Definition: solvers.hh:372
CGSolver(L &op, S &sp, P &prec, real_type reduction, int maxit, int verbose)
Set up conjugate gradient solver.
Definition: solvers.hh:400
virtual void apply(X &x, X &b, InverseOperatorResult &res)
Apply inverse operator.
Definition: solvers.hh:414
virtual void apply(X &x, X &b, double reduction, InverseOperatorResult &res)
Apply inverse operator with given reduction factor.
Definition: solvers.hh:517
CGSolver(L &op, P &prec, real_type reduction, int maxit, int verbose)
Set up conjugate gradient solver.
Definition: solvers.hh:386
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: solvers.hh:378
X::field_type field_type
The field type of the operator to be inverted.
Definition: solvers.hh:376
Generalized preconditioned conjugate gradient solver.
Definition: solvers.hh:1456
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: solvers.hh:1465
GeneralizedPCGSolver(L &op, P &prec, real_type reduction, int maxit, int verbose, int restart=10)
Set up nonlinear preconditioned conjugate gradient solver.
Definition: solvers.hh:1475
virtual void apply(X &x, X &b, double reduction, InverseOperatorResult &res)
Apply inverse operator with given reduction factor.
Definition: solvers.hh:1650
virtual void apply(X &x, X &b, InverseOperatorResult &res)
Apply inverse operator.
Definition: solvers.hh:1509
GeneralizedPCGSolver(L &op, S &sp, P &prec, real_type reduction, int maxit, int verbose, int restart=10)
Set up nonlinear preconditioned conjugate gradient solver.
Definition: solvers.hh:1494
X::field_type field_type
The field type of the operator to be inverted.
Definition: solvers.hh:1463
X range_type
The range type of the operator to be inverted.
Definition: solvers.hh:1461
X domain_type
The domain type of the operator to be inverted.
Definition: solvers.hh:1459
gradient method
Definition: solvers.hh:231
X domain_type
The domain type of the operator that we do the inverse for.
Definition: solvers.hh:234
virtual void apply(X &x, X &b, InverseOperatorResult &res)
Apply inverse operator.
Definition: solvers.hh:279
GradientSolver(L &op, P &prec, real_type reduction, int maxit, int verbose)
Set up solver.
Definition: solvers.hh:249
GradientSolver(L &op, S &sp, P &prec, real_type reduction, int maxit, int verbose)
Set up solver.
Definition: solvers.hh:264
X range_type
The range type of the operator that we do the inverse for.
Definition: solvers.hh:236
X::field_type field_type
The field type of the operator that we do the inverse for.
Definition: solvers.hh:238
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: solvers.hh:240
virtual void apply(X &x, X &b, double reduction, InverseOperatorResult &res)
Apply inverse operator with given reduction factor.
Definition: solvers.hh:347
derive error class from the base class in common
Definition: istlexception.hh:16
Abstract base class for all solvers.
Definition: solver.hh:79
void printHeader(std::ostream &s) const
helper function for printing header of solver output
Definition: solver.hh:121
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:130
virtual void applyscaleadd(field_type alpha, const X &x, Y &y) const =0
apply operator to x, scale and add:
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...
Preconditioned loop solver.
Definition: solvers.hh:57
X range_type
The range type of the operator that we do the inverse for.
Definition: solvers.hh:62
X domain_type
The domain type of the operator that we do the inverse for.
Definition: solvers.hh:60
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: solvers.hh:66
LoopSolver(L &op, P &prec, real_type reduction, int maxit, int verbose)
Set up Loop solver.
Definition: solvers.hh:88
virtual void apply(X &x, X &b, InverseOperatorResult &res)
Apply inverse operator,.
Definition: solvers.hh:131
virtual void apply(X &x, X &b, double reduction, InverseOperatorResult &res)
apply inverse operator, with given convergence criteria.
Definition: solvers.hh:209
LoopSolver(L &op, S &sp, P &prec, real_type reduction, int maxit, int verbose)
Set up loop solver.
Definition: solvers.hh:119
X::field_type field_type
The field type of the operator that we do the inverse for.
Definition: solvers.hh:64
Minimal Residual Method (MINRES)
Definition: solvers.hh:814
virtual void apply(X &x, X &b, InverseOperatorResult &res)
Apply inverse operator.
Definition: solvers.hh:859
X::field_type field_type
The field type of the operator to be inverted.
Definition: solvers.hh:821
MINRESSolver(L &op, S &sp, P &prec, real_type reduction, int maxit, int verbose)
Set up MINRES solver.
Definition: solvers.hh:845
MINRESSolver(L &op, P &prec, real_type reduction, int maxit, int verbose)
Set up MINRES solver.
Definition: solvers.hh:831
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: solvers.hh:823
X range_type
The range type of the operator to be inverted.
Definition: solvers.hh:819
virtual void apply(X &x, X &b, double reduction, InverseOperatorResult &res)
Apply inverse operator with given reduction factor.
Definition: solvers.hh:1044
X domain_type
The domain type of the operator to be inverted.
Definition: solvers.hh:817
virtual void post(X &x)=0
Clean up.
virtual void apply(X &v, const Y &d)=0
Apply one step of the preconditioner to the system A(v)=d.
virtual void pre(X &x, Y &b)=0
Prepare the preconditioner.
implements the Generalized Minimal Residual (GMRes) method
Definition: solvers.hh:1109
F basis_type
The field type of the basis vectors.
Definition: solvers.hh:1120
RestartedGMResSolver(L &op, P &prec, real_type reduction, int restart, int maxit, int verbose)
Set up solver.
Definition: solvers.hh:1148
Y range_type
The range type of the operator to be inverted.
Definition: solvers.hh:1114
virtual void apply(X &x, Y &b, double reduction, InverseOperatorResult &res)
Apply inverse operator.
Definition: solvers.hh:1205
X::field_type field_type
The field type of the operator to be inverted.
Definition: solvers.hh:1116
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: solvers.hh:1118
X domain_type
The domain type of the operator to be inverted.
Definition: solvers.hh:1112
virtual void apply(X &x, Y &b, InverseOperatorResult &res)
Apply inverse operator,.
Definition: solvers.hh:1195
RestartedGMResSolver(L &op, S &sp, P &prec, real_type reduction, int restart, int maxit, int verbose)
Set up solver.
Definition: solvers.hh:1183
Base class for scalar product and norm computation.
Definition: scalarproducts.hh:44
Default implementation for the scalar case.
Definition: scalarproducts.hh:96
A simple stop watch.
Definition: timer.hh:52
void reset()
Reset timer while keeping the running/stopped state.
Definition: timer.hh:66
double elapsed() const
Get elapsed user-time from last reset until now/last stop in seconds.
Definition: timer.hh:86
Definition of the DUNE_DEPRECATED macro for the case that config.h is not available.
Type traits to determine the type of reals (when working with complex numbers)
#define DUNE_DEPRECATED_MSG(text)
Mark some entity as deprecated.
Definition: deprecated.hh:169
#define DUNE_THROW(E, m)
Definition: exceptions.hh:243
Dune namespace.
Definition: alignment.hh:10
STL namespace.
Define general, extensible interface for operators. The available implementation wraps a matrix.
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:32
double elapsed
Elapsed time in seconds.
Definition: solver.hh:62
int iterations
Number of iterations.
Definition: solver.hh:50
double reduction
Reduction achieved: .
Definition: solver.hh:53
void clear()
Resets all data.
Definition: solver.hh:40
double conv_rate
Convergence rate (average reduction per step)
Definition: solver.hh:59
bool converged
True if convergence criterion has been met.
Definition: solver.hh:56
@ sequential
Category for sequential solvers.
Definition: solvercategory.hh:21
A simple timing class.
Traits for type conversions and type information.
Creative Commons License   |  Legal Statements / Impressum  |  Hosted by TU Dresden  |  generated with Hugo v0.111.3 (Nov 13, 23:29, 2024)