# Subversion Repositoriesslepc-dev

Rev
Rev Author Line No. Line
1843 jroman 1
%-------------------------------------------------------
2
% SLEPc Users Manual
3
%-------------------------------------------------------
4
\chapter{\label{cap:st}ST: Spectral Transformation}
5
%-------------------------------------------------------
6

7
\noindent The Spectral Transformation (\ident{ST}) is the \slepc object that encapsulates the functionality required for acceleration techniques based on the transformation of the spectrum. All the eigensolvers provided in \ident{EPS} work by applying an operator to a set of vectors and this operator can adopt different forms. The \ident{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 type of spectral transformation can be specified at run time, as well as several parameters such as the value of the shift.
2088 jroman 8

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

1843 jroman 11
%---------------------------------------------------
12
\section{General Description}
13

14
        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:
15
\begin{itemize}
16
\item Compute internal eigenvalues. In some applications, the eigenpairs of interest are not the extreme ones (largest magnitude, smallest magnitude, rightmost, leftmost), but those contained in a certain interval or those closest to a certain value of the complex plane.
17
\item 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.
18
\item Handle some special situations. For instance, in generalized problems when matrix $B$ is singular, it may be necessary to use a spectral transformation.
19
\end{itemize}
20

21
        \slepc separates spectral transformations from solution methods so that any combination of them can be specified by the user. To achieve this, all the eigensolvers contained in \ident{EPS} must be implemented in such a way that they are independent of which transformation has been selected by the user. 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.
22

23
        For technical details of the transformations described in this chapter, the interested user is referred to \citep{Ericsson:1980:STL}, \citep{Scott:1982:AIO}, \citep{Nour-Omid:1987:HIS}, and \citep{Meerbergen:1994:SCT}.
24

2088 jroman 25
\paragraph{Preconditioners.}
26

2093 jroman 27
As explained in the previous chapter, \ident{EPS} contains preconditioned eigensolvers. 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 \ident{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.
2088 jroman 28

1843 jroman 29
%---------------------------------------------------
30
\section{Basic Usage}
31

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

34
        Although the \ident{ST} context is hidden inside the \ident{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 \ident{ST} context from the \ident{EPS} object,
35
        \findex{EPSGetST}
36
        \begin{Verbatim}[fontsize=\small]
37
        EPSGetST(EPS eps,ST *st);
38
        \end{Verbatim}
39

40
        After this, one is able to set any options associated to the \ident{ST} object. For example, to set the value of the shift, the following function is available
41
        \findex{STSetShift}
42
        \begin{Verbatim}[fontsize=\small]
43
        STSetShift(ST st,PetscScalar shift);
44
        \end{Verbatim}
45
        This can also be done with the command line option \Verb!-st_shift <shift>!. Note that the argument \texttt{shift} is defined as a \texttt{PetscScalar}, and this means that complex shifts are not allowed unless the complex version of \slepc is used.
46

47
        Other object operations are available, which are not usually called by the user. The most important of such functions are \ident{STApply}, which applies the operator to a vector, and \ident{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.
48

49
%---------------------------------------------------
50
\section{Available Transformations}
51

52
        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
53
        \findex{STSetType}
54
        \begin{Verbatim}[fontsize=\small]
55
        STSetType(ST st,STType type);
56
        \end{Verbatim}
57
where \texttt{type} can be one of
58
\texttt{STSHIFT}, \texttt{STFOLD},
2092 jroman 59
\texttt{STSINVERT}, \texttt{STCAYLEY},
2093 jroman 60
\texttt{STPRECOND}, or \texttt{STSHELL}.
61
The \ident{ST} type can also be set with the command-line option \Verb!-st_type! followed by the name of the method (see Table \ref{tab:transforms}). The first five spectral transformations are described in detail in the rest of this section. The last possibility, \texttt{STSHELL}, uses a specific, application-provided spectral transformation. Section \ref{sec:shell} describes how to implement one of these transformations.
1843 jroman 62

63
\begin{table}
64
\centering
65
{\small \begin{tabular}{lllc}
66
                        &                   & {\footnotesize Options} &\\
67
Spectral Transformation & \ident{STType}    & {\footnotesize Name}    & Operator\\\hline
68
Shift of Origin         & \texttt{STSHIFT}  & \texttt{shift}   & $B^{-1}A+\sigma I$\\
69
Spectrum Folding        & \texttt{STFOLD}   & \texttt{fold}    & $(B^{-1}A-\sigma I)^2$\\
2092 jroman 70
Shift-and-invert        & \texttt{STSINVERT}& \texttt{sinvert} & $(A-\sigma B)^{-1}B$\\
71
Generalized Cayley      & \texttt{STCAYLEY} & \texttt{cayley}  & $(A-\sigma B)^{-1}(A+\nu B)$\\
2103 jroman 72
Preconditioner          & \texttt{STPRECOND}& \texttt{precond} & $K^{-1}\approx(A-\sigma B)^{-1}$\\\hline
1843 jroman 73
Shell Transformation    & \texttt{STSHELL}  & \texttt{shell}   & \emph{user-defined}\\\hline
74
\end{tabular} }
75
\caption{\label{tab:transforms}Spectral transformations available in the  \ident{ST} package.}
76
\end{table}
77

2091 jroman 78
        The last column of Table \ref{tab:transforms} 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 illustrated in Table \ref{tab:op}.
1843 jroman 79
        \begin{table}
80
        \centering
81
        {\small \begin{tabular}{llcc}
2091 jroman 82
        \ident{ST}     & Choice of $\sigma,\nu$ & Standard problem & Generalized problem \\\hline
1843 jroman 83
        \texttt{shift}
84
        & $\sigma=0$     & $A$           & $B^{-1}A$          \\
85
        & $\sigma\not=0$ & $A+\sigma I$  & $B^{-1}A+\sigma I$ \\ \hline
86
        \texttt{fold}
87
        & $\sigma=0$     & $A^2$         & $(B^{-1}A)^2$          \\
88
        & $\sigma\not=0$ & $(A-\sigma I)^2$ & $(B^{-1}A-\sigma I)^2$ \\ \hline
89
        \texttt{sinvert}
90
        & $\sigma=0$     & $A^{-1}$      & $A^{-1}B$          \\
91
        & $\sigma\not=0$ & $(A-\sigma I)^{-1}$  & $(A-\sigma B)^{-1}B$ \\ \hline
92
        \texttt{cayley}
2091 jroman 93
        & $\sigma\not=0,\nu=0$ & $(A-\sigma I)^{-1}A$  & $(A-\sigma B)^{-1}A$ \\
94
        & $\sigma=0,\nu\not=0$     & $I+\nu A^{-1}$      & $I+\nu A^{-1}B$ \\
95
        & $\sigma\not=0,\nu\not=0$ & $(A-\sigma I)^{-1}(A+\nu I)$  & $(A-\sigma B)^{-1}(A+\nu B)$ \\ \hline
2093 jroman 96
        \texttt{precond}
2103 jroman 97
        & $\sigma=0$     & $K^{-1}\approx A^{-1}$ & $K^{-1}\approx A^{-1}$          \\
98
        & $\sigma\not=0$ & $K^{-1}\approx(A-\sigma I)^{-1}$ & $K^{-1}\approx(A-\sigma B)^{-1}$ \\ \hline
1843 jroman 99
        \end{tabular} }
100
        \caption{\label{tab:op}Operators used in each spectral transformation mode.}
101
        \end{table}
102

103
        The expressions shown in Table \ref{tab:op} 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 \ident{KSP} object. The user can control the behavior of this object by adjusting the appropriate options, as will be illustrated with examples in \S\ref{sec:lin}.
104

2093 jroman 105
\paragraph{Relation between Target and Shift.}
106

107
        In all transformations except \ident{STSHIFT}, there is a direct connection between the target $\tau$ (described in \S\ref{sec:defprob}) 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).
108

1843 jroman 109
\subsection{Shift of Origin}
110

111
        By default, no spectral transformation is performed. This is equivalent to a shift of origin (\texttt{STSHIFT}) with $\sigma=0$, that is, the first line of Table \ref{tab:op}. The solver works with the original expressions of the eigenvalue problems,
112
\begin{equation}Ax=\lambda x\;\;,\end{equation}
113
for standard problems, and $Ax=\lambda Bx$ for generalized ones. Note that this last equation is actually treated internally as
114
\begin{equation}B^{-1}Ax=\lambda x\;\;.\end{equation}
115
When the eigensolver in \ident{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.
116

117
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 \emph{shift of origin}. To achieve this, the solver works with the shifted matrix, that is, the expressions it has to cope with are
118
\begin{equation}(A+\sigma I)x=\theta x\;\;,\end{equation}
119
for standard problems, and
120
\begin{equation}(B^{-1}A+\sigma I) x=\theta x\;\;,\end{equation}
121
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
122
\begin{equation}\theta=\lambda+\sigma\;\;.\end{equation}
123
This means that after the solution process, the value $\sigma$ has to be subtracted from the computed eigenvalues, $\theta$, in order to retrieve the solution of the original problem, $\lambda$. This is done by means of the function \ident{STBackTransform}, which does not need to be called directly by the user.
124

125
\subsection{Spectrum Folding}
126

127
Spectrum folding refers to a spectral transformation that involves squaring in addition to shifting. The transformed problems to be addressed are the following
128
\begin{equation}(A-\sigma I)^2x=\theta x\;\;,\end{equation}
129
for standard problems, and
130
\begin{equation}(B^{-1}A-\sigma I)^2 x=\theta x\;\;,\end{equation}
131
for generalized ones. In both cases, the following relation holds
132
\begin{equation}\theta=(\lambda-\sigma)^2\;\;.\end{equation}
2093 jroman 133
For real eigenvalues, the effect of this transformation 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 \ref{fig:fold}. 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, \citep{Canning:2000:PEP}. This transformation can be an effective, low-cost alternative to shift-and-invert (explained below).
1843 jroman 134

135
  \begin{figure}
136
    \centering
137
    % Requires gnuplot
138
    \begin{tikzpicture}[domain=-2:2]
139
      \draw[->] (-2.5,0) -- (2.5,0) node[right] {$\lambda$};
140
      \draw[->] (0,0) node[below] {$\sigma$} -- (0,4.2) node[above] {$\theta$};
141
      \draw[color=blue] plot[id=fold] function{x*x}
142
          node[right] {$\theta = (\lambda-\sigma)^2$};
143
      \draw[very thin,color=gray] (-1.3,0) |- (0,1.69);
144
      \draw[very thin,color=gray] (-0.5,0) |- (0,0.25);
145
      \draw[very thin,color=gray] (1.1,0) |- (0,1.21);
146
      \foreach \pos/\label in {1.69/$\theta_1$,0.25/$\theta_2$,1.21/$\theta_3$}
147
      \draw (-0.1,\pos) -- (0.1,\pos) (1mm,\pos cm) node[anchor=west,font=\footnotesize] {\label};
148
      \foreach \pos/\label in {-1.3/$\lambda_1$,-0.5/$\lambda_2$,1.1/$\lambda_3$}
149
      \draw (\pos,0.1) -- (\pos,-0.1) (\pos cm,-2.5ex) node[anchor=base,font=\footnotesize] {\label};
150
    \end{tikzpicture}
151
    \caption{\label{fig:fold}Illustration of the effect of spectrum folding.}
152
  \end{figure}
153

154
\subsection{Shift-and-invert}
155

2092 jroman 156
        The shift-and-invert spectral transformation (\texttt{STSINVERT}) is used to enhance convergence of eigenvalues in the neighborhood of a given value. In this case, the solver deals with the expressions
1843 jroman 157
\begin{eqnarray}
158
(A-\sigma I)^{-1}x=\theta x\;\;,\\
159
(A-\sigma B)^{-1}B x=\theta x\;\;,
160
\end{eqnarray}
161
for standard and generalized problems, respectively.
162
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 \ref{fig:sinvert} 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.
163
In this case, the relation between the eigenvalues of both problems is
164
\begin{equation}\theta=1/(\lambda-\sigma)\;\;.\end{equation}
165
Therefore, after the solution process, the operation to be performed in function \ident{STBackTransform} is $\lambda=\sigma+1/\theta$ for each of the computed eigenvalues.
166

167
  \begin{figure}
168
    \centering
169
    % Requires gnuplot
170
    \begin{tikzpicture}
171
      \draw[->] (-1.5,0) -- (5.5,0) node[right] {$\lambda$};
172
      \draw[->] (0,-4) -- (0,4) node[above] {$\theta$};
173
      \draw[dashed,color=gray] (2,-4) -- (2,4);
174
      \draw[color=blue] plot[id=sinv1,domain=-1.3:1.75] function{1/(x-2)};
175
      \draw[color=blue] plot[id=sinv2,domain=2.25:5.3] function{1/(x-2)};
176
      \node[below right] at (0,0) {$0$};
177
      \draw (2,0.1) -- (2,-0.1) node[below left] at (2,0) {$\sigma$};
178
      \node at (4,3.2) {$\theta = \frac{1}{\lambda-\sigma}$};
179
      \draw[very thin,color=gray] (2.5,0) |- (0,2);
180
      \draw[very thin,color=gray] (3.9,0) |- (0,0.5263);
181
      \foreach \pos/\label in {2/$\theta_1$,0.5263/$\theta_2$}
182
      \draw (-0.1,\pos) -- (0.1,\pos) (-0.8mm,\pos cm) node[anchor=east,font=\footnotesize] {\label};
183
      \foreach \pos/\label in {2.5/$\lambda_1$,3.9/$\lambda_2$}
184
      \draw (\pos,0.1) -- (\pos,-0.1) (\pos cm,-2.5ex) node[anchor=base,font=\footnotesize] {\label};
185
    \end{tikzpicture}
186
    \caption{\label{fig:sinvert}The shift-and-invert spectral transformation.}
187
  \end{figure}
188

189
\subsection{Cayley}
190
\label{sec:cayley}
191

192
        The generalized Cayley transform (\texttt{STCAYLEY}) is defined from the expressions
193
\begin{eqnarray}
2091 jroman 194
(A-\sigma I)^{-1}(A+\nu I)x=\theta x\;\;,\\
195
(A-\sigma B)^{-1}(A+\nu B)x=\theta x\;\;,
1843 jroman 196
\end{eqnarray}
2091 jroman 197
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
1843 jroman 198
        \findex{STCayleySetAntishift}
199
        \begin{Verbatim}[fontsize=\small]
2091 jroman 200
        STCayleySetAntishift(ST st,PetscScalar nu);
1843 jroman 201
        \end{Verbatim}
2090 jroman 202
or in the command line with \Verb!-st_cayley_antishift!.
1843 jroman 203

204
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 \citep[\S 11.2]{Bai:2000:TSA}, \citep{Lehoucq:2001:LEC}).
205

206
In this case, the relation between the eigenvalues of both problems is
2091 jroman 207
\begin{equation}\theta=(\lambda+\nu)/(\lambda-\sigma)\;\;.\end{equation}
208
Therefore, after the solution process, the operation to be performed in function \ident{STBackTransform} is $\lambda=(\theta\sigma+\nu)/(\theta-1)$ for each of the computed eigenvalues.
1843 jroman 209

2093 jroman 210
\subsection{Preconditioner}
2103 jroman 211
\label{sec:precond}
2093 jroman 212

213
        As mentioned in the introduction of this chapter, the special type \ident{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
214
\begin{equation}
215
d_i=K^{-1}P_i(A-\theta_i B)x_i,
216
\end{equation}
217
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
218
\begin{equation}
219
\label{eq:precon}
2103 jroman 220
K^{-1}\approx (A-\sigma B)^{-1}.
2093 jroman 221
\end{equation}
222

2103 jroman 223
        Similarly, in the JD eigensolver the expansion of the subspace is carried out by solving a correction equation similar to
2093 jroman 224
\begin{equation}
225
(I-x_ix_i^*)(A-\theta_i B)(I-x_ix_i^*)d_i=-(A-\theta_i B)x_i,
226
\end{equation}
2103 jroman 227
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 expression \ref{eq:precon} as well.
2093 jroman 228

229
        It should be clear from the previous discussion, that \ident{STPRECOND} does not work in the same way as the rest of spectral transformations. In particular, it is not intended to be used on the basis of operator applications with \ident{STApply}, and it does not rely on \ident{STBackTransform} either. It is rather a convenient mechanism for handling the preconditioner and linear solver (see examples in \S\ref{sec:lin}). The expressions shown in Tables \ref{tab:transforms} and \ref{tab:op} are just a reference to indicate from which matrix the preconditioner is built by default.
230

231
        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
232
        \findex{STPrecondSetMatForPC}
233
        \begin{Verbatim}[fontsize=\small]
234
        STPrecondSetMatForPC(ST st,Mat mat);
235
        \end{Verbatim}
236

237
        Note that preconditioned eigensolvers in \ident{EPS} select \ident{STPRECOND} by default, so the user does not need to specify it explicitly.
238

1843 jroman 239
%---------------------------------------------------
240
\section{Advanced Usage}
241

242
Using the \ident{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.
243

244
\subsection{Solution of Linear Systems}
245
\label{sec:lin}
246

247
        In many of the cases shown in Table \ref{tab:op}, 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 \ident{KSP} object.
248

249
        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:
250
\begin{Verbatim}[fontsize=\small]
251
        $program 252 \end{Verbatim} 2095 jroman 253 In this case, the \ident{ST} object associated to the \ident{EPS} solver creates a \ident{KSP} object whose coefficient matrix is$B$. By default, this \ident{KSP} object is set to use a direct solver\footnote{This is the default since \slepc 3.0.0.}, in particular an LU factorization. However, default settings can be changed, as illustrated below. 1843 jroman 254   2095 jroman 255  The following command-line is equivalent to the previous one\footnote{More precisely, the equivalent command would use the \texttt{redundant} \texttt{PC} in order to work in parallel.}: 1843 jroman 256 \begin{Verbatim}[fontsize=\small] 257 $ program -st_ksp_type preonly -st_pc_type lu
258
\end{Verbatim}
259
The two options specify the type of the linear solver and preconditioner to be used. The \Verb!-st_! prefix signifies that the option corresponds to the linear solver within \ident{ST}. The combination \texttt{preonly}$+$\texttt{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
260
\begin{Verbatim}[fontsize=\small]
261
        $program -st_ksp_type preonly -st_pc_type lu 262  -st_pc_factor_mat_solver_package mumps 263 \end{Verbatim} 2095 jroman 264 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 a truly parallel direct linear solver. 1843 jroman 265   2095 jroman 266  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 1843 jroman 267 \begin{Verbatim}[fontsize=\small] 268 $ program -st_ksp_type gmres -st_pc_type bjacobi -st_ksp_rtol 1e-8
269
\end{Verbatim}
2095 jroman 270
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^{-8}$ in the example, compared to $10^{-7}$ for the eigensolver).
1843 jroman 271

272
        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.
273

274
        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 \ident{KSP} object must be retrieved with the following function
275
        \findex{STGetKSP}
276
        \begin{Verbatim}[fontsize=\small]
277
        STGetKSP(ST st,KSP *ksp);
278
        \end{Verbatim}
279

2095 jroman 280
        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:
1843 jroman 281
\begin{Verbatim}[fontsize=\small]
2095 jroman 282
        $program -st_type sinvert -eps_target 10 -st_ksp_type bcgs -st_pc_type jacobi 1843 jroman 283 \end{Verbatim} 284  The shift-and-invert and Cayley transformations deserve special consideration. In these cases, 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 \S\ref{sec:explicit} below. 285   286 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 badly conditioned). In such cases, during the execution of the application, the user may get the following message: 287 \begin{Verbatim}[fontsize=\small] 288  [0]PETSC ERROR: Warning: KSP did not converge (-3)! 289 \end{Verbatim} 290 If this happens, chances are that the \ident{EPS} object fails to compute the eigensolution or that the retrieved solution is wrong whatsoever. In that situation, it is necessary to use a direct method for solving the linear systems, as explained above. 291   2095 jroman 292 \paragraph{The Case of Preconditioned Eigensolvers.} 293   294  The \texttt{KSP} object contained internally in \ident{ST} is also used for applying the preconditioner or solving the correction equation in preconditioned eigensolvers. 295   296  The GD eigensolver employs just a preconditioner. Therefore, by default it sets the \texttt{KSP} type to \texttt{preonly} (no other \texttt{KSP} is allowed) and the \texttt{PC} type to \texttt{jacobi}. The user may change the preconditioner, for example as 297 \begin{Verbatim}[fontsize=\small] 298 $ ./ex5 -eps_type gd -st_pc_type jacobi
299
\end{Verbatim}
300

301
        The JD eigensolver uses both an iterative linear solver and a preconditioner, so both \texttt{KSP} and \texttt{PC} are meaningful in this case. The default is \texttt{gmres}+\texttt{bjacobi}. 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:
302
\begin{Verbatim}[fontsize=\small]
303
        $./ex5 -eps_type jd -st_ksp_type gmres -st_pc_type jacobi -st_ksp_max_it 10 304 \end{Verbatim} 305   1843 jroman 306 \subsection{Explicit Computation of Coefficient Matrix} 307 \label{sec:explicit} 308   309  Three possibilities can be distinguished regarding the form of the coefficient matrix of the linear systems of equations associated to the different spectral transformations. The possible coefficient matrices are: 310  \begin{itemize} 311  \item Simple:$B$. 312  \item Shifted:$A-\sigma I$. 313  \item Axpy:$A-\sigma B$. 314  \end{itemize} 315  The first case has already been described and presents no difficulty. In the other two cases, there are three possible approaches: 316  \begin{description} 317  \item[\Verb!shell!''] To work with the corresponding expression without forming the matrix explicitly. This is achieved by internally setting a matrix-free matrix with \ident{MatCreateShell}. 318  \item[\Verb!inplace!''] To build the coefficient matrix explicitly. This is done by means of a \ident{MatShift} or a \ident{MatAXPY} operation, which overwrites matrix$A$with the corresponding expression. This alteration of matrix$A$is reversed after the eigensolution process has finished. 319  \item[\Verb!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$. 320  \end{description} 321  The default behavior is to build the coefficient matrix explicitly in a copy of$A$(option \Verb!copy!''). The user can change this as in the following example 322 \begin{Verbatim}[fontsize=\small] 323 $ program -st_type sinvert -st_shift 10 -st_ksp_type cg
324
                  -st_pc_type jacobi -st_matmode shell
325
\end{Verbatim}
326
        As always, the procedural equivalent is also available for specifying this option in the code of the program:
327
        \findex{STSetMatMode}
328
        \begin{Verbatim}[fontsize=\small]
329
        STSetMatMode(ST st,STMatMode mode);
330
        \end{Verbatim}
331

332
        The user must consider which approach is the most appropriate for the particular application. The different options have advantages and drawbacks. The first 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 \ident{MATOP\_GET\_DIAGONAL}) or a user-defined one.
333

334
        The second approach (\Verb!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 \ident{MatAXPY} can be very inefficient.
335
If the user knows that the pattern is the same (or a subset), then this can be specified with the function
336
        \findex{STSetMatStructure}
337
        \begin{Verbatim}[fontsize=\small]
338
        STSetMatStructure(ST st,MatStructure str);
339
        \end{Verbatim}
340

341
        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 and in that case it is not possible to use a symmetric definite factorization (Cholesky or ICC).
342

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

345
\subsection{Preserving the Symmetry in Generalized Eigenproblems}
346
\label{sec:symm}
347

2095 jroman 348
        As mentioned in \S\ref{sec:defprob}, 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).
1843 jroman 349

350
        The solution proposed in \slepc is based on selecting different kinds of inner products. Currently, we have the following choice of inner products:
351
  \begin{itemize}
352
    \item Standard Hermitian inner product: $\langle x,y\rangle=x^*y$.
353
    \item $B$-inner product: $\langle x,y\rangle_B=x^*\!B\,y$.
354
  \end{itemize}
355

356
        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 \S\ref{sec:purif}).
357

358
        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$.
359

360
\begin{figure}
361
\centering
362
\begin{tikzpicture}
363
  \begin{scope}[every node/.style={draw,thick,text ragged,font=\sffamily\scriptsize},every to/.style={draw,-stealth,very thick,dashed}]
364
  \node[above right,rounded corners,fill=black!10,text width=2.5cm,inner ysep=2mm] (n1c) at (0,0)
365
       {The user selects the solver};
366
  \node[below right,fill=black!5,text width=2.2cm] (n1) at (0.15,0.03)
367
       {\begin{minipage}{\textwidth}Power/RQI\\Subspace iteration\\Arnoldi\\Lanczos\\
368
       Krylov-Schur\\External solvers\end{minipage}};
369
  \node[above right,rounded corners,fill=black!10,text width=2.5cm,inner ysep=2mm] (n2c) at (6,1.4)
370
       {The user can specify the problem type};
371
  \node[below right,fill=black!5,text width=2.2cm] (n2) at (6.15,1.43)
372
       {\begin{minipage}{\textwidth}General\\Hermitian Definite\\Complex Symmetric
373
       \end{minipage}};
374
  \node[above right,rounded corners,fill=black!10,text width=2.5cm,inner ysep=2mm] (n3c) at (6,-1.0)
375
       {The user can specify the spectral transform};
376
  \node[below right,fill=black!5,text width=2.2cm] (n3) at (6.15,-0.97)
377
       {\begin{minipage}{\textwidth}Shift\\Shift-and-invert\\Cayley\\Fold
378
       \end{minipage}};
379
  \draw (n1c.east) to[out=65,in=120] (n2.west)
380
        (n2.west) to[out=-125,in=-40] (n1c.east);
381
  \draw (n1.east) to[out=25,in=120] (n3.west)
382
        (n3.west) to[out=-115,in=-70] (n1.east);
383
  \end{scope}
384
  \begin{scope}[every node/.style={text centered,text width=2cm,font=\sffamily\scriptsize}]
385
  \node at (4.5,.8) {Appropriate inner product is performed};
386
  \node at (4.5,-1.4) {Appropriate matrix-vector product is performed};
387
  \end{scope}
388
\end{tikzpicture}
389
\caption{\label{fig:abstr}Abstraction used by \slepc solvers.}
390
\end{figure}
391

392
        Internally, \slepc operates with the abstraction illustrated in Figure \ref{fig:abstr}. The operations indicated by dashed arrows are implemented as virtual functions: \ident{IPInnerProduct} and \ident{STApply}. 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 \ident{EPSSetProblemType} (see \S\ref{sec:defprob}).
2091 jroman 393
        In the case of the Cayley transform, \slepc is using $\langle x,y\rangle_{A+\nu B}$ as the inner product for preserving symmetry.
1843 jroman 394

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

2095 jroman 397
        Note that the above discussion is not directly applicable to \ident{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.
398

1843 jroman 399
\subsection{Purification of Eigenvectors}
400
\label{sec:purif}
401

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

404
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 function $\langle x,y\rangle_B$ introduced in \S\ref{sec:symm} 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)$.
405

406
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 \emph{eigenvector purification}.
407

408
Although more elaborate purification strategies have been proposed (usually trying to reduce the computational effort, see \citep{Nour-Omid:1987:HIS} and \citep{Meerbergen:1997:IRA}), 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$.
409

410
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 \ident{EPS\_PGNHEP}, see Table \ref{tab:ptype}.
411