The Bubka Project

Towards the ultimate AMG parallel preconditioner Part I: A massively parallel IBLU preconditioner

Jos van Kan

Department of Technical Applied Analysis, Delft University of Technology, Delft, Netherlands

1 Introduction

The Bubka project aims at constructing the ultimate parallel preconditioner with

1. Multigrid efficiency

2. Robustness against anisotropy

3. 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.

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

1. 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.

2. 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.

3. 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?

2 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:

\begin{aligned}A^{1} & = (D^{0} + L^{0} + U^{0})(D^{0})^{−1}A \\ & = D^{0} − L^{0}(D^{0})^{−1}U^{0} − U^{0}(D^{0})^{−1}L^{0} − L^{0}(D^{0})^{−1}L^{0} − U^{0}(D^{0})^{−1}U^{0} \\ & = D^{1} − L^{1} − U^{1} \end{aligned} (1)

where

D^{1} = D^{0} − L^{0}(D^{0})^{−1}U^{0} − U^{0}(D^{0})^{−1}L^{0}  (2)
L^{1} = L^{0}(D^{0})^{−1}L^{0}  (3)
U^{1} = U^{0}(D^{0})^{−1}U^{0}  (4)

In these formulae D^{1} is again a diagonal matrix, and L^{1} and U^{1} are sub and superdiagonals respectively, once removed; that is, of the form l_{i,i−2}^{1} and u_{i,i+2}^{1} respectively. We can continue in the same vein and define recursively:

\begin{aligned}A^{k+1} & = (D^{k} + L^{k} + U^{k})(D^{k})^{−1}A \\ & = D^{k} − L^{k}(D^{k})^{−1}U^{k} − U^{k}(D^{k})^{−1}L^{k} − L^{k}(D^{k})^{−1}L^{k} − U^{k}(D^{k})^{−1}U^{k} \\ & = D^{k+1} − L^{k+1} − U^{k+1} \end{aligned} (5)

where

D^{k+1} = D^{k} − L^{k}(D^{k})^{−1}U^{k} − U^{k}(D^{k})^{−1}L^{k}  (6)
L^{k+1} = L^{k}(D^{k})^{−1}L^{k}  (7)
U^{k+1} = U^{k}(D^{k})^{−1}U^{k}  (8)

The reader will have no trouble to verify, that L^{k} and U^{k} are of the form l_{i,i−2^{k}}^{k} and u_{i,i+2^{k}}^{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).

3 Incomplete Block LU Factorization

Various authors [1], [4], [2],[3] 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 block tridiagonal matrix. We write (following an idea of Eisenstadt’s)

A = (D̃ − L)D̃^{−1}(D̃ − U)  (9)

and we choose D̃ such, that D = D̃ + LD̃^{−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 D̃, D̃_{1} are identical and we have the recursion

D̃_{j+1} = D_{j+1} + L_{j+1,j}D̃_{j}^{−1}U_{j,j+1}  (10)

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 three diagonals in every D̃_{j}. Define the operator E_{m} as a diagonal extractor. E_{0} extracts the main diagonal, E_{1} three diagonals (sub, main and super), E_{2} five, etcetera. With this notation the incomplete version of 10 can be constructed as follows:

D̃_{j+1} = D_{j+1} + E_{1}(L_{j+1,j}D̃_{j}^{−1}U_{j,j+1})  (11)

The reader will have no trouble in verifying, that this is equivalent to

D̃_{j+1} = D_{j+1} + E_{1}(L_{j+1,j}E_{3}(D̃_{j}^{−1})U_{j,j+1}))  (12)

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 D̃_{j}^{−1}. It is explained a.o. in [3] how to calculate those diagonals. This algorithm, however, is strongly recursive and unsuitable for parallel computation.

4 Full Cyclic Block Reduction

To find a parallel IBLU, let us reinterpret the formulae of section 2 in the following way. Let A be a block tridiagonal matrix. Let the diagonal D^{0} of A consist of N_{x} blocks of size N_{y} × N_{y}, D_{j}^{0}, j = 1 … N_{x}, the superdiagonal U^{0} of blocks U_{j,j+1}^{0}, j = 1 … N_{x} − 1 and the subdiagonal L^{0} of blocks L_{j,j−1}^{0}, j = 2 … N_{x}. Now we may perform the exact same transformations as in section 2 to obtain a full block cyclic reduction. In fact, there is a complete block analogue to equations 1-8, using the same formulae but with the connotation that it concerns block matrices.

However, equations 6-8 may also be written down blockwise. Equation 6 then becomes:

D_{j}^{k+1} = D^{k} − U_{j,j+2^{k}}^{k}(D_{j+2^{k}}^{k})^{−1}L_{j+2^{k},j}^{k},  j = 1, …, 2^{k}  (13)
D_{j}^{k+1} = D^{k} − L_{j,j−2^{k}}^{k}(D_{j−2^{k}}^{k})^{−1}U_{j−2^{k},j}^{k} +  (14)
− U_{j,j+2^{k}}^{k}(D_{j+2^{k}}^{k})^{−1}L_{j+2^{k},j}^{k},  j = 2^{k} + 1, …, N_{x} − 2^{k}
D_{j}^{k+1} = D^{k} − L_{j,j−2^{k}}^{k}(D_{j−2^{k}}^{k})^{−1}U_{j−2^{k},j}^{k}, j = N_{x} − 2^{k} + 1, …, N_{x}  (15)

Equation 7 becomes

L_{j+2^{k+1},j}^{k+1} = L_{j+2^{k+1},j+2^{k}}^{k}(D_{j+2^{k}}^{k})^{−1}L_{j+2^{k},j}^{k}, j = 1, …, N_{x} − 2^{k+1}  (16)

and finally equation 8 becomes

U_{j,j+2^{k+1}}^{k+1} = U_{j,j+2^{k}}^{k}(D_{j+2^{k}}^{k})^{−1}U_{j+2^{k},j+2^{k+1}}^{k}  (17)

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 2 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_{j}^{k})^{−1}, which are size N_{y} × N_{y}. As we have shown in section 2 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 × N_{x} × N_{y}^{2} ×^{2} log N_{x} to store the intermediate matrices, because starting from A^{1} the blocks will be full.

5 Incomplete Cyclic Block Reduction

To reduce the memory requirements one has to do an incomplete cyclic block reduction. Instead of using equations 13-17 we only retain the main, super and subdiagonal in every nonzero block. This gives the following set of recursions: For the diagonal blocks:

D_{j}^{k+1} = D^{k} − E_{1}(U_{j,j+2^{k}}^{k}E_{3}(D_{j+2^{k}}^{k}^{−1})L_{j+2^{k},j}^{k}),  j = 1, …, 2^{k} (18)
D_{j}^{k+1} = D^{k} − E_{1}(L_{j,j−2^{k}}^{k}E_{3}(D_{j−2^{k}}^{k}^{−1})U_{j−2^{k},j}^{k})+  (19)
− E_{1}(U_{j,j+2^{k}}^{k}E_{3}(D_{j+2^{k}}^{k}^{−1})L_{j+2^{k},j}^{k}),  j = 2^{k} + 1, …, N_{x} − 2^{k}
D_{j}^{k+1} = D^{k} − E_{1}(L_{j,j−2^{k}}^{k}E_{3}(D_{j−2^{k}}^{k}^{−1})U_{j−2^{k},j}^{k}), j = N_{x} − 2^{k} + 1, …, N_{x} (20)

for the subdiagonal blocks

L_{j+2^{k+1},j}^{k+1} = E_{1}(L_{j+2^{k+1},j+2^{k}}^{k}E_{3}(D_{j+2^{k}}^{k}^{−1})L_{j+2^{k},j}^{k}), j = 1, …, N_{x} − 2^{k+1}  (21)

and for the super diagonal blocks

U_{j,j+2^{k+1}}^{k+1} = E_{1}(U_{j,j+2^{k}}^{k}E_{3}(D_{j+2^{k}}^{k}^{−1})U_{j+2^{k},j+2^{k+1}}^{k}), j = 1, …, N_{x} − 2^{k+1} (22)

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

6 The Magnificent Seven

The problem remains, however, how to calculate E_{3}(A^{−1}) if A is tridiagonal. We recall from section 2 how the inverse was constructed:

A^{−1} = (D^{k+1})^{−1}E^{k}E^{k−1} … E^{1}E^{0}  (23)

with E^{j} = (L^{j} + D^{j} + U^{j})(D^{j})^{−1} and 2^{k−1} < N ≤ 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.

Lemma 1. The pattern of nonzero diagonals of the product P^{j} = E^{k} ·E^{k−1} ·· · ··E^{k−j} is given by d_{m2^{k−j}}, m ∈ ZZ, −N < m2^{k−j} < N. All other diagonals are 0.

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 ∈ ZZ 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}}. ⊓ ⊔

The Lemma is instrumental in the following theorem that enables us to calculate E_{3}(A^{−1}) in a parallel way.

Theorem 1. 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}.

Proof. Following the reasoning of Lemma 1 it will be clear, that to calculate the seven diagonals d_{−3}, d_{−2}, …, 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}. ⊓ ⊔

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.

References

1. O. Axelsson, S. Brinkkemper, and V.P. Il’in. On some versions of incomplete block matrix factorization iterative methods. Lin. Algebra Appl. 59, pages 3–15, 1984.

2. P. Concus, G.H. Golub, and G. Meurant. Block preconditioning for the conjugate gradient method. SIAM J. Sci. Stat. Comp. 6, pages 220–252, 1985.

3. Robert Kettler. Linear Multigrid Methods for Numerical Reservoir Simulation. University of Technology, Delft, december 1987. Thesis.

4. G. Meurant. The block preconditioned conjugate gradient method on vector computers. BIT, 24:623– 633, 1984.