Using Arpack#

The arpackppSupport library contains the interface of XLiFE++ with Arpack library through Arpack++. Because of commonality and stability, Arpack library is chosen as the external eigen solver.

Arpack is a well known collection of FORTRAN subroutines designed to compute a few eigenvalues and eigenvectors of large scale sparse matrices and pencils. It implements a variant of the Arnoldi process for finding eigenvalues called Implicit restarted Arnoldi method (IRAM) and is capable of solving a great variety of problems from single precision positive definite symmetric problems to double precision complex non-Hermitian generalized problems.

Arpack++ is the C++ an interface of Arpack.

XLiFE++ offers an interface to Arpack++ through a well-organized collection of classes which hide all the complexity of Arpack++. User can easily call an eigen solver like a call to an iterative solver. Moreover, this interface offers advanced users a flexible way to make use of all features of Arpack++ with a minimum effort by modifying some classes.

This package is organized into 2 layers:

Users Interface wraps all Arpack++ interface and provides simple call prototype to user
Arpack++ Interface is collection of classes interacting directly with Arpack++

In order to use Arpack++, Arpack must be already installed.

Arpack++ package can be found at Arpack++. However, because of the deprecation of this version, a patch needs to be applied to make sure a correct compilation.

The following is organized in “bottom-up”, the Arpack++ interface is described then Users Interface.

Arpack++ interface layer#

Arpack++ interface is a group of class which is an interface between XLiFE++ and Arpack++. Each class of this group corresponds to a specific eigen problem with a computational mode.

Figure made with TikZ

LargeMatrixArpackpp is an interface class between XLiFE++ and Arpack++. This class is an abstract base class so that all classes to be used in the program are derived from this one. This class defines some basic properties and methods which are used in all derived classes.

template<typename K, class Mat>
class LargeMatrixArpackpp
{
  public:
   LargeMatrixArpackpp(const Mat& matA, const Number sizeProblem);
   LargeMatrixArpackpp(const Mat& matA, const Number sizeProblem, const String name);
   virtual ~LargeMatrixArpackpp();
   virtual void multOPx (K* x, K* y) = 0;
   String kindOfFactorization() const {return factName_;}

  protected:
   const Mat* matA_p;     //< Left hand-side matrix A, generally stiffness matrix (given data)
   Number size_;    //< Size of problem
   String name_;     //< Name of problem
   String factName_;     //< Name of the algorithm used to compute the factorization pointed to by fact_p
   Mat* fact_p;  //<Factorization of matB, matAsI or matAsB in order to speed up linear  systems solving.

  protected:
   void (Mat::*solver_)(std::vector<K>& b, std::vector<K>& x);
   void checkDims(const Mat* mat);
   void factorize(const Mat* mat);
   void bzMultAx(K* x, K* y);     //< Matrix-vector product y <- A * x
   void array2Vector(const K*, std::vector<K>&);     //< Convert array to Vector
   void vector2Array(std::vector<K>&, K*);     //< Convert Vector to Array
}; // end of Class LargeMatrixArpackpp

See also

library=eigenSolvers,
header=
implementation=,
test=test_EigenSolverArpackpp.cpp,
header dep=

The `LargeMatrixArpackppStdReg` class#

LargeMatrixArpackppStdReg class is an interface between XLiFE++ and Arpack++ designed to enable the use of Arpack++’s own classes for standard problems in regular mode. This class allows XLiFE++ users to solve real symmetric, non-symmetric and complex standard eigen problem in regular mode. This class holds informations of matrices and some methods related to product matrix-vector. In general, only matrix-vector product \(y\leftarrow Ax\) is needed.

template<typename K, class Mat>
class LargeMatrixArpackppStdReg: public LargeMatrixArpackpp<K, Mat>
{
  public:
    LargeMatrixArpackppStdReg(const Mat& matA, const Number sizeProblem);
    virtual ~LargeMatrixArpackppStdReg() {}
    virtual void multOPx(K* x, K* y);     //< Matrix-vector product y <- A * x

  private:
    //! Duplicating such object is disabled
    LargeMatrixArpackppStdReg(const LargeMatrixArpackppStdReg&);
    LargeMatrixArpackppStdReg& operator=(const LargeMatrixArpackppStdReg&);
}; // end of Class LargeMatrixArpackppStdReg

See also

library=eigenSolvers,
header=
implementation=,
test=test_EigenSolverArpackpp.cpp,
header dep=LargeMatrixArpackpp.hpp

The `LargeMatrixArpackppStdShf` class#

LargeMatrixArpackppStdShf class is an interface between XLiFE++ and Arpack++ designed to enable the use of Arpack++’s own classes for standard problems in shift and inverse mode. This class allows XLiFE++ users to solve real symmetric, non-symmetric and complex standard eigen problem in shift and inverse mode. This class holds informations of matrices and some methods related to product matrix-vector. In general, we need to provide matrix-vector product \(y \leftarrow (A-\sigma I)^{-1}x\). However, thanks to some built-in special solvers of XLiFE++, only \((A-\sigma I)\) is enough.

template<typename K, class Mat>
class LargeMatrixArpackppStdShf: public LargeMatrixArpackpp<K, Mat>
{
  public:
    LargeMatrixArpackppStdShf(const Mat& matA, const Number sizeProblem, const K sigma);
    virtual ~LargeMatrixArpackppStdShf();
    virtual void multOPx (K* x, K* y);     //< Matrix-vector product y <- inv(A-sigma*Id) * x
  private:
    const Mat* matAsI_p; //<         A - sigma Id matrix for shift and invert mode in standard problems
    bool localAlloc_;     //< flag to indicate if dynamic memory is allocated to store the matrix (A - sigma*Id)

    //! Duplicating such object is disabled
    LargeMatrixArpackppStdShf(const LargeMatrixArpackppStdShf&);
    LargeMatrixArpackppStdShf& operator=(const LargeMatrixArpackppStdShf&);
}; // end of Class LargeMatrixArpackppStdShf

See also

library=eigenSolvers,
header=
implementation=,
test=test_EigenSolverArpackpp.cpp,
header dep=LargeMatrixArpackpp.hpp

The `LargeMatrixArpackppGen` class#

LargeMatrixArpackppGen is an abstract class which plays a role of an interface between XLiFE++ and Arpack++ designed to enable the use of Arpack++’s own classes for generalized eigen problem \(Ax=\lambda Bx\). All classes for generalized problem are derived from this one.

template<typename K, class Mat>
class LargeMatrixArpackppGen: public LargeMatrixArpackpp<K, Mat>
{
  public:
    LargeMatrixArpackppGen(const Mat& matA, const Mat& matB, const Number sizeProblem);
    LargeMatrixArpackppGen(const Mat& matA, const Mat& matB, const Number sizeProblem, const String name);
    virtual ~LargeMatrixArpackppGen() {}
    virtual void multBx (K* x, K* y); //<   Matrix-vector product y <- B * x required for all derived classes

  protected:
    const Mat* matB_p; //< Right hand-side matrix B, e.g. mass matrix (given data)
    Mat* matAsB_p; //!< Matrix (A - sigma*B)
    bool localAlloc_; //<     flag to indicate if dynamic memory is allocated to store A - sigma B,

  private:
    //! Duplicating such object is disabled
    LargeMatrixArpackppGen(const LargeMatrixArpackppGen&);
    LargeMatrixArpackppGen& operator=(const LargeMatrixArpackppGen&);
}; // end of Class LargeMatrixArpackppGen

See also

library=eigenSolvers,
header=
implementation=,
test=test_EigenSolverArpackpp.cpp,
header dep=LargeMatrixArpackpp.hpp

The `LargeMatrixArpackppGenReg` class#

LargeMatrixArpackppGenReg class is interface classes between XLiFE++ and Arpack++ designed to enable the use of Arpack++’s own classes for generalized problems in regular mode except for real symmtric generalize problem. This class holds informations of matrices and some methods related to product matrix-vector.

In general, we need to provide matrix-vector product \(w\leftarrow (B)^{-1}Ax\) and \(z\leftarrow Bx\)

template<typename K, class Mat>
class LargeMatrixArpackppGenReg: public LargeMatrixArpackppGen<K, Mat>
{
  public:
   LargeMatrixArpackppGenReg(const Mat& matA, const Mat& matB, const Number sizeProblem);
   virtual ~LargeMatrixArpackppGenReg() {};
   virtual void multOPx (K* x, K* y);  //< Matrix-vector product y <- inv(B)*A * x

  private:
   //! Duplicating such object is disabled
   LargeMatrixArpackppGenReg(const LargeMatrixArpackppGenReg&);
   LargeMatrixArpackppGenReg& operator=(const LargeMatrixArpackppGenReg&);
}; // end of Class LargeMatrixArpackppGenReg

See also

library=eigenSolvers,
header=
implementation=,
test=test_EigenSolverArpackpp.cpp,
header dep=LargeMatrixArpackppGen.hpp

The `LargeMatrixArpackppSymGenReg` class#

LargeMatrixArpackppGenReg is basically not different from class LargeMatrixArpackppGenReg except on how the calculation results are returned. Maybe in the future, we don’t need this class anymore.

template<typename K, class Mat>
class LargeMatrixArpackppSymGenReg: public LargeMatrixArpackppGen<K, Mat>
{
  public:
   LargeMatrixArpackppSymGenReg(const Mat& matA, const Mat& matB, const Number sizeProblem);
   virtual ~LargeMatrixArpackppSymGenReg() {}
   virtual void multOPx (K* x, K* y);  //< Matrix-vector product y <- inv(B)*A * x

  private:
   //! Duplicating such object is disabled
   LargeMatrixArpackppSymGenReg(const LargeMatrixArpackppSymGenReg&);
   LargeMatrixArpackppSymGenReg& operator=(const LargeMatrixArpackppSymGenReg&);
}; // end of Class LargeMatrixArpackppSymGenReg

See also

library=eigenSolvers,
header=
implementation=,
test=test_EigenSolverArpackpp.cpp,
header dep=LargeMatrixArpackppGen.hpp

The `LargeMatrixArpackppGenShf` class#

Class LargeMatrixArpackppGenShf is interface between XLiFE++ and Arpack++ designed to enable the use of Arpack++’s own classes for generalized problems in shift and invert mode. All classes for generalized problem in relevant shift and invert mode (including Buckling and Cayley) are derived from this one.

template<typename K, class Mat>
class LargeMatrixArpackppGenShf: public LargeMatrixArpackppGen<K, Mat>
{
  public:
   LargeMatrixArpackppGenShf(const Mat& matA, Mat& matB, const Number sizeProblem, const K sigma, const String name = "Generalized Shift and Invert");
   virtual ~LargeMatrixArpackppGenShf() {}
   virtual void multOPx (K* x, K* y); //< Matrix-vector product y <- inv(A-sigma B)*x

  private:
   //! Duplicating such object is disabled
   LargeMatrixArpackppGenShf(const LargeMatrixArpackppGenShf&);
   LargeMatrixArpackppGenShf& operator=(const LargeMatrixArpackppGenShf&);
}; // end of Class LargeMatrixArpackppGenShf

See also

library=eigenSolvers,
header=
implementation=,
test=test_EigenSolverArpackpp.cpp,
header dep=LargeMatrixArpackppGen.hpp

The `LargeMatrixArpackppGenBuc` class#

Class LargeMatrixArpackppGenBuc is derived from class LargeMatrixArpackppGenShf to solve generalized problems Buckling mode. Two products matrix-vector need provided: \(w\leftarrow (A-\sigma B)^{-1}x\) and \(z\leftarrow Ax\)

template<typename K, class Mat>
class LargeMatrixArpackppGenBuc: public LargeMatrixArpackppGenShf<K, Mat>
{
 public:
   LargeMatrixArpackppGenBuc(const Mat& matA, Mat& matB, const Number sizeProblem, const K sigma);
   virtual ~LargeMatrixArpackppGenBuc() {}
   virtual void multBx (K* x, K* y);//<  Matrix-vector product y <- A * x

 private:
   //! Duplicating such object is disabled
   LargeMatrixArpackppGenBuc(const LargeMatrixArpackppGenBuc&);
   LargeMatrixArpackppGenBuc& operator=(const LargeMatrixArpackppGenBuc&);
}; // end of Class LargeMatrixArpackppGenBuc

See also

library=eigenSolvers,
header=
implementation=,
test=test_EigenSolverArpackpp.cpp,
header dep=LargeMatrixArpackppGen.hpp

The `LargeMatrixArpackppGenCay` class#

Class LargeMatrixArpackppGenCay is derived from class LargeMatrixArpackppGenShf to solve generalized problems in Cayley mode. Three products matrix-vector need provided: \(y\leftarrow (A-\sigma B)^{-1}z\), \(w\leftarrow Az\) and \(u\leftarrow Bz\)

template<typename K, class Mat>
class LargeMatrixArpackppGenBuc: public LargeMatrixArpackppGenShf<K, Mat>
{
 public:
   LargeMatrixArpackppGenBuc(const Mat& matA, Mat& matB, const Number sizeProblem, const K sigma);
   virtual ~LargeMatrixArpackppGenBuc() {}
   virtual void multBx (K* x, K* y);//<  Matrix-vector product y <- A * x

 private:
   //! Duplicating such object is disabled
   LargeMatrixArpackppGenBuc(const LargeMatrixArpackppGenBuc&);
   LargeMatrixArpackppGenBuc& operator=(const LargeMatrixArpackppGenBuc&);
}; // end of Class LargeMatrixArpackppGenBuc

See also

library=eigenSolvers,
header=
implementation=,
test=test_EigenSolverArpackpp.cpp,
header dep=LargeMatrixArpackppGen.hpp

The `LargeMatrixArpackppGenCsh` class#

Class LargeMatrixArpackppGenCsh is derived from class LargeMatrixArpackppGen to solve real non-symmetric generalized problem in complex shift and invert mode. Three products matrix-vector need provided: \(y\leftarrow (A-\sigma B)^{-1}z\), \(w\leftarrow Az\) and \(u\leftarrow Bz\)

template<typename K, class Mat>
class LargeMatrixArpackppGenCsh : public LargeMatrixArpackppGen<K, Mat>
{
 public:
  LargeMatrixArpackppGenCsh(const Mat& matA, Mat& matB, const Number sizeProblem, Complex sigma, const char part);
  virtual ~LargeMatrixArpackppGenCsh() {};
  /*!
      Matrix-vector product y <- real(inv(A - sigma*B)) if part = 'R'
                            y <- imag(inv(A - sigma*B)) if part = 'I'
  */
  void virtual multOPx (K* x, K* y);
  void multAx (K* x, K* y); //< Matrix-vector product y <- A * x
 private:
  //! Duplicating such object is disabled
  LargeMatrixArpackppGenCsh(const LargeMatrixArpackppGenCsh&);
  LargeMatrixArpackppGenCsh& operator=(const LargeMatrixArpackppGenCsh&);
}; // end of Class LargeMatrixArpackppGenCsh

See also

library=eigenSolvers,
header=
implementation=,
test=test_EigenSolverArpackpp.cpp,
header dep=LargeMatrixArpackppGen.hpp

Users interface layer#

User interface contains classes which wrap all Arpack++ interface layer and provide user a simple prototype. All eigen solvers are invoked by operator(), which is similar to other iterative solvers of solvers library. There are two principal classes, each of which solves a specific eigen problem

EigenSolverStdArpackpp for standard eigen problem \(Ax=\lambda x\)
EigenSolverGenArpackpp for generalized eigen problem \(Ax=\lambda Bx\)

For each type of eigen problem, input matrix \(A\) decides the value type (Real or Complex) of eigenvalue and/or eigenvector found. There are three cases corresponding to input matrix \(A\): real symmetric, complex and real non-symmetric. For the first two cases, the eigenvalue and/or eigenvector found have the same value type of matrix \(A\) (i.e: Real if matrix \(A\) is real symmetric and Complex if matrix \(A\) is complex). For the real non-symmetric case, eigenvalue and/or eigenvector can be complex. Based on value type of input matrix and output eigenvalue (and/or eigenvector), an operator () can be invoked to solve an eigen problem correctly. There are three prototypes of operator() corresponding to three cases:

real symmetric case: input matrix \(A\) is real symmetric and output eigenvalues/eigenvectors are real
real non-symmetric case: input matrix \(A\) is real non-symmetric and output eigenvalues/eigenvectors are complex
complex case : input matrix \(A\) is complex and output eigenvalues/eigenvectors are complex

The `EigenSolverBaseArpackpp` class#

The EigenSolverBaseArpackpp is the base class which provides some methods to construct object and retrieve eigenvalue/eigenvector returned by Arpack++. All other wrapper classes are derived from this one.

class EigenSolverBaseArpackpp
{
 public:
  EigenSolverBaseArpackpp(const String& nam);
  EigenSolverBaseArpackpp(const String& nam, const Number maxOfIt, const Real epsilon);
  virtual ~EigenSolverBaseArpackpp();
  // Iterative solver name for documentation purposes
  String name() { return name_; }
  // Some functions to extract Eigenvalue and EigenVector returned from Arpack
  // Extract EigenValue and EigenVector
  template<class ArProblem, class EigValue, class EigVector>
  Number extractEigenValueVector(ArProblem& arProb, EigValue& eigValue, EigVector& eigVector, Number nev, Number size, bool isReturnVector);

  // Extract EigenValue and EigenVector for non-symmetric problem
  template<class ArProblem, class EigValue, class EigVector>
  Number extractEigenValueVectorNonSym(ArProblem& arProb, EigValue& eigValue, EigVector& eigVector, Number nev, Number size, bool isReturnVector);

  // Set some parameters with user-defined value
  template<class ArProblem>
  void setParameters(ArProblem& arProb);
 protected:
  // Extract result EigenVector
  template<typename K, class ArProblem>
  void extractEigenVector(ArProblem& arProb, std::vector<K>& result, Number size, Number idx);

  // Extract result EigenVector of non-symmetric problem
  template<class ArProblem>
  void extractEigenVectorNonSym(ArProblem& arProb, std::vector<Complex>& result, Number size, Number idx, bool isNegative, bool isNegativeLeft);
  ...
};

The `EigenSolverStdArpackpp` class#

The EigenSolverStdArpackpp class enables the use of Arpack++ to solve the standard eigen problem \(Ax=\lambda x\) in all computational modes. Moreover, this class simplifies the post-processing stage, the found eigenvalues/eigenvectors are returned directly to user. However, advanced user can easily customize the some lines of code to get more features of Arpack++ (e.g Schur vector, etc, …). There are three prototypes of operator () corresponding to three cases: real symmetric, complex and real non-symmetric.

Real symmetric:

template<class Mat>
Number operator()(Mat& matA, Vector<Real>& eigValue, std::vector<Vector<Real> >& eigVector, Number nev, EigenComputationalMode compMode, Real sigma, const char* calMode = "LM")

Real non-symmetric:

template<class Mat>
Number operator()(Mat& matA, Vector<Complex>& eigValue, std::vector<Vector<Complex> >& eigVector, Number nev, EigenComputationalMode compMode, Complex sigma, const char* calMode = "LM")

Complex:

template<class Mat>
Number operator()(Mat& matA, Vector<Complex>& eigValue, std::vector<Vector<Complex> >& eigVector, Number nev, EigenComputationalMode compMode, Real sigma, const char* calMode = "LM")

Hint

If compMode isn’t _regular, only largeMatrix with storage of skylineStorage can be used as matrix input.

The operator () returns a value specifying the number of eigenvalues/eigenvectors being have found by Arpack++. This value is equal or less than the input \(nev\), which determines the number of eigenvalues/eigenvectors which user needs. sigma determines shift value in Shift mode (also in Buckling and Cayley mode). By default, \(nev\) eigenvalues with largest magnitude are turned; however, the function can return eigenvalues with smalles magnitude by changing \(calMode\) to “SM”.

The EigenComputationalMode is enumeration:

enum EigenComputationalMode {_regular, _buckling, _cayley, _shiftInv, _cplxShiftInv};

In order to solve eigen problem with Arpack++, it is necessary to define a Arpack++ problem then link this problem to XLiFE++ through Arpack++ interface layer (i.e LargeMatrixArpackpp, etc, …). There are three Arpack++ classes which are used to define a standard eigen problem: ARSymStdEig, ARNonSymStdEig, ARCompStdEig.

The following functions define Arpack++ problem and make link to XLiFE++ via Arpack++ interface layer.

template<typename K, class MatArpackpp, class Arpackpp, class Mat, class EigValue, class EigVector>
Number regMode(Mat& matA, Number size, Number nev, EigValue& eigValue, EigVector& eigVector, const char* calMode);
Number shfInvMode(Mat& matA, Number size, Number nev, K sigma, EigValue& eigValue, EigVector& eigVector, const char* calMode);
Number regModeNonSym(Mat& matA, Number size, Number nev, EigValue& eigValue, EigVector& eigVector, const char* calMode);
Number shfInvModeNonSym(Mat& matA, Number size, Number nev, K sigma, EigValue& eigValue, EigVector& eigVector, const char* calMode);

Hint

Users can easily customize the parameters of Arpack++ problem corresponding to their need.

See also

library=eigenSolvers,
header=
implementation=,
test=test_EigenSolverArpackpp.cpp,
header dep=EigenSolverBaseArpackpp.hpp, LargeMatrixArpackppStdReg.hpp, LargeMatrixArpackppStdShf.hpp

The `EigenSolverGenArpackpp` class#

The EigenSolverStdArpackpp class enables the use of Arpack++ to solve the generalized eigen problem \(Ax=\lambda Bx\) in all computational modes. Moreover, this class simplifies the post-processing stage, the found eigenvalues/eigenvectors are returned directly to user. However, advanced user can easily customize several lines of code to get more features of Arpack++ (e.g Schur vector, etc, …). There are three prototypes of operator () corresponding to three cases: real symmetric, complex and real non-symmetric.

Real symmetric:

template<class Mat>
Number operator()(const Mat& matA, Mat& matB, Vector<Real>& eigValue, std::vector<Vector<Real> >& eigVector, Number nev, EigenComputationalMode compMode, Real sigmaR, const char* calMode = "LM");

Real non-symmetric:

template<class Mat>
Number operator()(const Mat& matA, Mat& matB, Vector<Complex>& eigValue, std::vector<Vector<Complex> >& eigVector, Number nev, EigenComputationalMode compMode, Complex sigma, const char* calMode = "LM");

Complex:

template<class Mat>
Number operator()(const Mat& matA, Mat& matB, Vector<Complex>& eigValue, std::vector<Vector<Complex> >& eigVector, Number nev, EigenComputationalMode compMode, Real sigmaR, const char* calMode = "LM");

Real non-symmetric with complex shift:

template<class Mat>
Number operator()(const Mat& matA, Mat& matB, Vector<Complex>& eigValue, std::vector<Vector<Complex> >& eigVector, Number nev, EigenComputationalMode compMode, const char mode, Complex sigma, const char* calMode = "LM");

Hint

If compMode isn’t _regular, only largeMatrix with storage of skylineStorage can be used as matrix input.

Having matrix {B} as an additional input, the operator() of EigenSolverGenArpackpp class is called with the same parameters as the case of EigenSolverStdArpackpp. Moreover, it has one more prototype which is served for real non-symmetric case with complex shift value. Besides the same parameters like other prototypes, it is necessary for user to specify which part of complex shift value:real or imaginary, to use in the shift mode by parameter mode. This input can be R for real part or I for imaginary part of complex shift value sigma.

Same as EigenSolverStdArpackpp, it needs to define Arpack++ problems and link them to XLiFE++. ARSymGenEig, ARNonSymGenEig and ARCompGenEig class are used to define Arpack++ generalized problem.

The following functions define Arpack++ problem and make link to XLiFE++ via Arpack++ interface layer.

template<typename K, class MatArpackpp, class Arpackpp, class Mat, class EigValue, class EigVector>
Number regMode(const Mat& matA, const Mat& matB, Number size, Number nev, EigValue& eigValue, EigVector& eigVector, const char* calMode);
Number regModeNonSym(const Mat& matA, const Mat& matB, Number size, Number nev, EigValue& eigValue, EigVector& eigVector, const char* calMode);
Number symShfBucMode(const Mat& matA, Mat& matB, Number size, Number nev, K sigma, const char mode, EigValue& eigValue, EigVector& eigVector, const char* calMode);
Number symCaleyMode(const Mat& matA, Mat& matB, Number size, Number nev, K sigma, EigValue& eigValue, EigVector& eigVector, const char* calMode);
Number shfMode(const Mat& matA, Mat& matB, Number size, Number nev, K sigma, EigValue& eigValue, EigVector& eigVector, const char* calMode);
Number shfModeNonSym(const Mat& matA, Mat& matB, Number size, Number nev, K sigma, EigValue& eigValue, EigVector& eigVector, const char* calMode);
Number complexShfInvMode(const Mat& matA, Mat& matB, Number size, Number nev, Real sigmaR, Real sigmaI, const char mode, EigValue& eigValue, EigVector& eigVector, const char* calMode);

Hint

Users can easily customize the parameters of Arpack++ problem corresponding to their need.

See also

library=eigenSolvers,
header=
implementation=,
test=test_EigenSolverArpackpp.cpp,
header dep=EigenSolverBaseArpackpp.hpp, LargeMatrixArpackppGenReg.hpp, LargeMatrixArpackppGenShf.hpp

Examples#

The following matrices and vectors are used for all examples below.

const Number rowNum = 20;
const Number colNum = 20;
const Number nev = 15;
const String rMatrixDataSym("rSym20.data");
const String rMatrixDataSymPosDef("rSymPos20.data");
const String rMatrixDataNoSym("rnonSym20.data");
const String cMatrixDataSym("c20.data");
const String cMatrixData("c20b.data");
Vector<Real> rEigValue(colNum);
Vector<Complex> cEigValue(colNum);
std::vector<Vector<Real> > rEigVector(nev, rEigValue);
std::vector<Vector<Complex> > cEigVector(nev, cEigValue);

LargeMatrix<Real> rMatSkylineSym(dataPathTo(rMatrixDataNoSym), _dense, rowNum, colNum, _skyline, _sym);
LargeMatrix<Real> rMatSkylineSymSym(dataPathTo(rMatrixDataSym), _dense, rowNum, _skyline, _symmetric);

// A trick to make matrix A and matrix B have the same storage (just for testing purpose)
LargeMatrix<Real> rMatPosDefTemp(dataPathTo(rMatrixDataSymPosDef), _dense, rowNum, _skyline, _symmetric);
LargeMatrix<Real> rMatPosDefSym(rMatSkylineSymSym);
  for (Number i = 1; i < rowNum + 1; i++)
    for (Number j = 1; j < colNum + 1; j++)
    { rMatPosDefSym(i, j) = rMatPosDefTemp(i, j); }

// Positive definite in real non-symmetric case
LargeMatrix<Real> rMatPosDefNoSymTemp(dataPathTo(rMatrixDataSymPosDef), _dense, rowNum, colNum, _skyline, _sym);
LargeMatrix<Real> rMatPosDefNoSym(rMatSkylineSym);
  for (Number i = 1; i < rowNum + 1; i++)
    for (Number j = 1; j < colNum + 1; j++)
    { rMatPosDefNoSym(i, j) = rMatPosDefNoSymTemp(i, j); }

//  Complex Large Matrix
LargeMatrix<Complex> cMatSkylineSymSym(dataPathTo(cMatrixDataSym), _dense, rowNum, _skyline, _symmetric);

Real sigma = 10.0;
Complex csigma(-700, 20.0);
EigenSolverStdArpackpp eigenSolverStd;
EigenSolverGenArpackpp eigenSolverGen;

Standard eigenvalue problem#

Real symmetric case

// Regular mode
if ( 0 != eigenSolverStd(rMatSkylineSymSym, rEigValue, rEigVector, nev, _regular, sigma) ) {
  out << "Eigenvalue regular mode (Skyline Sym) : " << rEigValue << std::endl;
  out << "Eigenvector regular mode (Skyline Sym) : " << rEigVector << std::endl;
}
// Shift and invert mode
if ( 0 != eigenSolverStd(rMatSkylineSymSym, rEigValue, rEigVector, nev, _shiftInv, sigma) ) {
  out << "Eigenvalue shift and inverse mode :" << rEigValue << std::endl;
  out << "Eigenvalue shift and inverse mode: " << rEigVector << std::endl;
}

Real non -symmetric case

// Regular mode
if ( 0 != eigenSolverStd(rMatSkylineSym, cEigValue, cEigVector, nev, _regular, sigma)) {
  out << "Eigenvalue regular mode (Skyline Sym) : " << cEigValue << std::endl;
  out << "Eigenvector regular mode (Skyline Sym) : " << cEigVector << std::endl;
}
//Shift and invert mode
if ( 0 != eigenSolverStd(rMatSkylineSym, cEigValue, cEigVector, nev, _shiftInv, sigma)) {
  out << "Eigenvalue shift and inverse mode (Skyline Sym) :" << cEigValue << std::endl;
  out << "Eigenvalue shift and inverse mode (Skyline Sym): " << cEigVector << std::endl;
}

Complex case

// Regular mode
if ( 0 != eigenSolverStd(cMatSkylineSymSym, cEigValue, cEigVector, nev, _regular, csigma)) {
  out << "Eigenvalue regular mode (Skyline Sym) : " << cEigValue << std::endl;
  out << "Eigenvector regular mode (Skyline Sym) : " << cEigVector << std::endl;
}
// Shift and invert mode
if ( 0 != eigenSolverStd(cMatSkylineSymSym, cEigValue, cEigVector, nev, _shiftInv, csigma)) {
  out << "Eigenvalue shift and inverse mode (Skyline Sym) : " << cEigValue << std::endl;
  out << "Eigenvector shift and inverse mode (Skyline Sym) : " << cEigVector << std::endl;
}

Generalized eigenvalue problem#

Real symmetric case

// Regular mode
if ( 0 != eigenSolverGen(rMatSkylineSymSym, rMatPosDefSym, rEigValue, rEigVector, nev, _regular, sigma) ) {
out << "Eigenvalue regular mode (Skyline Sym):" << rEigValue << std::endl;
out << "Eigenvalue regular mode (Skyline Sym): " << rEigVector << std::endl;
}
// Shift and invert mode
if ( 0 != eigenSolverGen(rMatSkylineSymSym, rMatPosDefSym, rEigValue, rEigVector, nev, _shiftInv, sigma) ) {
  out << "Eigenvalue shift and inverse mode (Skyline Sym) :" << rEigValue << std::endl;
  out << "Eigenvalue shift and inverse mode (Skyline Sym) : " << rEigVector << std::endl;
}
// Buckling mode
if ( 0 != eigenSolverGen(rMatPosDefSym, rMatSkylineSymSym, rEigValue, rEigVector, nev, _buckling, sigma) ) {
  out << "Eigenvalue buckling mode (Skyline Sym) (positive definite):" << rEigValue << std::endl;
  out << "Eigenvalue buckling mode (Skyline Sym) (positive definite): " << rEigVector << std::endl;
}
// Cayley mode
if ( 0 != eigenSolverGen(rMatSkylineSymSym, rMatPosDefSym, rEigValue, rEigVector, nev, _cayley, sigma) ) {
  out << "Eigenvalue cayley mode (Skyline Sym):" << rEigValue << std::endl;
  out << "Eigenvalue cayley mode (Skyline Sym): " << rEigVector << std::endl;
}

Real non -symmetric case

// Regular mode
if ( 0 != eigenSolverGen(rMatSkylineSym, rMatPosDefSym, cEigValue, cEigVector, nev, _regular, sigma) ) {
  out << "Eigenvalue regular mode  (Skyline Dual) : " << cEigValue << std::endl;
  out << "Eigenvector regular mode (Skyline Dual) : " << cEigVector << std::endl;
}
// Shift and invert mode
if ( 0 != eigenSolverGen(rMatSkylineSym, rMatPosDefNoSym, cEigValue, cEigVector, nev, _shiftInv, sigma) ) {
  out << "Eigenvalue shift and inverse mode :" << cEigValue << std::endl;
  out << "Eigenvalue shift and inverse mode: " << cEigVector << std::endl;
}
// Complex shift with real part
if ( 0 != eigenSolverGen(rMatSkylineSym, rMatPosDefNoSym, cEigValue, cEigVector, nev, _cplxShiftInv, 'R', csigma) ) {
  out << "Eigenvalue complex shift and inverse mode (Real part) :" << cEigValue << std::endl;
  out << "Eigenvalue complex shift and inverse mode (Real part): " << cEigVector << std::endl;
}
// Complex shift with imaginary part
if ( 0 != eigenSolverGen(rMatSkylineSym, rMatPosDefNoSym, cEigValue, cEigVector, nev, _cplxShiftInv, 'I', csigma) ) {
  out << "Eigenvalue complex shift and inverse mode (Imaginary part):" << cEigValue << std::endl;
  out << "Eigenvalue complex shift and inverse mode (Imaginary part): " << cEigVector << std::endl;
}

Complex case

// Regular mode
if ( 0 != eigenSolverGen(cMatSkylineSymSym, cMatPosDef, cEigValue, cEigVector, nev, _regular, csigma)) {
  out << "Eigenvalue regular mode (Skyline Sym) : " << cEigValue << std::endl;
  out << "Eigenvector regular mode (Skyline Sym) : " << cEigVector << std::endl;
}
// Shift and invert mode
if ( 0 != eigenSolverGen(cMatSkylineSymSym, cMatPosDef, cEigValue, cEigVector, nev, _shiftInv, csigma, "LM")) {
  out << "Eigenvalue shift and inverse mode (Skyline Sym) : " << cEigValue << std::endl;
  out << "Eigenvector shift and inverse mode (Skyline Sym) : " << cEigVector << std::endl;
}

Besides making use of class EigenSolverStdArpackpp and class EigenSolverGenArpackpp, we can call Arpack++ problem object directly. However, it’s only recommended for advanced users!!

The following examples describe how to use Arpack++ object directly.

//1. Create interface object
LargeMatrixArpackppStdReg<Real, LargeMatrix<Real> > intMat(rMatSkylineSymSym, size);

//2. Create Arpack++ problem object
ARSymStdEig<Real, LargeMatrixArpackppStdReg<Real, LargeMatrix<Real> > >  arProb(size, nev, &intMat, &LargeMatrixArpackppStdReg<Real, LargeMatrix<Real> >::multOPx, (char *)"LM");

//3. Compute eigenvalues
Int nconv = arProb.FindEigenvalues();

for (Int i = 0; i < nconv; ++i) {
    out << arProb.Eigenvalue(i) << "  ";
}
out << "\n";
//4. Change some parameters
arProb.ChangeNev(nev-2);
arProb.ChangeTol(1.e-6);

nconv = arProb.FindEigenvectors();

for (Int i = 0; i < nconv; ++i) {
  out << arProb.Eigenvalue(i) << "  ";
}
out << "\n";

//4. Retrieve eigenVector
Real* eigenVec_p = 0;
for (Int i = 0; i < nconv; ++i) {
    eigenVec_p = arProb.RawEigenvector(i);
    out << "Eigenvector [" << i << "] :" ;
    for (Int j = 0; j < size; j++) {
        out <<   *eigenVec_p << "  ";
        eigenVec_p++;
    }
    out << "\n";
}

out << "\n";

Using Arpack#

Arpack++ interface layer#

The LargeMatrixArpackppStdReg class#

The LargeMatrixArpackppStdShf class#

The LargeMatrixArpackppGen class#

The LargeMatrixArpackppGenReg class#

The LargeMatrixArpackppSymGenReg class#

The LargeMatrixArpackppGenShf class#

The LargeMatrixArpackppGenBuc class#

The LargeMatrixArpackppGenCay class#

The LargeMatrixArpackppGenCsh class#

Users interface layer#

The EigenSolverBaseArpackpp class#

The EigenSolverStdArpackpp class#

The EigenSolverGenArpackpp class#

Examples#

Standard eigenvalue problem#

Generalized eigenvalue problem#

The `LargeMatrixArpackppStdReg` class#

The `LargeMatrixArpackppStdShf` class#

The `LargeMatrixArpackppGen` class#

The `LargeMatrixArpackppGenReg` class#

The `LargeMatrixArpackppSymGenReg` class#

The `LargeMatrixArpackppGenShf` class#

The `LargeMatrixArpackppGenBuc` class#

The `LargeMatrixArpackppGenCay` class#

The `LargeMatrixArpackppGenCsh` class#

The `EigenSolverBaseArpackpp` class#

The `EigenSolverStdArpackpp` class#

The `EigenSolverGenArpackpp` class#