Algebraic solutions
Definitions
Matrices
Matrix Definition A matrix is a linear transformation between finite dimensional vector spaces.
Assembling a matrix Assembling a matrix means defining its action as entries stored in a sparse or dense format. For example, in the finite element context, the storage format is sparse to take advantage of the many zero entries.
 Symmetric matrix

\$A = A^T\$
 Definite (resp. semidefinite) positive matrix

All eigenvalue are

\$>0\$ (resp \$\geq 0\$) or

\$x^T A x > 0, \forall\ x\$ (resp. \$x^T\ A\ x\geq 0\, \forall\ x\$)

 Definite (resp. seminegative) matrix

All eigenvalue are

\$<0\$ (resp. \$\leq 0\$) or

\$x^T\ A\ x < 0\ \forall\ x\$ (resp. \$x^T\ A\ x \leq 0\, \forall\ x\$)

 Indefinite matrix

There exists

positive and negative eigenvalue (Stokes, Helmholtz) or

there exists \$x,y\$ such that \$x^TAx > 0 > y^T A y\$

Preconditioners
Definition
Let \$A\$ be a \$\mathbb{R}^{n\times n}\$ matrix, \$x\$ and \$b\$ be \$\mathbb{R}^n\$ vectors, we wish to solve \$A x = b\$.
 Definition

A preconditioner \$\mathcal{P}\$ is a method for constructing a matrix (just a linear function, not assembled!) P^{1} = \mathcal{P}(A,A_p) using a matrix \$A\$ and extra information \$A_p\$, such that the spectrum of P^{1}A (left preconditioning) or A P^{1} (right preconditioning) is wellbehaved. The action of preconditioning improves the conditioning of the previous linear system.
Left preconditioning: We solve for (P^{1} A) x = P^{1} b and we build the Krylov space \{ P^{1} b, (P^{1}A) P^{1} b, (P^{1}A)^2 P^{1} b, \dots\}
Right preconditioning: We solve for (A P^{1}) P x = b and we build the Krylov space \{ b, (P^{1}A)b, (P^{1}A)^2b, \dotsc \}
Note that the product P^{1}A or A P^{1} is never assembled.
Properties
Let us now describe some properties of preconditioners

P^{1} is dense, P is often not available and is not needed

A is rarely used by \mathcal{P}, but A_p = A is common

A_p is often a sparse matrix, the \e preconditioning \e matrix
Here are some numerical methods to solve the system A x = b

Matrixbased: Jacobi, GaussSeidel, SOR, ILU(k), LU

Parallel: BlockJacobi, Schwarz, Multigrid, FETIDP, BDDC

Indefinite: Schurcomplement, Domain Decomposition, Multigrid
Preconditioner strategies
 Relaxation

Split into lower, diagonal, upper parts: \$A = L + D + U\$.
 Jacobi

Cheapest preconditioner: \$P^{1}=D^{1}\$.
# sequential
pctype=jacobi
# parallel
pctype=block_jacobi
 Successive overrelaxation (SOR)

Implemented as a sweep.

\$\omega = 1\$ corresponds to GaussSeidel.

Very effective at removing highfrequency components of residual.
# sequential
pctype=sor
Factorization
Two phases

symbolic factorization: find where fill occurs, only uses sparsity pattern.

numeric factorization: compute factors.
 LU decomposition


preconditioner.

Expensive, for \$m\times m\$ sparse matrix with bandwidth \$b\$, traditionally requires \$\mathcal{O}(mb^2)\$ time and \$\mathcal{O}(mb)\$ space.

Bandwidth scales as \$m^{\frac{d1}{d}}\$ in \$d\$dimensions.

Optimal in 2D: \$\mathcal{O}(m \cdot \log m)\$ space, \$\mathcal{O}(m^{3/2})\$ time.

Optimal in 3D: \$\mathcal{O}(m^{4/3})\$ space, \$\mathcal{O}(m^2)\$ time.



Symbolic factorization is problematic in parallel.

Incomplete LU

Allow a limited number of levels of fill: ILU(\$k\$).

Only allow fill for entries that exceed threshold: ILUT.

Usually poor scaling in parallel.

No guarantees.
1level Domain decomposition
Domain size \$L\$, subdomain size \$H\$, element size \$h\$

Overlapping/Schwarz

Solve Dirichlet problems on overlapping subdomains.

No overlap: \$\textit{its} \in \mathcal{O}\big( \frac{L}{\sqrt{Hh}} \big)\$.

Overlap \delta: \textit{its} \in \big( \frac L {\sqrt{H\delta}} \big).

pctype=gasm # has a coarse grid preconditioner
pctype=asm

NeumannNeumann

Solve Neumann problems on nonoverlapping subdomains.

\textit{its} \in \mathcal{O}\big( \frac{L}{H}(1+\log\frac H h) \big).

Tricky null space issues (floating subdomains).

Need subdomain matrices, not globally assembled matrix.

Notes: Multilevel variants knock off the leading \frac L H.
Both overlapping and nonoverlapping with this bound.

BDDC and FETIDP

Neumann problems on subdomains with coarse grid correction.

\textit{its} \in \mathcal{O}\big(1 + \log\frac H h \big).

Multigrid
Hierarchy: Interpolation and restriction operators \Pi^\uparrow : X_{\text{coarse}} \to X_{\text{fine}} \qquad \Pi^\downarrow : X_{\text{fine}} \to X_{\text{coarse}}

Geometric: define problem on multiple levels, use grid to compute hierarchy.

Algebraic: define problem only on finest level, use matrix structure to build hierarchy.
Galerkin approximation
Assemble this matrix: A_{\text{coarse}} = \Pi^\downarrow A_{\text{fine}} \Pi^\uparrow
Application of multigrid preconditioner ( V cycle)

Apply presmoother on fine level (any preconditioner).

Restrict residual to coarse level with \Pi^\downarrow.

Solve on coarse level A_{\text{coarse}} x = r.

Interpolate result back to fine level with \Pi^\uparrow.

Apply postsmoother on fine level (any preconditioner).
Multigrid convergence properties

Textbook: P^{1}A is spectrally equivalent to identity

Constant number of iterations to converge up to discretization error.


Most theory applies to SPD systems

variable coefficients (e.g. discontinuous): low energy interpolants.

mesh and/or physicsinduced anisotropy: semicoarsening/line smoothers.

complex geometry: difficult to have meaningful coarse levels.


Deeper algorithmic difficulties

nonsymmetric (e.g. advection, shallow water, Euler).

indefinite (e.g. incompressible flow, Helmholtz).


Performance considerations

Aggressive coarsening is critical in parallel.

Most theory uses SOR smoothers, ILU often more robust.

Coarsest level usually solved semiredundantly with direct solver.


Multilevel Schwarz is essentially the same with different language

assume strong smoothers, emphasize aggressive coarsening.

List of PETSc Preconditioners
See this PETSc page for a complete list.
PETSc 
Description 
Parallel 
none 
No preconditioner 
yes 
jacobi 
diagonal preconditioner 
yes 
bjacobi 
block diagonal preconditioner 
yes 
sor 
SOR preconditioner 
yes 
lu 
Direct solver as preconditioner 
depends on the factorization package (e.g.mumps,pastix: OK) 
shell 
User defined preconditioner 
depends on the user preconditioner 
mg 
multigrid prec 
yes 
ilu 
incomplete lu 

icc 
incomplete cholesky 

cholesky 
Cholesky factorisation 
yes 
asm 
Additive Schwarz Method 
yes 
gasm 
Scalable Additive Schwarz Method 
yes 
ksp 
Krylov subspace preconditioner 
yes 
fieldsplit 
block preconditioner framework 
yes 
lsc 
Least Square Commutator 
yes 
gamg 
Scalable Algebraic Multigrid 
yes 
hypre 
Hypre framework (multigrid…) 

bddc 
balancing domain decomposition by constraints preconditioner 
yes 
Principles
Feel++ abstracts the PETSc library and provides a subset (sufficient in most cases) to the PETSc features. It interfaces with the following PETSc libraries: Mat
, Vec
, KSP
, PC
, SNES.

Vec
Vector handling library 
Mat
Matrix handling library 
KSP
Krylov SubSpace library implements various iterative solvers 
PC
Preconditioner library implements various preconditioning strategies 
SNES
Nonlinear solver library implements various nonlinear solve strategies
All linear algebra are encapsulated within backends using the command line option backend=<backend>
or config file option backend=<backend>
which provide interface to several libraries
Library 
Format 
Backend 
PETSc 
sparse 

Eigen 
sparse 

Eigen 
dense 

The default backend
is petsc.
Somes generic examples
The configuration files .cfg
allow for a wide range of options to solve a linear or nonlinear system.
We consider now the following example [import:"marker1"](../../codes/mylaplacian.cpp)
To execute this example
# sequential
./feelpp_tut_laplacian
# parallel on 4 cores
mpirun np 4 ./feelpp_tut_laplacian
As described in section
Direct solver
cholesky
and lu
factorisation are available. However the parallel implementation depends on the availability of MUMPS. The configuration is very simple.
# no iterative solver
ksptype=preonly
#
pctype=cholesky
Using the PETSc backend allows to choose different packages to compute the factorization.
Package 
Description 
Parallel 

PETSc own implementation 
yes 

MUltifrontal Massively Parallel sparse direct Solver 
yes 

Unsymmetric MultiFrontal package 
no 

Parallel Sparse matriX package 
yes 
To choose between these factorization package
# choose mumps
pcfactormatsolverpackage=mumps
# choose umfpack (sequential)
pcfactormatsolverpackage=umfpack
In order to perform a cholesky type of factorisation, it is required to set the underlying matrix to be SPD.
// matrix
auto A = backend>newMatrix(_test=...,_trial=...,_properties=SPD);
// bilinear form
auto a = form2( _test=..., _trial=..., _properties=SPD );
Using iterative solvers
Using CG and ICC(3)
with a relative tolerance of 1e12:
ksprtol=1.e12
ksptype=cg
pctype=icc
pcfactorlevels=3
Using GMRES and ILU(3)
with a relative tolerance of 1e12 and a restart of 300:
ksprtol=1.e12
ksptype=gmres
kspgmresrestart=300
pctype=ilu
pcfactorlevels=3
Using GMRES and Jacobi
With a relative tolerance of 1e12 and a restart of 100:
ksprtol=1.e12
ksptype=gmres
kspgmresrestart 100
pctype=jacobi
Monitoring linear nonlinear and eigen problem solver residuals
# linear
ksp_monitor=1
# nonlinear
snesmonitor=1
# eigen value problem
epsmonitor=1
Solving the Laplace problem
We start with the quickstart Laplacian example, recall that we wish to, given a domain \Omega, find u such that
\nabla \cdot (k \nabla u) = f \mbox{ in } \Omega \subset \mathbb{R}^{2},\\ u = g \mbox{ on } \partial \Omega
Monitoring KSP solvers
feelpp_qs_laplacian kspmonitor=true
Viewing KSP solvers
shell> mpirun np 2 feelpp_qs_laplacian kspmonitor=1 kspview=1
0 KSP Residual norm 8.953261456448e01
1 KSP Residual norm 7.204431786960e16
KSP Object: 2 MPI processes
type: gmres
GMRES: restart=30, using Classical (unmodified) GramSchmidt
Orthogonalization with no iterative refinement
GMRES: happy breakdown tolerance 1e30
maximum iterations=1000
tolerances: relative=1e13, absolute=1e50, divergence=100000
left preconditioning
using nonzero initial guess
using PRECONDITIONED norm type for convergence test
PC Object: 2 MPI processes
type: shell
Shell:
linear system matrix = precond matrix:
Matrix Object: 2 MPI processes
type: mpiaij
rows=525, cols=525
total: nonzeros=5727, allocated nonzeros=5727
total number of mallocs used during MatSetValues calls =0
not using Inode (on process 0) routines
Solvers and preconditioners
You can now change the Krylov subspace solver using the ksptype
option and the preconditioner using pcptype
option.
For example,

to solve use the conjugate gradient,
cg
, solver and the default preconditioner use the following
./feelpp_qs_laplacian ksptype=cg kspview=1 kspmonitor=1

to solve using the algebraic multigrid preconditioner,
gamg
, withcg
as a solver use the following
./feelpp_qs_laplacian ksptype=cg kspview=1 kspmonitor=1 pctype=gamg
Block factorisation
Stokes
We now turn to the quickstart Stokes example, recall that we wish to, given a domain \Omega, find (\mathbf{u},p) such that
This problem is indefinite. Possible solution strategies are

Uzawa,

penalty(techniques from optimisation),

augmented lagrangian approach (Glowinski,Le Tallec)
Note that The Infsup condition must be satisfied. In particular for a multigrid strategy, the smoother needs to preserve it.
General approach for saddle point problems
The Krylov subspace solvers for indefinite problems are MINRES, GMRES. As to preconditioning, we look first at the saddle point matrix M and its block factorization M = LDL^T, indeed we have :

Elman, Silvester and Wathen propose 3 preconditioners:
where \$\tilde{S} \approx S^{1} = B^T A^{1} B\$ and \$\tilde{A}^{1} \approx A^{1}\$
Preconditioner strategies
Relaxation
Split into lower, diagonal, upper parts: \$A = L + D + U\$.
Jacobi
Cheapest preconditioner: \$P^{1}=D^{1}\$.
# sequential
pctype=jacobi
# parallel
pctype=block_jacobi
Successive overrelaxation (SOR)

Implemented as a sweep.

\$\omega = 1\$ corresponds to GaussSeidel.

Very effective at removing highfrequency components of residual.
# sequential
pctype=sor
Factorization
Two phases

symbolic factorization: find where fill occurs, only uses sparsity pattern.

numeric factorization: compute factors.
LU decomposition

preconditioner.

Expensive, for \$m\times m\$ sparse matrix with bandwidth \$b\$, traditionally requires \$\mathcal{O}(mb^2)\$ time and \$\mathcal{O}(mb)\$ space.

Bandwidth scales as \$m^{\frac{d1}{d}}\$ in ddimensions.

Optimal in 2D: \$\mathcal{O}(m \cdot \log m)\$ space, \$\mathcal{O}(m^{3/2})\$ time.

Optimal in 3D: \$\mathcal{O}(m^{4/3})\$ space, \$\mathcal{O}(m^2)\$ time.


Symbolic factorization is problematic in parallel.
Incomplete LU

Allow a limited number of levels of fill: ILU(\$k\$).

Only allow fill for entries that exceed threshold: ILUT.

Usually poor scaling in parallel.

No guarantees.
1level Domain decomposition
Domain size $$L$$, subdomain size $$H$$, element size $$h$$

Overlapping/Schwarz

Solve Dirichlet problems on overlapping subdomains.

No overlap: \$\textit{its} \in \mathcal{O}\big( \frac{L}{\sqrt{Hh}} \big)\$.

Overlap \$\delta: \textit{its} \in \big( \frac L {\sqrt{H\delta}} \big)\$.


NeumannNeumann

Solve Neumann problems on nonoverlapping subdomains.

\$\textit{its} \in \mathcal{O}\big( \frac{L}{H}(1+\log\frac H h) \big)\$.

Tricky null space issues (floating subdomains).

Need subdomain matrices, not globally assembled matrix.

Multilevel variants knock off the leading \$\frac L H\$. Both overlapping and nonoverlapping with this bound. 

BDDC and FETIDP

Neumann problems on subdomains with coarse grid correction.

\$\textit{its} \in \mathcal{O}\big(1 + \log\frac H h \big)\$.

Multigrid
Introduction
Hierarchy: Interpolation and restriction operators \$\Pi^\uparrow : X_{\text{coarse}} \to X_{\text{fine}} \qquad \Pi^\downarrow : X_{\text{fine}} \to X_{\text{coarse}} \$

Geometric: define problem on multiple levels, use grid to compute hierarchy.

Algebraic: define problem only on finest level, use matrix structure to build hierarchy.
Galerkin approximation
Assemble this matrix: \$A_{\text{coarse}} = \Pi^\downarrow A_{\text{fine}} \Pi^\uparrow\$
Application of multigrid preconditioner (\$V\$cycle)

Apply presmoother on fine level (any preconditioner).

Restrict residual to coarse level with \$\Pi^\downarrow\$.

Solve on coarse level \$A_{\text{coarse}} x = r\$.

Interpolate result back to fine level with \Pi^\uparrow.

Apply postsmoother on fine level (any preconditioner).
Multigrid convergence properties

Textbook: \$P^{1}A\$ is spectrally equivalent to identity

Constant number of iterations to converge up to discretization error.


Most theory applies to SPD systems

variable coefficients (e.g. discontinuous): low energy interpolants.

mesh and/or physicsinduced anisotropy: semicoarsening/line smoothers.

complex geometry: difficult to have meaningful coarse levels.


Deeper algorithmic difficulties

nonsymmetric (e.g. advection, shallow water, Euler).

indefinite (e.g. incompressible flow, Helmholtz).


Performance considerations

Aggressive coarsening is critical in parallel.

Most theory uses SOR smoothers, ILU often more robust.

Coarsest level usually solved semiredundantly with direct solver.


Multilevel Schwarz is essentially the same with different language

assume strong smoothers, emphasize aggressive coarsening.

List of PETSc Preconditioners
See this PETSc page for a complete list.
PETSc 
Description 
Parallel 
none 
No preconditioner 
yes 
jacobi 
diagonal preconditioner 
yes 
bjacobi 
block diagonal preconditioner 
yes 
sor 
SOR preconditioner 
yes 
lu 
Direct solver as preconditioner 
depends on the factorization package (e.g.mumps,pastix: OK) 
shell 
User defined preconditioner 
depends on the user preconditioner 
mg 
multigrid prec 
yes 
ilu 
incomplete lu 

icc 
incomplete cholesky 

cholesky 
Cholesky factorisation 
yes 
asm 
Additive Schwarz Method 
yes 
gasm 
Scalable Additive Schwarz Method 
yes 
ksp 
Krylov subspace preconditioner 
yes 
fieldsplit 
block preconditioner framework 
yes 
lsc 
Least Square Commutator 
yes 
gamg 
Scalable Algebraic Multigrid 
yes 
hypre 
Hypre framework (multigrid…) 

bddc 
balancing domain decomposition by constraints preconditioner 
yes 
Algebra Backends
NonLinear algebra backends are crucial components of Feel++. They provide a uniform interface between Feel++ data structures and underlying the linear algebra libraries used by Feel++.
Libraries
Feel++ interfaces the following libraries:

PETSc : Portable, Extensible Toolkit for Scientific Computation

SLEPc : Eigen value solver framework based on PETSc

Eigen3
Backend
Backend is a template class parametrized by the numerical type providing access to

vector

matrix : dense and sparse

algorithms : solvers, preconditioners, …
PETSc provides sequential and parallel data structures whereas Eigen3 is only sequential.
To create a Backend, use the free function backend(…)
which has the following interface:
backend(_name="name_of_backend",
_rebuild=... /// truefalse,
_kind=..., // type of backend,
_worldcomm=... // communicator
)
All these parameters are optional which means that the simplest call reads for example:
auto b = backend();
They take default values as described in the following table:
parameter 
type 
default value 

string 
"" (empty string ) 

boolean 
false 

string 
"petsc" 

WorldComm 
Environment::worldComm() 
_name
Backends are store in a Backend factory and handled by a manager that
allows to keep track of allocated backends. They a registered with
respect to their name and kind. The default name value is en empty
string (""
) which names the default Backend. The _name parameter is
important because it provides also the name for the command
line/config file option section associated to the associated Backend.
Only the command line/config file options for the default backend are registered. Developers have to register the option for each new Backend they define otherwise failures at runtime are to be expected whenever a Backend command line option is accessed.
Consider that you create a Backend name projection
in your code like
this
auto b = backend(_name="projection");
to register the command line options of this Backend
Environment env( _argc=argc, _argv=argv,
_desc=backend_options("projection") );
_kind
Feel++ supports three kind of Backends:

petsc : PETSC Backend

eigen_dense

eigen_sparse
SLEPc uses the PETSc Backend since it is based on PETSc. 
The kind of Backend can be changed from the command line or configuration file thanks to the "backend" option.
auto b = backend(_name="name",
_kind=soption(_prefix="name",_name="backend"))
and in the config file
[name]
backend=petsc
backend=eigen_sparse
_rebuild
If you want to reuse a Backend and not allocate a new one plus add the corresponding option to the command line/configuration file, you can rebuild the Backend in order to delete the data structures already associated to this Backend and in particular the preconditioner. It is mandatory to do that when you solve say a linear system first with dimensions \$m\times m\$ and then solve another one with different dimension \$n \times n\$ because in that case the Backend will throw an error saying that the dimensions are incompatible. To avoid that you have to rebuild the Backend.
auto b = backend(_name="mybackend");
// solve A x = f
b>solve(_solution=x,_matrix=A,_rhs=f);
// rebuild: clean up the internal Backend data structure
b = backend(_name="mybackend",_rebuild=true);
// solve A1 x1 = f1
b>solve(_solution=x1,_matrix=A1,_rhs=f1);
Although this feature might appear useful, you have to make sure that the solving strategy applies to all problems because you won’t be able to customize the solver/preconditioner for each problem. If you have different problems to solve and need to have custom solver/preconditioner it would be best to have different Backends. 
_worldComm
One of the strength of Feel++ is to be able to change the communicator and in the case of Feel++ the WorldComm. This allows for example to

solve sequential problems

solve a problem on a subset of MPI processes
For example passing a sequential WorldComm to backend()
would then in the subsequent use of the Backend generate sequential data structures (e.g. IndexSet, Vector and Matrix) and algorithms (e.g. Krylov Solvers).
// create a sequential Backend
auto b = backend(_name="seq",
_worldComm=Environment::worldCommSeq());
auto A = b>newMatrix(); // sequential Matrix
auto f = b>newVector(); // sequential Vector
Info The default WorldComm is provided by
Environment::worldComm()
and it conresponds to the default MPI communicatorMPI_COMM_WORLD
.
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 lowdimensional subspace, this is called a RayleighRitz 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.
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
.
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 :
Problem type  EigenProblemType  command line key 

Standard Hermitian 
HEP 
"hep" 
Standard nonHermitian 
NHEP 
"nhep" 
Generalized Hermitian 
GHEP 
"ghep" 
Generalized nonHermitian 
GNHEP 
"gnhep" 
Positive definite Generalized nonHermitian 
PGNHEP 
"pgnhep" 
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 :
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" 
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
.
Spectral transformation  SpectralTransformationType  command line key 

Shift 
SHIFT 
shift 
Shift and invert 
SINVERT 
shift_invert 
Cayley 
CAYLEY 
cayley 
Eigensolvers
The details of the implementation of the different solvers can be found in the SLEPc Technical Reports.
The default solver is KrylovSchur, but can be modified using
EigenSolverType
or the option solvereigen.solver
.
Solver  EigenSolverType  command line key 

Power 
POWER 
power 
Lapack 
LAPACK 
lapack 
Subspace 
SUBSPACE 
subspace 
Arnoldi 
Arnoldi 
arnoldi 
Lanczos 
LANCZOS 
lanczos 
KrylovSchur 
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.
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 
KrylovSchur 
any 
any 
Arpack 
any 
any 
Special cases of spectrum
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>
.
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.intervala
to set the beginning
of the interval and solvereigen.intervalb
to set the end.
In this case, be aware that the problem need to be generalized and
hermitian. The solver will be set to KrylovSchur 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.