4 #ifndef DUNE_ISTL_SOLVERS_HH 5 #define DUNE_ISTL_SOLVERS_HH 21 #include <dune/common/deprecated.hh> 22 #include <dune/common/exceptions.hh> 23 #include <dune/common/timer.hh> 24 #include <dune/common/ftraits.hh> 25 #include <dune/common/typetraits.hh> 67 typedef typename FieldTraits<field_type>::real_type
real_type;
88 template<
class L,
class P>
90 real_type reduction,
int maxit,
int verbose) :
91 ssp(), _op(op), _prec(prec), _sp(ssp), _reduction(reduction), _maxit(maxit), _verbose(verbose)
93 static_assert(static_cast<int>(L::category) == static_cast<int>(P::category),
94 "L and P have to have the same category!");
96 "L has to be sequential!");
119 template<
class L,
class S,
class P>
121 real_type reduction,
int maxit,
int verbose) :
122 _op(op), _prec(prec), _sp(sp), _reduction(reduction), _maxit(maxit), _verbose(verbose)
124 static_assert(static_cast<int>(L::category) == static_cast<int>(P::category),
125 "L and P must have the same category!");
126 static_assert(static_cast<int>(L::category) == static_cast<int>(S::category),
127 "L and S must have the same category!");
147 real_type def0 = _sp.norm(b);
152 std::cout <<
"=== LoopSolver" << std::endl;
164 int i=1; real_type def=def0;
165 for ( ; i<=_maxit; i++ )
171 real_type defnew=_sp.norm(b);
176 if (def<def0*_reduction || def<1E-30)
184 i=std::min(_maxit,i);
202 std::cout <<
"=== rate=" << res.
conv_rate 205 <<
", IT=" << i << std::endl;
212 real_type saved_reduction = _reduction;
213 _reduction = reduction;
214 (*this).apply(x,b,res);
215 _reduction = saved_reduction;
223 real_type _reduction;
241 typedef typename FieldTraits<field_type>::real_type
real_type;
249 template<
class L,
class P>
251 real_type reduction,
int maxit,
int verbose) :
252 ssp(), _op(op), _prec(prec), _sp(ssp), _reduction(reduction), _maxit(maxit), _verbose(verbose)
254 static_assert(static_cast<int>(L::category) == static_cast<int>(P::category),
255 "L and P have to have the same category!");
257 "L has to be sequential!");
264 template<
class L,
class S,
class P>
266 real_type reduction,
int maxit,
int verbose) :
267 _op(op), _prec(prec), _sp(sp), _reduction(reduction), _maxit(maxit), _verbose(verbose)
269 static_assert(static_cast<int>(L::category) == static_cast<int>(P::category),
270 "L and P have to have the same category!");
271 static_assert(static_cast<int>(L::category) == static_cast<int>(S::category),
272 "L and S have to have the same category!");
285 _op.applyscaleadd(-1,x,b);
290 real_type def0 = _sp.norm(b);
294 std::cout <<
"=== GradientSolver" << std::endl;
302 int i=1; real_type def=def0;
304 for ( ; i<=_maxit; i++ )
309 lambda = _sp.dot(p,b)/_sp.dot(q,p);
313 real_type defnew=_sp.norm(b);
318 if (def<def0*_reduction || def<1E-30)
326 i=std::min(_maxit,i);
333 res.
reduction =
static_cast<double>(def/def0);
337 std::cout <<
"=== rate=" << res.
conv_rate 340 <<
", IT=" << i << std::endl;
350 real_type saved_reduction = _reduction;
351 _reduction = reduction;
352 (*this).apply(x,b,res);
353 _reduction = saved_reduction;
361 real_type _reduction;
379 typedef typename FieldTraits<field_type>::real_type
real_type;
386 template<
class L,
class P>
387 CGSolver (L& op, P& prec, real_type reduction,
int maxit,
int verbose) :
388 ssp(), _op(op), _prec(prec), _sp(ssp), _reduction(reduction), _maxit(maxit), _verbose(verbose)
390 static_assert(static_cast<int>(L::category) == static_cast<int>(P::category),
391 "L and P must have the same category!");
393 "L must be sequential!");
400 template<
class L,
class S,
class P>
401 CGSolver (L& op, S& sp, P& prec, real_type reduction,
int maxit,
int verbose) :
402 _op(op), _prec(prec), _sp(sp), _reduction(reduction), _maxit(maxit), _verbose(verbose)
404 static_assert(static_cast<int>(L::category) == static_cast<int>(P::category),
405 "L and P must have the same category!");
406 static_assert(static_cast<int>(L::category) == static_cast<int>(S::category),
407 "L and S must have the same category!");
428 _op.applyscaleadd(-1,x,b);
433 real_type def0 = _sp.norm(b);
438 std::cout <<
"=== CGSolver: abort due to infinite or NaN initial defect" 440 DUNE_THROW(
SolverAbort,
"CGSolver: initial defect=" << def0
441 <<
" is infinite or NaN");
452 std::cout <<
"=== rate=" << res.
conv_rate 454 <<
", IT=0" << std::endl;
460 std::cout <<
"=== CGSolver" << std::endl;
469 field_type rho,rholast,lambda,alpha,beta;
474 rholast = _sp.dot(p,b);
478 for ( ; i<=_maxit; i++ )
482 alpha = _sp.dot(p,q);
483 lambda = rholast/alpha;
488 real_type defnew=_sp.norm(b);
497 std::cout <<
"=== CGSolver: abort due to infinite or NaN defect" 500 "CGSolver: defect=" << def <<
" is infinite or NaN");
503 if (def<def0*_reduction || def<1E-30)
520 i=std::min(_maxit,i);
527 res.
reduction =
static_cast<double>(def/def0);
533 std::cout <<
"=== rate=" << res.
conv_rate 536 <<
", IT=" << i << std::endl;
551 virtual void apply (X& x, X& b,
double reduction,
554 real_type saved_reduction = _reduction;
555 _reduction = reduction;
556 (*this).apply(x,b,res);
557 _reduction = saved_reduction;
565 real_type _reduction;
583 typedef typename FieldTraits<field_type>::real_type
real_type;
590 template<
class L,
class P>
592 real_type reduction,
int maxit,
int verbose) :
593 ssp(), _op(op), _prec(prec), _sp(ssp), _reduction(reduction), _maxit(maxit), _verbose(verbose)
595 static_assert(static_cast<int>(L::category) == static_cast<int>(P::category),
596 "L and P must be of the same category!");
598 "L must be sequential!");
605 template<
class L,
class S,
class P>
607 real_type reduction,
int maxit,
int verbose) :
608 _op(op), _prec(prec), _sp(sp), _reduction(reduction), _maxit(maxit), _verbose(verbose)
610 static_assert(static_cast<int>(L::category) == static_cast<int>(P::category),
611 "L and P must have the same category!");
612 static_assert(static_cast<int>(L::category) == static_cast<int>(S::category),
613 "L and S must have the same category!");
626 const real_type EPSILON=1e-80;
628 field_type rho, rho_new, alpha, beta, h, omega;
629 real_type norm, norm_old, norm_0;
649 _op.applyscaleadd(-1,x,r);
653 norm = norm_old = norm_0 = _sp.norm(r);
664 std::cout <<
"=== BiCGSTABSolver" << std::endl;
674 if ( norm < (_reduction * norm_0) || norm<1E-30)
689 for (it = 0.5; it < _maxit; it+=.5)
696 rho_new = _sp.dot(rt,r);
699 if (abs(rho) <= EPSILON)
700 DUNE_THROW(
SolverAbort,
"breakdown in BiCGSTAB - rho " 701 << rho <<
" <= EPSILON " << EPSILON
702 <<
" after " << it <<
" iterations");
703 if (abs(omega) <= EPSILON)
704 DUNE_THROW(
SolverAbort,
"breakdown in BiCGSTAB - omega " 705 << omega <<
" <= EPSILON " << EPSILON
706 <<
" after " << it <<
" iterations");
713 beta = ( rho_new / rho ) * ( alpha / omega );
729 if (abs(h) < EPSILON)
730 DUNE_THROW(
SolverAbort,
"abs(h) < EPSILON in BiCGSTAB - abs(h) " 731 << abs(h) <<
" < EPSILON " << EPSILON
732 <<
" after " << it <<
" iterations");
754 if ( norm < (_reduction * norm_0) )
771 omega = _sp.dot(t,r)/_sp.dot(t,t);
793 if ( norm < (_reduction * norm_0) || norm<1E-30)
803 it=std::min(static_cast<double>(_maxit),it);
809 res.
iterations =
static_cast<int>(std::ceil(it));
810 res.
reduction =
static_cast<double>(norm/norm_0);
814 std::cout <<
"=== rate=" << res.
conv_rate 817 <<
", IT=" << it << std::endl;
829 real_type saved_reduction = _reduction;
830 _reduction = reduction;
831 (*this).apply(x,b,res);
832 _reduction = saved_reduction;
840 real_type _reduction;
861 typedef typename FieldTraits<field_type>::real_type
real_type;
868 template<
class L,
class P>
869 MINRESSolver (L& op, P& prec, real_type reduction,
int maxit,
int verbose) :
870 ssp(), _op(op), _prec(prec), _sp(ssp), _reduction(reduction), _maxit(maxit), _verbose(verbose)
872 static_assert(static_cast<int>(L::category) == static_cast<int>(P::category),
873 "L and P must have the same category!");
875 "L must be sequential!");
882 template<
class L,
class S,
class P>
883 MINRESSolver (L& op, S& sp, P& prec, real_type reduction,
int maxit,
int verbose) :
884 _op(op), _prec(prec), _sp(sp), _reduction(reduction), _maxit(maxit), _verbose(verbose)
886 static_assert(static_cast<int>(L::category) == static_cast<int>(P::category),
887 "L and P must have the same category!");
888 static_assert(static_cast<int>(L::category) == static_cast<int>(S::category),
889 "L and S must have the same category!");
909 _op.applyscaleadd(-1,x,b);
912 real_type def0 = _sp.norm(b);
916 std::cout <<
"=== MINRESSolver" << std::endl;
932 std::cout <<
"=== rate=" << res.
conv_rate 941 real_type def = def0;
943 field_type alpha, beta;
945 std::array<real_type,2> c{{0.0,0.0}};
947 std::array<field_type,2> s{{0.0,0.0}};
950 std::array<field_type,3> T{{0.0,0.0,0.0}};
953 std::array<field_type,2> xi{{1.0,0.0}};
964 beta = sqrt(_sp.dot(b,z));
965 field_type beta0 = beta;
968 std::array<X,3> p{{b,b,b}};
974 std::array<X,3> q{{b,b,b}};
983 for( ; i<=_maxit; i++) {
992 q[i2].axpy(-beta,q[i0]);
996 alpha = _sp.dot(z,q[i2]);
997 q[i2].axpy(-alpha,q[i1]);
1000 _prec.apply(z,q[i2]);
1004 beta = sqrt(_sp.dot(q[i2],z));
1017 T[2] = c[(i+1)%2]*alpha - s[(i+1)%2]*T[1];
1018 T[1] = c[(i+1)%2]*T[1] + s[(i+1)%2]*alpha;
1024 generateGivensRotation(T[2],beta,c[i%2],s[i%2]);
1026 T[2] = c[i%2]*T[2] + s[i%2]*beta;
1028 xi[i%2] = -s[i%2]*xi[(i+1)%2];
1029 xi[(i+1)%2] *= c[i%2];
1033 p[i2].axpy(-T[1],p[i1]);
1034 p[i2].axpy(-T[0],p[i0]);
1038 x.axpy(beta0*xi[(i+1)%2],p[i2]);
1045 real_type defnew = abs(beta0*xi[i%2]);
1051 if(def < def0*_reduction || def < 1e-30 || i == _maxit ) {
1064 res.
reduction =
static_cast<double>(def/def0);
1066 res.
elapsed = watch.elapsed();
1070 std::cout <<
"=== rate=" << res.
conv_rate 1073 <<
", IT=" << i << std::endl;
1084 real_type saved_reduction = _reduction;
1085 _reduction = reduction;
1086 (*this).apply(x,b,res);
1087 _reduction = saved_reduction;
1092 void generateGivensRotation(field_type &dx, field_type &dy, real_type &cs, field_type &sn)
1096 real_type norm_dx = abs(dx);
1097 real_type norm_dy = abs(dy);
1098 if(norm_dy < 1e-15) {
1101 }
else if(norm_dx < 1e-15) {
1104 }
else if(norm_dy > norm_dx) {
1105 real_type temp = norm_dx/norm_dy;
1106 cs = 1.0/sqrt(1.0 + temp*temp);
1114 real_type temp = norm_dy/norm_dx;
1115 cs = 1.0/sqrt(1.0 + temp*temp);
1127 real_type _reduction;
1145 template<
class X,
class Y=X,
class F = Y>
1156 typedef typename FieldTraits<field_type>::real_type
real_type;
1160 template<
class L,
class P>
1161 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")
1162 RestartedGMResSolver (L& op, P& prec, real_type reduction,
int restart,
int maxit,
int verbose,
bool recalc_defect)
1168 , _reduction(reduction)
1172 static_assert(static_cast<int>(P::category) == static_cast<int>(L::category),
1173 "P and L must be the same category!");
1175 "L must be sequential!");
1185 template<
class L,
class P>
1188 ssp(), _sp(ssp), _restart(restart),
1189 _reduction(reduction), _maxit(maxit), _verbose(verbose)
1191 static_assert(static_cast<int>(P::category) == static_cast<int>(L::category),
1192 "P and L must be the same category!");
1194 "L must be sequential!");
1197 template<
class L,
class S,
class P>
1198 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")
1199 RestartedGMResSolver(L& op, S& sp, P& prec, real_type reduction,
int restart,
int maxit,
int verbose,
bool recalc_defect)
1204 , _reduction(reduction)
1208 static_assert(static_cast<int>(P::category) == static_cast<int>(L::category),
1209 " P and L must have the same category!");
1210 static_assert(static_cast<int>(P::category) == static_cast<int>(S::category),
1211 "P and S must have the same category!");
1220 template<
class L,
class S,
class P>
1223 _sp(sp), _restart(restart),
1224 _reduction(reduction), _maxit(maxit), _verbose(verbose)
1226 static_assert(static_cast<int>(P::category) == static_cast<int>(L::category),
1227 "P and L must have the same category!");
1228 static_assert(static_cast<int>(P::category) == static_cast<int>(S::category),
1229 "P and S must have the same category!");
1242 apply(x,b,_reduction,res);
1257 const int m = _restart;
1260 std::vector<field_type> s(m+1), sn(m);
1261 std::vector<real_type> cs(m);
1266 std::vector< std::vector<field_type> > H(m+1,s);
1267 std::vector<F> v(m+1,b);
1278 _A.applyscaleadd(-1.0,x,b);
1280 v[0] = 0.0; _W.apply(v[0],b);
1281 norm_0 = _sp.norm(v[0]);
1288 std::cout <<
"=== RestartedGMResSolver" << std::endl;
1295 if(norm_0 < EPSILON) {
1302 while(j <= _maxit && res.
converged !=
true) {
1307 for(i=1; i<m+1; i++)
1310 for(i=0; i < m && j <= _maxit && res.
converged !=
true; i++, j++) {
1315 _A.apply(v[i],v[i+1]);
1317 for(
int k=0; k<i+1; k++) {
1322 H[k][i] = _sp.dot(v[k],w);
1324 w.axpy(-H[k][i],v[k]);
1326 H[i+1][i] = _sp.norm(w);
1327 if(abs(H[i+1][i]) < EPSILON)
1329 "breakdown in GMRes - |w| == 0.0 after " << j <<
" iterations");
1332 v[i+1] = w; v[i+1] *= 1.0/H[i+1][i];
1335 for(
int k=0; k<i; k++)
1336 applyPlaneRotation(H[k][i],H[k+1][i],cs[k],sn[k]);
1339 generatePlaneRotation(H[i][i],H[i+1][i],cs[i],sn[i]);
1341 applyPlaneRotation(H[i][i],H[i+1][i],cs[i],sn[i]);
1342 applyPlaneRotation(s[i],s[i+1],cs[i],sn[i]);
1355 if(norm < reduction * norm_0)
1369 if( res.
converged !=
true && j <= _maxit ) {
1372 std::cout <<
"=== GMRes::restart" << std::endl;
1376 _A.applyscaleadd(-1.0,x,b);
1380 norm = _sp.norm(v[0]);
1391 res.
reduction =
static_cast<double>(norm/norm_0);
1393 res.
elapsed = watch.elapsed();
1404 std::cout <<
"=== rate=" << res.
conv_rate 1411 void update(X& w,
int i,
1412 const std::vector<std::vector<field_type> >& H,
1413 const std::vector<field_type>& s,
1414 const std::vector<X>& v) {
1416 std::vector<field_type> y(s);
1419 for(
int a=i-1; a>=0; a--) {
1421 for(
int b=a+1; b<i; b++)
1422 rhs -= H[a][b]*y[b];
1431 template<
typename T>
1432 typename std::enable_if<std::is_same<field_type,real_type>::value,T>::type conjugate(
const T& t) {
1436 template<
typename T>
1437 typename std::enable_if<!std::is_same<field_type,real_type>::value,T>::type conjugate(
const T& t) {
1448 if(norm_dy < 1e-15) {
1451 }
else if(norm_dx < 1e-15) {
1454 }
else if(norm_dy > norm_dx) {
1456 cs = 1.0/sqrt(1.0 + temp*temp);
1460 sn *= conjugate(dy)/norm_dy;
1463 cs = 1.0/sqrt(1.0 + temp*temp);
1465 sn *= conjugate(dy/dx);
1474 dy = -conjugate(sn) * dx + cs * dy;
1513 typedef typename FieldTraits<field_type>::real_type
real_type;
1522 template<
class L,
class P>
1525 ssp(), _op(op), _prec(prec), _sp(ssp), _reduction(reduction), _maxit(maxit),
1526 _verbose(verbose), _restart(
std::min(maxit,restart))
1528 static_assert(static_cast<int>(L::category) == static_cast<int>(P::category),
1529 "L and P have to have the same category!");
1530 static_assert(static_cast<int>(L::category) ==
1532 "L has to be sequential!");
1541 template<
class L,
class P,
class S>
1543 real_type reduction,
int maxit,
int verbose,
int restart=10) :
1544 _op(op), _prec(prec), _sp(sp), _reduction(reduction), _maxit(maxit), _verbose(verbose),
1545 _restart(
std::min(maxit,restart))
1547 static_assert(static_cast<int>(L::category) == static_cast<int>(P::category),
1548 "L and P must have the same category!");
1549 static_assert(static_cast<int>(L::category) == static_cast<int>(S::category),
1550 "L and S must have the same category!");
1562 _op.applyscaleadd(-1,x,b);
1564 std::vector<std::shared_ptr<X> > p(_restart);
1565 std::vector<typename X::field_type> pp(_restart);
1569 p[0].reset(
new X(x));
1571 real_type def0 = _sp.norm(b);
1580 std::cout <<
"=== rate=" << res.
conv_rate 1582 <<
", IT=0" << std::endl;
1588 std::cout <<
"=== GeneralizedPCGSolver" << std::endl;
1596 field_type rho, lambda;
1602 _prec.apply(*(p[0]),b);
1603 rho = _sp.dot(*(p[0]),b);
1604 _op.apply(*(p[0]),q);
1605 pp[0] = _sp.dot(*(p[0]),q);
1607 x.axpy(lambda,*(p[0]));
1611 real_type defnew=_sp.norm(b);
1615 if (def<def0*_reduction || def<1E-30)
1620 std::cout <<
"=== rate=" << res.
conv_rate 1623 <<
", IT=" << 1 << std::endl;
1630 int end=std::min(_restart, _maxit-i+1);
1631 for (ii=1; ii<end; ++ii )
1636 _prec.apply(prec_res,b);
1638 p[ii].reset(
new X(prec_res));
1639 _op.apply(prec_res, q);
1641 for(
int j=0; j<ii; ++j) {
1642 rho =_sp.dot(q,*(p[j]))/pp[j];
1643 p[ii]->axpy(-rho, *(p[j]));
1647 _op.apply(*(p[ii]),q);
1648 pp[ii] = _sp.dot(*(p[ii]),q);
1649 rho = _sp.dot(*(p[ii]),b);
1650 lambda = rho/pp[ii];
1651 x.axpy(lambda,*(p[ii]));
1655 real_type defNew=_sp.norm(b);
1661 if (def<def0*_reduction || def<1E-30)
1670 *(p[0])=*(p[_restart-1]);
1671 pp[0]=pp[_restart-1];
1682 res.
elapsed = watch.elapsed();
1686 std::cout <<
"=== rate=" << res.
conv_rate 1689 <<
", IT=" << i+1 << std::endl;
1698 virtual void apply (X& x, X& b,
double reduction,
1701 real_type saved_reduction = _reduction;
1702 _reduction = reduction;
1703 (*this).apply(x,b,res);
1704 _reduction = saved_reduction;
1711 real_type _reduction;
CGSolver(L &op, S &sp, P &prec, real_type reduction, int maxit, int verbose)
Set up conjugate gradient solver.
Definition: solvers.hh:401
X::field_type field_type
The field type of the operator to be inverted.
Definition: solvers.hh:581
implements the Generalized Minimal Residual (GMRes) method
Definition: solvers.hh:1146
virtual void apply(X &x, X &b, InverseOperatorResult &res)
Apply inverse operator.
Definition: solvers.hh:1557
X range_type
The range type of the operator to be inverted.
Definition: solvers.hh:579
Default implementation for the scalar case.
Definition: scalarproducts.hh:95
Thrown when a solver aborts due to some problem.
Definition: istlexception.hh:43
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:583
virtual void apply(X &x, X &b, double reduction, InverseOperatorResult &res)
Apply inverse operator with given reduction factor.
Definition: solvers.hh:551
conjugate gradient method
Definition: solvers.hh:370
Y range_type
The range type of the operator to be inverted.
Definition: solvers.hh:1152
X domain_type
The domain type of the operator to be inverted.
Definition: solvers.hh:855
virtual void apply(X &x, X &b, double reduction, InverseOperatorResult &res)
Apply inverse operator with given reduction factor.
Definition: solvers.hh:1698
virtual void apply(X &x, X &b, InverseOperatorResult &res)
Apply inverse operator.
Definition: solvers.hh:623
double reduction
Reduction achieved: .
Definition: solver.hh:53
X::field_type field_type
The field type of the operator that we do the inverse for.
Definition: solvers.hh:65
virtual void apply(X &x, X &b, double reduction, InverseOperatorResult &res)
apply inverse operator, with given convergence criteria.
Definition: solvers.hh:210
X::field_type field_type
The field type of the operator to be inverted.
Definition: solvers.hh:1511
Definition: example.cc:35
X domain_type
The domain type of the operator to be inverted.
Definition: solvers.hh:373
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:379
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:861
X domain_type
The domain type of the operator that we do the inverse for.
Definition: solvers.hh:61
RestartedGMResSolver(L &op, S &sp, P &prec, real_type reduction, int restart, int maxit, int verbose)
Set up solver.
Definition: solvers.hh:1221
double elapsed
Elapsed time in seconds.
Definition: solver.hh:62
virtual void pre(X &x, Y &b)=0
Prepare the preconditioner.
F basis_type
The field type of the basis vectors.
Definition: solvers.hh:1158
X::field_type field_type
The field type of the operator to be inverted.
Definition: solvers.hh:859
virtual void apply(X &x, Y &b, double reduction, InverseOperatorResult &res)
Apply inverse operator.
Definition: solvers.hh:1253
virtual void apply(X &v, const Y &d)=0
Apply one step of the preconditioner to the system A(v)=d.
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:1513
virtual void apply(X &x, X &b, InverseOperatorResult &res)
Apply inverse operator.
Definition: solvers.hh:280
MINRESSolver(L &op, P &prec, real_type reduction, int maxit, int verbose)
Set up MINRES solver.
Definition: solvers.hh:869
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:67
BiCGSTABSolver(L &op, S &sp, P &prec, real_type reduction, int maxit, int verbose)
Set up solver.
Definition: solvers.hh:606
X::field_type field_type
The field type of the operator to be inverted.
Definition: solvers.hh:1154
Category for sequential solvers.
Definition: solvercategory.hh:21
BiCGSTABSolver(L &op, P &prec, real_type reduction, int maxit, int verbose)
Set up solver.
Definition: solvers.hh:591
Minimal Residual Method (MINRES)
Definition: solvers.hh:852
virtual void apply(X &x, X &b, InverseOperatorResult &res)
Apply inverse operator.
Definition: solvers.hh:421
gradient method
Definition: solvers.hh:232
Define general, extensible interface for inverse operators.
Abstract base class for all solvers.
Definition: solver.hh:79
X domain_type
The domain type of the operator to be inverted.
Definition: solvers.hh:577
CGSolver(L &op, P &prec, real_type reduction, int maxit, int verbose)
Set up conjugate gradient solver.
Definition: solvers.hh:387
double conv_rate
Convergence rate (average reduction per step)
Definition: solver.hh:59
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:136
virtual void applyscaleadd(field_type alpha, const X &x, Y &y) const =0
apply operator to x, scale and add:
void clear()
Resets all data.
Definition: solver.hh:40
X::field_type field_type
The field type of the operator to be inverted.
Definition: solvers.hh:377
Generalized preconditioned conjugate gradient solver.
Definition: solvers.hh:1503
Define base class for scalar product and norm.
LoopSolver(L &op, S &sp, P &prec, real_type reduction, int maxit, int verbose)
Set up loop solver.
Definition: solvers.hh:120
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:241
Bi-conjugate Gradient Stabilized (BiCG-STAB)
Definition: solvers.hh:574
virtual void apply(X &x, Y &b, InverseOperatorResult &res)
Apply inverse operator.
Definition: solvers.hh:1240
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:1156
Definition: basearray.hh:19
bool converged
True if convergence criterion has been met.
Definition: solver.hh:56
GradientSolver(L &op, P &prec, real_type reduction, int maxit, int verbose)
Set up solver.
Definition: solvers.hh:250
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:1523
X domain_type
The domain type of the operator to be inverted.
Definition: solvers.hh:1507
virtual void apply(X &x, X &b, double reduction, InverseOperatorResult &res)
Apply inverse operator with given reduction factor.
Definition: solvers.hh:827
X::field_type field_type
The field type of the operator that we do the inverse for.
Definition: solvers.hh:239
virtual void apply(X &x, X &b, InverseOperatorResult &res)
Apply inverse operator.
Definition: solvers.hh:897
X domain_type
The domain type of the operator that we do the inverse for.
Definition: solvers.hh:235
GradientSolver(L &op, S &sp, P &prec, real_type reduction, int maxit, int verbose)
Set up solver.
Definition: solvers.hh:265
X domain_type
The domain type of the operator to be inverted.
Definition: solvers.hh:1150
virtual void apply(X &x, X &b, double reduction, InverseOperatorResult &res)
Apply inverse operator with given reduction factor.
Definition: solvers.hh:348
LoopSolver(L &op, P &prec, real_type reduction, int maxit, int verbose)
Set up Loop solver.
Definition: solvers.hh:89
void printHeader(std::ostream &s) const
helper function for printing header of solver output
Definition: solver.hh:127
RestartedGMResSolver(L &op, P &prec, real_type reduction, int restart, int maxit, int verbose)
Set up solver.
Definition: solvers.hh:1186
int iterations
Number of iterations.
Definition: solver.hh:50
Base class for scalar product and norm computation.
Definition: scalarproducts.hh:44
X range_type
The range type of the operator to be inverted.
Definition: solvers.hh:857
X::field_type field_type
The field type of the operator.
Definition: solver.hh:88
X range_type
The range type of the operator to be inverted.
Definition: solvers.hh:375
Define general, extensible interface for operators. The available implementation wraps a matrix...
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:1542
X range_type
The range type of the operator that we do the inverse for.
Definition: solvers.hh:63
MINRESSolver(L &op, S &sp, P &prec, real_type reduction, int maxit, int verbose)
Set up MINRES solver.
Definition: solvers.hh:883
X range_type
The range type of the operator to be inverted.
Definition: solvers.hh:1509
Preconditioned loop solver.
Definition: solvers.hh:58
X range_type
The range type of the operator that we do the inverse for.
Definition: solvers.hh:237
virtual void post(X &x)=0
Clean up.
Statistics about the application of an inverse operator.
Definition: solver.hh:31
virtual void apply(X &x, X &b, double reduction, InverseOperatorResult &res)
Apply inverse operator with given reduction factor.
Definition: solvers.hh:1082
virtual void apply(X &x, X &b, InverseOperatorResult &res)
Apply inverse operator,.
Definition: solvers.hh:132