4 #ifndef DUNE_ISTL_SOLVERS_HH
5 #define DUNE_ISTL_SOLVERS_HH
20 #include <dune/common/array.hh>
21 #include <dune/common/deprecated.hh>
22 #include <dune/common/timer.hh>
23 #include <dune/common/ftraits.hh>
24 #include <dune/common/typetraits.hh>
66 typedef typename FieldTraits<field_type>::real_type
real_type;
87 template<
class L,
class P>
89 real_type reduction,
int maxit,
int verbose) :
90 ssp(), _op(op), _prec(prec), _sp(ssp), _reduction(reduction), _maxit(maxit), _verbose(verbose)
92 static_assert(static_cast<int>(L::category) == static_cast<int>(P::category),
93 "L and P have to have the same category!");
95 "L has to be sequential!");
118 template<
class L,
class S,
class P>
120 real_type reduction,
int maxit,
int verbose) :
121 _op(op), _prec(prec), _sp(sp), _reduction(reduction), _maxit(maxit), _verbose(verbose)
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!");
146 real_type def0 = _sp.norm(b);
151 std::cout <<
"=== LoopSolver" << std::endl;
163 int i=1; real_type def=def0;
164 for ( ; i<=_maxit; i++ )
170 real_type defnew=_sp.norm(b);
175 if (def<def0*_reduction || def<1E-30)
183 i=std::min(_maxit,i);
201 std::cout <<
"=== rate=" << res.
conv_rate
204 <<
", IT=" << i << std::endl;
211 real_type saved_reduction = _reduction;
212 _reduction = reduction;
213 (*this).apply(x,b,res);
214 _reduction = saved_reduction;
222 real_type _reduction;
240 typedef typename FieldTraits<field_type>::real_type
real_type;
248 template<
class L,
class P>
250 real_type reduction,
int maxit,
int verbose) :
251 ssp(), _op(op), _prec(prec), _sp(ssp), _reduction(reduction), _maxit(maxit), _verbose(verbose)
253 static_assert(static_cast<int>(L::category) == static_cast<int>(P::category),
254 "L and P have to have the same category!");
256 "L has to be sequential!");
263 template<
class L,
class S,
class P>
265 real_type reduction,
int maxit,
int verbose) :
266 _op(op), _prec(prec), _sp(sp), _reduction(reduction), _maxit(maxit), _verbose(verbose)
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!");
289 real_type def0 = _sp.norm(b);
293 std::cout <<
"=== GradientSolver" << std::endl;
301 int i=1; real_type def=def0;
303 for ( ; i<=_maxit; i++ )
308 lambda = _sp.dot(p,b)/_sp.dot(q,p);
312 real_type defnew=_sp.norm(b);
317 if (def<def0*_reduction || def<1E-30)
325 i=std::min(_maxit,i);
332 res.
reduction =
static_cast<double>(def/def0);
336 std::cout <<
"=== rate=" << res.
conv_rate
339 <<
", IT=" << i << std::endl;
349 real_type saved_reduction = _reduction;
350 _reduction = reduction;
351 (*this).apply(x,b,res);
352 _reduction = saved_reduction;
360 real_type _reduction;
378 typedef typename FieldTraits<field_type>::real_type
real_type;
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)
389 static_assert(static_cast<int>(L::category) == static_cast<int>(P::category),
390 "L and P must have the same category!");
392 "L must be sequential!");
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)
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!");
424 real_type def0 = _sp.norm(b);
433 std::cout <<
"=== rate=" << res.
conv_rate
435 <<
", IT=0" << std::endl;
441 std::cout <<
"=== CGSolver" << std::endl;
450 field_type rho,rholast,lambda,alpha,beta;
455 rholast = _sp.dot(p,b);
459 for ( ; i<=_maxit; i++ )
463 alpha = _sp.dot(p,q);
464 lambda = rholast/alpha;
469 real_type defnew=_sp.norm(b);
475 if (def<def0*_reduction || def<1E-30)
492 i=std::min(_maxit,i);
499 res.
reduction =
static_cast<double>(def/def0);
505 std::cout <<
"=== rate=" << res.
conv_rate
508 <<
", IT=" << i << std::endl;
517 virtual void apply (X& x, X& b,
double reduction,
520 real_type saved_reduction = _reduction;
521 _reduction = reduction;
522 (*this).apply(x,b,res);
523 _reduction = saved_reduction;
531 real_type _reduction;
549 typedef typename FieldTraits<field_type>::real_type
real_type;
556 template<
class L,
class P>
558 real_type reduction,
int maxit,
int verbose) :
559 ssp(), _op(op), _prec(prec), _sp(ssp), _reduction(reduction), _maxit(maxit), _verbose(verbose)
561 static_assert(static_cast<int>(L::category) == static_cast<int>(P::category),
562 "L and P must be of the same category!");
564 "L must be sequential!");
571 template<
class L,
class S,
class P>
573 real_type reduction,
int maxit,
int verbose) :
574 _op(op), _prec(prec), _sp(sp), _reduction(reduction), _maxit(maxit), _verbose(verbose)
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!");
589 const real_type EPSILON=1e-80;
591 field_type rho, rho_new, alpha, beta, h, omega;
592 real_type norm, norm_old, norm_0;
616 norm = norm_old = norm_0 = _sp.norm(r);
627 std::cout <<
"=== BiCGSTABSolver" << std::endl;
637 if ( norm < (_reduction * norm_0) || norm<1E-30)
652 for (it = 0.5; it < _maxit; it+=.5)
659 rho_new = _sp.dot(rt,r);
662 if (std::abs(rho) <= EPSILON)
663 DUNE_THROW(
ISTLError,
"breakdown in BiCGSTAB - rho "
664 << rho <<
" <= EPSILON " << EPSILON
665 <<
" after " << it <<
" iterations");
666 if (std::abs(omega) <= EPSILON)
667 DUNE_THROW(
ISTLError,
"breakdown in BiCGSTAB - omega "
668 << omega <<
" <= EPSILON " << EPSILON
669 <<
" after " << it <<
" iterations");
676 beta = ( rho_new / rho ) * ( alpha / omega );
692 if ( std::abs(h) < EPSILON )
715 if ( norm < (_reduction * norm_0) )
732 omega = _sp.dot(t,r)/_sp.dot(t,t);
754 if ( norm < (_reduction * norm_0) || norm<1E-30)
764 it=std::min(static_cast<double>(_maxit),it);
770 res.
iterations =
static_cast<int>(std::ceil(it));
771 res.
reduction =
static_cast<double>(norm/norm_0);
775 std::cout <<
"=== rate=" << res.
conv_rate
778 <<
", IT=" << it << std::endl;
788 real_type saved_reduction = _reduction;
789 _reduction = reduction;
790 (*this).apply(x,b,res);
791 _reduction = saved_reduction;
799 real_type _reduction;
820 typedef typename FieldTraits<field_type>::real_type
real_type;
827 template<
class L,
class P>
828 MINRESSolver (L& op, P& prec, real_type reduction,
int maxit,
int verbose) :
829 ssp(), _op(op), _prec(prec), _sp(ssp), _reduction(reduction), _maxit(maxit), _verbose(verbose)
831 static_assert(static_cast<int>(L::category) == static_cast<int>(P::category),
832 "L and P must have the same category!");
834 "L must be sequential!");
841 template<
class L,
class S,
class P>
842 MINRESSolver (L& op, S& sp, P& prec, real_type reduction,
int maxit,
int verbose) :
843 _op(op), _prec(prec), _sp(sp), _reduction(reduction), _maxit(maxit), _verbose(verbose)
845 static_assert(static_cast<int>(L::category) == static_cast<int>(P::category),
846 "L and P must have the same category!");
847 static_assert(static_cast<int>(L::category) == static_cast<int>(S::category),
848 "L and S must have the same category!");
869 real_type def0 = _sp.norm(b);
873 std::cout <<
"=== MINRESSolver" << std::endl;
889 std::cout <<
"=== rate=" << res.
conv_rate
898 real_type def = def0;
900 field_type alpha, beta;
902 Dune::array<real_type,2> c{{0.0,0.0}};
904 Dune::array<field_type,2> s{{0.0,0.0}};
907 Dune::array<field_type,3> T{{0.0,0.0,0.0}};
910 Dune::array<field_type,2> xi{{1.0,0.0}};
921 beta = std::sqrt(_sp.dot(b,z));
922 field_type beta0 = beta;
925 Dune::array<X,3> p{{b,b,b}};
931 Dune::array<X,3> q{{b,b,b}};
940 for( ; i<=_maxit; i++) {
949 q[i2].axpy(-beta,q[i0]);
953 alpha = _sp.dot(z,q[i2]);
954 q[i2].axpy(-alpha,q[i1]);
957 _prec.
apply(z,q[i2]);
961 beta = std::sqrt(_sp.dot(q[i2],z));
974 T[2] = c[(i+1)%2]*alpha - s[(i+1)%2]*T[1];
975 T[1] = c[(i+1)%2]*T[1] + s[(i+1)%2]*alpha;
981 generateGivensRotation(T[2],beta,c[i%2],s[i%2]);
983 T[2] = c[i%2]*T[2] + s[i%2]*beta;
985 xi[i%2] = -s[i%2]*xi[(i+1)%2];
986 xi[(i+1)%2] *= c[i%2];
990 p[i2].axpy(-T[1],p[i1]);
991 p[i2].axpy(-T[0],p[i0]);
995 x.axpy(beta0*xi[(i+1)%2],p[i2]);
1002 real_type defnew = std::abs(beta0*xi[i%2]);
1008 if(def < def0*_reduction || def < 1e-30 || i == _maxit ) {
1021 res.
reduction =
static_cast<double>(def/def0);
1023 res.
elapsed = watch.elapsed();
1027 std::cout <<
"=== rate=" << res.
conv_rate
1030 <<
", IT=" << i << std::endl;
1041 real_type saved_reduction = _reduction;
1042 _reduction = reduction;
1043 (*this).apply(x,b,res);
1044 _reduction = saved_reduction;
1049 void generateGivensRotation(field_type &dx, field_type &dy, real_type &cs, field_type &sn)
1051 real_type norm_dx = std::abs(dx);
1052 real_type norm_dy = std::abs(dy);
1053 if(norm_dy < 1e-15) {
1056 }
else if(norm_dx < 1e-15) {
1059 }
else if(norm_dy > norm_dx) {
1060 real_type temp = norm_dx/norm_dy;
1061 cs = 1.0/std::sqrt(1.0 + temp*temp);
1069 real_type temp = norm_dy/norm_dx;
1070 cs = 1.0/std::sqrt(1.0 + temp*temp);
1078 SeqScalarProduct<X> ssp;
1079 LinearOperator<X,X>& _op;
1080 Preconditioner<X,X>& _prec;
1081 ScalarProduct<X>& _sp;
1082 real_type _reduction;
1100 template<
class X,
class Y=X,
class F = Y>
1111 typedef typename FieldTraits<field_type>::real_type
real_type;
1115 template<
class L,
class P>
1116 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")
1117 RestartedGMResSolver (L& op, P& prec, real_type reduction,
int restart,
int maxit,
int verbose,
bool recalc_defect)
1123 , _reduction(reduction)
1127 static_assert(static_cast<int>(P::category) == static_cast<int>(L::category),
1128 "P and L must be the same category!");
1130 "L must be sequential!");
1140 template<
class L,
class P>
1143 ssp(), _sp(ssp), _restart(restart),
1144 _reduction(reduction), _maxit(maxit), _verbose(verbose)
1146 static_assert(static_cast<int>(P::category) == static_cast<int>(L::category),
1147 "P and L must be the same category!");
1149 "L must be sequential!");
1152 template<
class L,
class S,
class P>
1153 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")
1154 RestartedGMResSolver(L& op, S& sp, P& prec, real_type reduction,
int restart,
int maxit,
int verbose,
bool recalc_defect)
1159 , _reduction(reduction)
1163 static_assert(static_cast<int>(P::category) == static_cast<int>(L::category),
1164 " P and L must have the same category!");
1165 static_assert(static_cast<int>(P::category) == static_cast<int>(S::category),
1166 "P and S must have the same category!");
1175 template<
class L,
class S,
class P>
1178 _sp(sp), _restart(restart),
1179 _reduction(reduction), _maxit(maxit), _verbose(verbose)
1181 static_assert(static_cast<int>(P::category) == static_cast<int>(L::category),
1182 "P and L must have the same category!");
1183 static_assert(static_cast<int>(P::category) == static_cast<int>(S::category),
1184 "P and S must have the same category!");
1190 apply(x,b,_reduction,res);
1200 const real_type EPSILON = 1e-80;
1201 const int m = _restart;
1202 real_type norm, norm_old = 0.0, norm_0;
1204 std::vector<field_type> s(m+1), sn(m);
1205 std::vector<real_type> cs(m);
1210 std::vector< std::vector<field_type> > H(m+1,s);
1211 std::vector<F> v(m+1,b);
1222 _A.applyscaleadd(-1.0,x,b);
1224 v[0] = 0.0; _W.apply(v[0],b);
1225 norm_0 = _sp.norm(v[0]);
1232 std::cout <<
"=== RestartedGMResSolver" << std::endl;
1239 if(norm_0 < EPSILON) {
1246 while(j <= _maxit && res.
converged !=
true) {
1251 for(i=1; i<m+1; i++)
1254 for(i=0; i < m && j <= _maxit && res.
converged !=
true; i++, j++) {
1259 _A.apply(v[i],v[i+1]);
1261 for(
int k=0; k<i+1; k++) {
1266 H[k][i] = _sp.dot(v[k],w);
1268 w.axpy(-H[k][i],v[k]);
1270 H[i+1][i] = _sp.norm(w);
1271 if(std::abs(H[i+1][i]) < EPSILON)
1273 "breakdown in GMRes - |w| == 0.0 after " << j <<
" iterations");
1276 v[i+1] = w; v[i+1] *= 1.0/H[i+1][i];
1279 for(
int k=0; k<i; k++)
1280 applyPlaneRotation(H[k][i],H[k+1][i],cs[k],sn[k]);
1283 generatePlaneRotation(H[i][i],H[i+1][i],cs[i],sn[i]);
1285 applyPlaneRotation(H[i][i],H[i+1][i],cs[i],sn[i]);
1286 applyPlaneRotation(s[i],s[i+1],cs[i],sn[i]);
1289 norm = std::abs(s[i+1]);
1299 if(norm < reduction * norm_0)
1313 if( res.
converged !=
true && j <= _maxit ) {
1316 std::cout <<
"=== GMRes::restart" << std::endl;
1320 _A.applyscaleadd(-1.0,x,b);
1324 norm = _sp.norm(v[0]);
1335 res.
reduction =
static_cast<double>(norm/norm_0);
1337 res.
elapsed = watch.elapsed();
1348 std::cout <<
"=== rate=" << res.
conv_rate
1355 void update(X& w,
int i,
1356 const std::vector<std::vector<field_type> >& H,
1357 const std::vector<field_type>& s,
1358 const std::vector<X>& v) {
1360 std::vector<field_type> y(s);
1363 for(
int a=i-1; a>=0; a--) {
1364 field_type rhs(s[a]);
1365 for(
int b=a+1; b<i; b++)
1366 rhs -= H[a][b]*y[b];
1375 template<
typename T>
1376 typename enable_if<is_same<field_type,real_type>::value,T>::type conjugate(
const T& t) {
1380 template<
typename T>
1381 typename enable_if<!is_same<field_type,real_type>::value,T>::type conjugate(
const T& t) {
1386 generatePlaneRotation(field_type &dx, field_type &dy, real_type &cs, field_type &sn)
1388 real_type norm_dx = std::abs(dx);
1389 real_type norm_dy = std::abs(dy);
1390 if(norm_dy < 1e-15) {
1393 }
else if(norm_dx < 1e-15) {
1396 }
else if(norm_dy > norm_dx) {
1397 real_type temp = norm_dx/norm_dy;
1398 cs = 1.0/std::sqrt(1.0 + temp*temp);
1402 sn *= conjugate(dy)/norm_dy;
1404 real_type temp = norm_dy/norm_dx;
1405 cs = 1.0/std::sqrt(1.0 + temp*temp);
1407 sn *= conjugate(dy/dx);
1413 applyPlaneRotation(field_type &dx, field_type &dy, real_type &cs, field_type &sn)
1415 field_type temp = cs * dx + sn * dy;
1416 dy = -conjugate(sn) * dx + cs * dy;
1420 LinearOperator<X,Y>& _A;
1421 Preconditioner<X,Y>& _W;
1422 SeqScalarProduct<X> ssp;
1423 ScalarProduct<X>& _sp;
1425 real_type _reduction;
1455 typedef typename FieldTraits<field_type>::real_type
real_type;
1464 template<
class L,
class P>
1467 ssp(), _op(op), _prec(prec), _sp(ssp), _reduction(reduction), _maxit(maxit),
1468 _verbose(verbose), _restart(
std::min(maxit,restart))
1470 static_assert(static_cast<int>(L::category) == static_cast<int>(P::category),
1471 "L and P have to have the same category!");
1472 static_assert(static_cast<int>(L::category) ==
1474 "L has to be sequential!");
1483 template<
class L,
class P,
class S>
1485 real_type reduction,
int maxit,
int verbose,
int restart=10) :
1486 _op(op), _prec(prec), _sp(sp), _reduction(reduction), _maxit(maxit), _verbose(verbose),
1487 _restart(
std::min(maxit,restart))
1489 static_assert(static_cast<int>(L::category) == static_cast<int>(P::category),
1490 "L and P must have the same category!");
1491 static_assert(static_cast<int>(L::category) == static_cast<int>(S::category),
1492 "L and S must have the same category!");
1506 std::vector<std::shared_ptr<X> > p(_restart);
1507 std::vector<typename X::field_type> pp(_restart);
1511 p[0].reset(
new X(x));
1513 real_type def0 = _sp.norm(b);
1522 std::cout <<
"=== rate=" << res.
conv_rate
1524 <<
", IT=0" << std::endl;
1530 std::cout <<
"=== GeneralizedPCGSolver" << std::endl;
1538 field_type rho, lambda;
1544 _prec.
apply(*(p[0]),b);
1545 rho = _sp.dot(*(p[0]),b);
1546 _op.
apply(*(p[0]),q);
1547 pp[0] = _sp.dot(*(p[0]),q);
1549 x.axpy(lambda,*(p[0]));
1553 real_type defnew=_sp.norm(b);
1557 if (def<def0*_reduction || def<1E-30)
1562 std::cout <<
"=== rate=" << res.
conv_rate
1565 <<
", IT=" << 1 << std::endl;
1572 int end=std::min(_restart, _maxit-i+1);
1573 for (ii=1; ii<end; ++ii )
1578 _prec.
apply(prec_res,b);
1580 p[ii].reset(
new X(prec_res));
1581 _op.
apply(prec_res, q);
1583 for(
int j=0; j<ii; ++j) {
1584 rho =_sp.dot(q,*(p[j]))/pp[j];
1585 p[ii]->axpy(-rho, *(p[j]));
1589 _op.
apply(*(p[ii]),q);
1590 pp[ii] = _sp.dot(*(p[ii]),q);
1591 rho = _sp.dot(*(p[ii]),b);
1592 lambda = rho/pp[ii];
1593 x.axpy(lambda,*(p[ii]));
1597 real_type defnew=_sp.norm(b);
1603 if (def<def0*_reduction || def<1E-30)
1612 *(p[0])=*(p[_restart-1]);
1613 pp[0]=pp[_restart-1];
1624 res.
elapsed = watch.elapsed();
1628 std::cout <<
"=== rate=" << res.
conv_rate
1631 <<
", IT=" << i+1 << std::endl;
1640 virtual void apply (X& x, X& b,
double reduction,
1643 real_type saved_reduction = _reduction;
1644 _reduction = reduction;
1645 (*this).apply(x,b,res);
1646 _reduction = saved_reduction;
1653 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:400
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
Y range_type
The range type of the operator to be inverted.
Definition: solvers.hh:1107
virtual void apply(X &x, Y &b, InverseOperatorResult &res)
Apply inverse operator,.
Definition: solvers.hh:1188
Base class for scalar product and norm computation.
Definition: scalarproducts.hh:43
Preconditioned loop solver.
Definition: solvers.hh:57
X::field_type field_type
The field type of the operator to be inverted.
Definition: solvers.hh:1109
Definition: basearray.hh:19
MINRESSolver(L &op, S &sp, P &prec, real_type reduction, int maxit, int verbose)
Set up MINRES solver.
Definition: solvers.hh:842
RestartedGMResSolver(L &op, P &prec, real_type reduction, int restart, int maxit, int verbose)
Set up solver.
Definition: solvers.hh:1141
double conv_rate
Convergence rate (average reduction per step)
Definition: solver.hh:59
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
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:1484
Define base class for scalar product and norm.
X range_type
The range type of the operator to be inverted.
Definition: solvers.hh:1451
derive error class from the base class in common
Definition: istlexception.hh:16
virtual void apply(X &x, X &b, InverseOperatorResult &res)
Apply inverse operator.
Definition: solvers.hh:1499
RestartedGMResSolver(L &op, S &sp, P &prec, real_type reduction, int restart, int maxit, int verbose)
Set up solver.
Definition: solvers.hh:1176
implements the Generalized Minimal Residual (GMRes) method
Definition: solvers.hh:1101
gradient method
Definition: solvers.hh:231
F basis_type
The field type of the basis vectors.
Definition: solvers.hh:1113
virtual void apply(X &x, X &b, double reduction, InverseOperatorResult &res)
Apply inverse operator with given reduction factor.
Definition: solvers.hh:1039
LoopSolver(L &op, S &sp, P &prec, real_type reduction, int maxit, int verbose)
Set up loop solver.
Definition: solvers.hh:119
BiCGSTABSolver(L &op, P &prec, real_type reduction, int maxit, int verbose)
Set up solver.
Definition: solvers.hh:557
double elapsed
Elapsed time in seconds.
Definition: solver.hh:62
int iterations
Number of iterations.
Definition: solver.hh:50
X domain_type
The domain type of the operator to be inverted.
Definition: solvers.hh:1449
virtual void post(X &x)=0
Clean up.
X::field_type field_type
The field type of the operator to be inverted.
Definition: solvers.hh:818
X domain_type
The domain type of the operator to be inverted.
Definition: solvers.hh:543
virtual void apply(X &x, X &b, InverseOperatorResult &res)
Apply inverse operator.
Definition: solvers.hh:587
virtual void apply(X &x, X &b, double reduction, InverseOperatorResult &res)
Apply inverse operator with given reduction factor.
Definition: solvers.hh:517
X range_type
The range type of the operator to be inverted.
Definition: solvers.hh:374
X::field_type field_type
The field type of the operator that we do the inverse for.
Definition: solvers.hh:238
Define general, extensible interface for operators. The available implementation wraps a matrix...
Category for sequential solvers.
Definition: solvercategory.hh:21
virtual void apply(X &x, X &b, double reduction, InverseOperatorResult &res)
Apply inverse operator with given reduction factor.
Definition: solvers.hh:1640
virtual void apply(X &x, X &b, InverseOperatorResult &res)
Apply inverse operator,.
Definition: solvers.hh:131
Minimal Residual Method (MINRES)
Definition: solvers.hh:811
virtual void applyscaleadd(field_type alpha, const X &x, Y &y) const =0
apply operator to x, scale and add:
Bi-conjugate Gradient Stabilized (BiCG-STAB)
Definition: solvers.hh:540
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
Define general, extensible interface for inverse operators.
X::field_type field_type
The field type of the operator that we do the inverse for.
Definition: solvers.hh:64
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:1455
MINRESSolver(L &op, P &prec, real_type reduction, int maxit, int verbose)
Set up MINRES solver.
Definition: solvers.hh:828
virtual void pre(X &x, Y &b)=0
Prepare the preconditioner.
double reduction
Reduction achieved: .
Definition: solver.hh:53
X domain_type
The domain type of the operator that we do the inverse for.
Definition: solvers.hh:60
void printHeader(std::ostream &s) const
helper function for printing header of solver output
Definition: solver.hh:121
Abstract base class for all solvers.
Definition: solver.hh:79
GradientSolver(L &op, P &prec, real_type reduction, int maxit, int verbose)
Set up solver.
Definition: solvers.hh:249
X domain_type
The domain type of the operator to be inverted.
Definition: solvers.hh:814
LoopSolver(L &op, P &prec, real_type reduction, int maxit, int verbose)
Set up Loop solver.
Definition: solvers.hh:88
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...
GradientSolver(L &op, S &sp, P &prec, real_type reduction, int maxit, int verbose)
Set up solver.
Definition: solvers.hh:264
conjugate gradient method
Definition: solvers.hh:369
virtual void apply(X &x, Y &b, real_type reduction, InverseOperatorResult &res)
Apply inverse operator.
Definition: solvers.hh:1198
X::field_type field_type
The field type of the operator to be inverted.
Definition: solvers.hh:376
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:820
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 to be inverted.
Definition: solvers.hh:1105
virtual void apply(X &x, X &b, InverseOperatorResult &res)
Apply inverse operator.
Definition: solvers.hh:279
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:1111
X domain_type
The domain type of the operator to be inverted.
Definition: solvers.hh:372
X range_type
The range type of the operator to be inverted.
Definition: solvers.hh:816
virtual void apply(X &x, X &b, double reduction, InverseOperatorResult &res)
Apply inverse operator with given reduction factor.
Definition: solvers.hh:347
virtual void apply(X &x, X &b, double reduction, InverseOperatorResult &res)
Apply inverse operator with given reduction factor.
Definition: solvers.hh:786
void printOutput(std::ostream &s, const DataType &iter, const DataType &norm, const DataType &norm_old) const
helper function for printing solver output
Definition: solver.hh:130
Definition: example.cc:34
virtual void apply(X &x, X &b, InverseOperatorResult &res)
Apply inverse operator.
Definition: solvers.hh:856
virtual void apply(X &x, X &b, InverseOperatorResult &res)
Apply inverse operator.
Definition: solvers.hh:414
X range_type
The range type of the operator to be inverted.
Definition: solvers.hh:545
X domain_type
The domain type of the operator that we do the inverse for.
Definition: solvers.hh:234
X range_type
The range type of the operator that we do the inverse for.
Definition: solvers.hh:236
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:1465
X::field_type field_type
The field type of the operator to be inverted.
Definition: solvers.hh:1453
bool converged
True if convergence criterion has been met.
Definition: solver.hh:56
Default implementation for the scalar case.
Definition: scalarproducts.hh:94
Generalized preconditioned conjugate gradient solver.
Definition: solvers.hh:1445
void clear()
Resets all data.
Definition: solver.hh:40
Statistics about the application of an inverse operator.
Definition: solver.hh:31
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
BiCGSTABSolver(L &op, S &sp, P &prec, real_type reduction, int maxit, int verbose)
Set up solver.
Definition: solvers.hh:572
CGSolver(L &op, P &prec, real_type reduction, int maxit, int verbose)
Set up conjugate gradient solver.
Definition: solvers.hh:386
X::field_type field_type
The field type of the operator to be inverted.
Definition: solvers.hh:547
virtual void apply(X &v, const Y &d)=0
Apply one step of the preconditioner to the system A(v)=d.
virtual void apply(X &x, X &b, double reduction, InverseOperatorResult &res)
apply inverse operator, with given convergence criteria.
Definition: solvers.hh:209