arnoldi.hh

Go to the documentation of this file.

#ifndef __RVL_ALG_LINEAR_SOLVER_H_
#define __RVL_ALG_LINEAR_SOLVER_H_

#include "linop.hh"
#include "alg.hh"

namespace RVLUmin {

  using namespace RVL;

  template<class Scalar>
  class ArnoldiSolverStep : public Algorithm {
  protected: 
    LinearOp<Scalar> & A;
    int k;
    Vector<Scalar> v;
    Vector<Scalar> f;
    CartesianPowerSpace<Scalar> vspc;
    ProductVector<Scalar> V; // length k.  Will push_back new column DC each iter.
    valarray<Scalar> H; // arranged columnwise.  Will add k+1 elements each iter.
    int hcount;
    int iter;

    ArnoldiSolverStep();

  public:
    ArnoldiSolverStep(LinearOp<Scalar> & _A, int num_eigs, const Vector<Scalar> & init_guess)
      :A(_A), k(num_eigs), v(init_guess), f(A.getRange()), 
       vspc(k, A.getDomain()), V(vspc), H(k*(k+1)/2+k), hcount(0), iter(0) {
      Vector<Scalar> w(A.getRange());
      v.scale(Scalar(1.0/v.norm()));
      A.applyOp(v,w); // w <- Av
      Scalar alpha = v.inner(w); // alpha <- v'*Av
      f.copy(w); f.linComb(-alpha, v); // f<- Av - alpha v
      Scalar c = v.inner(f); // c <- v'*(Av-alpha v)
      f.linComb(-c, v); // f<- f-c v
      alpha+=c;
      V[iter].copy(v);
      H[0] = alpha;
      hcount++;
    }

    void run() {
      // note: no-op if you're done - needs coupling to terminator in driver
      if( iter >= k-1 ) return;
      int i = 0;
      Vector<Scalar> w(A.getRange());
      iter++;
      Scalar beta = f.norm();
      H[hcount] = beta; // fill in subdiagonal entry on previous column;
      hcount++;
      v.scale(Scalar(1.0)/beta, f);
      V[iter].copy(v);
      A.applyOp(v,w);
      valarray<Scalar> h(iter+1), c(iter+1);
      // h = V'*w
      for(i=0; i <= iter; i++)
    h[i] = V[i].inner(w);
      // f = w-V*h
      f.copy(w);
      for(i=0; i<= iter; i++)
    f.linComb(-h[i], V[i]);
      // c = V'*f
      for(i=0; i <= iter; i++)
    c[i] = V[i].inner(f);
      // f = f-V*c
      for(i=0; i<= iter; i++)
    f.linComb(-c[i], V[i]);
 
      // fill in upper triangle
      for(i = 0; i < c.size(); i++, hcount++) {
    H[hcount] = c[i]+h[i];
      }
    }

    const valarray<Scalar> & getHessenberg() { return H; }
  
    int getHNumEntries() { return hcount; }

    const ProductVector<Scalar> & getOrthogonalVecs() { return V; }

    int getIterationCount() { return iter; }
    int getMaxIter() { return k; }
  };

  template<typename Scalar>
  class Arnoldi: public Algorithm {



}
#endif

---