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!");
284 _op.applyscaleadd(-1,x,b);
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!");
419 _op.applyscaleadd(-1,x,b);
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!");
590 const real_type EPSILON=1e-80;
592 field_type rho, rho_new, alpha, beta, h, omega;
593 real_type norm, norm_old, norm_0;
613 _op.applyscaleadd(-1,x,r);
617 norm = norm_old = norm_0 = _sp.norm(r);
628 std::cout <<
"=== BiCGSTABSolver" << std::endl;
638 if ( norm < (_reduction * norm_0) || norm<1E-30)
653 for (it = 0.5; it < _maxit; it+=.5)
660 rho_new = _sp.dot(rt,r);
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");
677 beta = ( rho_new / rho ) * ( alpha / omega );
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");
718 if ( norm < (_reduction * norm_0) )
735 omega = _sp.dot(t,r)/_sp.dot(t,t);
757 if ( norm < (_reduction * norm_0) || norm<1E-30)
767 it=std::min(static_cast<double>(_maxit),it);
773 res.
iterations =
static_cast<int>(std::ceil(it));
774 res.
reduction =
static_cast<double>(norm/norm_0);
778 std::cout <<
"=== rate=" << res.
conv_rate 781 <<
", IT=" << it << std::endl;
791 real_type saved_reduction = _reduction;
792 _reduction = reduction;
793 (*this).apply(x,b,res);
794 _reduction = saved_reduction;
802 real_type _reduction;
823 typedef typename FieldTraits<field_type>::real_type
real_type;
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)
834 static_assert(static_cast<int>(L::category) == static_cast<int>(P::category),
835 "L and P must have the same category!");
837 "L must be sequential!");
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)
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!");
871 _op.applyscaleadd(-1,x,b);
874 real_type def0 = _sp.norm(b);
878 std::cout <<
"=== MINRESSolver" << std::endl;
894 std::cout <<
"=== rate=" << res.
conv_rate 903 real_type def = def0;
905 field_type alpha, beta;
907 Dune::array<real_type,2> c{{0.0,0.0}};
909 Dune::array<field_type,2> s{{0.0,0.0}};
912 Dune::array<field_type,3> T{{0.0,0.0,0.0}};
915 Dune::array<field_type,2> xi{{1.0,0.0}};
926 beta = sqrt(_sp.dot(b,z));
927 field_type beta0 = beta;
930 Dune::array<X,3> p{{b,b,b}};
936 Dune::array<X,3> q{{b,b,b}};
945 for( ; i<=_maxit; i++) {
954 q[i2].axpy(-beta,q[i0]);
958 alpha = _sp.dot(z,q[i2]);
959 q[i2].axpy(-alpha,q[i1]);
962 _prec.apply(z,q[i2]);
966 beta = sqrt(_sp.dot(q[i2],z));
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;
986 generateGivensRotation(T[2],beta,c[i%2],s[i%2]);
988 T[2] = c[i%2]*T[2] + s[i%2]*beta;
990 xi[i%2] = -s[i%2]*xi[(i+1)%2];
991 xi[(i+1)%2] *= c[i%2];
995 p[i2].axpy(-T[1],p[i1]);
996 p[i2].axpy(-T[0],p[i0]);
1000 x.axpy(beta0*xi[(i+1)%2],p[i2]);
1007 real_type defnew = abs(beta0*xi[i%2]);
1013 if(def < def0*_reduction || def < 1e-30 || i == _maxit ) {
1026 res.
reduction =
static_cast<double>(def/def0);
1028 res.
elapsed = watch.elapsed();
1032 std::cout <<
"=== rate=" << res.
conv_rate 1035 <<
", IT=" << i << std::endl;
1046 real_type saved_reduction = _reduction;
1047 _reduction = reduction;
1048 (*this).apply(x,b,res);
1049 _reduction = saved_reduction;
1054 void generateGivensRotation(field_type &dx, field_type &dy, real_type &cs, field_type &sn)
1058 real_type norm_dx = abs(dx);
1059 real_type norm_dy = abs(dy);
1060 if(norm_dy < 1e-15) {
1063 }
else if(norm_dx < 1e-15) {
1066 }
else if(norm_dy > norm_dx) {
1067 real_type temp = norm_dx/norm_dy;
1068 cs = 1.0/sqrt(1.0 + temp*temp);
1076 real_type temp = norm_dy/norm_dx;
1077 cs = 1.0/sqrt(1.0 + temp*temp);
1089 real_type _reduction;
1107 template<
class X,
class Y=X,
class F = Y>
1118 typedef typename FieldTraits<field_type>::real_type
real_type;
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)
1130 , _reduction(reduction)
1134 static_assert(static_cast<int>(P::category) == static_cast<int>(L::category),
1135 "P and L must be the same category!");
1137 "L must be sequential!");
1147 template<
class L,
class P>
1150 ssp(), _sp(ssp), _restart(restart),
1151 _reduction(reduction), _maxit(maxit), _verbose(verbose)
1153 static_assert(static_cast<int>(P::category) == static_cast<int>(L::category),
1154 "P and L must be the same category!");
1156 "L must be sequential!");
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)
1166 , _reduction(reduction)
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!");
1182 template<
class L,
class S,
class P>
1185 _sp(sp), _restart(restart),
1186 _reduction(reduction), _maxit(maxit), _verbose(verbose)
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!");
1197 apply(x,b,_reduction,res);
1209 const int m = _restart;
1212 std::vector<field_type> s(m+1), sn(m);
1213 std::vector<real_type> cs(m);
1218 std::vector< std::vector<field_type> > H(m+1,s);
1219 std::vector<F> v(m+1,b);
1230 _A.applyscaleadd(-1.0,x,b);
1232 v[0] = 0.0; _W.apply(v[0],b);
1233 norm_0 = _sp.norm(v[0]);
1240 std::cout <<
"=== RestartedGMResSolver" << std::endl;
1247 if(norm_0 < EPSILON) {
1254 while(j <= _maxit && res.
converged !=
true) {
1259 for(i=1; i<m+1; i++)
1262 for(i=0; i < m && j <= _maxit && res.
converged !=
true; i++, j++) {
1267 _A.apply(v[i],v[i+1]);
1269 for(
int k=0; k<i+1; k++) {
1274 H[k][i] = _sp.dot(v[k],w);
1276 w.axpy(-H[k][i],v[k]);
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");
1284 v[i+1] = w; v[i+1] *= 1.0/H[i+1][i];
1287 for(
int k=0; k<i; k++)
1288 applyPlaneRotation(H[k][i],H[k+1][i],cs[k],sn[k]);
1291 generatePlaneRotation(H[i][i],H[i+1][i],cs[i],sn[i]);
1293 applyPlaneRotation(H[i][i],H[i+1][i],cs[i],sn[i]);
1294 applyPlaneRotation(s[i],s[i+1],cs[i],sn[i]);
1307 if(norm < reduction * norm_0)
1321 if( res.
converged !=
true && j <= _maxit ) {
1324 std::cout <<
"=== GMRes::restart" << std::endl;
1328 _A.applyscaleadd(-1.0,x,b);
1332 norm = _sp.norm(v[0]);
1343 res.
reduction =
static_cast<double>(norm/norm_0);
1345 res.
elapsed = watch.elapsed();
1356 std::cout <<
"=== rate=" << res.
conv_rate 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) {
1368 std::vector<field_type> y(s);
1371 for(
int a=i-1; a>=0; a--) {
1373 for(
int b=a+1; b<i; b++)
1374 rhs -= H[a][b]*y[b];
1383 template<
typename T>
1384 typename enable_if<is_same<field_type,real_type>::value,T>::type conjugate(
const T& t) {
1388 template<
typename T>
1389 typename enable_if<!is_same<field_type,real_type>::value,T>::type conjugate(
const T& t) {
1400 if(norm_dy < 1e-15) {
1403 }
else if(norm_dx < 1e-15) {
1406 }
else if(norm_dy > norm_dx) {
1408 cs = 1.0/sqrt(1.0 + temp*temp);
1412 sn *= conjugate(dy)/norm_dy;
1415 cs = 1.0/sqrt(1.0 + temp*temp);
1417 sn *= conjugate(dy/dx);
1426 dy = -conjugate(sn) * dx + cs * dy;
1465 typedef typename FieldTraits<field_type>::real_type
real_type;
1474 template<
class L,
class P>
1477 ssp(), _op(op), _prec(prec), _sp(ssp), _reduction(reduction), _maxit(maxit),
1478 _verbose(verbose), _restart(
std::min(maxit,restart))
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) ==
1484 "L has to be sequential!");
1493 template<
class L,
class P,
class S>
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))
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!");
1514 _op.applyscaleadd(-1,x,b);
1516 std::vector<std::shared_ptr<X> > p(_restart);
1517 std::vector<typename X::field_type> pp(_restart);
1521 p[0].reset(
new X(x));
1523 real_type def0 = _sp.norm(b);
1532 std::cout <<
"=== rate=" << res.
conv_rate 1534 <<
", IT=0" << std::endl;
1540 std::cout <<
"=== GeneralizedPCGSolver" << std::endl;
1548 field_type rho, lambda;
1554 _prec.apply(*(p[0]),b);
1555 rho = _sp.dot(*(p[0]),b);
1556 _op.apply(*(p[0]),q);
1557 pp[0] = _sp.dot(*(p[0]),q);
1559 x.axpy(lambda,*(p[0]));
1563 real_type defnew=_sp.norm(b);
1567 if (def<def0*_reduction || def<1E-30)
1572 std::cout <<
"=== rate=" << res.
conv_rate 1575 <<
", IT=" << 1 << std::endl;
1582 int end=std::min(_restart, _maxit-i+1);
1583 for (ii=1; ii<end; ++ii )
1588 _prec.apply(prec_res,b);
1590 p[ii].reset(
new X(prec_res));
1591 _op.apply(prec_res, q);
1593 for(
int j=0; j<ii; ++j) {
1594 rho =_sp.dot(q,*(p[j]))/pp[j];
1595 p[ii]->axpy(-rho, *(p[j]));
1599 _op.apply(*(p[ii]),q);
1600 pp[ii] = _sp.dot(*(p[ii]),q);
1601 rho = _sp.dot(*(p[ii]),b);
1602 lambda = rho/pp[ii];
1603 x.axpy(lambda,*(p[ii]));
1607 real_type defNew=_sp.norm(b);
1613 if (def<def0*_reduction || def<1E-30)
1622 *(p[0])=*(p[_restart-1]);
1623 pp[0]=pp[_restart-1];
1634 res.
elapsed = watch.elapsed();
1638 std::cout <<
"=== rate=" << res.
conv_rate 1641 <<
", IT=" << i+1 << std::endl;
1650 virtual void apply (X& x, X& b,
double reduction,
1653 real_type saved_reduction = _reduction;
1654 _reduction = reduction;
1655 (*this).apply(x,b,res);
1656 _reduction = saved_reduction;
1663 real_type _reduction;
X range_type
The range type of the operator to be inverted.
Definition: solvers.hh:1461
void clear()
Resets all data.
Definition: solver.hh:40
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:1112
GradientSolver(L &op, S &sp, P &prec, real_type reduction, int maxit, int verbose)
Set up solver.
Definition: solvers.hh:264
implements the Generalized Minimal Residual (GMRes) method
Definition: solvers.hh:1108
X::field_type field_type
The field type of the operator to be inverted.
Definition: solvers.hh:1463
void printHeader(std::ostream &s) const
helper function for printing header of solver output
Definition: solver.hh:121
X::field_type field_type
The field type of the operator to be inverted.
Definition: solvers.hh:1116
double elapsed
Elapsed time in seconds.
Definition: solver.hh:62
X range_type
The range type of the operator that we do the inverse for.
Definition: solvers.hh:62
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
X domain_type
The domain type of the operator to be inverted.
Definition: solvers.hh:372
CGSolver(L &op, P &prec, real_type reduction, int maxit, int verbose)
Set up conjugate gradient solver.
Definition: solvers.hh:386
virtual void apply(X &x, X &b, double reduction, InverseOperatorResult &res)
Apply inverse operator with given reduction factor.
Definition: solvers.hh:1650
Definition: example.cc:34
LoopSolver(L &op, P &prec, real_type reduction, int maxit, int verbose)
Set up Loop solver.
Definition: solvers.hh:88
derive error class from the base class in common
Definition: istlexception.hh:16
X domain_type
The domain type of the operator that we do the inverse for.
Definition: solvers.hh:234
int iterations
Number of iterations.
Definition: solver.hh:50
Default implementation for the scalar case.
Definition: scalarproducts.hh:95
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:378
virtual void apply(X &x, X &b, InverseOperatorResult &res)
Apply inverse operator.
Definition: solvers.hh:279
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:823
F basis_type
The field type of the basis vectors.
Definition: solvers.hh:1120
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:347
Minimal Residual Method (MINRES)
Definition: solvers.hh:814
virtual void apply(X &x, X &b, InverseOperatorResult &res)
Apply inverse operator.
Definition: solvers.hh:1509
X range_type
The range type of the operator to be inverted.
Definition: solvers.hh:819
LoopSolver(L &op, S &sp, P &prec, real_type reduction, int maxit, int verbose)
Set up loop solver.
Definition: solvers.hh:119
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
Define general, extensible interface for inverse operators.
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
gradient method
Definition: solvers.hh:231
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:1465
X range_type
The range type of the operator that we do the inverse for.
Definition: solvers.hh:236
BiCGSTABSolver(L &op, P &prec, real_type reduction, int maxit, int verbose)
Set up solver.
Definition: solvers.hh:557
virtual void apply(X &x, X &b, InverseOperatorResult &res)
Apply inverse operator.
Definition: solvers.hh:859
Y range_type
The range type of the operator to be inverted.
Definition: solvers.hh:1114
X::field_type field_type
The field type of the operator.
Definition: solver.hh:88
MINRESSolver(L &op, S &sp, P &prec, real_type reduction, int maxit, int verbose)
Set up MINRES solver.
Definition: solvers.hh:845
virtual void post(X &x)=0
Clean up.
RestartedGMResSolver(L &op, P &prec, real_type reduction, int restart, int maxit, int verbose)
Set up solver.
Definition: solvers.hh:1148
X::field_type field_type
The field type of the operator to be inverted.
Definition: solvers.hh:376
double reduction
Reduction achieved: .
Definition: solver.hh:53
Define base class for scalar product and norm.
virtual void apply(X &x, X &b, double reduction, InverseOperatorResult &res)
Apply inverse operator with given reduction factor.
Definition: solvers.hh:1044
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
Bi-conjugate Gradient Stabilized (BiCG-STAB)
Definition: solvers.hh:540
virtual void pre(X &x, Y &b)=0
Prepare the preconditioner.
virtual void apply(X &x, X &b, double reduction, InverseOperatorResult &res)
apply inverse operator, with given convergence criteria.
Definition: solvers.hh:209
Statistics about the application of an inverse operator.
Definition: solver.hh:31
X range_type
The range type of the operator to be inverted.
Definition: solvers.hh:374
Abstract base class for all solvers.
Definition: solver.hh:79
virtual void apply(X &x, Y &b, double reduction, InverseOperatorResult &res)
Apply inverse operator.
Definition: solvers.hh:1205
conjugate gradient method
Definition: solvers.hh:369
virtual void apply(X &v, const Y &d)=0
Apply one step of the preconditioner to the system A(v)=d.
RestartedGMResSolver(L &op, S &sp, P &prec, real_type reduction, int restart, int maxit, int verbose)
Set up solver.
Definition: solvers.hh:1183
virtual void apply(X &x, X &b, InverseOperatorResult &res)
Apply inverse operator.
Definition: solvers.hh:587
Base class for scalar product and norm computation.
Definition: scalarproducts.hh:44
Category for sequential solvers.
Definition: solvercategory.hh:21
Generalized preconditioned conjugate gradient solver.
Definition: solvers.hh:1455
virtual void apply(X &x, X &b, double reduction, InverseOperatorResult &res)
Apply inverse operator with given reduction factor.
Definition: solvers.hh:789
bool converged
True if convergence criterion has been met.
Definition: solver.hh:56
virtual void applyscaleadd(field_type alpha, const X &x, Y &y) const =0
apply operator to x, scale and add:
virtual void apply(X &x, X &b, double reduction, InverseOperatorResult &res)
Apply inverse operator with given reduction factor.
Definition: solvers.hh:517
Define general, extensible interface for operators. The available implementation wraps a matrix...
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
X::field_type field_type
The field type of the operator to be inverted.
Definition: solvers.hh:547
X::field_type field_type
The field type of the operator to be inverted.
Definition: solvers.hh:821
virtual void apply(X &x, X &b, InverseOperatorResult &res)
Apply inverse operator,.
Definition: solvers.hh:131
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 apply(X &x, Y &b, InverseOperatorResult &res)
Apply inverse operator,.
Definition: solvers.hh:1195
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
MINRESSolver(L &op, P &prec, real_type reduction, int maxit, int verbose)
Set up MINRES solver.
Definition: solvers.hh:831
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:1459
CGSolver(L &op, S &sp, P &prec, real_type reduction, int maxit, int verbose)
Set up conjugate gradient solver.
Definition: solvers.hh:400
Definition: basearray.hh:19
X range_type
The range type of the operator to be inverted.
Definition: solvers.hh:545
Preconditioned loop solver.
Definition: solvers.hh:57
X domain_type
The domain type of the operator to be inverted.
Definition: solvers.hh:817
X::field_type field_type
The field type of the operator that we do the inverse for.
Definition: solvers.hh:238
double conv_rate
Convergence rate (average reduction per step)
Definition: solver.hh:59