arpackpp.hh

Go to the documentation of this file.

 
#ifndef ARPACKPP_HH_
#define ARPACKPP_HH_
 
/* The embedded Arpack includes have virtual overloaded functions
   For this file, we shut this warning
 */
#include "basics/warnings/push.h"
#include "basics/warnings/ignore_warning_overloaded_virtual.h"
 
#include "ARPACK.hh"
//classes for complex (general) eigenvalue problems
#include <eigensolver/arpackpp/arscomp.h>
#include <eigensolver/arpackpp/argcomp.h>
//classes for real symmetric (general) eigenvalue problems
#include <eigensolver/arpackpp/argsym.h>
 
// Restore warning settings
#include "basics/warnings/pop.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) {
      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 = 1, 
                concepts::Cmplx sigma = 0.0, 
                Real tol = 0.0, 
                int maxiter = 300) :
      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 < 1 || kmax > ((int)OP.dimX() - 2))
          throw conceptsException(concepts::ExceptionBase("number of eigenpairs to be computed has to be larger than or equal to one and smaller than or equal to the dimension of the matrix minus 2"));
        ;
      }
    ;
 
    ArPackppStd(concepts::Operator<T>& OP, 
                int kmax = 1, 
                char* which = (char*) "LM",
                Real tol = 0.0, 
                int maxiter = 300) :
      whichEV_(which), 
      modus_(ArPack<Real>::SHIFTINV), 
      aow_(OP), 
      arpackSolver_(
        OP.dimX(), kmax, &aow_, &ArpackStdOperatorWrapper<T>::multOPx, 
        whichEV_, 0, tol, maxiter), 
      computed_(false) {
        if (kmax < 1 || kmax > ((int)OP.dimX() - 2))
          throw conceptsException(concepts::ExceptionBase("number of eigenpairs to be computed has to be larger than or equal to one and smaller than or equal to the dimension of the matrix minus 2"));
        ;
      }
    ;
 
    virtual ~ArPackppStd() {
      for (uint i = 0; i < this->arpackSolver_.EigenvectorsFound(); ++i)
        delete eigenVectors_[i];
    }
 
    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 << concepts::typeOf(*this) << "\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_;
    
    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);
 
    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);
 
//     /** Regular mode constructor
//      @param OP multiplication operator, OP * x = inv(B) * x
//      @param A matrix of the left hand side
//      @param B matrix of the right hand side
//      @param kmax number of eigenpairs to be computed, default 1
//      @param which defines which sort of eigenvalues should be computed, i.e. 
//        eigenvalues of largest magnitude ("LM"), smallest magnitude ("SM"),
//        largest real part ("LR"), smallest real part ("SR"),
//        largest imaginary part ("LI") or smallest imaginary part ("SI"), default 
//        "LM"
//      @param tol convergence tolerance for the eigenpairs, default 0.0 (is 
//        replaced by LAPACK's \c DLAMCH('EPS'))
//      @param maxiter maximum number of Arnoldi iterations, default 300
//      @warning \c A has to be real and symmetric.
//      @warning \c B has to be real, symmetric and positive definite.
//      @warning \c kmax has to be larger than or equal to 1 and smaller than 
//        or equal to the dimension of the matrix minus 2.
//     */
//     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);
 
    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 << concepts::typeOf(*this) << std::endl ;
      switch (modus_) {
        case ArPack<concepts::Real>::REGINV:
          os << "regular inverse mode" << std::endl;
          break;
        case ArPack<concepts::Real>::SHIFTINV:
          os << "shift-invert mode" << std::endl;
          break;
        case ArPack<concepts::Real>::NORMAL:
          os << "normal mode" << std::endl;
          break;
        case ArPack<concepts::Real>::SHIFTINVIMAG:
          os << "shift-invert mode" << std::endl;
          break;
      }
      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.0, 
                int maxiter = 300) :
      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 < 1 || kmax > ((int)OP.dimX() - 2))
          throw conceptsException(concepts::ExceptionBase("number of eigenpairs to be computed has to be larger than or equal to one and smaller than or equal to the dimension of the matrix minus 2"));
        ;
      }
    ;
 
    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) :
      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 < 1 || kmax > ((int)OP.dimX() - 2))
          throw conceptsException(concepts::ExceptionBase("number of eigenpairs to be computed has to be larger than or equal to one and smaller than or equal to the dimension of the matrix minus 2"));
        ;
      }
    ;
 
    virtual ~ArPackppGen() {
      for (uint i = 0; i < this->arpackSolver_.EigenvectorsFound(); ++i)
        delete eigenVectors_[i];
    }
 
    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 << concepts::typeOf(*this) << std::endl;
      switch (modus_) {
        case ArPack<concepts::Real>::REGINV:
          os << "regular inverse mode" << std::endl;
          break;
        case ArPack<concepts::Real>::SHIFTINV:
          os << "shift-invert mode" << std::endl;
          break;
        case ArPack<concepts::Real>::NORMAL:
          os << "normal mode" << std::endl;
          break;
        case ArPack<concepts::Real>::SHIFTINVIMAG:
          os << "shift-invert mode" << std::endl;
          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_ */