Skip to content
Snippets Groups Projects
Commit 7fb6aa4d authored by Christian Engwer's avatar Christian Engwer
Browse files

Merge branch 'feature/smaller-cleanups' into 'master' 

Feature/smaller cleanups

as the title suggests, these are a couple of smaller fixes, which are not directly related and aren't important enough to manage them in separate branches.

See merge request !13
parents 1299ddb1 37b668cd
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
dune/istl/bcrsmatrix.hh
+ 2
− 2
View file @ 7fb6aa4d
......@@ -1623,8 +1623,8 @@ namespace Dune {
}
//! y += alpha A x
template<class X, class Y>
void usmv (const field_type& alpha, const X& x, Y& y) const
template<typename F, class X, class Y>
void usmv (F&& alpha, const X& x, Y& y) const
{
#ifdef DUNE_ISTL_WITH_CHECKING
if (ready != built)
......
dune/istl/solver.hh
+ 4
− 4
View file @ 7fb6aa4d
......@@ -126,9 +126,9 @@ namespace Dune
}
//! helper function for printing solver output
template <class DataType>
template <typename CountType, typename DataType>
void printOutput(std::ostream& s,
const DataType& iter,
const CountType& iter,
const DataType& norm,
const DataType& norm_old) const
{
......@@ -139,9 +139,9 @@ namespace Dune
}
//! helper function for printing solver output
template <class DataType>
template <typename CountType, typename DataType>
void printOutput(std::ostream& s,
const DataType& iter,
const CountType& iter,
const DataType& norm) const
{
s << std::setw(iterationSpacing) << iter << " ";
......
dune/istl/solvers.hh
+ 48
− 38
View file @ 7fb6aa4d
......@@ -152,7 +152,7 @@ namespace Dune {
if (_verbose>1)
{
this->printHeader(std::cout);
this->printOutput(std::cout,real_type(0),def0);
this->printOutput(std::cout,0,def0);
}
}
......@@ -169,7 +169,7 @@ namespace Dune {
_op.applyscaleadd(-1,v,b); // update defect
real_type defnew=_sp.norm(b); // comp defect norm
if (_verbose>1) // print
this->printOutput(std::cout,real_type(i),defnew,def);
this->printOutput(std::cout,i,defnew,def);
//std::cout << i << " " << defnew << " " << defnew/def << std::endl;
def = defnew; // update norm
if (def<def0*_reduction || def<1E-30) // convergence check
......@@ -184,7 +184,7 @@ namespace Dune {
// print
if (_verbose==1)
this->printOutput(std::cout,real_type(i),def);
this->printOutput(std::cout,i,def);
// postprocess preconditioner
_prec.post(x);
......@@ -294,7 +294,7 @@ namespace Dune {
if (_verbose>1)
{
this->printHeader(std::cout);
this->printOutput(std::cout,real_type(0),def0);
this->printOutput(std::cout,0,def0);
}
}
......@@ -311,7 +311,7 @@ namespace Dune {
real_type defnew=_sp.norm(b); // comp defect norm
if (_verbose>1) // print
this->printOutput(std::cout,real_type(i),defnew,def);
this->printOutput(std::cout,i,defnew,def);
def = defnew; // update norm
if (def<def0*_reduction || def<1E-30) // convergence check
......@@ -325,7 +325,7 @@ namespace Dune {
i=std::min(_maxit,i);
if (_verbose==1) // printing for non verbose
this->printOutput(std::cout,real_type(i),def);
this->printOutput(std::cout,i,def);
_prec.post(x); // postprocess preconditioner
res.iterations = i; // fill statistics
......@@ -586,6 +586,7 @@ namespace Dune {
*/
virtual void apply (X& x, X& b, InverseOperatorResult& res)
{
using std::abs;
const real_type EPSILON=1e-80;
double it;
field_type rho, rho_new, alpha, beta, h, omega;
......@@ -628,7 +629,7 @@ namespace Dune {
if (_verbose>1)
{
this->printHeader(std::cout);
this->printOutput(std::cout,real_type(0),norm_0);
this->printOutput(std::cout,0,norm_0);
//std::cout << " Iter Defect Rate" << std::endl;
//std::cout << " 0" << std::setw(14) << norm_0 << std::endl;
}
......@@ -659,11 +660,11 @@ namespace Dune {
rho_new = _sp.dot(rt,r);
// look if breakdown occured
if (std::abs(rho) <= EPSILON)
if (abs(rho) <= EPSILON)
DUNE_THROW(ISTLError,"breakdown in BiCGSTAB - rho "
<< rho << " <= EPSILON " << EPSILON
<< " after " << it << " iterations");
if (std::abs(omega) <= EPSILON)
if (abs(omega) <= EPSILON)
DUNE_THROW(ISTLError,"breakdown in BiCGSTAB - omega "
<< omega << " <= EPSILON " << EPSILON
<< " after " << it << " iterations");
......@@ -689,8 +690,10 @@ namespace Dune {
// alpha = rho_new / < rt, v >
h = _sp.dot(rt,v);
if ( std::abs(h) < EPSILON )
DUNE_THROW(ISTLError,"h=0 in BiCGSTAB");
if (abs(h) < EPSILON)
DUNE_THROW(ISTLError,"abs(h) < EPSILON in BiCGSTAB - abs(h) "
<< abs(h) << " < EPSILON " << EPSILON
<< " after " << it << " iterations");
alpha = rho_new / h;
......@@ -709,7 +712,7 @@ namespace Dune {
if (_verbose>1) // print
{
this->printOutput(std::cout,real_type(it),norm,norm_old);
this->printOutput(std::cout,it,norm,norm_old);
}
if ( norm < (_reduction * norm_0) )
......@@ -748,7 +751,7 @@ namespace Dune {
if (_verbose > 1) // print
{
this->printOutput(std::cout,real_type(it),norm,norm_old);
this->printOutput(std::cout,it,norm,norm_old);
}
if ( norm < (_reduction * norm_0) || norm<1E-30)
......@@ -764,7 +767,7 @@ namespace Dune {
it=std::min(static_cast<double>(_maxit),it);
if (_verbose==1) // printing for non verbose
this->printOutput(std::cout,real_type(it),norm);
this->printOutput(std::cout,it,norm);
_prec.post(x); // postprocess preconditioner
res.iterations = static_cast<int>(std::ceil(it)); // fill statistics
......@@ -855,6 +858,8 @@ namespace Dune {
*/
virtual void apply (X& x, X& b, InverseOperatorResult& res)
{
using std::sqrt;
using std::abs;
// clear solver statistics
res.clear();
// start a timer
......@@ -873,7 +878,7 @@ namespace Dune {
std::cout << "=== MINRESSolver" << std::endl;
if(_verbose > 1) {
this->printHeader(std::cout);
this->printOutput(std::cout,real_type(0),def0);
this->printOutput(std::cout,0,def0);
}
}
......@@ -898,9 +903,9 @@ namespace Dune {
real_type def = def0;
// recurrence coefficients as computed in Lanczos algorithm
field_type alpha, beta;
// diagonal entries of givens rotation
// diagonal entries of givens rotation
std::array<real_type,2> c{{0.0,0.0}};
// off-diagonal entries of givens rotation
// off-diagonal entries of givens rotation
std::array<field_type,2> s{{0.0,0.0}};
// recurrence coefficients (column k of tridiag matrix T_k)
......@@ -918,7 +923,7 @@ namespace Dune {
// beta is real and positive in exact arithmetic
// since it is the norm of the basis vectors (in unpreconditioned case)
beta = std::sqrt(_sp.dot(b,z));
beta = sqrt(_sp.dot(b,z));
field_type beta0 = beta;
// the search directions
......@@ -958,7 +963,7 @@ namespace Dune {
// beta is real and positive in exact arithmetic
// since it is the norm of the basis vectors (in unpreconditioned case)
beta = std::sqrt(_sp.dot(q[i2],z));
beta = sqrt(_sp.dot(q[i2],z));
q[i2] *= 1.0/beta;
z *= 1.0/beta;
......@@ -999,10 +1004,10 @@ namespace Dune {
// check for convergence
// the last entry in the rhs of the min-problem is the residual
real_type defnew = std::abs(beta0*xi[i%2]);
real_type defnew = abs(beta0*xi[i%2]);
if(_verbose > 1)
this->printOutput(std::cout,real_type(i),defnew,def);
this->printOutput(std::cout,i,defnew,def);
def = defnew;
if(def < def0*_reduction || def < 1e-30 || i == _maxit ) {
......@@ -1012,7 +1017,7 @@ namespace Dune {
} // end for
if(_verbose == 1)
this->printOutput(std::cout,real_type(i),def);
this->printOutput(std::cout,i,def);
// postprocess preconditioner
_prec.post(x);
......@@ -1048,8 +1053,10 @@ namespace Dune {
void generateGivensRotation(field_type &dx, field_type &dy, real_type &cs, field_type &sn)
{
real_type norm_dx = std::abs(dx);
real_type norm_dy = std::abs(dy);
using std::sqrt;
using std::abs;
real_type norm_dx = abs(dx);
real_type norm_dy = abs(dy);
if(norm_dy < 1e-15) {
cs = 1.0;
sn = 0.0;
......@@ -1058,7 +1065,7 @@ namespace Dune {
sn = 1.0;
} else if(norm_dy > norm_dx) {
real_type temp = norm_dx/norm_dy;
cs = 1.0/std::sqrt(1.0 + temp*temp);
cs = 1.0/sqrt(1.0 + temp*temp);
sn = cs;
cs *= temp;
sn *= dx/norm_dx;
......@@ -1067,7 +1074,7 @@ namespace Dune {
sn *= dy/norm_dy;
} else {
real_type temp = norm_dy/norm_dx;
cs = 1.0/std::sqrt(1.0 + temp*temp);
cs = 1.0/sqrt(1.0 + temp*temp);
sn = cs;
sn *= dy/dx;
// dy and dx is real in exact arithmetic
......@@ -1197,6 +1204,7 @@ namespace Dune {
*/
virtual void apply (X& x, Y& b, double reduction, InverseOperatorResult& res)
{
using std::abs;
const real_type EPSILON = 1e-80;
const int m = _restart;
real_type norm, norm_old = 0.0, norm_0;
......@@ -1232,7 +1240,7 @@ namespace Dune {
std::cout << "=== RestartedGMResSolver" << std::endl;
if(_verbose > 1) {
this->printHeader(std::cout);
this->printOutput(std::cout,real_type(0),norm_0);
this->printOutput(std::cout,0,norm_0);
}
}
......@@ -1261,14 +1269,14 @@ namespace Dune {
for(int k=0; k<i+1; k++) {
// notice that _sp.dot(v[k],w) = v[k]\adjoint w
// so one has to pay attention to the order
// the in scalar product for the complex case
// in the scalar product for the complex case
// doing the modified Gram-Schmidt algorithm
H[k][i] = _sp.dot(v[k],w);
// w -= H[k][i] * v[k]
w.axpy(-H[k][i],v[k]);
}
H[i+1][i] = _sp.norm(w);
if(std::abs(H[i+1][i]) < EPSILON)
if(abs(H[i+1][i]) < EPSILON)
DUNE_THROW(ISTLError,
"breakdown in GMRes - |w| == 0.0 after " << j << " iterations");
......@@ -1286,11 +1294,11 @@ namespace Dune {
applyPlaneRotation(s[i],s[i+1],cs[i],sn[i]);
// norm of the defect is the last component the vector s
norm = std::abs(s[i+1]);
norm = abs(s[i+1]);
// print current iteration statistics
if(_verbose > 1) {
this->printOutput(std::cout,real_type(j),norm,norm_old);
this->printOutput(std::cout,j,norm,norm_old);
}
norm_old = norm;
......@@ -1385,8 +1393,10 @@ namespace Dune {
void
generatePlaneRotation(field_type &dx, field_type &dy, real_type &cs, field_type &sn)
{
real_type norm_dx = std::abs(dx);
real_type norm_dy = std::abs(dy);
using std::sqrt;
using std::abs;
real_type norm_dx = abs(dx);
real_type norm_dy = abs(dy);
if(norm_dy < 1e-15) {
cs = 1.0;
sn = 0.0;
......@@ -1395,14 +1405,14 @@ namespace Dune {
sn = 1.0;
} else if(norm_dy > norm_dx) {
real_type temp = norm_dx/norm_dy;
cs = 1.0/std::sqrt(1.0 + temp*temp);
cs = 1.0/sqrt(1.0 + temp*temp);
sn = cs;
cs *= temp;
sn *= dx/norm_dx;
sn *= conjugate(dy)/norm_dy;
} else {
real_type temp = norm_dy/norm_dx;
cs = 1.0/std::sqrt(1.0 + temp*temp);
cs = 1.0/sqrt(1.0 + temp*temp);
sn = cs;
sn *= conjugate(dy/dx);
}
......@@ -1530,7 +1540,7 @@ namespace Dune {
std::cout << "=== GeneralizedPCGSolver" << std::endl;
if (_verbose>1) {
this->printHeader(std::cout);
this->printOutput(std::cout,real_type(0),def0);
this->printOutput(std::cout,0,def0);
}
}
// some local variables
......@@ -1552,7 +1562,7 @@ namespace Dune {
// convergence test
real_type defnew=_sp.norm(b); // comp defect norm
if (_verbose>1) // print
this->printOutput(std::cout,real_type(++i),defnew,def);
this->printOutput(std::cout,++i,defnew,def);
def = defnew; // update norm
if (def<def0*_reduction || def<1E-30) // convergence check
{
......@@ -1597,7 +1607,7 @@ namespace Dune {
real_type defnew=_sp.norm(b); // comp defect norm
if (_verbose>1) // print
this->printOutput(std::cout,real_type(++i),defnew,def);
this->printOutput(std::cout,++i,defnew,def);
def = defnew; // update norm
if (def<def0*_reduction || def<1E-30) // convergence check
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment