LAPACK FOR DISTRIBUTED MEMORY ARCHITECTURES: PROGRESS REPORT

ED ANDERSON ^{∗}, ANNAMARIA BENZONI †, JACK DONGARRA ^{‡ ∗},

STEVE MOULTON ^{∗}, SUSAN OSTROUCHOV ^{∗},

BERNARD TOURANCHEAU §, AND ROBERT VAN DE GEIJN ¶ ^{∗}

Abstract. We describe an ongoing project that investigates the reusability of LAPACK code for distributed memory MIMD architectures. Related projects are described and experience with the Intel iPSC/860 is reported.

1. Introduction. In this paper, we discuss the current status of a project to implement a subset of routines from the LAPACK linear algebra subroutine library on distributed memory architectures. Parallel implementations of block algorithms for QR, LU, and Cholesky decomposition as well as Hessenberg and tridiagonal reduction for the algebraic eigenvalue problem are described. Results from experiments on the Intel iPSC/860 are reported.

The LAPACK project [1] is an effort to update the classical linear algebra software packages LIN-PACK and EISPACK to allow more efficient use of shared memory and traditional supercomputers. In addition, some of the algorithms have been made more robust.

The routines in LINPACK were written in terms of a library of routines known as the Basic Linear Algebra Subprograms (BLAS) [14]. A typical (Level 1) BLAS routine performs a vector operation. To improve performance on machines with hierarchical memories and/or vector processors, the BLAS were extended to include matrix-vector (Level 2 BLAS) and matrix-matrix (Level 3 BLAS) [8, 9] operations. In addition, these routines allow limited encapsulation of parallel operations, mostly on shared memory multiprocessors. LAPACK routines achieve efficiency by maximizing the amount of time spent in higher level BLAS.

In implementing LAPACK routines on distributed memory MIMD architectures, we decided to stress minimal change to the algorithms and codes, as well as portability of the final product. The codes are written in FORTRAN in the SPMD (Single Program Multiple Data) paradigm: all processor nodes of the multicomputer are driven by the same program. The path through the code

∗ Department of Computer Science, University of Tennessee, Knoxville. †IBM-ECSEC, Rome, Italy. ‡Mathematical Sciences Section, Oak Ridge National Laboratory. §LIP at ENS-Lyon, France. ¶On leave from the University of Texas at Austin.

This project was supported in part by the National Science Foundation Science and Technology Center Cooperative Agreement No. CCR-8809615 and the Applied Mathematical Sciences subprogram of the Office of Energy Research, U.S. Department of Energy, under Contract DE-AC05-84OR21400.

To appear in the proceedings of the Fifth SIAM Conference on Parallel Processing for Scientific Computing.

is determined by the relative position of the node. In addition, we have so far restricted ourselves to column-oriented mappings of matrices onto the nodes and one-dimensional communication. Portability of the codes is ensured by performing all communication between nodes through a proposed set of Basic Linear Algebra Communication Subprograms (BLACS).

This paper is organized as follows: We start by describing the model of parallel computing in Section 2. In Section 3 we describe the storage scheme used to date. A few of the algorithms that have been implemented are described in Section 4. Results from experiments on the Intel iPSC/860 are given in Section 5. Future directions are given in the conclusion.

2. Parallel Computing Model. We assume that our multicomputer consists of p nodes, indexed 0, …, p − 1. Each is equipped with CPU and local memory. The nodes are connected by some communication network that allows broadcasting of messages and combining of global data (in form of global summation), in addition to point-to-point communication. Messages that arrive over the network are buffered until absorbed by an appropriate call to a receive routine.

In our approach, it is actually not necessary to restrict the model to one process per node. Indeed, whenever we refer to a node, this could just as easily be a process, where many processes could be assigned to a single node.

3. Mapping Matrices to Nodes. There are many ways that matrices can be distributed among nodes. When p << n, the dimension of the matrix, the most common class of mappings requires the matrix to be “wrapped” onto the nodes. Assume n = mn_{b}, where m for convenience equals an integer. Partition

A = ( \begin{matrix}A_{11} & ⋯ & A_{1m} \\ ⋮ & & ⋮ \\ A_{m1} & ⋯ & A_{mm} \end{matrix} )

where each subblock A_{ij} is n_{b} × n_{b}. Then we define the (column) panel-wrapped storage scheme by the mapping that assigns A_{ij} to node (j − 1) mod p. Similarly, the row panel-wrapped scheme assigns A_{ij} to node (i − 1) mod p. Finally, if the nodes are indexed by (I, J), where 1 ≤ I ≤ M and 1 ≤ J ≤ N (NM = p), then the blocks can be wrapped onto this grid by assigning A_{ij} to node ((i − 1) mod M, (j − 1) mod N).

For the current generation of MIMD distributed memory computers, where p << n and communication is relatively expensive, column oriented assignment strategies are quite satisfactory and simplify the programming effort. Our effort to date has concentrated on the case where M = 1, i.e., (column) panel-wrapped storage. In this case, n_{b} represents the width of the panel. Panel-wrapped storage is a generalization of the familiar column-wrapped storage scheme, which is equivalent to panel-wrapped storage with n_{b} = 1 [13].

4. Parallel Implementation of the Algorithms. In this section, we outline the parallel implementation of the various algorithms. All are based on right-looking variants of the algorithms. Details can be found in [5, 7].

4.1. Factorization Algorithms.

LU Factorization: Let

A = (A_{1}, A_{2}, …, A_{m})

where A_{i} is a panel of width n_{b} (except for the last panel, which may be of smaller width). Assume these have been assigned to nodes in a panel-wrapped fashion.

Assume the algorithm has proceeded to the point where panels A_{i}, i < k, have been factored and the remaining panels have been updated accordingly. At the next step (step k), node (k − 1) mod p factors panel A_{k} and broadcasts the pivot and multiplier information to all other nodes. Next, all nodes update their portion of panels A_{(k+1)}, …, A_{m}. This algorithm pipelines nicely around a ring, allowing step k + 1 to start as soon as enough information is available. Moreover, a simple rearrangement of the computation known as compute-ahead allows some communication and computation to be overlapped.

Cholesky Factorization: This algorithm proceeds similarly, except that due to symmetry only the lower triangular portion is stored and updated. Moreover, no pivoting is necessary since the matrix is assumed to be positive definite.

QR Factorization: Again, this algorithm proceeds as the LU factorization, except that it uses Householder transformations to factor the panels. Again, no pivoting is performed.

4.2. Reduction Algorithms. The only communication required by the factorization algorithms consists of broadcasting information from the node that performs the factorization of the panel to all other nodes. By contrast, the reduction algorithms require not only a broadcast, but also one or more global summations of vectors across nodes. This means that all processors are involved in reducing a panel. Moreover, the global summations act as synchronization points.

While details of these algorithms go beyond the scope of this paper, the crucial step that sets them apart from the factorization algorithms is the necessity to repeatedly form matrix vector products that require data from different nodes to be combined: Let A be partitioned and distributed as before. Assume x^{T} = (x_{1}^{T} ⋯ x_{m}^{T}), is distributed so that x_{i} ∈ R^{n_{b}} is assigned to node (i−1) mod p. Then y = Ax can be formed in two steps, leaving the result on all nodes: First form

y^{i} = ∑_{A_{j} onnode i} A_{j}x_{j}

on each node i. Next, all nodes collaborate to form y = ∑_{i=0}^{p} y^{i}.

Our reduction to Hessenberg form requires at least n_{b} global summations of vectors of length n and n_{b} − 1 global summations of length ranging from 1 to n_{b} − 1, in addition to n_{b} broadcasts, per panel factored.

5. Experiments on the Intel iPSC/860. We have implemented the above algorithms on the Intel iPSC/860 at Oak Ridge National Laboratory, which consists of 128 nodes connected by a hypercube network. Each node consists of an Intel i860 microprocessor and 8 Mbytes of memory.

Our algorithms are implemented using the SPMD (Single Program, Multiple Data) paradigm and written in FORTRAN. Communication is implemented using a proposed set of Basic Linear Algebra Communication Subprograms (BLACS), which makes the code portable to a range of multicomputers. The algorithms use Level 1-3 BLAS (Basic Linear Algebra Subprograms) implemented by Kuck and Associates. Since these BLAS are assembly coded, reasonable performance can be attained, in the 50-65 MFLOPS per node range, compared to a peak performance of 80 MFLOPS. All experiments were performed in single precision arithmetic.

For typical implementations of the BLAS, the performance increases with the size of the data being manipulated. For the QR factorization and the reduction algorithms, this is offset to some degree by an increase in the operation count as the panel width increases. Given a problem size, there is some optimal choice of panel width for a given algorithm.

These same issues arise for the parallel implementation. However, the panel width also influences communication overhead and processor idle time. The performance of the algorithms as a function of the panel width n_{b} is illustrated in Fig. 1. In this figure, we also report results from profiling the Hessenberg reduction as a function of the panel width. The space between two curves equals

Fig. 1. LEFT: Performance of the algorithms as a function of n_{b} when p = 64 and n = 5000. RIGHT: Results from profiling the Hessenberg reduction as a function of the panel width is plotted. The area between two curves indicates the time spent in the indicated operation.

Fig. 2. Performance and execution time of the algorithms when p = 128, n_{b} equals the optimal panel width, and n is varied.

Fig. 3. Profile of the QR factorization (left) and Hessenberg reduction (right) when p = 128, n_{b} equals the optimal panel width, and n varies. The graphs report percentage of time on average spent in the various operations.

the time spent in the indicated operation. The erratic behavior of the BLAS3 curve is due to the implementation of the matrix-matrix multiply.

Fig. 2 shows the performance and execution time on 128 nodes when the problem size is varied and the panel width is fixed at the optimal for each of the algorithms, as determined from Fig.

1. Naturally, the performance is related to the ratio of computation to communication. For a given problem size, the communication overhead the factorization algorithms is essential the same. However, the operation counts for the sequential algorithms are about 4/3n^{3}, 2/3n^{3}, and 1/3n^{3}, for the QR, LU, and Cholesky factorizations, respectively. This accounts for the diminished performance of the latter two parallel algorithms.

In addition to broadcast overhead, the reduction algorithms incur the overhead due to numerous global summations. The difference is illustrated for the QR factorization and Hessenberg reduction in Fig. 3. For the Hessenberg reduction, most of the idle time is due to the nodes waiting for the broadcast arrive, before participating in relaying the information. Communication is similar for the parallel implementations of the two reduction algorithm, although the time spent in global summation is about half for the tridiagonal reduction. Operation counts for the sequential Hessenberg and tridiagonal reduction algorithms equal about 10/3n^{3} and 4/3n^{3}, respectively.

It should be noted that the broadcast for the parallel factorization algorithms consists of the information being passed around the ring, allowing some communication and computation to be overlapped. By contrast, the reduction algorithms are very synchronous in nature, due to the global summations, and a more appropriate minimum spanning tree broadcast is used.

6. Conclusion. For networks of relatively few powerful nodes (p << n), the approach we have taken for implementing the algorithms discussed in the previous sections appears to be both relatively straight forward and reasonably efficient. Performance on the individual nodes is enhanced by the use of higher level BLAS and portability is ensured through the use of the BLACS. Communication overhead is considerable. This can be attributed largely to the limited bandwidth between nodes on the iPSC/860: while each node has a peak performance of 80 MFLOPS, the communication bandwidth between nodes is only about 2 Mbytes per second. This is expected to improve by an order of magnitude for the next generation iPSC/860.

We are in the process of implementing algorithms for solving triangular systems arising from the discussed factorizations, as well as matrix inversion. A real challenge will be the implementation of algebraic eigenvalue algorithms. There are a lot of possibilities for the symmetric eigenvalue problem, like divide and conquer strategies [6, 15] in addition to Jacobi methods [3, 4, 18]. Parallel algorithms for nonsymmetric matrices appear to be more difficult to implement [11, 19, 2].

REFERENCES

[1] E. Anderson, Z. Bai, C. Bischof, J. Demmel, J. Dongarra, J. Du Croz, A. Greenbaum, S. Hammarling, A. McKenney, and D. Sorensen, LAPACK: A portable linear algebra library for high-performance computers, in Proceedings Supercomputing ’90, IEEE Computer Society Press, Los Alamitos, California, 1990, pp. 2–11.

[2] D. Boley and R. Maier, A parallel QR algorithm for the nonsymmetric eigenvalue problem, Computer Science report TR-88-12, Univ. of Minn., 1988.

[3] R. Brent and F. Luk, The solution of singular value and symmetric eigenvalue problems on multi-processor arrays, SIAM J. Sci. Stat. Comput., 6 (1985), pp. 69–84.

[4] J. Demmel and K. Veselic, Jacobi’s method is more accurate than QR, LAPACK Working Note 15, Technical Report CS-89-88, University of Tennessee, Oct. 1989.

[5] J. Dongarra and S. Ostrouchov, LAPACK block factorization algorithms on the Intel iPSC/860, LAPACK Working Note 24, Technical Report CS-90-115, University of Tennessee, Oct. 1990.

[6] J. Dongarra and D. Sorensen, A fully parallel algorithm for the symmetric eigenvalue problem, SIAM J. Sci. Stat. Comput., 8 (1987).

[7] J. Dongarra and R. van de Geijn, Reduction to condensed form for the eigenvalue problem on distributed memory architectures, LAPACK Working Note 30, Technical Report CS-91-130, University of Tennessee, 1991.

[8] J. J. Dongarra, J. Du Croz, S. Hammarling, and I. Duff, A set of level 3 basic linear algebra subprograms, ACM Trans. Math. Soft., 16 (1990), pp. 1–17.

[9] J. J. Dongarra, J. Du Croz, S. Hammarling, and R. J. Hanson, An extended set of FORTRAN basic linear algebra subprograms, ACM Trans. Math. Soft., 14 (1988), pp. 1–17.

[10] T. Dunigan, Performance of the Intel iPSC/860 and NCUBE 6400 hypercubes, Tech. Report ORNL/TM-11790, Oak Ridge National Laboratory, Oak Ridge, Tennessee, 1991.

[11] G. Geist and G. Davis, Finding eigenvalues and eigenvectors of unsymmetric matrices using a distributedmemory multiprocessor, Parallel Computing, 13 (1990), pp. 199–209.

[12] G. H. Golub and C. F. V. Loan, Matrix Computations, The Johns Hopkins University Press, Baltimore, second ed., 1989.

[13] I. Ipsen, Y. Saad, and M. Schultz, Complexity of dense-linear-system solution on a multiprocessor ring, Lin. Alg. Appl., 77 (1986), pp. 205–239.

[14] C. L. Lawson, R. J. Hanson, D. R. Kincaid, and F. T. Krogh, Basic linear algebra subprograms for Fortran usage, ACM Trans. Math. Soft., 5 (1979), pp. 308–323.

[15] S. Lo, B. Philippe, and A. Sameh, A multiprocessor algorithm for the symmetric tridiagonal eigenvalue problem, SIAM J. Sci. Stat. Comput., 8 (1987).

[16] C. Moler, J. Little, and S. Bangert, Pro-Matlab, User’s Guide, The Mathworks, Inc., 1987.

[17] R. van de Geijn, On global combine operations, LAPACK Working Note 29, Technical Report CS-91-129, University of Tennessee, 1991.

[18] , Storage schemes for parallel eigenvalue algorithms, in Numerical Linear Algebra, Digital Signal Processing and Parallel Algorithms, G. Golub and P. van Dooren, eds., NATO ASI Series, Springer Verlag, 1991, pp. 639–648.

[19] R. van de Geijn and D. Hudson, An efficient parallel implementation of the nonsymmetric QR algorithm, in Proceedings of the Fourth Conference on Hypercube Concurrent Computers and Applications, 1989.