4 #ifndef DUNE_ISTL_SOLVERS_HH 5 #define DUNE_ISTL_SOLVERS_HH 12 #include <type_traits> 15 #include <dune/common/conditional.hh> 16 #include <dune/common/exceptions.hh> 17 #include <dune/common/hybridutilities.hh> 18 #include <dune/common/rangeutilities.hh> 19 #include <dune/common/std/type_traits.hh> 20 #include <dune/common/timer.hh> 82 _op->applyscaleadd(-1,x,b);
90 std::cout <<
"=== LoopSolver" << std::endl;
108 _op->applyscaleadd(-1,v,b);
114 if (all_true(def<def0*
_reduction) || max_value(def)<1E-30)
133 res.
reduction =
static_cast<double>(max_value(def/def0));
140 std::cout <<
"=== rate=" << res.
conv_rate 143 <<
", IT=" << i << std::endl;
183 _op->applyscaleadd(-1,x,b);
192 std::cout <<
"=== GradientSolver" << std::endl;
207 lambda =
_sp->dot(p,b)/
_sp->dot(q,p);
216 if (all_true(def<def0*
_reduction) || max_value(def)<1E-30)
231 res.
reduction =
static_cast<double>(max_value(def/def0));
235 std::cout <<
"=== rate=" << res.
conv_rate 238 <<
", IT=" << i << std::endl;
269 using enableConditionEstimate_t = Dune::Std::bool_constant<(std::is_same<field_type,float>::value || std::is_same<field_type,double>::value)>;
285 condition_estimate_(condition_estimate)
288 condition_estimate_ =
false;
289 std::cerr <<
"WARNING: Condition estimate was disabled. It is only available for double and float field types!" << std::endl;
301 scalar_real_type reduction,
int maxit,
int verbose,
bool condition_estimate) :
IterativeSolver<X,X>(op, sp, prec, reduction, maxit, verbose),
302 condition_estimate_(condition_estimate)
304 if (condition_estimate && !(std::is_same<field_type,float>::value || std::is_same<field_type,double>::value)) {
305 condition_estimate_ =
false;
306 std::cerr <<
"WARNING: Condition estimate was disabled. It is only available for double and float field types!" << std::endl;
329 _op->applyscaleadd(-1,x,b);
336 if (!all_true(isfinite(def0)))
339 std::cout <<
"=== CGSolver: abort due to infinite or NaN initial defect" 341 DUNE_THROW(
SolverAbort,
"CGSolver: initial defect=" << def0
342 <<
" is infinite or NaN");
345 if (max_value(def0)<1E-30)
353 std::cout <<
"=== rate=" << res.
conv_rate 355 <<
", IT=0" << std::endl;
361 std::cout <<
"=== CGSolver" << std::endl;
369 std::vector<real_type> lambdas(0);
370 std::vector<real_type> betas(0);
379 rholast =
_sp->dot(p,b);
387 alpha =
_sp->dot(p,q);
388 lambda = rholast/alpha;
390 if (condition_estimate_)
391 lambdas.push_back(std::real(
id(lambda)));
403 if (!all_true(isfinite(def)))
406 std::cout <<
"=== CGSolver: abort due to infinite or NaN defect" 409 "CGSolver: defect=" << def <<
" is infinite or NaN");
412 if (all_true(def<def0*
_reduction) || max_value(def)<1E-30)
424 if (condition_estimate_)
425 betas.push_back(std::real(
id(beta)));
440 res.
reduction =
static_cast<double>(max_value(max_value(def/def0)));
446 std::cout <<
"=== rate=" << res.
conv_rate 449 <<
", IT=" << i << std::endl;
453 if (condition_estimate_) {
461 COND_MAT T(i, i, COND_MAT::row_wise);
465 row.insert(row.index()-1);
466 row.insert(row.index());
467 if (row.index() < T.
N() - 1)
468 row.insert(row.index()+1);
470 for (
int row = 0; row < i; ++row) {
472 T[row][row-1] = std::sqrt(
id(betas[row-1])) / lambdas[row-1];
475 T[row][row] = 1.0 / id(lambdas[row]);
477 T[row][row] += betas[row-1] / lambdas[row-1];
481 T[row][row+1] = std::sqrt(
id(betas[row])) / lambdas[row];
497 std::cout <<
"Min eigv estimate: " << min_eigv << std::endl;
498 std::cout <<
"Max eigv estimate: " << max_eigv << std::endl;
499 std::cout <<
"Condition estimate: " << max_eigv / min_eigv << std::endl;
503 std::cerr <<
"WARNING: Condition estimate was requested. This requires ARPACK, but ARPACK was not found!" << std::endl;
509 bool condition_estimate_ =
false;
554 field_type rho, rho_new, alpha, beta, h, omega;
575 _op->applyscaleadd(-1,x,r);
579 norm = norm_old = norm_0 =
_sp->norm(r);
590 std::cout <<
"=== BiCGSTABSolver" << std::endl;
600 if ( all_true(norm < (
_reduction * norm_0)) || max_value(norm)<1E-30)
615 for (it = 0.5; it <
_maxit; it+=.5)
622 rho_new =
_sp->dot(rt,r);
625 if (all_true(abs(rho) <= EPSILON))
626 DUNE_THROW(
SolverAbort,
"breakdown in BiCGSTAB - rho " 627 << rho <<
" <= EPSILON " << max_value(EPSILON)
628 <<
" after " << it <<
" iterations");
629 if (all_true(abs(omega) <= EPSILON))
630 DUNE_THROW(
SolverAbort,
"breakdown in BiCGSTAB - omega " 631 << omega <<
" <= EPSILON " << max_value(EPSILON)
632 <<
" after " << it <<
" iterations");
639 beta = ( rho_new / rho ) * ( alpha / omega );
655 if ( all_true(abs(h) < EPSILON) )
656 DUNE_THROW(
SolverAbort,
"abs(h) < EPSILON in BiCGSTAB - abs(h) " 657 << abs(h) <<
" < EPSILON " << max_value(EPSILON)
658 <<
" after " << it <<
" iterations");
697 omega =
_sp->dot(t,r)/
_sp->dot(t,t);
719 if ( all_true(norm < (
_reduction * norm_0)) || max_value(norm)<1E-30)
729 it=std::min(static_cast<double>(_maxit),it);
735 res.
iterations =
static_cast<int>(std::ceil(it));
736 res.
reduction =
static_cast<double>(max_value(norm/norm_0));
740 std::cout <<
"=== rate=" << res.
conv_rate 743 <<
", IT=" << it << std::endl;
792 _op->applyscaleadd(-1,x,b);
799 std::cout <<
"=== MINRESSolver" << std::endl;
807 if( max_value(def0) < 1e-30 ) {
815 std::cout <<
"=== rate=" << res.
conv_rate 828 std::array<real_type,2> c{{0.0,0.0}};
830 std::array<field_type,2> s{{0.0,0.0}};
833 std::array<field_type,3> T{{0.0,0.0,0.0}};
836 std::array<field_type,2> xi{{1.0,0.0}};
847 beta = sqrt(
_sp->dot(b,z));
851 std::array<X,3> p{{b,b,b}};
857 std::array<X,3> q{{b,b,b}};
875 q[i2].axpy(-beta,q[i0]);
879 alpha =
_sp->dot(z,q[i2]);
880 q[i2].axpy(-alpha,q[i1]);
883 _prec->apply(z,q[i2]);
887 beta = sqrt(
_sp->dot(q[i2],z));
900 T[2] = c[(i+1)%2]*alpha - s[(i+1)%2]*T[1];
901 T[1] = c[(i+1)%2]*T[1] + s[(i+1)%2]*alpha;
907 generateGivensRotation(T[2],beta,c[i%2],s[i%2]);
909 T[2] = c[i%2]*T[2] + s[i%2]*beta;
911 xi[i%2] = -s[i%2]*xi[(i+1)%2];
912 xi[(i+1)%2] *= c[i%2];
916 p[i2].axpy(-T[1],p[i1]);
917 p[i2].axpy(-T[0],p[i0]);
921 x.axpy(beta0*xi[(i+1)%2],p[i2]);
935 || max_value(def) < 1e-30 || i ==
_maxit ) {
948 res.
reduction =
static_cast<double>(max_value(def/def0));
954 std::cout <<
"=== rate=" << res.
conv_rate 957 <<
", IT=" << i << std::endl;
971 real_type to_real(
const std::complex<real_type> & v)
985 real_type norm_max = max(norm_dx, norm_dy);
986 real_type norm_min = min(norm_dx, norm_dy);
989 cs = cond(norm_dy < eps,
993 cond(norm_dy > norm_dx,
997 sn = cond(norm_dy < eps,
1001 cond(norm_dy > norm_dx,
1033 template<
class X,
class Y=X,
class F = Y>
1099 const int m = _restart;
1102 std::vector<field_type,fAlloc> s(m+1), sn(m);
1103 std::vector<real_type,rAlloc> cs(m);
1108 std::vector< std::vector<field_type,fAlloc> > H(m+1,s);
1109 std::vector<F> v(m+1,b);
1120 _op->applyscaleadd(-1.0,x,b);
1122 v[0] = 0.0;
_prec->apply(v[0],b);
1123 norm_0 =
_sp->norm(v[0]);
1130 std::cout <<
"=== RestartedGMResSolver" << std::endl;
1137 if(all_true(norm_0 < EPSILON)) {
1149 for(i=1; i<m+1; i++)
1157 _op->apply(v[i],v[i+1]);
1158 _prec->apply(w,v[i+1]);
1159 for(
int k=0; k<i+1; k++) {
1164 H[k][i] =
_sp->dot(v[k],w);
1166 w.axpy(-H[k][i],v[k]);
1168 H[i+1][i] =
_sp->norm(w);
1169 if(all_true(abs(H[i+1][i]) < EPSILON))
1171 "breakdown in GMRes - |w| == 0.0 after " << j <<
" iterations");
1174 v[i+1] = w; v[i+1] *=
real_type(1.0)/H[i+1][i];
1177 for(
int k=0; k<i; k++)
1178 applyPlaneRotation(H[k][i],H[k+1][i],cs[k],sn[k]);
1181 generatePlaneRotation(H[i][i],H[i+1][i],cs[i],sn[i]);
1183 applyPlaneRotation(H[i][i],H[i+1][i],cs[i],sn[i]);
1184 applyPlaneRotation(s[i],s[i+1],cs[i],sn[i]);
1197 if(all_true(norm <
real_type(reduction) * norm_0))
1214 std::cout <<
"=== GMRes::restart" << std::endl;
1218 _op->applyscaleadd(-1.0,x,b);
1221 _prec->apply(v[0],b);
1222 norm =
_sp->norm(v[0]);
1233 res.
reduction =
static_cast<double>(max_value(norm/norm_0));
1235 res.
elapsed = watch.elapsed();
1246 std::cout <<
"=== rate=" << res.
conv_rate 1253 void update(X& w,
int i,
1254 const std::vector<std::vector<field_type,fAlloc> >& H,
1255 const std::vector<field_type,fAlloc>& s,
1256 const std::vector<X>& v) {
1258 std::vector<field_type,fAlloc> y(s);
1261 for(
int a=i-1; a>=0; a--) {
1263 for(
int b=a+1; b<i; b++)
1264 rhs -= H[a][b]*y[b];
1273 template<
typename T>
1274 typename std::enable_if<std::is_same<field_type,real_type>::value,T>::type conjugate(
const T& t) {
1278 template<
typename T>
1279 typename std::enable_if<!std::is_same<field_type,real_type>::value,T>::type conjugate(
const T& t) {
1292 real_type to_real(
const std::complex<real_type> & v)
1307 real_type norm_max = max(norm_dx, norm_dy);
1308 real_type norm_min = min(norm_dx, norm_dy);
1311 cs = cond(norm_dy < eps,
1315 cond(norm_dy > norm_dx,
1319 sn = cond(norm_dy < eps,
1323 cond(norm_dy > norm_dx,
1334 dy = -conjugate(sn) * dx + cs * dy;
1414 _op->applyscaleadd(-1,x,b);
1416 std::vector<std::shared_ptr<X> > p(_restart);
1417 std::vector<field_type,fAlloc> pp(_restart);
1421 p[0].reset(
new X(x));
1424 if ( max_value(def0) < 1E-30 )
1432 std::cout <<
"=== rate=" << res.
conv_rate 1434 <<
", IT=0" << std::endl;
1440 std::cout <<
"=== GeneralizedPCGSolver" << std::endl;
1454 _prec->apply(*(p[0]),b);
1455 rho =
_sp->dot(*(p[0]),b);
1456 _op->apply(*(p[0]),q);
1457 pp[0] =
_sp->dot(*(p[0]),q);
1459 x.axpy(lambda,*(p[0]));
1468 if (all_true(def<def0*
_reduction) || max_value(def)<1E-30)
1473 std::cout <<
"=== rate=" << res.
conv_rate 1476 <<
", IT=" << 1 << std::endl;
1483 int end=std::min(_restart,
_maxit-i+1);
1484 for (ii=1; ii<end; ++ii )
1489 _prec->apply(prec_res,b);
1491 p[ii].reset(
new X(prec_res));
1492 _op->apply(prec_res, q);
1494 for(
int j=0; j<ii; ++j) {
1495 rho =
_sp->dot(q,*(p[j]))/pp[j];
1496 p[ii]->axpy(-rho, *(p[j]));
1500 _op->apply(*(p[ii]),q);
1501 pp[ii] =
_sp->dot(*(p[ii]),q);
1502 rho =
_sp->dot(*(p[ii]),b);
1503 lambda = rho/pp[ii];
1504 x.axpy(lambda,*(p[ii]));
1508 defnew=
_sp->norm(b);
1515 if (all_true(def<def0*_reduction) || max_value(def)<1E-30)
1524 *(p[0])=*(p[_restart-1]);
1525 pp[0]=pp[_restart-1];
1534 res.
reduction =
static_cast<double>(max_value(def/def0));
1536 res.
elapsed = watch.elapsed();
1540 std::cout <<
"=== rate=" << res.
conv_rate 1543 <<
", IT=" << i+1 << std::endl;
Base class for all implementations of iterative solvers.
Definition: solver.hh:195
GeneralizedPCGSolver(LinearOperator< X, X > &op, Preconditioner< X, X > &prec, scalar_real_type reduction, int maxit, int verbose, int restart=10)
Set up nonlinear preconditioned conjugate gradient solver.
Definition: solvers.hh:1387
X::field_type field_type
The field type of the operator.
Definition: solver.hh:100
Definition: allocator.hh:7
double reduction
Reduction achieved: .
Definition: solver.hh:62
A sparse block matrix with compressed row storage.
Definition: bcrsmatrix.hh:422
virtual void apply(X &x, X &b, InverseOperatorResult &res)
Apply inverse operator.
Definition: solvers.hh:548
typename AllocatorTraits< T >::type::template rebind< X >::other ReboundAllocatorType
Definition: allocator.hh:33
virtual void apply(X &x, X &b, InverseOperatorResult &res)
Apply inverse operator.
Definition: solvers.hh:780
gradient method
Definition: solvers.hh:160
Define general, extensible interface for inverse operators.
virtual void apply(X &x, Y &b, InverseOperatorResult &res)
Apply inverse operator.
Definition: solvers.hh:1082
SimdScalar< real_type > scalar_real_type
scalar type underlying the field_type
Definition: solver.hh:106
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: solver.hh:103
std::shared_ptr< ScalarProduct< X > > _sp
Definition: solver.hh:302
CreateIterator createbegin()
get initial create iterator
Definition: bcrsmatrix.hh:1056
size_type N() const
number of rows (counted in blocks)
Definition: bcrsmatrix.hh:1894
void printHeader(std::ostream &s) const
helper function for printing header of solver output
Definition: solver.hh:155
A linear operator.
Definition: operators.hh:64
Generalized preconditioned conjugate gradient solver.
Definition: solvers.hh:1362
Define base class for scalar product and norm.
Bi-conjugate Gradient Stabilized (BiCG-STAB)
Definition: solvers.hh:528
RestartedGMResSolver(LinearOperator< X, Y > &op, Preconditioner< X, Y > &prec, scalar_real_type reduction, int restart, int maxit, int verbose)
Set up RestartedGMResSolver solver.
Definition: solvers.hh:1058
std::shared_ptr< LinearOperator< X, X > > _op
Definition: solver.hh:300
scalar_real_type _reduction
Definition: solver.hh:303
virtual void apply(X &x, X &b, InverseOperatorResult &res)
Apply inverse operator.
Definition: solvers.hh:178
Dune::Std::bool_constant<(std::is_same< field_type, float >::value||std::is_same< field_type, double >::value)> enableConditionEstimate_t
Definition: solvers.hh:269
double elapsed
Elapsed time in seconds.
Definition: solver.hh:74
Implementation of the BCRSMatrix class.
Preconditioned loop solver.
Definition: solvers.hh:56
virtual void apply(X &x, X &b, InverseOperatorResult &res)
Apply inverse operator.
Definition: solvers.hh:1409
Define general, extensible interface for operators. The available implementation wraps a matrix...
Thrown when a solver aborts due to some problem.
Definition: istlexception.hh:43
Wrapper to use a range of ARPACK++ eigenvalue solvers.
Definition: arpackpp.hh:241
bool converged
True if convergence criterion has been met.
Definition: solver.hh:65
implements the Generalized Minimal Residual (GMRes) method
Definition: solvers.hh:1034
conjugate gradient method
Definition: solvers.hh:254
CGSolver(LinearOperator< X, X > &op, Preconditioner< X, X > &prec, scalar_real_type reduction, int maxit, int verbose, bool condition_estimate)
Constructor to initialize a CG solver.
Definition: solvers.hh:283
virtual void apply(X &x, X &b, InverseOperatorResult &res)
Apply inverse operator,.
Definition: solvers.hh:70
void computeSymMaxMagnitude(const Real &epsilon, BlockVector &x, Real &lambda) const
Assume the matrix to be square, symmetric and perform IRLM to compute an approximation lambda of its ...
Definition: arpackpp.hh:286
int _maxit
Definition: solver.hh:304
Minimal Residual Method (MINRES)
Definition: solvers.hh:762
CGSolver(LinearOperator< X, X > &op, ScalarProduct< X > &sp, Preconditioner< X, X > &prec, scalar_real_type reduction, int maxit, int verbose, bool condition_estimate)
Constructor to initialize a CG solver.
Definition: solvers.hh:300
X domain_type
Type of the domain of the operator to be inverted.
Definition: solver.hh:94
A vector of blocks with memory management.
Definition: bvector.hh:316
void clear()
Resets all data.
Definition: solver.hh:49
virtual void apply(X &x, X &b, InverseOperatorResult &res)
Apply inverse operator.
Definition: solvers.hh:322
void computeSymMinMagnitude(const Real &epsilon, BlockVector &x, Real &lambda) const
Assume the matrix to be square, symmetric and perform IRLM to compute an approximation lambda of its ...
Definition: arpackpp.hh:388
Base class for scalar product and norm computation.
Definition: scalarproducts.hh:46
virtual void apply(X &x, Y &b, double reduction, InverseOperatorResult &res)
Apply inverse operator.
Definition: solvers.hh:1095
Y range_type
Type of the range of the operator to be inverted.
Definition: solver.hh:97
GeneralizedPCGSolver(LinearOperator< X, X > &op, ScalarProduct< X > &sp, Preconditioner< X, X > &prec, scalar_real_type reduction, int maxit, int verbose, int restart=10)
Set up nonlinear preconditioned conjugate gradient solver.
Definition: solvers.hh:1399
Statistics about the application of an inverse operator.
Definition: solver.hh:40
int _verbose
Definition: solver.hh:305
CreateIterator createend()
get create iterator pointing to one after the last block
Definition: bcrsmatrix.hh:1062
RestartedGMResSolver(LinearOperator< X, Y > &op, ScalarProduct< X > &sp, Preconditioner< X, Y > &prec, scalar_real_type reduction, int restart, int maxit, int verbose)
Set up RestartedGMResSolver solver.
Definition: solvers.hh:1069
int iterations
Number of iterations.
Definition: solver.hh:59
std::shared_ptr< Preconditioner< X, X > > _prec
Definition: solver.hh:301
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:164
double conv_rate
Convergence rate (average reduction per step)
Definition: solver.hh:68
double condition_estimate
Estimate of condition number.
Definition: solver.hh:71