\documentclass [12pt] {llncs}
\usepackage {a4,amsmath,palatino}
\pagestyle{plain}
\title {A Parallel Incomplete Block LU Factorization}
\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}
\begin {document}
\maketitle
\begin {abstract}
A parallel algorithm is presented to construct an incomplete block factorization
for an $N \times N$ block tridiagonal matrix in which the blocks themselves are also tridiagonal matrices
of size $N_x \times N_x$. $N=N_x\times N_y.$ The
algorithm uses $O({}^2\log N_x \times {}^2\log N_y)$ recursive cycles, whereas 
a fully recursive algorithm would require $O(N)$ recursive cycles.

The factorization may be used as a preconditioner in a Krylov space method or as smoother in a
multigrid method.
\end {abstract}
\section {Introduction}
Let $A$ be an $N \times N$ tridiagonal block matrix consisting of blocks of size 
$N_x \times N_x$. The blocks themselves are also tridiagonal and $N = N_x
\times N_y.$ Such matrices occur for instance in the context of Laplace like
equations in general coordinates.

Incomplete block LU factorization (IBLU) is known to be a powerful preconditioner
for Krylov space methods like CG, BiCGStab~\cite{vdV92} and GMRes~\cite{Saa86S}.
It also is a powerful smoother for multigrid methods \cite{Hem83KWdZ}.
The conventional IBLU algorithm is strongly recursive on two levels: at the
block level to make the LU decomposition and within a block to construct the
seven diagonals of the inverse of a tridiagonal matrix. We shall show, that both
these strong recursions can be replaced by algorithms that take typically
$O({}^2\log M)$ recursive cycles on both levels. The outer recursion then 
will take $O({}^2\log N_y)$ recursive cycles, and on the blocklevel $O({}^2\log
N_x)$ recursive cycles. 

\section {\label {sec:FCR}Full Cyclic Reduction for Tridiagonal Matrices}
We shall start out by taking a different view on cyclic reduction. Consider an
$M \times M$ 
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
$M$,
the size of the matrix, $L^k$ and $U^k$ will be empty and hence it will take 
${}^2\log M$ 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 $M$ or more processors the  number of  cycles 
would effectively reduce to $O({}^2\log M).$

\section {Incomplete Block LU Factorization (recursive)}
Various authors \cite{Axe84BI}, \cite{Con85GM},\cite{Kett_87} and \cite{Meu84} 
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_y$ blocks of size $N_x \times N_x$, $D^0_j, j=1\dots N_y$, the
superdiagonal $U^0$ of blocks $U^0_{j,j+1}, j=1\dots N_y-1$ and the
subdiagonal $L^0$ of blocks $L^0_{j,j-1}, j=2\dots N_y$. 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_y - 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_y -
2^k + 1, \dots, N_y
\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_y - 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 parallel. 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_y$ cycles on a massively parallel 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_x \times N_x$. As we have shown in
section~\ref{sec:FCR} this will take ${}^2\log N_x$ recursive cycles on a 
massively
parallel machine. So the whole operation will take $O(\log N_x \times \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_y - 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_y
\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_y - 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_y - 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
$M \times M$ and 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} < M \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, -M < m2^{k-j} < M$. 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 are only 
${}^2\log M$ recursive cycles in this calculation and that within such a cycle 
all calculations can be executed in parallel.

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