\documentclass [12pt] {llncs}
\usepackage {a4,amsmath,palatino}
\pagestyle{plain}
\title {The Bubka Project}
\subtitle {\small Towards the ultimate AMG parallel preconditioner\\
Part I: A massively parallel IBLU preconditioner}
\author {Jos van Kan \institute {Department of Technical Applied Analysis, Delft
University of Technology, Delft, Netherlands}}
\newcommand \calE {{\cal E}}
\newcommand \IZ {\bbbz}
\newcommand \eop {\hfill$\Box$\par\medskip\noindent}
%\newtheorem {lemma}{Lemma}[section]
%\newtheorem {theorem}{Theorem}[section]
\begin {document}
\maketitle
\section {Introduction}
The Bubka project aims at constructing the ultimate parallel preconditioner with
\begin {enumerate}
\item Multigrid efficiency
\item Robustness against anisotropy
\item Typically only O(${}^2\log N$) recursive steps needed.  On a massively parallel
computer this would imply an amount of computational cycles of the same order 
of magnitude.
\end {enumerate}
It will not achieve this aim in one fell swoop, but rather like the famous
Russian athlete that lends his name to the project proceed in small steps,
everyone of which will achieve a subgoal.

These subgoals are
\begin {enumerate}
\item Parallellization of the incomplete block LU-decomposition. This in itself is
a worthwile aim: as it things stand now IBLU decomposition is completely
recursive, both globally (blocks in the matrix) and locally (within a single 
block). Since our new algorithm comes down to full cyclic reduction, we 'll 
call it: CRIBLU.
\item Alternating directions in CRIBLU. As we are about to discover, IBLU and 
also CRIBLU have
a preferred direction, that is not mirrored in the problem itself. The idea is,
that alternating directions in the algorithm may make the preconditioner more 
robust.
\item CRIBLU operates on a hierarchy of subsets that decrease in size
exponentially. The idea of exploiting these subsets in a multigrid like way is a
natural one. Viewed like this, CRIBLU solves a set of approximate problems on
a specific level, by using an incomplete block decomposition on the lower level.
If you'd use the solution in one direction as the `coarse grid correction' and
the solution in the other direction (but of the residual) as a `smoother' (or
vice versa) you'd have a perfectly symmetric AMG algorithm. Big question: does
this AMG algorithm have MG efficiency?
\end {enumerate}

\section {\label {sec:FCR}Full Cyclic Reduction for Tridiagonal Matrices}
We'll start out, by taking a different view on cyclic reduction. Consider a
tridiagonal matrix $A^0 = D^0-L^0-U^0$ in which $D^0$ is a a diagonal, $L^0$
a subdiagonal
and $U^0$ a superdiagonal matrix. Now consider the following transformation:  
\newcommand \dinv [1] {(D^{#1})^{-1}}
\begin {equation}
\begin {split} \label {eq:new_a1}
A^1 &= (D^0+L^0+U^0) \dinv0 A \\
    &= D^0 - L^0\dinv0 U^0 - U^0\dinv0 L^0 - L^0\dinv0 L^0
 - U^0\dinv0 U^0\\
    &= D^1 - L^1 - U^1\\
\end {split}
\end {equation}
where
\begin {align} \label {eq:new_d1}
D^1 &= D^0 - L^0\dinv0 U^0 - U^0 \dinv0 L^0\\
\label {eq:new_l1}
L^1 &= L^0 \dinv0 L^0\\
\label {eq:new_u1}
U^1 &= U^0 \dinv0 U^0
\end {align}
In these formulae $D^1$ is again a diagonal matrix, and $L^1$ and $U^1$ are sub and
superdiagonals respectively, {\em once removed}; that is, of the form $l^1_{i,
i-2}$ and $u^1_{i, i+2}$ respectively.
We can continue in the same vein and define recursively:
\begin {equation}
\begin {split} \label {eq:new_ak}
A^{k+1} &= (D^{k}+L^{k}+U^{k}) \dinv{k} A \\
&= D^{k} - L^{k}\dinv{k} U^{k} - U^{k}
        \dinv{k} L^{k} - L^{k}\dinv{k} L^{k} - U^{k}\dinv{k} U^{k}\\
    &= D^{k+1} - L^{k+1} - U^{k+1}\\
\end {split}
\end {equation}

where
\begin {align} 
\label {eq:new_dk}
D^{k+1} &= D^{k} - L^{k}\dinv{k} U^{k} - U^{k} \dinv{k} L^{k}\\
\label {eq:new_lk} 
L^{k+1} &= L^{k} \dinv{k} L^{k}\\
\label {eq:new_uk} 
U^{k+1} &= U^{k} \dinv{k} U^{k}
\end {align}

The reader will have no trouble to verify, that $L^k$ and $U^k$ are of the form
$l^k_{i, i-2^k}$ and $u^k_{i, i+2^k}$ respectively. As soon as $2^k$ exceeds
$N$,
the size of the matrix, $L^k$ and $U^k$ will be empty and hence it will take 
${}^2\log N$ steps to
reduce $A^0$
in this way to a diagonal matrix. This would make
little sense for a serial computer, or a parallel computer with few processors.
But all matrix multiplications are in principle parallel operations without
recursion, so on a computer with $N$ or more processors the  number of  cycles 
would effectively reduce to $O({}^2\log N).$

\section {Incomplete Block LU Factorization}
Various authors \cite{Axe84BI}, \cite{Meu84}, \cite{Con85GM},\cite{Kett_87} have considered incomplete block LU
factorization, called IBLU (incomplete block LU), both as a preconditioner for a
Krylov space method or as a smoother for a multi grid method. 
The principle is
as follows. Let $A = D - L - U$ be a {\em block} tridiagonal matrix. We write
(following an idea of Eisenstadt's) 
\newcommand \td {\tilde D}
\begin {equation}
A = (\td - L) \td^{-1} (\td - U)
\end {equation}
and we choose $\td$ such, that $D = \td + L \td^{-1} U$. It will be readily
seen, that the latter identity can be obtained by block recursion. The first
block of $D$, $D_1$ and of $\td$, $\td_1$ are identical and we have the recursion
\begin {equation} \label {eq:exact_recursion}
\td_{j+1} = D_{j+1} + L_{j+1,j} \td_j^{-1} U_{j,j+1}
\end {equation}
So far this is an exact relation. If the blocks of $D$ are themselves
tridiagonal, the inverses will be full matrices. An incomplete version is
obtained, by retaining only {\em three} diagonals in every $\td_{j}$. Define the
operator $\calE_m$ as a {\em diagonal extractor}. $\calE_0$ extracts the main
diagonal, $\calE_1$ three diagonals (sub, main and super), $\calE_2$ five,
etcetera. With this notation the incomplete version of \ref{eq:exact_recursion}
can be constructed as follows:
\begin {equation} \label {eq:incomplete_recursion}
\td_{j+1} = D_{j+1} + {\calE_1} (L_{j+1,j} \td_j^{-1} U_{j,j+1})
\end {equation}
The reader will have no trouble in verifying, that this is equivalent to  
\begin {equation} \label {eq:incomplete_recursion2}
\td_{j+1} = D_{j+1} + \calE_1 (L_{j+1,j} \calE_3 (\td_j^{-1}) U_{j,j+1}))
\end {equation}
since $L_{j+1,j}$ and $U_{j,j+1}$ are tridiagonal.
So in order to calculate this expression one needs the seven diagonals of
$\td_j^{-1}$. It is explained a.o. in \cite{Kett_87} how to calculate those
diagonals.
This algorithm, however, is strongly recursive and unsuitable for parallel 
computation.

\section {Full Cyclic Block Reduction}
To find a parallel IBLU, let us reinterpret the formulae of 
section~\ref{sec:FCR} in the following way.
Let $A$ be a {\em block} tridiagonal matrix. Let the diagonal $D^0$ of $A$ 
consist
of $N_x$ blocks of size $N_y \times N_y$, $D^0_j, j=1\dots N_x$, the
superdiagonal $U^0$ of blocks $U^0_{j,j+1}, j=1\dots N_x-1$ and the
subdiagonal $L^0$ of blocks $L^0_{j,j-1}, j=2\dots N_x$. Now we may perform the exact
same transformations as in section~\ref{sec:FCR} to obtain a full block cyclic
reduction. In fact, there is a complete block analogue to
equations~\ref{eq:new_a1}-\ref{eq:new_uk}, using the same formulae but with the
connotation that it concerns block matrices.

However, equations~\ref{eq:new_dk}-\ref{eq:new_uk} may also be written
down blockwise. Equation~\ref{eq:new_dk} then becomes:
\begin {alignat} {2} \label {eq:new_blockdk}
D^{k+1}_j &= D^k - U^k_{j,j+2^k}(D^k_{j+2^k})^{-1}L^k_{j+2^k,j}, &\quad
j=1,\dots,2^k\\
D^{k+1}_j &= D^k - L^k_{j,j-2^k}(D^k_{j-2^k})^{-1} U^k_{j-2^k,j}+ \\\nonumber
          &\qquad- U^k_{j,j+2^k}(D^k_{j+2^k})^{-1}L^k_{j+2^k,j}, &\quad j=2^k + 1,
	  \dots, N_x - 2^k\\
D^{k+1}_j &= D^k - L^k_{j,j-2^k}(D^k_{j-2^k})^{-1} U^k_{j-2^k,j}, &\quad j = N_x -
2^k + 1, \dots, N_x
\end {alignat}
Equation~\ref{eq:new_lk} becomes
\begin {equation} \label {eq:new_blocklk}
L^{k+1}_{j+2^{k+1},j} = L^k_{j+2^{k+1}, j+2^k} (D^k_{j+2^k})^{-1}L^k_{j+2^k,
j}, j = 1, \dots, N_x - 2^{k+1}
\end {equation}
and finally equation~\ref{eq:new_uk} becomes
\begin {equation} \label {eq:new_blockuk}
U^{k+1}_{j, j+2^{k+1}} = U^{k}_{j, j+2^k} (D^k_{j+2^k})^{-1}
U^{k}_{j+2^k,j+2^{k+1}}
\end {equation}

Note, that there is no recursion on j in these expressions. For any given $k$
all calculations can be done in parallell. The only recursion left is that on
$k$ and for the same reason as in section~\ref{sec:FCR} this will take 
${}^2\log N_x$ cycles on a massive parallell machine. 

There are a couple of ifs and buts about this. You have to calculate the
inverses $(D^k_j)^{-1}$, which are size $N_y \times N_y$. As we have shown in
section~\ref{sec:FCR} this will take ${}^2\log N_y$ cycles on a massively
parallel machine. So the whole operation will take $O(\log N_x + \log N_y)$
cycles. The memory requirements will be prohibitive though: 
$3 \times N_x \times N_y^2 \times {}^2\log N_x$ to store the intermediate
matrices, because starting from $A^1$ the blocks will be full. 

\section {Incomplete Cyclic Block Reduction}

To reduce the memory requirements one has to do an incomplete cyclic block 
reduction.
Instead of using equations~\ref{eq:new_blockdk}-\ref{eq:new_blockuk} we only retain
the main, super and subdiagonal in every nonzero block. This gives the following
set of recursions:
For the diagonal blocks:
\begin {alignat} {2} \label {eq:new_incblockdk}
D^{k+1}_j &= D^k - \calE_1(U^k_{j,j+2^k}\calE_3({D^k_{j+2^k}}^{-1})L^k_{j+2^k,j}), &\quad
j=1,\dots,2^k\\
D^{k+1}_j &= D^k - \calE_1(L^k_{j,j-2^k}\calE_3({D^k_{j-2^k}}^{-1}) U^k_{j-2^k,j}) +
&\\\nonumber
&\qquad-\calE_1(U^k_{j,j+2^k}\calE_3({D^k_{j+2^k}}^{-1}) L^k_{j+2^k,j}), &\quad j=2^k + 1,
\dots, N_x - 2^k\\
D^{k+1}_j &= D^k - \calE_1(L^k_{j,j-2^k}\calE_3({D^k_{j-2^k}}^{-1}) U^k_{j-2^k,j}), &\quad j = N_x -
2^k + 1,\dots, N_x
\end {alignat}
for the subdiagonal blocks
\begin {equation} \label {eq:new_incblocklk}
L^{k+1}_{j+2^{k+1},j} = \calE_1(L^k_{j+2^{k+1}, j+2^k} \calE_3
({D^k_{j+2^k}}^{-1}) L^k_{j+2^k,
j}), j = 1, \dots, N_x - 2^{k+1}
\end {equation}
and for the super diagonal blocks
\begin {equation} \label {eq:new_incblockuk}
U^{k+1}_{j, j+2^{k+1}} = \calE_1(U^{k}_{j, j+2^k} \calE_3({D^k_{j+2^k}}^{-1})
U^{k}_{j+2^k,j+2^{k+1}}), \quad j=1,\dots,N_x - 2^{k+1}
\end {equation}

It is easily verified, that the structure of all nonzero blocks remains
tridiagonal.

\section {The Magnificent Seven}
The problem remains, however, how to calculate $\calE_3 (A^{-1})$ if $A$ is
tridiagonal. We recall from section~\ref{sec:FCR} how the inverse was
constructed: 
\begin {equation}
A^{-1}= \dinv {k+1} E^k E^{k-1} \dots E^1 E^0
\end {equation}
with $E^j = (L^j + D^j + U^j)
(D^j)^{-1}$ and $2^{k-1} < N \leq 2^k$. The typical structure of $E^j$ is
tridiagonal, with a main diagonal and a super and subdiagonal $2^j$ places
removed. Let us denote the pattern of such a matix by $d_{2^j} + d_0 + d_{-2^j}$. 
The product of matrices of such a structure has a specific property,
that will be the subject of the next lemma.
\begin {lemma} \label {lem:pattern}
The pattern of nonzero diagonals of the product 
$P^j = E^k\cdot E^{k-1}\cdot\dots\cdot E^{k-j}$ is given by
$d_{m2^{k-j}}, m \in \IZ, -N < m2^{k-j} < N$. All other diagonals are 0.
\end {lemma}

%\noindent{\bf Proof}
\begin {proof}
The lemma holds apparently for $j=0$. Suppose it holds up to $j$. Since
$E^{k-j-1}$ has pattern $d_{2^{k-j-1}} + d_0 + d_{-2^{k-j-1}}$ the resulting
matrix has a diagonal pattern that is a result of multiplying any existing
diagonal by one of the three diagonals of $E^{k-j-1}$. One easily verifies, that
multiplying a matrix consisting of a single diagonal $d_n$ by a matrix
consisting of a single diagonal $d_m$ results in a matrix consisting of a single
diagonal $d_{n+m}$. Now let $d_{m2^{k-j}}, m \in \IZ$ be any diagonal of $P^j$.
Multiplication by $d_0$ results in the same pattern:
$d_{m2^{k-j}}=d_{(2m)2^{k-j-1}}$.
Multiplication by $d_{2^{k-j-1}}$ results in a diagonal $d_{(2m+1)2^{k-j-1}}$ and
multiplication by $d_{-2^{k-j-1}}$ in a diagonal $d_{(2m-1)2^{k-j-1}}.$\qed
\end {proof}

The Lemma is instrumental in the
following theorem that enables us to calculate $\calE_3 (A^{-1})$ in a parallel
way.
\begin {theorem}
To calculate 7 central nonzero diagonals of $P^k$, it is sufficient to 
retain 3 central nonzero diagonals of $P^j, j=1, k-2$ and 5 central nonzero
diagonals of $P^{k-1}$.   
\end {theorem}

%\noindent{\bf Proof} 
\begin {proof} Following the reasoning of Lemma~\ref{lem:pattern} it will
be clear, that to calculate the seven diagonals $d_{-3},d_{-2},\dots, d_2,d_3$
of $P^k$ one needs $d_{-4}, d_{-2}, d_0, d_2, d_4$ of $P^{k-1}$ and
$d_{-1},d_0,d_1$ of $E^0$. 
The five diagonals of $P^{k-1}$ 
can be calculated from $d_{-4}, d_0, d_4$ of $P^{k-2}$ and $d_{-2},d_0,d_2$ of 
$E^1$. And the 3 central diagonals $d_{-2^j},d_0,d_{2^j}$ of $P^{k-j}$ can be
calculated from the 3 central diagonals $d_{-2^{j+1}},d_0,d_{2^{j+1}}$ of
$P^{k-j-1}$ and $d_{-2^j},d_0,d_{2^j}$ of $E^j$.\qed 
\end {proof}
 
\noindent{\bf Remark} It will be evident, that again there is only a ${}^2\log N$
cycle recursion in this calculation and that the rest of the calculations can be
executed in parallel.

\bibliography {BubkaBib,biblioT}
\bibliographystyle {plain}
\end {document}
