6.4. DESCRIPTION OF THE PLAPACK INFRASTRUCTURE

6.3.3.Elemental

Both ScaLAPACK and PLAPACK are relatively old software packages. With the advent of new many-core architectures, close to distributed-memory clusters inside a chip [81, 103], the interest in investigating and developing novel message-passing libraries is renewed, as they will have to be mapped to single-chip environments in a near future.

Elemental [113] is a new framework for the development of message-passing dense matrix computations that incorporates many of the insights introduced by the scientific community since the appearance of ScaLAPACK and PLAPACK one decade ago. Besides being a complete rewrite, with many of the functionality o ered by other dense linear algebra libraries, for distributed-memory architectures, Elemental proposes many novelties from the lowest level of its development and philosophy. Many of these new contributions are beyond the scope of this thesis, but the layered approach of the framework makes it feasible to port it to accelerated hybrid architectures in a similar way as the port of the PLAPACK infrastructure has been carried out in the framework of this dissertation.

6.4.Description of the PLAPACK infrastructure

6.4.1.Layered approach of PLAPACK

Figure 6.9 o ers an overview of the layered structure of the PLAPACK infrastructure. The components in the bottom line of the table are machine-dependent. The rest of the infrastructure was designed to be machine-independent. In fact, this is the key decision that makes it possible to perform a clean port to accelerated clusters without a ecting independent modules of the library.

Machine-dependent layer: PLAPACK aims at using standardized components, namely MPI for communication, BLAS for local kernel execution, standard memory management systems (malloc/calloc), etc.

Machine-independent layer: PLAPACK o ers a number of higher-level interfaces to the machinedependent layer. These interfaces include wrappers for MPI communications, memory management routines, distribution primitives, and BLAS interfaces.

PLAPACK abstraction layer: It provides an abstraction that isolates users from details such as indexing, communication, and local computation.

Linear algebra objects manipulation. All information related to the description of a linear algebra object (vectors or matrices) is stored in an opaque object; this is the same approach adopted years later by the FLAME methodology. This component of the abstraction layer allows the programmer to create, initialize, manipulate and destroy objects. It also provides mechanisms to perform indexing within matrices, sub-matrices and vectors.

Copy/reduce: duplication and consolidation of data. PLAPACK does not provide explicit routines for the communication of data. Instead, the objects themselves describe how data has to be distributed or duplicated. Communication itself is reduced to copy or reduction from one object to another.

PLAPACK local BLAS. Rather than explicitly requiring that the user extracts object data and invokes BLAS calls, PLAPACK provides a full set of routines to correctly

181

CHAPTER 6. MATRIX COMPUTATIONS ON CLUSTERS OF GPUS

Figure 6.9: PLAPACK software hierarchy. The PLAPACK abstraction layer (in red) is the only part of the infrastructure that needs to be modified in order to adapt it to an accelerated cluster.

extract data and call the appropriate local BLAS routine to perform per-process local BLAS invocations.

Library layer: PLAPACK was created not only as a library but also as an infrastructure to create new application-level libraries on top of it.

PLAPACK global BLAS. PLAPACK provides a full global BLAS implementation that allows the library developer to build higher-level libraries based on it without exposing communication routines.

PLAPACK higher-level linear algebra routines. Although these routines can be built directly from the global BLAS provided by the infrastructure, higher performance is expected if the abstraction layer functionalities are directly used. PLAPACK o ers implementation of a wide variety of LAPACK-level routines for message-passing architectures.

PLAPACK application interface: Even though the abstraction layer o ers routines to manage linear algebra objects, the PLAPACK API allows the programmer to directly manipulate the internals of those objects directly from the application.

6.4.2.Usage of the PLAPACK infrastructure. Practical cases

We next propose three di erent implementations of dense linear algebra operations to illustrate the main functionality of the PLAPACK infrastructure. We show the necessary steps to implement the matrix-vector multiplication, matrix-matrix multiplication and the Cholesky factorization using PLAPACK. Each example aims at showing a given functionality of the infrastructure. In the matrix-vector multiplication, we focus on how the communication is performed in PLAPACK. In the matrix-matrix multiplication, we illustrate how PLAPACK deals with partitioning and repartitioning objects. In the Cholesky factorization, we show how LAPACK-level routines can be implemented using the distributed BLAS routines provided by PLAPACK, extracting parallelism inside those BLAS calls transparently for the programmer.

182

6.4. DESCRIPTION OF THE PLAPACK INFRASTRUCTURE

Figure 6.10: Parallel matrix-vector multiplication. (a) Sub-vectors of x are distributed among processes in column-major order. (b) x is spread within columns. After that, local matrix-vector products can commence. (c) Each process has a partial contribution to the result, which needs to be summed within rows of processes. (d) Sub-vectors of y end up being distributed between processes in column-major order.

Parallelizing matrix-vector multiplication

To illustrate the use of PLAPACK as an infrastructure for developing new dense linear algebra operations, we show the parallel implementation of the dense matrix-vector multiplication using PLAPACK. We consider the basic case

y = Ax

where A is an m × n matrix. Three di erent objects must be created in the driver program to handle vectors x and y and matrix A. We will refer to these objects as x, y, and a. In addition, a template is created with a r × c communicator and block size nb (templates are the abstraction introduced by PLAPACK to describe internally the distribution of objects within processes). With this template, vectors x and y are aligned with the first element of the template vector, and A with the upper-left element of the template matrix.

Consider that x and y are distributed as vectors. In this case, we assume that x and y are identically distributed across the processes. Notice how, by spreading (collecting) vector x within columns, all the necessary elements of x are duplicated, so that local matrix vector multiplication

183

CHAPTER 6. MATRIX COMPUTATIONS ON CLUSTERS OF GPUS

can start on each process. After this local computation is done, a reduction (a summation in this case) within rows of processes of the local partial results yields vector y.

Therefore, three di erent steps are necessary to perform the operation:

1.Spread elements of x within columns of processes.

2.Perform the local matrix-vector multiply.

3.Perform a distributed reduce of the local results within rows. The global result is left in vector y.

We next illustrate how PLAPACK perform these steps. PLAPACK uses a mechanism to communicate based on the description of the initial and final distribution as objects, performing then a copy or reduction between them. All communication patterns in PLAPACK follow this approach. Thus, the following code produces the necessary spread of the elements of x within columns of processes:

As in the PLA Copy routine, all lower-level details underlying the reduction are encapsulated inside the routine PLA Reduce.

Note how, beyond the implementation details explained and those that can be obtained in the literature, there is an important factor that will have capital importance for the GPU implementations. All communications are encapsulated in only two routines: PLA Copy and PLA Reduce. The

184

6.4. DESCRIPTION OF THE PLAPACK INFRASTRUCTURE

developer is not exposed to any other communication primitive or detail, neither at a lower nor at a higher level. The code pattern is common for all communications: objects are created according to the way data has to be duplicated or reduced, and the proper PLAPACK communication routine is invoked. The infrastructure is designed to deal with lower level details without intervention of the programmer.

Parallelizing matrix-matrix multiplication

We next describe a blocked algorithm to illustrate how to perform a scalable matrix-matrix multiplication on message-passing architectures using PLAPACK. This approach was introduced in 1995 by van de Geijn and Watts [135], and is known as SUMMA (Scalable Universal Matrix Multiplication Algorithm). We only illustrate here the formation of the matrix product C := αAB + βC, where A, B, and C are matrices, and α, β are scalars.

For simplicity, consider that all three matrices A, B, and C are n × n. We use the following partitionings:

X = x0 x1 . . . xn−1 ,

with X A, B, C and x a, b, c, where xj is the j-th column of matrix X. Similarly

where xˆTi represents the i-th row of matrix X. Observe that

and thus the concurrent execution of the matrix-matrix multiplication can be performed as a set of rank-1 updates, with vectors y and x being the appropriate column and row of A and B, respectively.

Higher performance can be obtained if one replaces the matrix-vector multiplication with multiplication of a matrix by a panel of vectors. In this case, rank-1 updates are substituted by rank-k updates. Let A and B be partitioned as:

where A1 and B1 are m × s and s × n matrices, respectively. In that case

C ← βC

C ← α[A1B1 + A2B2],

which leads to the algorithm:

1.Scale C ← βC.

2.Let Acur = A and Bcur = B.

185

CHAPTER 6. MATRIX COMPUTATIONS ON CLUSTERS OF GPUS

3. Repeat until done

Determine width of the split. If Acur and Bcur are empty, break the loop.

where A1 and B1 are m × s and s × n matrices, respectively. Update C ← C + αA1B1.

Let Acur = A2 and Bcur = B2.

In this case, the symbol ← indicates assignment, while the symbol = indicates a reference or a partitioning. Note that this approach is usually referred as the panel-panel variant, as the basic operation uses panels of each matrix to perform the rank-k update. The code that follows gives the PLAPACK implementation for the operation. Note how PLAPACK provides routines to perform the partitioning and repartitioning of the matrices necessary in this algorithm. The codes also introduces the concept of projected multi-vector. A multi-vector can be considered as a group of vectors that are treated as a whole and are operated (and communicated) simultaneously. All vectors in a multi-vector are of equal length and present identical alignment. Intuitively, a multivector can be seen as a matrix with a few columns, where all columns are equally distributed, like vectors.

186

6.4. DESCRIPTION OF THE PLAPACK INFRASTRUCTURE

187