[Dune-devel] [Dune-Commit] [Commit] dune-istl - 1f4cd6b: fixed GMRes
Christian Engwer
christian.engwer at uni-muenster.de
Thu May 8 18:55:05 CEST 2014
Hallo Steffen,
gibt es einen Grund für so viele NOP-changes. Da ist ganz viel
whitespace geändert und/oder Variablen umbenannt worden, was die
patches total kaputt macht. Bisher waren deine Patches da um einiges
aufgeräumter ;-)
Christian
On Thu, May 08, 2014 at 05:42:09PM +0200, Steffen Müthing wrote:
> New commit, appeared at Thu May 8 17:42:09 2014 +0200
> as part of the following ref changes:
>
> branch refs/heads/master updated from 7a48043 -> 42f157e
>
> Browsable version: http://cgit.dune-project.org/repositories/dune-istl/commit/?id=1f4cd6b4b7431e4aeceb828f8785637ea00c1962
>
> ======================================================================
>
> commit 1f4cd6b4b7431e4aeceb828f8785637ea00c1962
> Author: Marian Piatkowski <marian.piatkowski at iwr.uni-heidelberg.de>
> AuthorDate: Wed Dec 18 13:01:38 2013 +0100
> Commit: Steffen Müthing <steffen.muething at iwr.uni-heidelberg.de>
> CommitDate: Wed May 7 15:29:34 2014 +0200
>
> fixed GMRes
>
> removed bool recalc_defect
> reimplement the update function with respect to optimality
> fixed strange behaviour when maxiter is reached and GMRes wants to restart infinitely many times
> fixed strange behaviour when linear reduction is reached but GMRes wants to restart and continue calculation
> fixed the printing of the results
>
> dune/istl/solvers.hh | 306 ++++++++++++++++++++++-----------------------------
> 1 file changed, 133 insertions(+), 173 deletions(-)
>
>
>
> diff --git a/dune/istl/solvers.hh b/dune/istl/solvers.hh
> index d19e059..8632db2 100644
> --- a/dune/istl/solvers.hh
> +++ b/dune/istl/solvers.hh
> @@ -1046,11 +1046,9 @@ namespace Dune {
> Generalized Minimal Residual method as described the SIAM Templates
> book (http://www.netlib.org/templates/templates.pdf).
>
> - \todo construct F via rebind and an appropriate field_type
> -
> */
>
> - template<class X, class Y=X, class F = Y>
> + template<class X, class Y=X>
> class RestartedGMResSolver : public InverseOperator<X,Y>
> {
> public:
> @@ -1060,29 +1058,23 @@ namespace Dune {
> typedef Y range_type;
> //! \brief The field type of the operator to be inverted
> typedef typename X::field_type field_type;
> - //! \brief The real type of the field type (is the same of using real numbers, but differs for std::complex)
> - typedef typename FieldTraits<field_type>::real_type real_type;
> - //! \brief The field type of the basis vectors
> - typedef F basis_type;
>
> /*!
> \brief Set up solver.
>
> \copydoc LoopSolver::LoopSolver(L&,P&,double,int,int)
> \param restart number of GMRes cycles before restart
> - \param recalc_defect recalculate the defect after everey restart or not [default=false]
> */
> template<class L, class P>
> - RestartedGMResSolver (L& op, P& prec, double reduction, int restart, int maxit, int verbose, bool recalc_defect = false) :
> - _A_(op), _M(prec),
> + RestartedGMResSolver (L& op, P& prec, double reduction, int restart, int maxit, int verbose) :
> + _A(op), _W(prec),
> ssp(), _sp(ssp), _restart(restart),
> - _reduction(reduction), _maxit(maxit), _verbose(verbose),
> - _recalc_defect(recalc_defect)
> + _reduction(reduction), _maxit(maxit), _verbose(verbose)
> {
> dune_static_assert(static_cast<int>(P::category) == static_cast<int>(L::category),
> "P and L must be the same category!");
> - dune_static_assert( static_cast<int>(L::category) == static_cast<int>(SolverCategory::sequential),
> - "L must be sequential!");
> + dune_static_assert(static_cast<int>(L::category) == static_cast<int>(SolverCategory::sequential),
> + "L must be sequential!");
> }
>
> /*!
> @@ -1090,14 +1082,12 @@ namespace Dune {
>
> \copydoc LoopSolver::LoopSolver(L&,S&,P&,double,int,int)
> \param restart number of GMRes cycles before restart
> - \param recalc_defect recalculate the defect after everey restart or not [default=false]
> */
> template<class L, class S, class P>
> - RestartedGMResSolver (L& op, S& sp, P& prec, double reduction, int restart, int maxit, int verbose, bool recalc_defect = false) :
> - _A_(op), _M(prec),
> + RestartedGMResSolver (L& op, S& sp, P& prec, double reduction, int restart, int maxit, int verbose) :
> + _A(op), _W(prec),
> _sp(sp), _restart(restart),
> - _reduction(reduction), _maxit(maxit), _verbose(verbose),
> - _recalc_defect(recalc_defect)
> + _reduction(reduction), _maxit(maxit), _verbose(verbose)
> {
> dune_static_assert(static_cast<int>(P::category) == static_cast<int>(L::category),
> "P and L must have the same category!");
> @@ -1118,218 +1108,189 @@ namespace Dune {
> */
> virtual void apply (X& x, Y& b, double reduction, InverseOperatorResult& res)
> {
> - int m = _restart;
> - real_type norm;
> - real_type norm_old = 0.0;
> - real_type norm_0;
> - real_type beta;
> - int i, j = 1, k;
> + const double EPSILON = 1e-80;
> + const int m = _restart;
> + field_type norm, norm_old = 0.0, norm_0;
> + int j = 1;
> std::vector<field_type> s(m+1), cs(m), sn(m);
> + // need copy of rhs if GMRes has to be restarted
> + Y b2(b);
> // helper vector
> X w(b);
> std::vector< std::vector<field_type> > H(m+1,s);
> - std::vector<F> v(m+1,b);
> + std::vector<X> v(m+1,b);
>
> // start timer
> - Timer watch; // start a timer
> + Dune::Timer watch;
> + watch.reset();
>
> - // clear solver statistics
> + // clear solver statistics and set res.converged to false
> res.clear();
> - _M.pre(x,b);
> - if (_recalc_defect)
> - {
> - // norm_0 = norm(M^-1 b)
> - w = 0.0; _M.apply(w,b); // w = M^-1 b
> - norm_0 = _sp.norm(w); // use defect of preconditioned residual
> - // r = _M.solve(b - A * x);
> - w = b;
> - _A_.applyscaleadd(-1,x, /* => */ w); // w = b - Ax;
> - v[0] = 0.0; _M.apply(v[0],w); // r = M^-1 w
> - beta = _sp.norm(v[0]);
> - }
> - else
> - {
> - // norm_0 = norm(b-Ax)
> - _A_.applyscaleadd(-1,x, /* => */ b); // b = b - Ax;
> - v[0] = 0.0; _M.apply(v[0],b); // r = M^-1 b
> - beta = _sp.norm(v[0]);
> - norm_0 = beta; // use defect of preconditioned residual
> - }
> + _W.pre(x,b);
>
> - // avoid division by zero
> - if (norm_0 == 0.0)
> - norm_0 = 1.0;
> - norm = norm_old = beta;
> + // calculate defect and overwrite rhs with it
> + _A.applyscaleadd(-1.0,x,b); // b -= Ax
> + // calculate preconditioned defect
> + v[0] = 0.0; _W.apply(v[0],b); // r = W^-1 b
> + norm_0 = _sp.norm(v[0]);
> + norm = norm_0;
> + norm_old = norm;
>
> // print header
> - if (_verbose > 0)
> - {
> - std::cout << "=== RestartedGMResSolver" << std::endl;
> - if (_verbose > 1)
> + if(_verbose > 0)
> {
> - this->printHeader(std::cout);
> - this->printOutput(std::cout,0,norm_0);
> - this->printOutput(std::cout,0,norm, norm_0);
> + std::cout << "=== RestartedGMResSolver" << std::endl;
> + if(_verbose > 1) {
> + this->printHeader(std::cout);
> + this->printOutput(std::cout,0,norm_0);
> + }
> }
> - }
>
> - // check convergence
> - if (norm <= reduction * norm_0) {
> - _M.post(x); // postprocess preconditioner
> - res.converged = true;
> - if (_verbose > 0) // final print
> + if(norm_0 < EPSILON) {
> + _W.post(x);
> + res.converged = true;
> + if(_verbose > 0) // final print
> print_result(res);
> - return;
> }
>
> - while (j <= _maxit && res.converged != true) {
> - v[0] *= (1.0 / beta);
> - for (i=1; i<=m; i++) s[i] = 0.0;
> - s[0] = beta;
> + while(j <= _maxit && res.converged != true) {
>
> - for (i = 0; i < m && j <= _maxit && res.converged != true; i++, j++) {
> + int i = 0;
> + v[0] *= 1.0/norm;
> + s[0] = norm;
> + for(i=1; i<m+1; i++)
> + s[i] = 0.0;
> +
> + for(i=0; i < m && j <= _maxit && res.converged != true; i++, j++) {
> w = 0.0;
> - v[i+1] = 0.0; // use v[i+1] as temporary vector
> - _A_.apply(v[i], /* => */ v[i+1]);
> - _M.apply(w, v[i+1]);
> - for (k = 0; k <= i; k++) {
> - H[k][i] = _sp.dot(w, v[k]);
> - // w -= H[k][i] * v[k];
> - w.axpy(-H[k][i], v[k]);
> + // use v[i+1] as temporary vector
> + v[i+1] = 0.0;
> + // do Arnoldi algorithm
> + _A.apply(v[i],v[i+1]);
> + _W.apply(w,v[i+1]);
> + for(int k=0; k<i+1; k++) {
> + H[k][i] = _sp.dot(w,v[k]);
> + // w -= H[k][i] * v[k]
> + w.axpy(-H[k][i],v[k]);
> }
> H[i+1][i] = _sp.norm(w);
> - if (H[i+1][i] == 0.0)
> - DUNE_THROW(ISTLError,"breakdown in GMRes - |w| "
> - << " == 0.0 after " << j << " iterations");
> - // v[i+1] = w * (1.0 / H[i+1][i]);
> - v[i+1] = w; v[i+1] *= (1.0 / H[i+1][i]);
> + if(std::abs(H[i+1][i]) < EPSILON)
> + DUNE_THROW(ISTLError,
> + "breakdown in GMRes - |w| == 0.0 after " << j << " iterations");
> +
> + // normalize new vector
> + v[i+1] = w; v[i+1] *= 1.0/H[i+1][i];
>
> - for (k = 0; k < i; k++)
> - applyPlaneRotation(H[k][i], H[k+1][i], cs[k], sn[k]);
> + // update QR factorization
> + for(int k=0; k<i; k++)
> + applyPlaneRotation(H[k][i],H[k+1][i],cs[k],sn[k]);
>
> - generatePlaneRotation(H[i][i], H[i+1][i], cs[i], sn[i]);
> - applyPlaneRotation(H[i][i], H[i+1][i], cs[i], sn[i]);
> - applyPlaneRotation(s[i], s[i+1], cs[i], sn[i]);
> + // compute new givens rotation
> + generatePlaneRotation(H[i][i],H[i+1][i],cs[i],sn[i]);
> + // finish updating QR factorization
> + applyPlaneRotation(H[i][i],H[i+1][i],cs[i],sn[i]);
> + 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]);
>
> - if (_verbose > 1) // print
> - {
> + // print current iteration statistics
> + if(_verbose > 1) {
> this->printOutput(std::cout,j,norm,norm_old);
> }
>
> norm_old = norm;
>
> - if (norm < reduction * norm_0) {
> + // check convergence
> + if(norm < reduction * norm_0)
> res.converged = true;
> - }
> - }
>
> - if (_recalc_defect)
> - {
> - // update x
> - update(x, i - 1, H, s, v);
> -
> - // update residuum
> - // r = M^-1 (b - A * x);
> - w = b; _A_.applyscaleadd(-1,x, /* => */ w);
> - _M.apply(v[0], w);
> - beta = _sp.norm(v[0]);
> - norm = beta;
> - }
> - else
> - {
> - // calc update vector
> - w = 0;
> - update(w, i - 1, H, s, v);
> -
> - // update x
> - x += w;
> -
> - // r = M^-1 (b - A * x);
> - // update defect
> - _A_.applyscaleadd(-1,w, /* => */ b);
> - // r = M^-1 (b - A * x);
> - v[0] = 0.0; _M.apply(v[0],b); // r = M^-1 b
> - beta = _sp.norm(v[0]);
> - norm = beta;
> -
> - res.converged = false;
> - }
> -
> - //correct i which is wrong if convergence was not achieved.
> - j=std::min(_maxit,j);
> -
> - if (_verbose > 1) // print
> - {
> - this->printOutput(std::cout,j,norm,norm_old);
> + } // end for
> +
> + // calculate update vector
> + w = 0.0;
> + update(w,i,H,s,v);
> + // and current iterate
> + x += w;
> +
> + // restart GMRes if convergence was not achieved,
> + // i.e. linear defect has not reached desired reduction
> + // and if j < _maxit
> + if( res.converged != true && j <= _maxit ) {
> +
> + if(_verbose > 0)
> + std::cout << "=== GMRes::restart" << std::endl;
> + // get saved rhs
> + b = b2;
> + // calculate new defect
> + _A.applyscaleadd(-1.0,x,b); // b -= Ax;
> + // calculate preconditioned defect
> + v[0] = 0.0;
> + _W.apply(v[0],b);
> + norm = _sp.norm(v[0]);
> + norm_old = norm;
> }
>
> - norm_old = norm;
> -
> - if (norm < reduction * norm_0) {
> - // fill statistics
> - res.converged = true;
> - }
> + } //end while
>
> - if (res.converged != true && _verbose > 0)
> - std::cout << "=== GMRes::restart\n";
> - }
> -
> - _M.post(x); // postprocess preconditioner
> + // postprocess preconditioner
> + _W.post(x);
>
> - res.iterations = j;
> - res.reduction = norm / norm_0;
> - res.conv_rate = pow(res.reduction,1.0/j);
> + // save solver statistics
> + res.iterations = j-1; // it has to be j-1!!!
> + res.reduction = norm/norm_0;
> + res.conv_rate = pow(res.reduction,1.0/(j-1));
> res.elapsed = watch.elapsed();
>
> - if (_verbose>0)
> + if(_verbose>0)
> print_result(res);
> +
> }
> - private:
>
> - void
> - print_result (const InverseOperatorResult & res) const
> - {
> - int j = res.iterations>0 ? res.iterations : 1;
> + private :
> +
> + void print_result(const InverseOperatorResult& res) const {
> + int k = res.iterations>0 ? res.iterations : 1;
> std::cout << "=== rate=" << res.conv_rate
> << ", T=" << res.elapsed
> - << ", TIT=" << res.elapsed/j
> + << ", TIT=" << res.elapsed/k
> << ", IT=" << res.iterations
> << std::endl;
> }
>
> - static void
> - update(X &x, int k,
> - std::vector< std::vector<field_type> > & h,
> - std::vector<field_type> & s, std::vector<F> v)
> - {
> + void update(X& w, int i,
> + std::vector<std::vector<field_type> >& H,
> + std::vector<field_type>& s,
> + std::vector<X>& v) {
> + // solution vector of the upper triangular system
> std::vector<field_type> y(s);
>
> - // Backsolve:
> - for (int i = k; i >= 0; i--) {
> - y[i] /= h[i][i];
> - for (int j = i - 1; j >= 0; j--)
> - y[j] -= h[j][i] * y[i];
> - }
> + // backsolve
> + for(int a=i-1; a>=0; a--) {
> + field_type rhs(s[a]);
> + for(int b=a+1; b<i; b++)
> + rhs -= H[a][b]*y[b];
> + y[a] = rhs/H[a][a];
>
> - for (int j = 0; j <= k; j++)
> - // x += v[j] * y[j];
> - x.axpy(y[j],v[j]);
> + // compute update on the fly
> + // w += y[a]*v[a]
> + w.axpy(y[a],v[a]);
> + }
> }
>
> void
> generatePlaneRotation(field_type &dx, field_type &dy, field_type &cs, field_type &sn)
> {
> - if (dy == 0.0) {
> + field_type temp = std::abs(dy);
> + if (std::abs(dy) < 1e-15 ) {
> cs = 1.0;
> sn = 0.0;
> - } else if (std::abs(dy) > std::abs(dx)) {
> - field_type temp = dx / dy;
> + } else if (temp > std::abs(dx)) {
> + temp = dx / dy;
> sn = 1.0 / std::sqrt( 1.0 + temp*temp );
> cs = temp * sn;
> } else {
> - field_type temp = dy / dx;
> + temp = dy / dx;
> cs = 1.0 / std::sqrt( 1.0 + temp*temp );
> sn = temp * cs;
> }
> @@ -1344,15 +1305,14 @@ namespace Dune {
> dx = temp;
> }
>
> - LinearOperator<X,X>& _A_;
> - Preconditioner<X,X>& _M;
> + LinearOperator<X,X>& _A;
> + Preconditioner<X,X>& _W;
> SeqScalarProduct<X> ssp;
> ScalarProduct<X>& _sp;
> int _restart;
> double _reduction;
> int _maxit;
> int _verbose;
> - bool _recalc_defect;
> };
>
>
>
> _______________________________________________
> Dune-Commit mailing list
> Dune-Commit at dune-project.org
> http://lists.dune-project.org/mailman/listinfo/dune-commit
--
Prof. Dr. Christian Engwer
Institut für Numerische und Angewandte Mathematik
Fachbereich Mathematik und Informatik der Universität Münster
Einsteinstrasse 62
48149 Münster
E-Mail christian.engwer at uni-muenster.de
Telefon +49 251 83-35067
FAX +49 251 83-32729
More information about the Dune-devel
mailing list