Installation of external eigensolvers

InitMatrix

Syntax :

  void InitMatrix(Matrix& K);
  void InitMatrix(Matrix& K, Matrix& M);

This method allows the initialization of pointers for the stiffness matrix and mass matrix. It is mandatory to call it when using SparseEigenProblem, DenseEigenProblem or VirtualEigenProblem. K is the stiffness matrix and M the mass matrix, we search a scalar λ and a non-null vector x such that
K x = λ M x
If M is not provided, it is assumed to be equal to the identity matrix.

Example :

{

// declaration of the eigenproblem
DenseEigenProblem<double, double, General, RowMajor> var_eig;
Vector<double> eigen_values, eigen_imag;
Matrix<double> eigen_vec;

// you can set some parameters
var_eig.SetNbAskedEigenvalues(5);
var_eig.SetComputationalMode(var_eig.REGULAR_MODE);
var_eig.SetTypeSpectrum(var_eig.LARGE_EIGENVALUES, 0);

// then you can construct the stiffness matrix
int n = 20;
Matrix<double, General, RowMajor> K(n, n);
K.FillRand();

// and you provide this matrix to the eigenproblem
// standard eigenvalue problem K x = lambda x
var_eig.InitMatrix(K);

// then you can call GetEigenvaluesEigenvectors
GetEigenvaluesEigenvectors(var_eig, lambda, lambda_imag, eigen_vec);

// for a generalized eigenvalue problem, provide M
Matrix<double, Symmetric, RowSymPacked> M(n, n);
M.SetIdentity();
// searching K x = lambda M x
var_eig.InitMatrix(K, M);

// then you can call GetEigenvaluesEigenvectors
GetEigenvaluesEigenvectors(var_eig, lambda, lambda_imag, eigen_vec);

// if the type of M is different from Matrix<double, Symmetric, RowSymPacked>
// and is Matrix<Tm, PropM, StorageM>
// use DenseEigenProblem<double, double, General, RowMajor, Tm, PropM, StorageM> var_eig;
}

{
// for a sparse eigenproblem, similar stuff
// declaration of the eigenproblem
SparseEigenProblem<double, Matrix<double, General, ArrayRowSparse> > var_eig;
Vector<double> eigen_values, eigen_imag;
Matrix<double> eigen_vec;

// you can set some parameters
var_eig.SetNbAskedEigenvalues(5);
var_eig.SetComputationalMode(var_eig.REGULAR_MODE);
var_eig.SetTypeSpectrum(var_eig.LARGE_EIGENVALUES, 0);

// then you can construct the stiffness matrix
int n = 20;
Matrix<double, General, ArrayRowSparse> K(n, n);
K.FillRand();

// and you provide this matrix to the eigenproblem
// standard eigenvalue problem K x = lambda x
var_eig.InitMatrix(K);

// then you can call GetEigenvaluesEigenvectors
GetEigenvaluesEigenvectors(var_eig, lambda, lambda_imag, eigen_vec);

// for a generalized eigenvalue problem, provide M
Matrix<double, Symmetric, ArrayRowSymSparse> M(n, n);
M.SetIdentity();
// searching K x = lambda M x
var_eig.InitMatrix(K, M);

// then you can call GetEigenvaluesEigenvectors
GetEigenvaluesEigenvectors(var_eig, lambda, lambda_imag, eigen_vec);

// if the type of M is different from Matrix<double, Symmetric, ArrayRowSymSparse>
// and is Matrix<Tm, PropM, StorageM>
// use SparseEigenProblem<double, Matrix<double, General, ArrayRowSparse>,
//                                Matrix<Tm, PropM, StorageM> > var_eig;
}

{
  // matrix-free eigenproblem
  // stiffness matrix has type MyMatrixClass which derives from VirtualMatrix
  VirtualEigenProblem<double> var_eig;

  // you can set some parameters
  var_eig.SetNbAskedEigenvalues(5);
  var_eig.SetComputationalMode(var_eig.REGULAR_MODE);
  var_eig.SetTypeSpectrum(var_eig.LARGE_EIGENVALUES, 0);
 
  // then you can construct the stiffness matrix
  MyMatrixClass K;

  // and you provide this matrix to the eigenproblem
  // standard eigenvalue problem K x = lambda x
  var_eig.InitMatrix(K);

  // then you can call GetEigenvaluesEigenvectors
  GetEigenvaluesEigenvectors(var_eig, lambda, lambda_imag, eigen_vec);
}

Location :

GetComputationalMode, SetComputationalMode

Syntax :

  void SetComputationalMode(int);
  int GetComputationalMode();

SetComputationalMode sets the computational mode to use for the research of eigenvalues and eigenvectors, whereas GetComputationalMode returns the computational mode used. This computational mode is used only for Arpack or Slepc eigenvalue solver. The default computational mode is REGULAR_MODE. This mode is particularly well suited if you are looking for the largest eigenvalues of the matrix. However, it can induce a lot of iterations to converge if smallest or clustered eigenvalues are researched. Two major computational modes are available :

REGULAR_MODE : we solve K x = λ M x
This mode will perform only matrix-vector products with M or K. In that mode, the mass matrix M must be real symmetric positive.
INVERT_MODE : we compute eigenvalues μ of the following eigenvalue problem (K - σ M)^-1 M x = μ x
In that mode, the mass matrix M and stiffness matrix K can be any type of matrices. The eigenvalue λ of the original eigenvalue problem is equal to σ + 1/μ. This mode will perform linear solves and is well suited to find eigenvalues close to σ. If large eigenvalues are searched, this mode will compute the eigenvalues of M^-1 K, requiring an inversible mass matrix.

In regular use, the user should prefer REGULAR_MODE to compute large eigenvalues, whereas INVERT_MODE is more appropriate to compute clustered eigenvalues (around a shift σ). If you are using Slepc and you have a mass matrix M different from the identity, you will need to use INVERT_MODE to compute large eigenvalues.

A drawback of INVERT_MODE is that we have to compute eigenvalues of a non-symmetric matrix even though the original matrices K and M are symmetric. In order to overcome this issue, there are two solutions : using the Cholesky factorization of M or using other computational modes (provided only by Arpack). We describe these two solutions in the following paragraphs

Factorization of M

If the mass matrix is diagonal (M = D) with positive terms, we search the eigenvalues of
D^-1/2 K D^-1/2 x = λ x
By this way, we have to solve a standard symmetric eigenproblem. You have to call the method SetDiagonalMass in order to inform that the mass matrix is a positive diagonal and use this approach

If the mass matrix is definite positive, we can compute the Cholesky factorization M = L L^T to obtain the following eigenvalue problem
L^-1 K L^-T x = λ x
which is a standard symmetric eigenproblem. If you want to use this approach, you have to call the metho SetCholeskyFactoForMass.

Other modes for Arpack

If you are using Arpack, and if M is a real symmetric positive matrix, you can use the following modes (we put between parenthesis the mode number used in Arpack) :

SHIFTED_MODE (3 for Arpack) : operator (K - σ M)^-1 M is used
It uses real part of (K- σ M)^-1 M for complex shifts.
IMAG_SHIFTED_MODE (4 for Arpack, only for real unsymmetric matrices)
It uses imaginary part of (K- σ M)^-1 M for complex shifts.
BUCKLING_MODE (4 for Arpack) : operator (K - σ M)^-1 K is used
In that mode, the mass matrix M and stiffness matrix K must be real symmetric positive.
CAYLEY_MODE (5 for Arpack) : operator (K - σ M)^-1 (K + σ M) is used
In that mode, the mass matrix M must be real symmetric positive, the stiffness matrix K must be real symmetric.

Example :

// declaration of the eigenproblem
SparseEigenProblem<double, Matrix<double, General, ArrayRowSparse> > var_eig;
Vector<double> eigen_values, eigen_imag;
Matrix<double> eigen_vec;

// you can set some parameters
var_eig.SetNbAskedEigenvalues(5);
var_eig.SetTypeSpectrum(var_eig.LARGE_EIGENVALUES, 0);

// and choose the computational mode relevant with researched spectrum
var_eig.SetComputationalMode(var_eig.REGULAR_MODE);

// then you can construct the stiffness matrix
int n = 20;
Matrix<double, General, ArrayRowSparse> K(n, n);
K.FillRand();

// and you provide this matrix to the eigenproblem
// standard eigenvalue problem K x = lambda x
var_eig.InitMatrix(K);

// then you can call GetEigenvaluesEigenvectors
GetEigenvaluesEigenvectors(var_eig, lambda, lambda_imag, eigen_vec);

// for a generalized eigenvalue problem, provide M
Matrix<double, Symmetric, ArrayRowSymSparse> M(n, n);
M.SetIdentity();
// searching K x = lambda M x
var_eig.InitMatrix(K, M);

// searching eigenvalues close to a shift
var_eig.SetTypeSpectrum(var_eig.CENTERED_EIGENVALUES, 0.5);
var_eig.SetComputationalMode(var_eig.SHIFTED_MODE);

// then you can call GetEigenvaluesEigenvectors
GetEigenvaluesEigenvectors(var_eig, lambda, lambda_imag, eigen_vec);

Location :

GetCommunicator

Syntax :

  const MPI_Comm& GetCommunicator() const;
  void SetCommunicator(const MPI_Comm&);
  int GetRankCommunicator() const;

The method SetCommunicator sets the MPI communicator (corresponding to L3_comm for Feast) shared by rows of the matrix. The method GetCommunicator returns the MPI communicator associated with the eigenvalue problem. For specialized classes such as SparseEigenProblem, it is not necessary to call SetCommunicator, since the communicator will be copied from the distributed matrices. The method GetRankCommunicator returns the rank of the current processor in the MPI communicator.

Example :

// filling distributed matrices
DistributedMatrix<complex<double>, General, ArrayRowSparse> M, K;
// MPI communicator is given when filling these matrices

// for a parallel computation, DistributedMatrix are needed
SparseEigenProblem<complex<double>,
                   DistributedMatrix<complex<double>, General, ArrayRowSparse>,
                   DistributedMatrix<complex<double>, General, ArrayRowSparse> > var_eig;

var_eig.SetStoppingCriterion(1e-12);
var_eig.SetNbAskedEigenvalues(20);
var_eig.InitMatrix(K, M); // here the communicator is copied
var_eig.SetPrintLevel(4);

// you can retrieve the communicator
MPI_Comm comm = var_eig.GetCommunicator();
// and the rank
int rank = var_eig.GetRankCommunicator();

complex<double> center(0.1, 2.6);
var_eig.SetTypeSpectrum(var_eig.CENTERED_EIGENVALUES, center);
var_eig.SetComputationalMode(var_eig.INVERT_MODE);

Vector<complex<double> > lambda, lambda_imag;
Matrix<complex<double>, General, ColMajor> eigen_vec;

// we call Slepc to find eigenvalues
FindEigenvaluesSlepc(var_eig, lambda, lambda_imag, eigen_vec);

Location :

GetRciCommunicator

Syntax :

  const MPI_Comm& GetRciCommunicator() const;

This method returns the intermediate communicator needed by pfeast (L2_comm). This method does not need to be used in regular use.

Example :

// filling distributed matrices
DistributedMatrix<complex<double>, General, ArrayRowSparse> M, K;
// MPI communicator is given when filling these matrices

// for a parallel computation, DistributedMatrix are needed
SparseEigenProblem<complex<double>, DistributedMatrix<complex<double>, General, ArrayRowSparse>,
                   DistributedMatrix<complex<double>, General, ArrayRowSparse> > var_eig;

// SetGlobalCommunicator is mandatory when using pfeast : you provide the communicator
// grouping all the processors that will perform the eigenvalue computation (L1_comm for Feast)
var_eig.SetGlobalCommunicator(MPI_COMM_WORLD);
var_eig.SetStoppingCriterion(1e-12);
var_eig.SetNbAskedEigenvalues(20);
var_eig.InitMatrix(K, M); // here the communicator for matrices is copied (L3_comm for Feast))
var_eig.SetPrintLevel(4);

// you can retrieve the communicator
MPI_Comm comm = var_eig.GetCommunicator();

// you can retrieve the L2_comm that will be provided to Feast
MPI_Comm L2_comm = var_eig.GetRciCommunicator();

complex<double> center(0.1, 2.6);
FeastParam& param = var_eig.GetFeastParameters();
param.SetCircleSpectrum(center, 0.2);

Vector<complex<double> > lambda, lambda_imag;
Matrix<complex<double>, General, ColMajor> eigen_vec;

// we call Feast to find eigenvalues
FindEigenvaluesFeast(var_eig, lambda, lambda_imag, eigen_vec);

Location :

GetGlobalCommunicator, SetGlobalCommunicator

Syntax :

  const MPI_Comm& GetGlobalCommunicator() const;
  void SetGlobalCommunicator(const MPI_Comm& comm) const;
  int GetGlobalRankCommunicator() const;

The method GetGlobalCommunicator returns the global communicator needed by pfeast (L1_comm). The method SetGlobalCommunicator sets this communicator. It is mandatory to call this method before providing distributed matrices. The method GetRankCommunicator returns the rank of the current processor in the global MPI communicator.

Example :

// filling distributed matrices
DistributedMatrix<complex<double>, General, ArrayRowSparse> M, K;
// MPI communicator is given when filling these matrices

// for a parallel computation, DistributedMatrix are needed
SparseEigenProblem<complex<double>, DistributedMatrix<complex<double>, General, ArrayRowSparse>,
                   DistributedMatrix<complex<double>, General, ArrayRowSparse> > var_eig;

// SetGlobalCommunicator is mandatory when using pfeast : you provide the communicator
// grouping all the processors that will perform the eigenvalue computation (L1_comm for Feast)
var_eig.SetGlobalCommunicator(MPI_COMM_WORLD);
var_eig.SetStoppingCriterion(1e-12);
var_eig.SetNbAskedEigenvalues(20);
var_eig.InitMatrix(K, M); // here the communicator for matrices is copied (L3_comm for Feast))
var_eig.SetPrintLevel(4);

// you can retrieve the global communicator
MPI_Comm comm_glob = var_eig.GetGlobalCommunicator();
// and its associated rank
int rank_glob = var_eig.GetGlobalRankCommunicator();

complex<double> center(0.1, 2.6);
FeastParam& param = var_eig.GetFeastParameters();
param.SetCircleSpectrum(center, 0.2);

Vector<complex<double> > lambda, lambda_imag;
Matrix<complex<double>, General, ColMajor> eigen_vec;

// we call Feast to find eigenvalues
FindEigenvaluesFeast(var_eig, lambda, lambda_imag, eigen_vec);

Location :

GetNbAskedEigenvalues

Syntax :

  void SetNbAskedEigenvalues(int nev);
  int GetNbAskedEigenvalues();

SetNbAskedEigenvalues sets the number of converged eigenvalues you desire to know, whereas GetNbAskedEigenvalues returns the number of desired eigenvalues. This method is mandatory since the default number of eigenvalues is equal to 0.

Example :

// declaration of the eigenproblem
SparseEigenProblem<double, Matrix<double, General, ArrayRowSparse> > var_eig;
Vector<double> eigen_values, eigen_imag;
Matrix<double> eigen_vec;

// you can set some parameters
var_eig.SetNbAskedEigenvalues(5);
var_eig.SetTypeSpectrum(var_eig.LARGE_EIGENVALUES, 0);

// and choose the computational mode relevant with researched spectrum
var_eig.SetComputationalMode(var_eig.REGULAR_MODE);

// then you can construct the stiffness matrix
int n = 20;
Matrix<double, General, ArrayRowSparse> K(n, n);
K.FillRand();

// and you provide this matrix to the eigenproblem
// standard eigenvalue problem K x = lambda x
var_eig.InitMatrix(K);

// then you can call GetEigenvaluesEigenvectors
GetEigenvaluesEigenvectors(var_eig, lambda, lambda_imag, eigen_vec);

Location :

GetNbAdditionalEigenvalues, SetNbAdditionalEigenvalues

Syntax :

  void SetNbAdditionalEigenvalues(int nev);
  int GetNbAdditionalEigenvalues();

SetNbAdditionalEigenvalues sets the number of additional eigenvalues to avoid problems in Arpack, whereas GetNbAdditionalEigenvalues returns the number of additional eigenvalues. Usually Arpack works correctly, but if you see that the recover of eigenvalues generates a segmentation fault, you can provide additional eigenvalues, it may fix the problem..

Example :

// declaration of the eigenproblem
SparseEigenProblem<double, Matrix<double, General, ArrayRowSparse> > var_eig;
Vector<double> eigen_values, eigen_imag;
Matrix<double> eigen_vec;

// you can set some parameters
var_eig.SetNbAskedEigenvalues(5);
var_eig.SetNbAdditionalEigenvalues(1);
var_eig.SetTypeSpectrum(var_eig.LARGE_EIGENVALUES, 0);

// and choose the computational mode relevant with researched spectrum
var_eig.SetComputationalMode(var_eig.REGULAR_MODE);

// then you can construct the stiffness matrix
int n = 20;
Matrix<double, General, ArrayRowSparse> K(n, n);
K.FillRand();

// and you provide this matrix to the eigenproblem
// standard eigenvalue problem K x = lambda x
var_eig.InitMatrix(K);

// then you can call GetEigenvaluesEigenvectors
GetEigenvaluesEigenvectors(var_eig, lambda, lambda_imag, eigen_vec);

Location :

GetAnasaziParameters

Syntax :

  AnasaziParam& GetAnasaziParameters();

This method returns the parameters associated with Anasazi solver.

Example :

// declaration of the eigenproblem
SparseEigenProblem<double, Matrix<double, General, ArrayRowSparse> > var_eig;
Vector<double> eigen_values, eigen_imag;
Matrix<double> eigen_vec;

// you can set some parameters
var_eig.SetNbAskedEigenvalues(5);
var_eig.SetTypeSpectrum(var_eig.LARGE_EIGENVALUES, 0);

// and also parameters specific to Anasazi
AnasaziParam& param = var_eig.GetAnasaziParameters();
param.SetEigensolverType(param.SOLVER_BKS);

// and choose the computational mode relevant with researched spectrum
var_eig.SetComputationalMode(var_eig.REGULAR_MODE);

// then you can construct the stiffness matrix
int n = 20;
Matrix<double, General, ArrayRowSparse> K(n, n);
K.FillRand();

// and you provide this matrix to the eigenproblem
// standard eigenvalue problem K x = lambda x
var_eig.InitMatrix(K);

// then you can call GetEigenvaluesEigenvectors
GetEigenvaluesEigenvectors(var_eig, lambda, lambda_imag, eigen_vec, TypeEigenvalueSolver::ANASAZI);

Location :

GetSlepcParameters

Syntax :

  SlepcParam& GetSlepcParameters();

This method returns the parameters associated with SLEPc solver.

Example :

// declaration of the eigenproblem
SparseEigenProblem<double, Matrix<double, General, ArrayRowSparse> > var_eig;
Vector<double> eigen_values, eigen_imag;
Matrix<double> eigen_vec;

// you can set some parameters
var_eig.SetNbAskedEigenvalues(5);
var_eig.SetTypeSpectrum(var_eig.LARGE_EIGENVALUES, 0);

// and also parameters specific to SLEPc
SlepcParam& param = var_eig.GetSlepcParameters();
param.SetEigensolverType(param.KRYLOVSCHUR);

// and choose the computational mode relevant with researched spectrum
var_eig.SetComputationalMode(var_eig.REGULAR_MODE);

// then you can construct the stiffness matrix
int n = 20;
Matrix<double, General, ArrayRowSparse> K(n, n);
K.FillRand();

// and you provide this matrix to the eigenproblem
// standard eigenvalue problem K x = lambda x
var_eig.InitMatrix(K);

// then you can call GetEigenvaluesEigenvectors
GetEigenvaluesEigenvectors(var_eig, lambda, lambda_imag, eigen_vec, TypeEigenvalueSolver::SLEPC);

Location :

GetFeastParameters

Syntax :

  FeastParam& GetFeastParameters();

This method returns the parameters associated with FEAST solver.

Example :

// declaration of the eigenproblem
SparseEigenProblem<double, Matrix<double, General, ArrayRowSparse> > var_eig;
Vector<double> eigen_values, eigen_imag;
Matrix<double> eigen_vec;

// you can set some parameters
complex<double> center(0.4, 0.2);
var_eig.SetNbAskedEigenvalues(5);
var_eig.SetTypeSpectrum(var_eig.CENTERED_EIGENVALUES, center);

// and also parameters specific to Feast
FeastParam& param = var_eig.GetFeastParameters();
param.SetCircleSpectrum(center, 0.2);

// and you provide this matrix to the eigenproblem
// standard eigenvalue problem K x = lambda x
var_eig.InitMatrix(K);

// then you can call GetEigenvaluesEigenvectors
GetEigenvaluesEigenvectors(var_eig, lambda, lambda_imag, eigen_vec, TypeEigenvalueSolver::FEAST);

Location :

GetNbBlocks, SetNbBlocks

Syntax :

  void SetNbBlocks(int n);
  int GetNbBlocks();

SetNbBlocks sets the number of blocks to use in Anasazi, whereas GetNbBlocks returns the number of blocks.

Example :

  SparseEigenProblem<T, MatrixK, MatrixM> var_eig;

  var_eig.SetStoppingCriterion(1e-14);
  var_eig.SetNbMaximumIterations(20000);
  var_eig.SetNbAskedEigenvalues(nb_eigenval);
  var_eig.SetComputationalMode(var_eig.REGULAR_MODE);

  var_eig.InitMatrix(K, M);

  AnasaziParam& param = var_eig.GetAnasaziParameters();
  param.SetEigensolverType(param.SOLVER_BKS);

  param.SetNbMaximumRestarts(200);
  param.SetNbBlocks(10);
  var_eig.SetNbArnoldiVectors(nb_eigenval+1);

  // Large eigenvalues
  var_eig.SetTypeSpectrum(var_eig.LARGE_EIGENVALUES, 0, var_eig.SORTED_MODULUS);
  
  GetEigenvaluesEigenvectors(var_eig, lambda, lambda_imag, eigen_vec);
  DISP(lambda); DISP(lambda_imag);

Location :

GetNbMaximumRestarts, SetNbMaximumRestarts

Syntax :

  void SetNbMaximumRestarts(int n);
  int GetNbMaximumRestarts();

SetNbMaximumRestarts sets the maximal number of restarts allowed in Anasazi, whereas GetNbMaximumRestarts returns this number.

Example :

  SparseEigenProblem<T, MatrixK, MatrixM> var_eig;

  var_eig.SetStoppingCriterion(1e-14);
  var_eig.SetNbMaximumIterations(20000);
  var_eig.SetNbAskedEigenvalues(nb_eigenval);
  var_eig.SetComputationalMode(var_eig.REGULAR_MODE);

  var_eig.InitMatrix(K, M);

  AnasaziParam& param = var_eig.GetAnasaziParameters();
  param.SetEigensolverType(param.SOLVER_BKS);

  param.SetNbMaximumRestarts(200);
  param.SetNbBlocks(10);
  var_eig.SetNbArnoldiVectors(nb_eigenval+1);

  // Large eigenvalues
  var_eig.SetTypeSpectrum(var_eig.LARGE_EIGENVALUES, 0, var_eig.SORTED_MODULUS);
  
  GetEigenvaluesEigenvectors(var_eig, lambda, lambda_imag, eigen_vec);
  DISP(lambda); DISP(lambda_imag);

Location :

GetOrthoManager

Syntax :

    int GetOrthoManager();

This method returns the orthogonal manager to use in Anasazi solvers.

Location :

GetEigensolverType, SetEigensolverType

Syntax :

  void SetEigensolverType(int type);
  int GetEigensolverType();

SetEigensolverType sets the type of eigensolver to use in Anasazi, whereas GetEigensolverType returns this number. You can choose between SOLVER_LOBPCG, SOLVER_BKS and SOLVER_BD.

Example :

  SparseEigenProblem<T, MatrixK, MatrixM> var_eig;

  var_eig.SetStoppingCriterion(1e-14);
  var_eig.SetNbMaximumIterations(20000);
  var_eig.SetNbAskedEigenvalues(nb_eigenval);
  var_eig.SetComputationalMode(var_eig.REGULAR_MODE);

  var_eig.InitMatrix(K, M);
  var_eig.SetNbArnoldiVectors(nb_eigenval+1);

  AnasaziParam& param = var_eig.GetAnasaziParameters();
  param.SetEigensolverType(param.SOLVER_BKS);
  param.SetNbBlocks(10);
  param.SetNbMaximumRestarts(200);        

  // Large eigenvalues
  var_eig.SetTypeSpectrum(var_eig.LARGE_EIGENVALUES, 0, var_eig.SORTED_MODULUS);
  
  GetEigenvaluesEigenvectors(var_eig, lambda, lambda_imag, eigen_vec, TypeEigenvalueSolver::ANASAZI);
  DISP(lambda); DISP(lambda_imag);

Location :

GetTypeSpectrum, SetTypeSpectrum

Syntax :

  void SetTypeSpectrum(int type, T sigma);
  void SetTypeSpectrum(int type, T sigma, int type_sort);
  int GetTypeSpectrum();

SetTypeSpectrum sets the part of the spectrum you desire to know, whereas GetTypeSpectrum returns the type of researched spectrum. You can ask LARGE_EIGENVALUES, SMALL_EIGENVALUES and CENTERED_EIGENVALUES, to research largest eigenvalues, smallest eigenvalues or eigenvalues closest to the shift sigma. If you search the largest or smallest eigenvalues, the shift sigma will not be used, and you can only use REGULAR_MODE (for linear eigenproblems). In the case of small eigenvalues, this mode may induce a high number of iterations, it can be a good idea to search eigenvalues close to a small value. type_sort can be used to specify how values are sorted, you can search largest eigenvalues by their modulus, real part or imaginary part (SORTED_MODULUS, SORTED_REAL and SORTED_IMAG). SORTED_MODULUS is the default sorting strategy.

Example :

// declaration of the eigenproblem
SparseEigenProblem<double, Matrix<double, General, ArrayRowSparse> > var_eig;
Vector<double> eigen_values, eigen_imag;
Matrix<double> eigen_vec;

// you can set some parameters
var_eig.SetNbAskedEigenvalues(5);
var_eig.SetTypeSpectrum(var_eig.LARGE_EIGENVALUES, 0);

// for eigenvalues closest to 0.25+0.3i
var_eig.SetTypeSpectrum(var_eig.CENTERED_EIGENVALUES, complex<double>(0.25, 0.3), SORTED_IMAG);

if (var_eig.GetTypeSpectrum() != var_eig.CENTERED_EIGENVALUES)
  cout << "not possible" << endl;

// and choose the computational mode relevant with researched spectrum
var_eig.SetComputationalMode(var_eig.SHIFTED_MODE);

// then you can construct the stiffness matrix
int n = 20;
Matrix<double, General, ArrayRowSparse> K(n, n);
K.FillRand();

// and you provide this matrix to the eigenproblem
// standard eigenvalue problem K x = lambda x
var_eig.InitMatrix(K);

// then you can call GetEigenvaluesEigenvectors
GetEigenvaluesEigenvectors(var_eig, lambda, lambda_imag, eigen_vec);

Location :

GetTypeSorting

Syntax :

  int GetTypeSorting();

This methods returns the sorting strategy used (SORTED_MODULUS, SORTED_REAL or SORTED_IMAG).

Example :

// declaration of the eigenproblem
SparseEigenProblem<double, Matrix<double, General, ArrayRowSparse> > var_eig;
Vector<double> eigen_values, eigen_imag;
Matrix<double> eigen_vec;

// you can set some parameters
var_eig.SetNbAskedEigenvalues(5);
var_eig.SetTypeSpectrum(var_eig.LARGE_EIGENVALUES, 0);

if (var_eig.GetTypeSorting() != var_eig.SORTED_MODULUS)
  cout << "not possible" << endl;

// and choose the computational mode relevant with researched spectrum
var_eig.SetComputationalMode(var_eig.REGULAR_MODE);

// then you can construct the stiffness matrix
int n = 20;
Matrix<double, General, ArrayRowSparse> K(n, n);
K.FillRand();

// and you provide this matrix to the eigenproblem
// standard eigenvalue problem K x = lambda x
var_eig.InitMatrix(K);

// then you can call GetEigenvaluesEigenvectors
GetEigenvaluesEigenvectors(var_eig, lambda, lambda_imag, eigen_vec);

GetShiftValue, GetImagShiftValue

Syntax :

  T GetShiftValue();
  T GetImagShiftValue();
  void SetShiftValue(T sr);
  void SetImagShiftValue(T si);

The method GetShiftValue returns the shift (real or complex). For real unsymmetric matrices, you can also retrieve the imaginary part of the shift. For complex eigenproblems, the imaginary part is not relevant since the shift is already complex. You can modify the shift by calling SetShiftValue. Usually the shift is set by calling SetTypeSpectrum.

Example :

// declaration of the eigenproblem
SparseEigenProblem<double, Matrix<double, General, ArrayRowSparse> > var_eig;
Vector<double> eigen_values, eigen_imag;
Matrix<double> eigen_vec;

// you can set some parameters
var_eig.SetNbAskedEigenvalues(5);

// for eigenvalues closest to 0.25+0.3i
var_eig.SetTypeSpectrum(var_eig.CENTERED_EIGENVALUES, complex<double>(0.25, 0.3), SORTED_IMAG);

// shift_r = 0.25 and shift_i = 0.3
double shift_r = var_eig.GetShiftValue();
double shift_i = var_eig.GetImagShiftValue();

// if you want to find eigenvalues closes to 0.4+0.2i, you can modify the values
var_eig.SetShiftValue(0.4);
var_eig.SetImagShiftValue(0.2);

GetComplexShift

Syntax :

  void GetComplexShift(T sr, T si, T& s)

This methods sets s = sr + I si if sr and si are real, and s = sr if sr and si are complex. This method should not be called in regular use.

SetIntervalSpectrum

Syntax :

  void SetIntervalSpectrum(double Lmin, double Lmax);
  double GetLowerBoundInterval();
  double GetUpperBoundInterval();

The method SetIntervalSpectrum sets the interval where eigenvalues are searched (relevant only for an hermitian eigenproblem).

Example :

  Matrix<double, Symmetric, ArrayRowSymSparse> K; // K needs to be filled

  SparseEigenProblem<T, Matrix<double, Symmetric, ArrayRowSymSparse> > var_eig;
  
  var_eig.SetStoppingCriterion(1e-12);
  var_eig.SetNbAskedEigenvalues(nb_eigenval);  
  var_eig.InitMatrix(K);  
  
  // eigenvalues are searched within a given interval
  FeastParam& param = var_eig.GetFeastParameters();
  param.SetIntervalSpectrum(-1.7, -1.0);

  cout << "Interval = [" << param.GetLowerBoundInterval() << ", " << param.GetUpperBoundInterval() << "]" << endl;

  Vector<double> lambda, lambda_imag;
  Matrix<double, General, ColMajor> eigen_vec;
  
  FindEigenvaluesFeast(var_eig, lambda, lambda_imag, eigen_vec);

SetCholeskyFactoForMass, UseCholeskyFactoForMass

Syntax :

  void SetCholeskyFactoForMass(bool);
  void SetCholeskyFactoForMass();
  bool UseCholeskyFactoForMass();

When solving generalized eigenvalue problems K x = λ M x, if K and M are symmetric, it can be attractive to perform a Cholesky factorization of M = L L^T, and consider the standard problem L^-1 K L^-T x = λ x. If you want to use that strategy, SetCholeskyFactoForMass has to be called. UseCholeskyFactoForMass returns true if this strategy is used.

Example :

// declaration of the eigenproblem
SparseEigenProblem<double, Matrix<double, General, ArrayRowSparse> > var_eig;
Vector<double> eigen_values, eigen_imag;
Matrix<double> eigen_vec;

// you can set some parameters
var_eig.SetNbAskedEigenvalues(5);
var_eig.SetTypeSpectrum(var_eig.CENTERED_EIGENVALUES, 0.4);

// then use Cholesky factorization
var_eig.SetCholeskyFactoForMass();
if (var_eig.UseCholekyFactoForMass())
    cout << "Searching eigenvalues of L^-1  K L^-T" << endl;

// if you want to cancel that, or solve another problem without Cholesky
var_eig.SetCholeskyFactoForMass(false);

// then declare a generalized problem
var_eig.InitMatrix(K, M);

SetDiagonalMass, DiagonalMass

Syntax :

  void SetDiagonalMass();
  void SetDiagonalMass(bool);
  bool DiagonalMass();

When solving generalized eigenvalue problems K x = λ M x, if K and M are symmetric and M diagonal, it can be attractive to consider the standard problem M^-1/2 K M^-1/2 x = λ x. If you want to use that strategy, SetDiagonalMass has to be called. DiagonalMass returns true if this strategy is used.

Example :

// declaration of the eigenproblem
SparseEigenProblem<double, Matrix<double, General, ArrayRowSparse> > var_eig;
Vector<double> eigen_values, eigen_imag;
Matrix<double> eigen_vec;

// you can set some parameters
var_eig.SetNbAskedEigenvalues(5);
var_eig.SetTypeSpectrum(var_eig.CENTERED_EIGENVALUES, 0.4);

// then inform that mass matrix is diagonal
var_eig.SetDiagonalMass();
if (var_eig.DiagonalMass())
    cout << "Searching eigenvalues of M^-1/2  K M^-1/2" << endl;

// if you want to cancel that, or solve another problem without diagonal mass
var_eig.SetDiagonalMass(false);

// then declare a generalized problem
var_eig.InitMatrix(K, M);

Location :

SetStoppingCriterion, GetStoppingCriterion

Syntax :

  void SetStoppingCriterion(double);
  double GetStoppingCriterion();

With those methods, you can modify and retrieve the stopping criterion used by the iterative algorithm. The default value is equal to 1e-6.

Example :

// declaration of the eigenproblem
SparseEigenProblem<double, Matrix<double, General, ArrayRowSparse> > var_eig;
Vector<double> eigen_values, eigen_imag;
Matrix<double> eigen_vec;

// you can set some parameters
var_eig.SetNbAskedEigenvalues(5);
var_eig.SetStoppingCriterion(1e-12);
var_eig.SetTypeSpectrum(var_eig.CENTERED_EIGENVALUES, 0.4);

Location :

SetNbMaximumIterations, GetNbMaximumIterations

Syntax :

  void SetNbMaximumIterations(int);
  int GetNbMaximumIterations();

With those methods, you can modify and retrieve the maximal number of iterations completed by the iterative algorithm. If the iterative algorithm spends more iterations, it is stopped. The default value is equal to 1000.

Example :

// declaration of the eigenproblem
SparseEigenProblem<double, Matrix<double, General, ArrayRowSparse> > var_eig;
Vector<double> eigen_values, eigen_imag;
Matrix<double> eigen_vec;

// you can set some parameters
var_eig.SetNbAskedEigenvalues(5);
var_eig.SetStoppingCriterion(1e-12);
var_eig.SetNbMaximumIterations(10000);
var_eig.SetTypeSpectrum(var_eig.CENTERED_EIGENVALUES, 0.4);

Location :

GetNbMatrixVectorProducts

Syntax :

  int GetNbMatrixVectorProducts();

This method returns the number of matrix vector products performed by the eigenvalue solver.

Example :

// declaration of the eigenproblem
VirtualEigenProblem<double> var_eig;
Vector<double> eigen_values, eigen_imag;
Matrix<double> eigen_vec;

// you can set some parameters
var_eig.SetNbAskedEigenvalues(5);
var_eig.SetTypeSpectrum(var_eig.LARGE_EIGENVALUES, 0.0);

var_eig.InitMatrix(K);

// research of eigenvalues
GetEigenvaluesEigenvectors(var_eig, eigen_values, eigen_imag, eigen_vec);

// then displaying the number of matrix vector products performed :
cout << "Number of times K x has been done : " << var_eig.GetNbMatrixVectorProducts() << endl;

Location :

GetNbLinearSolves

Syntax :

  int GetNbLinearSolves();

This method returns the number of linear solves performed by the eigenvalue solver.

Example :

// declaration of the eigenproblem
VirtualEigenProblem<double> var_eig;
Vector<double> eigen_values, eigen_imag;
Matrix<double> eigen_vec;

// you can set some parameters
var_eig.SetNbAskedEigenvalues(5);
var_eig.SetTypeSpectrum(var_eig.CENTERED_EIGENVALUES, 1.0);

var_eig.InitMatrix(K);

// research of eigenvalues
GetEigenvaluesEigenvectors(var_eig, eigen_values, eigen_imag, eigen_vec);

// then displaying the number of linear solves performed :
cout << "Number of linear solves (Ax = b) computed : " << var_eig.GetNbLinearSolves() << endl;

Location :

GetNbArnoldiVectors, SetNbArnoldiVectors

Syntax :

  int GetNbArnoldiVectors();
  void SetNbArnoldiVectors(int n);

You can modify and retrieve the number of Arnoldi vectors used by the iterative algorithm. It can be a good idea to increase the number of Arnoldi vectors in order to improve the convergence. By default, the number of Arnoldi vectors is set to 2 nev + 2 where nev is the number of asked eigenvalues.

Example :

// declaration of the eigenproblem
VirtualEigenProblem<double> var_eig;
Vector<double> eigen_values, eigen_imag;
Matrix<double> eigen_vec;

// you can set some parameters
var_eig.SetNbAskedEigenvalues(5);
var_eig.SetTypeSpectrum(var_eig.LARGE_EIGENVALUES, 0.0);

// and require a larger number of Arnoldi vectors
var_eig.SetNbArnoldiVectors(20);

var_eig.InitMatrix(K);

// research of eigenvalues
GetEigenvaluesEigenvectors(var_eig, eigen_values, eigen_imag, eigen_vec);

cout << "Number of Arnoldi vectors used by simulation" << var_eig.GetNbArnoldiVectors() << endl;

Location :

GetM, GetGlobalM, GetN

Syntax :

  int GetM();
  int GetN();
  int GetGlobalM();

GetM and GetN return the local number of rows of the eigenproblem. GetGlobalM() returns the global number of rows, i.e. the sum of local numbers for the involved processors.

Location :

GetPrintLevel

Syntax :

  int GetPrintLevel();
  void SetPrintLevel(int n);

By increasing print level, more messages should be displayed on the standard output.

Location :

IncrementProdMatVect

Syntax :

void IncrementProdMatVect();

This method is used internally to increment the variable containing the number of matrix vector products. It should not be called by the user.

Location :

IncrementLinearSolves

Syntax :

void IncrementLinearSolves();

This method is used internally to increment the variable containing the number of linear solves. It should not be called by the user.

Location :

PrintErrorInit

Syntax :

void PrintErrorInit();

This method is used internally to display an error message and stop the program if InitMatrix has not been called.

Location :

IsSymmetricProblem

Syntax :

bool IsSymmetricProblem();

This method returns true if the eigenvalue problem involves symmetric matrices. This method is virtual such that you can overload it to specify if your own eigenproblem is symmetric or not.

Location :

IsHermitianProblem

Syntax :

bool IsHermitianProblem();

This method returns true if the eigenvalue problem involves a symmetric mass matrix and an hermitian stiffness matrix (for a linear eigenproblem). This method is virtual such that you can overload it to specify if your own eigenproblem is hermitian or not.

Location :

FactorizeDiagonalMass

Syntax :

void FactorizeDiagonalMass();

This method computes the square root of diagonal mass matrix. It is a virtual method that can be overloaded if you define your own eigenproblem.

Example:

template<class T, class MatStiff, class MatMass>
class MyOwnEigenProblem : public VirtualEigenProblem<T, typename MatStiff::entry_type, typename MatMass::entry_type>
{
  protected:
    // you can add members to handle diagonal mass matrix
    // for example, you could consider a diagonal with only
    // two coefficient D = [alpha I, 0; 0, beta I]
    double alpha, beta, sqrt_alpha, sqrt_beta;

  public :

  // computation of diagonal mass matrix
  void ComputeDiagonalMass()
  {
     // you may use Mh
     if (this->Mh != NULL)
       {
         MatMass& M = dynamic_cast<MatMass&>(*this->Mh);
         alpha = M(0, 0);
         beta = M(1, 1);
       }
     else
       {
         // or not
         alpha = 1.5; beta = 3.0;
       }
   }
                                                                              }

  // computation of sqrt(D)
  void FactorizeDiagonalMass()
  {
     sqrt_alpha = sqrt(alpha);
     sqrt_beta = sqrt(beta);
  }

  // sqrt(D) can be written in an output vector
  void GetSqrtDiagonal(Vector<T>& D)
  {
    D.Reallocate(this->n_);
    for (int i = 0; i < this->n_; i += 2)
      {
        D(i) = sqrt_alpha;
        D(i+1) = sqrt_beta;
      }
  }

  // application of sqrt(D)
  void MltSqrtDiagonalMass(Vector<double>& X)
  {
    for (int i = 0; i < X.GetM(); i+=2)
      {
        X(i) *= sqrt_alpha;
        X(i+1) *= sqrt_beta;
      }
  }

  // application of 1/sqrt(D)
  void MltInvSqrtDiagonalMass(Vector<double>& X)
  {
    for (int i = 0; i < X.GetM(); i+=2)
      {
         X(i) /= sqrt_alpha;
         X(i+1) /= sqrt_beta;
      }
   }

};


int main()
{

// then you can use this class :
MyOwnEigenProblem<double, MyStiffMatrix, MyMassMatrix> var_eig;
Vector<double> lambda, lambda_imag;
Matrix<double> eigen_vec;

// set some parameters
var_eig.SetNbAskedEigenvalues(5);
// specify that mass matrix is actually diagonal
var_eig.SetDiagonalMass();

// initialized pointers to matrices
var_eig.InitMatrix(K, M);

// then call computation of eigenvalues
GetEigenvaluesEigenvectors(var_eig, lambda, lambda_imag, eigen_vec);

}

Location :

MltSqrtDiagonalMass, MltInvSqrtDiagonalMass

Syntax :

void MltSqrtDiagonalMass(Vector&);
void MltInvSqrtDiagonalMass(Vector&);

These methods apply M^1/2 or M^-1/2 to a given vector. It is a virtual method that can be overloaded if you define your own eigenproblem. You can look at the example given in the method FactorizeDiagonalMass.

Location :

GetSqrtDiagonal

Syntax :

void GetSqrtDiagonal(Vector& D);

This method fills sqrt(M) in the vector D when M is the mass matrix and has been declared as diagonal. It is a virtual method that can be overloaded if you define your own eigenproblem. You can look at the example given in the method FactorizeDiagonalMass.

Location :

ComputeDiagonalMass

Syntax :

void ComputeDiagonalMass();

This method computes the diagonal of mass matrix. It is used internally, it should not be called directly by the user. However, if you derive a class for your own eigenproblem, and you want to change how diagonal mass matrices are handled, you can overload this method. You can look at the example given in the method FactorizeDiagonalMass.

Location :

ComputeMassForCholesky

Syntax :

void ComputeMassForCholesky();

This method is assumed to compute the mass matrix when a Cholesky factorisation is required. It can differ from the computation of the mass matrix when a matrix vector product is required, since that in that latter case, the matrix may be not stored. This method is used internally and should not be called directly by the user. However, if you need to change how Cholesky factorisation is handled, you can overload this function as shown in the example of member function FactorizeCholeskyMass.

Location :

ComputeMassMatrix

Syntax :

void ComputeMassMatrix();

This method is assumed to compute the mass matrix when a matrix vector product is required. This method is called internally, and the user should not call it directly. However, if you have your own storage and matrix-vector product with mass matrix, you can overload this function as shown in the example of member function MltMass.

Location :

MltMass

Syntax :

void MltMass(const Vector& X, Vector& Y);
void MltMass(const SeldonTranspose&, const Vector& X, Vector& Y);

This method is assumed to compute Y = M X, where M is the mass matrix. In the second syntax, you can multiply by the transpose of M or not. This method is called internally, and the user should not call it directly. However, if you have your own storage and matrix-vector product with mass matrix, you can overload this function as shown below :

Example :

template<class T, class MatStiff, class MatMass>
class MyOwnEigenProblem : public VirtualEigenProblem<T, typename MatStiff::entry_type, typename MatMass::entry_type>
{
  protected:
  // you can add members to handle mass matrix
  double coef;

  public :
  
  // computation of mass matrix
  void ComputeMassMatrix()
  {
    // here you compute all you need to perform matrix-vector products
    // with the mass matrix
    coef = 2.5;
  }  

  // application of mass matrix (tridiagonal here)
  void MltMass(const Vector<T>& X, Vector<T>& Y)
  {
    int n = X.GetM();
    Y(0) = 4.0*X(0) + X(1);
    Y(n-1) = 4.0*X(n-1) + X(n-2);
    for (int i = 1; i < n-1; i++)
      Y(i) = 4.0*X(i) + X(i-1) + X(i+1);
    
    Mlt(coef, Y);
  }

  // application of mass matrix (tridiagonal here)
  void MltMass(const SeldonTranspose& trans, const Vector<T>& X, Vector<T>& Y)
  {
  // symmetric matrix -> trans is not used
    int n = X.GetM();
    Y(0) = 4.0*X(0) + X(1);
    Y(n-1) = 4.0*X(n-1) + X(n-2);
    for (int i = 1; i < n-1; i++)
      Y(i) = 4.0*X(i) + X(i-1) + X(i+1);
    
    Mlt(coef, Y);
  }

};

int main()
{

// then you can use this class :
MyOwnEigenProblem<double, MyStiffMatrix, MyMassMatrix> var_eig;
Vector<double> lambda, lambda_imag;
Matrix<double> eigen_vec;

// set some parameters
var_eig.SetNbAskedEigenvalues(5);

// initializing pointers to matrices
var_eig.InitMatrix(K, M);

// then calling computation of eigenvalues
GetEigenvaluesEigenvectors(var_eig, lambda, lambda_imag, eigen_vec);

}

Location :

ComputeStiffnessMatrix

Syntax :

void ComputeStiffnessMatrix();
void ComputeStiffnessMatrix(const T& a, const T& b);

This method is assumed to compute the stiffness matrix when a matrix vector product is required. This method is called internally, and the user should not call it directly. However, if you have your own storage and matrix-vector product with stiffness matrix, you can overload this function as shown in the example of member function MltStiffness. One other option is to write your own class for stiffness matrix inheriting from the class VirtualMatrix and overload the methods MltVector and MltAddVector as shown in the section "Advanced Use", then you can use VirtualEigenProblem class. Nevertheless, this class is not appropriate for shift-invert mode (only regular mode can be used).

Location :

MltStiffness

Syntax :

void MltStiffness(const Vector& X, Vector& Y);
void MltStiffness(const T& a, const T& b, const Vector& X, Vector& Y);
void MltStiffness(const SeldonTranspose& trans, const Vector& X, Vector& Y);

This method is assumed to compute Y = K X (or Y = (a M + b K) X ) or Y = K^T x, where M is the mass matrix and K the stiffness matrix. This method is called internally, and the user should not call it directly. We assume here that either we need to compute Y = K X, either Y = (a M + b K) X (depending on the computational mode), but we do not need to perform both operations. However, if you have your own storage and matrix-vector product with stiffness matrix, you can overload this function as shown below :

Example :

template<class T, class MatStiff, class MatMass>
class MyOwnEigenProblem : public VirtualEigenProblem<T, typename MatStiff::entry_type, typename MatMass::entry_type>
{
  protected:
  // you can add members to handle stiffness matrix
  double alpha, beta, gamma;

  public :
  
  // computation of stiffness matrix
  void ComputeStiffnessMatrix()
  {
    // you compute here all you need 
    alpha = 1.4;
    beta = 2.5;
    gamma = 0.8;
  }

  // computation of a M + b K
  void ComputeStiffnessMatrix(const T& a, const T& b)
  {
    // we consider that M = I
    alpha = 1.4*b;
    beta = a + 2.5*b;
    gamma = 0.8*b;
  }
  

  // application of stiffness matrix (tridiagonal here)
  void MltStiffness(const Vector<T>& X, Vector<T>& Y)
  {
    int n = X.GetM();
    Y(0) = beta*X(0) + gamma*X(1);
    Y(n-1) = beta*X(n-1) + alpha*X(n-2);
    for (int i = 1; i < n-1; i++)
      Y(i) = beta*X(i) + alpha*X(i-1) + gamma*X(i+1);
    
  }


  // computation of Y = (a M + b K) X
  void MltStiffness(const T& a, const T& b, const Vector<T>& X, Vector<T>& Y)
  {
     // we use the fact that either we have to compute Y = K X 
     // either Y = (a M + b K) X, but not both
     MltStiffness(X, Y);
  }

  // application of stiffness matrix (tridiagonal here)
  void MltStiffness(const SeldonTranspose& trans, const Vector<T>& X, Vector<T>& Y)
  {
    if (trans.NoTrans())
      return MltStiffness(X, Y);

    // multiplication with transpose : alpha and gamma are exchanged
    int n = X.GetM();
    Y(0) = beta*X(0) + alpha*X(1);
    Y(n-1) = beta*X(n-1) + gamma*X(n-2);
    for (int i = 1; i < n-1; i++)
      Y(i) = beta*X(i) + gamma*X(i-1) + alpha*X(i+1);
    
  }

};

int main()
{

// then you can use this class :
MyOwnEigenProblem<double, MyStiffMatrix, MyMassMatrix> var_eig;
Vector<double> lambda, lambda_imag;
Matrix<double> eigen_vec;

// set some parameters
var_eig.SetNbAskedEigenvalues(5);

// initializing pointers to matrices
var_eig.InitMatrix(K, M);

// then calling computation of eigenvalues
GetEigenvaluesEigenvectors(var_eig, lambda, lambda_imag, eigen_vec);

}

Location :

ComputeAndFactorizeStiffnessMatrix

Syntax :

void ComputeAndFactorizeStiffnessMatrix(const T& a, const T& b, int real_p = COMPLEX_PART);
void ComputeAndFactorizeStiffnessMatrix(const complex<T>& a, const complex<T>& b, int real_p = COMPLEX_PART);

This method is assumed to compute matrix a M + b K and factorize this matrix. If this matrix is solved by an iterative algorithm, the factorization can consist of constructing a preconditioning. This method is called internally, and the user should not call it directly. However, if you have your own procedure to solve linear system (a M + b K) y = x , you can overload this function as shown in the example of member function ComputeSolution. For real unsymmetric matrices, coefficients a and b can be complex and we require complex part, real part or imaginary part of the factorization of a M + b K.

Location :

ComputeSolution

Syntax :

void ComputeSolution(const Vector& X, Vector& Y);
void ComputeSolution(const SeldonTranspose&, const Vector& X, Vector& Y);

This method is assumed to solve (a M + b K) Y = X, where M is the mass matrix and K the stiffness matrix. Coefficients a and b have been given when calling the function ComputeAndFactorizeStiffnessMatrix. This method is called internally, and the user should not call it directly. However, if you have your own procedure to solver the linear system (a M + b K) y = x , you can overload this function as shown below :

Example :

template<class T, class MatStiff, class MatMass>
class MyOwnEigenProblem : public VirtualEigenProblem<T, typename MatStiff::entry_type, typename MatMass::entry_type>,
       public VirtualMatrix<T>, public Preconditioner_Base<T>
{
  protected:
  // you can add members to handle factorization
  double alpha, beta;
  Vector<double> precond;

  public :
  
  // computation of  a M + b K and factorization
  void ComputeAndFactorizeStiffnessMatrix(const T& a, const T& b, int real_p = COMPLEX_PART)
  {
    const MatMass& M = dynamic_cast<MatMass&>(*this->Mh);
    const MatStiff& K = dynamic_cast<MatStiff&>(*this->Kh);
    alpha = a;
    beta = b;
    // iterative process, constructing preconditioning
    // example of diagonal preconditioning
    precond.Reallocate(this->n_);
    for (int i = 0; i < this->n_; i++)
      precond(i) = 1.0 / (a*M(i, i) + b*K(i, i));
  }  

  // computing solution of (a M + b K) y = x
  void ComputeSolution(const Vector<T>& X, Vector<T>& Y)
  {
    // use of conjugate gradient to solve the linear system
    Iteration<double> iter(10000, 1e-12); 
    Cg(*this, Y, X, *this, iter);
  }

  void MltVector(const Vector<T>& x, Vector<T>& y)
  {  
    Mlt(*this->Mh, x, y);
    MltAdd(beta, *this->Kh, x, alpha, y);
  }

  void MltAddVector(const T& a, const Vector<T>& x, const T& b, Vector<T>& y)
  {
    MltAdd(a*alpha, *this->Mh, x, b, y);
    MltAdd(a*beta, *this->Kh, x, T(1), y);
  }

  // application of preconditioning
  void Solve(const VirtualMatrix<T>& A, const Vector<T>& r, Vector<T>& z)
  {
    for (int i = 0; i < this->n_; i++)
      z(i) = r(i)*precond(i);
  }

  void ComputeSolution(const SeldonTranspose& trans, const Vector<T>& X, Vector<T>& Y)
  {
    // system is assumed to be real symmetric here
    // otherwise you have compute the solution of the transpose system or conjugate transpose
    ComputeSolution(X, Y);
  }

};


int main()
{

// then you can use this class :
MyOwnEigenProblem<double, MyStiffMatrix, MyMassMatrix> var_eig;
Vector<double> lambda, lambda_imag;
Matrix<double> eigen_vec;

// set some parameters
var_eig.SetNbAskedEigenvalues(5);
var_eig.SetTypeSpectrum(var_eig.CENTERED_EIGENVALUES, 0.1);
var_eig.SetComputationalMode(var_eig.SHIFTED_MODE);

// initializing pointers to matrices
var_eig.InitMatrix(K, M);

// then calling computation of eigenvalues
GetEigenvaluesEigenvectors(var_eig, lambda, lambda_imag, eigen_vec);

}

Location :

FactorizeCholeskyMass

Syntax :

void FactorizeCholeskyMass();

This method computes a Cholesky factorization of mass matrix. It is used internally, and should not be called directly. However, if you derive a class for your own eigenproblem, and you want to change how Cholesky factorisation is handled, you can overload this method. Then, look at the example detailed in MltCholeskyMass.

Location :

SolveCholeskyMass, MltCholeskyMass

Syntax :

void MltCholeskyMass(SeldonTranspose, Vector& x);
void SolveCholeskyMass(SeldonTranspose, Vector& x);

These methods apply L, L^T, L^-1 or L^-T to a given vector. It is used internally, and should not be called directly. However, if you derive a class for your own eigenproblem, and you want to change how Cholesky factorization is handled, you can overload this method as shown in the example below.

Example :

template<class T, class MatStiff, class MatMass>
class MyOwnEigenProblem : public VirtualEigenProblem<T, typename MatStiff::entry_type, typename MatMass::entry_type>
{
  protected:
  // you can add members to handle Cholesky mass matrix
  // for example, you could consider a diagonal with only
  // two coefficient D = [alpha I, 0; 0, beta I]
  double alpha, beta, sqrt_alpha, sqrt_beta;

  public :
  
  // computation of mass in preparation of Cholesky factorization
  void ComputeMassForCholesky()
  {
  }  

  // computation of Cholesky factor
  void FactorizeCholeskyMass()
  {
    // TO DO
  }
  
  // application of L or L^T
  void MltCholesyMass(const SeldonTranspose& TransA, Vector<double>& X)
  {
    // TO DO
  }

  // application of L^-1 or L^-T
  void SolveCholeskyMass(const SeldonTranspose& TransA, Vector<double>& X)
  {
    // TO DO
  }
};

int main()
{

// then you can use this class :
MyOwnEigenProblem<double, MyStiffMatrix, MyMassMatrix> var_eig;
Vector<double> lambda, lambda_imag;
Matrix<double> eigen_vec;

// set some parameters
var_eig.SetNbAskedEigenvalues(5);
// specify that you want to use Cholesky factors of mass matrix
var_eig.SetCholeskyFactoForMass();

// initialized pointers to matrices
var_eig.InitMatrix(K, M);

// then call computation of eigenvalues
GetEigenvaluesEigenvectors(var_eig, lambda, lambda_imag, eigen_vec);

}

Location :

Clear

Syntax :

void Clear();

This method releases memory used for the factorization of mass matrix and/or stiffness matrix.

Location :

VirtualEigenvalueSolver.cxx
PolynomialEigenvalueSolver.cxx
NonLinearEigenvalueSolver.cxx

GetEigenvaluesEigenvectors

Syntax :

void GetEigenvaluesEigenvectors(EigenProblem_Base& var, Vector& lambda, Vector& lambda_imag, Matrix& eigen_vec, int type_solver = 0);
void GetEigenvaluesEigenvectors(PolynomialEigenProblem_Base& var, Vector& lambda, Vector& lambda_imag, Matrix& eigen_vec, int type_solver = 0);
void GetEigenvaluesEigenvectors(NonLinearEigenProblem_Base& var, Vector& lambda, Vector& lambda_imag, Matrix& eigen_vec, int type_solver = 0);

The first function compute some eigenvalues and eigenvectors of a matrix (dense, sparse or defined by the user) : K x = λ x, or a generalized eigenvalue problem : K x = λ M x, K being called stiffness matrix and M mass matrix. For dense matrices, you can use DenseEigenProblem. For sparse matrices (Seldon format), you can use SparseEigenProblem, it will use direct solvers interfaced by Seldon to solve intermediary linear systems (if needed). For matrices defined by the user, you can evaluate largest eigenvalues with only the matrix vector product thanks to the class VirtualEigenProblem. Eigenvectors are placed in a dense matrix, each column being associated with the corresponding eigenvalue. In the case of real unsymmetric eigenvalue problems, the imaginary part of eigenvalues is placed in vector lambda_imag, and the column i of eigen_vec represents the real part u of the eigenvector, the column i+1 stores the imaginary part v (eigenvectors are then equal to u + i v and u - i v ). An optional last argument can be provided in order to select the eigenvalue solver you prefer. By default it will select (by order of preference) Arpack, Feast or Slepc. The second function handles polynomial eigenproblems while the last one handles non-linear eigenproblems.

// declaration of the eigenvalue problem
DenseEigenProblem<double, double, Symmetric, RowSymPacked> var_eig;

// setting parameters of eigenproblem
var_eig.SetStoppingCriterion(1e-12);
var_eig.SetNbAskedEigenvalues(10);
// you can ask largest eigenvalues of M^-1 K, smallest 
// or eigenvalues closest to a shift sigma
var_eig.SetTypeSpectrum(var_eig.LARGE_EIGENVALUES, 0);

// then you select the computational mode
var_eig.SetComputationalMode(var_eig.REGULAR_MODE);

// giving matrices K and M to the eigenvalue problem
int n = 20;
Matrix<double, Symmetric, RowSymPacked> K(n, n), M(n, n);
K.FillRand(); M.FillRand();
var_eig.InitMatrix(K, M);

// declaring arrays that will contains eigenvalues and eigenvectors
Vector<double> lambda, lambda_imag;
Matrix<double> eigen_vec;
GetEigenvaluesEigenvectors(var_eig, lambda, lambda_imag, eigen_vec);

// you can also impose a solver as a last argument
GetEigenvaluesEigenvectors(var_eig, lambda, lambda_imag, eigen_vec, TypeEigenvalueSolver::FEAST);

Location :

Init

Syntax :

void Init(int n, int nev, int ncv, int maxit, double tolerance, string solver_type,
          int mode, string which, char bmat, char HowMny, bool arpack_with_verbose)

This function inits parameters related to eigenvalue problem (see Arpack documentation).

Location :

CheckParameter

Syntax :

void CheckParameter();

This function checks if the parameters - given when calling Init - are correct.

Location :

Clear

Syntax :

void Clear()

This function clears memory used by internal arrays.

Location :

Allocate

Syntax :

void Allocate()

This function allocates memory used by internal arrays.

Location :

Deallocate

Syntax :

void Deallocate()

This function deallocates memory used by internal arrays.

Location :

SetArpackVerbose

Syntax :

void SetArpackVerbose()

This function sets verbose mode.

Location :

ClearArpackVerbose

Syntax :

void ClearArpackVerbose()

This function returns to silent mode.

Location :

GetFirstWorkVector

Syntax :

T* GetFirstWorkVector()

This function returns a pointer to the first work vector.

Location :

GetSecondWorkVector

Syntax :

T* GetSecondWorkVector()

This function returns a pointer to the second work vector.

Location :

GetEigenValue

Syntax :

T GetEigenValue(int i)

This function returns the i-th eigenvalue.

Location :

GetEigenVector

Syntax :

T* GetEigenVector(int i)

This function returns the i-th eigenvector.

Location :

GetReverseCommunicationFlag, SetReverseCommunicationFlag

Syntax :

void SetReverseCommunicationFlag(int flag)
int GetReverseCommunicationFlag()

This function gives access to the reverse-communication flag (ido).

Location :

GetInfoFlag, SetInfoFlag

Syntax :

void SetInfoFlag(int flag)
int GetInfoFlag()

This function gives access to the information flag (info).

Location :

GetConvergedNumber

Syntax :

int GetConvergedNumber() const

This function returns the number of converged eigenvalues.

Location :

Continue

Syntax :

bool Continue()

This function completes one step of the algorithm, it should return true if the iterative algorithm has ended.

Location :

Finish

Syntax :

bool Finish()

This function completes the computation of eigenvalues and eigenvectors, it is called once the iterative algorithm has ended.

Location :

GetEigensolverType, SetEigensolverType, GetEigensolverChar

Syntax :

  void SetEigensolverType(int type);
  int GetEigensolverType();
  const char* GetEigensolverChar();

SetEigensolverType sets the type of eigensolver to use in Slepc, whereas GetEigensolverType returns this number. You can choose between POWER, SUBSPACE, ARNOLDI, LANCZOS, KRYLOVSCHUR, GD, JD, RQCG, LOBPCG, CISS, LAPACK, ARPACK, BLZPACK, TRLAN, BLOPEX, PRIMME or FEAST. This list corresponds to the list provided by Slepc. The method GetEigensolverChar returns the solver type as a string.

Example :

  SparseEigenProblem<T, MatrixK, MatrixM> var_eig;

  var_eig.SetStoppingCriterion(1e-14);
  var_eig.SetNbAskedEigenvalues(nb_eigenval);
  var_eig.SetComputationalMode(var_eig.REGULAR_MODE);

  var_eig.InitMatrix(K, M);

  SlepcParam& param = var_eig.GetSlepcParameters();
  param.SetEigensolverType(param.ARNOLDI);

  // Large eigenvalues
  var_eig.SetTypeSpectrum(var_eig.LARGE_EIGENVALUES, 0, var_eig.SORTED_MODULUS);
  
  GetEigenvaluesEigenvectors(var_eig, lambda, lambda_imag, eigen_vec, TypeEigenvalueSolver::SLEPC);
  DISP(lambda); DISP(lambda_imag);

Location :

GetBlockSize, SetBlockSize

Syntax :

  void SetBlockSize(int n);
  int GetBlockSize();

SetBlockSize sets the block size to use in Slepc, whereas GetBlockSize returns this number. This block size is used when Seldon calls the following Slepc functions : EPSBLOPEXSetBlockSize, EPSCISSSetSizes, EPSGDSetBlockSize, EPSJDSetBlockSize, EPSLOBPCGSetBlockSize, EPSPRIMMESetBlockSize.

Example :

  SparseEigenProblem<T, MatrixK, MatrixM> var_eig;

  var_eig.SetStoppingCriterion(1e-14);
  var_eig.SetNbAskedEigenvalues(nb_eigenval);
  var_eig.SetComputationalMode(var_eig.REGULAR_MODE);

  var_eig.InitMatrix(K, M);

  SlepcParam& param = var_eig.GetSlepcParameters();
  param.SetBlockSize(20);

  // Large eigenvalues
  var_eig.SetTypeSpectrum(var_eig.LARGE_EIGENVALUES, 0, var_eig.SORTED_MODULUS);
  
  GetEigenvaluesEigenvectors(var_eig, lambda, lambda_imag, eigen_vec, TypeEigenvalueSolver::SLEPC);
  DISP(lambda); DISP(lambda_imag);

Location :

GetMaximumBlockSize, SetMaximumBlockSize

Syntax :

  void SetMaximumBlockSize(int n);
  int GetMaximumBlockSize();

SetMaximumBlockSize sets the maximum block size to use in Slepc, whereas GetMaximumBlockSize returns this number. This maximum block size is used when Seldon calls the following Slepc function : EPSCISSSetSizes.

Example :

  SparseEigenProblem<T, MatrixK, MatrixM> var_eig;

  var_eig.SetStoppingCriterion(1e-14);
  var_eig.SetNbAskedEigenvalues(nb_eigenval);
  var_eig.SetComputationalMode(var_eig.REGULAR_MODE);

  var_eig.InitMatrix(K, M);

  SlepcParam& param = var_eig.GetSlepcParameters();
  param.SetMaximumBlockSize(20);

  // Large eigenvalues
  var_eig.SetTypeSpectrum(var_eig.LARGE_EIGENVALUES, 0, var_eig.SORTED_MODULUS);
  
  GetEigenvaluesEigenvectors(var_eig, lambda, lambda_imag, eigen_vec, TypeEigenvalueSolver::SLEPC);
  DISP(lambda); DISP(lambda_imag);

Location :

GetNumberOfSteps, SetNumberOfSteps

Syntax :

  void SetNumberOfSteps(int n);
  int GetNumberOfSteps();

SetNumberOfSteps sets the number of steps to use in Slepc, whereas GetNumberOfSteps returns this number. This number of steps is used when Seldon calls the following Slepc functions : EPSBlzpackSetNSteps, EPSRQCGSetReset.

Example :

  SparseEigenProblem<T, MatrixK, MatrixM> var_eig;

  var_eig.SetStoppingCriterion(1e-14);
  var_eig.SetNbAskedEigenvalues(nb_eigenval);
  var_eig.SetComputationalMode(var_eig.REGULAR_MODE);

  var_eig.InitMatrix(K, M);

  SlepcParam& param = var_eig.GetSlepcParameters();
  param.SetNumberOfSteps(5);
  
  // Large eigenvalues
  var_eig.SetTypeSpectrum(var_eig.LARGE_EIGENVALUES, 0, var_eig.SORTED_MODULUS);
  
  GetEigenvaluesEigenvectors(var_eig, lambda, lambda_imag, eigen_vec, TypeEigenvalueSolver::SLEPC);
  DISP(lambda); DISP(lambda_imag);

Location :

GetExtractionType, SetExtractionType

Syntax :

  void SetExtractionType(int n);
  int GetExtractionType();

SetExtractionType sets the extraction method to use in the solver CISS of Slepc, whereas GetExtractionType returns this method. You can choose between EXTRACT_RITZ and EXTRACT_HANKEL (in class SlepcParam). This method used when Seldon calls the following Slepc function : EPSCISSSetExtraction.

Example :

  SparseEigenProblem<T, MatrixK, MatrixM> var_eig;

  var_eig.SetStoppingCriterion(1e-14);
  var_eig.SetNbAskedEigenvalues(nb_eigenval);
  var_eig.SetComputationalMode(var_eig.REGULAR_MODE);

  var_eig.InitMatrix(K, M);

  SlepcParam& param = var_eig.GetSlepcParameters();
  param.SetExtractionType(param.EXTRACT_RITZ);
  
  // Large eigenvalues
  var_eig.SetTypeSpectrum(var_eig.LARGE_EIGENVALUES, 0, var_eig.SORTED_MODULUS);
  
  GetEigenvaluesEigenvectors(var_eig, lambda, lambda_imag, eigen_vec, TypeEigenvalueSolver::SLEPC);
  DISP(lambda); DISP(lambda_imag);

Location :

GetQuadratureRuleType, SetQuadratureRuleType

Syntax :

  void SetQuadratureRuleType(int n);
  int GetQuadratureRuleType();

SetQuadratureRuleType sets the quadrature rule to use in the solver CISS of Slepc, whereas GetQuadratureRuleType returns this rule. You can choose between QUADRULE_TRAPEZE and QUADRULE_CHEBY (in class SlepcParam). This method is used when Seldon calls the following Slepc function : EPSCISSSetQuadRule.

Example :

  SparseEigenProblem<T, MatrixK, MatrixM> var_eig;

  var_eig.SetStoppingCriterion(1e-14);
  var_eig.SetNbAskedEigenvalues(nb_eigenval);
  var_eig.SetComputationalMode(var_eig.REGULAR_MODE);

  var_eig.InitMatrix(K, M);

  SlepcParam& param = var_eig.GetSlepcParameters();
  param.SetQuadratureRuleType(param.QUADRULE_CHEBY);
  
  // Large eigenvalues
  var_eig.SetTypeSpectrum(var_eig.LARGE_EIGENVALUES, 0, var_eig.SORTED_MODULUS);
  
  GetEigenvaluesEigenvectors(var_eig, lambda, lambda_imag, eigen_vec, TypeEigenvalueSolver::SLEPC);
  DISP(lambda); DISP(lambda_imag);

Location :

GetInnerSteps, SetInnerSteps

Syntax :

  void SetInnerSteps(int n);
  int GetInnerSteps();

SetInnerSteps sets the number of inner steps to use in the solver CISS of Slepc, whereas GetInnerSteps returns this number. This method is used when Seldon calls the following Slepc function : EPSCISSSetRefinement.

Example :

  SparseEigenProblem<T, MatrixK, MatrixM> var_eig;

  var_eig.SetStoppingCriterion(1e-14);
  var_eig.SetNbAskedEigenvalues(nb_eigenval);
  var_eig.SetComputationalMode(var_eig.REGULAR_MODE);

  var_eig.InitMatrix(K, M);

  SlepcParam& param = var_eig.GetSlepcParameters();
  param.SetInnerSteps(10);
  
  // Large eigenvalues
  var_eig.SetTypeSpectrum(var_eig.LARGE_EIGENVALUES, 0, var_eig.SORTED_MODULUS);
  
  GetEigenvaluesEigenvectors(var_eig, lambda, lambda_imag, eigen_vec, TypeEigenvalueSolver::SLEPC);
  DISP(lambda); DISP(lambda_imag);

Location :

GetOuterSteps, SetOuterSteps

Syntax :

  void SetOuterSteps(int n);
  int GetOuterSteps();

SetOuterSteps sets the number of outer steps to use in the solver CISS of Slepc, whereas GetOuterSteps returns this number. This method is used when Seldon calls the following Slepc function : EPSCISSSetRefinement.

Example :

  SparseEigenProblem<T, MatrixK, MatrixM> var_eig;

  var_eig.SetStoppingCriterion(1e-14);
  var_eig.SetNbAskedEigenvalues(nb_eigenval);
  var_eig.SetComputationalMode(var_eig.REGULAR_MODE);

  var_eig.InitMatrix(K, M);

  SlepcParam& param = var_eig.GetSlepcParameters();
  param.SetOuterSteps(10);
  
  // Large eigenvalues
  var_eig.SetTypeSpectrum(var_eig.LARGE_EIGENVALUES, 0, var_eig.SORTED_MODULUS);
  
  GetEigenvaluesEigenvectors(var_eig, lambda, lambda_imag, eigen_vec, TypeEigenvalueSolver::SLEPC);
  DISP(lambda); DISP(lambda_imag);

Location :

GetNumberIntegrationPoints, SetNumberIntegrationPoints

Syntax :

  void SetNumberIntegrationPoints(int n);
  int GetNumberIntegrationPoints();

SetNumberIntegrationPoints sets the number of integration points to use in the solver CISS (or Feast interface) of Slepc, whereas GetNumberIntegrationPoints returns this number. This method is used when Seldon calls the following Slepc functions : EPSCISSSetSizes, EPSFEASTSetNumPoints.

Example :

  SparseEigenProblem<T, MatrixK, MatrixM> var_eig;

  var_eig.SetStoppingCriterion(1e-14);
  var_eig.SetNbAskedEigenvalues(nb_eigenval);
  var_eig.SetComputationalMode(var_eig.REGULAR_MODE);

  var_eig.InitMatrix(K, M);

  SlepcParam& param = var_eig.GetSlepcParameters();
  param.SetNumberIntegrationPoints(10);
  
  // Large eigenvalues
  var_eig.SetTypeSpectrum(var_eig.LARGE_EIGENVALUES, 0, var_eig.SORTED_MODULUS);
  
  GetEigenvaluesEigenvectors(var_eig, lambda, lambda_imag, eigen_vec, TypeEigenvalueSolver::SLEPC);
  DISP(lambda); DISP(lambda_imag);

Location :

GetMomentSize, SetMomentSize

Syntax :

  void SetMomentSize(int n);
  int GetMomentSize();

SetMomentSize sets the moment size to use in the solver CISS of Slepc, whereas GetMomentSize returns this number. This method is used when Seldon calls the following Slepc function : EPSCISSSetSizes.

Example :

  SparseEigenProblem<T, MatrixK, MatrixM> var_eig;

  var_eig.SetStoppingCriterion(1e-14);
  var_eig.SetNbAskedEigenvalues(nb_eigenval);
  var_eig.SetComputationalMode(var_eig.REGULAR_MODE);

  var_eig.InitMatrix(K, M);

  SlepcParam& param = var_eig.GetSlepcParameters();
  param.SetMomentSize(10);
  
  // Large eigenvalues
  var_eig.SetTypeSpectrum(var_eig.LARGE_EIGENVALUES, 0, var_eig.SORTED_MODULUS);
  
  GetEigenvaluesEigenvectors(var_eig, lambda, lambda_imag, eigen_vec, TypeEigenvalueSolver::SLEPC);
  DISP(lambda); DISP(lambda_imag);

Location :

GetNumberPartitions, SetNumberPartitions

Syntax :

  void SetNumberPartitions(int n);
  int GetNumberPartitions();

SetNumberPartitions sets the number of partitions to use in the solver CISS of Slepc, whereas GetNumberPartitions returns this number. This method is used when Seldon calls the following Slepc function : EPSCISSSetSizes.

Example :

  SparseEigenProblem<T, MatrixK, MatrixM> var_eig;

  var_eig.SetStoppingCriterion(1e-14);
  var_eig.SetNbAskedEigenvalues(nb_eigenval);
  var_eig.SetComputationalMode(var_eig.REGULAR_MODE);

  var_eig.InitMatrix(K, M);

  SlepcParam& param = var_eig.GetSlepcParameters();
  param.SetNumberPartitions(8);
  
  // Large eigenvalues
  var_eig.SetTypeSpectrum(var_eig.LARGE_EIGENVALUES, 0, var_eig.SORTED_MODULUS);
  
  GetEigenvaluesEigenvectors(var_eig, lambda, lambda_imag, eigen_vec, TypeEigenvalueSolver::SLEPC);
  DISP(lambda); DISP(lambda_imag);

Location :

GetThresholdRank, SetThresholdRank

Syntax :

  void SetThresholdRank(double d);
  double GetThresholdRank();

SetThresholdRank sets the rank threshold to use in the solver CISS of Slepc, whereas GetThresholdRank returns this number. This method is used when Seldon calls the following Slepc function : EPSCISSSetThreshold.

Example :

  SparseEigenProblem<T, MatrixK, MatrixM> var_eig;

  var_eig.SetStoppingCriterion(1e-14);
  var_eig.SetNbAskedEigenvalues(nb_eigenval);
  var_eig.SetComputationalMode(var_eig.REGULAR_MODE);

  var_eig.InitMatrix(K, M);

  SlepcParam& param = var_eig.GetSlepcParameters();
  param.SetThresholdRank(1e-8);
  
  // Large eigenvalues
  var_eig.SetTypeSpectrum(var_eig.LARGE_EIGENVALUES, 0, var_eig.SORTED_MODULUS);
  
  GetEigenvaluesEigenvectors(var_eig, lambda, lambda_imag, eigen_vec, TypeEigenvalueSolver::SLEPC);
  DISP(lambda); DISP(lambda_imag);

Location :

GetThresholdSpurious, SetThresholdSpurious

Syntax :

  void SetThresholdSpurious(double d);
  double GetThresholdSpurious();

SetThresholdSpurious sets the spurious threshold to use in the solver CISS of Slepc, whereas GetThresholdRank returns this number. This method is used when Seldon calls the following Slepc function : EPSCISSSetThreshold.

Example :

  SparseEigenProblem<T, MatrixK, MatrixM> var_eig;

  var_eig.SetStoppingCriterion(1e-14);
  var_eig.SetNbAskedEigenvalues(nb_eigenval);
  var_eig.SetComputationalMode(var_eig.REGULAR_MODE);

  var_eig.InitMatrix(K, M);

  SlepcParam& param = var_eig.GetSlepcParameters();
  param.SetThresholdSpurious(1e-8);
  
  // Large eigenvalues
  var_eig.SetTypeSpectrum(var_eig.LARGE_EIGENVALUES, 0, var_eig.SORTED_MODULUS);
  
  GetEigenvaluesEigenvectors(var_eig, lambda, lambda_imag, eigen_vec, TypeEigenvalueSolver::SLEPC);
  DISP(lambda); DISP(lambda_imag);

Location :

GetBorthogonalization, SetBorthogonalization

Syntax :

  void SetBorthogonalization(int b);
  int GetBorthogonalization();

SetBorthogonalization selects the orthogonalization to use in the solver GD/JD of Slepc, whereas GetBorthogonalization returns this orthogonalization. This method is used when Seldon calls the following Slepc functions : EPSGDSetBOrth, EPSJDSetBOrth.

Example :

  SparseEigenProblem<T, MatrixK, MatrixM> var_eig;

  var_eig.SetStoppingCriterion(1e-14);
  var_eig.SetNbAskedEigenvalues(nb_eigenval);
  var_eig.SetComputationalMode(var_eig.REGULAR_MODE);

  var_eig.InitMatrix(K, M);

  SlepcParam& param = var_eig.GetSlepcParameters();
  param.SetBorthogonalization(true);
  
  // Large eigenvalues
  var_eig.SetTypeSpectrum(var_eig.LARGE_EIGENVALUES, 0, var_eig.SORTED_MODULUS);
  
  GetEigenvaluesEigenvectors(var_eig, lambda, lambda_imag, eigen_vec, TypeEigenvalueSolver::SLEPC);
  DISP(lambda); DISP(lambda_imag);

Location :

GetDoubleExpansion, SetDoubleExpansion

Syntax :

  void SetDoubleExpansion(int d);
  int GetDoubleExpansion();

SetDoubleExpansion selects the expansion to use in the solver GD/JD of Slepc, whereas GetDoubleExpansion returns this expansion. This method is used when Seldon calls the following Slepc functions : EPSGDSetDoubleExpansion, EPSJDSetDoubleExpansion.

Example :

  SparseEigenProblem<T, MatrixK, MatrixM> var_eig;

  var_eig.SetStoppingCriterion(1e-14);
  var_eig.SetNbAskedEigenvalues(nb_eigenval);
  var_eig.SetComputationalMode(var_eig.REGULAR_MODE);

  var_eig.InitMatrix(K, M);

  SlepcParam& param = var_eig.GetSlepcParameters();
  param.SetDoubleExpansion(true);
  
  // Large eigenvalues
  var_eig.SetTypeSpectrum(var_eig.LARGE_EIGENVALUES, 0, var_eig.SORTED_MODULUS);
  
  GetEigenvaluesEigenvectors(var_eig, lambda, lambda_imag, eigen_vec, TypeEigenvalueSolver::SLEPC);
  DISP(lambda); DISP(lambda_imag);

Location :

GetInitialSize, SetInitialSize

Syntax :

  void SetInitialSize(int d);
  int GetInitialSize();

SetInitialSize sets the initial size to use in the solver GD/JD of Slepc, whereas GetInitialSize returns this number. This method is used when Seldon calls the following Slepc functions : EPSGDSetInitialSize, EPSJDSetInitialSize.

Example :

  SparseEigenProblem<T, MatrixK, MatrixM> var_eig;

  var_eig.SetStoppingCriterion(1e-14);
  var_eig.SetNbAskedEigenvalues(nb_eigenval);
  var_eig.SetComputationalMode(var_eig.REGULAR_MODE);

  var_eig.InitMatrix(K, M);

  SlepcParam& param = var_eig.GetSlepcParameters();
  param.SetInitialSize(22);
  
  // Large eigenvalues
  var_eig.SetTypeSpectrum(var_eig.LARGE_EIGENVALUES, 0, var_eig.SORTED_MODULUS);
  
  GetEigenvaluesEigenvectors(var_eig, lambda, lambda_imag, eigen_vec, TypeEigenvalueSolver::SLEPC);
  DISP(lambda); DISP(lambda_imag);

Location :

GetKrylovRestart, SetKrylovRestart

Syntax :

  void SetKrylovRestart(int d);
  int GetKrylovRestart();

SetKrylovRestart activates the search with Krylov basis in the solver GD/JD of Slepc, whereas GetKrylovRestart returns this activation. This method is used when Seldon calls the following Slepc functions : EPSGDSetKrylovStart, EPSGDSetKrylovStart.

Example :

  SparseEigenProblem<T, MatrixK, MatrixM> var_eig;

  var_eig.SetStoppingCriterion(1e-14);
  var_eig.SetNbAskedEigenvalues(nb_eigenval);
  var_eig.SetComputationalMode(var_eig.REGULAR_MODE);

  var_eig.InitMatrix(K, M);

  SlepcParam& param = var_eig.GetSlepcParameters();
  param.SetKrylovRestart(true);
  
  // Large eigenvalues
  var_eig.SetTypeSpectrum(var_eig.LARGE_EIGENVALUES, 0, var_eig.SORTED_MODULUS);
  
  GetEigenvaluesEigenvectors(var_eig, lambda, lambda_imag, eigen_vec, TypeEigenvalueSolver::SLEPC);
  DISP(lambda); DISP(lambda_imag);

Location :

GetRestartNumber, SetRestartNumber

Syntax :

  void SetRestartNumber(int d);
  int GetRestartNumber();

SetRestartNumber sets the restart number in the solver GD/JD of Slepc, whereas GetRestartNumber returns this number. This method is used when Seldon calls the following Slepc functions : EPSGDSetRestart, EPSJDSetRestart.

Example :

  SparseEigenProblem<T, MatrixK, MatrixM> var_eig;

  var_eig.SetStoppingCriterion(1e-14);
  var_eig.SetNbAskedEigenvalues(nb_eigenval);
  var_eig.SetComputationalMode(var_eig.REGULAR_MODE);

  var_eig.InitMatrix(K, M);

  SlepcParam& param = var_eig.GetSlepcParameters();
  param.SetRestartNumber(20);
  
  // Large eigenvalues
  var_eig.SetTypeSpectrum(var_eig.LARGE_EIGENVALUES, 0, var_eig.SORTED_MODULUS);
  
  GetEigenvaluesEigenvectors(var_eig, lambda, lambda_imag, eigen_vec, TypeEigenvalueSolver::SLEPC);
  DISP(lambda); DISP(lambda_imag);

Location :

GetRestartNumberAdd, SetRestartNumberAdd

Syntax :

  void SetRestartNumberAdd(int d);
  int GetRestartNumberAdd();

SetRestartNumberAdd sets the number of vectors saved from previous iteration in the solver GD/JD of Slepc, whereas GetRestartNumberAdd returns this number. This method is used when Seldon calls the following Slepc functions : EPSGDSetRestart, EPSJDSetRestart.

Example :

  SparseEigenProblem<T, MatrixK, MatrixM> var_eig;

  var_eig.SetStoppingCriterion(1e-14);
  var_eig.SetNbAskedEigenvalues(nb_eigenval);
  var_eig.SetComputationalMode(var_eig.REGULAR_MODE);

  var_eig.InitMatrix(K, M);

  SlepcParam& param = var_eig.GetSlepcParameters();
  param.SetRestartNumberAdd(3);
  
  // Large eigenvalues
  var_eig.SetTypeSpectrum(var_eig.LARGE_EIGENVALUES, 0, var_eig.SORTED_MODULUS);
  
  GetEigenvaluesEigenvectors(var_eig, lambda, lambda_imag, eigen_vec, TypeEigenvalueSolver::SLEPC);
  DISP(lambda); DISP(lambda_imag);

Location :

GetNumberConvergedVectors, SetNumberConvergedVectors, GetNumberConvergedVectorsProjected, SetNumberConvergedVectorsProjected

Syntax :

  void SetNumberConvergedVectors(int d);
  int GetNumberConvergedVectors();
  void SetNumberConvergedVectorsProjected(int d);
  int GetNumberConvergedVectorsProjected();

SetNumberConvergedVectors sets the number of converged vectors in the projector for the solver GD/JD of Slepc, whereas GetNumberConvergedVectors returns this number. SetNumberConvergedVectorsProjected sets the number of converged vectors in the projected problem, whereas GetNumberConvergedVectorsProjected returns this number.

Example :

  SparseEigenProblem<T, MatrixK, MatrixM> var_eig;

  var_eig.SetStoppingCriterion(1e-14);
  var_eig.SetNbAskedEigenvalues(nb_eigenval);
  var_eig.SetComputationalMode(var_eig.REGULAR_MODE);

  var_eig.InitMatrix(K, M);

  SlepcParam& param = var_eig.GetSlepcParameters();
  param.SetNumberConvergedVectors(10);
  param.SetNumberConvergedVectorsProjected(12);
  
  // Large eigenvalues
  var_eig.SetTypeSpectrum(var_eig.LARGE_EIGENVALUES, 0, var_eig.SORTED_MODULUS);
  
  GetEigenvaluesEigenvectors(var_eig, lambda, lambda_imag, eigen_vec, TypeEigenvalueSolver::SLEPC);
  DISP(lambda); DISP(lambda_imag);

Location :

UseNonLockingVariant, SetNonLockingVariant

Syntax :

  void SetNonLockingVariant(bool d);
  bool UseNonLockingVariant();

SetNonLockingVariant enables the non-locking variant of Slepc solvers, whereas UseNonLockingVariant returns true if the non-locking variant is used. This method is used when Seldon calls the following Slepc functions : EPSKrylovSchurSetLocking, EPSLOBPCGSetLocking.

Example :

  SparseEigenProblem<T, MatrixK, MatrixM> var_eig;

  var_eig.SetStoppingCriterion(1e-14);
  var_eig.SetNbAskedEigenvalues(nb_eigenval);
  var_eig.SetComputationalMode(var_eig.REGULAR_MODE);

  var_eig.InitMatrix(K, M);

  SlepcParam& param = var_eig.GetSlepcParameters();
  param.SetNonLockingVariant(true);
  
  // Large eigenvalues
  var_eig.SetTypeSpectrum(var_eig.LARGE_EIGENVALUES, 0, var_eig.SORTED_MODULUS);
  
  GetEigenvaluesEigenvectors(var_eig, lambda, lambda_imag, eigen_vec, TypeEigenvalueSolver::SLEPC);
  DISP(lambda); DISP(lambda_imag);

Location :

GetRestartRatio, SetRestartRatio

Syntax :

  void SetRestartRatio(double d);
  double GetRestartRatio();

SetRestartRatio sets the restart ratio of Slepc solvers, whereas GetRestartRatio returns this ratio. This method is used when Seldon calls the following Slepc functions : EPSKrylovSchurSetRestart, EPSLOBPCGSetRestart.

Example :

  SparseEigenProblem<T, MatrixK, MatrixM> var_eig;

  var_eig.SetStoppingCriterion(1e-14);
  var_eig.SetNbAskedEigenvalues(nb_eigenval);
  var_eig.SetComputationalMode(var_eig.REGULAR_MODE);

  var_eig.InitMatrix(K, M);

  SlepcParam& param = var_eig.GetSlepcParameters();
  param.SetRestartRatio(0.4);
  
  // Large eigenvalues
  var_eig.SetTypeSpectrum(var_eig.LARGE_EIGENVALUES, 0, var_eig.SORTED_MODULUS);
  
  GetEigenvaluesEigenvectors(var_eig, lambda, lambda_imag, eigen_vec, TypeEigenvalueSolver::SLEPC);
  DISP(lambda); DISP(lambda_imag);

Location :

GetMethod, SetMethod

Syntax :

  void SetMethod(string s);
  string GetMethod();

SetMethod selects the method to use in solver PRIMME of Slepc, whereas GetMethod returns this method. You can choose between DYNAMIC, DEFAULT_MIN_TIME, DEFAULT_MIN_MATVECS, ARNOLDI, GD, GD_PLUSK, GD_OLSEN_PLUSK, JD_OLSEN_PLUSK, RQI, JDQR, JDQMR, JDQMR_ETOL, SUBSPACE_ITERATION, LOBPCG_ORTHOBASIS and LOBPCG_ORTHOBASIS. Those methods are used when Seldon calls the following Slepc function : EPSPRIMMESetMethod.

Example :

  SparseEigenProblem<T, MatrixK, MatrixM> var_eig;

  var_eig.SetStoppingCriterion(1e-14);
  var_eig.SetNbAskedEigenvalues(nb_eigenval);
  var_eig.SetComputationalMode(var_eig.REGULAR_MODE);

  var_eig.InitMatrix(K, M);

  SlepcParam& param = var_eig.GetSlepcParameters();
  param.SetMethod("DYNAMIC");
  
  // Large eigenvalues
  var_eig.SetTypeSpectrum(var_eig.LARGE_EIGENVALUES, 0, var_eig.SORTED_MODULUS);
  
  GetEigenvaluesEigenvectors(var_eig, lambda, lambda_imag, eigen_vec, TypeEigenvalueSolver::SLEPC);
  DISP(lambda); DISP(lambda_imag);

Location :

GetShiftType, SetShiftType

Syntax :

  void SetShiftType(int s);
  int GetShiftType();

SetShiftType selects the type of shift strategy to use in solver POWER of Slepc, whereas GetShiftType returns this type. You can choose between SHIFT_CONSTANT, SHIFT_RAYLEIGH and SHIFT_WILKINSON. Those methods are used when Seldon calls the following Slepc function : EPSPowerSetShiftType.

Example :

  SparseEigenProblem<T, MatrixK, MatrixM> var_eig;

  var_eig.SetStoppingCriterion(1e-14);
  var_eig.SetNbAskedEigenvalues(nb_eigenval);
  var_eig.SetComputationalMode(var_eig.REGULAR_MODE);

  var_eig.InitMatrix(K, M);

  SlepcParam& param = var_eig.GetSlepcParameters();
  param.SetShiftType(param.SHIFT_RAYLEIGH);
  
  // Large eigenvalues
  var_eig.SetTypeSpectrum(var_eig.LARGE_EIGENVALUES, 0, var_eig.SORTED_MODULUS);
  
  GetEigenvaluesEigenvectors(var_eig, lambda, lambda_imag, eigen_vec, TypeEigenvalueSolver::SLEPC);
  DISP(lambda); DISP(lambda_imag);

Location :

SetUserComparisonClass

Syntax :

  void SetUserComparisonClass(EigenvalueComparisonClass<T>* fct);

The user can specify its own ordering of eigenvalues by providing a class derived from EigenValueComparisonClass. This new ordering will be used only with SLEPc eigensolver which implements this functionality.

Example :

// you define a new ordering of eigenvalues
// double for real eigenproblems and complex<double> for complex eigenproblems
class MyOrdering : public EigenvalueComparisonClass<double>
{

// double has to be replaced by complex<double> for complex eigenproblems
// in that case Li and Li2 are not relevant
// result is 0 is two eigenvalues are "equal", 1 if z < z2 and -1 if z > z2
virtual int CompareEigenvalue(const double& Lr, const double& Li,
                              const double& Lr2, const double& Li2)
{
  // forming z and z2 from real and imaginary parts
  complex<double> z(Lr, Li), z2(Lr2, Li2);

  // we sort by ascending |1 + z + z^2|
  double zmod = abs(1 + z + z*z);
  double zmod2 = abs(1 + z2 + z2*z2);
  if (abs(zmod - zmod2) <= 1e-12)
    return 0;
  if (z > zmod2)
    return -1;

  return 1;
}

};

int main()
{
  SparseEigenProblem<double, Matrix<double, General, ArrayRowSparse> > var_eig;

  var_eig.SetStoppingCriterion(1e-14);
  var_eig.SetNbAskedEigenvalues(nb_eigenval);
  var_eig.SetComputationalMode(var_eig.REGULAR_MODE);

  var_eig.InitMatrix(K, M);

  MyOrdering fct;
  var_eig.SetUserComparisonClass(&fct);
  
  // we seek large eigenvalues with user-defined sorting function
  var_eig.SetTypeSpectrum(var_eig.LARGE_EIGENVALUES, 0, var_eig.SORTED_USER);

  // calling Slepc
  GetEigenvaluesEigenvectors(var_eig, lambda, lambda_imag, eigen_vec, TypeEigenvalueSolver::SLEPC);

Location :

GetComparisonEigenvalueSlepc

Syntax :

  void GetComparisonEigenvalueSlepc(T Lr, T Li, T Lr2, T Li2, int* res, void* ctx);

This method is used in the interface of Slepc in order to call the user-defined comparison of two eigenvalues (see method SetUserComparisonClass) for an example. This method should not be called solely.

Location :

FillComplexEigenvectors

Syntax :

  void FillComplexEigenvectors(int m, complex<T> Emid, double epsilon,
                               const Vector<complex<T> >& lambda_cplx,
                               const Matrix<complex<T>, General, ColMajor>& eigen_vec_cplx,
                               Vector<T>& Lr, Vector<T>& Li, Matrix<T, General, ColMajor>& eigen_vec);

This method is used to convert complex eigenvalues/eigenvectors to Lapack form (where a column is formed of real part and the other one of imaginary part). If T is complex, this method just copies lambda_cplx to Lr, and eigen_vec_cplx to eigen_vec. If T is real, it fills Lr and Li, and eigen_vec by using Lapack form since eigenvalues are complex conjugate. This function should not be used in a regular use.

Location :

DistributeEigenvectors

Syntax :

  void DistributeEigenvectors(Matrix<T, General, ColMajor>& eigen_vec);

This method is used to distribute eigenvectors between processors (especially for rows that are shared between processors) for a parallel execution. This function should not be used in a regular use.

Location :

to_complex_eigen

Syntax :

  void to_complex_eigen(const T1& x, T2 y);

This function copies the scalar x to the scalar y. If x is complex and y real, it takes the real part.

Example :

double x = 3; complex<double> y;
to_complex_eigen(x, y); // y should be equal to 3

y = complex<double>(2, -1);
to_complex_eigen(y, x); // x should be equal to 2

Location :

ApplyScalingEigenvec

Syntax :

  void ApplyScalingEigenvec(EigenProblem_Base& var, Vector& lambda, Vector& lambda_imag,
                            Matrix& eigen_vec, T shift_r, T shift_i);

This function modifies the eigenvectors in the two following cases.

mass matrix M is diagonal : the eigenvector has be multiplied by M^-1/2
because eigenvalues/eigenvectors have been computed for the matrix M^1/2 K M ^1/2
mass matrix M has been factorized (Cholesky) : the eigenvector has be multiplied by L^-T
because eigenvalues/eigenvectors have been computed for the matrix L K L ^T where L is the Cholesky factor

The eigenvalues are also modified to be equal to σ + 1/μ (in the case of INVERT_MODE) where σ stands for the shift (shift_r, shift_i). Eigenvectors are also distributed among processors. The user should not call this function, this function is already automatically called by functions like GetEigenvaluesEigenvectors, FindEigenvaluesArpack or FindEigenvaluesSlepc.

Location :

SortEigenvalues

Syntax :

  void SortEigenvalues(Vector& lambda, Vector& lambda_imag, Matrix& eigenvec_old,
                       Matrix& eigenvec_new, int type_spectrum, int type_sort,
                       T shift_r, T shift_i);

This function sorts eigenvalues by ascending order. The associated eigenvectors are stored in the matrix eigenvec_new. This function is actually already called by function GetEigenvaluesEigenvectors.

Location :

EnableEstimateNumberEigenval, EstimateNumberEigenval

Syntax :

  void EnableEstimateNumberEigenval();
  void EnableEstimateNumberEigenval(false);
  bool EstimateNumberEigenval();

The method EnableEstimateNumberEigenval enables the estimate of the number of eigenvalues located in the selected region. If the optional argument is false, this estimate is not performed. If the estimate is enabled, the call to FindEigenvaluesFeast will display on the screen the estimated number of eigenvalues located in the selected region of the spectrum. Finally, the method EstimateNumberEigenval returns true is an estimate is asked, false otherwise.

Example :

// filling matrices
Matrix<complex<double>, General, ArrayRowSparse> M, K;

// declaring the eigenproblem
SparseEigenProblem<complex<double>, Matrix<complex<double>, General, ArrayRowSparse> > var_eig;

var_eig.SetStoppingCriterion(1e-12);
var_eig.SetNbAskedEigenvalues(20); // here you can put a moderate number in order to estimate the actual number of eigenvalues to search
var_eig.InitMatrix(K, M); 
var_eig.SetPrintLevel(1);

// you define the selected region (here a disc) where you want to know how many eigenvalues lie in. 
complex<double> center(0.1, 2.6);
FeastParam& param = var_eig.GetFeastParameters();
param.SetCircleSpectrum(center, 0.2);

// you require the estimate
param.EnableEstimateNumberEigenval();

Vector<complex<double> > lambda, lambda_imag;
Matrix<complex<double>, General, ColMajor> eigen_vec;

// we call Feast to display the estimate
FindEigenvaluesFeast(var_eig, lambda, lambda_imag, eigen_vec);

Location :

GetNumOfQuadraturePoints, SetNumOfQuadraturePoints

Syntax :

  int GetNumOfQuadraturePoints();
  void SetNumOfQuadraturePoints(int);

The method GetNumOfQuadraturePoints returns the number of quadrature points to use when calling Feast. The method SetNumOfQuadraturePoints sets the number of quadrature points to use. By default (if this method is not called), eight points will be used for an hermitian eigenproblem and 16 points for a general eigenproblem.

Example :

// filling matrices
Matrix<complex<double>, General, ArrayRowSparse> M, K;

// declaring the eigenproblem
SparseEigenProblem<complex<double>, Matrix<complex<double>, General, ArrayRowSparse> > var_eig;

var_eig.SetStoppingCriterion(1e-12);
var_eig.SetNbAskedEigenvalues(20);
var_eig.InitMatrix(K, M); 
var_eig.SetPrintLevel(1);

complex<double> center(0.1, 2.6);
FeastParam& param = var_eig.GetFeastParameters();
param.SetCircleSpectrum(center, 0.2);

// if you want to use a different number of quadrature points
param.SetNumOfQuadraturePoints(12);

Vector<complex<double> > lambda, lambda_imag;
Matrix<complex<double>, General, ColMajor> eigen_vec;

// we call Feast to compute the eigenvalues
FindEigenvaluesFeast(var_eig, lambda, lambda_imag, eigen_vec);

Location :

GetTypeIntegration

Syntax :

  int GetTypeIntegration();
  void SetTypeIntegration(int);

The method GetTypeIntegration returns the integration method to use in Feast (0, 1 or 2). The method SetTypeIntegration sets the integration method to use. By default (if this method is not called), the default method will be selected by Feast.

Example :

// filling matrices
Matrix<double, Symmetric, ArrayRowSymSparse> M, K;

// declaring the eigenproblem
SparseEigenProblem<double, Matrix<double, Symmetric, ArrayRowSymSparse> > var_eig;

var_eig.SetStoppingCriterion(1e-12);
var_eig.SetNbAskedEigenvalues(20);
var_eig.InitMatrix(K, M); 
var_eig.SetPrintLevel(1);

complex<double> center(0.1, 2.6);
FeastParam& param = var_eig.GetFeastParameters();
param.SetCircleSpectrum(center, 0.2);

// if you want to use a different number of quadrature points
param.SetNumOfQuadraturePoints(12);
// and a different integration method
// 0 : Gauss, 1 : Trapezoidal, 2 : Zolotarev
param.SetTypeIntegration(2);

Vector<complex<double> > lambda, lambda_imag;
Matrix<complex<double>, General, ColMajor> eigen_vec;

// we call Feast to compute the eigenvalues
FindEigenvaluesFeast(var_eig, lambda, lambda_imag, eigen_vec);

Location :

SetCircleSpectrum, GetCircleCenterSpectrum, GetCircleRadiusSpectrum

Syntax :

  void SetCircleSpectrum(complex<double> z, double r);
  complex<double> GetCircleCenterSpectrum();
  double GetCircleRadiusSpectrum();

The method SetCircleSpectrum sets a circle where eigenvalues have to be searched by Feast. It is relevant only for non-hermitian eigenproblems. The methods GetCircleCenterSpectrum and GetCircleRadiusSpectrum return the center of the circle and the radius, respectively.

Example :

// filling matrices
Matrix<complex<double>, General, ArrayRowSparse> M, K;

// declaring the eigenproblem
SparseEigenProblem<complex<double>, Matrix<complex<double>, General, ArrayRowSparse> > var_eig;

var_eig.SetStoppingCriterion(1e-12);
var_eig.SetNbAskedEigenvalues(20);
var_eig.InitMatrix(K, M); 
var_eig.SetPrintLevel(1);

// here we want to find eigenvalues located inside a disc
// of center 0.1 + 2.6i and of radius 0.2
complex<double> center(0.1, 2.6);
FeastParam& param = var_eig.GetFeastParameters();
param.SetCircleSpectrum(center, 0.2);

Vector<complex<double> > lambda, lambda_imag;
Matrix<complex<double>, General, ColMajor> eigen_vec;

// we call Feast to compute the eigenvalues
FindEigenvaluesFeast(var_eig, lambda, lambda_imag, eigen_vec);

Location :

SetEllipseSpectrum, GetRatioEllipseSpectrum, GetAngleEllipseSpectrum

Syntax :

  void SetEllipseSpectrum(complex<double> z, double r, double ratio, double angle);
  double GetRatioEllipseSpectrum();
  double GetAngleEllipseSpectrum();

The method SetEllipseSpectrum sets an ellipse where eigenvalues have to be searched by Feast. It is relevant only for non-hermitian eigenproblems. The methods GetCircleCenterSpectrum, GetCircleRadiusSpectrum, GetRatioEllipseSpectrum and GetAngleEllipseSpectrum return the center, horizontal radius, ratio and angle of the ellipse respectively. The ratio of the ellipse is the ratio (in percents) of the vertical axis over the horizontal axis. A ratio equal to 100 corresponds to a circle. The angle of the ellipse is an angle in degrees between -180 and 180.

Example :

// filling matrices
Matrix<complex<double>, General, ArrayRowSparse> M, K;

// declaring the eigenproblem
SparseEigenProblem<complex<double>, Matrix<complex<double>, General, ArrayRowSparse> > var_eig;

var_eig.SetStoppingCriterion(1e-12);
var_eig.SetNbAskedEigenvalues(20);
var_eig.InitMatrix(K, M); 
var_eig.SetPrintLevel(1);

// here we want to find eigenvalues located inside an ellipse
// of center 0.1 + 2.6i and of horizontal radius of 0.2, vertical radius of 0.1, and
// oriented with an angle of 45 degrees
complex<double> center(0.1, 2.6);
FeastParam& param = var_eig.GetFeastParameters();
param.SetEllipseSpectrum(center, 0.2, 50.0, 45.0);

Vector<complex<double> > lambda, lambda_imag;
Matrix<complex<double>, General, ColMajor> eigen_vec;

// we call Feast to compute the eigenvalues
FindEigenvaluesFeast(var_eig, lambda, lambda_imag, eigen_vec);

Location :

FindEigenvaluesSlepc

Syntax :

  void FindEigenvaluesSlepc(EigenProblem_Base& var, Vector& lambda, Vector& lambda_imag, Matrix& eigen_vec);
  void FindEigenvaluesSlepc(PolynomialEigenProblem_Base& var, Vector& lambda, Vector& lambda_imag, Matrix& eigen_vec);
  void FindEigenvaluesSlepc(NonLinearEigenProblem_Base& var, Vector& lambda, Vector& lambda_imag, Matrix& eigen_vec);

These functions call Slepc in order to compute the eigenvalues (lambda_imag is the imaginary part for real eigenproblems) and eigenvectors. The eigenproblem can be linear (derived class of EigenProblem_Base), polynomial (derived class of PolynomialEigenProblem_Base) or non-linear (derived class of NonLinearEigenProblem_Base). The eigenvalues are placed in lambda and lambda_imag, the eigenvectors in eigen_vec.

Example :

// filling matrices
Matrix<complex<double>, General, ArrayRowSparse> M, K;

// declaring the eigenproblem
// here SparseEigenProblem derives from EigenProblem_Base
SparseEigenProblem<complex<double>, Matrix<complex<double>, General, ArrayRowSparse> > var_eig;

var_eig.SetStoppingCriterion(1e-12);
var_eig.SetNbAskedEigenvalues(20);
var_eig.InitMatrix(K, M); 
var_eig.SetPrintLevel(1);

complex<double> center(0.1, 2.6);
var_eig.SetTypeSpectrum(var_eig.CENTERED_EIGENVALUES, center);
var_eig.SetComputationalMode(var_eig.INVERT_MODE);

Vector<complex<double> > lambda, lambda_imag;
Matrix<complex<double>, General, ColMajor> eigen_vec;

// we call Slepc to find eigenvalues
FindEigenvaluesSlepc(var_eig, lambda, lambda_imag, eigen_vec);

Location :

Slepc.cxx

FindEigenvaluesArpack

Syntax :

  void FindEigenvaluesArpack(EigenProblem_Base& var, Vector& lambda, Vector& lambda_imag, Matrix& eigen_vec);

This function calls Arpack in order to compute the eigenvalues (lambda_imag is the imaginary part for real eigenproblems) and eigenvectors. The eigenproblem is linear, i.e. a derived class of EigenProblem_Base. The eigenvalues are placed in lambda and lambda_imag, the eigenvectors in eigen_vec.

Example :

// filling matrices
Matrix<complex<double>, General, ArrayRowSparse> M, K;

// declaring the eigenproblem
// here SparseEigenProblem derives from EigenProblem_Base
SparseEigenProblem<complex<double>, Matrix<complex<double>, General, ArrayRowSparse> > var_eig;

var_eig.SetStoppingCriterion(1e-12);
var_eig.SetNbAskedEigenvalues(20);
var_eig.InitMatrix(K, M); 
var_eig.SetPrintLevel(1);

complex<double> center(0.1, 2.6);
var_eig.SetTypeSpectrum(var_eig.CENTERED_EIGENVALUES, center);
var_eig.SetComputationalMode(var_eig.INVERT_MODE);

Vector<complex<double> > lambda, lambda_imag;
Matrix<complex<double>, General, ColMajor> eigen_vec;

// we call Arpack to find eigenvalues
FindEigenvaluesArpack(var_eig, lambda, lambda_imag, eigen_vec);

Location :

Arpack.cxx

FindEigenvaluesFeast

Syntax :

  void FindEigenvaluesArpack(EigenProblem_Base& var, Vector& lambda, Vector& lambda_imag, Matrix& eigen_vec);
  void FindEigenvaluesArpack(PolynomialEigenProblem_Base& var, Vector& lambda, Vector& lambda_imag, Matrix& eigen_vec);

This function calls Arpack in order to compute the eigenvalues (lambda_imag is the imaginary part for real eigenproblems) and eigenvectors. The eigenproblem can be linear (derived class of EigenProblem_Base) or polynomial (derived class of PolynomialEigenProblem_Base). The eigenvalues are placed in lambda and lambda_imag, the eigenvectors in eigen_vec.

Example :

// filling matrices
Matrix<complex<double>, General, ArrayRowSparse> M, K;

// declaring the eigenproblem
// here SparseEigenProblem derives from EigenProblem_Base
SparseEigenProblem<complex<double>, Matrix<complex<double>, General, ArrayRowSparse> > var_eig;

var_eig.SetStoppingCriterion(1e-12);
var_eig.SetNbAskedEigenvalues(20);
var_eig.InitMatrix(K, M); 
var_eig.SetPrintLevel(1);

// Feast needs to define precisely the region where eigenvalues are searched
complex<double> center(0.1, 2.6);
var_eig.SetTypeSpectrum(var_eig.CENTERED_EIGENVALUES, center);
FeastParam& param = var_eig.GetFeastParameters();
param.SetCircleSpectrum(center, 0.2);

Vector<complex<double> > lambda, lambda_imag;
Matrix<complex<double>, General, ColMajor> eigen_vec;

// we call Feast to find eigenvalues
FindEigenvaluesFeast(var_eig, lambda, lambda_imag, eigen_vec);

Location :

Feast.cxx

GetEigensolverType, SetEigensolverType

Syntax :

  void SetEigensolverType(int type);
  int GetEigensolverType();

SetEigensolverType sets the type of eigensolver to use in Slepc, whereas GetEigensolverType returns this number. You can choose between TOAR, STOAR, QARNOLDI, LINEAR or JD. This list corresponds to the list provided by Slepc. By default, the selected solver is TOAR.

Example :

  Matrix<double, General, ArrayRowSparse> K, S, M;
  Vector<Matrix<double, General, ArrayRowSparse>* > vecK(2);
  vecK(0) = &K; vecK(1) = &S;
  
  PolynomialSparseEigenProblem<double, Matrix<double, General, ArrayRowSparse> > var_eig;

  var_eig.SetStoppingCriterion(1e-14);
  var_eig.SetNbAskedEigenvalues(10);
  var_eig.InitMatrix(vecK, M);

  // you can select an alternative solver
  SlepcParamPep& param = var_eig.GetSlepcParameters();
  param.SetEigensolverType(param.QARNOLDI);

  // clustered eigenvalues
  var_eig.SetTypeSpectrum(var_eig.CENTERED_EIGENVALUES, -6.0);
  var_eig.SetSpectralTransformation(true);
    
  GetEigenvaluesEigenvectors(var_eig, lambda, lambda_imag, eigen_vec, TypeEigenvalueSolver::SLEPC);
  DISP(lambda); DISP(lambda_imag);

Location :

UseSpectralTransformation, SetSpectralTransformation

Syntax :

  void SetSpectralTransformation(bool t = true);
  bool UseSpectralTransformation();

SetSpectralTransformation enables the use of a spectral transformation, i.e. we search θ such that
λ = σ + 1/θ
where σ is the shift. With this spectral transformation, a linear system needs to be factorized and not only the mass matrix. However, this spectral transformation is strongly recommended if you want to compute eigenvalues close to a given shift, since it will require a small number of iterations. UseSpectralTransformation returns true if a spectral transformation has to be used.

Example :

  // polynomial eigenproblem (K + lambda S + lambda^2 M) x = 0
  Matrix<double, General, ArrayRowSparse> K, S, M;
  Vector<Matrix<double, General, ArrayRowSparse>* > vecK(2);
  vecK(0) = &K; vecK(1) = &S;
  
  PolynomialSparseEigenProblem<double, Matrix<double, General, ArrayRowSparse> > var_eig;

  var_eig.SetStoppingCriterion(1e-14);
  var_eig.SetNbAskedEigenvalues(10);
  var_eig.InitMatrix(vecK, M);

  // you can select an alternative solver
  SlepcParamPep& param = var_eig.GetSlepcParameters();
  param.SetEigensolverType(param.QARNOLDI);

  // clustered eigenvalues => spectral transformation for a fast computation
  var_eig.SetTypeSpectrum(var_eig.CENTERED_EIGENVALUES, -6.0);
  var_eig.SetSpectralTransformation(true);
    
  GetEigenvaluesEigenvectors(var_eig, lambda, lambda_imag, eigen_vec, TypeEigenvalueSolver::SLEPC);
  DISP(lambda); DISP(lambda_imag);

Location :

GetSlepcParameters

Syntax :

  SlepcParamPep& GetSlepcParameters();

This method returns the Slepc parameters associated with PEP solver.

Example :

  // polynomial eigenproblem (K + lambda S + lambda^2 M) x = 0
  Matrix<complex<double>, General, ArrayRowSparse> K, S, M;
  Vector<Matrix<complex<double>, General, ArrayRowSparse>* > vecK(2);
  vecK(0) = &K; vecK(1) = &S;
  
  PolynomialSparseEigenProblem<complex<double>, Matrix<complex<double>, General, ArrayRowSparse> > var_eig;

  var_eig.SetStoppingCriterion(1e-14);
  var_eig.SetNbAskedEigenvalues(10);
  var_eig.InitMatrix(vecK, M);

  complex<double> center(-6.0, 0.5);  
  var_eig.SetTypeSpectrum(var_eig.CENTERED_EIGENVALUES, center);

  // you can select an alternative solver
  SlepcParamPep& param = var_eig.GetSlepcParameters();
  param.SetEigensolverType(param.QARNOLDI);

  // Slepc is called    
  FindEigenvaluesSlepc(var_eig, lambda, lambda_imag, eigen_vec);

Location :

GetFeastParameters

Syntax :

  FeastParam& GetFeastParameters();

This method returns the Feast parameters associated with the eigenproblem.

Example :

  // polynomial eigenproblem (K + lambda S + lambda^2 M) x = 0
  Matrix<complex<double>, General, ArrayRowSparse> K, S, M;
  Vector<Matrix<complex<double>, General, ArrayRowSparse>* > vecK(2);
  vecK(0) = &K; vecK(1) = &S;
  
  PolynomialSparseEigenProblem<complex<double>, Matrix<complex<double>, General, ArrayRowSparse> > var_eig;

  var_eig.SetStoppingCriterion(1e-14);
  var_eig.SetNbAskedEigenvalues(10);
  var_eig.InitMatrix(vecK, M);

  // With Feast, you need to define the region in the complex plane
  // where eigenvalues are researched
  complex<double> center(-6.0, 0.5);  
  var_eig.SetTypeSpectrum(var_eig.CENTERED_EIGENVALUES, center);
  FeastParam& param = var_eig.GetFeastParameters();
  param.SetCircleSpectrum(center, 0.2);

  // Feast is called    
  FindEigenvaluesFeast(var_eig, lambda, lambda_imag, eigen_vec);

Location :

GetPolynomialDegree

Syntax :

  int GetPolynomialDegree()

This method returns the polynomial degree of the eigenproblem.

Example :

  // polynomial eigenproblem (K + lambda S + lambda^2 M) x = 0
  Matrix<complex<double>, General, ArrayRowSparse> K, S, M;
  Vector<Matrix<complex<double>, General, ArrayRowSparse>* > vecK(2);
  vecK(0) = &K; vecK(1) = &S;
  
  PolynomialSparseEigenProblem<complex<double>, Matrix<complex<double>, General, ArrayRowSparse> > var_eig;

  var_eig.SetStoppingCriterion(1e-14);
  var_eig.SetNbAskedEigenvalues(10);
  var_eig.InitMatrix(vecK, M);

  // degree should be equal to 2 for this case
  cout << "Degree = " << var_eig.GetPolynomialDegree() << endl;

Location :

SetDiagonalMass, DiagonalMass

Syntax :

void SetDiagonalMass(bool diag = true);
bool DiagonalMass();

The method SetDiagonalMass informs that the mass matrix is actually diagonal. This information is useful, since it will induce fast solutions with the mass matrix (completed by the function SolveMass). The method DiagonalMass returns true if the mass matrix is diagonal.

Example :

  // polynomial eigenproblem (K + lambda S + lambda^2 M) x = 0
  Matrix<complex<double>, General, ArrayRowSparse> K, S, M;
  Vector<Matrix<complex<double>, General, ArrayRowSparse>* > vecK(2);
  vecK(0) = &K; vecK(1) = &S;
  
  PolynomialSparseEigenProblem<complex<double>, Matrix<complex<double>, General, ArrayRowSparse> > var_eig;

  var_eig.SetStoppingCriterion(1e-14);
  var_eig.SetNbAskedEigenvalues(10);
  var_eig.InitMatrix(vecK, M);

  // you know that M is diagonal, you provide that information
  var_eig.SetDiagonalMass();
  
  // large eigenvalues
  var_eig.SetTypeSpectrum(var_eig.LARGE_EIGENVALUES, 0.0);
  
  // Slepc is called    
  FindEigenvaluesSlepc(var_eig, lambda, lambda_imag, eigen_vec);

Location :

ComputeOperator, MltOperator

Syntax :

void ComputeOperator(int num, const Vector& coef);
void MltOperator(int num, const SeldonTranspose& trans, const Vector& X, Vector& Y);

The method ComputeOperator is called only if a spectral transformation is enabled (see SetSpectralTransformation). It requires that the operator num is computed as
T_num = coef₀ A₀ + coef₁ A₁ + coef_d-1 A_d-1 + coef_d M
Then a call to MltOperator will compute
Y = T_num X or Y = T^T_num X
If a spectral transformation is not required, MltOperator is supposed to compute
Y = A_num X or Y = A_num^T X
These two virtual methods can be overloaded if you want to define your own class defining the eigenproblem.

Location :

FactorizeMass, SolveMass

Syntax :

void FactorizeMass();
void SolveMass(const SeldonTranspose& trans, const Vector& X, Vector& Y);

The method FactorizeMass is called to factorize the mass matrix M. Then a call to SolveMass will solve the linear system
M Y = X or M^T Y = X
These two virtual methods can be overloaded if you want to define your own class defining the eigenproblem.

Location :

FactorizeOperator, SolveOperator

Syntax :

void FactorizeOperator(const Vector& coef);
void SolveOperator(const SeldonTranspose& trans, const Vector& X, Vector& Y);

The method ComputeOperator is called only if a spectral transformation is enabled (see SetSpectralTransformation). It requires that the following operator
OP = coef₀ A₀ + coef₁ A₁ + coef_d-1 A_d-1 + coef_d M
is factorized. Then a call to SolveOperator will solve
OP Y = X or OP^T Y = X
These two virtual methods can be overloaded if you want to define your own class defining the eigenproblem.

Location :

InitMatrix

Syntax :

for PolynomialDenseEigenProblem
 void InitMatrix(const Vector<Matrix<T, Prop, Storage>* >& Ai);

for PolynomialEigenProblem
void InitMatrix(const Vector<VirtualMatrix<T>* >& Ai, int n = -1);

for PolynomialSparseEigenProblem
void InitMatrix(const Vector<MatStiff*>& Ai, MatMass& K);

This method initializes the matrices A_i. The syntax differ depending on the chosen eigenproblem. For PolynomialEigenProblem, the optional argument is the local number of rows.

Example :

  // polynomial eigenproblem (K + lambda S + lambda^2 M) x = 0
  Matrix<complex<double>, General, ArrayRowSparse> K, S, M;
  Vector<Matrix<complex<double>, General, ArrayRowSparse>* > vecAi(3);
  vecAi(0) = &K; vecAi(1) = &S; vecAi(2) = &M;

  // Here we have stored matrices K, S, M, but you can implement matrix-free structures
  // if you want to use matrix-free eigenproblem
  PolynomialEigenProblem<complex<double> > var_eig;

  var_eig.SetStoppingCriterion(1e-14);
  var_eig.SetNbAskedEigenvalues(10);
  var_eig.InitMatrix(vecAi); // you provide the matrices Ai

  // large eigenvalues are searched
  var_eig.SetTypeSpectrum(var_eig.LARGE_EIGENVALUES, 0.0);
  
  // Slepc is called    
  FindEigenvaluesSlepc(var_eig, lambda, lambda_imag, eigen_vec);

Location :

GetEigensolverType, SetEigensolverType

Syntax :

  void SetEigensolverType(int type);
  int GetEigensolverType();

SetEigensolverType sets the type of eigensolver to use in Slepc, whereas GetEigensolverType returns this number. You can choose between RII, SLP, NARNOLDI, CISS, INTERPOL or NLEIGS. This list corresponds to the list provided by Slepc. By default, the selected solver is RII.

Example :

// non-linear eigenproblem
SplitSparseNonLinearEigenProblem<Petsc_Scalar, Symmetric, ArrayRowSymSparse> var_eig;
var_eig.SetStoppingCriterion(1e-12);
var_eig.SetNbAskedEigenvalues(nb_eigenval);
// matrices A_i and polynomial P_i and Q_i are provided
var_eig.InitMatrix(vec_Ai, coef_Pi, coef_Qi);
      
// eigenvalues close to a shift sigma are searched
var_eig.SetTypeSpectrum(var_eig.CENTERED_EIGENVALUES, complex<double>(1.0, 0.2));

// if you want to use Nleigs solver
SlepcParamNep& param = var_eig.GetSlepcParameters();
param.SetEigensolverType(param.NLEIGS);
// for this solver a region [a, b] x [c, d] is required
param.SetIntervalRegion(0.3, 0.6, -0.3, -1e-4);

// eigenvalues and eigenvectors are computed by calling Slepc
Vector<Petsc_Scalar> lambda, lambda_imag;
Matrix<Petsc_Scalar, General, ColMajor> eigen_vec;
      
FindEigenvaluesSlepc(var_eig, lambda, lambda_imag, eigen_vec);

Location :

GetLrMin, GetLrmax, GetLiMin, GetLiMax

Syntax :

  double GetLrMin();
  double GetLrMax();
  double GetLiMin();
  double GetLiMax();

These methods return the extremities of the region where eigenvalues are searched.

Example :

// non-linear eigenproblem
SplitSparseNonLinearEigenProblem<Petsc_Scalar, Symmetric, ArrayRowSymSparse> var_eig;
var_eig.SetStoppingCriterion(1e-12);
var_eig.SetNbAskedEigenvalues(nb_eigenval);
// matrices A_i and polynomial P_i and Q_i are provided
var_eig.InitMatrix(vec_Ai, coef_Pi, coef_Qi);
      
// eigenvalues close to a shift sigma are searched
var_eig.SetTypeSpectrum(var_eig.CENTERED_EIGENVALUES, complex<double>(1.0, 0.2));

// if you want to use Nleigs solver
SlepcParamNep& param = var_eig.GetSlepcParameters();
param.SetEigensolverType(param.NLEIGS);
// for this solver a region [a, b] x [c, d] is required
param.SetIntervalRegion(0.3, 0.6, -0.3, -1e-4);
// param.GetLrMin() should return 0.3, param.GetLrMax() 0.6, param.GetLiMin() -0.3 and param.GetLiMax() -1e-4

// eigenvalues and eigenvectors are computed by calling Slepc
Vector<Petsc_Scalar> lambda, lambda_imag;
Matrix<Petsc_Scalar, General, ColMajor> eigen_vec;
      
FindEigenvaluesSlepc(var_eig, lambda, lambda_imag, eigen_vec);

Location :

SetIntervalRegion

Syntax :

  void SetIntervalRegion(double a, double b, double c, double d);
  bool InsideRegion(complex<double> z);

The method SetInterval sets the region where eigenvalues are searched. Real part of eigenvalues should lie in interval [a, b]. Imaginary part of eigenvalues should lie in interval [c, d]. The method InsideRegion returns true if the argument z lies inside the region.

Example :

// non-linear eigenproblem
SplitSparseNonLinearEigenProblem<Petsc_Scalar, Symmetric, ArrayRowSymSparse> var_eig;
var_eig.SetStoppingCriterion(1e-12);
var_eig.SetNbAskedEigenvalues(nb_eigenval);
// matrices A_i and polynomial P_i and Q_i are provided
var_eig.InitMatrix(vec_Ai, coef_Pi, coef_Qi);
      
// eigenvalues close to a shift sigma are searched
var_eig.SetTypeSpectrum(var_eig.CENTERED_EIGENVALUES, complex<double>(0.4, -0.1));

// if you want to use Nleigs solver
SlepcParamNep& param = var_eig.GetSlepcParameters();
param.SetEigensolverType(param.NLEIGS);
// for this solver a region [a, b] x [c, d] is required
param.SetIntervalRegion(0.3, 0.6, -0.3, -1e-4);

// should return true
bool res = param.InsideRegion(complex<double>(0.45, -0.1));

// eigenvalues and eigenvectors are computed by calling Slepc
Vector<Petsc_Scalar> lambda, lambda_imag;
Matrix<Petsc_Scalar, General, ColMajor> eigen_vec;
      
FindEigenvaluesSlepc(var_eig, lambda, lambda_imag, eigen_vec);

Location :

EnableCommandLineOptions, UseCommandLineOptions

Syntax :

  void EnableCommandLineOptions(bool yes = true);
  bool UseCommandLineOptions();

The method EnableCommandLineOptions enables the use of options given in the command line. The method UseCommandLineOptions returns true if it is enabled.

Example :

// non-linear eigenproblem
SplitSparseNonLinearEigenProblem<Petsc_Scalar, Symmetric, ArrayRowSymSparse> var_eig;
var_eig.SetStoppingCriterion(1e-12);
var_eig.SetNbAskedEigenvalues(nb_eigenval);
// matrices A_i and polynomial P_i and Q_i are provided
var_eig.InitMatrix(vec_Ai, coef_Pi, coef_Qi);
      
// eigenvalues close to a shift sigma are searched
var_eig.SetTypeSpectrum(var_eig.CENTERED_EIGENVALUES, complex<double>(0.4, -0.1));

// if you want to use Nleigs solver
SlepcParamNep& param = var_eig.GetSlepcParameters();
param.SetEigensolverType(param.NLEIGS);
// for this solver a region [a, b] x [c, d] is required
param.SetIntervalRegion(0.3, 0.6, -0.3, -1e-4);

// you can autorize the user to provide parameters in the command line
param.EnableCommandLineOptions();


// eigenvalues and eigenvectors are computed by calling Slepc
Vector<Petsc_Scalar> lambda, lambda_imag;
Matrix<Petsc_Scalar, General, ColMajor> eigen_vec;
      
FindEigenvaluesSlepc(var_eig, lambda, lambda_imag, eigen_vec);

Location :

SetDefaultPetscSolver, UseDefaultPetscSolver

Syntax :

  void SetDefaultPetscSolver(bool yes = true);
  bool UseDefaultPetscSolver();

The method SetDefaultPetscSolver informs that the linear systems will be solved by the default Petsc solver. The method UseDefaultPetscSolver returns true if the default Petsc solver will be used.

Example :

// non-linear eigenproblem
SplitSparseNonLinearEigenProblem<Petsc_Scalar, Symmetric, ArrayRowSymSparse> var_eig;
var_eig.SetStoppingCriterion(1e-12);
var_eig.SetNbAskedEigenvalues(nb_eigenval);
// matrices A_i and polynomial P_i and Q_i are provided
var_eig.InitMatrix(vec_Ai, coef_Pi, coef_Qi);
      
// eigenvalues close to a shift sigma are searched
var_eig.SetTypeSpectrum(var_eig.CENTERED_EIGENVALUES, complex<double>(0.4, -0.1));

// if you want to use Nleigs solver
SlepcParamNep& param = var_eig.GetSlepcParameters();
param.SetEigensolverType(param.NLEIGS);
// for this solver a region [a, b] x [c, d] is required
param.SetIntervalRegion(0.3, 0.6, -0.3, -1e-4);

// you can use Petsc solver instead of Seldon solver
param.SetDefaultPetscSolver();

// eigenvalues and eigenvectors are computed by calling Slepc
Vector<Petsc_Scalar> lambda, lambda_imag;
Matrix<Petsc_Scalar, General, ColMajor> eigen_vec;
      
FindEigenvaluesSlepc(var_eig, lambda, lambda_imag, eigen_vec);

Location :

SetFullBasis, FullBasis

Syntax :

  void SetFullBasis(bool yes = true);
  bool FullBasis();

The method SetFullBasis informs that a full basis has to be used (it corresponds to NEPNLEIGSSetFullBasis Slepc function). The method FullBasis returns true if the full basis will be used.

Example :

// non-linear eigenproblem
SplitSparseNonLinearEigenProblem<Petsc_Scalar, Symmetric, ArrayRowSymSparse> var_eig;
var_eig.SetStoppingCriterion(1e-12);
var_eig.SetNbAskedEigenvalues(nb_eigenval);
// matrices A_i and polynomial P_i and Q_i are provided
var_eig.InitMatrix(vec_Ai, coef_Pi, coef_Qi);
      
// eigenvalues close to a shift sigma are searched
var_eig.SetTypeSpectrum(var_eig.CENTERED_EIGENVALUES, complex<double>(0.4, -0.1));

// if you want to use Nleigs solver
SlepcParamNep& param = var_eig.GetSlepcParameters();
param.SetEigensolverType(param.NLEIGS);
// for this solver a region [a, b] x [c, d] is required
param.SetIntervalRegion(0.3, 0.6, -0.3, -1e-4);

// you can use a full basis
param.SetFullBasis();

// eigenvalues and eigenvectors are computed by calling Slepc
Vector<Petsc_Scalar> lambda, lambda_imag;
Matrix<Petsc_Scalar, General, ColMajor> eigen_vec;
      
FindEigenvaluesSlepc(var_eig, lambda, lambda_imag, eigen_vec);

Location :

SetLockingVariant, LockingVariant

Syntax :

  void SetLockingVariant(bool yes = true);
  bool LockingVariant();

The method SetLockingVariant informs that a locking variant has to be used (it corresponds to NEPNLEIGSSetLocking Slepc function). The method LockingVariant returns true if the locking variant will be used.

Example :

// non-linear eigenproblem
SplitSparseNonLinearEigenProblem<Petsc_Scalar, Symmetric, ArrayRowSymSparse> var_eig;
var_eig.SetStoppingCriterion(1e-12);
var_eig.SetNbAskedEigenvalues(nb_eigenval);
// matrices A_i and polynomial P_i and Q_i are provided
var_eig.InitMatrix(vec_Ai, coef_Pi, coef_Qi);
      
// eigenvalues close to a shift sigma are searched
var_eig.SetTypeSpectrum(var_eig.CENTERED_EIGENVALUES, complex<double>(0.4, -0.1));

// if you want to use Nleigs solver
SlepcParamNep& param = var_eig.GetSlepcParameters();
param.SetEigensolverType(param.NLEIGS);
// for this solver a region [a, b] x [c, d] is required
param.SetIntervalRegion(0.3, 0.6, -0.3, -1e-4);

// you can use a locking variant
param.SetLockingVariant();

// eigenvalues and eigenvectors are computed by calling Slepc
Vector<Petsc_Scalar> lambda, lambda_imag;
Matrix<Petsc_Scalar, General, ColMajor> eigen_vec;
      
FindEigenvaluesSlepc(var_eig, lambda, lambda_imag, eigen_vec);

Location :

SetInterpolationDegree, GetInterpolationDegree

Syntax :

  void SetInterpolationDegree(int p);
  int GetInterpolationDegree();

The method SetInterpolationDegree provides the polynomial degree to be used. The method GetInterpolationDegree returns the polynomial degree to use.

Example :

// non-linear eigenproblem
SplitSparseNonLinearEigenProblem<Petsc_Scalar, Symmetric, ArrayRowSymSparse> var_eig;
var_eig.SetStoppingCriterion(1e-12);
var_eig.SetNbAskedEigenvalues(nb_eigenval);
// matrices A_i and polynomial P_i and Q_i are provided
var_eig.InitMatrix(vec_Ai, coef_Pi, coef_Qi);
      
// eigenvalues close to a shift sigma are searched
var_eig.SetTypeSpectrum(var_eig.CENTERED_EIGENVALUES, complex<double>(0.4, -0.1));

// if you want to use Nleigs solver
SlepcParamNep& param = var_eig.GetSlepcParameters();
param.SetEigensolverType(param.NLEIGS);
// for this solver a region [a, b] x [c, d] is required
param.SetIntervalRegion(0.3, 0.6, -0.3, -1e-4);

// you can enforce a polynomial degree
param.SetInterpolationDegree(6);

// eigenvalues and eigenvectors are computed by calling Slepc
Vector<Petsc_Scalar> lambda, lambda_imag;
Matrix<Petsc_Scalar, General, ColMajor> eigen_vec;
      
FindEigenvaluesSlepc(var_eig, lambda, lambda_imag, eigen_vec);

Location :

SetInterpolationTolerance, GetInterpolationTolerance

Syntax :

  void SetInterpolationTolerance(double eps);
  double GetInterpolationTolerance();

The method SetInterpolationTolerance provides the threshold to be used to compute the interpolation. The method GetInterpolationTolerance returns this threshold.

Example :

// non-linear eigenproblem
SplitSparseNonLinearEigenProblem<Petsc_Scalar, Symmetric, ArrayRowSymSparse> var_eig;
var_eig.SetStoppingCriterion(1e-12);
var_eig.SetNbAskedEigenvalues(nb_eigenval);
// matrices A_i and polynomial P_i and Q_i are provided
var_eig.InitMatrix(vec_Ai, coef_Pi, coef_Qi);
      
// eigenvalues close to a shift sigma are searched
var_eig.SetTypeSpectrum(var_eig.CENTERED_EIGENVALUES, complex<double>(0.4, -0.1));

// if you want to use Nleigs solver
SlepcParamNep& param = var_eig.GetSlepcParameters();
param.SetEigensolverType(param.NLEIGS);
// for this solver a region [a, b] x [c, d] is required
param.SetIntervalRegion(0.3, 0.6, -0.3, -1e-4);

// you can specify a tolerance between interpolation and true function
param.SetInterpolationTolerance(1e-12);

// eigenvalues and eigenvectors are computed by calling Slepc
Vector<Petsc_Scalar> lambda, lambda_imag;
Matrix<Petsc_Scalar, General, ColMajor> eigen_vec;
      
FindEigenvaluesSlepc(var_eig, lambda, lambda_imag, eigen_vec);

Location :

SetRestartNleigs, GetRestartNleigs

Syntax :

  void SetRestartNleigs(double eps);
  double GetRestartNleigs();

The method SetRestartNleigs provides the restart parameter to be used for Nleigs (it corresponds to the NEPNLEIGSSetRestart Slepc function). The method GetRestartNleigs returns the restart parameter.

Example :

// non-linear eigenproblem
SplitSparseNonLinearEigenProblem<Petsc_Scalar, Symmetric, ArrayRowSymSparse> var_eig;
var_eig.SetStoppingCriterion(1e-12);
var_eig.SetNbAskedEigenvalues(nb_eigenval);
// matrices A_i and polynomial P_i and Q_i are provided
var_eig.InitMatrix(vec_Ai, coef_Pi, coef_Qi);
      
// eigenvalues close to a shift sigma are searched
var_eig.SetTypeSpectrum(var_eig.CENTERED_EIGENVALUES, complex<double>(0.4, -0.1));

// if you want to use Nleigs solver
SlepcParamNep& param = var_eig.GetSlepcParameters();
param.SetEigensolverType(param.NLEIGS);
// for this solver a region [a, b] x [c, d] is required
param.SetIntervalRegion(0.3, 0.6, -0.3, -1e-4);

// you can specify a restart parameter
param.SetRestartNleigs(0.3);

// eigenvalues and eigenvectors are computed by calling Slepc
Vector<Petsc_Scalar> lambda, lambda_imag;
Matrix<Petsc_Scalar, General, ColMajor> eigen_vec;
      
FindEigenvaluesSlepc(var_eig, lambda, lambda_imag, eigen_vec);

Location :

SetRKShifts, GetRKShifts

Syntax :

  void SetRKShifts(const Vector& s);
  const Vector& GetRKShifts();

The method SetRKShifts provides the shifts to be used for Nleigs. The method GetRKShifts returns these shifts..

Location :

GetSlepcParameters

Syntax :

  SlepcParamNep& GetSlepcParameters();

This method returns the Slepc parameters associated with NEP solver.

Example :

  // non-linear eigenproblem
  SplitSparseNonLinearEigenProblem<Petsc_Scalar, Symmetric  , ArrayRowSymSparse> var_eig;

  var_eig.SetStoppingCriterion(1e-14);
  var_eig.SetNbAskedEigenvalues(10);
  // matrices A_i and polynomial P_i and Q_i are provided
  var_eig.InitMatrix(vec_Ai, coef_Pi, coef_Qi);

  complex<double> center(-6.0, 0.5);  
  var_eig.SetTypeSpectrum(var_eig.CENTERED_EIGENVALUES, center);

  // you can change Slepc parameters
  SlepcParamNep& param = var_eig.GetSlepcParameters();
  param.SetEigensolverType(param.SLP);

  // Slepc is called    
  FindEigenvaluesSlepc(var_eig, lambda, lambda_imag, eigen_vec);

Location :

SetExactPreconditioning, ExactPreconditioning

Syntax :

  void SetExactPreconditioning(bool yes = true);
  bool ExactPreconditioning();

This method SetExactPreconditioning informs that the preconditioning used to solve T(λ) is actually exact. This flag is already enabled for the class SplitSparseNonLinearEigenProblem since this classes uses a direct solver as preconditioning. The method ExactPreconditioning returns true if the preconditioning is assumed to be exact.

Location :

GetSingularities, SetSingularities

Syntax :

  void SetSingularities(const Vector& s);
  const Vector& GetSingularities();

The method SetSingularities sets the list of singular points of the operator T (i.e. poles of T). It is not needed to call this function for the class SplitSparseNonLinearEigenProblem since this class is based on the split formulation with rational functions. Slepc computes the singular points directly from these rational functions. The method GetSingularities returns the stored singular points.

Location :

ExplicitMatrix, SetExplicitMatrix

Syntax :

  void SetExplicitMatrix(bool yes = true);
  bool ExplicitMatrix();

The method SetExplicitMatrix tells to store explicitely (with Petsc format) all the involved matrices (instead of using shell matrices). If you activate this flag, you will be able to perform sequential computations only, since this functionality has not been implemented in parallel. The method ExplicitMatrix returns true if you activated this flag.

Location :

SetSplitMatrices, UseSplitMatrices, GetNbSplitMatrices

Syntax :

  void SetSplitMatrices(int n);
  bool UseSplitMatrices();
  int GetNbSplitMatrices();

The method SetSplitMatrices tells to use the split form for defining Slepc eigenrproblem:
T(λ) = Σ f_i(λ) A_i
The argument n is the number of matrices A_i. This method is already automatically called if you are instantiating the class SplitSparseNonLinearEigenProblem. The method UseSplitMatrices returns true if the split form is used. The method GetNbSplitMatrices returns the number of matrices A_i.

Location :

SetNumeratorSplitFct, SetDenominatorSplitFct, GetNumeratorSplitFct, GetDenominatorSplitFct

Syntax :

  void SetNumeratorSplitFct(int i, const Vector&);
  void SetDenominatorSplitFct(int i, const Vector&);
  const Vector& GetNumeratorSplitFct(int i);
  const Vector& GetDenominatorSplitFct(int i);

The method SetNumeratorSplitFct sets the numerator P_i in the split form:
T(λ) = Σ P_i(λ) / Q_i(λ) A_i
while the method SetDenominatorSplitFct sets the denominator Q_i. For the class SplitSparseNonLinearEigenProblem, the numerators and denominators are already set when you call InitMatrix but you can modify them with those methods. The methods GetNumeratorSplitFct and GetDenominatorSplitFct return the numerator P_i and the denominator Q_i respectively.

Example :

  // non-linear eigenproblem
  SplitSparseNonLinearEigenProblem<Petsc_Scalar, Symmetric, ArrayRowSymSparse> var_eig;

  var_eig.SetStoppingCriterion(1e-14);
  var_eig.SetNbAskedEigenvalues(10);
  // matrices A_i and polynomial P_i and Q_i are provided
  var_eig.InitMatrix(vec_Ai, coef_Pi, coef_Qi);

  // if you want to change P_2 for instance
  Vector<Petsc_Scalar> numer(4);
  numer(0) = 1.0; numer(1) = -2.0; numer(2) = 3.0; numer(3) = 4.0; // P_2 = x^3 - 2 x^2 + 3 x + 4
  var_eig.SetNumeratorSplitFct(2, numer);
  
  complex<double> center(-6.0, 0.5);  
  var_eig.SetTypeSpectrum(var_eig.CENTERED_EIGENVALUES, center);

  // Slepc is called    
  FindEigenvaluesSlepc(var_eig, lambda, lambda_imag, eigen_vec);

Location :

CheckValue

Syntax :

  void CheckValue(const T& L);

This method checks if the given value of L is valid (not infinite, or nan). This method is already automatically called when you call FindEigenvaluesSlepc.

Location :

ComputeOperator, MltOperator

Syntax :

  void ComputeOperator(const T& L);
  void MltOperator(const T& L, const SeldonTranspose&, const Vector& X, Vector& Y);

The method ComputeOperator computes the operator T(λ) such that the computation of Y = T(λ) X will be fast (performed by the method MltOperator). These two methods are virtual, such that you can overload them for your own non-linear eigenproblem.

Location :

ComputeOperatorExplicit

Syntax :

  void ComputeOperatorExplicit(const T& L, DistributedMatrix<T, General, ArrayRowSparse>& A);

This method computes explicitly the operator T(λ) and stores it in the matrix A. This method is virtual, such that you can overload them for your own non-linear eigenproblem. This method is called only if matrices are explicitely computed (see method SetExplicitMatrix).

Location :

ComputeJacobian, MltJacobian

Syntax :

  void ComputeJacobian(const T& L);
  void MltJacobian(const T& L, const SeldonTranspose&, const Vector& X, Vector& Y);

The method ComputeJacobian computes the operator T'(λ) such that the computation of Y = T'(λ) X will be fast (performed by the method MltJacobian). These two methods are virtual, such that you can overload them for your own non-linear eigenproblem. T' stands here for the derivative of T with respect to λ, that's why it is labelled as a jacobian.

Location :

ComputeJacobianExplicit

Syntax :

  void ComputeJacobianExplicit(const T& L, DistributedMatrix<T, General, ArrayRowSparse>& A);

This method computes explicitly the operator T'(λ) and stores it in the matrix A. This method is virtual, such that you can overload them for your own non-linear eigenproblem. This method is called only if matrices are explicitely computed (see method SetExplicitMatrix).

Location :

ComputePreconditioning, ApplyPreconditioning

Syntax :

  void ComputePreconditioning(const T& L);
  void ComputePreconditioning(const Vector& vecL, const Vector& coef);
  void ApplyPreconditioning(const SeldonTranspose&, const Vector& X, Vector& Y);

The method ComputePreconditioning computes a preconditioning for solving the operator T(λ) such that the computation of Y = M^-1 X will be fast (performed by the method ApplyPreconditioning). M^-1 stands for the preconditioning matrix. The second syntax is used to compute the preconditioning for solving the operator:
A = Σ coef_i T(λ_i)
These three methods are virtual, such that you can overload them for your own non-linear eigenproblem.

Location :

ComputeSplitPreconditioning

Syntax :

  void ComputeSplitPreconditioning(const Vector& numL, const Vector& coef);

The method ComputeSplitPreconditioning computes a preconditioning for solving the operator:
A = Σ coef_i A_numL(i)
This method is virtual, such that you can overload it for your own non-linear eigenproblem.

Location :

ComputeExplicitPreconditioning

Syntax :

  void ComputeExplicitPreconditioning(DistributedMatrix<T, General, ArrayRowSparse>& A);

This method computes a preconditioning for the provided matrix A. This method is called only if matrices are explicitely computed (see method SetExplicitMatrix). This method is virtual, such that you can overload it for your own non-linear eigenproblem.

Location :

ComputeOperatorSplitExplicit

Syntax :

  void ComputeOperatorSplitExplicit(int i, DistributedMatrix<T, General, ArrayRowSparse>& A);

This method computes explicitely the matrix A_i and stores it in matrix A. This method is called only if matrices are explicitely computed (see method SetExplicitMatrix). This method is virtual, such that you can overload it for your own non-linear eigenproblem.

Location :

ComputeOperatorSplitExplicit

Syntax :

  void MltOperatorSplit(int i, const SeldonTranspose&, const Vector& X, Vector& Y);

This method computes Y = A_i X if a split form of T has been set. This method is virtual, such that you can overload it for your own non-linear eigenproblem.

Location :

InitMatrix

Syntax :

  void InitMatrix(Vector<DistributedMatrix<T, Prop, Storage> >&,
                  const Vector<Vector<T> >& numer, const Vector<Vector<T>& denom);

This method initializes the matrices A_i and polynomials P_i and Q_i. We remind that the following split form is used in this class
T(λ) = Σ P_i(λ) / Q_i(λ) A_i

Example:

// we want to solve a rational eigenvalue problem
// T(lambda) x = 0
// where T(lambda) = \sum_i P_i(lambda) / Q_i(lambda) A_i
// matrices A_i are read from a file
Vector<DistributedMatrix<Petsc_Scalar, Symmetric, ArrayRowSymSparse> > vec_Ai(3);
vec_Ai(0).ReadText("A0.dat");
vec_Ai(1).ReadText("A1.dat");
vec_Ai(2).ReadText("A2.dat");

// numerators P_i and denominators Q_i can be set manually
Petsc_Scalar a = 1.0, b = 2.0; c = 3.0;
Vector<Vector<Petsc_Scalar> > coef_Pi(3), coef_Qi(3);
// here P_0(lambda) = a lambda^2 + b lambda + c
coef_Pi(0).Reallocate(3); coef_Pi(0)(0) = a; coef_Pi(0)(1) = b; coef_Pi(0)(2) = c;
// Q_0(lambda) = 2 lambda - 5;
coef_Qi(0).Reallocate(2); coef_Qi(0)(0) = 2.0; coef_Qi(0)(1) = -5.0;
// ...

// main object handling sparse rational eigenvalue problem
SplitSparseNonLinearEigenProblem<Petsc_Scalar, Symmetric, ArrayRowSymSparse> var_eig;
var_eig.SetStoppingCriterion(1e-12);
var_eig.SetNbAskedEigenvalues(nb_eigenval);
// matrices A_i and polynomial P_i and Q_i are provided
var_eig.InitMatrix(vec_Ai, coef_Pi, coef_Qi);
      
// eigenvalues close to a shift sigma are searched
var_eig.SetTypeSpectrum(var_eig.CENTERED_EIGENVALUES, complex<double>(-6.0, 2.0));

// if you want to use Nleigs solver
SlepcParamNep& param = var_eig.GetSlepcParameters();
param.SetEigensolverType(param.NLEIGS);
// for this solver a region [a, b] x [c, d] is required
param.SetIntervalRegion(0.3, 0.6, -0.3, -1e-4);
  
// eigenvalues and eigenvectors are computed by calling Slepc
Vector<Petsc_Scalar> lambda, lambda_imag;
Matrix<Petsc_Scalar, General, ColMajor> eigen_vec;
      
FindEigenvaluesSlepc(var_eig, lambda, lambda_imag, eigen_vec);

Location :

GetDirectSolver

Syntax :

  SparseDistributedSolver<T>& GetDirectSolver();

This method gives access to the direct solver stored in the object.

Example:

// main object handling sparse rational eigenvalue problem
SplitSparseNonLinearEigenProblem<Petsc_Scalar, Symmetric, ArrayRowSymSparse> var_eig;
var_eig.SetStoppingCriterion(1e-12);
var_eig.SetNbAskedEigenvalues(nb_eigenval);
// matrices A_i and polynomial P_i and Q_i are provided
var_eig.InitMatrix(vec_Ai, coef_Pi, coef_Qi);
      
// eigenvalues close to a shift sigma are searched
var_eig.SetTypeSpectrum(var_eig.CENTERED_EIGENVALUES, complex<double>(-6.0, 2.0));

// if you want to use Nleigs solver
SlepcParamNep& param = var_eig.GetSlepcParameters();
param.SetEigensolverType(param.NLEIGS);
// for this solver a region [a, b] x [c, d] is required
param.SetIntervalRegion(0.3, 0.6, -0.3, -1e-4);

// if you want to parameter the direct solver
SparseDistributedSolver<Petsc_Scalar>& solver = var_eig.GetDirectSolver();
solver.RefineSolution();

// eigenvalues and eigenvectors are computed by calling Slepc
Vector<Petsc_Scalar> lambda, lambda_imag;
Matrix<Petsc_Scalar, General, ColMajor> eigen_vec;
      
FindEigenvaluesSlepc(var_eig, lambda, lambda_imag, eigen_vec);

Location :