STK_CG.h

Go to the documentation of this file.

 /*--------------------------------------------------------------------*/
/*     Copyright (C) 2013-2015  Serge Iovleff, Quentin Grimonprez
 
    This program is free software; you can redistribute it and/or modify
    it under the terms of the GNU Lesser General Public License as
    published by the Free Software Foundation; either version 2 of the
    License, or (at your option) any later version.
 
    This program is distributed in the hope that it will be useful,
    but WITHOUT ANY WARRANTY; without even the implied warranty of
    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
    GNU Lesser General Public License for more details.
 
    You should have received a copy of the GNU Lesser General Public
    License along with this program; if not, write to the
    Free Software Foundation, Inc.,
    59 Temple Place,
    Suite 330,
    Boston, MA 02111-1307
    USA
 
    Contact : S..._DOT_I..._AT_stkpp.org (see copyright for ...)
*/
 
/*
 * Project:  stkpp::Algebra
 * created on: 24 mai 2013
 * Author:   Quentin Grimonprez, Serge Iovleff
 **/
 
#ifndef STK_CG_H
#define STK_CG_H
 
#include <Sdk.h>
 
namespace STK
{
 
template<class ColVector>
struct DefaultFunctor
{
    ColVector operator()() const
    { return ColVector();}
};
 
template<class MultFunctor, class ColVector, class InitFunctor = DefaultFunctor<ColVector> >
class CG
{
  public:
    typedef typename ColVector::Type Type;
    CG(): x_(), r_(), eps_(0.), iter_(0), nbStart_(0), p_mult_(0), p_init_(0), p_b_(0) {}
    CG( MultFunctor const& mult
      , ColVector const& b
      , InitFunctor* const& p_init =0
      , Type eps=Arithmetic<Type>::epsilon())
      : x_(), r_()
      , eps_(eps), iter_(0), nbStart_(0)
      , p_mult_(&mult)
      , p_init_(p_init)
      , p_b_(&b)
    {}
 
    CG( CG const& cg)
      : x_(cg.x_)
      , r_(cg.r_)
      , eps_(cg.eps_), iter_(cg.iter_), nbStart_(cg.nbStart_)
      , p_mult_(cg.p_mult_)
      , p_init_(cg.p_init_)
      , p_b_(cg.p_b_)
    {}
    ~CG() {}
    CG* clone() const { return new CG(*this);}
 
    inline ColVector const& x() const { return x_;}
    inline Real const& x(int const& i) const { return x_[i];}
    inline int const& iter() const { return iter_;}
    inline int const& nbStart() const { return nbStart_;}
    inline ColVector const& r() const { return r_;}
    inline void setEps(Type const& eps) {eps_ = eps;}
    inline void setB(ColVector const& b) { p_b_=&b;}
    inline void setMultFunctor(MultFunctor const& mult) { p_mult_= &mult; }
    inline void setInitFunctor(InitFunctor* const& p_init) { p_init_=p_init; }
    int run() { return cg();}
    inline String const& error() const { return msg_error_;}
 
  protected:
    String msg_error_;
    int cg()
    {
      iter_ = 0;
      int nbStart = 0;
      ColVector xOld, z, p_;
 
      // initialization
      if(!p_init_) {x_ = *p_b_;}
      else { x_ = (*p_init_)();}
      while(nbStart<3)
      {
        int step = 0; //number of step
        Real bnorm2 = p_b_->norm2(), alpha, beta; //
        if (bnorm2 == 0.) bnorm2 = 1.;
        //compute the residuals
        r_ = *p_b_ - (*p_mult_)(x_);
        if (r_.norm2()/bnorm2 <eps_) { break;}
        //initialization of the conjugate direction
        p_= r_;
        while(1)
        {
          Real rnorm2 = r_.norm2();
          //compute z=A p
          z.move((*p_mult_)(p_));
          //compute alpha
          alpha = rnorm2/p_.dot(z);
          //update x_
          xOld.exchange(x_);
          x_ = xOld + alpha * p_;
          iter_++;
          //update residuals
          r_ = r_ - (alpha * z);
          //compute beta
          beta = 1/rnorm2;
          if ((rnorm2=r_.norm2())/bnorm2 <eps_) { nbStart = 2; break;}
          beta *= rnorm2;
          //update p_
          p_ = (p_ * beta) + r_;
          step++;
          if( step > 50 ) { break;}
        }
        nbStart++;
      }
      // return an error
      return iter_;
    }
 
  private:
    ColVector x_;
    ColVector r_;
    Type eps_;
    int iter_;
    int nbStart_;
    MultFunctor const*  p_mult_;
    InitFunctor const*  p_init_;
    ColVector const* p_b_;
};
 
template<class MultFunctor, class CondFunctor, class ColVector, class InitFunctor = DefaultFunctor<ColVector> >
class PCG
{
  public:
    typedef typename ColVector::Type Type;
    PCG( MultFunctor const& mult, CondFunctor const& cond, ColVector const& b, InitFunctor* const& p_init =0, Type eps=Arithmetic<Type>::epsilon())
       : x_(), r_()
       , iter_(0), eps_(eps)
       , p_mult_(&mult)
       , p_cond_(&cond)
       , p_init_(p_init)
       , p_b_(&b)
    {};
    PCG( PCG const& pcg)
      : x_(pcg.x_), r_(pcg.r_)
      , eps_(pcg.eps_), iter_(0)
      , p_mult_(pcg.p_mult_)
      , p_cond_(pcg.p_cond_)
      , p_init_(pcg.p_init_)
      , p_b_(pcg.p_b_)
    {};
    ~PCG() {};
    PCG* clone() const { return new PCG(*this);}
 
    inline ColVector const& x() const { return x_;}
    inline Real const& x(int const& i) const { return x_[i];}
    inline ColVector const& r() const { return r_;}
    inline void setEps(Type const& eps) {eps_ = eps;}
    inline void setB(ColVector const& b) { p_b_=&b;}
    inline void setInitFunctor(InitFunctor const& init) { p_init_= &init; }
    inline void setMultFunctor(MultFunctor const& mult) { p_mult_= &mult; }
    inline void setCondFunctor(CondFunctor const& cond) { p_cond_= &cond; }
    inline int run() { return pcg();}
    inline String const& error() const { return msg_error_;}
 
  protected:
    String msg_error_;
    int pcg()
    {
      int nbStart = 0;
      ColVector xOld, y, z, p;
 
      Real bnorm2 = p_b_->norm2(), alpha, beta; //
      iter_= 0;
      // initialization
      if(!p_init_) {x_ = *p_b_;}
      else { x_ = (*p_init_)();}
      // first loop -> allow to restart algorithm in case of divergence.
      while(nbStart<2)
      {
        if (bnorm2 == 0.) bnorm2 = 1.;
        //compute the residuals
        r_ = *p_b_ - (*p_mult_)(x_);
        if (r_.norm2()/bnorm2 <eps_) { break;}
        //initialization of the conjugate direction
        y = (*p_cond_)(r_);
        p = y;
        Real rty=r_.dot(y);
        while(1)
        {
          Real rnorm2 = r_.norm2();
          //compute z=A p
          z.move((*p_mult_)(p));
          //compute alpha
          alpha = rty/p.dot(z);
          //update x_
          xOld.exchange(x_);
          x_ = xOld + alpha * p;
          iter_++;
          //update residuals
          r_ = r_ - (alpha * z);
          //update y
          y = (*p_cond_)(r_);
          //compute beta
          beta = 1/rty;
          if ((rnorm2=r_.norm2())/bnorm2 <eps_) { nbStart = 2; break;}
          rty = r_.dot(y);
          beta *= rty;
          //update p_
          p = (p * beta) + y;
        }
        nbStart++;
      }
      return iter_;
    }
 
  private:
    ColVector x_;
    ColVector r_;
    int iter_;
    Type eps_;
    MultFunctor const*  p_mult_;
    CondFunctor const*  p_cond_;
    InitFunctor const*  p_init_;
    ColVector const* p_b_;
};
 
} // namespace STK
 
#undef MAXITER
 
#endif /* STK_CG_H_ */