A declared two-level MGRIT coarse-correction cycle

Input: a linear fine-grid system $A_{t}u=b$ with $N\ge 1$ fine intervals as in The all-at-once bidiagonal time operator and C/F points, a coarsening factor $m\ge 2$ dividing $N$ and its C/F partition as in Eliminating F-points gives the exact coarse time operator, practical coarse propagators ${\mathrm{\Psi }}_j\in {\mathbb{R}}^{d\times d}$ for $1\le j\le N_c:=N\mathrm{/}m$, an initial fine-grid iterate $u_{\mathrm{i}\mathrm{n}}$, pre- and post-relaxation choices from $\{F,FCF\}$, a real-valued norm on the fine-grid residual space as in A norm on a real vector space, the induced metric, and the dictionary with the metric axioms, a nonnegative integer cycle budget $J_{\mathrm{c}\mathrm{y}\mathrm{c}}$, and a residual tolerance $\epsilon \ge 0$. The symbol $J_{\mathrm{c}\mathrm{y}\mathrm{c}}$ is the number of permitted cycles; the coarse-interval count is $N_c$ (the letter $J$ in Eliminating F-points gives the exact coarse time operator denotes that lemma's coarse-interval count). The output is the last iterate, its residual norm, the number of completed cycles, and status converged or budget-exhausted.

In F/C ordering define ideal interpolation $P$ and restriction $R$ by $P=\left[\begin{array}{c}-A_{FF}^{-1}{A}_{FC}\\ I_C\end{array}\right],R=\left[\begin{array}{cc}-A_{CF}{A}_{FF}^{-1}& I_C\end{array}\right]\mathrm{.}$

The exact coarse operator is the Schur complement $R{A}_{t}P=A_{CC}-A_{CF}{A}_{FF}^{-1}{A}_{FC}$ from Eliminating F-points gives the exact coarse time operator. The practical correction matrix $A_c$ has identity diagonal blocks and subdiagonal blocks $-{\mathrm{\Psi }}_j$ for $2\le j\le N_c$; its first row corresponds to the fixed zero initial correction $z_0=0$. For any right-hand side $\rho$, forward substitution uniquely solves $A_{c}z=\rho$ by setting $z_0=0$ and, in increasing order for $j=1,\dots ,N_c$, setting $z_j={\rho }_j+{\mathrm{\Psi }}_j{z}_{j-1}$. Thus $A_c$ is invertible even when an individual ${\mathrm{\Psi }}_j$ is singular.

Define the relaxations on $A_{t}u=b$ as follows. An F-relaxation holds every C-value fixed and, for each coarse interval $j=1,\dots ,N_c$, updates its F-indices $i=(j-1)m+1,\dots ,jm-1$ in increasing order: set $u_1=b_1$ if $i=1$, and otherwise set $u_i=b_i+{\mathrm{\Phi }}_i{u}_{i-1}$. Different coarse intervals have independent F-updates and may be processed in parallel. A C-relaxation holds every F-value fixed and, for $j=1,\dots ,N_c$ in increasing order, sets $u_{jm}=b_{jm}+{\mathrm{\Phi }}_{jm}{u}_{jm-1}$; these C-updates are independent because $m\ge 2$ and each predecessor is an F-index. An FCF-relaxation applies F-relaxation, then C-relaxation, then F-relaxation. Each relaxation is a finite linear update; write $S_{\mathrm{p}\mathrm{r}\mathrm{e}}$ and $S_{\mathrm{p}\mathrm{o}\mathrm{s}\mathrm{t}}$ for the error maps of the selected pre- and post-relaxations.

For one cycle from a current iterate $u$, apply the selected pre-relaxation; form $r=b-A_{t}u$ and $\rho =Rr$; solve $A_{c}z=\rho$ by the forward substitution above; correct $u\leftarrow u+Pz$; then apply the selected post-relaxation. The exact error-propagation matrix of this fixed cycle is $E=S_{\mathrm{p}\mathrm{o}\mathrm{s}\mathrm{t}}(I-P{A}_c^{-1}R{A}_t)S_{\mathrm{p}\mathrm{r}\mathrm{e}}\mathrm{.}$

Initialize $u^{(0)}=u_{\mathrm{i}\mathrm{n}}$ and compute $\mathrm{\parallel }b-A_t{u}^{(0)}\mathrm{\parallel }$. If this is at most $\epsilon$, return $u^{(0)}$, zero completed cycles, and status converged. Otherwise, for $j_{\mathrm{c}\mathrm{y}\mathrm{c}}=1,\dots ,J_{\mathrm{c}\mathrm{y}\mathrm{c}}$, perform one cycle, call the resulting iterate $u^{(j_{\mathrm{c}\mathrm{y}\mathrm{c}})}$, and compute $\mathrm{\parallel }b-A_t{u}^{(j_{\mathrm{c}\mathrm{y}\mathrm{c}})}\mathrm{\parallel }$. If this norm is at most $\epsilon$, return that iterate, $j_{\mathrm{c}\mathrm{y}\mathrm{c}}$ completed cycles, and status converged; if not, continue. If all $J_{\mathrm{c}\mathrm{y}\mathrm{c}}$ cycles finish without this condition (including when $J_{\mathrm{c}\mathrm{y}\mathrm{c}}=0$), return the last iterate, its residual norm, and status budget-exhausted. The matrix $E$ describes exact arithmetic and exact coarse solves; inexact solves contribute an update defect that needs its own bound.