Massively Parallel LINPACK Benchmark on the Intel Touchstone DELTA and iPSC/860 Systems

Progress Report

Robert A. van de Geijn ^{∗}

Department of Computer Sciences The University of Texas Austin, Texas 78712

December 5, 1991

Abstract

We describe an effort to implement the LINPACK Benchmark on two massively parallel distributed memory MIMD computers, the Intel iPSC/860 and DELTA Systems.

1 Introduction

For over a decade, the LINPACK benchmark has provided a measure of comparison between computers [6]. This benchmark reflects the performance of computers ranging from the home-used PC to the most powerful supercomputers when solving a dense system of linear equations, which is a component of many technical applications.

The original LINPACK benchmark is reported as three numbers: the performance, in MFLOPS, of the standard LINPACK code when applied to a 100 × 100 problem; the performance for a 1000 × 1000 problem is solved using an equivalent method; and the theoretical peak performance. The MFLOPS attained is computed by dividing the number of floating point operations required, 2n^{3}/3+ 2n^{2}, by the execution time.

More recently, the benchmark has been extended for massively parallel architectures to allow arbitrarily large problems to be used. The rationale is that such machines are meant to target larger problems than conventional computers. Applications that lead to large dense systems include problems employing method of moments or integral equation/Green’s function techniques. Again, the code must be equivalent to that of the original LINPACK benchmark.

∗This project was supported in part by Intel Corporation. Additional support came from the National Science Foundation Science and Technology Center for Research on Parallel Computation, Cooperative Agreement No. CCR-8809615, while the author was on leave at the University of Tennessee.

Several features of the Intel iPSC/860 and Touchstone Delta systems need to be considered in order to optimize the implementation of the LINPACK Benchmark. First, the i860 processor has a peak performance rate of 60 MFLOPS. However, for the kind of computation that is at the core of the dense linear solve, namely dense matrix-matrix multiplication, the theoretical peak is closer to 40 MFLOPS. A high quality implementation of the Level 1, 2, and 3 Basic Linear Algebra Subprograms (BLAS), developed by Kuck and Associates for the i860, is available from Intel and attains performance in the 30-35 MFLOPS range for matrix-matrix operations. ^{1}. Next, the Intel systems use a distributed memory MIMD architecture. Communication is required when more than one node is used to perform a given task. The cost of communication is such that the algorithm has to be reasonably coarse grained if high performance is to be achieved.

The above considerations led us to choose a right-looking variant of the LAPACK LU factorization [1] for the decomposition component of the code and a variant of the Li-Coleman (cyclic) algorithm for the triangular solves [16, 12, 10].

2 Mapping Matrices to Nodes

There are many ways that matrices can be distributed among nodes. When the number of nodes, p, is much less than the dimension of the matrix, n, 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). We will refer to this last mapping as the block-torus-wrapped mapping. In this case, n_{b} refers to the size of the block.

3 Parallel Implementation of the Algorithms

In this section, we briefly outline the parallel implementation of the primary components of the LINPACK Benchmark. All algorithms assume the matrix is distributed using the block-torus-wrapped mapping.

1More recently, new implementations of the BLAS have attained 36-37 Mflops per node. Unfortuntely, we did not get a chance to update our graphs for this paper.

3.1 LU Factorization

The LU factorization is a variant of the right-looking LAPACK LU factorization routine. Below, we give a brief sketch of the parallel implementation of this portion of the benchmark code.

Assume the LU factorization has proceeded so that all but the labeled portions of the matrix have been updated:

[figure: C; B; E]

where B ∈ R^{m×n_{b}}, C ∈ R^{n_{b}×(m−n_{b})}, and E ∈ R^{(m−n_{b})×(m−n_{b})}. During the next step, the rightlooking algorithm factors panel B, pivoting if necessary. Next, the pivots are applied to the remainder of the matrix. Blocks C and E now become blocks C̄ and Ē, a triangular solve updates submatrix C̄, and a rank n_{b} update updates submatrix Ē. This process continues with the updated matrix.

Turning now to the distributed memory implementation, assume the matrix is block-torus-wrapped onto a grid of nodes of dimension M × N, with block size n_{b}. The above described process proceeds as follows:

• The column of nodes that holds B colaborating to factor this panel.

• Pivot information is distributed to all other columns of nodes.

• Columns of nodes colaborate to pivot the remainders of the matrix rows.

• Factored panel B is distributed within rows of nodes.

• The row that holds C̄ performs the triangular solve, the results of which are distributed within columns of nodes.

• The update of Ē is performed in parallel.

Several implementation details are important to keep in mind:

• Broadcasting the factored panel can be done by sending the panel around the ring of node columns. This has proved to be efficient for column panel-wrapped algorithms, allowing communication and computation to be pipelined and is efficient for block-torus-wrapped approaches for similar reasons [2, 5].

• The broadcast of the updated C̄ must be done in a manner that requires the least total time. We have implemented this as a minimum spanning tree broadcast.

• During the factoring of panel B, determining the pivot row and the distribution of this row has been combined to avoid excessive startup overhead.

3.2 Triangular Solve

The upper and lower triangular solve routines are variants on the Li-Coleman (cyclic) algorithm. The major difference is that the Li-Coleman algorithm assumes that the matrix is column-wrapped onto an embedded ring of nodes, while our algorithm assumes a block-torus-wrapped mapping. Indeed, our algorithm can be seen as a cross between the Li-Coleman algorithm for a ring of nodes and an algorithm proposed by Johnsson [14] for a square grid of nodes, modified to act on blocks. Details go beyond the scope of this paper.

4 Experiments on the iPSC/860 and DELTA Systems

We have implemented the above algorithms on the Intel iPSC/860 and the Intel Touchstone DELTA systems. The iPSC/860 consists of 128 nodes connected by a hypercube network. Each node includes an Intel i860 microprocessor and 8 Mbytes of memory. While the DELTA is a mesh of 528 nodes, at this time only 512 can be used for a single application. The nodes of the DELTA include an Intel i860 and 16 Mbytes of memory. Our algorithms are implemented using the SPMD (Single Program, Multiple Data) paradigm and written in FORTRAN. Communication is implemented using a library that resembles a 2D version of the proposed Basic Linear Algebra Communication Subprograms (BLACS) [3], which makes the code portable to a range of multicomputers. The algorithms use Level 1-3 BLAS and are highly portable among distributed memory MIMD architectures. All experiments were performed in double precision arithmetic. We start by comparing the iPSC/860 and DELTA systems.

4.1 iPSC/860 vs. DELTA

For typical implementations of the BLAS, the performance changes with the size of the data being manipulated. Given a problem size, there is some optimal choice of panel width for a given algorithm. This is illustrated in Figure 1 for the benchmark code on 128 nodes of both the iPSC/860 and DELTA systems, configured as a 8 × 16 grid. In this figure, the space between two curves equals the time spent in the indicated operation. The performance of the BLAS matrix-matrix multiply routine dgemm dominates the computation. We conclude that n_{b} = 5 is a reasonable choice for the rest of our experiments.

A second way of tuning the performance is to vary the grid dimensions, given a fixed number of nodes. This is illustrated in Figure 2. Our current code restricts the row dimension to equal a power of 2 (due to the minimum spanning-tree broadcast). Notice that for p = 128, the optimal grid sizes are 8 × 16 and 4 × 32 for the iPSC/860 and DELTA, respectively. This can be attributed to the diffence in ratio between latency and bandwidth for the two machines, as well as the fact that the DELTA is physically a mesh, and hence the broadcasts can generate network conflicts.

Figure 1: Total computation time for 128 nodes, configured as a 8 × 16 grid, when n = 5000 and the block size n_{b} is varied. The space between two curves equals the time spent in the indicated operation.

Figure 2: Total computation time for 128 nodes, as a function of the row dimension of the grid of nodes, when n = 5000 and n_{b} = 5.

Figure 3: Proportion of time spent in various parts of the code, when the 128 nodes are configured as a 8 × 16 grid for the iPSC/860 and a 4 × 32 grid for the DELTA, n_{b} = 5, and n is varied.

For the optimal grid size and block size, the proportion of time spent in the various part of the code is reported in Figure 3, where the matrix size n is allowed to vary. In this figure, dgemm indicates time spent in the matrix-matrix multiply, bcast row panel indicates time spent in broadcasting the updated submatrix C̄, bcast col panel indicates time spent in broadcasting the factored panel B, and pivot indicates time spent pivoting matrix rows across rows of nodes. Finally, ltr solve and utr solve indicate time spent in performing the lower and upper triangular solve after the matrix has been factored.

Since communication overhead is considerably lower on the DELTA, better utilization of the 128 nodes can be realized, as illustrated in Figure 4.

4.2 Performance of the DELTA

Similar experiments were performed on the maximum number of nodes available to us on the DELTA, p = 512. We used the same block size, n_{b} = 5 for these experiments.

In Figure 5, the time is reported as a function of the row dimension of the node grid. This figure shows that configuring the nodes as a 8×64 logical grid is a good choice. In Figure 5, we report what proportion of time is spent in the various portions of the code as n is varied. The performance of the 512 mesh is contrasted to that of the 128 node DELTA and iPSC/860 in Figure 4. The DELTA completes the 25000 × 25000 problem in less than 15 minutes.

Figure 4: GFLOPS attained as the problem size is varied.

More recently, 13.9 GFLOPS was attained by this code on the Intel Touchstone Delta.

4.3 DELTA vs. the Rest

The final question that needs to be answered is how the DELTA measures up to other supercomputers according to the LINPACK benchmark. Performance of various architectures is reported in [6]. In this report, four numbers are reported for massively parallel computers.

• p, the number of nodes.

• r_{max}, the maximum GFLOPS attained for the benchmark.

• n_{max}, the size of the problem that yields r_{max}.

• n_{1/2}, the size of the problem that yields a performance of r_{max}/2.

• r_{peak}, the theoretical peak performance. The current ranking (as of June 20, 1991) of the most powerful parallel machines is given by:

Figure 5: LEFT: Total computation time for 512 nodes of the DELTA, as a function of the row dimension of the grid of nodes, when n = 10000 and n_{b} = 5. RIGHT: Proportion of time spent in various parts of the code, when the 512 nodes of the DELTA are configured as a 8 × 64 grid, n_{b} = 5, and n is varied.

Computer  p  \begin{matrix}r_{max} \\ GFLOPS \end{matrix}  \begin{matrix}n_{max} \\ order \end{matrix}  \begin{matrix}n_{1/2} \\ order \end{matrix}  \begin{matrix}r_{peak} \\ GFLOPS \end{matrix}
Intel Delta (40 MHz)  512  11.9  25000  7000  20 ^{2}
Thinking Machines CM-200 (10 MHz) 2048  9  28672  11264  20
Intel Delta (40 MHz)  256  5.9  18000  5000  10
Thinking Machines CM-2 (7 MHz)  2048  5.2  26624  11000  14
Intel Delta (40 MHz)  128  3.0  12500  3500  5
Intel iPSC/860 (40 MHz)  128  2.4  13000  4500  5
Intel Delta (40 MHz)  64  1.5  8000  3000  2.5

This compares very favorably with the fastest vector computers and shared memory multiprocessors. For example, a one processor NEC SX-3/14 attains 4.1 GFLOPS. However, such computers can attain this performance for much smaller problems (n = 1000).

5 Conclusion

In this paper, we have reported preliminary results of the performance of the LINPACK benchmark on the Intel machines. The figures show that there is still some room for improvement. In particular, it is still possible to fit larger problems in memory, which will result in better performance. The current (preliminary) implementation of the communication library does not make an effort to minimize network conflict on the DELTA. Finally, eventually it may be possible to utilize all 528 nodes of the DELTA.

2It should be noted that the theoretical peak for the Delta is 60 MFLOPS per node. However, the theoretical peak for matrix-matrix multiplication on a single node is considerably less, resulting in r_{peak} ≈ 20GFLOPS.

We will continue to update this report in the future. Eventually, we intend to include a rough time complexity analysis. This will give us further insight into whether communication overhead is mainly due to network latency or bandwidth.

6 Acknowledgements

Special thanks are given to the members of the Concurrent Supercomputing Consortium for their cooperation and contribution of time and access to the Touchstone Delta System. The author is also indebted to Sharon Brunett at Caltech for her help while performing the experiments on this machine. Further experiments were conducted on the Intel iPSC/860 with 8 MBytes of memory per node located at Oak Ridge National Laboratory. Addition experiments on a Intel iPSC/860 with 16 MBytes per node were performed by Jerry Baugh at Intel. Finally, we would like to acknowledge Jack Dongarra for numerous suggestions on how to correctly implement the LINPACK benchmark.

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] E. Anderson, A. Benzoni, J. Dongarra, S. Moulton, S. Ostrouchov, B. Tourancheau, and R. van de Geijn, LAPACK for distributed memory architectures: Progress report., in Parallel Processing for Scientific Computing, Fifth SIAM Conference, SIAM, 1991.

[3] , Basic Linear Algebra Communication Subprograms, in Sixth Distributed Memory Computing Conference Proceedings, IEEE Computer Society Press, 1991, pp. 287–290.

[4] A. Benzoni, V. S. Sunderam, and R. van de Geijn, Matrix factorization on a IBM RISC System/6000 workstation network, in Proceedings of the Second Symposium on High Performance Computing, Montpellier, France, 1991.

[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. J. Dongarra, Performance of various computers using standard linear equations software, Technical Report CS-89-85, University of Tennessee, Oct. 1989.

[7] 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.

[8] 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.

[9] 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.

[10] S. C. Eisenstat, M. T. Heath, C. S. Henkel, and C. H. Romine, Modified cyclic algorithms for solving triangular systems on distributed-memory multiprocessors, SIAM J. Sci. Stat. Comput., 9 (1988), pp. 589–600.

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

[12] M. T. Heath and C. H. Romine, Parallel solution of triangular systems on distributed memory multiprocessors, SIAM J. Sci. Stat. Comput., 9 (1988).

[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] S. L. Johnsson, Communication efficient basic linear algebra computations on hypercube architectures, J. of Par. Distr. Comput., 4 (1987), pp. 133–172.

[15] 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.

[16] G. Li and T. F. Coleman, A parallel triangular solver for a distributed-memory multiprocessor, SIAM J. Sci. Stat. Comput., 9 (1988), pp. 558–588.