Jack Dongarra

Jean-Franc ¸ois Pineau

Yves Robert

Zhiao Shi
email: shi@cs.utk.edu

F R É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 master-worker platforms. While matrix-product is well-understood for homogeneous 2D-arrays of processors, 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. 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. 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: we want to achieve the smallest possible execution time while using as few processors as possible. 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. In Section 3, we start with a theoretical study of the simplest version of the problem, without memory limitation, which is intended to show the intrinsic difficulty of the scheduling problem. Next, in Section 4, 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 5, we propose a scheduling algorithm that includes resource selection. We report several MPI experiments in Section 6. Section 7 is dedicated to 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. 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 linear 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 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 nonblocking communications in MPI[END_REF], who report that asynchronous MPI sends get serialized as soon as message sizes exceed a hundred kilobytes. Their result holds 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 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.

Combinatorial complexity of a simple version of the problem

This section is almost a digression; it is devoted to the study of the simplest variant of the problem. It is intended to show the intrinsic combinatorial difficulty of the problem. We make the following simplifications:

• We consider only rank-one block updates; in other words, and with previous notations, we focus on the case where t = 1.

• Results need not be returned to the master.

• Workers have no memory limitation; they receive each stripe only once and can re-use them for other computations.

There are five parameters in the problem; three platform parameters (c, w, and the number of workers p) and two application parameters (r and s). The scheduling problem amounts to deciding which files should be sent to which workers and in which order. A given file may well be sent several times, to further distribute computations. For instance, a simple strategy is to partition A and to duplicate B, i.e., send each block A i only once and each block B j p times; all workers would then be able to work fully in parallel. The dependence graph of the problem is depicted in Figure 2. It suggests a natural strategy for enabling workers to start computing as soon as possible. Indeed, the master should alternate sending A-blocks and B-blocks. Of course it must be decided how many workers to enroll and in which order to send the blocks to the enrolled workers. But with a single worker, we can show that the alternating greedy algorithm is optimal: Proposition 1 With a single worker, the alternating greedy algorithm is optimal.

C 31 A 2 A 1 A 3 C 11 C 21 C 12 C 22 C 32 B 2 B 1
Proof. In this algorithm, the master sends blocks as soon as possible, alternating a block of type A and a block of type B (and proceeds with the remaining blocks when one type is exhausted). This strategy maximizes at each step the total number of tasks that can be processed by the worker. To see this, after x communication steps, with y files of type A sent, and z files of type B sent, where y + z = x, the worker can process at most y × z tasks. The greedy algorithm enforces y = x 2 and z = x 2 (as long as max(x, y) ≤ min(r, s), and then sends the remaining files), hence its optimality.

2 Unfortunately, for more than one worker, we did not succeed in determining an optimal algorithm. The difficulty lies in the selection of the next worker to enroll. One can see this on the following example. There are (at least) two greedy algorithms that can be devised for p workers:

Min-Min Thrifty A2 B1 A1 B2 B3 C11 C12 C21 C22 C13 C23 B1 A3 B2 B3 C31 C32 C33 B3 A2 B2 A1 B1 C11 C12 C21 C22 C13 C23 C33 C32 C31 A3 B1 B2 B3 A2 P1 P2 P1 P2 Figure 3:
Example showing that Thrifty is not optimal: with p = 2, c = 4, w = 7, and r = s = 3, Min-min has a lower makespan. For each algorithm, the first line presents the communications for P 1 , the second the computations of P 1 , the third the communications for P 2 , and the fourth the computations of P 2 . Thrifty: This algorithm "spares" resources as it aims at keeping each enrolled worker fully active. It works as follows:

A1 B2 C11 C12 A2 B3 C21 C22 B1 A6 C13 A3 C61 B2 C23 C31 C32 C33 A4 C41 C42 C43 A5 B3 C62 C51 C52 C63 C53 B1 A1 B2 A2 C11 C12 C21 C22 B3 B1 A6 C61 B2 A5 C23 C13 C62 A3 A4 C31 B3 C51 C52 C63 C53 C43 A4 C32 C33 C41 C42 B1 P1 P2 P2 P1 Thrifty Min-Min
• Send enough blocks to the first worker so that it is never idle,

• Send blocks to a second worker during spare communication slots, and

• Enroll a new worker (and send blocks to it) only if this does not delay previously enrolled workers.

Min-min: This algorithm is based on the well-known min-min heuristic [START_REF] Maheswaran | Dynamic matching and scheduling of a class of independent tasks onto heterogeneous computing systems[END_REF]. At each step, all tasks are considered. For each of them, we compute their possible starting date on each worker, given the files that have already been sent to this worker and all decisions taken previously; we select the best worker, hence the first min in the heuristic. We take the minimum of starting dates over all tasks, hence the second min.

It turns out that neither greedy algorithm is optimal. See Figure 3 for an example where Min-min is better than Thrifty, and Figure 4 for an example of the opposite situation.

We now go back to our original model.

C C C C C C C C C C C C C C C C A B B B B 1 µ µ 2
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

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 satisfying our hypotheses (centralized data and limited memory). 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 5. Four consecutive execution steps are shown in Figure 6. 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 5, 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 need 2µ 2 communications to send and retrieve C blocks. 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

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.

How can we assess the performance of the maximum re-use 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 and Winograd's algorithms, 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 .
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 [START_REF] Li | Parallel matrix multiplication on a linear array with a reconfigurable pipelined bus system[END_REF] . 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 [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 matrixmultiply algorithm of [START_REF] Toledo | A survey of out-of-core algorithms in numerical linear algebra[END_REF].

K = √ N A N B N C . Here K = √ αβγ × m √ mq 3 . We obtain α = β = γ = 2

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 , where τ c and τ a respectively represent the speed of the communication link and the speed of the processor. 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 send- 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 ; ing, 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 technically be 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 get much more complicated.

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 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. In order to Algorithm 2: 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; build a master-worker platform, we arbitrarily choose one processor as the master, and the other processors will serve as workers. Finally we used MPI WTime as timer in all experiments.

Algorithms

We choose four different algorithms from the general literature, and adapt them into six algorithms designed for matrix multiplication. The objective is to use these six algorithms as the basis for comparison with our new algorithm. The four original algorithms are Round-Robin, Minimum Completion Time (or Min-Min) and Demand-Driven, which are well known algorithms in the scheduling literature, and Toledo's algorithm [START_REF] Toledo | A survey of out-of-core algorithms in numerical linear algebra[END_REF], which is a dedicated algorithm for matrix multiplication.

As the first three chosen algorithms are designed for job distribution, and not for matrix multiplication, we adapt them in the following way: they all use our memory layout to divide matrices into chunks and to determine in which order chunks have to be sent to participating workers. The only difference between these algorithms and ours is the order in which the master serves workers.

All the algorithms used during the MPI experiments are now described. First, our algorithm:

• Homogeneous algorithm: HoLM is our homogeneous algorithm, described in section 5. It makes resource selection, and sends blocks to the selected workers in a round-robin fashion.

The other four algorithms using our memory layout:

• 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 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 to send blocks to the others, HoLM and OMMOML should have similar behavior on homogeneous platforms. • 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 compare to ODDOML.

Finally, the two adapted algorithms which do not use our memory allocation:

• Block Matrix Multiply: BMM is Toledo's algorithm. 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, each square block having the same size. Then it sends blocks to the workers in a demand-driven 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 does 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. We restricted the number of workers to eight after an initial experiment using the whole platform. This experiment showed that, because of platform parameters and of memory limitations, if the master serves five processors or more, on average one processor is idle at any time. 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 algorithms 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 7 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 only uses three workers. For instance, HoLM uses four workers, and is as competitive as the other algorithms which all use the eight available workers. This difference can be explained by the resource selection process. Before computing the schedules, we measure the communication and computation times by sending the same task to all slaves. The different measures obtained are not absolutely identical, for instance because of experimental measurement errors. OMMOML uses all the exact measures. As HoLM has to consider an absolutely homogeneous platform, it works on a more pessimistic platform using the maximum of the computation times, and the maximum of the communication times.

In Figure 8 we have the impact of the worker memory size on the performance of the algorithms. As expected, the performance increases (as the execution time needed to perform the multiplication decreases) 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.

In Figure 9, we see the impact of q on the performance of our algorithms. BMM and OBMM have constant execution times in the 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 execution 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.

Finally, figure 10 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.

Related work

In this section, we provide a brief overview of the related work. Even if, in this article, we only focused on homogeneous platforms, we consider the related work in the more general scope of heterogeneous platforms. We classify the related papers along the following six main lines:

Load balancing on heterogeneous platforms -Load balancing strategies for heterogeneous platforms have been widely studied. Distributing the computations (together with the associated data) can be performed either dynamically or statically, or a mixture of both. Some simple schedulers are available, but they use naive mapping strategies such as master-worker techniques or paradigms based upon the idea "use the past to predict the future", i.e. use the currently observed speed of computation of each machine to decide for the next distribution of work [START_REF] Cierniak | Compile-time scheduling algorithms for heterogeneous network of workstations[END_REF][START_REF] Cierniak | Customized dynamic load balancing for a network of workstations[END_REF][START_REF] Berman | High-performance schedulers[END_REF]. Dynamic strategies such as self-guided scheduling [START_REF] Polychronopoulos | Compiler optimization for enhancing parallelism and their impact on architecture design[END_REF] could be useful too. There is a challenge in determining a trade-off between the data distribution parameters and the process spawning and possible migration policies. Redundant computations might also be necessary to use a heterogeneous cluster at its best capabilities. However, dynamic strategies are outside the scope of this paper (but mentioned here for the sake of completeness). Because we have a library designer's perspective, we concentrate on static allocation schemes that are less general and more difficult to design than dynamic approaches, but which are better suited for the implementation of fixed algorithms such as linear algebra kernels from the ScaLAPACK library [START_REF] Blackford | ScaLAPACK Users' Guide[END_REF].

Out-of-core linear algebra routines -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].

Matrix product on reconfigurable architectures-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 [37]. 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.

Linear algebra algorithms on heterogeneous clusters -Several authors have dealt with the static implementation of matrix-multiplication algorithms on heterogeneous platforms. One simple approach is given by Kalinov and Lastovetsky [START_REF] Kalinov | Heterogeneous distribution of computations while solving linear algebra problems on networks of heterogeneous computers[END_REF]. Their idea is to achieve a perfect load-balance as follows: first they take a fixed layout of processors arranged as a collection of processor columns; then the load is evenly balanced within each processor column independently; next the load is balanced between columns; this is the "heterogeneous block cyclic distribution" of [START_REF] Kalinov | Heterogeneous distribution of computations while solving linear algebra problems on networks of heterogeneous computers[END_REF]. Another approach is proposed by Crandall and Quinn [START_REF] Crandall | Block data decomposition for data-parallel programming on a heterogeneous workstation network[END_REF], who propose a recursive partitioning algorithm, and by Kaddoura, Ranka, and Wang [START_REF] Kaddoura | Array decomposition for nonuniform computational environments[END_REF], who refine the latter algorithm and provide several variations. They report several numerical simulations. As pointed out in the introduction, theoretical results for matrix multiplication and LU decomposition on 2D-grids of heterogeneous processors are reported in [START_REF] Beaumont | A proposal for a heterogeneous cluster ScaLAPACK (dense linear solvers)[END_REF], while extensions to general 2D partitioning are considered in [START_REF] Beaumont | Matrix multiplication on heterogeneous platforms[END_REF]. See also Lastovetsky and Reddy [START_REF] Lastovetsky | Data partitioning with a realistic performance model of networks of heterogeneous computers[END_REF] for another partitioning approach.

Recent papers aim at making easier the process of tuning linear algebra kernels on heterogeneous systems. Self-optimization methodologies are described by Cuenca et al. [START_REF] Cuenca | Processes distribution of homogeneous parallel linear algebra routines on heterogeneous clusters[END_REF] and by Chen et al. [START_REF] Chen | Self adapting software for numerical linear algebra and lapack for clusters[END_REF]. Along the same line, Chakravarti et al. [START_REF] Chakravarti | Self-organizing scheduling on the organic grid[END_REF] describe an implementation of Cannon's algorithm using self-organizing agents on a peer-to-peer network.

Models for heterogeneous platforms -In the literature, one-port models come in two variants. In the unidirectional variant, a processor cannot be involved in more than one communication at a given time-step, either a send or a receive. This is the model that we have used throughout the paper. In the bidirectional model, a processor can send and receive in parallel, but at most to a given neighbor in each direction. In both variants, if P u sends a message to P v , both P u and P v are blocked throughout the communication.

The bidirectional one-port model is used by Bhat et al. [START_REF] Bhat | Efficient collective communication in distributed heterogeneous systems[END_REF][START_REF] Bhat | Efficient collective communication in distributed heterogeneous systems[END_REF] for fixed-size messages. They 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 nonblocking communications in MPI[END_REF], who report that asynchronous MPI sends get serialized as soon as message sizes exceed a few megabytes. Their results hold for two popular MPI implementations, MPICH on Linux clusters and IBM MPI on the SP2.

The one-port model fully accounts for the heterogeneity of the platform, as each link has a different bandwidth. It generalizes a simpler model studied by Banikazemi et al. [START_REF] Banikazemi | Efficient collective communication on heterogeneous networks of workstations[END_REF], Liu [START_REF] Liu | Broadcast scheduling optimization for heterogeneous cluster systems[END_REF], and Khuller and Kim [START_REF] Khuller | On broadcasting in heterogenous networks[END_REF]. In this simpler model, the communication time only depends on the sender, not on the receiver. In other words, the communication speed from a processor to all its neighbors is the same.

Finally, we note that some papers [START_REF] Banikazemi | Communication modeling of heterogeneous networks of workstations for performance characterization of collective operations[END_REF][START_REF] Bar-Noy | Message multicasting in heterogeneous networks[END_REF] depart form the one-port model as they allow a sending processor to initiate another communication while a previous one is still on-going on the network. However, such models insist that there is an overhead time to pay before being engaged in another operation, so they are not allowing for fully simultaneous communications.

Master-worker on the computational grid -Master-worker scheduling on the grid can be based on a network-flow approach [START_REF] Shao | Master/slave computing on the grid[END_REF][START_REF] Shao | Adaptive scheduling of master/worker applications on distributed computational resources[END_REF] or on an adaptive strategy [START_REF] Heymann | Adaptive scheduling for masterworker applications on the computational grid[END_REF]. Note that the network-flow approach of [START_REF] Shao | Master/slave computing on the grid[END_REF][START_REF] Shao | Adaptive scheduling of master/worker applications on distributed computational resources[END_REF] is possible only when using a full multiple-port model, where the number of simultaneous communications for a given node is not bounded. This approach has also been studied in [START_REF] Hong | Bandwidth-aware resource allocation for heterogeneous computing systems to maximize throughput[END_REF]. Enabling frameworks to facilitate the implementation of master-worker tasking are described in [START_REF] Goux | An enabling framework for master-worker applications on the computational grid[END_REF][START_REF] Weissman | Scheduling multi-component applications in heterogeneous widearea networks[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., a way 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 theoretical 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.

Figure 1 :

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

Figure 2 :

 2 Figure 2: Dependence graph of the problem (with r = 3 and s = 2).

Figure 4 :

 4 Figure 4: Example showing that Min-min is not optimal: with p = 2, c = 8, w = 9, r = 6, and s = 3, Thrifty has a lower makespan.

Figure 5 :

 5 Figure 5: 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 6 :

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

3 , and k = 8 27= 32 8m 1 8m

 38321 , so that the lower bound for the communicationto-computation ratio becomes: CCR opt = 27 8m . The maximum re-use algorithm does not achieve the lower bound: CCR ∞ = 2 √ m but it is quite close! Finally, we point out that the bound CCR opt improves upon the best-known value derived in

Algorithm 1 :

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

Figure 7 :

 7 Figure 7: Performance of the algorithms on different matrices.

Figure 8 :

 8 Figure 8: Impact of memory size on algorithm performance.

Figure 9 :Figure 10 :

 910 Figure 9: Impact of block size q on algorithm performance.

Acknowledgments. We thank the reviewers for their numerous comments and suggestions, which greatly improved the final version of the paper.