4 #ifndef DUNE_SOLVERS_HH
5 #define DUNE_SOLVERS_HH
19 #include <dune/common/timer.hh>
20 #include <dune/common/ftraits.hh>
21 #include <dune/common/shared_ptr.hh>
22 #include <dune/common/static_assert.hh>
83 template<
class L,
class P>
85 double reduction,
int maxit,
int verbose) :
86 ssp(), _op(op), _prec(prec), _sp(ssp), _reduction(reduction), _maxit(maxit), _verbose(verbose)
88 dune_static_assert(static_cast<int>(L::category) == static_cast<int>(P::category),
89 "L and P have to have the same category!");
91 "L has to be sequential!");
114 template<
class L,
class S,
class P>
116 double reduction,
int maxit,
int verbose) :
117 _op(op), _prec(prec), _sp(sp), _reduction(reduction), _maxit(maxit), _verbose(verbose)
119 dune_static_assert( static_cast<int>(L::category) == static_cast<int>(P::category),
120 "L and P must have the same category!");
121 dune_static_assert( static_cast<int>(L::category) == static_cast<int>(S::category),
122 "L and S must have the same category!");
142 double def0 = _sp.norm(b);
147 std::cout <<
"=== LoopSolver" << std::endl;
159 int i=1;
double def=def0;
160 for ( ; i<=_maxit; i++ )
166 double defnew=_sp.norm(b);
171 if (def<def0*_reduction || def<1E-30)
179 i=std::min(_maxit,i);
197 std::cout <<
"=== rate=" << res.
conv_rate
200 <<
", IT=" << i << std::endl;
207 std::swap(_reduction,reduction);
208 (*this).apply(x,b,res);
209 std::swap(_reduction,reduction);
240 template<
class L,
class P>
242 double reduction,
int maxit,
int verbose) :
243 ssp(), _op(op), _prec(prec), _sp(ssp), _reduction(reduction), _maxit(maxit), _verbose(verbose)
245 dune_static_assert(static_cast<int>(L::category) == static_cast<int>(P::category),
246 "L and P have to have the same category!");
248 "L has to be sequential!");
255 template<
class L,
class S,
class P>
257 double reduction,
int maxit,
int verbose) :
258 _op(op), _prec(prec), _sp(sp), _reduction(reduction), _maxit(maxit), _verbose(verbose)
260 dune_static_assert(static_cast<int>(L::category) == static_cast<int>(P::category),
261 "L and P have to have the same category!");
262 dune_static_assert(static_cast<int>(L::category) == static_cast<int>(S::category),
263 "L and S have to have the same category!");
281 double def0 = _sp.norm(b);
285 std::cout <<
"=== GradientSolver" << std::endl;
293 int i=1;
double def=def0;
295 for ( ; i<=_maxit; i++ )
300 lambda = _sp.dot(p,b)/_sp.dot(q,p);
304 double defnew=_sp.norm(b);
309 if (def<def0*_reduction || def<1E-30)
317 i=std::min(_maxit,i);
328 std::cout <<
"=== rate=" << res.
conv_rate
331 <<
", IT=" << i << std::endl;
341 std::swap(_reduction,reduction);
342 (*this).apply(x,b,res);
343 std::swap(_reduction,reduction);
374 template<
class L,
class P>
375 CGSolver (L& op, P& prec,
double reduction,
int maxit,
int verbose) :
376 ssp(), _op(op), _prec(prec), _sp(ssp), _reduction(reduction), _maxit(maxit), _verbose(verbose)
378 dune_static_assert( static_cast<int>(L::category) == static_cast<int>(P::category),
379 "L and P must have the same category!");
381 "L must be sequential!");
388 template<
class L,
class S,
class P>
389 CGSolver (L& op, S& sp, P& prec,
double reduction,
int maxit,
int verbose) :
390 _op(op), _prec(prec), _sp(sp), _reduction(reduction), _maxit(maxit), _verbose(verbose)
392 dune_static_assert( static_cast<int>(L::category) == static_cast<int>(P::category),
393 "L and P must have the same category!");
394 dune_static_assert( static_cast<int>(L::category) == static_cast<int>(S::category),
395 "L and S must have the same category!");
413 double def0 = _sp.norm(b);
422 std::cout <<
"=== rate=" << res.
conv_rate
424 <<
", IT=0" << std::endl;
430 std::cout <<
"=== CGSolver" << std::endl;
444 rholast = _sp.dot(p,b);
448 for ( ; i<=_maxit; i++ )
452 alpha = _sp.dot(p,q);
453 lambda = rholast/alpha;
458 double defnew=_sp.norm(b);
464 if (def<def0*_reduction || def<1E-30)
481 i=std::min(_maxit,i);
494 std::cout <<
"=== rate=" << res.
conv_rate
497 <<
", IT=" << i << std::endl;
506 virtual void apply (X& x, X& b,
double reduction,
509 std::swap(_reduction,reduction);
510 (*this).apply(x,b,res);
511 std::swap(_reduction,reduction);
537 typedef typename FieldTraits<field_type>::real_type
real_type;
544 template<
class L,
class P>
546 double reduction,
int maxit,
int verbose) :
547 ssp(), _op(op), _prec(prec), _sp(ssp), _reduction(reduction), _maxit(maxit), _verbose(verbose)
549 dune_static_assert(static_cast<int>(L::category) == static_cast<int>(P::category),
"L and P must be of the same category!");
557 template<
class L,
class S,
class P>
559 double reduction,
int maxit,
int verbose) :
560 _op(op), _prec(prec), _sp(sp), _reduction(reduction), _maxit(maxit), _verbose(verbose)
562 dune_static_assert( static_cast<int>(L::category) == static_cast<int>(P::category),
563 "L and P must have the same category!");
564 dune_static_assert( static_cast<int>(L::category) == static_cast<int>(S::category),
565 "L and S must have the same category!");
575 const double EPSILON=1e-80;
577 field_type rho, rho_new, alpha, beta, h, omega;
602 norm = norm_old = norm_0 = _sp.norm(r);
613 std::cout <<
"=== BiCGSTABSolver" << std::endl;
623 if ( norm < (_reduction * norm_0) || norm<1E-30)
638 for (it = 0.5; it < _maxit; it+=.5)
645 rho_new = _sp.dot(rt,r);
648 if (std::abs(rho) <= EPSILON)
649 DUNE_THROW(
ISTLError,
"breakdown in BiCGSTAB - rho "
650 << rho <<
" <= EPSILON " << EPSILON
651 <<
" after " << it <<
" iterations");
652 if (std::abs(omega) <= EPSILON)
653 DUNE_THROW(
ISTLError,
"breakdown in BiCGSTAB - omega "
654 << omega <<
" <= EPSILON " << EPSILON
655 <<
" after " << it <<
" iterations");
662 beta = ( rho_new / rho ) * ( alpha / omega );
678 if ( std::abs(h) < EPSILON )
701 if ( norm < (_reduction * norm_0) )
718 omega = _sp.dot(t,r)/_sp.dot(t,t);
740 if ( norm < (_reduction * norm_0) || norm<1E-30)
750 it=std::min((
double)_maxit,it);
756 res.
iterations =
static_cast<int>(std::ceil(it));
761 std::cout <<
"=== rate=" << res.
conv_rate
764 <<
", IT=" << it << std::endl;
774 std::swap(_reduction,reduction);
775 (*this).apply(x,b,res);
776 std::swap(_reduction,reduction);
805 typedef typename FieldTraits<field_type>::real_type
real_type;
812 template<
class L,
class P>
813 MINRESSolver (L& op, P& prec,
double reduction,
int maxit,
int verbose) :
814 ssp(), _op(op), _prec(prec), _sp(ssp), _reduction(reduction), _maxit(maxit), _verbose(verbose)
816 dune_static_assert( static_cast<int>(L::category) == static_cast<int>(P::category),
817 "L and P must have the same category!");
819 "L must be sequential!");
826 template<
class L,
class S,
class P>
827 MINRESSolver (L& op, S& sp, P& prec,
double reduction,
int maxit,
int verbose) :
828 _op(op), _prec(prec), _sp(sp), _reduction(reduction), _maxit(maxit), _verbose(verbose)
830 dune_static_assert( static_cast<int>(L::category) == static_cast<int>(P::category),
831 "L and P must have the same category!");
832 dune_static_assert( static_cast<int>(L::category) == static_cast<int>(S::category),
833 "L and S must have the same category!");
858 std::cout <<
"=== rate=" << res.
conv_rate <<
", T=" << res.
elapsed <<
", TIT=" << res.
elapsed <<
", IT=0" << std::endl;
864 std::cout <<
"=== MINRESSolver" << std::endl;
890 beta = std::sqrt(std::abs(_sp.dot(z,b)));
896 q[0].resize(b.size());
897 q[1].resize(b.size());
898 q[2].resize(b.size());
904 p[0].resize(b.size());
905 p[1].resize(b.size());
906 p[2].resize(b.size());
916 for ( ; i<=_maxit; i++)
926 q[i2].axpy(-beta, q[i0]);
927 alpha = _sp.dot(q[i2],z);
928 q[i2].axpy(-alpha, q[i1]);
931 _prec.
apply(z,q[i2]);
933 beta = std::sqrt(std::abs(_sp.dot(q[i2],z)));
948 T[2] = c[(i+1)%2]*alpha - s[(i+1)%2]*T[1];
949 T[1] = c[(i+1)%2]*T[1] + s[(i+1)%2]*alpha;
956 c[i%2] = 1.0/std::sqrt(T[2]*T[2] + beta*beta);
957 s[i%2] = beta*c[i%2];
961 T[2] = c[i%2]*T[2] + s[i%2]*beta;
964 xi[i%2] = -s[i%2]*xi[(i+1)%2];
965 xi[(i+1)%2] *= c[i%2];
969 p[i2].axpy(-T[1],p[i1]);
970 p[i2].axpy(-T[0],p[i0]);
974 x.axpy(beta0*xi[(i+1)%2], p[i2]);
985 real_type defnew = std::abs(beta0*xi[i%2]);
991 if (def<def0*_reduction || def<1E-30 || i==_maxit)
999 i=std::min(_maxit,i);
1008 res.
elapsed = watch.elapsed();
1012 std::cout <<
"=== rate=" << res.
conv_rate
1015 <<
", IT=" << i << std::endl;
1027 std::swap(_reduction,reduction);
1028 (*this).apply(x,b,res);
1029 std::swap(_reduction,reduction);
1053 template<
class X,
class Y=X,
class F = Y>
1064 typedef typename FieldTraits<field_type>::real_type
real_type;
1075 template<
class L,
class P>
1076 RestartedGMResSolver (L& op, P& prec,
double reduction,
int restart,
int maxit,
int verbose,
bool recalc_defect =
false) :
1078 ssp(), _sp(ssp), _restart(restart),
1079 _reduction(reduction), _maxit(maxit), _verbose(verbose),
1080 _recalc_defect(recalc_defect)
1082 dune_static_assert(static_cast<int>(P::category) == static_cast<int>(L::category),
1083 "P and L must be the same category!");
1085 "L must be sequential!");
1095 template<
class L,
class S,
class P>
1096 RestartedGMResSolver (L& op, S& sp, P& prec,
double reduction,
int restart,
int maxit,
int verbose,
bool recalc_defect =
false) :
1098 _sp(sp), _restart(restart),
1099 _reduction(reduction), _maxit(maxit), _verbose(verbose),
1100 _recalc_defect(recalc_defect)
1102 dune_static_assert(static_cast<int>(P::category) == static_cast<int>(L::category),
1103 "P and L must have the same category!");
1104 dune_static_assert(static_cast<int>(P::category) == static_cast<int>(S::category),
1105 "P and S must have the same category!");
1111 apply(x,b,_reduction,res);
1127 std::vector<field_type> s(m+1), cs(m), sn(m);
1130 std::vector< std::vector<field_type> > H(m+1,s);
1131 std::vector<F> v(m+1,b);
1142 w = 0.0; _M.
apply(w,b);
1143 norm_0 = _sp.norm(w);
1147 v[0] = 0.0; _M.
apply(v[0],w);
1148 beta = _sp.norm(v[0]);
1154 v[0] = 0.0; _M.
apply(v[0],b);
1155 beta = _sp.norm(v[0]);
1162 norm = norm_old = beta;
1167 std::cout <<
"=== RestartedGMResSolver" << std::endl;
1176 if (norm <= reduction * norm_0) {
1184 while (j <= _maxit && res.
converged !=
true) {
1185 v[0] *= (1.0 / beta);
1186 for (i=1; i<=m; i++) s[i] = 0.0;
1188 int end=std::min(m, _maxit-j+1);
1189 for (i = 0; i < end && res.
converged !=
true; i++, j++) {
1192 _A_.
apply(v[i], v[i+1]);
1193 _M.
apply(w, v[i+1]);
1194 for (k = 0; k <= i; k++) {
1195 H[k][i] = _sp.dot(w, v[k]);
1197 w.axpy(-H[k][i], v[k]);
1199 H[i+1][i] = _sp.norm(w);
1200 if (H[i+1][i] == 0.0)
1201 DUNE_THROW(
ISTLError,
"breakdown in GMRes - |w| "
1202 <<
" == 0.0 after " << j <<
" iterations");
1204 v[i+1] = w; v[i+1] *= (1.0 / H[i+1][i]);
1206 for (k = 0; k < i; k++)
1207 applyPlaneRotation(H[k][i], H[k+1][i], cs[k], sn[k]);
1209 generatePlaneRotation(H[i][i], H[i+1][i], cs[i], sn[i]);
1210 applyPlaneRotation(H[i][i], H[i+1][i], cs[i], sn[i]);
1211 applyPlaneRotation(s[i], s[i+1], cs[i], sn[i]);
1213 norm = std::abs(s[i+1]);
1222 if (norm < reduction * norm_0) {
1230 update(x, i - 1, H, s, v);
1236 beta = _sp.norm(v[0]);
1243 update(w, i - 1, H, s, v);
1252 v[0] = 0.0; _M.
apply(v[0],b);
1253 beta = _sp.norm(v[0]);
1261 if (norm < reduction * norm_0) {
1266 if (res.
converged !=
true && _verbose > 0)
1267 std::cout <<
"=== GMRes::restart\n";
1275 res.
elapsed = watch.elapsed();
1286 std::cout <<
"=== rate=" << res.
conv_rate
1295 const std::vector< std::vector<field_type> > & h,
1296 const std::vector<field_type> & s,
const std::vector<F> &v)
1298 std::vector<field_type> y(s);
1301 for (
int i = k; i >= 0; i--) {
1303 for (
int j = i - 1; j >= 0; j--)
1304 y[j] -= h[j][i] * y[i];
1307 for (
int j = 0; j <= k; j++)
1318 }
else if (std::abs(dy) > std::abs(dx)) {
1320 sn = 1.0 / std::sqrt( 1.0 + temp*temp );
1324 cs = 1.0 / std::sqrt( 1.0 + temp*temp );
1334 dy = -sn * dx + cs * dy;
1338 LinearOperator<X,X>& _A_;
1339 Preconditioner<X,X>& _M;
1340 SeqScalarProduct<X> ssp;
1341 ScalarProduct<X>& _sp;
1346 bool _recalc_defect;
1381 template<
class L,
class P>
1384 ssp(), _op(op), _prec(prec), _sp(ssp), _reduction(reduction), _maxit(maxit),
1385 _verbose(verbose), _restart(std::min(maxit,restart))
1387 dune_static_assert(static_cast<int>(L::category) == static_cast<int>(P::category),
1388 "L and P have to have the same category!");
1389 dune_static_assert(static_cast<int>(L::category) ==
1391 "L has to be sequential!");
1400 template<
class L,
class P,
class S>
1402 double reduction,
int maxit,
int verbose,
int restart=10) :
1403 _op(op), _prec(prec), _sp(sp), _reduction(reduction), _maxit(maxit), _verbose(verbose),
1404 _restart(std::min(maxit,restart))
1406 dune_static_assert( static_cast<int>(L::category) == static_cast<int>(P::category),
1407 "L and P must have the same category!");
1408 dune_static_assert( static_cast<int>(L::category) == static_cast<int>(S::category),
1409 "L and S must have the same category!");
1423 std::vector<shared_ptr<X> > p(_restart);
1424 std::vector<typename X::field_type> pp(_restart);
1428 p[0].reset(
new X(x));
1430 double def0 = _sp.norm(b);
1439 std::cout <<
"=== rate=" << res.
conv_rate
1441 <<
", IT=0" << std::endl;
1447 std::cout <<
"=== GeneralizedPCGSolver" << std::endl;
1461 _prec.
apply(*(p[0]),b);
1462 rho = _sp.dot(*(p[0]),b);
1463 _op.
apply(*(p[0]),q);
1464 pp[0] = _sp.dot(*(p[0]),q);
1466 x.axpy(lambda,*(p[0]));
1470 double defnew=_sp.norm(b);
1474 if (def<def0*_reduction || def<1E-30)
1479 std::cout <<
"=== rate=" << res.
conv_rate
1482 <<
", IT=" << 1 << std::endl;
1489 int end=std::min(_restart, _maxit-i+1);
1490 for (ii=1; ii<end; ++ii )
1495 _prec.
apply(prec_res,b);
1497 p[ii].reset(
new X(prec_res));
1498 _op.
apply(prec_res, q);
1500 for(
int j=0; j<ii; ++j) {
1501 rho =_sp.dot(q,*(p[j]))/pp[j];
1502 p[ii]->axpy(-rho, *(p[j]));
1506 _op.
apply(*(p[ii]),q);
1507 pp[ii] = _sp.dot(*(p[ii]),q);
1508 rho = _sp.dot(*(p[ii]),b);
1509 lambda = rho/pp[ii];
1510 x.axpy(lambda,*(p[ii]));
1514 double defnew=_sp.norm(b);
1520 if (def<def0*_reduction || def<1E-30)
1529 *(p[0])=*(p[_restart-1]);
1530 pp[0]=pp[_restart-1];
1541 res.
elapsed = watch.elapsed();
1545 std::cout <<
"=== rate=" << res.
conv_rate
1548 <<
", IT=" << i+1 << std::endl;
1557 virtual void apply (X& x, X& b,
double reduction,
1560 std::swap(_reduction,reduction);
1561 (*this).apply(x,b,res);
1562 std::swap(_reduction,reduction);