chebalg.hh

Go to the documentation of this file.

// chebalg.hh
// created by Yin 01/11/13
// last modified 01/25/13

// Much code is shamelessly stolen from the umin Cheb.H
// originally written by WWS, with his permission

/*************************************************************************

Copyright Rice University, 2004.
All rights reserved.

Permission is hereby granted, free of charge, to any person obtaining a
copy of this software and associated documentation files (the "Software"),
to deal in the Software without restriction, including without limitation
the rights to use, copy, modify, merge, publish, distribute, and/or sell
copies of the Software, and to permit persons to whom the Software is
furnished to do so, provided that the above copyright notice(s) and this
permission notice appear in all copies of the Software and that both the
above copyright notice(s) and this permission notice appear in supporting
documentation.

THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT OF THIRD PARTY
RIGHTS. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR HOLDERS INCLUDED IN THIS
NOTICE BE LIABLE FOR ANY CLAIM, OR ANY SPECIAL INDIRECT OR CONSEQUENTIAL
DAMAGES, OR ANY DAMAGES WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR
PROFITS, WHETHER IN AN ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS
ACTION, ARISING OUT OF OR IN CONNECTION WITH THE USE OR PERFORMANCE OF
THIS SOFTWARE.

Except as contained in this notice, the name of a copyright holder shall
not be used in advertising or otherwise to promote the sale, use or other
dealings in this Software without prior written authorization of the
copyright holder.

**************************************************************************/


#ifndef __RVLALG_UMIN_Cheb_H
#define __RVLALG_UMIN_Cheb_H

#include "alg.hh"
#include "terminator.hh"
#include "linop.hh"
#include "table.hh"

using namespace RVLAlg;

namespace RVLUmin {

  using namespace RVL;
  using namespace RVLAlg;    

  template<typename Scalar>
  class ChebStep : public Algorithm {

    typedef typename ScalarFieldTraits<Scalar>::AbsType atype;

  public:

    ChebStep(LinearOp<Scalar> const & _A,
         Vector<Scalar> & _x,
         Vector<Scalar> const & _b,
         atype & _rnorm, 
         atype & _nrnorm,
         std::vector<atype> &_coeff,
         int & _kc,
         int & _nrt,
         atype _gamma = 0.04,     // inversion level
         atype _epsilon = 0.001,  // error reduction
         atype _alpha = 1.1,      // 'fudge factor'
         atype _lbd_est = 0.0,    // input operator spectrum in case it is know prior
         ostream & _str = cout)
      : A(_A), x(_x), b(_b), rnorm(_rnorm), nrnorm(_nrnorm), coeff(_coeff), gamma(_gamma), epsilon(_epsilon), alpha(_alpha), lbd_est(_lbd_est),kc(_kc), nrt(_nrt), str(_str), dx(A.getDomain()), r(A.getDomain()), ndx(A.getDomain()), tmp_vector(A.getRange()), x_init(_x),dx_init(A.getDomain()),r_init(A.getDomain()),ndx_init(A.getDomain()){
        
      atype one = ScalarFieldTraits<atype>::One();
      atype tmp = one;
      nrt = 0;
      // NOTE: set up dx_init, r_init, ndx_int
      // The real step of initialization is in the run() method
      rnorm = b.norm();
      A.applyAdjOp(b,r_init);
      dx_init.zero();
      nrnorm=r_init.norm();
      ndx_init.zero();
      A.applyOp(r_init,tmp_vector);
      A.applyAdjOp(tmp_vector,ndx_init);
      //nrnorm=ndx_init.norm();
      
      tmp = abs(r_init.inner(ndx_init));
      if (ProtectedDivision<atype>(tmp,r_init.normsq(),RQ)) {
          RVLException e;
          e<<"Error: ChebStep::ChebStep() from ProtectedDivision: RQ\n";
          throw e;
      }
      if(RQ > lbd_est){
          lbd_est = alpha * RQ;
          nrt = nrt + 1;
      }

      ProtectedDivision<atype>(2*one,(one+gamma)*lbd_est,s);
    }
      
    void run() {
      try {
          atype one = ScalarFieldTraits<atype>::One();
          if (kc == 0){
              x.copy(x_init);
              r.copy(r_init);   
              dx.copy(dx_init);
              ndx.copy(ndx_init);
              str<<"Estimated spectrum bound at iter ["<<nrt<<"] = "<< lbd_est << endl;
              str << "NOTE: The  " << nrt << "-th restart\n";
          }
      
          atype tmp;
          dx.linComb(s*coeff[kc+1],r,coeff[kc+1]-one); // this = a*x+b*this
          
          x.linComb(one,dx);
          A.applyOp(dx,tmp_vector);
          A.applyAdjOp(tmp_vector,ndx);
          A.applyOp(x,tmp_vector);
          tmp_vector.linComb(-1*one,b);
          rnorm = tmp_vector.norm();
          
          r.linComb(-1*one,ndx);
          nrnorm=r.norm();
          //nrnorm = ndx.norm();
          tmp = abs(dx.inner(ndx));
          if (ProtectedDivision<atype>(tmp,dx.normsq(),RQ)) {
              RVLException e;
              e<<"Error: ChebStep::run() from ProtectedDivision: RQ\n";
              throw e;
          } 
          if(RQ > lbd_est){
              lbd_est = alpha * RQ;
              ProtectedDivision<atype>(2*one,(one+gamma)*lbd_est,s);
              kc = -1;
              nrt = nrt + 1;
          }
          kc = kc + 1;
      }
      catch (RVLException & e) {
    e<<"\ncalled from ChebStep::run()\n";
    throw e;
      }
     
    }
    
    atype getSpectrumBound() const { return lbd_est; }

    ~ChebStep() {}

  private:

    // references to external objects
    LinearOp<Scalar> const & A;
    Vector<Scalar> & x;
    Vector<Scalar> const & b;
    atype & rnorm;
    atype & nrnorm;
    
    // added for Chebyshev
    std::vector<atype> &coeff;
 
    atype gamma;           // inversion level
    atype epsilon;         // error reduction factor
    atype alpha;           // 'fudge factor'
    atype beta;
    atype RQ;
    atype lbd_est;         // estimated spectrum bound
    atype s;
    int &kc;
    int &nrt;
    ostream & str;
    // end of added
    
    // need four work vectors and one scalar as persistent object data
    Vector<Scalar> dx;
    Vector<Scalar> r;
    Vector<Scalar> ndx;
    Vector<Scalar> tmp_vector;
      
    // need vectors to store initial values for restarting purpose
    Vector<Scalar> x_init;
    Vector<Scalar> dx_init;
    Vector<Scalar> r_init;
    Vector<Scalar> ndx_init;
  };
  
 
  template<typename Scalar>
  class ChebAlg: public Algorithm {

    typedef typename ScalarFieldTraits<Scalar>::AbsType atype;

  public:

    ChebAlg(RVL::Vector<Scalar> & _x, 
        LinearOp<Scalar> const & _inA, 
        Vector<Scalar> const & _rhs, 
        atype & _rnorm,
        atype & _nrnorm,
        atype _gamma = 0.04,     // inversion level
        atype _epsilon = 0.001,  // error reduction
        atype _alpha = 1.1,      // 'fudge factor'
            atype _lbd_est=0.0,
        int _maxcount = 10,      // upper bound of iterations
        ostream & _str = cout)  
      : inA(_inA), x(_x), rhs(_rhs), rnorm(_rnorm), nrnorm(_nrnorm), gamma(_gamma),epsilon(_epsilon),alpha(_alpha),lbd_est(_lbd_est),maxcount(_maxcount), kmax(_maxcount), kc(0), ktot(0), str(_str), step(inA,x,rhs,rnorm,nrnorm,coeff,kc,nrt,gamma,epsilon,alpha,lbd_est,str)
    // NOTE: reference to coeff has been passed to ChebStep object step, however
    // coeff has not been initialized.
    { x.zero();
      // NOTE: estimating the necessary number of iterations
      atype one = ScalarFieldTraits<atype>::One();
      atype beta=0;
      atype epsest=0;
      ProtectedDivision<atype>(sqrt(alpha)*epsilon,one+sqrt(alpha),epsest);
      atype epsk;
      epsk= epsest + one;
      atype tmp = one;
      ProtectedDivision<atype>(one - gamma,one + gamma,beta);
      atype q=0;
      ProtectedDivision<atype>(beta,one + sqrt(one-beta*beta),q);
      int k = 1;
      while (epsk > epsest && k<=kmax) {
    tmp = tmp * q;
    ProtectedDivision<atype>(2*tmp,one + tmp*tmp,epsk);
      k = k + 1;
    }
    kmax = k-1;
    str << "NOTE: computed number of iterations needed:  " << kmax << "\n";
    // NOTE: compute Chebyshev coefficients
    atype ckm = one;
    atype ck = one / beta;
    atype ckp=0;
    coeff.reserve(kmax+1);   // allocate memory for coeff
    coeff[0] = ScalarFieldTraits<atype>::Zero();
    coeff[1] = one;
    for (int i=2; i<kmax+1; i++) {
      ProtectedDivision<atype>(2*ck,beta,ckp);
      ckp = ckp - ckm;
      coeff[i] = one + ckm / ckp;
      ckm = ck;
      ck = ckp;
      //str<<"coeff["<<i<<"] = "<< coeff[i] << endl;
    }
      }
  

    ~ChebAlg() {}

    void run() { 
      
      // terminator for Cheb iteration
      vector<string> names(2);
      vector<atype *> nums(2);
      vector<atype> tols(2);
      atype rnorm0=rnorm;
      atype nrnorm0=nrnorm;
      names[0]="Residual Norm"; nums[0]=&rnorm; tols[0]=ScalarFieldTraits<atype>::Zero();
      names[1]="Normal Residual Norm"; nums[1]=&nrnorm; tols[1]=tols[0]=ScalarFieldTraits<atype>::Zero();
      str<<"========================== BEGIN Cheb =========================\n";
      VectorCountingThresholdIterationTable<atype> stop1(maxcount,names,nums,tols,str);
      MaxTerminator<int> stop2(kc,kmax-1);
      //MaxTerminator<int> stop3(nrt,maxrt);
      stop1.init();
      // terminate if either
      OrTerminator stop(stop1,stop2);
      // loop
      LoopAlg doit(step,stop);
      doit.run();
      ktot = stop1.getCount();

      str<<"=========================== END Cheb ==========================\n";
        // display results
        str<<"\n ******* summary ********  "<<endl;
        str<<"initial residual norm      = "<<rnorm0<<endl;
        str<<"residual norm              = "<<rnorm<<endl;
        str<<"residual redn              = "<<rnorm/rnorm0<<endl;
        str<<"initial gradient norm      = "<<nrnorm0<<endl;
        str<<"gradient norm              = "<<nrnorm<<endl;
        str<<"gradient redn              = "<<nrnorm/nrnorm0<<endl;
    }

    int getCount() const { return kc; }
    int getTotalCount() const { return ktot; }
    int getRestartCount() const { return nrt+1; }
    atype getSpectrumBound() const { return step.getSpectrumBound(); }


  private:

    LinearOp<Scalar> const & inA;  // operator
    Vector<Scalar> & x;            // state - solution vector
    Vector<Scalar> const & rhs;    // reference to rhs
    atype & rnorm;                 // residual norm
    atype & nrnorm;                // gradient norm

    // added for Chebyshev
    std::vector<atype> coeff;
    atype gamma;           // inversion level
    atype epsilon;         // error reduction factor
    atype alpha;           // 'fudge factor'
    atype lbd_est;         // spectrum bound of op
    int maxcount;                  // upper bound for iteration count
    int nrt;                       // upper bound for restarting
    // end of added 

    int kmax;
    int kc;
    int ktot;
    ostream & str;                 // stream for report output
    ChebStep<Scalar> step;         // single step of Cheb

    // disable default, copy constructors
    ChebAlg();
    ChebAlg(ChebAlg<Scalar> const &);

  };

  template<typename Scalar>
  class ChebPolicyData {
    
    typedef typename ScalarFieldTraits<Scalar>::AbsType atype;
    
  public:
    
    int maxcount;
    atype gamma;     // inversion level
    atype epsilon;   // error reduction
    atype alpha;     // 'fudge factor'
    atype lbd_est;   // estimated spectrum bound
    bool verbose;
    
    ChebPolicyData(atype _gamma = 0.04,
           atype _epsilon = 0.01,
           atype _alpha = 1.001,
           atype _lbd_est = 0.0,
           int _maxcount = 0,
           bool _verbose = false)
      : maxcount(_maxcount), gamma(_gamma), epsilon(_epsilon), alpha(_alpha),lbd_est(_lbd_est),verbose(_verbose) {}
    
    ChebPolicyData(ChebPolicyData<Scalar> const & a)
      : maxcount(a.maxcount), gamma(a.gamma), epsilon(a.epsilon), alpha(a.alpha), lbd_est(a.lbd_est), verbose(a.verbose) {}
      
    ostream & write(ostream & str) const {
          str<<"\n";
          str<<"==============================================\n";
          str<<"ChebPolicyData: \n";
          str<<"gamma     = "<<gamma<<"\n";
          str<<"epsilon   = "<<epsilon<<"\n";
          str<<"alpha     = "<<alpha<<"\n";
          str<<"lbd_est   = "<<lbd_est<<"\n";
          str<<"maxcount  = "<<maxcount<<"\n";
          str<<"verbose   = "<<verbose<<"\n";
          str<<"==============================================\n";
          return str;
      }
    
  };
  
  template<typename Scalar> 
  class ChebPolicy {

    typedef typename ScalarFieldTraits<Scalar>::AbsType atype;

  public:
    ChebAlg<Scalar> * build(Vector<Scalar> & x, 
                LinearOp<Scalar> const & A,
                Vector<Scalar> const & d,
                atype & rnorm,
                atype & nrnorm,
                ostream & str) const {
        if (verbose)
      return new ChebAlg<Scalar>(x, A, d,rnorm,nrnorm,gamma,epsilon,alpha,lbd_est,maxcount,str);
        else
      return new ChebAlg<Scalar>(x, A, d,rnorm,nrnorm,gamma,epsilon,alpha,lbd_est,maxcount,nullstr);
    }

    void assign(atype _rtol, atype _nrtol, atype _gamma, atype _epsilon, atype _alpha, atype _lbd_est, int _maxcount, bool _verbose) {
        gamma = _gamma; epsilon=_epsilon; alpha=_alpha; lbd_est=_lbd_est; maxcount=_maxcount; verbose=_verbose;
    }

    void assign(Table const & t) {
      gamma=getValueFromTable<atype>(t,"Cheb_gamma");
      epsilon=getValueFromTable<atype>(t,"Cheb_epsilon");
      alpha=getValueFromTable<atype>(t,"Cheb_alpha");
      lbd_est=getValueFromTable<atype>(t,"Cheb_lbd_est");
      maxcount=getValueFromTable<int>(t,"Cheb_MaxItn");
      verbose=getValueFromTable<bool>(t,"Cheb_Verbose");
    }

    void assign(ChebPolicyData<Scalar> const & s) {
      gamma=s.gamma;
      epsilon=s.epsilon;
      alpha=s.alpha;
      lbd_est=s.lbd_est;
      maxcount=s.maxcount;
      verbose=s.verbose;
    }

    ChebPolicy(atype _gamma = 0.04,
               atype _epsilon = 0.01,
               atype _alpha = 1.001,
               atype _lbd_est = 0.0,
           int _maxcount = 0,
            bool _verbose = true)
      : gamma(_gamma), epsilon(_epsilon), alpha(_alpha), lbd_est(_lbd_est),maxcount(_maxcount), verbose(_verbose), nullstr(0){}
      
    ChebPolicy(ChebPolicy<Scalar> const & p)
      : gamma(p.gamma),
    epsilon(p.epsilon),
    alpha(p.alpha),
    lbd_est(p.lbd_est),
    maxcount(p.maxcount),
    verbose(p.verbose),
    nullstr(0) {}
      
  private:
    mutable atype gamma;
    mutable atype epsilon;
    mutable atype alpha;
    mutable atype lbd_est;
    mutable int maxcount;
    mutable bool verbose;
    mutable std::ostream nullstr;
  };
}

#endif

---