eigensolver/arpackpp.hh

Go to the documentation of this file.

#ifndef ARPACKPP_HH_
#define ARPACKPP_HH_

#include "ARPACK.hh"
//classes for complex (general) eigenvalue problems
#include <arpack++/arscomp.h>
#include <arpack++/argcomp.h>
//classes for real symmetric (general) eigenvalue problems
#include <arpack++/argsym.h>
#include "operator.hh"

namespace eigensolver {

  using concepts::Real;
  using concepts::Cmplx;

  template<class T>
  class ArpackStdOperatorWrapper {
  public:

    ArpackStdOperatorWrapper(concepts::Operator<T>& OP) :
      OP_(OP) {
    }

    template<class F>
    void multOPx(F* v, F* w) {
      concepts::Vector<F> vec1(OP_.dimX(), v);
      concepts::Vector<F> vec2(OP_.dimY());
      OP_(vec1, vec2);
      for (uint i = 0; i < OP_.dimY(); ++i)
        w[i] = vec2[i];
    }

  private:
    concepts::Operator<T>& OP_;
  };

  template<class T, class U, class V>
  class ArpackOperatorWrapper {
  public:

    ArpackOperatorWrapper(concepts::Operator<T>& OP, concepts::Operator<U>& A,
        concepts::Operator<V>& B) :
      OP_(OP), A_(A), B_(B) {
    }

    //FIXME replace for-loop with memcpy
    //TODO what happens if dimensions don't agree? Exception??
    template<class S>
    void multAxp(S* v, S* w) {
      concepts::Vector<S> vec1(A_.dimX(), v);
      concepts::Vector<S> vec2(A_.dimY());
      A_(vec1, vec2);
      for (uint i = 0; i < A_.dimY(); ++i)
        w[i] = vec2[i];
    }

    template<class S>
    void multBxp(S* v, S* w) {
      //std::cout << "multBx" << std::endl;
      concepts::Vector<S> vec1(B_.dimX(), v);
      concepts::Vector<S> vec2(B_.dimY());
      B_(vec1, vec2);
      for (uint i = 0; i < B_.dimY(); ++i)
        w[i] = vec2[i];
    }

    template<class S>
    void multOPx(S* v, S* w) {
      concepts::Vector<S> vec1(OP_.dimX(), v);
      concepts::Vector<S> vec2(OP_.dimY());
      OP_(vec1, vec2);
      for (uint i = 0; i < OP_.dimY(); ++i)
        w[i] = vec2[i];
    }

    template<class S>
    void multOPxRegular(S* v, S* w) {

      concepts::Vector<S> vec1(OP_.dimX(), v);
      concepts::Vector<S> vec2(A_.dimY());
      //vec2 = A*v
      A_(vec1, vec2);
      //vec1 = OP * vec2 = OP * A * v
      OP_(vec2, vec1);
      for (uint i = 0; i < OP_.dimY(); ++i)
        w[i] = vec1[i];
    }

  private:
    concepts::Operator<T>& OP_;
    concepts::Operator<U>& A_;
    concepts::Operator<V>& B_;
  };

  template<class T>
  class ArPackppStd: public EigenSolver<concepts::Cmplx> {
  public:

    ArPackppStd(concepts::Operator<T>& OP, int kmax, concepts::Cmplx sigma,
        Real tol = 0.0, int maxiter = 300,
        concepts::Vector<concepts::Cmplx>* start = 0, concepts::Array<
            concepts::Cmplx>* resid = 0, bool schur = false) :
      whichEV_((char *) "LM"), modus_(ArPack<Real>::SHIFTINV), aow_(OP),
          arpackSolver_(OP.dimX(), kmax, &aow_,
              &ArpackStdOperatorWrapper<T>::multOPx, sigma, whichEV_, 0, tol,
              maxiter), computed_(false) {
      if (kmax <= 0 || kmax > (((int)OP.dimX()) - 1))
        throw conceptsException(concepts::ExceptionBase("Dimension must be greater than zero and smaller then the dimension of the problem -1"));
      ;
    }
    ;

    ArPackppStd(concepts::Operator<T>& OP, int kmax = 1, char* which = "LM",
        Real tol = 0.0, int maxiter = 300,
        concepts::Vector<concepts::Cmplx>* start = 0, concepts::Array<
            concepts::Cmplx>* resid = 0, bool schur = false) :
      whichEV_(which), modus_(ArPack<Real>::SHIFTINV), aow_(OP), arpackSolver_(
          OP.dimX(), kmax, &aow_, &ArpackStdOperatorWrapper<T>::multOPx,
          whichEV_, 0, tol, maxiter), computed_(false) {
      if (kmax <= 0 || kmax > (((int)OP.dimX()) - 1))
        throw conceptsException(concepts::ExceptionBase("Dimension must be greater than zero and smaller then the dimension of the problem -1"));
      ;
    }
    ;

    virtual ~ArPackppStd() {
      for (uint i = 0; i < this->arpackSolver_.EigenvectorsFound(); ++i)
        delete eigenVectors_[i];
    }

    /*
     * Getter for the computed eigenvalues. Computes Eigenvalus and Eigenvectors if
     * they were not computed.
     */
    virtual const concepts::Array<Cmplx>&
    getEV() {
      if (!computed_)
        compute_();
      return eigenValues_;
    }

    virtual concepts::Array<concepts::Vector<Cmplx>*>&
    getEF() {
      if (!computed_)
        compute_();
      return eigenVectors_;
    }

    virtual uint iterations() const {
      return iter_;
    }

    virtual uint converged() const {
      return convEigenvalues_;
    }

    inline concepts::Array<Cmplx> getRESID() {
      concepts::Array<Cmplx> res(this->arpackSolver_.RawResidualVector(),
          this->arpackSolver_.GetN());
      return res;
    }
  protected:
    virtual std::ostream&
    info(std::ostream& os) const {
      os << "ArPackppSymmGen \n";
      switch (modus_) {
        case ArPack<concepts::Real>::REGINV:
          os << "regular inverse mode \n";
          break;
        case ArPack<concepts::Real>::SHIFTINV:
          os << "shift-invert mode \n";
          break;
        case ArPack<concepts::Real>::NORMAL:
          os << "normal mode \n";
          break;
        case ArPack<concepts::Real>::SHIFTINVIMAG:
          os << "shift-invert mode \n";
          break;
      }
      return os;
    }
  private:
    void compute_() {
      int numOfVecs = this->arpackSolver_.FindEigenvectors();
      //std::cout << numOfVecs << std::endl;
      eigenValues_.resize(this->arpackSolver_.ConvergedEigenvalues());

      for (int i = 0; i < this->arpackSolver_.ConvergedEigenvalues(); ++i) {
        eigenValues_[i] = this->arpackSolver_.Eigenvalue(i);
      }
      eigenVectors_.resize(this->arpackSolver_.GetN());
      for (int i = 0; i < numOfVecs; ++i)
        eigenVectors_[i] = new concepts::Vector<concepts::Cmplx>(
            this->arpackSolver_.GetN(), this->arpackSolver_.RawEigenvectors()
                + i * this->arpackSolver_.GetN());

      //for (int i = 0; i < 8; ++i)
      //std::cout << "AHA: " << this->arpackSolver_.RawEigenvectors()[i]
      //    << std::endl;

      computed_ = true;
      iter_ = arpackSolver_.GetIter();
      //std::cout << "Num of iterations: " << iter_ << std::endl;
      convEigenvalues_ = arpackSolver_.ConvergedEigenvalues();
    }
    ;
    char* whichEV_;
    ArPack<Real>::modus modus_;
    uint n_;
    ArpackStdOperatorWrapper<T> aow_;
    ARCompStdEig<Real, ArpackStdOperatorWrapper<T> > arpackSolver_;

    uint iter_;
    uint convEigenvalues_;
    bool computed_;

    concepts::Array<Cmplx> eigenValues_;
    concepts::Array<concepts::Vector<Cmplx>*> eigenVectors_;

  };

  class ArPackppSymGen: public EigenSolver<Real> {

  public:

    ArPackppSymGen(concepts::Operator<Real>& OP, concepts::Operator<Real>& A,
        concepts::Operator<Real>& B, int kmax = 1, Real sigma = 0.1, Real tol =
            0.0, int maxiter = 300, concepts::Vector<Real>* start = 0,
        concepts::Array<Real>* resid = 0, bool schur = false);

    ArPackppSymGen(char invertType, concepts::Operator<Real>& OP,
        concepts::Operator<Real>& A, int kmax = 1, Real sigma = 0.1, Real tol =
            0.0, int maxiter = 300, concepts::Vector<Real>* start = 0,
        concepts::Array<Real>* resid = 0, bool schur = false);

    //    /** Constructor for SHIFTandINVERT mode and Buckling mode
    //     @warning B have to be a symmetric positive definite matrix
    //     @param OP is the multiplication operator B^-1
    //     @param A Stiffness matrix
    //     @param B is the Matrix on the right hand side
    //     @param kmax Number of eigenpairs to be computed,
    //     must be lower then the dimension of the matrices -1
    //     @warning kmax have to be lower than the dimension of the matrix -1!
    //     @param tol Convergence tolerance for the eigenpairs. The default
    //     value 0.0 is replaced by \c DLAMCH('EPS') from LAPACK.
    //     @param maxiter Maximum number of Arnoldi iterations allowed
    //     @param target What sort of eigenvalues to compute
    //     @param sigma Shift for the shift-invert modes
    //     @warning In the buckling mode sigma should be a real number that is NOT ZERO!
    //     @param start Initial vector for iteration
    //     @param schur Calculate Schur vector basis instead of eigenvectors
    //     */
    //    ArPackppSymGen(concepts::Operator<Real>& OP, concepts::Operator<Real>& A,
    //        concepts::Operator<Real>& B, int kmax = 1, char* which = (char*) "LM",
    //        Real tol = 0.0, int maxiter = 300, concepts::Vector<Real>* start = 0,
    //        concepts::Array<Real>* resid = 0, bool schur = false);

    virtual ~ArPackppSymGen() {
      for (uint i = 0; i < this->arpackSolver_.EigenvectorsFound(); ++i)
        delete eigenVectors_[i];
    }

    virtual const concepts::Array<Real>&
    getEV() {
      if (!computed_)
        compute_();
      return eigenValues_;
    }

    virtual concepts::Array<concepts::Vector<Real>*>&
    getEF() {
      if (!computed_)
        compute_();
      return eigenVectors_;
    }

    virtual uint iterations() const {
      return iter_;
    }

    virtual uint converged() const {
      return convEigenvalues_;
    }

    inline concepts::Array<Real> getRESID() {
      concepts::Array<Real> res(this->arpackSolver_.RawResidualVector(),
          this->arpackSolver_.GetN());
      return res;
    }

  protected:
    virtual std::ostream& info(std::ostream& os) const {
      os << "ArPackppSymmGen \n";
      return os;
    }

  private:
    void compute_();
    char* whichEV_;
    ArPack<Real>::modus modus_;
    ArpackOperatorWrapper<concepts::Real, concepts::Real, concepts::Real> aow_;
    ARSymGenEig<Real, ArpackOperatorWrapper<concepts::Real, concepts::Real,
        concepts::Real> , ArpackOperatorWrapper<concepts::Real, concepts::Real,
        concepts::Real> > arpackSolver_;

    uint n_;
    uint iter_;
    uint convEigenvalues_;
    bool computed_;

    concepts::Array<Real> eigenValues_;
    concepts::Array<concepts::Vector<Real>*> eigenVectors_;

  };
  template<class F, class G, class H = concepts::Real>
  class ArPackppGen: public EigenSolver<concepts::Cmplx> {
  public:

    ArPackppGen(concepts::Operator<F>& OP, concepts::Operator<G>& A,
        concepts::Operator<H>& B, int kmax = 1, concepts::Cmplx sigma = 0.0,
        Real tol = 0.01, int maxiter = 300,
        concepts::Vector<concepts::Cmplx>* start = 0, concepts::Array<
            concepts::Cmplx>* resid = 0, bool schur = false) :
      whichEV_((char*) "LM"), modus_(ArPack<Real>::SHIFTINV), aow_(OP, A, B),
          arpackSolver_(A.dimX(), kmax, &aow_,
              &ArpackOperatorWrapper<F, G, H>::multOPx, &aow_,
              &ArpackOperatorWrapper<F, G, H>::multBxp, sigma, whichEV_, 0,
              tol, maxiter), computed_(false) {
      if (kmax <= 0 || kmax > (((int)OP.dimX()) - 1))
        throw conceptsException(concepts::ExceptionBase("Dimension must be greater than zero and smaller then the dimension of the problem -1"));
      ;
    }
    ;

    ArPackppGen(concepts::Operator<F>& OP, concepts::Operator<G>& A,
        concepts::Operator<H>& B, int kmax = 1, char* which = (char*) "LM",
        Real tol = 0.0, int maxiter = 300,
        concepts::Vector<concepts::Cmplx>* start = 0, concepts::Array<
            concepts::Cmplx>* resid = 0, bool schur = false) :
      whichEV_(which), modus_(ArPack<Real>::SHIFTINV), aow_(OP, A, B),
          arpackSolver_(A.dimX(), kmax, &aow_,
              &ArpackOperatorWrapper<F, G, H>::multOPxRegular, &aow_,
              &ArpackOperatorWrapper<F, G, H>::multBxp, whichEV_, 0, tol,
              maxiter), computed_(false) {
      if (kmax <= 0 || kmax > (((int)OP.dimX()) - 1))
        throw conceptsException(concepts::ExceptionBase("Dimension must be greater than zero and smaller then the dimension of the problem -1"));
      ;
    }
    ;

    virtual ~ArPackppGen() {
      for (uint i = 0; i < this->arpackSolver_.EigenvectorsFound(); ++i)
        delete eigenVectors_[i];
    }

    /*
     * Getter for the computed eigenvalues. Computes Eigenvalus and Eigenvectors if
     * they were not computed.
     */
    virtual const concepts::Array<Cmplx>&
    getEV() {
      if (!computed_)
        compute_();
      return eigenValues_;
    }

    virtual concepts::Array<concepts::Vector<Cmplx>*>&
    getEF() {
      if (!computed_)
        compute_();
      return eigenVectors_;
    }

    virtual uint iterations() const {
      return iter_;
    }

    virtual uint converged() const {
      return convEigenvalues_;
    }

    inline concepts::Array<Cmplx> getRESID() {
      concepts::Array<Cmplx> res(this->arpackSolver_.RawResidualVector(),
          this->arpackSolver_.GetN());
      return res;
    }
  protected:
    virtual std::ostream&
    info(std::ostream& os) const {
      os << "ArPackppSymmGen \n";
      switch (modus_) {
        case ArPack<concepts::Real>::REGINV:
          os << "regular inverse mode \n";
          break;
        case ArPack<concepts::Real>::SHIFTINV:
          os << "shift-invert mode \n";
          break;
        case ArPack<concepts::Real>::NORMAL:
          os << "normal mode \n";
          break;
        case ArPack<concepts::Real>::SHIFTINVIMAG:
          os << "shift-invert mode \n";
          break;
      }
      return os;
    }
  private:
    void compute_() {
      int numOfVecs = this->arpackSolver_.FindEigenvectors();
      eigenValues_.resize(this->arpackSolver_.ConvergedEigenvalues());

      for (int i = 0; i < this->arpackSolver_.ConvergedEigenvalues(); ++i) {
        eigenValues_[i] = this->arpackSolver_.Eigenvalue(i);
      }
      eigenVectors_.resize(this->arpackSolver_.GetN());
      for (int i = 0; i < numOfVecs; ++i)
        eigenVectors_[i] = new concepts::Vector<concepts::Cmplx>(
            this->arpackSolver_.GetN(), this->arpackSolver_.RawEigenvectors()
                + i * this->arpackSolver_.GetN());

      computed_ = true;
      iter_ = arpackSolver_.GetIter();
      convEigenvalues_ = arpackSolver_.ConvergedEigenvalues();
    }
    ;
    char* whichEV_;
    ArPack<Real>::modus modus_;
    uint n_;
    ArpackOperatorWrapper<F, G, H> aow_;
    ARCompGenEig<Real, ArpackOperatorWrapper<F, G, H> , ArpackOperatorWrapper<
        F, G, H> > arpackSolver_;

    uint iter_;
    uint convEigenvalues_;
    bool computed_;

    concepts::Array<Cmplx> eigenValues_;
    concepts::Array<concepts::Vector<Cmplx>*> eigenVectors_;

  };

} // namespace eigensolver

#endif /* ARPACKPP_HH_ */