ST: Spectral Transformation

The Spectral Transformation (ST) is the object that encapsulates the functionality required for acceleration techniques based on the transformation of the spectrum. Most eigensolvers in EPS work by applying an operator to a set of vectors and this operator can adopt different forms. The ST object handles all the different possibilities in a uniform way, so that the solver can proceed without knowing which transformation has been selected. The spectral transformation can be specified at run time, as well as related options such as which linear solver to use.

Despite being a rather unrelated concept, ST is also used to handle the preconditioners and correction-equation solvers used in preconditioned eigensolvers such as GD and JD.

The description in this chapter focuses on the use of ST in the context of EPS. For usage within other solver classes, we will provide further details, e.g., shift-and-invert for polynomial eigenproblems in section Spectral Transformation.

General Description

Spectral transformations are powerful tools for adjusting the way in which eigensolvers behave when coping with a problem. The general strategy consists in transforming the original problem into a new one in which eigenvalues are mapped to a new position while eigenvectors remain unchanged. These transformations can be used with several goals in mind:

• Compute internal eigenvalues. In some applications, the eigenpairs of interest are not the extreme ones (largest magnitude, rightmost, leftmost), but those contained in a certain interval or those closest to a certain value of the complex plane.

• Accelerate convergence. Convergence properties typically depend on how close the eigenvalues are from each other. With some spectral transformations, difficult eigenvalue distributions can be remapped in a more favorable way in terms of convergence.

• Handle some special situations. For instance, in generalized problems when matrix \(B\) is singular, it may be necessary to use a spectral transformation.

SLEPc separates spectral transformations from solution methods so that any combination of them can be specified by the user. To achieve this, most eigensolvers contained in EPS are implemented in such a way that they are independent of which transformation has been selected by the user (the exception are preconditioned solvers, see below). That is, the solver algorithm has to work with a generic operator, whose actual form depends on the transformation used. After convergence, eigenvalues are transformed back appropriately.

For technical details of the transformations described in this chapter, the interested user is referred to [Ericsson and Ruhe, 1980], [Scott, 1982], [Nour-Omid et al., 1987], and [Meerbergen et al., 1994].

Preconditioners

As explained in the previous chapter, EPS contains preconditioned eigensolvers such as GD or JD. These solvers either apply a preconditioner at a certain step of the computation, or need to solve a correction equation with a preconditioned linear solver. One of the main goals of these solvers is to achieve a similar effect as an inverse-based spectral transformation such as shift-and-invert, but with less computational cost. For this reason, a “preconditioner” spectral transformation has been included in the ST object. However, this is just a convenient way of organizing the functionality, since this fake spectral transform cannot be used with non-preconditioned eigensolvers, and conversely preconditioned eigensolvers cannot be used with conventional spectral transformations.

Basic Usage

The ST module is the analog of some PETSc modules such as PC. The user does not usually need to create a stand-alone ST object explicitly. Instead, every EPS object internally sets up an associated ST. Therefore, the usual object management methods such as STCreate, STDestroy, STView, STSetFromOptions, are not usually called by the user.

Although the ST context is hidden inside the EPS object, the user still has control over all the options, by means of the command line, or also inside the program. To allow application programmers to set any of the spectral transformation options directly within the code, the following routine is provided to extract the ST context from the EPS object, EPSGetST

EPSGetST(EPS eps,ST *st);

After this, one is able to set any options associated with the ST object. For example, to set the value of the shift, the following function is available STSetShift

STSetShift(ST st,PetscScalar shift);

This can also be done with the command line option -st_shift <shift>. Note that the argument shift is defined as a PetscScalar, and this means that complex shifts are not allowed unless the complex version of SLEPc is used.

Other object operations are available, which are not usually called by the user. The most important of such functions are STApply, which applies the operator to a vector, and STSetUp, which prepares all the necessary data structures before the solution process starts. The term “operator” refers to one of \(A\), \(B^{-1}\!A\), \(A-\sigma I\), … depending on which kind of spectral transformation is being used.

Available Transformations

This section describes the spectral transformations that are provided in SLEPc. As in the case of eigensolvers, the spectral transformation to be used can be specified procedurally or via the command line. The application programmer can set it by means of the command STSetType

STSetType(ST st,STType type);

where type can be one of STSHIFT, STSINVERT, STCAYLEY, STPRECOND, STFILTER, or STSHELL. The ST type can also be set with the command-line option -st_type followed by the name of the method (see table Spectral transformations available in the ST package.). The first five spectral transformations are described in detail in the rest of this section. The last possibility, STSHELL, uses a specific, application-provided spectral transformation. Section Extending SLEPc describes how to implement one of these transformations.

Spectral transformations available in the ST package.

Spectral Transformation	STType	Options Name	Operator
Shift of Origin	STSHIFT	shift	\(B^{-1}A-\sigma I\)
Shift-and-invert	STSINVERT	sinvert	\((A-\sigma B)^{-1}B\)
Generalized Cayley	STCAYLEY	cayley	\((A-\sigma B)^{-1}(A+\nu B)\)
Preconditioner	STPRECOND	precond	\(K^{-1}\approx(A-\sigma B)^{-1}\)
————————-	————-	————-	——————————–
Polynomial Filter	STFILTER	filter	\(p(A)\)
Shell Transformation	STSHELL	shell	user-defined

The last column of table Spectral transformations available in the ST package. shows a general form of the operator used in each case. This generic operator can adopt different particular forms depending on whether the eigenproblem is standard or generalized, or whether the value of the shift (\(\sigma\)) and anti-shift (\(\nu\)) is zero or not. All the possible combinations are listed in table Operators used in each spectral transformation mode..

Operators used in each spectral transformation mode.

ST	Choice of \(\sigma,\nu\)	Standard problem	Generalized problem
shift	\(\sigma=0\)	\(A\)	\(B^{-1}A\)
	\(\sigma\not=0\)	\(A-\sigma I\)	\(B^{-1}A-\sigma I\)
sinvert	\(\sigma=0\)	\(A^{-1}\)	\(A^{-1}B\)
	\(\sigma\not=0\)	\((A-\sigma I)^{-1}\)	\((A-\sigma B)^{-1}B\)
cayley	\(\sigma\not=0,\nu=0\)	\((A-\sigma I)^{-1}A\)	\((A-\sigma B)^{-1}A\)
	\(\sigma=0,\nu\not=0\)	\(I+\nu A^{-1}\)	\(I+\nu A^{-1}B\)
	\(\sigma\not=0,\nu\not=0\)	\((A-\sigma I)^{-1}(A+\nu I)\)	\((A-\sigma B)^{-1}(A+\nu B)\)
precond	\(\sigma=0\)	\(K^{-1}\approx A^{-1}\)	\(K^{-1}\approx A^{-1}\)
	\(\sigma\not=0\)	\(K^{-1}\approx(A-\sigma I)^{-1}\)	\(K^{-1}\approx(A-\sigma B)^{-1}\)

The expressions shown in table Operators used in each spectral transformation mode. are not built explicitly. Instead, the appropriate operations are carried out when applying the operator to a certain vector. The inverses imply the solution of a linear system of equations that is managed by setting up an associated KSP object. The user can control the behavior of this object by adjusting the appropriate options, as will be illustrated with examples in section Solution of Linear Systems.

Relation between Target and Shift

In all transformations except STSHIFT, there is a direct connection between the target \(\tau\) (described in section Defining the Problem) and the shift \(\sigma\), as will be discussed below. The normal usage is that the user sets the target and then \(\sigma\) is set to \(\tau\) automatically (though it is still possible for the user to set a different value of the shift).

Shift of Origin

By default, no spectral transformation is performed. This is equivalent to a shift of origin (STSHIFT) with \(\sigma=0\), that is, the first line of table Operators used in each spectral transformation mode.. The solver works with the original expressions of the eigenvalue problems,

(1)\[Ax=\lambda x,\]

for standard problems, and \(Ax=\lambda Bx\) for generalized ones. Note that this last equation is actually treated internally as

(2)\[B^{-1}Ax=\lambda x.\]

When the eigensolver in EPS requests the application of the operator to a vector, a matrix-vector multiplication by matrix \(A\) is carried out (in the standard case) or a matrix-vector multiplication by matrix \(A\) followed by a linear system solve with coefficient matrix \(B\) (in the generalized case). Note that in the last case, the operation will fail if matrix \(B\) is singular.

When the shift, \(\sigma\), is given a value different from the default, 0, the effect is to move the whole spectrum by that exact quantity, \(\sigma\), which is called shift of origin. To achieve this, the solver works with the shifted matrix, that is, the expressions it has to cope with are

(3)\[(A-\sigma I)x=\theta x,\]

for standard problems, and

(4)\[(B^{-1}A-\sigma I) x=\theta x,\]

for generalized ones. The important property that is used is that shifting does not alter the eigenvectors and that it does change the eigenvalues in a simple known way, it shifts them by \(\sigma\). In both the standard and the generalized problems, the following relation holds

(5)\[\theta=\lambda-\sigma.\]

This means that after the solution process, the value \(\sigma\) has to be added[1] to the computed eigenvalues, \(\theta\), in order to retrieve the solution of the original problem, \(\lambda\). This is done by means of the function STBackTransform, which does not need to be called directly by the user.

Shift-and-invert

The shift-and-invert spectral transformation.

The shift-and-invert spectral transformation (STSINVERT) is used to enhance convergence of eigenvalues in the neighborhood of a given value. In this case, the solver deals with the expressions

(6)\[(A-\sigma I)^{-1}x=\theta x,\]

(7)\[(A-\sigma B)^{-1}B x=\theta x,\]

for standard and generalized problems, respectively. This transformation is effective for finding eigenvalues near \(\sigma\) since the eigenvalues \(\theta\) of the operator that are largest in magnitude correspond to the eigenvalues \(\lambda\) of the original problem that are closest to the shift \(\sigma\) in absolute value, as illustrated in figure The shift-and-invert spectral transformation. for an example with real eigenvalues. Once the wanted eigenvalues have been found, they may be transformed back to eigenvalues of the original problem. Again, the eigenvectors remain unchanged. In this case, the relation between the eigenvalues of both problems is

(8)\[\theta=1/(\lambda-\sigma).\]

Therefore, after the solution process, the operation to be performed in function STBackTransform is \(\lambda=\sigma+1/\theta\) for each of the computed eigenvalues.

This spectral transformation is used in the spectrum slicing technique, see section Spectrum Slicing.

Cayley

The generalized Cayley transform (STCAYLEY) is defined from the expressions

(9)\[(A-\sigma I)^{-1}(A+\nu I)x=\theta x,\]

(10)\[(A-\sigma B)^{-1}(A+\nu B)x=\theta x,\]

for standard and generalized problems, respectively. Sometimes, the term Cayley transform is applied for the particular case in which \(\nu=\sigma\). This is the default if \(\nu\) is not given a value explicitly. The value of \(\nu\) (the anti-shift) can be set with the following function STCayleySetAntishift

STCayleySetAntishift(ST st,PetscScalar nu);

or in the command line with -st_cayley_antishift.

This transformation is mathematically equivalent to shift-and-invert and, therefore, it is effective for finding eigenvalues near \(\sigma\) as well. However, in some situations it is numerically advantageous with respect to shift-and-invert (see [Bai et al., 2000, 11.2], [Lehoucq and Salinger, 2001]).

In this case, the relation between the eigenvalues of both problems is

(11)\[\theta=(\lambda+\nu)/(\lambda-\sigma).\]

Therefore, after the solution process, the operation to be performed in function STBackTransform is \(\lambda=(\theta\sigma+\nu)/(\theta-1)\) for each of the computed eigenvalues.

Preconditioner

As mentioned in the introduction of this chapter, the special type STPRECOND is used for handling preconditioners or preconditioned iterative linear solvers, which are used in the context of preconditioned eigensolvers for expanding the subspace. For instance, in the GD solver the so-called correction vector \(d_i\) to be added to the subspace in each iteration is computed as

(12)\[d_i=K^{-1}P_i(A-\theta_i B)x_i,\]

where \((\theta_i,x_i)\) is the current approximation of the sought-after eigenpair, and \(P_i\) is a projector involving \(x_i\) and \(K^{-1}x_i\). In the above expressions, \(K\) is a preconditioner matrix that is built from \(A-\theta_i B\). However, since \(\theta_i\) changes at each iteration, which would force recomputation of the preconditioner, we opt for using

(13)\[K^{-1}\approx (A-\sigma B)^{-1}.\]

Similarly, in the JD eigensolver the expansion of the subspace is carried out by solving a correction equation similar to

(14)\[(I-x_ix_i^*)(A-\theta_i B)(I-x_ix_i^*)d_i=-(A-\theta_i B)x_i,\]

where the system is solved approximately with a preconditioned iterative linear solver. For building the preconditioner of this linear system, the projectors \(I-x_ix_i^*\) are ignored, and again it is not recomputed in each iteration. Therefore, the preconditioner is built as in equation (13) as well.

It should be clear from the previous discussion, that STPRECOND does not work in the same way as the rest of spectral transformations. In particular, it does not rely on STBackTransform. It is rather a convenient mechanism for handling the preconditioner and linear solver (see examples in section Solution of Linear Systems). The expressions shown in tables Spectral transformations available in the ST package. and Operators used in each spectral transformation mode. are just a reference to indicate from which matrix the preconditioner is built by default.

There is the possibility that the user overrides the default behaviour, that is, to explicitly supply a matrix from which the preconditioner is to be built, with STSetPreconditionerMat

STSetPreconditionerMat(ST st,Mat mat);

The above function can also be used in other spectral transformations such as shift-and-invert in case the user has a cheap approximation \(K\) of the coefficient matrix \(A-\sigma B\). An alternative is to pass approximations of both \(A\) and \(B\) so that ST builds the preconditioner matrix internally, with STSetSplitPreconditioner

STSetSplitPreconditioner(ST st,PetscInt n,Mat Psplit[],MatStructure strp);

Note that preconditioned eigensolvers in EPS select STPRECOND by default, so the user does not need to specify it explicitly.

Polynomial Filtering

The type STFILTER is also special. It is used in the case of standard symmetric (or Hermitian) eigenvalue problems when the eigenvalues of interest are interior to the spectrum and we want to avoid the high cost associated with the matrix factorization of the shift-and-invert spectral transformation. The techniques generically known as polynomial filtering aim at this goal.

The polynomial filtering methods address the eigenvalue problem

(15)\[p(A)x=\theta x,\]

where \(p(\cdot)\) is a suitable high-degree polynomial. Once the polynomial is built, the eigensolver relies on STApply to compute approximations of the eigenvalues \(\theta\) of the transformed problem. These approximations must be processed in some way in order to recover the \(\lambda\) eigenvalues. Note that in this case there is no STBackTransform operation. Details of the method can be found in [Fang and Saad, 2012].

Polynomial filtering techniques are still under development in SLEPc, and will be improved in future versions. Note that the external package EVSL also implements polynomial filters to compute all eigenvalues in an interval.

Advanced Usage

Using the ST object is very straightforward. However, when using spectral transformations many things are happening behind the scenes, mainly the solution of linear systems of equations. The user must be aware of what is going on in each case, so that it is possible to guide the solution process in the most beneficial way. This section describes several advanced aspects that can have a considerable impact on efficiency.

Solution of Linear Systems

In many of the cases shown in table Operators used in each spectral transformation mode., the operator contains an inverted matrix, which means that a linear system of equations must be solved whenever the application of the operator to a vector is required. These cases are handled internally by means of a KSP object.

In the simplest case, a generalized problem is to be solved with a zero shift. Suppose you run a program that solves a generalized eigenproblem, with default options:

$ ./program

In this case, the ST object associated with the EPS solver creates a KSP object whose coefficient matrix is \(B\). By default, this KSP object is set to use a direct solver, in particular an LU factorization. However, default settings can be changed, as illustrated below.

The following command-line is equivalent to the previous one:

$ ./program -st_ksp_type preonly -st_pc_type lu

The two options specify the type of the linear solver and preconditioner to be used. The -st_ prefix indicates that the option corresponds to the linear solver within ST. The combination preonly\(+\)lu instructs to use a direct solver (LU factorization, see PETSc’s documentation for details), so this is the same as the default. Adding a new option changes the default behaviour, for instance

$ ./program -st_ksp_type preonly -st_pc_type lu
            -st_pc_factor_mat_solver_type mumps

In this case, an external linear solver package is used (MUMPS, see PETSc’s documentation for other available packages). Note that an external package is required for computing a matrix factorization in parallel, since PETSc itself only provides sequential direct linear solvers.

Instead of a direct linear solver, it is possible to use an iterative solver. This may be necessary in some cases, specially for very large problems. However, the user is warned that using an iterative linear solver makes the overall solution process less robust (see also the discussion of preconditioned eigensolvers below). As an example, the command-line

$ ./program -st_ksp_type gmres -st_pc_type bjacobi -st_ksp_rtol 1e-9

selects the GMRES solver with block Jacobi preconditioning. In the case of iterative solvers, it is important to use an appropriate tolerance, usually slightly more stringent for the linear solves relative to the desired accuracy of the eigenvalue calculation (\(10^{-9}\) in the example, compared to \(10^{-8}\) for the eigensolver).

Although the direct solver approach may seem too costly, note that the factorization is only carried out at the beginning of the eigenvalue calculation and this cost is amortized in each subsequent application of the operator. This is also the case for iterative methods with preconditioners with high-cost set-up such as ILU.

The application programmer is able to set the desired linear systems solver options also from within the code. In order to do this, first the context of the KSP object must be retrieved with the following function STGetKSP

STGetKSP(ST st,KSP *ksp);

The above functionality is also applicable to the other spectral transformations. For instance, for the shift-and-invert technique with \(\tau=10\) using BiCGStab+Jacobi:

$ ./program -st_type sinvert -eps_target 10 -st_ksp_type bcgs -st_pc_type jacobi

In shift-and-invert and Cayley, unless \(\sigma=0\), the coefficient matrix is not a simple matrix but an expression that can be explicitly constructed or not, depending on the user’s choice. This issue is examined in detail in section Explicit Computation of Coefficient Matrix below.

In many cases, especially if a shift-and-invert or Cayley transformation is being used, iterative methods may not be well suited for solving linear systems (because of the properties of the coefficient matrix that can be indefinite and ill-conditioned). When using an iterative linear solver, it may be helpful to run with the option -st_ksp_converged_reason, which will display the number of iterations required in each operator application. In extreme cases, the iterative solver fails, so EPSSolve aborts with an error

[0]PETSC ERROR: KSP did not converge (reason=DIVERGED_ITS)!

If this happens, it is necessary to use a direct method for solving the linear systems, as explained above.

The Case of Preconditioned Eigensolvers

The KSP object contained internally in ST is also used for applying the preconditioner or solving the correction equation in preconditioned eigensolvers.

The GD eigensolver employs just a preconditioner. Therefore, by default it sets the KSP type to preonly (no other KSP is allowed) and the PC type to jacobi. The user may change the preconditioner, for example as

$ ./ex5 -eps_type gd -st_pc_type asm

The JD eigensolver uses both an iterative linear solver and a preconditioner, so both KSP and PC are meaningful in this case. It is important to note that, contrary to the ordinary spectral transformations where a direct linear solver is recommended, in JD using an iterative linear solver is usually better than a direct solver. Indeed, the best performance may be achieved with a few iterations of the linear solver (or a large tolerance). For instance, the next example uses JD with GMRES+Jacobi limiting to 10 the number of allowed iterations for the linear solver:

$ ./ex5 -eps_type jd -st_ksp_type gmres -st_pc_type jacobi -st_ksp_max_it 10

A discussion on the different options available for the Davidson solvers can be found in [Romero and Roman, 2014].

Explicit Computation of Coefficient Matrix

Three possibilities can be distinguished regarding the form of the coefficient matrix of the linear systems of equations associated with the different spectral transformations. The possible coefficient matrices are:

• Simple: \(B\).

• Shifted: \(A-\sigma I\).

• Axpy: \(A-\sigma B\).

The first case has already been described and presents no difficulty. In the other two cases, there are three possible approaches:

“shell”

To work with the corresponding expression without forming the matrix explicitly. This is achieved by internally setting a matrix-free matrix with MatCreateShell.

“inplace”

To build the coefficient matrix explicitly. This is done by means of a MatShift or a MatAXPY operation, which overwrites matrix \(A\) with the corresponding expression. This alteration of matrix \(A\) is reversed after the eigensolution process has finished.

“copy”

To build the matrix explicitly, as in the previous option, but using a working copy of the matrix, that is, without modifying the original matrix \(A\).

The default behavior is to build the coefficient matrix explicitly in a copy of \(A\) (option “copy”). The user can change this as in the following example

$ ./program -st_type sinvert -eps_target 10 -st_ksp_type cg
            -st_pc_type jacobi -st_matmode shell

As always, the procedural equivalent is also available for specifying this option in the code of the program: STSetMatMode

STSetMatMode(ST st,STMatMode mode);

The user must consider which approach is the most appropriate for the particular application. The different options have advantages and drawbacks. The “shell” approach is the simplest one but severely restricts the number of possibilities available for solving the system, in particular most of the PETSc preconditioners would not be available, including direct methods. The only preconditioners that can be used in this case are Jacobi (only if matrices \(A\) and \(B\) have the operation MATOP_GET_DIAGONAL) or a user-defined one.

The second approach (“inplace”) can be much faster, specially in the generalized case. A more important advantage of this approach is that, in this case, the linear system solver can be combined with any of the preconditioners available in PETSc, including those which need to access internal matrix data-structures such as ILU. The main drawback is that, in the generalized problem, this approach probably makes sense only in the case that \(A\) and \(B\) have the same sparse pattern, because otherwise the function MatAXPY might be inefficient. If the user knows that the pattern is the same (or a subset), then this can be specified with the function STSetMatStructure

STSetMatStructure(ST st,MatStructure str);

Note that when the value of the shift \(\sigma\) is very close to an eigenvalue, then the linear system will be ill-conditioned and using iterative methods may be problematic. On the other hand, in symmetric definite problems, the coefficient matrix will be indefinite whenever \(\sigma\) is a point in the interior of the spectrum.

The third approach (“copy”) uses more memory but avoids a potential problem that could appear in the “inplace” approach: the recovered matrix might be slightly different from the original one (due to roundoff).

Preserving the Symmetry in Generalized Eigenproblems

As mentioned in section Defining the Problem, some eigensolvers can exploit symmetry and compute a solution for Hermitian problems with less storage and/or computational cost than other methods. Also, symmetric solvers can be more accurate in some cases. However, in the case of generalized eigenvalue problems in which both \(A\) and \(B\) are symmetric, it happens that, due to the spectral transformation, symmetry is lost since none of the transformed operators \(B^{-1}\!A-\sigma I\), \((A-\sigma B)^{-1}B\), etc. is symmetric (the same applies in the Hermitian case for complex matrices).

The solution proposed in SLEPc is based on selecting different kinds of inner products. Currently, we have the following choice of inner products:

• Standard Hermitian inner product: \(\langle x,y\rangle=x^*y\).

• \(B\)-inner product: \(\langle x,y\rangle_B=x^*\!B\,y\).

The second one can be used for preserving the symmetry in symmetric definite generalized problems, as described below. Note that \(\langle x,y\rangle_B\) is a genuine inner product only if \(B\) is symmetric positive definite (for the case of symmetric positive semi-definite \(B\) see section Purification of Eigenvectors).

It can be shown that \(\mathbb{R}^n\) with the \(\langle x,y\rangle_B\) inner product is isomorphic to the Euclidean \(n\)-space \(\mathbb{R}^n\) with the standard Hermitian inner product. This means that if we use \(\langle x,y\rangle_B\) instead of the standard inner product, we are just changing the way lengths and angles are measured, but otherwise all the algebraic properties are maintained and therefore algorithms remain correct. What is interesting to observe is that the transformed operators such as \(B^{-1}\!A\) or \((A-\sigma B)^{-1}B\) are self-adjoint with respect to \(\langle x,y\rangle_B\).

Abstraction used by SLEPc solvers.

Internally, SLEPc operates with the abstraction illustrated in figure Abstraction used by SLEPc solvers.. The operations indicated by dashed arrows are implemented as virtual functions. From the user point of view, all the above explanation is transparent. The only thing he/she has to care about is to set the problem type appropriately with EPSSetProblemType (see section Defining the Problem). In the case of the Cayley transform, SLEPc is using \(\langle x,y\rangle_{A+\nu B}\) as the inner product for preserving symmetry.

Using the \(B\)-inner product may be attractive also in the non-symmetric case (\(A\) non-symmetric) as described in the next subsection.

Note that the above discussion is not directly applicable to STPRECOND and the preconditioned eigensolvers, in the sense that the goal is not to recover the symmetry of the operator. Still, the \(B\)-inner product is also used in generalized symmetric-definite problems.

Purification of Eigenvectors

In generalized eigenproblems, the case of singular \(B\) deserves especial consideration. In this case the default spectral transformation (STSHIFT) cannot be used since \(B^{-1}\) does not exist.

In shift-and-invert with operator matrix \(T=(A-\sigma B)^{-1}B\), when \(B\) is singular all the eigenvectors that belong to finite eigenvalues are also eigenvectors of \(T\) and belong to the range of \(T\), \(\mathcal{R}(T)\). In this case, the bilinear form \(\langle x,y\rangle_B\) introduced in section Preserving the Symmetry in Generalized Eigenproblems is a semi-inner product and \(\|x\|_B=\sqrt{\langle x,x\rangle_B}\) is a semi-norm. As before, \(T\) is self-adjoint with respect to this inner product since \(B\,T=T^*B\). Also, \(\langle x,y\rangle_B\) is a true inner product on \(\mathcal{R}(T)\).

The implication of all this is that, for singular \(B\), if the \(B\)-inner product is used throughout the eigensolver then, assuming that the initial vector has been forced to lie in \(\mathcal{R}(T)\), the computed eigenvectors should be correct, i.e., they should belong to \(\mathcal{R}(T)\) as well. Nevertheless, finite precision arithmetic spoils this nice picture, and computed eigenvectors are easily corrupted by components of vectors in the null-space of \(B\). Additional computation is required for achieving the desired property. This is usually referred to as eigenvector purification.

Although more elaborate purification strategies have been proposed (usually trying to reduce the computational effort, see [Nour-Omid et al., 1987] and [Meerbergen and Spence, 1997]), the approach in SLEPc is simply to explicitly force the initial vector in the range of \(T\), with \(v_0\leftarrow Tv_0\), as well as the computed eigenvectors at the end, \(x_i\leftarrow Tx_i\). Since this computation can be costly, it can be deactivated if the user knows that \(B\) is non-singular, with EPSSetPurify

EPSSetPurify(EPS eps,PetscBool purify);

A final comment is that eigenvector corruption happens also in the non-symmetric case. If \(A\) is non-symmetric but \(B\) is symmetric positive semi-definite, then the scheme presented above (\(B\)-inner product together with purification) can still be applied and is generally more successful than the straightforward approach with the standard inner product. For using this scheme in SLEPc, the user has to specify the special problem type EPS_PGNHEP, see table Problem types considered in EPS..

Spectrum Slicing

In the context of symmetric-definite generalized eigenvalue problems (EPS_GHEP) it is often required to compute all eigenvalues contained in a given interval \([a,b]\). This poses some difficulties, such as:

• The number of eigenvalues in the interval is not known a priori.

• There might be many eigenvalues, in some applications a significant percentage of the spectrum (20%, say).

• We must make certain that no eigenvalues are missed, and in particular all eigenvalues must be computed with their correct multiplicity.

• In some applications, the interval is open in one end, i.e., either \(a\) or \(b\) can be infinite.

One possible strategy to solve this problem is to sweep the interval from one end to the other, computing chunks of eigenvalues with a spectral transformation that updates the shift dynamically. This is generally referred to as spectrum slicing. The method implemented in SLEPc is similar to that proposed by Grimes et al. [1994], where inertia information is used to validate sub-intervals. Given a symmetric-indefinite triangular factorization

(16)\[A-\sigma B=LDL^T,\]

by Sylvester’s law of inertia we know that the number of eigenvalues on the left of \(\sigma\) is equal to the number of negative eigenvalues of \(D\),

(17)\[\nu(A-\sigma B)=\nu(D).\]

A detailed description of the method available in SLEPc can be found in [Campos and Roman, 2012]. The SLEPc interface hides all the complications of the algorithm. However, the user must be aware of all the restrictions for this technique to be employed:

• This is currently implemented only in Krylov-Schur.

• The method is based on shift-and-invert, so STSINVERT must be used. Furthermore, direct linear solvers are required.

• The direct linear solver must provide the matrix inertia (see PETSc’s MatGetInertia).

An example command-line that sets up all the required options is:

$ ./ex2 -n 50 -eps_interval 0.4,0.8 -st_type sinvert
              -st_ksp_type preonly -st_pc_type cholesky

Note that PETSc’s Cholesky factorization is not parallel, so for doing spectrum slicing in parallel it is required to use an external solver that supports inertia. For example, with MUMPS (see section Solution of Linear Systems on how to use external linear solvers) we would do:

$ ./ex2 -n 50 -eps_interval 0.4,0.8 -st_type sinvert
              -st_ksp_type preonly -st_pc_type cholesky
              -st_pc_factor_mat_solver_type mumps -st_mat_mumps_icntl_13 1

The last option is required by MUMPS to compute the inertia. An alternative is to use SuperLU_DIST, in which case the required options would be:

$ ./ex2 -n 50 -eps_interval 0.4,0.8 -st_type sinvert
              -st_ksp_type preonly -st_pc_type cholesky
              -st_pc_factor_mat_solver_type superlu_dist
              -st_mat_superlu_dist_rowperm NOROWPERM

In the latter example, MatSetOption must be used in both matrices to explicitly state that they are symmetric (or Hermitian in the complex case).

Apart from the above recommendations, the following must be taken into account:

• The parameters nev and ncv from EPSSetDimensions are determined internally (user-provided values are ignored, and set to the number of eigenvalues in the interval). It is possible for the user to specify the “local” nev and ncv parameters that will be used when computing eigenvalues around each shift, with EPSKrylovSchurSetDimensions.

• The user can also tune the computation by setting the value of max_it.

Usage with Complex Scalars

Some external packages that provide inertia information (MUMPS, Pardiso) do so only in real scalars, but not in the case of complex scalars. Hence, with complex scalars spectrum slicing is available only sequentially (with PETSc’s Cholesky factorization) or via SuperLU_DIST (as in the last example above). An alternative to spectrum slicing is to use the CISS solver with a region enclosing an interval on the real axis, see Maeda et al. [2016] for details.

Use of Multiple Communicators

Since spectrum slicing requires direct linear solves, parallel runs may suffer from bad scalability in the sense that increasing the number of MPI processes does not imply a performance gain. For this reason, SLEPc provides the option of using multiple communicators, that is, splitting the initial MPI communicator in several groups, each of them in charge of processing part of the interval.

The multi-communicator setting is activated with a value of npart>1 in EPSKrylovSchurSetPartitions

EPSKrylovSchurSetPartitions(EPS eps,PetscInt npart);

The interval \([a,b]\) is then divided in npart subintervals of equal size, and the problem of computing all eigenvalues in \([a,b]\) is divided in npart independent subproblems. Each subproblem is solved using only a subset of the initial \(p\) processes, with \(\lceil p/\texttt{npart}\rceil\) processes at most. A final step will gather all computed solutions so that they are available in the whole EPS communicator.

The division of the interval in subintervals is done blindly, and this may result in load imbalance if some subintervals contain much more eigenvalues than others. This can be prevented by passing a list of subinterval boundaries, provided that the user has a priori information to roughly determine the eigenvalue distribution: EPSKrylovSchurSetSubintervals

EPSKrylovSchurSetSubintervals(EPS eps,PetscReal *subint);

An additional benefit of multi-communicator support is that it enables parallel spectrum slicing runs without the need to install a parallel direct solver (MUMPS), by setting the number of partitions equal to the number of MPI processes. The following command-line example uses sequential linear solves in 4 partitions, one process each:

$ mpiexec -n 4 ./ex25 -eps_interval 0.4,0.8 -eps_krylovschur_partitions 4
          -st_type sinvert -st_ksp_type preonly -st_pc_type cholesky

The analog example using MUMPS with 5 processes in each partition:

$ mpiexec -n 20 ./ex25 -eps_interval 0.4,0.8 -eps_krylovschur_partitions 4
          -st_type sinvert -st_ksp_type preonly -st_pc_type cholesky
          -st_pc_factor_mat_solver_type mumps -st_mat_mumps_icntl_13 1

Spectrum Folding

In SLEPc versions prior to 3.5, ST had another type intended to perform the spectrum folding technique described below. It is no longer available with ST, but it can be implemented directly in application code as illustrated in example ex24.c.

Spectrum folding involves squaring in addition to shifting. This makes sense for standard Hermitian eigenvalue problems, where the transformed problem to be addressed is

(18)\[(A-\sigma I)^2x=\theta x.\]

The following relation holds

(19)\[\theta=(\lambda-\sigma)^2.\]

Note that the mapping between \(\lambda\) and \(\theta\) is not injective, and hence this cannot be considered a true spectral transformation.

The effect is that the spectrum is folded around the value of \(\sigma\). Thus, eigenvalues that are closest to the shift become the smallest eigenvalues in the folded spectrum, as illustrated in figure Illustration of the effect of spectrum folding.. For this reason, spectrum folding is commonly used in combination with eigensolvers that compute the smallest eigenvalues, for instance in the context of electronic structure calculations [Canning et al., 2000]. This transformation can be an effective, low-cost alternative to shift-and-invert.

Illustration of the effect of spectrum folding.

Footnotes

[Bai:2000:TSA]

Z. Bai, J. Demmel, J. Dongarra, A. Ruhe, and H. van der Vorst, editors. Templates for the Solution of Algebraic Eigenvalue Problems: A Practical Guide. Society for Industrial and Applied Mathematics, Philadelphia, PA, 2000.

[Campos:2012:SSS]

C. Campos and J. E. Roman. Strategies for spectrum slicing based on restarted Lanczos methods. Numerical Algorithms, 60(2):279–295, 2012.

[Canning:2000:PEP]

A. Canning, L. W. Wang, A. Williamson, and A. Zunger. Parallel empirical pseudopotential electronic structure calculations for million atom systems. Journal of Computational Physics, 160(1):29–41, 2000.

[Ericsson:1980:STL]

T. Ericsson and A. Ruhe. The spectral transformation Lanczos method for the numerical solution of large sparse generalized symmetric eigenvalue problems. Mathematics of Computation, 35(152):1251–1268, 1980.

[Fang:2012:FLP]

H. Fang and Y. Saad. A filtered Lanczos procedure for extreme and interior eigenvalue problems. SIAM Journal on Scientific Computing, 34(4):A2220–A2246, 2012.

[Grimes:1994:SBL]

Roger G. Grimes, John G. Lewis, and Horst D. Simon. A shifted block Lanczos algorithm for solving sparse symmetric generalized eigenproblems. SIAM Journal on Matrix Analysis and Applications, 15(1):228–272, 1994.

[Lehoucq:2001:LEC]

R. B. Lehoucq and A. G. Salinger. Large-scale eigenvalue calculations for stability analysis of steady flows on massively parallel computers. International Journal for Numerical Methods in Fluids, 36:309–327, 2001.

[str-11]

Y. Maeda, T. Sakurai, and J. E. Roman. Contour Integral Spectrum Slicing Method in SLEPc. Technical Report STR-11, Universitat Politècnica de València, 2016. Available at https://slepc.upv.es.

[Meerbergen:1994:SCT]

K. Meerbergen, A. Spence, and D. Roose. Shift-invert and Cayley transforms for detection of rightmost eigenvalues of nonsymmetric matrices. BIT Numerical Mathematics, 34(3):409–423, 1994.

[Meerbergen:1997:IRA]

Karl Meerbergen and Alastair Spence. Implicitly restarted Arnoldi with purification for the shift-invert transformation. Mathematics of Computation, 66(218):667–689, 1997.

[Nour-Omid:1987:HIS] (1,2)

Bahram Nour-Omid, Beresford N. Parlett, Thomas Ericsson, and Paul S. Jensen. How to implement the spectral transformation. Mathematics of Computation, 48(178):663–673, 1987.

[Romero:2014:PID]

E. Romero and J. E. Roman. A parallel implementation of Davidson methods for large-scale eigenvalue problems in SLEPc. ACM Transactions on Mathematical Software, 40(2):13:1–13:29, 2014. doi:10.1145/2543696.

[Scott:1982:AIO]

D. S. Scott. The advantages of inverted operators in Rayleigh-Ritz approximations. SIAM Journal on Scientific and Statistical Computing, 3(1):68–75, 1982.

---

[1]

Note that the sign changed in SLEPc 3.5 with respect to previous versions.