Eigen Problem

To solve standard and generalized eigenvalue problems, Feel++ interfaces SLEPc. SLEPc is a library which extends PETSc to provide the functionality necessary for the solution of eigenvalue problems. It comes with many strategies for both standard and generalized problems, Hermitian or not.

We want to find \((\lambda_i,x_i)\) such that \(Ax = \lambda x\). To do that, most eigensolvers project the problem onto a low-dimensional subspace, this is called a Rayleigh-Ritz projection. + Let \(V_j=[v_1,v_2,...,v_j\)] be an orthogonal basis of this subspace, then the projected problem reads:

Find \((\theta_i,s_i)\) for \(i=1,\ldots,j\) such that \(B_j s_i=\theta_i s_i\) where \(B_j=V_j^T A V_j\).

Then the approximate eigenpairs \((\lambda_i,x_i)\) of the original problem are obtained as: \(\lambda_i=\theta_i\) and \(x_i=V_j s_i\).

The eigensolvers differ from each other in the way the subspace is built.

1. Code

In Feel++, there is two functions that can be used to solve this type of problems, eigs and veigs.

Here is an example of how to use veigs.

auto Vh = Pch<Order>( mesh );
auto a = form2( _test=Vh, _trial=Vh );
// fill a
auto b = form2( _test=Vh, _trial=Vh );
// fill b
auto eigenmodes = veigs( _formA=a, _formB=b );

where eigenmodes is a std::vector<std::pair<value_type, element_type> > with value_type the type of the eigenvalue, and element_type the type of the eigenvector, a function of the space Vh.

The eigs function does not take the bilinear forms but two matrices. Also, the solver used, the type of the problem, the position of the spectrum and the spectral transformation are not read from the options.

auto Vh = Pch<Order>( mesh );
auto a = form2( _test=Vh, _trial=Vh );
// fill a
auto matA = a.matrixPtr();
auto b = form2( _test=Vh, _trial=Vh );
// fill b
auto matB = b.matrixPtr();
auto eigenmodes = eigs( _matrixA=aHat,
                        _matrixB=bHat,
                        _solver=(EigenSolverType)EigenMap[soption("solvereigen.solver")],
                        _problem=(EigenProblemType)EigenMap[soption("solvereigen.problem")],
                        _transform=(SpectralTransformType)EigenMap[soption("solvereigen.transform")],
                        _spectrum=(PositionOfSpectrum)EigenMap[soption("solvereigen.spectrum")]
                         );
auto femodes = std::vector<decltype(Vh->element())>( eigenmodes.size(), Vh->element() );
int i = 0;
for( auto const& mode : modes )
    femodes[i++] = *mode.second.get<2>();

where eigenmodes is a std::map<real_type, eigenmodes_type> with real_type of the magnitude of the eigenvalue. And eigenmodes_type is a boost::tuple<real_type, real_type, vector_ptrtype> with the first real_type representing the real part of the eigenvalue, the second real_type the imaginary part and the vector_ptrtype is a vector but not an element of a functionspace.

The two functions take a parameter _nev that tel how many eigenpair to compute. This can be set from the command line option --solvereigen.nev. + Another important parameter is _ncv which is the size of the subspace, j above. This parameter should always be greater than nev. SLEPc recommends to set it to max(nev+15, 2*nev). This can be set from the command line option --solvereigen.ncv.

2. Problem type

The standard formulation reads :

Find \(\lambda\in \mathbb{R}\) such that \(Ax = \lambda x\)

where \(\lambda\) is an eigenvalue and \(x\) an eigenvector.

But in the case of the finite element method, we will deal with the generalized form :

Find \(\lambda\in\mathbb{R}\) such that \(Ax = \lambda Bx\)

A standard problem is Hermitian if the matrix A is Hermitian (\(A=A^*\)). + A generalized problem is Hermitian if the matrices \(A\) and \(B\) are Hermitian and if \(B\) is positive definite. + If the problem is Hermitian, then the eigenvalues are real. A special case of the generalized problem is when the matrices are not Hermitian but \(B\) is positive definite.

The type of the problem can be specified using the EigenProblemType, or at run time with the command line option --solvereigen.problem and the following value :

Table 1. Table of problem type
Problem type	EigenProblemType	command line key
Standard Hermitian	HEP	"hep"
Standard non-Hermitian	NHEP	"nhep"
Generalized Hermitian	GHEP	"ghep"
Generalized non-Hermitian	GNHEP	"gnhep"
Positive definite Generalized non-Hermitian	PGNHEP	"pgnhep"

3. Position of spectrum

You can choose which eigenpairs will be computed. The user can set it programmatically with PositionOfSpectrum or at run time with the command line option --solvereigen.spectrum and the following value :

Table 2. Table of position of spectrum
Position of spectrum	PositionOfSpectrum	command line key
Largest magnitude	LARGEST_MAGNITUDE	"largest_magnitude"
Smallest magnitude	SMALLEST_MAGNITUDE	"smallest_magnitude"
Largest real	LARGEST_REAL	"largest_real"
Smallest real	SMALLEST_REAL	"smallest_real"
Largest imaginary	LARGEST_IMAGINARY	"largest_imaginary"
Smallest imaginary	SMALLEST_IMAGINARY	"smallest_imaginary"

4. Spectral transformation

It is observed that the algorithms used to solve the eigenvalue problems find solutions at the extremities of the spectrum. To improve the convergence, one need to compute the eigenpairs of a transformed operator. Those spectral transformations allow to compute solutions that are not on the boundary of the spectrum.

There are 3 types of spectral transformation:

Shift: \(A-\sigma I\) or \(B^{-1}A-\sigma I\)
Shift and invert: \((A-\sigma I)^{-1}\) or \((A-\sigma B)^{-1}B\)
Cayley: \((A-\sigma I)^{-1}(A+\nu I)\) or \((A-\sigma B)^{-1}(A+\nu B)\)

By default, shift and invert is used. You can change it with --solvereigen.transform.

Table 3. Table of spectral transformation
Spectral transformation	SpectralTransformationType	command line key
Shift	SHIFT	shift
Shift and invert	SINVERT	shift_invert
Cayley	CAYLEY	cayley

5. Eigensolvers

The details of the implementation of the different solvers can be found in the SLEPc Technical Reports.

The default solver is Krylov-Schur, but can be modified using EigenSolverType or the option --solvereigen.solver.

Table 4. Table of eigensolver
Solver	EigenSolverType	command line key
Power	POWER	power
Lapack	LAPACK	lapack
Subspace	SUBSPACE	subspace
Arnoldi	Arnoldi	arnoldi
Lanczos	LANCZOS	lanczos
Krylov-Schur	KRYLOVSCHUR	krylovschur
Arpack	ARPACK	arpack

Be careful that all solvers can not compute all the problem types and positions of the spectrum. The possibilities are summarize in the following table.

Table 5. Supported problem type for the eigensolvers
Solver	Position of spectrum	Problem type
Power	Largest magnitude	any
Lapack	any	any
Subspace	Largest magnitude	any
Arnoldi	any	any
Lanczos	any	standard and generalized Hermitian
Krylov-Schur	any	any
Arpack	any	any

6. Special cases of spectrum

6.1. Computing a large portion of the spectrum

In the case where you want compute a large number of eigenpairs, the rule for ncv implies a huge amount of memory to be used. To improve the performance, you can set the mpd parameter, which will limit the dimension of the projected problem.

You can set it via the command line with --solvereigen.mpd <mpd>.

6.2. Computing all the eigenpairs in an interval

If you want to compute all the eigenpairs in a given interval, you need to use the option --solvereigen.interval-a to set the beginning of the interval and --solvereigen.interval-b to set the end.

In this case, be aware that the problem need to be generalized and hermitian. The solver will be set to Krylov-Schur and the transformation to shift and invert. Beside, you’ll need to use a linear solver that will compute the inertia of the matrix, this is set to Cholesky, with mumps if you can use it. + For now, this method is only implemented in the eigs function.