Parallel-in-time MKY. Erlangga

# Parallel-in-time Multilevel Krylov Methods: A Prototype

Yogi A. Erlangga Zayed University, Mathematics and Statistics Department, Abu
Dhabi Campus, P.O. Box 144534, Abu Dhabi, United Arab Emirates ()
<EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

This paper presents a parallel-in-time multilevel iterative method for solving
differential algebraic equation, arising from a discretization of linear time-
dependent partial differential equation. The core of the method is the
multilevel Krylov method, introduced by Erlangga and Nabben [SIAM J. Sci.
Comput., 30(2008), pp. 1572–1595]. In the method, special time restriction and
interpolation operators are proposed to coarsen the time grid and to map
functions between fine and coarse time grids. The resulting Galerkin coarse-
grid system can be interpreted as time integration of an equivalent
differential algebraic equation associated with a larger time step and a
modified $\theta$-scheme. A perturbed coarse time-grid matrix is used on the
coarsest level to decouple the coarsest-level system, allowing full
parallelization of the method. Within this framework, spatial coarsening can
be included in a natural way, reducing further the size of the coarsest grid
problem to solve. Numerical results are presented for the 1- and 2-dimensional
heat equation using simulated parallel implementation, suggesting the
potential computational speed-up of up to 9 relative to the single-processor
implementation and the speed-up of about 3 compared to the sequential
$\theta$-scheme.

###### keywords:

Time-dependent PDEs, heat equation, parallel-in-time, multilevel, Krylov
subspace iterative methods.

65F10, 65M55, 65Y05, 68W10

##  1 Introduction

In the early 2000’s, Lions et al. [17, 18] introduced a time integration
method, called parareal, for numerically solving the linear initial-boundary
value problem:

(1) $\displaystyle\frac{\partial\widehat{u}}{\partial
t}=f\left(\widehat{u},\frac{\partial\widehat{u}}{\partial
x},\frac{\partial^{2}\widehat{u}}{\partial x^{2}}\right),\text{ in
}(0,T]\times\Omega,$

with $x\in\Omega\subset\mathbb{R}^{d}$, the initial condition
$\widehat{u}(0,x)=\widehat{u}_{0}$, and some boundary conditions. In constrast
to the classical approaches that require sequential, time marching solves
starting from the initial condition/solution, parareal computes the solution
over the time interval $\mathcal{T}=(0,T]$ in parallel, in an iterative
fashion.

To describe parareal, consider a spatially discretized version of (1):

(2) $\displaystyle\frac{\partial u(t)}{\partial t}=A_{n}u(t)+g(t),$

where $A_{n}\in\mathbb{R}^{n\times n}$, with $n$ the number of spatial grid
points (or nodal points), $u(t)\in\mathbb{R}^{n}$ is the solution vector at an
instant time $t\in\mathcal{T}$, and $g$ is the boundary condition vector.
Applying the $\theta$-scheme for time integration results in the recurrence

$u_{k}-u_{k-1}=(1-\theta)\Delta tA_{n}u_{k-1}+\theta\Delta tA_{n}u_{k}+\Delta
t((1-\theta)g_{k-1}+\theta g_{k}),$

for $k=1,2,\dots,n_{t}$ and $\theta\in[0,1]$. Here, $u_{k}:=u(t_{k})=u(k\Delta
t)$, where $\Delta t=T/n_{t}$ is the fine time step, used to compute the
highly accurate solution in $\mathcal{T}$, with $n_{t}$ be the number of fine
time steps. Note that $\theta=0,\frac{1}{2},1$ corresponds to the explicit,
Crank-Nicolson, and implicit scheme, respectively. Rearranging the above
equation results in the recurrence

(3) $\displaystyle\psi_{n}u_{k}=\phi_{n}u_{k-1}+\Delta t\bar{g}_{k},$

where $\phi_{n}=I_{n}+(1-\theta)\Delta tA_{n}$, $\psi_{n}=I_{n}-\theta\Delta
tA_{n}$, with $I_{n}$ be the identity matrix of size $n$, and
$\bar{g}_{k}=(1-\theta)g_{k-1}+\theta g_{k}$. Starting from the initial
solution $u_{0}=u(0)$, the solution at any time $t_{k}=k\Delta t$ is computed
sequentially using the recipe

(4) $\displaystyle u_{k}=\psi_{n}^{-1}(\phi_{n}u_{k-1}+\Delta t\bar{g}_{k}).$

Let $\mathcal{T}$ be partitioned into $n_{T}$ subintervals
$\mathcal{T}_{k}=[T_{k-1},T_{k}]$, $k=1,\dots,n_{T}$, with
$|\mathcal{T}_{k}|=\Delta T=\nu\Delta t$, $1<\nu\in\mathbb{N}$, $T_{0}=0$ and
$T_{n_{T}}=T$. The exact111Here, the term“exact” is referred to as the
solution of the fully discretized counterpart of the initial-value problem,
and not that of (1). solution at the time $T_{k}=t_{\nu(k-1)}+\nu\Delta t$ is
given by

(5) $\displaystyle U_{k}:=U(T_{k})=u_{\nu
k}=\Upsilon_{\nu}U_{k-1}+\tilde{g}_{\nu},\quad k=1,2,\dots,n_{T},$

where $\Upsilon_{\nu}=(\psi^{-1}_{n}\phi_{n})^{\nu}$, $\tilde{g}_{\nu}=\Delta
t\psi^{-1}_{n}\displaystyle\sum_{m=0}^{\nu-1}(\phi_{n}\psi^{-1}_{n})^{m-\nu+1}g_{m+\nu(k-1)}$.
This approach performs less time stepping, but requires access to the operator
$\Upsilon_{\nu}$, which is not easily available. This computation can not also
be performed in parallel because of the lack of initial condition at $T_{k-1}$
in the subinterval $\mathcal{T}_{k}$, except for $k=1$.

To allow for parallelization, parareal replaces (5) by the coarse time-step
analogue to (3):

(6) $\displaystyle U_{k}=\Psi^{-1}_{n}(\Phi_{n}U_{k-1}+\Delta T\bar{G}_{k}),$

where $\Phi_{n}=I_{n}+(1-\theta)\Delta TA_{n}$, $\Psi_{n}=I_{n}-\theta\Delta
TA_{n}$, and $\bar{G}_{k}=(1-\theta)G_{k-1}+\theta G_{k}$, and uses it as an
estimator for the initial condition in each subinterval $\mathcal{T}_{k}$.
With an initial-value problem now defined in each subinterval
$\mathcal{T}_{k}$, time integration using the fine time step $\Delta t$ based
on (3) can be computed in parallel. This sequence of computation on the coarse
and fine time grid however only produces the exact solution in
$\mathcal{T}_{1}$, where the exact initial condition is used. A more accurate
solution in the entire computational domain $\mathcal{T}$ is computed by
recalculating the coarse time-grid solution via a sequential correction
procedure, which expectedly produces a more accurate initial condition for
each subinterval $\mathcal{T}_{k}$. The new fine time solution in each time
subinterval $\mathcal{T}_{i}$ is then updated using this improved initial
condition, in parallel. This recipe, summarized in Algorithm 1, is repeated
until a convergence on the fine time grid is achieved.

Algorithm 1 Parareal

1: $U^{0}_{0}=u_{0}$;

2: for $k=1,\dots,n_{T}$ do

3: $U^{0}_{k}=\Psi^{-1}_{n}(\Phi_{n}U^{0}_{k-1}+\Delta T\bar{G}_{k})$;

4: for $\ell=1,2,\dots$ until convergence do

5: for $k=1,\dots,n_{T}$ do

6: for $j=1,\dots,\nu$ do

7:
$u^{\ell}_{\nu(k-1)+j}=\psi^{-1}_{n}(\phi_{n}u^{\ell-1}_{\nu(k-1)+j-1}+\Delta
t\bar{g}_{\nu(k-1)+j-1})$;

8: $U^{\ell}_{0}=U_{\ell-1}^{0}=u_{0}$;

9: for $k=1,\dots,n_{T}$ do

10: $U^{\ell}_{k}=u_{\nu
k}^{0}+\Psi^{-1}\Phi(U^{\ell}_{k-1}-U^{\ell-1}_{k-1})$;

Since the introduction of parareal, interests in time parallelization have
arisen in the last two decades. Parareal, however, may not historically be
regarded as the first work on time parallelization. The earliest work on this
can be traced back to the work by Nievergelt [22] in the 1960s, an era where
parallel computers did not even exist. While parareal and Nievergelt’s method
may look like two distant methods, they can in fact be classified under the
same family: shooting methods; see [10]. There are also methods which were
developed in the context time-space domain decomposition and multigrid methods
[9, 11, 14, 13, 15, 20]. We refer the reader to [10] for a literature account
on time parallelization methods prior to parareal.

There are some interpretations of the initial version of parareal, which
algebraically links parareal with well-known methods for statiationary
problems, and which has led to several alternative parallel-in-time approaches
for solving initial-value problems [12, 10]. One interpretation relates
parareal with a two-(time)grid method, with a nonstandard interpolation and
restriction operator. Such an interpretation has allowed the construction of a
parallel-in-time multigrid method, called MGRIT (multigrid reduction in time),
where time parallelization is done in a multilevel framework [5]. At present,
MGRIT is one of the most heavily tested and popular methods for time
parallelization. Further development of MGRIT includes the inclusion of
spatial coarsening in the method and extension to nonlinear evolution problems
[1, 7, 6, 16].

In this paper, we propose another approach to time parallelization based on
the multilevel Krylov (MK) framework, introduced in [3]. We note here that, as
a multilevel method, MK can not be considered as another multigrid method, due
to, e.g., the nonexistence of smoothing process required for an effective
multigrid method. In fact, MK as outlined in [3] uses GMRES recursively to
solve the coarse-grid system instead of multigrid smoothers such as Jacobi or
Gauss-Seidel. Nevertheless, MK shares components that also are present in
multigrid: fine-coarse grid maps and coarse-grid system. This allows the use
of standard multigrid components in MK [4]. We shall demonstrate in this paper
through numerical tests that a fast convergence of MK can still be achieved by
using a simple technique for the fine-coarse time-grid mapping. Even though
GMRES used on each coarse level increases the computational complexity of the
method, we shall also show that for the model problem considered in the paper
(the heat equation), the parallel-in-time MK method can be made a faster
method than the serial $\theta$-scheme, when space coarsening is also
performed. Due to the author’s limited computational resources, the numerical
tests were, however, only performed in a simulated parallel setting, where a
parallel code is run sequentially on a single processor.

The remainder of this paper is organized as follows. In Section 2, to motivate
the construction of our parallel-in-time multilevel Krylov method, we discuss
the relation between parareal and two-grid methods. Section 3 presents the
construction of the two-level Krylov method for time parallelization and
discusses numerical tests and performance results using simulated parallel
implementation. Extension to the multilevel Krylov method and numerical
results are detailed in Section 4 together with numerical results. In Section
5, we draw some concluding remarks and directions for future’s work.

##  2 Parareal as a two-grid method

To facilitate a link between parareal and a two-grid method, we first rewrite
the recurrence (3) as a large system of equations

(7) $\mathcal{A}_{h}\bm{u}_{h}=\bm{f}_{h},$

where

(8) $\mathcal{A}_{h}=\begin{bmatrix}I_{n}&&&\\\ -\phi_{n}&\psi_{n}&&\\\
&\ddots&\ddots&\\\
&&-\phi_{n}&\psi_{n}\end{bmatrix},\bm{u}_{h}=\begin{bmatrix}u_{0}\\\ u_{1}\\\
\vdots\\\ u_{n_{t}}\end{bmatrix},\text{ and
}\bm{f}_{h}=\begin{bmatrix}u_{0}\\\ \Delta t\bar{g}_{1}\\\ \vdots\\\ \Delta
t\bar{g}_{n_{t}}\end{bmatrix}.$

Here, $\mathcal{A}_{h}$ is an $(n_{t}+1)\times(n_{t}+1)$ block matrix. The
serial solve of (3) is equivalent to solving (7) using block forward
substitution. We shall refer this system to as the fine time-grid system
throughout.

Let the $(n_{T}+1)\times(n_{T}+1)$ block matrix

(9) $\displaystyle\mathcal{A}_{H}=\begin{bmatrix}I_{n}&&&\\\
-\Phi_{n}&\Psi_{n}&&\\\ &\ddots&\ddots&\\\ &&-\Phi_{n}&\Psi_{n}\end{bmatrix}$

be the coarse time-grid matrix analogue to $\mathcal{A}_{h}$. Furthermore, let
$\bm{u}_{H}=[U_{k}]_{k=0}^{n_{T}}$, with $U_{k}=u_{\nu k}$, and
$\bm{f}_{H}=[u_{0}^{\top},\Delta T\bar{G}_{1}^{\top},\dots,\Delta
T\bar{G}_{n_{T}}^{\top}]^{\top}$, with the superscript “⊤” denoting the
transposition of a matrix of vector.

To reformulate parareal as a two-grid method, we define the following specific
choice of components:

  1. 1.

the prolongation matrix $\mathcal{P}$, defined as a $(n_{t}+1)\times(n_{T}+1)$
block matrix where the $(i,k)$-block is zero, except when $i=1+\nu k$,
$k=0,\dots,n_{T}$, where the block $\mathcal{P}_{i,k}=I_{n}$;

  2. 2.

the restriction matrix $\mathcal{R}$, defined as a $(n_{T}+1)\times(n_{t}+1)$
block matrix where the $(k,i)$-block is zero, except when $i=1+\nu k$,
$k=0,\dots,n_{T}$, where the block $\mathcal{R}_{i,k}=\Phi_{n}^{-1}\psi_{n}$;

  3. 3.

the coarse time-grid matrix $\mathcal{A}_{H}$, given by (9). Note that this
coarse time-grid matrix is not of Galerkin type, i.e.,
$\mathcal{A}_{H}\neq\mathcal{R}\mathcal{A}_{h}\mathcal{P}$;

  4. 4.

smoother, which is represented algebraically by the $n_{T}\times n_{T}$
modified block Jacobi matrix

(10) $\mathcal{M}_{J}=\begin{bmatrix}\bar{M}_{J}&&&\\\ &M_{J}&&\\\
&&\ddots&\\\ &&&M_{J}\end{bmatrix},$

where the diagonal blocks

(11) $\bar{M}_{J}=\begin{bmatrix}I_{n}&&&\\\ -\phi_{n}&\psi_{n}&&\\\
&\ddots&\ddots&\\\ &&-\phi_{n}&\psi_{n}\end{bmatrix}\text{ and
}M_{J}=\begin{bmatrix}\psi_{n}&&&\\\ -\phi_{n}&\ddots&&\\\ &\ddots&\ddots&\\\
&&-\phi_{n}&\psi_{n}\end{bmatrix}$

are an $(\nu+1)\times(\nu+1)$ and $\nu\times\nu$ block matrix, respectively.
Furthermore, let $\mathcal{I}_{J}$ be the $n_{t}\times n_{t}$ modified block
identity matrix with zero matrix on the block that corresponds to the $u_{\nu
k}$ unknown, for $k=0,1,\dots,n_{T}$.

It can be shown that the following two-grid procedure produces parareal
iterands $u_{k}^{\ell}$ and $U_{k}^{\ell}$:

  1. 1.

Set the initial approximate solution $\bm{u}^{0}_{h}$ using the following
steps:

    1. 1.1.

Solve the coarse time-grid problem:

$\mathcal{A}_{H}\bm{u}_{H}^{0}=\bm{f}_{H}.$

    2. 1.2.

Prolongate the coarse time-grid solution to the fine time-grid:

$\bm{u}_{h}^{0}=\mathcal{P}\bm{u}^{0}_{H}.$

  2. 2.

For $\ell=1,2,\dots$ until convergence:

    1. 2.1.

Pre-smoothing:

$\widetilde{\bm{u}}_{h}^{\ell-1}=\bm{u}^{\ell-1}_{h}+\mathcal{I}_{J}\mathcal{M}^{-1}_{J}(\bm{f}_{h}-\mathcal{A}_{h}\bm{u}_{h}^{\ell-1}).$

    2. 2.2.

Restriction:

$\displaystyle\widetilde{\bm{r}}_{H}^{\ell-1}$
$\displaystyle=\mathcal{R}\left(\bm{f}_{h}-\mathcal{A}_{h}\widetilde{\bm{u}}_{h}^{\ell-1}\right).$

    3. 2.3.

Coarse time-grid solve:

$\bm{d}_{H}^{\ell-1}=\mathcal{A}_{H}^{-1}\widetilde{\bm{r}}_{H}^{\ell-1}.$

    4. 2.4.

Prolongation and correction:

$\bm{u}_{h}^{\ell}=\widetilde{\bm{u}}_{h}^{\ell-1}+\mathcal{P}\bm{d}_{H}^{\ell-1}.$

The above two-grid method produces a sequence of error vectors
$\bm{e}_{h}^{\ell}$ that satisfies the relation

(12) $\bm{e}_{h}^{\ell}=\mathcal{K}_{MG}\bm{e}_{h}^{\ell-1},$

where
$\mathcal{K}_{MG}:=(\mathcal{I}_{h}-\mathcal{P}\mathcal{A}_{H}^{-1}\mathcal{R}\mathcal{A}_{h})(\mathcal{I}_{h}-\mathcal{I}_{J}\mathcal{M}_{J}^{-1}\mathcal{A}_{h})$.

In $\mathcal{K}_{MG}$, the term
$\mathcal{I}_{h}-\mathcal{P}\mathcal{A}_{H}^{-1}\mathcal{R}\mathcal{A}_{h}=:Q_{D,h}$,
which represents the coarse-grid correction, is similar to the deflation
operator [21, 8, 2], used to accelerate a Krylov iterative method. However,
like in standard multigrid methods, deflation uses Galerkin coarse-grid matrix
and some conditions for $\mathcal{R}$ and $\mathcal{P}$ such that the Galerkin
coarse-grid matrix is nonsingular. In this case, applying $\mathcal{Q}_{D,h}$
on $\mathcal{A}_{h}$ deflates or shifts some eigenvalues of $\mathcal{A}_{h}$
to zero, resulting in a convergence acceleration of a Krylov method [19]. The
number of eigenvalues of $\mathcal{A}_{h}$ deflated to zeros typically depends
on the rank condition of the matrix $\mathcal{R}$ and $\mathcal{P}$.

Due to a non-Galerkin coarse-grid approach, $Q_{D,h}$ in this case does not
play the role one typically expects from deflation: that near-zero eigenvalues
in the spectrum of $\mathcal{A}_{h}$ are deflated so that they are no longer
present in the spectrum of $\mathcal{A}_{h}Q_{D,h}$, leading to improved
convergence.

In [24], several other variants of shifting-based preconditioning operators
derived from multigrid and domain decomposition methods are derived and
analyzed.

##  3 Two-level Krylov for time parallelization

###  3.1 General approach

Consider again the partition of the time interval $\mathcal{T}=[0,T]$ into
$n_{T}$ uniform subintervals, $\mathcal{T}_{k}=[T_{k-1},T_{k}]$,
$k=1,\dots,n_{T}$, with $T_{0}=t_{0}=0$ and $T_{n_{T}}=t_{n_{t}}=T$, such that
$\Delta T=T/n_{T}=\nu T/n_{t}=\nu\Delta t$, $1<\nu\in\mathbb{N}$. Using this
partitioning, the linear system (7) can be rewritten as

(13) $\displaystyle\begin{bmatrix}\mathcal{A}^{h}_{0,0}&&&\\\
\mathcal{A}^{h}_{1,0}&\mathcal{A}^{h}_{1,1}&&\\\ &\ddots&\ddots&\\\
&&\mathcal{A}^{h}_{n_{T},n_{T}-1}&\mathcal{A}^{h}_{n_{T},n_{T}}\end{bmatrix}\begin{bmatrix}\bm{u}^{h}_{0}\\\
\bm{u}^{h}_{1}\\\ \vdots\\\
\bm{u}^{h}_{n_{T}}\end{bmatrix}=\begin{bmatrix}\bm{f}^{h}_{0}\\\
\bm{f}^{h}_{1}\\\ \vdots\\\ \bm{f}^{h}_{n_{T}}\end{bmatrix},$

where $\mathcal{A}_{0,0}^{h}=I_{n}$, $\bm{u}_{0}^{h}=u_{0}$, and, for
$k=1,\dots,n_{T}$, $\bm{u}^{h}_{k}=[u^{h}_{\nu(k-1)+i}]_{i=1,\dots,\nu}$ and
$\mathcal{A}^{h}_{k,k}$ is a $\nu\times\nu$-block matrix.

The two-level Krylov method [3] for (13) is a Krylov method applied to the
preconditioned system

(14)
$\displaystyle\mathcal{A}_{h}\mathcal{Q}_{h}\widetilde{\bm{u}}_{h}=\bm{f}_{h},$

where

(15)
$\displaystyle\mathcal{Q}_{h}=\mathcal{I}_{h}-\mathcal{Z}\mathcal{A}_{H}^{-1}\mathcal{Y}^{\top}\mathcal{A}_{h}+\mu\mathcal{Z}\mathcal{A}_{H}^{-1}\mathcal{Y}^{\top},\quad\mathcal{A}_{H}=\mathcal{Y}^{\top}\mathcal{A}_{h}\mathcal{Z},$

with $\mu\in\mathbb{R}$. Here, both $\mathcal{A}_{H}$ and $\mathcal{Q}_{h}$
are assumed to be nonsingular so that
$\widetilde{\bm{u}}_{h}=\mathcal{Q}_{h}^{-1}\bm{u}^{h}$ is defined.

The choice $\mu=0$ in (15) corresponds to the deflation operator [8, 24]. This
choice however leads to a singular $\mathcal{Q}_{h}$ and is known to be not a
good choice for a multilevel implementation, due to “reappearance” of near-
zero eigenvalues in the spectrum of $\mathcal{A}_{h}\mathcal{Q}_{h}$ under an
inaccurate inversion of the coarse time-grid matrix $\mathcal{A}_{H}$. Since
near-zero eigenvalues generally are responsible for slow convergence, an
effective deflation method requires an exact solve of the coarse time-grid
system. Such a requirement limits the deflation approach not only to a two-
level setting – unless the coarse time-level system is solved very accurately
using a multilevel technique – but also the size of the coarse time-level
system to solve.

Alternatively, as used in [3], the parameter $\mu$ is set equal to the largest
eigenvalues of $\mathcal{A}_{h}$; $\mu=\lambda_{\max}(\mathcal{A}_{h})$. With
this choice, the action of $\mathcal{Q}_{h}$ “shifts” some near-zero
eigenvalues of $\mathcal{A}_{h}$ to the $\lambda_{\max}(\mathcal{A}_{h})$.
Under an inaccurate inversion of $\mathcal{A}_{h}$, some eigenvalues of
$\mathcal{A}_{h}\mathcal{Q}_{h}$ are expected to appear slightly away from
$\lambda_{\max}(\mathcal{A}_{h})$, but remain far from zero. This property
allows the use of an inaccurate inversion of $\mathcal{A}_{H}$ using, for
instance, an iterative method, making the method practical. Furthermore, with
the spectrum of $\mathcal{A}_{h}\mathcal{Q}_{h}$ distributed away from zero,
we may expect an improved convergence of a Krylov method.

The matrices $\mathcal{Y},\mathcal{Z}\in\mathbb{R}^{n_{t}n\times r}$ in (15)
are called the deflation matrices, assumed to be of rank $r<n_{t}n$ and such
that $\mathcal{A}_{H}=\mathcal{Y}^{T}\mathcal{A}_{h}\mathcal{Z}$ is
nonsingular. From a theoretical view point, $\mathcal{Y}$ and $\mathcal{Z}$
are chosen such that the spectrum of $\mathcal{A}_{h}\mathcal{Q}_{h}$ has no
near-zero eigenvalues. An example of such a choice is when the columns of
$\mathcal{Y}$ and $\mathcal{Z}$ are left and, respectively, right eigenvectors
associated with the $r$ near-zero eigenvalues of $\mathcal{A}_{h}$. In this
case, $r$ near-zero eigenvalues are shifted to $\lambda_{\max}$, with the
other $n_{t}n-r$ eigenvalues remaining untouched. Furthermore, a larger $r$
means more near-zero eigenvalues are shifted, resulting in a more clustered
eigenvalues of $\mathcal{A}_{h}\mathcal{Q}_{h}$, and expectedly faster
convergence. Computing a set of $r$ eigenvectors is however very costly,
especially when the size of $\mathcal{A}_{h}$ is large. Since $\mathcal{Y}$
and $\mathcal{Z}$ algebraically resemble integrid operators in multigrid (cf.
$\mathcal{K}_{MG}$ below Eq. (12)), in practice, we shall use the easier-to-
construct restriction and interpolation (or integrid) matrix known in
multigrid.

A two-level Krylov implementation is presented in Algorithm 2, where flexible
GMRES (FGMRES) [23] is used as the basis.

Algorithm 2 FGMRES-based Multilevel Krylov implementation

1: Choose $\bm{u}_{h}^{0}$ and $\mu$.

2: Compute $\bm{r}_{h}^{0}=\bm{f}_{h}-\mathcal{A}_{h}\bm{u}_{h}^{0}$,
$\beta=\|\bm{r}_{h}^{0}\|_{2}$, and $\bm{v}_{h}^{1}=\bm{r}_{h}^{0}/\beta$.

3: for $j=1,2,\dots,J$ do

4: $\bm{x}_{h}^{j}:=\mathcal{Q}_{h}\bm{v}_{h}^{j}$;

5: $\bm{w}_{h}:=\mathcal{A}_{h}\bm{x}_{h}^{j}$;

6: for $l=1,\dots,j$ do

7: $h_{l,j}:=(\bm{w}_{h},\bm{v}_{h}^{j})$;

8: $\bm{w}_{h}:=\bm{w}_{h}-h_{l,j}\bm{v}_{h}^{j}$;

9: $h_{j+1,j}:=\|\bm{w}_{h}\|_{2}$;

10: $\bm{v}_{h}^{j+1}:=\bm{w}_{h}/h_{j+1,j}$;

11: Set $\mathcal{X}_{J}=[\bm{x}_{h}^{1}\dots\bm{x}_{h}^{J}]$ and
$\hat{H}_{J}=[h_{i,j}]_{1\leq i\leq j+1;1\leq j\leq J}$;

12: Compute
$\bm{y}_{J}=\arg\min_{\bm{y}}\|\beta\bm{e}_{1}-\hat{H}_{J}\bm{y}\|_{2}$ and
$\bm{u}_{h}^{J}=\bm{u}_{h}^{0}+\mathcal{X}_{J}\bm{y}_{J}$.

In Algorithm 2, Step 4 corresponds to the (right-preconditioned) two-level
Krylov step. In this step, the matrix $\mathcal{Q}_{h}$ is never explicitly
formed. Instead, its action on the vector $\bm{v}_{h}^{j}$ is performed
implicitly using the following substeps:

$\displaystyle\text{\small 4.1:}~{}\bm{s}_{h}:=\mathcal{A}_{h}\bm{v}_{h}^{j};$
$\displaystyle\text{\small
4.2:}~{}\widehat{\bm{v}}_{H}:=\mathcal{Y}^{\top}(\bm{s}_{h}-\mu\bm{v}_{h}^{j});$
$\displaystyle\text{\small
4.3:}~{}\widetilde{\bm{v}}_{H}:=\mathcal{A}_{H}^{-1}\widehat{\bm{v}}_{H};$
$\displaystyle\text{\small
4.4:}~{}\bm{t}_{h}:=\mathcal{Z}\widetilde{\bm{v}}_{H};$
$\displaystyle\text{\small 4.5:}~{}\bm{x}_{h}^{j}:=\bm{v}_{h}^{j}-\bm{t}_{h};$

###  3.2 Intergrid matrices

To construct the matrix $\mathcal{Y}$ and $\mathcal{Z}$, we shall repurpose
the multigrid intergrid matrices. For instance, following MGRIT [5], we can
set $\mathcal{Y}$ as a matrix that corresponds to a direct injection process
of a fine time-grid quantity to the coarse time-grid, i.e.,

(16) $\displaystyle\mathcal{Y}=\begin{bmatrix}0\\\
\mathcal{I}_{c}\end{bmatrix}$

and $\mathcal{Z}$ the ideal interpolation:

(17)
$\displaystyle\mathcal{Z}=\begin{bmatrix}-(\mathcal{A}^{h}_{ff})^{-1}\mathcal{A}^{h}_{fc}\\\
\mathcal{I}_{c}\end{bmatrix},$

where $\mathcal{A}^{h}_{ff}$ and $\mathcal{A}^{h}_{fc}$ are submatrices of
$\mathcal{A}_{h}$ associated with partitioning of unknowns into the fine time-
grid and course time-grid unknowns:

$\mathcal{A}_{h}=\begin{bmatrix}\mathcal{A}_{ff}^{h}&\mathcal{A}_{fc}^{h}\\\
\mathcal{A}_{cf}^{h}&\mathcal{A}_{cc}^{h}\end{bmatrix}.$

Consequently,
$\mathcal{A}_{H}=\mathcal{Y}^{\top}\mathcal{A}_{h}\mathcal{Z}=\mathcal{A}_{cc}^{h}-\mathcal{A}_{cf}^{h}(\mathcal{A}_{ff}^{h})^{-1}\mathcal{A}_{fc}^{h}$.
As in parareal, this expensive-to-invert $\mathcal{A}_{H}$ can be replaced by
the coarse time-grid analogue of $\mathcal{A}_{h}$, associated with a
sequential solve of time marching problem with the larger time step $\Delta
T$, given by (9).

For the two-level Krylov method, we however consider the specific construction
of the intergrid operators based on the $(n_{T}+1)\times(n_{T}+1)$-block
diagonal matrix

(18)
$\displaystyle\quad\mathcal{Y}=blockdiag(\mathcal{Y}_{k}),\quad\mathcal{Z}=blockdiag(\mathcal{Z}_{k}),$

where the submatrices $\mathcal{Y}_{k}$ and $\mathcal{Z}_{k}$,
$k=0,\dots,n_{T}$ are full-rank. The size of $\mathcal{Y}_{k}$ and
$\mathcal{Z}_{k}$ is assumed to be compatible under left- and right-
multiplication with $\mathcal{A}^{h}_{k,k}$, so that the coarse time-grid
matrix $\mathcal{A}_{H}$ is defined. In this case,

(19) $\displaystyle\mathcal{A}_{H}=\begin{bmatrix}\mathcal{A}^{H}_{0,0}&&&\\\
\mathcal{A}^{H}_{1,0}&\mathcal{A}^{H}_{1,1}&&\\\ &\ddots&\ddots&\\\
&&\mathcal{A}^{H}_{n_{T},n_{T}-1}&\mathcal{A}^{H}_{n_{T},n_{T}}\end{bmatrix},$

where the diagonal blocks
$\mathcal{A}^{H}_{k,k}=\mathcal{Y}_{k}^{\top}\mathcal{A}^{h}_{k,k}\mathcal{Z}_{k}$,
$k=0,\dots,n_{T}$, are nonsingular, and
$\mathcal{A}^{H}_{k,k-1}=\mathcal{Y}_{k}^{\top}\mathcal{A}^{h}_{k,k-1}\mathcal{Z}_{k-1}$,
$k=1,\dots,n_{T}$.

In the sequel, we discuss two constructions, which lead to time coarsening and
time-space coarsening two-level Krylov methods.

####  3.2.1 Case 1: Time coarsening

One possible setting for the intergrid matrices, which corresponds to the
(pure) time (or T-coarsening, throughout), is based on $\mathcal{Z}_{0}=I_{n}$
and

(20) $\displaystyle\mathcal{Z}_{k}=[I_{n}\,I_{n}\,\dots\,I_{n}]^{\top},\quad
k=1,\dots,n_{T},$

where $\mathcal{Z}_{k}$ is a $\nu\times 1$ block matrix, and the setting
$\mathcal{Y}_{k}=\mathcal{Z}_{k}$. Submatrices in the coarse-grid matrix
$\mathcal{A}_{H}$ (19) can be computed explicitly and are given by

(21) $\displaystyle\mathcal{A}^{H}_{0,0}$ $\displaystyle=I_{n},$ (22)
$\displaystyle\mathcal{A}^{H}_{kk}$
$\displaystyle=\nu\psi_{n}-(\nu-1)\phi_{n}=I_{n}-(\nu-\theta)\Delta
tA_{n}=I_{n}-(1-\tilde{\theta})\Delta TA_{n}=:\widetilde{\Psi}_{n},$ (23)
$\displaystyle\mathcal{A}^{H}_{k,k-1}$
$\displaystyle=-\phi_{n}=-(I_{n}+\theta\Delta
tA_{n})=-(I_{n}+\tilde{\theta}\Delta TA_{n})=:-\widetilde{\Phi}_{n},$

for $k=1,\dots,n_{T}$, where $\tilde{\theta}=\theta/\nu\in[0,1/2]$.

Note that the equations (21)–(23) differ slightly from parareal’s and MGRIT’s
coarse time-grid problem in the parameter $\theta$ (cf. Eq. (9)). With
(21)–(23), the coarse time-grid system can be interpreted as a modfied
$\theta$-scheme with a larger time step $\Delta T=\nu\Delta t$ (see Figure 1).
A meaningful time coarsening requires $\nu\leq 2$, resulting in
$\tilde{\theta}<\theta$, which may lead to stability issues on the coarse grid
for large $\Delta T$. The case $\tilde{\theta}=\theta$ (as in parareal) only
occurs when $\theta=0$ (the explicit method).

Figure 1: Time coarsening, with $\Delta T=\nu\Delta t$, with $\nu=3$.

####  3.2.2 Case 2: Time-space coarsening

With the computational work dominated by the inversion of the implicit term
$\psi_{n}$, the total computational work of the sequential $\theta$-scheme is
of order $n_{t}C_{inv}(n)$, with $C_{inv}$ the complexity associated with the
action of $\psi_{n}^{-1}$ on a vector (for 2D cases, typically
$C^{h}_{inv}(n)\sim n^{2}$). For the GMRES-based two-level Krylov method, the
order of computational work is

$W_{2L}\sim N_{iter}\frac{n_{t}}{\nu}C^{H}_{inv},$

where $N_{iter}$ is the number of GMRES iterations needed to reach convergence
and $C^{H}_{inv}$ is the complexity associated with the action of
$\tilde{\Phi}_{n}^{-1}$ in the coarse time-grid system, which is the same as
$\psi_{n}^{-1}$. Therefore, for a two-level method to be competitive, the
ratio $N_{iter}/\nu$ has to be made less than 1, via an aggressive time
coarsening with large $\nu$. This is however not easy to attain solely using
an aggresive time coarsening. As we shall see in the numerical example, a
large $\nu$ used in an aggressive time coarsening slows down the convergence
of GMRES (so, $N_{iter}/\nu>1$) that even makes the method impractical, due to
a significant increase in the number of Arnoldi vectors to store.

A reduction in the computational complexity can be attained via a combination
of time and space (TS) coarsening, resulting in a much smaller coarse time-
space-grid system. For the construction of the time-space coarsening, we
consider the intergrid matrices (18), but now with

(24)
$\displaystyle\mathcal{Z}_{k}=[Z\,Z\,\dots\,Z]^{\top}\quad\text{and}\quad\mathcal{Y}_{k}=[Y\,Y\,\dots\,Y]^{\top},$

where $\mathcal{Z}_{k},\mathcal{Y}_{k}$ are a $\nu\times 1$-block diagonal
matrix, with a full rank $Z\in\mathbb{R}^{n\times m}$, $m<n$. In this case,
the submatrices in $\mathcal{A}_{H}$ are given by

(25) $\displaystyle\mathcal{A}^{H}_{0,0}$
$\displaystyle=Y^{\top}Z\in\mathbb{R}^{m\times m},$ (26)
$\displaystyle\mathcal{A}^{H}_{kk}$
$\displaystyle=Y^{\top}Z-(1-\tilde{\theta})\Delta
TY^{\top}AZ\in\mathbb{R}^{m\times m},$ (27)
$\displaystyle\mathcal{A}^{H}_{k,k-1}$
$\displaystyle=-(Y^{\top}Z+\tilde{\theta}\Delta
TY^{\top}AZ)\in\mathbb{R}^{m\times m},$

for $k=1,\dots,n_{T}$.

For $Y^{\top}$ and $Z$, we use simple-to-construct agglomeration matrices,
which are not a typical choice for multigrid but work very well for multilevel
Krylov [3]. We consider a partition of the spatial domain
$\Omega\subset\mathbb{R}^{d}$ into $N$ nonoverlapping subdomains $\Omega_{j}$,
such that $\Omega=\bigcup_{j=1}^{N}\Omega_{j}$. Let $\bm{x}_{i}\in\Omega$ be a
spatial gridpoint, and let $\mathcal{I}=\\{i|\bm{x}_{i}\in\Omega\\}$ be the
index set. Furthermore, let
$\mathcal{I}_{j}=\\{i\in\mathcal{I}|\bm{x}_{i}\in\Omega_{j}\\}$. Figure 2
illustrates this partitioning for $\Omega\subset\mathbb{R}^{2}$, with
$\Omega_{j}\subset\Omega$ containing 12 gridpoints with indices in
$\mathcal{I}_{j}$.

Figure 2: Agglomeration of 12 fine spatial gridpoints in
$\Omega_{j}\subset\Omega\subset\mathbb{R}^{2}$ into 1 coarse spatial
gridpoint.

The matrix $Z=[z_{ij}]$ is defined as follows:

(28) $\displaystyle z_{ij}=\begin{cases}1,&i\in\mathcal{I}_{j}\\\
0,&i\notin\mathcal{I}_{j}\end{cases}.$

The matrix $Y=[y_{ji}]$ is given by

(29) $\displaystyle
y_{ji}=\begin{cases}1/|\mathcal{I}_{j}|,&i\in\mathcal{I}_{j}\\\
0,&i\notin\mathcal{I}_{j}\end{cases}.$

A specific example in 1D is shown in Figure 3, where standard spatial
coarsening, i.e., where the coarse space grid is generated by agglomerating 2
fine spatial gridpoints, is used. One can show that $Y^{\top}Z=I_{m}$, with
$m=4$ for this particular example. Furthermore, $Y^{\top}AZ$ is associated
with a direct discretization with the coarse grid size $\Delta
X=\sqrt{2}\Delta x$.

Figure 3: 1D agglomeration using standard spatial coarsening. Left: the
construction of $Z$. Right: the construction of $Y$.

###  3.3 Two-level convergence

To have insights on the convergence of the two-level Krylov method, we shall
consider the heat equation

(30) $\displaystyle\frac{\partial\hat{u}}{\partial
t}=\tau\sum_{i=1}^{d}\frac{\partial^{2}\hat{u}}{\partial x_{i}^{2}},\text{ in
}\Omega\subset\mathbb{R}^{d}\quad(d\in\\{1,2\\}),$

with the initial condition $u(0,\cdot)=u_{0}$ and the boundary conditions
$u(\cdot,\bm{x})=0$, for $\bm{x}\in\partial\Omega$. The spatial derivatives
are discretized by the second-order finite difference scheme. The resultant
differential algebraic equation is integrated using the $\theta$-scheme, with
$\theta=0.5$ (Crank-Nicolson).

We shall first discuss spectra of the two-level matrix
$\mathcal{A}_{h}\mathcal{Q}_{h}$ for the 1-dimensional case. By a direct
computation it can be shown that

(31)
$\displaystyle\mathcal{A}_{h}\mathcal{Q}_{h}=\begin{bmatrix}\mathcal{A}^{h}_{0,0}\mathcal{Q}^{h}_{0}&&&\\\
\star&\mathcal{A}^{h}_{1,1}\mathcal{Q}^{h}_{1}&&\\\ \star&\star&\ddots&\\\
\star&\star&\star&\mathcal{A}^{h}_{n_{T},n_{T}}\mathcal{Q}^{h}_{n_{T}}\end{bmatrix},$

where
$\mathcal{Q}^{h}_{k}=\mathcal{I}^{h}_{k}-\mathcal{Z}_{k}(\mathcal{A}^{H}_{k,k})^{-1}\mathcal{Y}^{\top}_{k}\mathcal{A}^{h}_{k,k}+\mu\mathcal{Z}_{k}(\mathcal{A}^{H}_{k,k})^{-1}\mathcal{Y}^{\top}_{k}$.
Consequently,

(32)
$\displaystyle\sigma(\mathcal{A}_{h}\mathcal{Q}_{h})=\bigcup_{k=0}^{n_{T}}\sigma(\mathcal{A}^{h}_{k,k}\mathcal{Q}^{h}_{k}),$

where $\sigma(\cdot)$ denotes the spectrum of the matrix in the argument.
Thus, the eigenvalues can be computed block-wise. Furthermore, for uniform
coarsening (with fixed $\Delta T$), the block matrix
$\mathcal{A}^{h}_{k,k}\mathcal{Q}^{h}_{k,k}$, $k=1,\dots,n_{T}$, is constant.
Thus, the complete spectrum can be obtained by computing eigenvalues of
$\mathcal{A}^{h}_{k,k}\mathcal{Q}^{h}_{k,k}$, $k=0,1$.

####  3.3.1 1-dimensional case

Figure 4 shows computed eigenvalues (in black) of the two-level matrix
$\mathcal{A}_{h}\mathcal{Q}_{h}$ using the T-coarsening for the time
coarsening parameter $\nu\in\\{2,4,8\\}$ and the shift parameter $\mu=1$. The
spectrum of the original matrix $\mathcal{A}_{h}$, shown in blue in Figure 4,
lies along the real axis between 1 and 1.5. For $\nu=2$, the spectrum of
$\mathcal{A}_{h}\mathcal{Q}_{h}$ splits into two clusters: one at 1, which
corresponds to the shift $\mu=1$, and the other cluster spread between 1.5 and
2 in the real axis. For $\nu=4$ and $8$, the spectra are generally less
clustered with eigenvalues split in 4 and 8 clusters in the complex plane,
respectively, with one cluster exactly at 1. With less clustered eigenvalues,
we should expect slower convergence rate of the two-level method, when a more
aggresive time coarsening approach is used.

(a) $\nu=2$

(b) $\nu=4$

(c) $\nu=8$

Figure 4: Spectrum of $\mathcal{A}_{h}\mathcal{Q}_{h}$ (in black), from the
1-dimensional heat equation. The two-level setting uses time coarsening
intergrid operators (20) and $\mathcal{Y}_{k}=\mathcal{Z}_{k}$. The spectrum
of $\mathcal{A}_{h}$ is plotted in blue. For the discretization, $n_{x}=63$,
$n_{t}=32$, the Courant number $C=2k\Delta t/\Delta x^{2}=0.32$, and
$\theta=0.5$.

Spectra of the two-level matrix with the time-space (TS) coarsening and MGRIT
integrid components are shown in Figure 5 and 6, respectively. For the TS-
coarsening, the spectra are similar to the one with the T-coarsening,
characterized by clusters of eigenvalues in the complex plane. Some
eigenvalues are however spread near one. This cluster of eigenvalues tend to
spread farther from one with an increase in $\nu$. The MGRIT-based two-level
matrix leads to a completely different eigenvalue distribution from the T- and
TS-coarsening. In this case, the eigenvalues lie inside a convex hull that
contains the eigenvalues of the original matrix $\mathcal{A}_{h}$ (recall that
in this case, the matrix $\mathcal{A}_{h}$ is rearranged according the fine-
coarse partition).

(a) $\nu=2$

(b) $\nu=4$

(c) $\nu=8$

Figure 5: Spectrum of $\mathcal{A}_{h}\mathcal{Q}_{h}$ (in black), from
1-dimensional heat equation. The two-level setting uses time-space coarsening
integrid operators (24), (28) and (29). The spectrum of $\mathcal{A}_{h}$ is
plotted in blue. For the discretization, $n_{x}=63$ and $n_{t}=32$,
corresponding to the Courant number $C=2k\Delta t/\Delta x^{2}=0.32$, and
$\theta=0.5$.

(a) $\nu=2$

(b) $\nu=4$

(c) $\nu=8$

Figure 6: Spectrum of $\mathcal{A}_{h}\mathcal{Q}_{h}$ (in black), based on
1-dimensional heat equation. The two-level setting uses MGRIT intergrid
operators (16) and (17), and the coarse-grid matrix (9). The spectrum of
$\mathcal{A}_{h}$ is plotted in blue. For the discrerization, $n_{x}=64$ and
$n_{t}=32$, corresponding to the Courant number $C=2\tau\Delta t/\Delta
x^{2}=0.32$, and $\theta=0.5$.

With this type of spectra, it is not easy to deduce if a two-level Krylov
approach will bring a significant convergence improvement, especially when
$\nu>2$, where the spectra become even much less clustered. Despite this
spread of eigenvalues, the computed condition numbers $\kappa$ of the two-
level matrix, presented in Table 1, are significantly lower than that of
$\mathcal{A}_{h}$. Furthermore, the table suggests that a better conditioning
can be achieved by the T-coarsening strategy. The T-coarsening also leads to
the condition number of the two-level matrix which is practically independent
of the Courant number $C=2\tau\Delta t/\Delta x^{2}$. For all coarsening
strategies, the conditioning deteriorates with increasing $\nu$, with the
T-coarsening showing the least deterioration. This suggests that aggresive
time coarsening (large $\nu$) may not be a good scenario for achieving a low
number of iterations to reach convergence. For the TS-coarsening, because
$\text{rank}Z<\text{rank}I_{n}=n$, we expect less near-zero eigenvalues being
shifted to 1, leading to slower convergence than with the T-coarsening.

Table 1: Condition numbers of $\mathcal{A}_{h}\mathcal{Q}_{h}$, based on the 1-dimensional heat equation. For the discretization, $n_{x}=64$, with $\Delta t$ chosen according to the Courant number $C=2\tau\Delta t/\Delta x^{2}$, $\theta=0.5$. | $C=\frac{2\tau\Delta t}{\Delta x^{2}}$ | $\kappa(\mathcal{A}_{h}\mathcal{Q}_{h})$ | $\kappa(\mathcal{A}_{h})$
---|---|---|---
| $\nu=2$ | $\nu=4$ | $\nu=8$ |
T-coarsening | 0.16 | 2.62 | 4.05 | 8.69 | 82.95
0.32 | 2.62 | 4.05 | 8.69 | 42.39
0.64 | 2.62 | 4.05 | 8.67 | 22.11
TS-coarsening | 0.16 | 17.20 | 19.13 | 24.05 | 82.95
0.32 | 8.82 | 9.85 | 12.87 | 42.39
0.64 | 4.84 | 5.53 | 8.87 | 22.11
MGRIT-based intergrid operators | 0.16 | 2.64 | 4.05 | 10.76 | 82.95
0.32 | 2.69 | 5.52 | 10.73 | 42.39
0.64 | 2.85 | 5.59 | 10.74 | 22.11

To see how the above results on condition numbers are translated to the actual
convergence, we performed a numerical test by solving the two-level system
(14) using GMRES, with the fine (level 1) time-space discretization
corresponding to the Courant number $C=2\tau\Delta t/\Delta x^{2}=0.64$. In
this and all other numerical tests, we set $\bm{u}^{0}_{h}=0$ and declare
convergence at the $J$-th iteration if the following criteria is met:

$\|\bm{f}_{h}-\mathcal{A}_{h}\bm{u}^{J}_{h}\|_{2}/\|\bm{f}_{h}\|_{2}<10^{-6}.$

We note here in passing that in all numerical tests presented here, at
convergence, the norm of difference between the solutions of the two- or
multilevel Krylov method and the $\theta$-scheme is of order of $10^{-5}$.

The convergence results in terms of iteration counts are shown in Table 2. As
expected from its low condition numbers, the T-coarsening consistently leads
to a more superior approach than the other two approaches. The plain use of
MGRIT components does not generally lead to a better approach than the other
options, except when the standard time coarsening ($\nu=2$) is used. For the
latter case, MGRIT-based approach converges faster than the TS-coarsening
approach.

Table 2: Numbers of GMRES iterations for solving the 1-dimensional heat equation with two-level Krylov: $C=2\tau\Delta t/\Delta x^{2}=0.64$, $\theta=0.5$ |  | T-coarsening | TS-coarsening | MGRIT-based
---|---|---|---|---
|  | $\nu:$ | $\nu:$ | $\nu:$
$n_{x}$ | $n_{t}$ | 2 | 4 | 8 | 2 | 4 | 8 | 2 | 4 | 8
128 | 128 | 6 | 16 | 35 | 19 | 24 | 38 | 12 | 25 | 50
256 | 512 | 6 | 15 | 32 | 18 | 22 | 36 | 11 | 23 | 49
512 | 2048 | 5 | 14 | 30 | 18 | 22 | 34 | 11 | 22 | 43
1024 | 8192 | 5 | 13 | 27 | 17 | 21 | 32 | 10 | 22 | 42

Figure 7 shows iteration counts to reach convergence as a function of the
Courant number $C$, with $n_{x}=512$ and $1024$. Again, the T-coarsening
approach results in the fastest two-level method, with iteration counts to
convergence practically independent of the Courant number considered in this
test. The MGRIT-based two-level method shows a slight increase in the
iteration count as $C$ increases. The TS-coarsening, while does not result in
the fastest method in terms of iteration counts, shows a significantly
improved convergence in cases where $C$ is high. As suggested by the increased
condition numbers with larger $\nu$, aggressive coarsening affects the
convergence negatively.

Figure 7: The effect of the Courant number $C=2\tau\Delta t/\Delta x^{2}$ on
the number of iterations. The test problem is the 1-dimensional heat equation.

####  3.3.2 2-dimensional case

In this section, we discuss computational performance of the 2-level Krylov
method for solving the 2-dimensional heat equation. We use the standard time
coarsening $(\nu=2)$ in the T-coarsening and MGRIT-based approach and
additional standard space coarsening in the TS-coarsening. The computation is
performed using Matlab on an Intel Core i7-7th Gen machine with 24 Gbytes of
RAM. For the tests, the numbers of spatial gridpoints in the $x$\- and
$y$-direction, denoted by $n_{x}$ and $n_{y}$ respectively, are set to be
equal, with $n_{x}\in\\{128,192\\}$. Even though the method is implemented
sequentially, all parallelizable parts in the method are coded in a way that
mimics parallel computations. The computational time thus is measured as if
such a task is performed in parallel in multiple processors. This simulated
parallel computation only gives us an approximation to the actual parallel
computational time, as practically there is no communication time.

Convergence results of the sequential (single processor) two-level Krylov
method in terms of numbers of iterations for varying values of $C=3\tau\Delta
t/\Delta x^{2}$ are shown in Figure 8:top-left. As seen from the figure, the
T-coarsening outperforms the other two approaches in terms of iteration
counts. For the T-coarsening and MGRIT-based approach, the convergence rate
deteriorates as the Courant number $C$ increases, with the latter at a much
higher rate. The convergence of the TS-coarsening approach improves when $C$
increases towards 1 but deteriorates as $C$ increases further (not shown in
the figure). Overall, this convergence behavior is consistent with the
behavior we observed earlier in the 1-dimensional case; cf. Section 3.3.1.

Figure 8: Convergence performance of the sequential two-level Krylov on single
processor for solving the 2-dimensional heat equation with the T- and TS-
coarsening, and MGRIT-based approach for various values $C=3\tau\Delta
t/\Delta x^{2}$. Top-left figure: numbers of iterations; Top-right figure:
total time relative to the sequential $\theta$-scheme; Bottom figure: Coarse
time-grid solve (CGS) time relative to total time.

The T-coarsening strategy also results in the fastest method in terms of the
total computational time; see Figure 8:top-right. Note that, in this figure,
the time is presented as the time relative to the computational time of the
sequential $\theta$-scheme. All methods, however, require more computational
time than the $\theta$-scheme, with the T-coarsening two-level Krylov method
be slower by a factor of at least 2. Due to the inversion needed in the
interpolation step, the MGRIT-based two-level method becomes much slower than
the T-coarsening. With $\nu=2$, for instance, the MGRIT-based two-level method
has to perform $n_{t}/2$ inversions of a matrix of size of $n^{2}$, with
$n=n_{x}n_{y}$ in the interpolation step per iteration in addition to
$n_{t}/2$ inversion of a matrix of the same size in the coarse time-grid
solve. In constrast, the simple interpolation used in the T-coarsening does
not involve inversion and costs only of order $n=n_{x}n_{y}$ computational
work, which is negligible compared to the cost of inversion. Consequently, for
the T-coarsening method, the computational work is dominated (about 90%) by
the coarse time-grid solve (see Figure 8:bottom).

Computational time can be reduced via parallel implementation of the method.
Note that in the two-level approach, apart from the matrix inversion, the
other operations such as matrix-vector multiplications, restriction, and
interpolation can naturally be performed in parallel. Theoretically, taking
into account only the matrix-vector multiplication and the coarse time-grid
solve and assuming the inversion work of order of $n_{x}^{2}n_{y}^{2}$ for the
2-dimensional heat equation, the work a GMRES-based two-level method requires
to reach convergence is

$W_{\text{2-Level}}\sim
2n_{\text{iter}}^{(1)}n_{t}n_{x}n_{y}+n_{\text{iter}}^{(1)}\frac{n_{t}}{\nu}n_{x}^{2}n_{y}^{2},$

with $n_{\text{iter}}^{(1)}$ the number of GMRES iteration on the first
(finest) time level. Comparing with the $\theta$-scheme, whose total work is
given by

$W_{\text{$\theta$-scheme}}\sim n_{t}n_{x}n_{y}+n_{t}n_{x}^{2}n_{y}^{2},$

it is practically impossible to make a two-level multilevel Krylov faster than
the $\theta$-scheme without parallelization. With parallelization of matrix-
vector operations, the computational work on a processor is given by

$W_{\text{2-Level}}\sim\frac{2n_{\text{iter}}^{(1)}}{N_{\text{proc}}}n_{t}n_{x}n_{y}+n_{\text{iter}}^{(1)}\frac{n_{t}}{\nu}n_{x}^{2}n_{y}^{2},$

where $N_{\text{Proc}}$ is the number of processors. In this case, the
computational work on matrix-vector operations can be made much less than that
of the $\theta$-scheme using a large number of processors. The effect of
parallelization of the matrix/vector operations on the computational time is
shown in Figure 9. As seen in the figure, parallel implementation of the
T-coarsening approach reduces the computational time to a point that the
contribution of the matrix/vector operations to the overall computational time
becomes negligible. In the MGRIT-based approach, even though the interpolation
can also be performed in parallel, matrix inversions involved in the
interpolation process has a non-negligible contribution to the total
computational time.

Figure 9: Time spent on the coarse time-grid solve in the simulated parallel
setting. The underlying problem is the 2-dimensional heat equation,
discretized with $C=3\tau\Delta t/\Delta x^{2}=0.16$.

Increased parallelism in the coarse time-grid solve can be attained by zeroing
some off-diagonal parts of the coarse time-grid matrix $\mathcal{A}_{H}$,
hence decoupling the coarse time-grid system. For example, in the extreme
case, we may drop all off-diagonal blocks and use the two-level Krylov matrix

(33)
$\displaystyle\widetilde{\mathcal{Q}}_{h}=\mathcal{I}_{h}-\mathcal{Z}\widetilde{\mathcal{A}}_{H}^{-1}\mathcal{Y}^{\top}\mathcal{A}_{h}+\mu\mathcal{Z}\widetilde{\mathcal{A}}_{H}^{-1}\mathcal{Y}^{\top},$

where $\widetilde{\mathcal{A}}_{H}=\text{blockdiag}(\mathcal{A}^{H}_{k,k})$.
In this case, the coarse time-grid system is now completely decoupled and the
solution procedure becomes parallelizable. For a choice of $N_{\text{CGS}}$
off-diagonal blocks set to zero, the course-grid system is decoupled into
$N_{\text{CGS}}$ systems, which can be solved in parallel. In this case, with
the approximate coarse-grid matrx $\widetilde{\mathcal{A}}_{H}$, one can show
that

(34)
$\displaystyle\mathcal{A}_{h}\widetilde{\mathcal{Q}}_{h}=\begin{bmatrix}\mathcal{A}^{h}_{0,0}\mathcal{Q}^{h}_{0}&&&\\\
\star&\mathcal{A}^{h}_{1,1}\mathcal{Q}^{h}_{1}&&\\\ \star&\star&\ddots&\\\
&\star&\star&\mathcal{A}^{h}_{n_{T},n_{T}}\mathcal{Q}^{h}_{n_{T}}\end{bmatrix}.$

Hence,
$\text{blockdiag}(\mathcal{A}_{h}\widetilde{\mathcal{Q}}_{h})=\text{blockdiag}(\mathcal{A}_{h}\mathcal{Q}_{h})$,
and consequently,

(35)
$\displaystyle\sigma(\mathcal{A}_{h}\widetilde{\mathcal{Q}}_{h})=\sigma(\mathcal{A}_{h}\mathcal{Q}_{h}).$

Figure 10:left shows numbers of iterations to reach convergence using the
T-coarsening and the approximate coarse time-grid matrix
$\widetilde{\mathcal{A}}_{H}$. All matrix-vector operations are performed in
parallel using 32 simulated processors. The coarse time-grid system is
decoupled into $N_{\text{CGS}}$ systems by zeroing
$N_{\text{CGS}}=N_{\text{Proc}}$ off-diagonal blocks, such that each simulated
processor has an equal load, equivalent to inverting an
$\frac{n_{T}}{N_{\text{CGS}}}\times\frac{n_{T}}{N_{\text{CGS}}}$ block matrix.
As shown in Figure 10:left, the number of iterations increases linearly in the
number of processors, but is almost indendent of $C$. Thus, even though
$\mathcal{Q}_{h}$ and $\widetilde{Q}_{h}$ lead to spectrally equivalent two-
level Krylov matrices, this fact does not necessarily translate to an idential
convergence behavior, because eigenvalues are not the only factor that
determines the convergence. The parallel coarse-grid solve however reduces the
total computational time down significantly, even though there is a tendency
to increase again with a larger number of used processors. In particular, with
$n_{x}=n_{y}=192$, the two-level Krylov method outperforms the $\theta$-scheme
(Figure 10:mid). We note here that the total work with full parallelization of
the two-level method is

$W^{\text{full
parallel}}_{\text{2-level}}\sim\frac{2n_{\text{iter}}^{(1)}}{N_{\text{proc}}}n_{t}n_{x}n_{y}+\frac{n_{\text{iter}}^{(1)}}{\nu
N_{\text{CGS}}}n_{t}n_{x}^{2}n_{y}^{2}.$

With $N_{\text{CGS}}=N_{\text{proc}}$ and
$n^{(1)}_{\text{iter}}/N_{\text{CGS}}<1.5$ as observed in Figure 10:left, we
have the constant $n^{(1)}_{\text{iter}}/(\nu N_{\text{proc}})<1$, justifying
this observed speed-up. The propotion of work on the coarse time-grid solve
also decreases (Figure 10:right; compare this with Figure 9) when more
processors are used.

Figure 10: Effect of parallel coarse time-grid solves on the performance of
two-level Krylov method with the T-coarsening and approximate coarse time-grid
matrix $\widetilde{A}_{H}$. All matrix-vector operations are performed in
parallel using 32 simulated processors. The coarse time-grid system is
decoupled into $N_{\text{CGS}}=N_{\text{proc}}$ systems. Left figure:
iteration numbers; Middle figure: total time relative to the $\theta$-scheme;
Right-figure: Time spent on the coarse time-grid solve relative to the total
computational time.

##  4 Parallel-in-Time Multilevel Krylov and numerical results

In this section, we discuss the extension of the two-level Krylov method
discussed in Section 3.3 to the multilevel Krylov method. We note that the
multilevel Krylov method is flexible in the way the coarse-grid system is
solved. What is crucial in the multilevel Krylov method is the accuracy and
the effeciency of the coarse-grid solve. One can employ, e.g., MGRIT, to solve
the coarse-grid system. A standard multilevel Krylov method, however, is built
upon a recursive application of the two-level Krylov method. For elliptic
PDEs, as the spatial grid becomes coarser, the underlying coarse space-grid
matrix tends to become better conditioned. Hence, the number of iterations
needed to approximately solve the coarse-grid system on the coarser levels can
be reduced. As reported in [3, 4], a multilevel Krylov method needs only a few
GMRES iterations on the coarser levels.

###  4.1 Time coarsening

For the T-coarsening multilevel Krylov method with $L$ levels, it can be shown
the computational work is

(36) $\displaystyle W^{\text{full
parallel}}_{\text{$L$-level}}\sim\frac{2n_{\text{iter}}^{(1)}}{N_{\text{proc}}}\left(1+\sum_{\ell=2}^{L-1}\frac{1}{\nu^{\ell-1}}\prod_{j=2}^{\ell}n_{\text{iter}}^{(j)}\right)n_{t}n_{x}n_{y}+\frac{n_{\text{iter}}^{(1)}}{N_{\text{CGS}}}\left(\prod_{\ell=2}^{L-1}\frac{n_{\text{iter}}^{(\ell)}}{\nu^{L-2}}\right)n_{t}n_{x}^{2}n_{y}^{2}.$

The above estimate suggests that the multilevel implementation adds extra cost
in matrix/vector operations. The cost of coarse time-grid solve can however
potentially be reduced if the number of iterations used on the coarse levels
is less than $\nu$. For standard coarsening, this will require
$n_{\text{iter}}^{(\ell)}\leq\nu=2$. As observed in the experiments, this
requirement is quite stringent if applied on every level, and results in slow
convergence due to inaccuracy in the coarse-grid solve. Because the accuracy
of the second-level coarse-grid solve seems to be crucial, we shall allow more
GMRES iterations on the second level and compensate the extra cost they bring
by performing only one or two iterations on the coarser levels.

Figure 11 shows performance of 2-, 3-, 4-level Krylov method with T-coarsening
and with standard coarsening ($\nu=2$) applied to generate a sequence of
coarse time-grid systems. Matrix-vector operations at all levels are performed
on 32 simulated processors. On the coarsest level, the coarse time-grid system
is solved in parallel using the perturbed coarse-grid matrix
$\widetilde{A}_{H}$, with $N_{\text{CGS}}=N_{\text{proc}}$. At any level
between the finest and the coarsest level, 2 GMRES iterations are used.

As shown in Figure 11:left, with T-coarsening, the number of iterations tends
to increase as more coarse levels are used. This in turn increases the total
CPU time needed to reach convergence; see Figure 11:mid. In this case, the
parallel multilevel approach does not lead to a faster method than the
sequential $\theta$-scheme. Even though more iteration numbers have to be
spent when more processors are used on the coarsest level, the smaller coarse
time-grid system each processor has to solve results in reduced CPU times. As
shown in Figure 11:right, the coarse time-grid solve remains the most time
consuming operation, taking more than 90% of the total time.

Figure 11: Performance comparison of multilevel Krylov with time coarsening
and various numbers of processors used to solve the coarse-grid system in
parallel, with $n_{x}=n_{y}=128$, $n_{t}=768$. Other than the coarse grid
solve, 32 processors are used to perform all related matrix/vector operations.

We remark here that for time-dependent problems, coarser time grid may cause
instability, leading to highly inaccurate solve of the coarse time-grid
system. This inaccurate coarse time-grid solve can probably be compensated by
allowing more GMRES iterations on the fine levels. This however may in the end
render the method inefficient.

###  4.2 Alternate time-space coarsening

Based on the estimate (36), the computational time can be reduced further via
the inclusion of spatial coarsening. For the 2d case, the use of standard
coarsening in space (roughly doubling the spatial mesh size) at a level
reduces the first and the second term in the estimate by a factor of 2 and 4,
respectively. One strategy for spatial coarsening is by performing time-space
coarsening simultaneously. However, as suggested by the test results in
Section 3.3.1, this approach leads to a slow-to-converge two-level Krylov
method. As the coarse-grid solve needs to be done as accurately as possible,
applying a time-space coarsening at a particular level will require a
significant number of iterations on the level the simultaneous time-space
coarsening is performed.

An alternative to the simultaneous time-space coarsening is by applying
spatial coarsening on some coarse levels whenever the time grid has been
sufficiently coarsened. This strategy is illustrated in Figure 12, which
involves 5 coarsening. Starting from the fine time-space grid ($\ell=1$), the
next three coarse-level systems are constructed via time coarsening, indicated
by the black squares. After this sequence of time coarsening, the remaining
coarse-level systems are constructed via space and time coarsening, applied
alternatingly (the space coarsening is indicated by the white square). This
coarsening strategy results in a coarse-level system associated with $n_{t}/8$
time steps $(\Delta T=8\Delta t)$. In each time step, a coarse matrix of size
$n_{x}n_{y}/4\times n_{x}n_{y}/4$ has to be inverted. For the numerical
results presented here, the first space coarsening is performed whenever the
sequence of time coarsening in the finer level results in time-space grids
which violates the stability condition.

Figure 12: An example of time-space coarsening, resulting in 6 levels.

We show computational performance of the multilevel Krylov method with the
time-space coarsening outlined above in Figure 13. The computations were
performed for various Courant numbers $C$, with two spatial gridpoint
settings: $n_{x}=n_{y}=128$ and $192$. In this case, with $C=0.16$ for
example, the initial time coarsening is performed until the 4th coarse level,
followed by alternating space-time coarsening. On each level, computations are
done in parallel using a fixed number of simulated processors, except on the
coarsest level, where in all cases, the coarsest-grid system is decoupled into
16 systems, solved in parallel on 16 simulated processors. On the coarse
levels, 2 FGMRES iterations are performed, except on the level before the
coarsest one, where 1 FGMRES iteration is used.

As shown in Figure 13:left, iteration numbers to converge depend on the
Courant number $C$, with lower $C$ resulting in a faster convergence. The
results also suggests that the number of iterations to converge is mildly
dependent on the grid settings, so long as $C$ is kept constant. The total CPU
time relative to the sequential $\theta$-scheme is shown in Figure 13:mid. The
graph demonstrates the gain in CPU time when a large number of processors is
used (in this case, $n_{\text{proc}}>8$), outperforming the serial
$\theta$-scheme for all $C$ used. Furthemore, the computational time on the
coarse-grid solve is reduced to below 90% of the total computational time.

Figure 13: Performance comparison of multilevel Krylov with time-space
coarsening and various numbers of processors used to solve the coarse-grid
system in parallel. On the coarsest level, 16 processors is used to solve the
coarse-grid system.

Finally, to investigate the effect of number of levels used in the method, we
consider a case with $n_{t}=2^{11}$ and $n_{x}=n_{y}=192$ on the finest level,
associated with $C=0.16$. The time-space grid is coarsened in the similar way
as in the previous test, with the coarsest level system solved in a decoupled
way on 16 simulated processors. Apart from the coarsest-grid solves,
matrix/vector operations are performed in parallel using either $32$ or $64$
simulated processors. The convergence results are shown in Table 3. We note
here that for the case with 32 processors and 7 levels, the data in the table
corresponds to the right-end of the curve in Figure 13, with $n_{x}=n_{y}=192$
and $C=0.16$. Thus, Table 3 suggests that the CPU time relative to the
$\theta$-scheme can be further reduced via the use of more levels to decrease
the size of the coarsest-grid system, without practically affecting the number
of iterations. In this case, about 90% of the computational time is spent on
the coarse-grid solves. With 11 levels, the method achieves a computational
speed-up of about 9.39 from the single-processor implementation.

Using 64 simulated processors on the matrix/vector operations reduces the
total CPU time further, as in this case, the computational time on the
coarsest-grid solves increases to around 95%; see Table 3. Under this setting,
we however stop the coarsening when reaching 64 time grids and after
performing the necessary spatial coarsening (on Level 9). It is indeed
possible that after Level 9, the coarsening is continued and the resulting
coarser-grid systems are handled on, for instance, 32 simulated processors.
This implementation scenario is, however, not considered here. Even with 9
levels, the method gains a computational speed up of about 9.85 relative to
the single-processor implementation.

Table 3: Performance comparison using different number of levels, using an alternate TS-coarsening, $n_{x}=n_{y}=192$, $n_{t}=2^{11}$. On the the coarsest level, 16 processors are used. Other than the coarsest level, all matrix/vectors operations are performed on 32 or 64 processors. The symbol “–” means coarser levels that can be fit in the assigned number of processors could not be constructed further. | coarsest level | No. of Processors
---|---|---
| time-space | 32 | 64
Level | grid | #Iter. | $T_{\text{CPU}}$ | $T_{\text{CGS}}$ | #Iter. | $T_{\text{CPU}}$ | $T_{\text{CGS}}$
7 | $128\times 48\times 48$ | 9 | 0.436 | 0.927 | 10 | 0.322 | 0.948
8 | $64\times 48\times 48$ | 10 | 0.428 | 0.913 | 10 | 0.365 | 0.947
9 | $64\times 24\times 24$ | 10 | 0.416 | 0.918 | 10 | 0.365 | 0.951
10 | $32\times 24\times 24$ | 10 | 0.422 | 0.916 | – | – | –
11 | $32\times 12\times 12$ | 9 | 0.384 | 0.916 | – | – | –

##  5 Conclusion

This paper discussed and extended the multilevel Krylov method for stationary
problems to the time-dependent problems, which allows parallel-in-time
computational strategies. Likewise parareal- and MGRIT-based methods, the
method is iterative in nature and involves a coarse-level problem on which
parallelization can be exploited in a multilevel fashion. As the coarsest
level is typically the bottleneck for full time parallelization, we proposed
an approximation to the coarsest-grid system by dropping some off-diagonal
blocks in the coarsest-grid matrix, effectively decoupling the systems and
allowing full time parallelization on all levels. This approximation was
observed to be quite effective, especially when a sufficient number of coarse
time levels was used.

Numerical tests suggested that a significant gain in computational time could
not be achieved solely with time coarsening. The proposed method incorporated
the space coarsening in an alternating way, following a sequence of first few
time coarsening until a stability condition (dictated by the Courant number)
is violated. This approached led to the total computational time as low as
0.32 relative to the sequential $\theta$-scheme.

The current implementation was limited by the author’s computational
resources. It is therefore expected to perform an actual test of the method in
the future on a parallel computer. Furthermore, as this paper was meant to
present only the prototypical implementation of the method, future’s work will
include, among others, the effective construction of the deflation
(prolongation and restriction) matrices, even though the current simple
construction seemed to work very well, time-space coarsening and as well the
multilevel strategies, which would lead to further gain in computational time.

## References

  * [1] V. Dobrev, T. Kolev, N. Petersson, and J. Schroder, Two-level convergence theory for multigrid reduction in time (mgrit), SIAM Journal on Scientific Computing, 39 (2017), pp. S501–S527.
  * [2] Y. A. Erlangga and R. Nabben, Deflation and balancing preconditioners for krylov subspace methods applied to nonsymmetric matrices, SIAM Journal on Matrix Analysis and Applications, 30 (2008), pp. 684–699.
  * [3] Y. A. Erlangga and R. Nabben, Multilevel projection-based nested Krylov iteration for boundary value problems, SIAM Journal on Scientific Computing, 30 (2008), pp. 1572–1595.
  * [4] Y. A. Erlangga and R. Nabben, Algebraic mutilevel krylov methods, SIAM Journal on Scientific Computing, 31 (2009), pp. 3417–3437.
  * [5] R. D. Falgout, S. Friedhoff, T. V. Kolev, S. P. MacLachlan, and J. B. Schroder, Parallel time integration with multigrid, SIAM Journal on Scientific Computing, 36 (2014), pp. C635–C661.
  * [6] R. D. Falgout, M. Lecouvez, and C. S. Woodward, Parallel time integration with multigrid, SIAM Journal on Scientific Computing, 37 (2019), p. 101029.
  * [7] R. D. Falgout, T. A. Manteuffel, B. O’Neill, and J. B. Schroder, Multigrid reduction intime for nonlinear parabolic problems: a case study, SIAM Journal on Scientific Computing, 39 (2017), pp. S298–S322.
  * [8] J. Frank and C. Vuik, On the construction of deflation-based preconditioners, SIAM Journal on Scientific Computing, 23 (2001), pp. 442–462.
  * [9] M. J. Gander, Overlapping schwarz for linear and nonlinear parabolic systems, in Proceedings of the 9th International Conference on Domain Decomposition, 1996, pp. 97–104.
  * [10] M. J. Gander, 50 years of time parallel time integration, in Multiple Shooting and Time Domain Decomposition Methods, T. Carraro, M. Geiger, S. Körkel, and R. Rannache, eds., Cham, 2015, Springer International Publishing, pp. 69–113.
  * [11] M. J. Gander, L. Halpern, and F. Nataf, Optimal schwarz waveform relaxation for the one dimensional wave equation, SIAM Journal on Numerical Analysis, 41 (2003), pp. 1643–1681.
  * [12] M. J. Gander and S. Vandewalle, Analysis of the parareal time-parallel timeintegration method, SIAM Journal on Scientific Computing, 29 (2007), pp. 556–578.
  * [13] W. Hackbusch, Parabolic multigrid methods, in Computing Methods in Applied Sciences and Engineering VI, R. Glowinski and J.-L. Lions, eds., North-Holland, 1984, pp. 189–197.
  * [14] T.-P. Hoang, J. Jaffré, C. Japhet, M. Kern, and J. E. Roberts, Space-time domain decomposition methods for diffusion problems in mixed formulation, SIAM Journal on Numerical Analysis, 51 (2013), pp. 3532–3559.
  * [15] G. Horton and S. Vandewalle, A space-time multigrid method for parabolic partial differential equations, SIAM Journal on Scientific Computing, 16 (1995), pp. 848–864.
  * [16] A. J. Howse, H. de Sterck, R. D. Falgout, S. Maclachlan, and J. Schroder, Parallel-in-time multigrid with adaptive spatial coarsening for the linear advection and inviscid burger’s equation, SIAM Journal on Scientific Computing, 41 (2019), pp. A538–A565.
  * [17] J.-L. Lions, Y. Maday, and G. Turinici, A parareal in time discretization of PDEs, Acad. Sci. Paris, Serie I, 332 (2021), pp. 661–668.
  * [18] Y. Maday and G. Turinici, The parareal in time iterative solver: a further direction to parallel implementation, in Domain Decomposition Methods in Science and Engineering, Springer, 2005, pp. 441–448.
  * [19] R. Nabben and C. Vuik, A comparison of deflation and coarse grid correction applied to porous media flow, SIAM Journal on Numerical Analysis, 42 (2004), pp. 1631–1647.
  * [20] M. Neümuller, Space-Time Methods: Fast sover and Applications, PhD thesis, University of Graz, 2014.
  * [21] R. A. Nicolaides, Deflation of conjugate gradients with applications to boundary value problems, SIAM J. Numer. Anal., 24 (1987), pp. 355–365.
  * [22] J. Nievergelt, Parallel methods for intergrating ordinary differential equations, Communications of the ACM, 7 (1964), pp. 731–733.
  * [23] Y. Saad, A flexible inner-outer preconditioned gmres algorithm, SIAM Journal on Scientific Computing, 14 (1993), pp. 461–469.
  * [24] J. M. Tang, R. Nabben, C. Vuik, and Y. A. Erlangga, Comparison of two-level preconditioners derived from deflation, domain decomposition and multigrid methods., Journal of Scientific Computing, 39 (2009), pp. 340–370.