[Dune-devel] [Dune-Commit] [Commit] dune-istl - 1f4cd6b: fixed GMRes

Christian Engwer christian.engwer at uni-muenster.de
Thu May 8 18:57:39 CEST 2014


Sorry all others...

the **** reply-to "feature" hit me, in the first place the mail was
intended for Steffen only... As it is out now... the question was why
there are so many white-space and variable-rename changes which make
it really hard to spot the actual changes.

Christian

On Thu, May 08, 2014 at 06:55:05PM +0200, Christian Engwer wrote:
> 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
> 
> _______________________________________________
> Dune-devel mailing list
> Dune-devel at dune-project.org
> http://lists.dune-project.org/mailman/listinfo/dune-devel
> 

-- 
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