Jack Dongarra
email: dongarra@cs.utk.edu

Jean-François Pineau
email: jean-francois.pineau@ens-lyon.fr

Yves Robert
email: yves.robert@ens-lyon.fr

Frédéric Vivien
email: frederic.vivien@ens-lyon.fr

Matrix Product on Heterogeneous Master-Worker Platforms

Keywords: F.2.1 [Analysis of algorithms and problem complexity]: Numerical Algorithms and Problems -Computations on matrices; F.2.2 [Analysis of algorithms and problem complexity]: Nonnumerical Algorithms and Problems -Scheduling Algorithms, Experimentation, Theory Matrix product, limited memory, communication

This paper is focused on designing efficient parallel matrix-product algorithms for heterogeneous master-worker platforms. While matrix-product is well-understood for homogeneous 2D-arrays of processors (e.g., Cannon algorithm and ScaLAPACK outer product algorithm), there are three key hypotheses that render our work original and innovative: -Centralized data. We assume that all matrix files originate from, and must be returned to, the master. The master distributes data and computations to the workers (while in ScaLAPACK, input and output matrices are supposed to be equally distributed among participating resources beforehand). Typically, our approach is useful in the context of speeding up MATLAB or SCILAB clients running on a server (which acts as the master and initial repository of files).

-Heterogeneous star-shaped platforms. We target fully heterogeneous platforms, where computational resources have different computing powers. Also, the workers are connected to the master by links of different capacities. This framework is realistic when deploying the application from the server, which is responsible for enrolling authorized resources.

-Limited memory. As we investigate the parallelization of large problems, we cannot assume that full matrix column blocks can be stored in the worker memories and be re-used for subsequent updates (as in ScaLAPACK).

We have devised efficient algorithms for resource selection (deciding which workers to enroll) and communication ordering (both for input and result messages), and we report a set of numerical experiments on a platform at our site. The experiments show that our matrix-product algorithm has smaller execution times than existing ones, while it also uses fewer resources.

Introduction

Matrix product is a key computational kernel in many scientific applications and it has been extensively studied on parallel architectures. Two well-known parallel versions are Cannon's algorithm [START_REF] Cannon | A cellular computer to implement the Kalman filter algorithm[END_REF] and the ScaLAPACK outer product algorithm [START_REF] Blackford | ScaLAPACK Users' Guide[END_REF]. Typically, parallel implementations work well on 2D processor grids, because the input matrices are sliced horizontally and vertically into square blocks that are mapped one-to-one onto the physical resources; several communications can take place in parallel, both horizontally and vertically. Even better, most of these communications can be overlapped with (independent) computations. All these characteristics render the matrix product kernel quite amenable to an efficient parallel implementation on 2D processor grids.

However, current architectures typically take the form of heterogeneous clusters, which are composed of heterogeneous computing resources, interconnected by a sparse network; there are no direct links between any pair of workers. Therefore, an accurate estimation of the communication cost requires precise knowledge of the underlying target platform. In addition, it becomes necessary to include the cost of both the initial distribution of the matrices to the processors and of collecting back the results. These input/output operations have always been neglected in the analysis of the conventional algorithms. This is because only Θ(n 2 ) coefficients need to be distributed in the beginning, and gathered at the end, as opposed to the Θ(n 3 ) computations 1 to be performed (where n is the problem size). The assumption that these communications can be ignored could have made sense on dedicated processor grids like, say, the Intel Paragon, but it is no longer reasonable on heterogeneous platforms. Furthermore, when processors cannot store all the matrices in their memory, the total volume of communication required can be larger than Θ(n 2 ), as a same matrix element may have to be sent several times to a same processor.

In this paper, we are not interested in adapting 2D processor grid strategies to heterogeneous clusters, as proposed in [START_REF] Kalinov | Heterogeneous distribution of computations solving linear algebra problems on networks of heterogeneous computers[END_REF][START_REF] Beaumont | A proposal for a heterogeneous cluster ScaLAPACK (dense linear solvers)[END_REF][START_REF] Beaumont | Matrix multiplication on heterogeneous platforms[END_REF]. Instead, we adopt a realistic application scenario, where input files are read from a fixed repository (such as a disk on a data server). Computations are delegated to available computational resources, and results will be returned to the repository. This calls for a masterworker paradigm, or more precisely for a computational scheme where the master (the processor holding the input data) assigns computations to other resources, the workers. In this centralized approach, all matrix files originate from, and must be returned to, the master. The master distributes both data and computations to the workers (while in ScaLAPACK, input and output matrices are supposed to be equally distributed among participating resources beforehand). Typically, our approach is useful in the context of speeding up MATLAB or SCILAB clients running on a server (which acts as the master and initial repository of files).

We target fully heterogeneous master-worker platforms, where computational resources have different computing powers, and workers are connected to the master by links of different capacities. This framework is realistic when deploying the application from a server, which is responsible for enrolling authorized resources.

Finally, because we investigate the parallelization of large problems, we cannot assume that full matrix panels, or column blocks, can be stored in worker memories and re-used for subsequent updates (as in ScaLAPACK). The amount of memory available in each worker is expressed as a given number mi of buffers, where a buffer can store a square block of matrix elements. The size q of these square blocks is chosen so as to harness the power of Level 3 BLAS routines [START_REF] Blackford | ScaLAPACK Users' Guide[END_REF]: q = 80 or 100 on most platforms.

To summarize, the target platform is composed of several workers with different computing powers, different bandwidth links to/from the master, and different, limited, memory capacities. The first problem is resource selection. Which workers should be enrolled in the execution: all of them? only the fastest computing ones? only the fastest-communicating ones? Once participating resources have been selected, there remain several scheduling problems: how to minimize the number of communications? in which order workers should receive input data and return results? what amount of communications can be overlapped with (independent) computations? The main contributions of this paper are twofold:

• On the theoretical side, we have refined the existing bounds on the volume of communications needed to perform a matrixproduct on a platform with insufficient memory to simultaneously store the whole three matrices. • On the practical side, we have designed an algorithm for heterogeneous platforms which is quicker than the existing ones, and which uses less computational resources, according to our MPI experiments. The rest of the paper is organized as follows. In Section 2, we state the scheduling problem precisely, and we introduce notations. Next, in Section 3, we proceed with the analysis of the total communication volume that is needed in the presence of memory constraints, and we improve a well-known bound by Toledo [START_REF] Toledo | A survey of out-of-core algorithms in numerical linear algebra[END_REF][START_REF] Ironya | Communication lower bounds for distributed-memory matrix multiplication[END_REF]. In order to help the reader apprehend the solution for heterogeneous platforms, we first deal with homogeneous platforms in Section 4, and we propose a scheduling algorithm that includes resource selection. Section 5 addresses the design of algorithms for heterogeneous platforms, which turns out to be a truly challenging algorithmic task. We report several MPI experiments in Section 6. Section 7 is devoted to an overview of related work. Finally, we state some concluding remarks in Section 8.

Framework Application

We deal with the computational kernel C ← C+A×B. We partition the three matrices A, B, and C as illustrated in Figure 1. Formally:

• We use a block-oriented approach. The atomic elements that we manipulate are not matrix coefficients but instead square blocks of size q × q (hence with q 2 coefficients). This is to harness the power of Level 3 BLAS routines [START_REF] Blackford | ScaLAPACK Users' Guide[END_REF]. Typically, q = 80 or 100 when using ATLAS-generated routines [START_REF] Whaley | Automated empirical optimizations of software and the atlas project[END_REF]. • The input matrix A is of size nA × nAB:

-we split A into r horizontal stripes Ai, 1 ≤ i ≤ r, where r = nA/q; -we split each stripe Ai into t square q × q blocks A i,k , 1 ≤ k ≤ t, where t = nAB/q. • The input matrix B is of size nAB × nB:

-we split B into s vertical stripes Bj, 1 ≤ j ≤ s, where s = nB/q; -we split stripe Bj into t square q × q blocks B k,j , 1 ≤ k ≤ t.

• We compute C = C + A × B. Matrix C is accessed (both for
input and output) by square q × q blocks Ci,j, 1 ≤ i ≤ r, 1 ≤ j ≤ s. There are r × s such blocks. We point out that with such a decomposition, all stripes and blocks have the same size. This will greatly simplify the analysis of communication costs.

Platform

We target a star network S = {P0, P1, P2, . . . , Pp}, composed of a master P0 and of p workers Pi, 1 ≤ i ≤ p.In practice, this star network can be a logical overlay built upon a different physical interconnection network. Because we manipulate large data blocks, we adopt a linear cost model, both for computations and communications (i.e., we neglect start-up overheads). We have the following notations:

• It takes X.wi time-units to execute a task of size X on Pi;

• It takes X.ci time units for the master P0 to send a message of size X to Pi or to receive a message of size X from Pi.

Our star platforms are thus fully heterogeneous, both in terms of computations and of communications. A fully homogeneous star platform would be a star platform with identical workers and identical communication links: wi = w and ci = c for each worker Pi, 1 ≤ i ≤ p. Without loss of generality, we assume that the master has no processing capability (otherwise, add a fictitious extra worker paying no communication cost to simulate computation at the master).

For the communication model, we adopt the one-port model [START_REF] Bhat | Efficient collective communication in distributed heterogeneous systems[END_REF], which is defined as follows: 1) the master can only send data to, and receive data from, a single worker at a given time-step; 2) a given worker cannot start execution before it has terminated the reception of the message from the master; similarly, it cannot start sending the results back to the master before finishing the computation; 3) a given worker can overlap computation and communication of independent tasks.

The one-port model is realistic. Bhat, Raghavendra, and Prasanna [START_REF] Bhat | Efficient collective communication in distributed heterogeneous systems[END_REF] advocate its use because "current hardware and software do not easily enable multiple messages to be transmitted simultaneously." Even if non-blocking, multi-threaded communication libraries allow for initiating multiple send and receive operations, they claim that all these operations "are eventually serialized by the single hardware port to the network." Experimental evidence of this fact has recently been reported by Saif and Parashar [START_REF] Saif | Understanding the behavior and performance of non-blocking communications in MPI[END_REF], who report that asynchronous MPI sends get serialized as soon as message sizes exceed a hundred kilobytes. Their results hold for two popular MPI implementations, MPICH on Linux clusters and IBM MPI on the SP2. Note that all the MPI experiments in Section 6 obey the one-port model.

Our final assumption is related to memory capacity; we assume that a worker Pi can only store mi blocks (either from A, B, or C). For large problems, this memory limitation will considerably impact the design of the algorithms, as data re-use will be greatly dependent on the amount of available buffers.

Minimization of the communication volume

In this section, we derive a lower bound on the total number of communications (sent from, or received by, the master) that are needed to execute any matrix multiplication algorithm. Any parallel algorithm can always be simulated on a single worker. Therefore, since we are not interested in optimizing the execution time, but only in minimizing the total communication volume, we only need to consider the one-worker case. We then deal with the following formulation of the problem:

• The master sends blocks A ik , B kj , and Cij,

• The master retrieves final values of blocks Cij, and

• We enforce limited memory on the worker; only m buffers are available, which means that at most m blocks of A, B, and/or C can simultaneously be stored on the worker. First, we improve the lower bound on the communication volume established by Toledo et al. [START_REF] Toledo | A survey of out-of-core algorithms in numerical linear algebra[END_REF][START_REF] Ironya | Communication lower bounds for distributed-memory matrix multiplication[END_REF]. Then, we describe an algorithm that aims at re-using C blocks as much as possible after they have been loaded, and we assess its performance.

Lower bound on the communication volume

To derive the lower bound, we refine an analysis due to Toledo [START_REF] Toledo | A survey of out-of-core algorithms in numerical linear algebra[END_REF]. The idea is to estimate the number of computations made thanks to m consecutive communication steps (once again, the unit here is a matrix block). We need some notations:

• We let α old , β old , and γ old be the number of buffers used by blocks of A, B, and C right before the beginning of the m communication steps; • We let αrecv, βrecv, and γrecv be the number of A, B, and C blocks sent by the master during the m communication steps; • Finally, we let γ send be the number of C blocks returned to the master during these m steps. Initially, the memory contains at most m blocks, and we consider a sequence of m communications. Therefore:

α old + β old + γ old ≤ m αrecv + βrecv + γrecv + γ send = m
We then use Loomis-Whitney inequality [START_REF] Ironya | Communication lower bounds for distributed-memory matrix multiplication[END_REF]: in any algorithm that uses the standard way of multiplying matrices, if NA elements of A, NB elements of B, and NC elements of C are accessed, no more than K computations can be done, where K = √ NANBNC . Here K = (α old + αrecv)(β old + βrecv)(γ old + γrecv). K is maximized when α old + αrecv = β old + βrecv = γ old + γrecv = The maximum re-use algorithm

In the above study, the lower-bound on the communication volume is obtained when the three matrices A, B, and C are equally accessed during a sequence of computations. This may suggest allocating one third of the memory to each of these matrices. In fact, Toledo [START_REF] Toledo | A survey of out-of-core algorithms in numerical linear algebra[END_REF] uses this memory layout. He even proves that, in the context of multiplication of square matrices of size r, his algorithm is "asymptotically optimal" as soon as the processor cannot store in its memory more than one sixth of one of the matrices: such an algorithm must have a communication-per-computation ratio which is

C C C C C C C C C C C C C C C C A B B B B 1 µ µ 2
Ω r 3 √ m
when his algorithm has a communication to computa-

tion ratio of O r 3 √ m .
We can, however, still significantly improve the performance of matrix multiplication in our context by reducing the constant hidden in the order of complexity of this ratio, as our experiments will show in Section 6.

A closer look at our problem shows that the multiplied matrices A and B have the same behavior, which differs from the behavior of the result matrix C. Indeed, if an element of C is no longer used, it cannot be simply discarded from the memory as the elements of A and B are, but it must be sent back to the master. Intuitively, sending an element of C to a worker also costs the communication needed to retrieve it from the worker, and is thus twice as expensive as sending an element of A or B. Therefore, we designed an algorithm that reuses as much as possible the elements of C.

Cannon's algorithm [START_REF] Cannon | A cellular computer to implement the Kalman filter algorithm[END_REF] and the ScaLAPACK outer product algorithm [START_REF] Blackford | ScaLAPACK Users' Guide[END_REF] both distribute square blocks of C to the processors. Intuitively, squares are better than elongated rectangles because their perimeter (which is proportional to the data needed by a worker to compute the area) is smaller for the same area. We use the same approach here.

The maximum re-use algorithm uses the memory layout illustrated in Figure 2. Four consecutive execution steps are shown in Figure 3. Assume that there are m available buffers. First we find µ as the largest integer such that 1 + µ + µ 2 ≤ m. The idea is to use one buffer to store A blocks, µ buffers to store B blocks, and µ 2 buffers to store C blocks. In the outer loop of the algorithm, a µ × µ square of C blocks is loaded. Once these µ 2 blocks have been loaded, they are repeatedly updated in the inner loop of the algorithm until their final value is computed. Then the blocks are returned to the master, and µ 2 new C blocks are sent by the master and stored by the worker. As illustrated in Figure 2, we need µ buffers to store a row of B blocks, but only one buffer for A blocks: A blocks are sent in sequence, each of them is used in combination with a row of µ B blocks to update the corresponding row of C blocks. This leads to the following sketch of the algorithm: Outer loop: while there remain C blocks to be computed

• Store µ 2 blocks of C in worker's memory: {Ci,j | i0 ≤ i < i0 + µ, j0 ≤ j < j0 + µ}. • Inner loop: For each k from 1 to t: 1. Send a row of µ elements {B k,j | j0 ≤ j < j0 + µ}; 2. Sequentially send µ elements of column {A i,k | i0 ≤ i < i0 + µ}. For each A i,k , update µ elements of C • Return results to master.
The performance of one iteration of the outer loop of the maximum re-use algorithm can readily be determined:

• We need 2µ 2 communications to send and retrieve C blocks.

• For each value of t:

-we need µ elements of A and µ elements of B; -we update µ 2 blocks. In terms of block operations, the communication-to-computation ratio achieved by the algorithm is thus

C14 C13 C12 C11 C41 C42 C43 C44 C34 C32 C33 C31 C24 C21 C22 C23 B14 B13 B12 B11 A11 C22 C21 C24 C23 C41 C42 C43 C44 C34 C32 C33 C31 C14 C13 C12 C11 B14 B13 B12 B11 A21 C34 C33 C32 C31 C41 C42 C43 C44 C24 C21 C22 C23 C14 C13 C12 C11 B14 B13 B12 B11 A31 C44 C43 C41 C42 C34 C32 C33 C31 C24 C21 C22 C23 C14 C13 C12 C11 B14 B13 B12 B11 A41
CCR = 2µ 2 + 2µt µ 2 t = 2 t + 2 µ .
For large problems, i.e., large values of t, we see that CCR is asymptotically close to the value CCR∞ = 2 √ m . We point out that, in terms of data elements, the communication-to-computation ratio is divided by a factor q. Indeed, a block consists of q 2 coefficients but an update requires q 3 floating-point operations. Also, the ratio CCR∞ achieved by the maximum re-use algorithm is lower by a factor √ 3 than the ratio achieved by the blocked matrix-multiply algorithm of [START_REF] Toledo | A survey of out-of-core algorithms in numerical linear algebra[END_REF]. Finally, we remark that the performance of the maximum re-use algorithm is quite close to the lower bound derived

earlier: CCR∞ = 2 √ m = 32 8m .

Algorithms for homogeneous platforms

We now adapt the maximum re-use algorithm to fully homogeneous platforms. We have a limitation on the memory capacity. So we must first decide which part of the memory will be used to stock which part of the original matrices, in order to maximize the total number of computations per time unit.

Principle of the algorithm

We load into the memory of each worker µ blocks of A and µ blocks of B to compute µ 2 blocks of C. In addition, we need 2µ extra buffers, split into µ buffers for A and µ for B, in order to overlap computation and communication steps. In fact, µ buffers for A and µ for B would suffice for each update, but we need to prepare for the next update while computing. Overall, the number µ 2 of C blocks that we can simultaneously load into memory is defined by the largest integer µ such that µ 2 + 4µ ≤ m.

We have to determine the number of participating workers P. On the communication side, we know that in a round (computing a C block entirely), the master exchanges with each worker 2µ 2 blocks of C (µ 2 sent and µ 2 received), and sends µt blocks of A and µt blocks of B. During this round, on the computation side, each worker computes µ 2 t block updates. If we enroll too many processors, the communication capacity of the master will be exceeded: there is a limit on the number of blocks that it can send per time unit. On the contrary, if we enroll too few processors, Split matrix C into squares C i ,j of µ 2 q × q blocks : C i ,j = {Ci,j | (i -1)µ+1 ≤ i ≤ i µ, (j -1)µ+1 ≤ j ≤ j µ} for j ← 0 to s Pµ by Step P do for i ← 1 to r µ do for id worker ← 1 to P do j ← j + id worker Send block C i ,j to worker id worker for k ← 1 to t do for id worker ← 1 to P do j ← j + id worker for j ← (j -

1)µ + 1 to j µ do Send B k,j for i ← (i -1)µ + 1 to i µ do Send A i,k
for id worker ← 1 to P do j ← j + id worker Receive C i ,j from worker id worker Algorithm 2: Homogeneous version, worker program.

for all blocks do Receive C i ,j from master

for k ← 1 to t do for j ← (j -1)µ + 1 to j µ do Receive B k,j for i ← (i -1)µ + 1 to i µ do Receive A i,k for j ← (j -1)µ + 1 to j µ do Ci,j ← Ci,j + A i,k .B k,j
Return C i ,j to master they may be overloaded. We can compute the number of processors P so that the time needed to send blocks to P processors will be roughly equal to (or slightly greater than) the time spent by one processor for its computations. P is the smallest integer such that 2µtc × P ≥ µ 2 tw. Indeed, this is the smallest value to saturate the communication capacity of the master while sustaining the corresponding computations. Finally, we cannot use more processors than available. In the context of matrix multiplication, c and w are of the form c = q 2 τc and w = q 3 τa, where τc and τa respectively represent the elementary communication and computation times. Hence:

P = min p, µq 2
τa τc .
For the sake of simplicity, we suppose that r is divisible by µ, and s by Pµ. We allocate µ block columns (i.e., qµ consecutive columns of the original matrix) of C to each processor. The algorithm is made of two parts: Algorithm 1 outlines the program of the master, while Algorithm 2 is the program of each worker.

Impact of the start-up overhead

If we follow the execution of the homogeneous algorithm, we may wonder whether we can really neglect the input/output of C blocks. We decide to sequentialize the sending, computing, and receiving of the C blocks, so that each worker loses 2c time-units per block, i.e., per tw time-units. As there are P ≤ µw 2c + 1 workers, the total loss would be of 2cP time-units every tw time-units, which is less than µ t + 2c tw . For example, with c = 2, w = 4.5, µ = 4 and t = 100, we enroll P = 5 workers, and the total lost is at most 4%, which is small enough to be neglected. Note that it would be technically possible to design an algorithm where the sending of the next block is overlapped with the last computations of the current block, but the whole procedure would become much more complicated.

Algorithms for heterogeneous platforms

We now consider the general problem, i.e., when processors are heterogeneous in term of memory size as well as computation or communication time. As in the previous section, mi is the number of q × q blocks that fit in the memory of worker Pi, and we need to load into the memory of Pi 2µi blocks of A, 2µi blocks of B, and µ 2 i blocks of C. This number of blocks loaded into memory changes from worker to worker, as it depends on their memory capacities: µi is the largest integer such that µ 2 i + 4µi ≤ mi. To adapt our maximum re-use algorithm to heterogeneous platforms, we first present a steady-state-like approach and discuss its limitations. We then introduce our algorithm for heterogeneous platforms.

Bandwidth-centric resource selection

Each worker Pi has parameters ci, wi, and µi, and each participating Pi needs a time 2µitci to receive its blocks and a time tµ 2 i wi to perform its computations. Once again, we neglect I/O for C blocks. Consider the steady-state of a schedule. During one time-unit, Pi receives a certain amount yi of blocks, both of A and B, and computes xi C blocks. We express the constraints, in terms of communication -the master has limited bandwidth-and of computation -a worker has limited computing power and cannot perform more work than it receives. The objective is to maximize the amount of work performed per time-unit. Altogether, we gather the linear program presented in Figure 1. The optimal solution for this system is a bandwidth-centric strategy [START_REF] Banino | Scheduling strategies for master-slave tasking on heterogeneous processor platforms[END_REF]: we sort workers by non-decreasing values of 2c i µ i and we enroll them as long as

2c i µ i w i ≤ 1.
In this way, we can achieve the throughput ρ ≈ i enrolled 1 w i . This solution seems to be close to the optimal. However, the problem is that workers may not have enough memory to execute it! Consider the example described by Table 2. Using the bandwidth-centric strategy, every 8x seconds:

• P1 receives 4x blocks (x µ1 × µ1 chunks) in 4x seconds, and computes 4x blocks in 8x seconds; • P2 receives 4 blocks (1 µ2 × µ2 chunk) in 4x seconds, and computes 4 blocks in 8x seconds. But P1 computes too quickly: during the time x needed to send a block to P2, P1 updates x 2 blocks, which requires at least √ 2x blocks and as many buffers. As x can be arbitrary large, the bandwidth-centric solution cannot always be realized in practice, and we turn to another algorithm, described below. To avoid the previous buffer problems, resource selection will be performed through a step-by-step simulation. However, the steady-state solution is an upper bound on the performance that can be achieved.

Incremental resource selection

The different memory capacities of the workers imply that we assign them chunks of different sizes. This requirement complicates the global partitioning of the C matrix among the workers. To take this into account, while simplifying the implementation, we decide to assign only full matrix column blocks in the algorithm. This is done in a two-phase approach.

In the first phase we pre-compute the allocation of blocks to processors, using a processor selection algorithm we will describe later. We start as if we had a huge matrix of size ∞ × n i=1 µi. Each time Pi is chosen by the processor selection algorithm, it

               MAXIMIZE i xi SUBJECT TO i yici ≤ 1 ∀i xiwi ≤ 1 ∀i xi µ 2 i ≤ yi 2µi
Table 1. Linear program maximizing the amount of work performed per time-unit.

P1 P2 ci 1 x wi 2 2x µi 2 2 2c i µ i w i 1 2 1 2
Table 2. Platform for which the bandwidth centric solution is not feasible.

is assigned a square chunk of µ 2 i C blocks. As soon as Pi has enough blocks to fill up µi block columns of the initial matrix, we decide that Pi will indeed execute these columns during the parallel execution. Therefore we maintain a panel of p i=1 µi block columns and fill them out by assigning blocks to processors. We stop this phase as soon as all the r × s blocks of the initial matrix have been allocated columnwise by this process. Pi will be assigned a block column after it has been selected r µ i times by the algorithm.

In the second phase we perform the actual execution. Messages will be sent to workers according to the previous selection process. The first time a processor Pi is selected, it receives a square chunk of µ 2 i C blocks, which initializes its repeated pattern of operation: the following t times, Pi receives µi A and µi B blocks, which requires 2µici time-units.

There remains to decide which processor to select at each step of the first phase. We have no closed-form formula for the allocation of blocks to processors. Instead, we use an incremental algorithm to compute which worker the next blocks will be assigned to. We have two variants of the incremental algorithm, a global one that aims at optimizing the overall communication-to-computation ratio, and a local one that selects the best processor for the next stage.

Global selection algorithm. The intuitive idea, here, is to select the processor that maximizes the ratio of the total work achieved so far (in terms of block updates) over the completion time of the last communication. The latter represents the time spent by the master so far, either sending data to workers or staying idle waiting for the workers to finish their current computations. We have:

ratio ← total work achieved completion time of last communication

Estimating computations is easy: Pi executes µ 2 i block updates per assignment. Communications are slightly more complicated to deal with; we cannot just use the communication time 2µici of Pi for the A and B blocks because we need to take ready times into account: if Pi is currently busy executing work, it cannot receive additional data too much in advance as its memory is limited.

Local selection algorithm. The global selection algorithm picks, as the next processor, the one that maximizes the ratio of the total amount of work assigned over the time needed to send all the required data. On the other hand, the local selection algorithm chooses, as destination of the i-th communication, the processor that maximizes the ratio of the amount of work assigned by this communication over the time during which the communication link is used to performed this communication (i.e., the elapsed time between the end of (i -1)-th communication and the end of the i-th communication). As previously, if processor Pj is the target of the i-th communication, the i-th communication is the sending of µj blocks of A and µj blocks of B to processor Pj, which enables it to perform µ 2 j updates.

Variants. At this point, we can have a global or a local selection process. Each process can either make a decision only looking at the next communication, or in a look-ahead approach. Another variant would be to take into account the cost of initially sending µ 2 i blocks of matrix C to processor Pi the first time its receives blocks. Overall, we thus define eight different selection algorithms (global or local, look-ahead or not, µ 2 i C costs or not). There is no reason for one of these heuristics to always dominate the others. We will thus consider the eight of them in our experiments.

MPI experiments

In this section, we aim at validating the previous theoretical results and algorithms. We conduct a variety of MPI experiments to compare our new schemes with several algorithms from the literature. We target heterogeneous platforms, and we assess the impact of the degree of heterogeneity (in processor speed, link bandwidth, and memory capacity) on the performance of the various algorithms.

The code and the experimental results can be downloaded from: http://graal.ens-lyon.fr/~jfpineau/Downloads/ppopp.tgz.

Platforms

We used a heterogeneous cluster composed of twenty-seven processors located in Lyon. It is composed of four different homogeneous sets of machines. The different sets are composed of: 1) 8 SuperMicro servers 5013-GM, with processors P4 2.4 GHz; 2) 5 SuperMicro servers 6013PI, with processors P4 Xeon 2.4 GHz; 3) 7 SuperMicro servers 5013SI, with processors P4 Xeon 2.6 GHz; 4) 7 SuperMicro servers IDE250W, with processors P4 2.8 GHz.

All nodes have 1 GB of memory and are running the Linux operating system. The nodes are connected with a switched 10 Mbps Fast Ethernet network. As this platform may not be as heterogeneous as we would like, we sometimes artificially modify its heterogeneity. In order to artificially slow down a communication link, we send the same message several times to one worker. The same idea works for processor speeds: we ask a worker to compute a given matrix-product several times in order to slow down its computation capability. In all experiments, except the last batch, we used nine processors: one master and eight workers.

Algorithms

We choose four different algorithms from the literature which we compare our algorithms to. The closest work addressing our problem is Toledo's out-of-core algorithm [START_REF] Toledo | A survey of out-of-core algorithms in numerical linear algebra[END_REF]. Hence, this work will serve as the baseline reference. Then we will study hybrid algorithms, i.e., algorithms which use our memory layout and are based on classical principles such as round-robin, min-min [START_REF] Maheswaran | Dynamic matching and scheduling of a class of independent tasks onto heterogeneous computing systems[END_REF], or a dynamic demand-driven approach. The first six algorithms below use our memory allocation, the only difference between them is the order in which the master sends blocks to workers. Homogeneous algorithm (Hom) is our homogeneous algorithm.

It makes resource selection and sends blocks to the selected workers in a round-robin fashion. When run on a heterogeneous platform, it tries to build a very simple homogeneous platform.

As the algorithm's only constraint is to send same size blocks to all participating workers, for a given memory size, we consider the homogeneous virtual platform composed of those workers having at least that amount of memory, and we estimate the total execution time of our homogeneous algorithm, for the targeted matrix-product, on that virtual platform (the apparent processor speed is the minimum of the processor speeds, the apparent communication bandwidth is the minimum of the communication bandwidths). We do this process for all the different memory sizes present in the actual platform, and we pick the virtual platform that minimizes the total estimated execution time.

Homogeneous algorithm improved (HomI) is our homogeneous algorithm running on a more carefully chosen homogeneous platform. For each memory size, communication speed, and computation speed present in an heterogeneous platform, we consider the homogeneous virtual platform composed of those workers having at least that performance. Then, we compute the total execution time of our homogeneous algorithm, for the targeted matrix-product, on that virtual platform (the apparent processor speed is the considered processor speed, and so on). We do this process for all the existing values, and we pick the virtual platform that minimizes the total execution time. Heterogeneous algorithm (Het) is our heterogeneous algorithm.

As we can have eight different versions of the resource selection, in a first step we simulate the eight versions, and then we pick and run the best one. Overlapped Round-Robin, Optimized Memory Layout (OR-ROML) sends tasks to all available workers in a round-robin fashion. It does not make any resource selection. Overlapped Min-Min, Optimized Memory Layout (OMMO-ML) is a static scheduling heuristic, which sends the next block to the first worker that will finish it. As it is looking for potential workers in a given order, this algorithm performs some resource selection too. Theoretically, as our homogeneous resource selection ensures that the first worker is free to compute when we finish sending blocks to the others, OMMOML and Hom should have a similar behavior on homogeneous platforms. Overlapped Demand-Driven, Optimized Memory Layout (O-DDOML) is a demand-driven algorithm. In our memory layout, two buffers of size µi are reserved for matrix A, and two for matrix B. In order to use the two available extra buffers (the second for A and the second for B), one sends the next block to the first worker which can receive it. This would be a dynamic version of our algorithm, if it took worker selection into account. Block Matrix Multiply (BMM) is Toledo's algorithm [START_REF] Toledo | A survey of out-of-core algorithms in numerical linear algebra[END_REF]. It splits each worker memory equally into three parts and allocates one slot for a square block of A, another for a square block of B, and the last one for a square block of C, with the square blocks having the same size. It sends blocks to the workers in a demand-driven fashion. First a worker receives a block of C, then it receives corresponding blocks of A and B in order to update C, until C is fully computed. Note that the six algorithms using our optimized memory layout are considering matrices as composed of square blocks of size q × q = 80 × 80, while BMM loads three panels, each of size one third of the available memory, for A, B and C.

When launching an algorithm on the platform, the very first step we do is to determine the platform's parameters. For that purpose, we launch a benchmark on it, in order to get the memory size, the communication speed, and the computation speed. The different speeds are determined by sending and computing a square block of size q × q ten times on each worker, and computing the median of the times obtained. This step takes between 20 and 80 seconds, depending on the speed of the workers, and is made before each algorithm, even ORROML, ODDOML, and BMM, which only need the memory size. This step represents at most 2% of the total time of execution.

In the following section, the times given takes into account the decision process of the algorithms, i.e., the simulation of the eight different versions of the resource selection for Het, the construction of an homogeneous platform for Hom and HomI, etc.

Experimental results

In the first three sets of experiments, we only have one parameter of heterogeneity, either the amount of memory, the communication speed, or the computation speed. We test the algorithms on such platforms with five matrices of increasing sizes. As we do not want to change several parameters at a time, we only change the value of parameter s (rather than, for instance, always consider square matrices). Matrix A is of size 8000 × 8000 whereas B is of increasing sizes 8000 × 64000, 8000 × 80000, 8000 × 96000, 8000 × 112000, and 8000 × 128000. For all other experiments, A is of size 8000 × 8000 and B is of size 8000 × 80000. The heterogeneous workers have different memory capacities, which implies that each algorithm, even BMM, assigns them chunks of different sizes. In order to simplify the global partitioning of matrix C, we decide to only assign workers full matrix column blocks.

As we want to assess whether the performance of any studied algorithm depends on the matrix size, we look at the relative cost of the algorithms rather than at their absolute execution times. The relative cost of a given algorithm on a particular instance is equal to the execution time, or makespan, achieved on that instance by the algorithm, divided by the minimum makespan achieved on that instance. Using relative cost also enables us to build statistics on the performance of algorithms.

Besides relative cost, we take into account the number of processors used. To assess the efficiency of a given algorithm, we look at its relative work, which is equal, for a given instance, to its makespan times the number of enrolled processors, divided by the minimum of this value over all studied algorithms.

Heterogeneous memory size

Here we assess the performance of our algorithms with respect to memory heterogeneity. We launch the algorithms on a homogeneous platform in terms of communication and computation capabilities, but where workers have different memory capacities. We suppose that two workers only have 256 MB of memory, four of them 512 MB, and the last two 1024 MB. 4(a) presents the relative cost of the algorithms, whose general shape is very similar for all five matrix sizes. ODDOML and our heterogeneous algorithm Het have the best makespans. At the other end of the spectrum, OMMOML is twice as bad. In between, Hom, HomI, ORROML, and BMM are roughly twenty percent slower. To give an idea of execution times, Het needs about 2000 seconds to compute the product of the smallest matrices, and about 3500 seconds for the largest.

The variations in the cost of BMM can easily be justified. The memory layout used in BMM is different from the other algorithms. Therefore the size of the matrix chunks used by BMM are different. The matrices are rather small (to be able to evaluate a significant number of algorithms on a significant number of platforms). Hence we observe some non negligible side effects (matrix size divided by chunk size not being a multiple of the number of processors used). Therefore, for a given platform, some memory sizes are more favorable for BMM than for the algorithms using our memory layout. We will see throughout our experiments on heterogeneous platforms that even if sometimes these side effects help BMM achieve reasonable cost, it may also have a very bad cost in some other situations..

The ranking of the algorithms is quite different when we look at the relative work (Figure 4(b)). First we have OMMOML, which only uses two workers, and is thus very thrifty, at the expense of its absolute cost. Then we have HomI, Het, Hom, and ODDOML. Hom relative work is always better than that of ODDOML: Hom is performing some resource selection contrarily to ODDOML, which always uses all the processors. And HomI is even better than Hom, thanks to a better platform selection. This is also the reason why the gap between the relative work of Het and that of ODDOML is significantly larger than the gap between their relative costs. Finally, we have ORROML and BMM, which do not achieve an efficient makespan and make no resource selection, and thus achieve very bad relative work.

Heterogeneous communication links

We now assess the performance of our algorithms when communication links have heterogeneous capabilities. The target platform is composed of two workers with a 10Mbps communication link, four workers with a 5Mbps communication link, and the last two have a 1Mbps communication link.

Figure 5(a) shows the relative cost. The superiority of our heterogeneous algorithm over BMM is clear.

Het, HomI, and OMMOML have excellent makespans, and make a good resource selection, as seen on Figure 5(b)). The first figure also shows the gap between HomI and Hom: Hom performs close to ODDOML, while HomI achieves a close to best makespan. This figure underlines the importance of carefully choosing the processors on which launching the algorithm: Hom only uses two processors because of the platform parameters and the way it extracts a homogeneous platform. BMM has the worst makespan and makes no resource selection, which explains its worse relative work. BMM achieves a makespan that is 70 to 90 percent worse than the best one. Concerning execution times, Het needs about 2500 seconds to compute the product of the smallest matrices, and about 5000 seconds for the largest.

Heterogeneous computations

Here we assess the performance of our algorithms when computation capabilities are heterogeneous. Workers have homogeneous communications and memory capacities, but different computation speeds. The platform is composed of two fast workers of speed S, four workers of speed S/2 and two workers of speed S/4.

In Figure 6(a), we see the relative cost obtained during this set of experiments. BMM performs rather well, but its makespan is larger than that of Het. Moreover, looking at the relative work in Figure 6(b), we see that the gap between the algorithms becomes larger, as our algorithms enroll fewer resources during execution. Among the other algorithms, we see that ODDOML performs well. If we look at the relative work, we also see that Het uses more and more processors as matrix size increases. During these experiments, Het needs about 2000 seconds to compute the product of the smallest matrices, and about 4000 seconds for the largest.

Fully heterogeneous platforms

We now consider fully heterogeneous platforms. Communication links, computation capabilities, and memory capacities can take two different values, which leads to eight possibilities, one per worker. We build that way two different platforms by fixing the ratio between the small and large values for each characteristics to either 2 in the first setting or 4 in the second one (first two columns in Figures 7(a We see that Het achieves the best makespan for all but two of the 12 platforms, and in the remaining cases is no more than 9% and 2% away from the best studied algorithm. All the other algorithms are, at least once, more than 41% away from the best solution. For example, ORROML can be up to 88% worse than the best achieved makespan. Only ODDOML achieves reasonable makespans on average but, as it does not select resources, its relative work is far worse. The relative work of Het is the best among all algorithms except HomI, and the unusable OMMOML, whose makespan can be 215% away from the best solution. But even if our improved homogeneous algorithm performs a good resource selection, the makespan it achieves can be up to 80% larger than the best makespan, and is 34% larger on average. According to the experiments, Het needs between 2700 seconds and 6000 seconds.

Real platform

In this set of experiments, we use almost all the processors of our platform. We do not modify the communication speed nor the com-putation speed of the workers. We take five processors of each of the four sets of machines, which gives us a rather homogeneous platform. We either use this platform "as is" (August 2007 configuration of Figure 8(a)) or we limit the amount of memory available on each processor to its value before the last memory upgrade (November 2006 configuration of Figure 8(b)). The actual platform was then:

• 5 SuperMicro servers 5013-GM, with processors P4 2.4 GHz with 256 MB of memory; • 5 SuperMicro servers 6013PI, with processors P4 Xeon 2.4

GHz with 1 GB of memory; • 5 SuperMicro servers 5013SI, with processors P4 Xeon 2.6

GHz with 1 GB of memory; • 5 SuperMicro servers IDE250W, with processors P4 2.8 GHz with 256 MB of memory. We use an extra processor as the master. The matrix is of size 8000 × 8000 for A and 8000 × 320000 for B. The results of these experiments are summarized in Figure 8. The results on the actual platform are similar to those obtained on homogeneous platforms [START_REF] Pineau | Revisiting matrix product on master-worker platforms[END_REF]. All the algorithms but BMM have similar makespan. All algorithms making resource selection use eleven workers among the twenty available, which explains why they achieve similar relative work.

The results of the experiments on the older version of the platform are very similar to the ones previously obtained on memory heterogeneous platforms. ODDOML and our heterogeneous algorithm Het achieve the best makespans. Then we have OR-ROML, our homogeneous algorithms, BMM, and finally OM-MOML which is 60% worse than Het. The execution time is around 7800 seconds for Het. If we look at resource selection, Het uses only the ten workers which have 1 GB of memory, and achieves a makespan close to ODDOML's, which uses the whole platform. On another side, Hom, HomI, OMMOML use six workers with small memories. All other algorithms use the whole platform. We can thus see the impact of the resource selection on the relative work of the algorithms. respectively present the relative cost and the relative work obtained for each experiment by our heterogeneous algorithm (Het), by Toledo's algorithm (BMM), and by the best of the dynamic heuristics using our memory layout (ODDOML). The results show the superiority of our memory allocation. Furthermore, if we add the resource selection of Het, not only do we achieve, most of the time, the best makespan, but also the best relative work as we also spare resources. Using our memory layout (ODDOML) rather than Toledo's (BMM) enables us to gain 19% of execution time on average. When this is combined with resource selection, this enables us to gain additionally 10%, which is 27% against Toledo's running time. We achieve this significant gain while sparing resources. Our Het algorithm is on average 1% away from the best achieved makespan. At worst Het is 14% away from the best makespan, ODDOML 61%, and BMM 128%. Moreover, we have seen that 80% of the time, the performance of Het was in fact obtained thanks to a global resource selection.

Summary

The steady-state approach described in Section 5 gives us an upper-bound on the best achievable throughput. This upper-bound is very optimistic as it assumes unbounded memories and does not take into account the communication costs due to the elements of matrix C. This upper bound is nevertheless on average only 2.29 times greater than the throughput achieved by Het (and at worst is 3.42 times greater). Therefore, considering this upper-bound tells us that our Het algorithm not only has good relative cost when compared to the other algorithms, but also has very good absolute cost.

Altogether, we have thus been able to design an efficient, thrifty, and reliable algorithm.

Related work

As already mentioned, the design of parallel algorithms for limited memory processors is very similar to the design of out-of-core routines for classical parallel machines. On the theoretical side, Hong and Kung [START_REF] Hong | I/O complexity: the red-blue pebble game[END_REF] investigate the I/O complexity of several computational kernels in their pioneering paper. Toledo [START_REF] Toledo | A survey of out-of-core algorithms in numerical linear algebra[END_REF] proposes a nice survey on the design of out-of-core algorithms for linear algebra, including dense and sparse computations. We refer to [START_REF] Toledo | A survey of out-of-core algorithms in numerical linear algebra[END_REF] for a complete list of implementations. The design principles followed by most implementations are introduced and analyzed by Dongarra et al. [START_REF] Dongarra | Key concepts for parallel out-of-core LU factorization[END_REF].

A similar thread of work, although in a different context, deals with reconfigurable architectures, either pipelined bus systems [START_REF] Li | Parallel matrix multiplication on a linear array with a reconfigurable pipelined bus system[END_REF], or FPGAs [START_REF] Zhuo | Scalable and modular algorithms for floating-point matrix multiplication on reconfigurable computing systems[END_REF]. In the latter approach, tradeoffs must be found to optimize the size of the on-chip memory and the available memory bandwidth, leading to partitioned algorithms that re-use data intensively.

Please see the companion research report of this paper [START_REF] Pineau | Revisiting matrix product on master-worker platforms[END_REF] for an overview or relevant literature dealing with (i) load balancing techniques on heterogeneous platforms; (ii) linear algebra algorithms on heterogeneous clusters; (iii) models for heterogeneous platforms; and (iv) master-worker scheduling on grid platforms. This research report also subsumes our previous work focusing on homogeneous platforms [START_REF] Pineau | Revisiting matrix product on master-worker platforms[END_REF].

Conclusion

The main contributions of this paper are the following: 1. On the theoretical side, we have derived a new, tighter, bound on the minimal volume of communications needed to multiply two matrices. From this lower bound, we have defined an efficient memory layout, i.e., an algorithm to share the memory available on the workers among the three matrices. 2. On the practical side, starting from our memory layout, we have designed an algorithm for homogeneous platforms whose performance is quite close to the communication volume lower bound. We have extended this algorithm to deal with heterogeneous platforms. Please refer to the companion research report [START_REF] Pineau | Revisiting matrix product on master-worker platforms[END_REF] for a discussion on how to adapt the approach for LU factorization. 3. Through MPI experiments, we have shown that our algorithm for heterogeneous platforms has far better performance than solutions using the memory layout proposed in [START_REF] Toledo | A survey of out-of-core algorithms in numerical linear algebra[END_REF]. Furthermore, this static heterogeneous algorithm has slightly better performance than dynamic algorithms using the same memory layout, but uses fewer processors, and has a far better worst case. It is therefore a very good candidate for deploying applications on heterogeneous platforms. It would be interesting to assess whether our memory layout could prove useful in the context of out-of-core algorithms for heterogeneous platforms.

Figure 1 .

 1 Figure 1. Partitioning of the three matrices A, B, and C.

2 3 m 1 8m

 31 , and the lower bound for the communication-to-computation ratio is: CCRopt ≥ 27 8m . This bound improves upon the bestknown value derived by Ironya, Toledo, and Tiskin[START_REF] Ironya | Communication lower bounds for distributed-memory matrix multiplication[END_REF].

Figure 2 .

 2 Figure 2. Memory layout for the maximum re-use algorithm when m = 21: µ = 4; 1 block is used for A, µ for B, and µ 2 for C.

Figure 3 .

 3 Figure 3. Four steps of the maximum re-use algorithm, with m = 21 and µ = 4. Updated elements of C are written white on black.

Algorithm 1 : 4 + m - 2 P

 142 Homogeneous version, master program. µ ← √ ← min p, µw 2c

  (b) Efficiency of algorithms.

Figure 4 .

 4 Figure 4. Heterogeneous memory.

Figure

  Figure4(a) presents the relative cost of the algorithms, whose general shape is very similar for all five matrix sizes. ODDOML and our heterogeneous algorithm Het have the best makespans. At the other end of the spectrum, OMMOML is twice as bad. In between, Hom, HomI, ORROML, and BMM are roughly twenty percent slower. To give an idea of execution times, Het needs about 2000 seconds to compute the product of the smallest matrices, and about 3500 seconds for the largest.The variations in the cost of BMM can easily be justified. The memory layout used in BMM is different from the other algorithms. Therefore the size of the matrix chunks used by BMM are different. The matrices are rather small (to be able to evaluate a significant number of algorithms on a significant number of platforms). Hence we observe some non negligible side effects (matrix size divided by chunk size not being a multiple of the number of processors used). Therefore, for a given platform, some memory sizes are more favorable for BMM than for the algorithms using our memory layout. We will see throughout our experiments on heterogeneous platforms that even if sometimes these side effects help BMM achieve reasonable cost, it may also have a very bad cost in some other situations..The ranking of the algorithms is quite different when we look at the relative work (Figure4(b)). First we have OMMOML, which only uses two workers, and is thus very thrifty, at the expense of its absolute cost. Then we have HomI, Het, Hom, and ODDOML. Hom relative work is always better than that of ODDOML: Hom is performing some resource selection contrarily to ODDOML, which always uses all the processors. And HomI is even better than Hom, thanks to a better platform selection. This is also the reason why the gap between the relative work of Het and that of ODDOML is significantly larger than the gap between their relative costs. Finally, we have ORROML and BMM, which do not achieve an efficient makespan and make no resource selection, and thus achieve very bad relative work.

  (a) Cost of algorithms. (b) Efficiency of algorithms.

Figure 5 .

 5 Figure 5. Heterogeneous communication links.

  (a) Cost of algorithms. (b) Efficiency of algorithms.

Figure 6 .

 6 Figure 6. Heterogeneous computations.

  ) and 7(b)). In order to show that our heterogeneous algorithm works on any heterogeneous platform, we also randomly create ten different platforms (last ten columns on the same figures). The ratio between minimum and maximum values of communication links, computation capacities, and memory size is up to four. Matrix A is of size 8000 × 8000 and B of size 8000 × 80000.

  (a) Cost of algorithms.

  Efficiency of algorithms.

Figure 7 .

 7 Figure 7. Fully heterogeneous platforms. The results of these experiments are summarized on Figures 7(a) and 7(b).We see that Het achieves the best makespan for all but two of the 12 platforms, and in the remaining cases is no more than 9% and 2% away from the best studied algorithm. All the other algorithms are, at least once, more than 41% away from the best solution. For example, ORROML can be up to 88% worse than the best achieved makespan. Only ODDOML achieves reasonable makespans on average but, as it does not select resources, its relative work is far worse. The relative work of Het is the best among all algorithms except HomI, and the unusable OMMOML, whose makespan can be 215% away from the best solution. But even if our improved homogeneous algorithm performs a good resource selection, the makespan it achieves can be up to 80% larger than the best makespan, and is 34% larger on average. According to the experiments, Het needs between 2700 seconds and 6000 seconds.

Figure 8 .

 8 Figure 8. Real platform.

Figure 9

 9 Figure 9 summarizes all our MPI experiments. Figures 9(a) and 9(b)respectively present the relative cost and the relative work obtained for each experiment by our heterogeneous algorithm (Het), by Toledo's algorithm (BMM), and by the best of the dynamic heuristics using our memory layout (ODDOML). The results show the superiority of our memory allocation. Furthermore, if we add the resource selection of Het, not only do we achieve, most of the time, the best makespan, but also the best relative work as we also spare resources. Using our memory layout (ODDOML) rather than Toledo's (BMM) enables us to gain 19% of execution time on average. When this is combined with resource selection, this enables us to gain additionally 10%, which is 27% against Toledo's running time. We achieve this significant gain while sparing resources. Our Het algorithm is on average 1% away from the best achieved makespan. At worst Het is 14% away from the best makespan, ODDOML 61%, and BMM 128%. Moreover, we have seen that 80% of the time, the performance of Het was in fact obtained thanks to a global resource selection.The steady-state approach described in Section 5 gives us an upper-bound on the best achievable throughput. This upper-bound is very optimistic as it assumes unbounded memories and does not take into account the communication costs due to the elements of matrix C. This upper bound is nevertheless on average only 2.29 times greater than the throughput achieved by Het (and at worst is 3.42 times greater). Therefore, considering this upper-bound tells us that our Het algorithm not only has good relative cost when compared to the other algorithms, but also has very good absolute cost.Altogether, we have thus been able to design an efficient, thrifty, and reliable algorithm.

Figure 9 .

 9 Figure 9. Summary of experiments.

Of course, there are Θ(n 3 ) computations if we only consider algorithms that uses the standard way of multiplying matrices; e.g., this excludes Strassen's and Winograd's algorithms.


Jack Dongarra

Jean-Franc ¸ois Pineau

Yves Robert

Zhiao Shi
email: shi@cs.utk.edu

Frédéric Vivien
email: frederic.vivien@ens-lyon.frdongarra

Revisiting Matrix Product on Master-Worker Platforms

This paper is aimed at designing efficient parallel matrix-product algorithms for homogeneous masterworker platforms. While matrix-product is well-understood for homogeneous 2D-arrays of processors (e.g., Cannon algorithm and ScaLAPACK outer product algorithm), there are two key hypotheses that render our work original and innovative: -Centralized data. We assume that all matrix files originate from, and must be returned to, the master. The master distributes both data and computations to the workers (while in ScaLAPACK, input and output matrices are initially distributed among participating resources). Typically, our approach is useful in the context of speeding up MATLAB or SCILAB clients running on a server (which acts as the master and initial repository of files).

-Limited memory. Because we investigate the parallelization of large problems, we cannot assume that full matrix panels can be stored in the worker memories and re-used for subsequent updates (as in ScaLAPACK). The amount of memory available in each worker is expressed as a given number of buffers, where a buffer can store a square block of matrix elements. These square blocks are chosen so as to harness the power of Level 3 BLAS routines; they are of size 80 or 100 on most platforms.

We have devised efficient algorithms for resource selection (deciding which workers to enroll) and communication ordering (both for input and result messages), and we report a set of MPI experiments conducted on a platform at the University of Tennessee.

Introduction

Matrix product is a key computational kernel in many scientific applications, and it has been extensively studied on parallel architectures. Two well-known parallel versions are Cannon's algorithm [START_REF] Cannon | A cellular computer to implement the Kalman filter algorithm[END_REF] and the ScaLAPACK outer product algorithm [START_REF] Blackford | ScaLAPACK Users' Guide[END_REF]. Typically, parallel implementations work well on 2D processor grids, because the input matrices are sliced horizontally and vertically into square blocks that are mapped one-to-one onto the physical resources; several communications can take place in parallel, both horizontally and vertically. Even better, most of these communications can be overlapped with (independent) computations. All these characteristics render the matrix product kernel quite amenable to an efficient parallel implementation on 2D processor grids.

However, current architectures typically take the form of clusters, which are composed of computing resources interconnected by a sparse network: there are no direct links between any pair of processors. Instead, messages from one processor to another are routed via several links, likely to have different capacities. Worse, congestion will occur when two messages, involving two different sender/receiver pairs, collide because a same physical link happens to belong to the two routing paths. Therefore, an accurate estimation of the communication cost requires a precise knowledge of the underlying target platform. In addition, it becomes necessary to include the cost of both the initial distribution of the matrices to the processors and of collecting back the results. These input/output operations have always been neglected in the analysis of the conventional algorithms. This is because only O(n 2 ) coefficients need to be distributed in the beginning, and gathered at the end, as opposed to the O(n 3 ) computations to be performed (where n is the problem size). The assumption that these communications can be ignored could have made sense on dedicated processor grids like, say, the Intel Paragon, but it is no longer reasonable on networks of workstations.

In this paper, we do not try to adapt the 2D processor grid strategy to networks of workstations. Instead, we adopt a realistic application scenario, where input files are read from a fixed repository (disk on a data server). Computations will be delegated to available resources in the target architecture, and results will be returned to the repository. This calls for a master-worker paradigm, or more precisely for a computational scheme where the master (the processor holding the input data) assigns computations to other resources, the workers. In this centralized approach, all matrix files originate from, and must be returned to, the master. The master distributes both data and computations to the workers (while in ScaLAPACK, input and output matrices are supposed to be equally distributed among participating resources beforehand). Typically, our approach is useful in the context of speeding up MATLAB or SCILAB clients running on a server (which acts as the master and initial repository of files).

Because we investigate the parallelization of large problems, we cannot assume that full matrix panels can be stored in worker memories and re-used for subsequent updates (as in ScaLAPACK). The amount of memory available in each worker is expressed as a given number m of buffers, where a buffer can store a square block of matrix elements. The size q of these square blocks is chosen so as to harness the power of Level 3 BLAS routines: q = 80 or 100 on most platforms.

To summarize, the target platform is composed of several workers with limited memory capacities. The first problem is resource selection. How many workers should be enrolled in the execution? All of them, or maybe only a fraction? Once participating resources have been selected, there remain several scheduling decisions to take: how to minimize the number of communications? in which order workers should receive input data and return results? what amount of communications can be overlapped with (independent) computations? The goal of this paper is to design efficient algorithms for resource selection and communication ordering. In addition, we report MPI experiments on platforms at the University of Tennessee.

The rest of the paper is organized as follows. In Section 2, we state the scheduling problem precisely, and we introduce some notations. Next, in Section 3, we proceed with the analysis of the total communication volume that is needed in the presence of memory constraints, and we improve a well-known bound by Toledo [START_REF] Toledo | A survey of out-of-core algorithms in numerical linear algebra[END_REF][START_REF] Ironya | Communication lower bounds for distributed-memory matrix multiplication[END_REF]. In Section 4, we propose a scheduling algorithm that includes resource selection. We report several MPI experiments in Section 5. Finally, we state some concluding remarks in Section 6. Due to lack of space, related work is not discussed in this paper: please refer to the extended version [START_REF] Dongarra | Revisiting matrix product on master-worker platforms[END_REF] for an overview of relevant literature.

Framework Application

We deal with the computational kernel C ← C + A × B. We partition the three matrices A, B, and C as illustrated in Figure 1. More precisely:

• We use a block-oriented approach. The atomic elements that we manipulate are not matrix coefficients but instead square blocks of size q × q (hence with q 2 coefficients). This is to harness the power of Level 3 BLAS routines [START_REF] Blackford | ScaLAPACK Users' Guide[END_REF]. Typically, q = 80 or 100 when using ATLAS-generated routines [START_REF] Whaley | Automatically tuned linear algebra software[END_REF]. • The input matrix A is of size n A × n AB :

-we split A into r horizontal stripes A i , 1 ≤ i ≤ r, where r = n A /q; -we split each stripe A i into t square q×q blocks A i,k , 1 ≤ k ≤ t, where t = n AB /q. • The input matrix B is of size n AB × n B :

-we split B into s vertical stripes B j , 1 ≤ j ≤ s, where s = n B /q; -we split stripe B j into t square q × q blocks B k,j ,

1 ≤ k ≤ t. • We compute C = C + A × B. Matrix C is accessed
(both for input and output) by square q × q blocks C i,j , 1 ≤ i ≤ r, 1 ≤ j ≤ s. There are r × s such blocks. We point out that with such a decomposition all stripes and blocks have same size. This will greatly simplify the analysis of communication costs.

Platform

We target a star network S = {P 0 , P 1 , P 2 , . . . , P p }, composed of a master P 0 and of p identical workers P i , 1 ≤ i ≤ p. Because we manipulate large data blocks, we adopt a lin-ear cost model, both for computations and communications (i.e., we neglect start-up overheads). We have the following notations:

• It takes X.w time-units to execute a task of size X on P i ; • It takes X.c time units for the master P 0 to send a message of size X to P i or to receive a message of size X from P i . Without loss of generality, we assume that the master has no processing capability (otherwise, add a fictitious extra worker paying no communication cost to simulate computation at the master).

Next, we need to define the communication model. We adopt the one-port model [START_REF] Bhat | Efficient collective communication in distributed heterogeneous systems[END_REF][START_REF] Bhat | Efficient collective communication in distributed heterogeneous systems[END_REF], which is defined as follows: (i) the master can only send data to, and receive data from, a single worker at a given time-step; (ii) a given worker cannot start execution before it has terminated the reception of the message from the master; similarly, it cannot start sending the results back to the master before finishing the computation. In fact, this one-port model naturally comes in two flavors, depending upon whether we allow the master to simultaneously send and receive messages or not. If we do allow for simultaneous sends and receives, we have actually the two-port model. Here we concentrate on the true one-port model, where the master cannot be enrolled in more than one communication at any time-step.

The one-port model is realistic. Bhat, Raghavendra, and Prasanna [START_REF] Bhat | Efficient collective communication in distributed heterogeneous systems[END_REF] advocate its use because "current hardware and software do not easily enable multiple messages to be transmitted simultaneously." Even if non-blocking multithreaded communication libraries allow for initiating multiple send and receive operations, they claim that all these operations "are eventually serialized by the single hardware port to the network." Experimental evidence of this fact has recently been reported by Saif and Parashar [START_REF] Saif | Understanding the behavior and performance of non-blocking communications in MPI[END_REF], who report that asynchronous MPI sends get serialized as soon as message sizes exceed a hundred kilobytes. Their result hold for two popular MPI implementations, MPICH on Linux clusters and IBM MPI on the SP2. Note that all the MPI experiments in Section 5 obey the one-port model.

Our final assumption is related to memory capacity; we assume that a worker P i can only store m blocks (either from A, B, or C). For large problems, this memory limitation will considerably impact the design of the algorithms, as data re-use will be greatly dependent on the amount of available buffers.

Minimization of the communication volume

In this section, we derive a lower bound on the total number of communications (sent from, or received by, the master) that are needed to execute any matrix multiplication al-gorithm. Since we aim at minimizing the total communication volume, we can simulate any parallel algorithm on a single worker. Therefore, we only need to consider the one-worker case. We deal with the original, and realistic, formulation of the problem as follows:

• The master sends blocks A ik , B kj , and C ij ,

• The master retrieves final values of blocks C ij , and • We enforce limited memory on the worker; only m buffers are available, which means that at most m blocks of A, B, and/or C can simultaneously be stored on the worker. First, we describe an algorithm that aims at re-using C blocks as much as possible after they have been loaded. Next, we assess the performance of this algorithm. Finally, we improve a lower bound previously established by Toledo [START_REF] Toledo | A survey of out-of-core algorithms in numerical linear algebra[END_REF][START_REF] Ironya | Communication lower bounds for distributed-memory matrix multiplication[END_REF].

The maximum re-use algorithm

Below we introduce and analyze the performance of the maximum re-use algorithm, whose memory management is illustrated in Figure 2. Four consecutive execution steps are shown in Figure 3. Assume that there are m available buffers. First we find µ as the largest integer such that 1 + µ+µ 2 ≤ m. The idea is to use one buffer to store A blocks, µ buffers to store B blocks, and µ 2 buffers to store C blocks. In the outer loop of the algorithm, a µ×µ square of C blocks is loaded. Once these µ 2 blocks have been loaded, they are repeatedly updated in the inner loop of the algorithm until their final value is computed. Then the blocks are returned to the master, and µ 2 new C blocks are sent by the master and stored by the worker. As illustrated in Figure 2, we need µ buffers to store a row of B blocks, but only one buffer for A blocks: A blocks are sent in sequence, each of them is used in combination with a row of µ B blocks to update the corresponding row of C blocks. This leads to the following sketch of the algorithm: Outer loop: while there remain C blocks to be computed

• Store µ 2 blocks of C in worker's memory: send a µ × µ square {C i,j / i 0 ≤ i < i 0 + µ, j 0 ≤ j < j 0 + µ} • Inner loop: For each k from 1 to t: 1. Send a row of µ elements {B k,j / j 0 ≤ j < j 0 + µ}; 2. Sequentially send µ elements of column {A i,k / i 0 ≤ i < i 0 + µ}. For each A i,k , update µ elements of C
• Return results to master.

Performance and lower bound

The performance of one iteration of the outer loop of the maximum re-use algorithm can readily be determined. We For each value of t, we need µ elements of A and µ elements of B, and we update µ 2 blocks. In terms of block operations, the communication-to-computation ratio achieved by the algorithm is thus

C C C C C C C C C C C C C C C C A B B B B 1 µ µ 2
CCR = 2µ 2 + 2µt µ 2 t = 2 t + 2 µ .
For large problems, i.e., large values of t, we see that CCR is asymptotically close to the value CCR ∞ = 2 √ m . We point out that, in terms of data elements, the communication-tocomputation ratio is divided by a factor q. Indeed, a block consists of q 2 coefficients but an update requires q 3 floatingpoint operations.

How can we assess the performance of the maximum reuse algorithm? How good is the value of CCR? To see this, we refine an analysis due to Toledo [START_REF] Toledo | A survey of out-of-core algorithms in numerical linear algebra[END_REF]. The idea is to estimate the number of computations made thanks to m consecutive communication steps (again, the unit is a matrix block here). We need some notations:

• We let α old , β old , and γ old be the number of buffers dedicated to A, B, and C at the beginning of the m communication steps; • We let α recv , β recv , and γ recv be the number of A, B, and C blocks sent by the master during the m communication steps; • Finally, we let γ send be the number of C blocks returned to the master during these m steps. Obviously, the following equations must hold true:

α old + β old + γ old ≤ m α recv + β recv + γ recv + γ send = m
The following lemma is given in [START_REF] Toledo | A survey of out-of-core algorithms in numerical linear algebra[END_REF]: consider any algorithm that uses the standard way of multiplying matrices (this excludes Strassen's or Winograd's algorithm, for instance). If N A elements of A, N B elements of B and N C elements of C are accessed, then no more than K computations can be done, where

K = min (N A + N B ) N C , (N A + N C ) N B , (N B + N C ) N A . (1)
To use this result here, we see that no more than α old +α recv blocks of A are accessed, hence N A = (α old + α recv )q 2 . Similarly, N B = (β old + β recv )q 2 and N C = (γ old + γ recv )q 2 (the C blocks returned are already counted). We simplify notations by writing α old + α recv = αm, β old + β recv = βm, and γ old + γ recv = γm. Then we obtain

K = min (α+β) √ γ, (β +γ) √ α, (γ +α) β × m √ mq 3
Writing K = km √ mq 3 , we obtain the following system of equations

             MAXIMIZE k S.T. k ≤ (α + β) √ γ k ≤ (β + γ) √ α k ≤ (γ + α) β α + β + γ ≤ 2
whose solution is easily found to be α = β = γ = 2 3 , AND k = 32 27 . This gives a lower bound for the communication-to-computation ratio (in terms of blocks) of any algorithm:

CCR opt = m km √ m = 27 32m .
In fact, it is possible to refine this bound. Instead of using the lemma given in [START_REF] Toledo | A survey of out-of-core algorithms in numerical linear algebra[END_REF], we use Loomis-Whitney inequality [START_REF] Ironya | Communication lower bounds for distributed-memory matrix multiplication[END_REF]: if N A elements of A, N B elements of B, and N C elements of C are accessed, then no more than K computations can be done, where

K = √ N A N B N C . Here K = √ αβγ × m √ mq 3 . We obtain α = β = γ = 2 3
, and k = 8 27 , so that the lower bound for the communication-tocomputation ratio becomes: CCR opt = 27 8m . The maximum re-use algorithm does not achieve the lower bound:

CCR ∞ = 2 √ m = 32 8m
but it is quite close! Finally, we point out that the bound CCR opt improves upon the best-known value 1 8m derived in [START_REF] Ironya | Communication lower bounds for distributed-memory matrix multiplication[END_REF]. Also, the ratio CCR ∞ achieved by the maximum re-use algorithm is lower by a factor √ 3 than the ratio achieved by the blocked matrix-multiply algorithm of [START_REF] Toledo | A survey of out-of-core algorithms in numerical linear algebra[END_REF].

Algorithms for homogeneous platforms

In this section, we adapt the maximum re-use algorithm to fully homogeneous platforms. We must first decide which part of the memory will be used to stock which part of the original matrices, in order to maximize the total number of computations per time unit. Cannon's algorithm [START_REF] Cannon | A cellular computer to implement the Kalman filter algorithm[END_REF] and the ScaLAPACK outer product algorithm [START_REF] Blackford | ScaLAPACK Users' Guide[END_REF] both distribute square blocks of C to the processors. Intuitively, squares are better than elongated rectangles because their perimeter (which is proportional to the communication volume) is smaller for the same area. We use the same approach here, but we have not been able to assess any optimal result.

Principle of the algorithm

We load into the memory of each worker µ q × q blocks of A and µ q × q blocks of B to compute µ 2 q × q blocks of C. In addition, we need 2µ extra buffers, split into µ buffers for A and µ for B, in order to overlap computation and communication steps. In fact, µ buffers for A and µ for B would suffice for each update, but we need to prepare for the next update while computing. Overall, the number of C blocks that we can simultaneously load into memory is the largest integer µ such that µ 2 + 4µ ≤ m.

We have to determine the number of participating workers P. For that purpose, we proceed as follows. On the communication side, we know that in a round (computing a C block entirely), the master exchanges with each worker 2µ 2 blocks of C (µ 2 sent and µ 2 received), and sends µt blocks of A and µt blocks of B. Also during this round, on the computation side, each worker computes µ 2 t block updates.

If we enroll too many processors, the communication capacity of the master will be exceeded. There is a limit on the number of blocks sent per time unit, hence on the maximal processor number P, which we compute as follows: P is the smallest integer such that

2µtc × P ≥ µ 2 tw.
Indeed, this is the smallest value to saturate the communication capacity of the master required to sustain the corresponding computations. We derive that

P = µ 2 tw 2µtc = µw 2c .
In the context of matrix multiplication, we have c = q 2 τ c and w = q 3 τ a , hence P = µq 2 τa τc . Moreover, we need to enforce that P ≤ p, hence we finally obtain P = min p, µq 2 τa τc . For the sake of simplicity, we suppose that r is divisible by µ, and that s is divisible by Pµ. We allocate µ block columns (i.e., qµ consecutive columns of the original matrix) of C to each processor. The algorithm is decomposed into two parts. Algorithm 1 outlines the program of the master, while Algorithm 2 is the program of each worker.

Impact of the start-up overhead

If we follow the execution of the homogeneous algorithm, we may wonder whether we can really neglect the input/output of C blocks. Here we sequentialize the sending, computing, and receiving of the C blocks, so that each worker loses 2c time-units per block, i.e., per tw time-units. As there are P ≤ µw 2c + 1 workers, the total loss would be of 2cP time-units every tw time-units, which is less than µ t + 2c tw . For example, with c = 2, w = 4.5, µ = 4 and t = 100, we enroll P = 5 workers, and the total lost is at most 4%, which is small enough to be neglected. Note that it would be technically possible to design an algorithm where the sending of the next block is overlapped with the last computations of the current block, but the whole procedure gets much more complicated.

MPI experiments

In this section, we aim at validating the previous theoretical results and algorithms. We conduct a variety of MPI Split the matrix into squares C i ,j of µ 2 blocks (of size q × q): C i ,j = {C i,j \ (i -1)µ + 1 ≤ i ≤ i µ, (j -1)µ + 1 ≤ j ≤ j µ}; for j ← 0 to s Pµ by Step P do for i ← 1 to r µ do for id worker ← 1 to P do j ← j + id worker ; Send block C i ,j to worker id worker ; for k ← 1 to t do for id worker ← 1 to P do j ← j + id worker ; for j ← (j -1)µ + 1 to j µ do Send B k,j ; for i ← (i -1)µ + 1 to i µ do Send A i,k ;

for id worker ← 1 to P do j ← j + id worker ; Receive C i ,j from worker id worker ; experiments to compare our new schemes with several other algorithms from the literature.

Platform

For our experiments we are using a platform at the University of Tennessee. All experiments are performed on a cluster of 64 Xeon 3.2GHz dual-processor nodes running the Linux operating system. Each node has four Gigabytes of memory, but we only use 512 MB of memory to further stress the impact of limited memories. The nodes are connected with a switched 100Mbps Fast Ethernet network.

Algorithm 2: Homogeneous version, worker program.

for all blocks do Receive C i ,j from master; for k ← 1 to t do for j ← (j -1)µ + 1 to j µ do Receive B k,j ;

for i ← (i -1)µ + 1 to i µ do Receive A i,k ; for j ← (j -1)µ + 1 to j µ do C i,j ← C i,j + A i,k .B k,j ;
Return C i ,j to master;

In order to build a master-worker platform, we arbitrarily choose one processor as the master, and the other processors become the workers. Finally we used MPI WTime as timer in all experiments.

Algorithms

We choose and adapt four different algorithms from the general literature to compare our algorithm to. We partition these algorithms into two sets. The first set is composed of algorithms which use the same memory allocation than ours. The only difference between the algorithms is the order in which the master sends blocks to workers.

• Homogeneous algorithm: HoLM is our homogeneous algorithm. It makes resource selection, and sends blocks to the selected workers in a round-robin fashion. • Overlapped Round-Robin, Optimized Memory Layout:

ORROML is very similar to our homogeneous algorithm. The only difference between them is that it does not make any resource selection, and so sends tasks to all available workers in a round-robin fashion. • Overlapped Min-Min, Optimized Memory Layout:

OMMOML is a static scheduling heuristic, which sends the next block to the first worker that will be available to compute it. As it is looking for potential workers in a given order, this algorithm performs some resource selection too. Theoretically, as our homogeneous resource selection ensures that the first worker is free to compute when we finish to send blocks to the others, they should have similar behavior. • Overlapped Demand-Driven, Optimized Memory Layout: ODDOML is a demand-driven algorithm. In order to use the extra buffers available in the worker memories, it will send the next block to the first worker which can receive it. This would be a dynamic version of our algorithm, if it took worker selection into account. • Demand-Driven, Optimized Memory Layout: DDOML is a very simple dynamic demand-driven algorithm, close to ODDOML. It sends the next block to the first worker which is free for computation. As workers never have to receive and compute at the same time, the algorithm has no extra buffer, so the memory available to store A, B, and C is greater. This may change the value of µ and so the behavior of the algorithm.

In the second set we have algorithms which do not use our memory allocation:

• Block Matrix Multiply: BMM is Toledo's algorithm [START_REF] Toledo | A survey of out-of-core algorithms in numerical linear algebra[END_REF]. It splits each worker memory equally into three parts, and allocate one slot for a square block of A, another for a square block of B, and the last one for a square block of C, each square block having the same size. Then it sends blocks to the workers in a demanddriven fashion, when a worker is free for computation. First a worker receives a block of C, then it receives corresponding blocks of A and B in order to update C, until C is fully computed. In this version, a worker do not overlap computation with the receiving of the next blocks. • Overlapped Block Matrix Multiply: OBMM is our attempt to improve the previous algorithm. We try to overlap the communications and the computations of the workers. To that purpose, we split each worker memory into five parts, so as to receive one block of A and one block of B while previous ones are used to update C.

Experiments

We have built several experimental protocols in order to assess the performance of the various algorithms. In the following experiments we use nine processors, one master and eight workers. In all experiments we compare the execution time needed by the algorithms which use our memory allocation to the execution time of the other algorithms. We also point out the number of processors used by each algorithm, an important parameter when comparing execution times.

In the first set of experiments, we test the different algorithms on matrices of different sizes and shapes. The matrices we are multiplying are of actual size -8000 × 8000 for A and 8000 × 64000 for B, -16000 × 16000 for A and 16000 × 128000 for B, and -8000 × 64000 for A and 64000 × 64000 for B. All the algorithms using our optimized memory layout consider these matrices as composed of square blocks of size q × q = 80 × 80. For instance in the first case we have r = t = 100 and s = 800.

In the second set of experiments we check whether the choice of q was wise. For that purpose, we launch the algo-rithms on matrices of size 8000 × 8000 and 8000 × 64000, changing from one experiment to another the size of the elementary square blocks. Then q will be respectively equal to 40 and 80. As the global matrix size is the same in both experiments, we expect both results to be the same.

In the third set of experiments we investigate the impact of the worker memory size onto the performance of the algorithms. In order to have reasonable execution times, we use matrices of size 16000 × 16000 and 16000 × 64000, and the memory size will vary from 132MB to 512MB. We choose these values to reduce side effects due to the partition of the matrices into blocks of size µq × µq.

In the fourth and last set of experiments we check the stability of the previous results. To that purpose we launch the same execution five times, in order to determine the maximum gap between two runs.

Results and discussion

We see in Figure 4 the results of the first set of experiments, where algorithms are computing different matrices. The first remark is that the shape of the two experiments is the same for all matrix sizes. We also underline the superiority of most of the algorithms which use our memory allocation against BMM: HoLM, ORROML, ODDOML, and DDOML are the best algorithms and have similar performance. Only OMMOML needs more time to complete its execution. This delay comes from its resource selection: it uses only two workers. For instance, HoLM uses four workers, and is as competitive as the other algorithms which all use the eight available workers.

In Figure 6, we see the impact of q on the performance of our algorithms. BMM and OBMM have same execution times in the three experiments as these algorithms do not split matrices into elementary square blocks of size q × q but, instead, call the Level 3 BLAS routines directly on the whole m 3 × m 3 matrices. In the two cases we see that the time of the algorithms are similar. We point out that this experiment shows that the choice of q has little impact on the algorithms performance.

In Figure 5 we have the impact of the worker memory size on the performance of the algorithms. As expected, the performance increases with the amount of memory available. It is interesting to underline that our resource selection always performs in the best possible way. HoLM will use respectively two and four workers when the memory available increases, compared to the other algorithms which will use all eight available workers on each test. OMMOML also makes some resource selection, but it performs worse.

Finally, figure 7 shows the difference that we can have between two runs. This difference is around 6%. Thus if two algorithms have less than 6% of difference in execution time, they should be considered as similar. To conclude, these experiments stress the superiority of our memory allocation. Furthermore, our homogeneous algorithm is as competitive as the others but uses fewer resources.

Conclusion

The main contributions of this paper are the following:

1. On the theoretical side, we have derived a new, tighter, bound on the minimal volume of communications needed to multiply two matrices. From this lower bound, we have defined an efficient memory layout, i.e., an algorithm to share the memory available on the workers among the three matrices.

2. On the practical side, starting from our memory layout, we have designed an algorithm for homogeneous platforms whose performance is quite close to the communication volume lower bound.

3. Through MPI experiments, we have shown that our algorithm for homogeneous platforms has far better performance than solutions using the memory layout proposed in [START_REF] Toledo | A survey of out-of-core algorithms in numerical linear algebra[END_REF]. Furthermore, this static homogeneous algorithm has similar performance as dynamic algorithms using the same memory layout, but uses fewer processors. It is therefore a very good candidate for deploying applications on regular, homogeneous platforms.

Future work is devoted to extending the memory management strategy and the corresponding algorithms to heterogeneous platforms. Preliminary results are available in [START_REF] Dongarra | Revisiting matrix product on master-worker platforms[END_REF].

Figure 1 .

 1 Figure 1. Partition of the three matrices A, B, and C.

Figure 2 .

 2 Figure 2. Memory usage for the maximum re-use algorithm when m = 21: µ = 4; 1 block is used for A, µ for B, and µ 2 for C.

Figure 3 .

 3 Figure 3. Four steps of the maximum re-use algorithm, with m = 21 and µ = 4. The elements of C updated are displayed on white on black.

Algorithm 1 :

 1 Homogeneous version, master program. µ ← √ 4 + m -2 ; P ← min p, µw 2c ;

Figure 6 .

 6 Figure 4. Performance of the algorithms on different matrices.

Figure 7 .

 7 Figure 7. Variation of algorithm execution times.