Guillaume Aupy

Anne Benoit

Thomas Hérault

Yves Robert

Jack Dongarra

Optimal Checkpointing Period: Time vs. Energy

Keywords: Fault-tolerance, checkpoint, energy, algorithm, model, exascale

This short paper deals with parallel scientific applications using non-blocking and periodic coordinated checkpointing to enforce resilience. We provide a model and detailed formulas for total execution time and consumed energy. We characterize the optimal period for both objectives, and we assess the range of time/energy trade-offs to be made by instantiating the model with a set of realistic scenarios for Exascale systems. We give a particular emphasis to I/O transfers, because the relative cost of communication is expected to dramatically increase, both in terms of latency and consumed energy, for future Exascale platforms.

Périodes optimales de checkpoint : temps vs énergie

Résumé : Dans ce court papier, nous considérons des applications scientifiques parallèles, utilisant des protocoles de sauvegarde de points de reprise (checkpoints) coordonnés et périodiques afin de se préserver des fautes. Nous fournissons un modèle et des formules détaillées pour le temps total d'exécution et la consommation d'énergie de ces applications. Nous caractérisons pour les deux objectifs, les périodes optimales, et par de nombreuses simulations réalistes pour des systèmes exascales, montrons le compromis temps/énergie induit par ces périodes. Finalement, nous insistons particulièrement dans cette étude sur le coût des tranferts I/O : il est attendu que le coût relatif des communications augmente drastiquement, à la fois en termes de latence, et en termes énergétiques, pour les futures plateformes exascales.

Mots-clés : Tolérance aux pannes, checkpoint, énergie, algorithmes, modèle, exascale

Introduction

A significant research effort is focusing on the characteristics, features, and challenges of High Performance Computing (HPC) systems capable of reaching the Exaflop performance mark [START_REF] Dongarra | The international exascale software project: a call to cooperative action by the global high-performance community[END_REF][START_REF] Sarkar | Exascale software study: Software challenges in extreme scale systems[END_REF]. The portrayed Exascale systems will necessitate billion way parallelism, resulting not only in a massive increase in the number of processing units (cores), but also in terms of computing nodes. Considering the relative slopes describing the evolution of the reliability of individual components on one side, and the evolution of the number of components on the other side, the reliability of the entire platform is expected to decrease, due to probabilistic amplification. Even if each independent component is quite reliable, the Mean Time Between Failures (MTBF) is expected to drop drastically. Executions of large parallel applications on these systems will have to tolerate a higher degree of errors and failures than in current systems. The defacto general-purpose error recovery technique in high performance computing is checkpoint and rollback recovery. Such protocols employ checkpoints to periodically save the state of a parallel application, so that when an error strikes some process, the application can be restored into one of its former states. The most widely used protocol is coordinated checkpointing, where all processes periodically stop computing and synchronize to write critical application data onto stable storage. Coordinated checkpointing is well understood, at least in its blocking form (when no computing activity takes place during checkpoints), and good approximations of the optimal checkpoint interval exist; they are known as Young's and Daly's formula [START_REF] Young | A first order approximation to the optimum checkpoint interval[END_REF][START_REF] Daly | A higher order estimate of the optimum checkpoint interval for restart dumps[END_REF].

While reliability is a major concern for Exascale, another key challenge is to minimize energy consumption, both for economic and environmental reasons. One of the most power-consuming components of today's systems is the processor: even when idle, it dissipates a significant fraction of the total power. However, for future Exascale systems, the power dissipated to execute I/O transfers is likely to play an even more important role, because the relative cost of communication is expected to dramatically increase, both in terms of latency and consumed energy [START_REF] Shalf | Exascale computing technology challenges[END_REF].

In this short paper, we investigate trade-offs between execution time and energy consumption for the execution of parallel applications on future Exascale systems. The optimal period T opt Time given by Young's and Daly's formula [START_REF] Young | A first order approximation to the optimum checkpoint interval[END_REF][START_REF] Daly | A higher order estimate of the optimum checkpoint interval for restart dumps[END_REF] will minimize (expected) execution time. But will it minimize energy consumption? The answer is negative, mainly because the fraction of power P Cal spent when computing (by the CPUs) is not the same as the fraction of power P I/O spent when checkpointing. In particular, we revisit the work of Meneses, Sarood and Kalé [START_REF] Meneses | Assessing Energy Efficiency of Fault Tolerance Protocols for HPC Systems[END_REF] for checkpoint/restart, where formulas are given to compute the time-optimum and energy-optimum periods. However, our model is more precise: (i) we carefully assess the impact of the power consumption required for I/O activity, which is likely to play a key role at the Exascale; (ii) we consider non-blocking checkpointing that can be partially overlapped with computations; (iii) we give a more accurate analysis of the consumed energy.

Altogether, this short paper provides the following main contributions:

• We provide a refined analytical model to compute both the execution time and the consumed energy with a given checkpoint period. The model handles the case where checkpointing activity can be non-blocking, i.e., partially overlapped with computations.

• We provide analytical formulas to approximate the optimal period for time T opt Time as well as the optimal period for energy T opt Energy , thereby refining and extending Daly [START_REF] Daly | A higher order estimate of the optimum checkpoint interval for restart dumps[END_REF] and Meneses, Sarood and Kalé [START_REF] Meneses | Assessing Energy Efficiency of Fault Tolerance Protocols for HPC Systems[END_REF] results to nonblocking checkpoints.

• We assess the range of time/energy trade-offs to be made by instantiating the model with a set of realistic scenarios for Exascale systems.

Model

In this section, we introduce all the model parameters. We start with parameters related to resilience (checkpointing) before moving to parameters related to energy consumption.

Checkpointing

We model coordinated checkpointing [START_REF] Chandy | Distributed snapshots: Determining global states of distributed systems[END_REF] where checkpoints are taken at regular intervals, after some fixed amount of work units have been performed. This corresponds to an execution partitioned into periods of duration T . Every period, a checkpoint of length C is taken. An important question is whether checkpoints are blocking or not. On some architectures, we may have to stop executing the application before writing to the stable storage where the checkpoint data is saved; in that case checkpoint is fully blocking. On other architectures, checkpoint data can be saved on the fly into a local memory before the checkpoint is sent to the stable storage, while computation can resume progress; in that case, checkpoints can be fully overlapped with computations. To deal with all situations, we introduce a slowdown factor ω: during a checkpoint of duration C, the work that is performed is ωC work units. In other words, (1 -ω)C work units are wasted due to checkpoint jitter disrupting the progress of computation. Here, 0 ≤ ω ≤ 1 is an arbitrary parameter. The case ω = 0 corresponds to a fully blocking checkpoint, while ω = 1 corresponds to a checkpoint totally overlapped with computations. All intermediate situations can be represented.

Next we have to account for failures. During t time units of execution, the expectation of the number of failures is t µ , where µ is the MTBF (Mean Time Between Failures) of the platform. Note that if the platform if made of N identical resources whose individual mean time between failures is µ ind , then µ = µ ind N . This relation is agnostic of the granularity of the resources, which can be anything from a single CPU to a complex multi-core socket. When a failure strikes, there is a downtime of length D (time to reboot the resource or set up a spare), and then a recovery of length R (time to read the last stored checkpoint). The work executed by the application since the last checkpoint and before the failure needs to be re-executed. Clearly, the shorter the period T , the less work to re-execute, but also the more overhead due to frequent checkpoints in a failure-free execution. The best trade-off when ω = 0 (blocking checkpoint) is achieved for T = √ 2Cµ+C (Young's formula [START_REF] Young | A first order approximation to the optimum checkpoint interval[END_REF]) or T = 2C(µ + D + R)+C (Daly's formula [START_REF] Daly | A higher order estimate of the optimum checkpoint interval for restart dumps[END_REF]). Both formulas are first-order approximations and valid only if all checkpoint parameters C, D and R are small in front of µ (and these

RR n°8387

formulas collapse if they become negligible). In Section 3, we show how to extend these formulas to the case of non-blocking checkpoints (see also [START_REF] Bosilca | Unified model for assessing checkpointing protocols at extreme-scale[END_REF] for more details).

Energy

To compute the energy consumption of the application, we need to consider the energy consumption of the different phases, and hence the power consumption at each time-step. To this purpose, we define:

• P Static : this is the base power consumed when the platform is switched on.

• P Cal : when the platform is active, we have to consider the CPU overhead in addition to the static power P Static .

• P I/O : similarly, this is the power overhead due to file I/O. This supplementary power consumption is induced by checkpointing, or when recovering from a failure.

• P Down : for coordinated checkpointing, when one processor fails, the rest of the machine stays idle. P Down is the power consumption overhead when one machine is down, that may be incurred for instance by rebooting the machine. In general, we let P Down = 0.

Meneses, Sarood and Kalé [START_REF] Meneses | Assessing Energy Efficiency of Fault Tolerance Protocols for HPC Systems[END_REF] have a simpler model with two parameters, namely L, the base power (corresponding to P Static with our notations), and H, the maximum power (corresponding to P Static + P Cal with our notations). They use P I/O = P Down = 0.

In Section 3, we show how to compute the optimal period that minimizes the energy consumption. In Section 4, we instantiate the model with expected values for power consumption of Exascale platforms.

Optimal checkpointing period

We consider a parallel application whose execution time is T base without any overhead due to the resilience method or the occurrence of failures. We compute the expectation T final of the total execution time (accounting both for checkpointing and for failures) in Section 3.1, and the expectation E final of the total energy consumed during this execution of length T final in Section 3.2. We will compute the optimal period T that minimizes the objective, either T final or E final .

Execution time

The total execution time T final of the application depends on two sources of overhead. We first compute T ff , the time taken by a fault-free execution, thereby accounting only for the overhead due to periodic checkpointing. Then we compute T fails , the time lost due to failures. Finally, T final = T ff + T fails . We detail here both computations:

• The reasoning to derive T ff is simple. We need to execute a total amount of work equal to T base . During each period of length T , there is an amount of time T -C where only computations take place, and an amount of time C of checkpointing, where only a work ωC is done. Therefore, the total number of work units executed during a period of length

T is T -C + ωC = T -(1 -ω)C, and
T ff = T base T T -(1 -ω)C .
• The reasoning to compute T fails is the following. Since the mean time between two failures is µ, the average number of failures during execution is T final µ . For each failure, the time lost is expressed as: -D + R for downtime and recovery; -a time ωC for the work that was done during the previous checkpoint and that has to be redone because it was not checkpointed (because of the failure);

with probability T -C T , the failure happens while we are not checkpointing, and the time lost is on average A = T -C 2 ; -otherwise, with probability C T , the failure happens while we are checkpointing, and the time lost is on average

B = T -C + C 2 = T -C 2 .
The time lost for each failure is

D + R + ωC + T -C T A + C T B = D + R + ωC + T 2 .
Finally,

T fails = T final µ D + R + ωC + T 2 .
We are now ready to express the total execution time:

T final = T ff + T fails = T base T T -(1 -ω)C + T final µ D + R + ωC + T 2 = T (T -(1 -ω)C) 1 -D+R+ωC+T /2 µ T base = T (T -a) b -T 2µ T base ,
where a = (1 -ω)C and b = 1 -D+R+ωC µ . This equation is minimized for

T opt Time = 2(1 -ω)C(µ -(D + R + ωC)). (1)
When ω = 0, we obtain an expression close to that of Young and Daly, but slightly different because they have less accurately approximated the total execution time. In the following, we let AlgoT be the checkpointing strategy that checkpoints with period T opt Time .

Energy consumption

In order to compute the total energy consumption of the execution, we consider the different phases during which the different powers introduced in Section 2.2 are used:

• First, we consume P Static during each time-step of the execution. Indeed, even when a node fails and is shutdown, we still pay for the power of all the other nodes, for the cooling system, etc. The corresponding energy cost is T final P Static .

• Next, let T Cal be the time during which the CPU is used, inducing a power overhead P Cal . T Cal includes the base work T base , and T re-exec , the work that must be re-executed after each failure (which we multiply by the number of failures T final /µ):

with probability T -C T , the failure does not happen during a checkpoint, and the work to re-execute is A = ωC + T -C 2 ; -with probability C T , the failure happens during the execution of a checkpoint, and the work to re-execute is

B = ωC + T -C + ωC 2 .
We derive

T re-exec = T -C T A + C T B, hence T re-exec = ωC + T 2 -C 2 2T + ωC 2 2T .
Finally, we have:

T Cal = T base + T final µ ωC + T 2 -C 2 2T + ωC 2 2T .
The corresponding energy consumption is T Cal P Cal .

• Let T I/O be the time during which the I/O system is used, inducing a power overhead P I/O . This time corresponds to checkpointing and recovery from failures.

-The total number of checkpoints that are taken in a fault-free execution is equal to the number of periods,

T base T -(1-ω)C
, and the time taken by checkpoints is therefore

T base C T -(1-ω)C .
-For each failure, there is an additional overhead:

1. the system needs to recover, which lasts R time-steps; 2. with probability T -C T , the failure does not happen during a checkpoint, and there is no additional I/O overhead; 3. however, with probability C T , the failure happens during a checkpoint, and the I/O time wasted is (in average) C 2 . Altogether, we obtain

T I/O = T base C T -(1 -ω)C + T final µ R + C 2 2T .
The corresponding energy consumption is T I/O P I/O .

RR n°8387

• Finally, let T Down be the total down time, incurring a power overhead P Down . We have

T Down = T final µ D,
and the corresponding energy cost is T Down P Down . This term is only included for full generality, as we expect to have P Down = 0 in most scenarios.

The final expression for the total energy consumed is

E final = T Cal P Cal + T I/O P I/O + T Down P Down + T final P Static = T base + T final µ ωC + T 2 -C 2 2T + ωC 2 2T P Cal + T final µ R + C 2 2T + C T base T -(1 -ω)C P I/O + T final µ DP Down + T final P Static .
It is important to understand that T final = T Cal + T I/O + T Down , unless ω = 0. Indeed, CPU and I/O activities are overlapped (and both consumed) when checkpointing. To ease the derivation of the optimal period that minimizes E final , we introduce some notations and let P Cal = αP Static , P I/O = βP Static , and P Down = γP Static . Re-using parameters a = (1 -ω)C and b = 1 -D+R+ωC µ from Section 3.1, we obtain:

T ′ final T base = -ab + T 2 2µ (T -a) 2 b -T 2µ 2
, and

E ′ final P Static = T ′ final µ αωC +βR+γD+ αT 2 - α(1 -ω)C 2 2T + βC 2 2T +µ + T final 2µ α+ α(1 -ω)C 2 T 2 - βC 2 T 2 - βCT base (T -(1-ω)C) 2 .
Then, letting

K = (T -a) 2 (b-T 2µ ) 2
PStaticT base , we have:

RR n°8387 KE ′ final = -ab+ T 2 2µ µ (αωC +βR+γD+µ)+ αT 2 + α(1 -ω)C 2 2T + βC 2 2T + (T -a)(b-T 2µ ) 2µ α+ α(1 -ω)C 2 -βC 2 T -βC b- T 2µ 2 = T 3 1 4µ - 1 4µ +T 2 αωC +βR+γD 2µ 2 + b+ a 2µ 2µ - βC 4µ 2 + 1 2µ + T - ab 2µ - ab 2µ + βCb µ -2 (α(1 -ω)-β)C 2 4µ 2 -βCb 2 - ab (αωC +βR+γD+µ) µ - b 2µ - a 4µ 2 (α(1 -ω)-β)C 2 + 1 T (α(1 -ω)-β) C 2µ -(α(1 -ω)-β) C 2µ = T 2 αωC +βR+γD 2µ 2 + b 2µ + a-βC 4µ 2 + 1 2µ + T (βC -a)b µ -2 (α(1 -ω)-β)C 2 4µ 2 - ab (αωC +βR+γD+µ) µ -βCb 2 + b 2µ + a 4µ 2 (α(1 -ω)-β)C 2 .
Let T opt Energy be the only positive root of this quadratic polynomial in T :

T opt
Energy is the value that minimizes E final . In the following, we let AlgoE be the checkpointing strategy that checkpoints with period T opt Energy . As a side note, let us emphasize the differences with the approach of Meneses, Sarood and Kalé [START_REF] Meneses | Assessing Energy Efficiency of Fault Tolerance Protocols for HPC Systems[END_REF] when restricting to the case ω = 0 (because they only consider the blocking variant). For each failure, they consider that:

• energy lost due to re-execution is T -2C 2 P Cal , while we have

T -C T T -C 2 + C T (T -C) P Cal = T 2 -C 2 2T P Cal
• energy lost due to I/O is CP I/O , while we have C 2 2T P I/O . Theses differences come from our more detailed analysis of the impact of the failure location, which can strike either during the computation phase, or during the checkpointing phase, of the whole period.

Experiments

In this section, we instantiate the previous model with scenarios taken from current projections for Exascale platforms [START_REF] Dongarra | The international exascale software project: a call to cooperative action by the global high-performance community[END_REF][START_REF] Sarkar | Exascale software study: Software challenges in extreme scale systems[END_REF][START_REF] Shalf | Exascale computing technology challenges[END_REF][START_REF] Ferreira | Evaluating the Viability of Process Replication Reliability for Exascale Systems[END_REF]. We choose realistic values for all model parameters: this includes all types of power consumption (P Static , P Cal , P I/O and P Down ), all checkpoint parameters (C, R, D and ω), and the platform MTBF µ. We start with a word of caution: our choices for these parameters may be somewhat arbitrary, and do not cover the whole range of scenarios that can be investigated. However, a key feature of our model is its robustness: as long as µ is reasonably large in front of checkpoint times, the model is able to accurately predict the best period for execution time and for energy consumption.

The power consumption of an Exascale machine is capped to 20 Mega-watts. With 10 6 nodes, this represents a nominal power of 20 milli-watts per node. Let us express all power values in milli-watts. A reasonable scenario is to assume that half this power is used for operating the platform, hence to let P Static = 10. The overhead due to computing would represent the other half, hence P Cal = 10. As for communications and I/Os, which are expected to cost an order of magnitude more than computing [START_REF] Shalf | Exascale computing technology challenges[END_REF], we take an overhead of 100, hence P I/O = 100. A key parameter for the experimental study is the ratio

ρ = P Static + P I/O P Static + P Cal = 1 + β 1 + α . (2)
With our values, we get ρ = 5.5. Note that if we used P Static = 5 and kept the same overheads 10 and 100 for computing and I/O respectively, we would get P Cal = 10, P I/O = 100, and ρ = 7. These two representative values of ρ (ρ = 5.5 and ρ = 7) are emphasized by vertical arrows in the plots below on Figure 1. As for P Down , the power during downtime, we use P Down = 0, meaning that during downtime we only account for the static power P Static of the processors that are idle.

The Jaguar platform, with N = 45, 208 processors, is reported to have experienced about one fault per day [START_REF] Zheng | A scalable double in-memory checkpoint and restart scheme towards exascale[END_REF], which leads to an individual (processor) MTBF µ ind equal to 45,208 365 ≈ 125 years. Therefore, we set the individual (processor) MTBF to µ ind = 125 years. Letting the total number of processors N vary from N = 219, 150 to N = 2, 191, 500 (future exascale platforms), the platform MTBF µ varies from µ = 300 min (5 hours) down to µ = 30 min. The experiments use resilience parameters that are representative of current and forthcoming large-scale platforms [START_REF] Ferreira | Evaluating the Viability of Process Replication Reliability for Exascale Systems[END_REF][START_REF] Cappello | Preventive migration vs. preventive checkpointing for extreme scale supercomputers[END_REF]. We take C = R = 10 min, D = 1 min, and ω = 1/2.

On Figures 1 and2, we evaluate the impact of the ratio ρ (see Equation ( 2)) on the gain in energy and loss in time of AlgoE with respect to AlgoT. The general trend is that using AlgoE can lead to significant gains in energy at the price of a small increase in execution time.

We then study in Figure 3 the scalability of the approach on forthcoming platforms. We set the duration of the complete checkpoint and rollback (C and R, respectively) to 1 minute, independently of the number of processors, and we let the downtime D equal to 0.1 minutes. It is reasonable to consider that checkpoint storage time will not increase with the number of nodes in the future, but on the contrary will remain constant. Indeed, system designers are studying a couple of alternative approaches. One consists in featuring each computing node with local storage capability, ensuring through the hardware that this storage will remain available during a failure of the node. Another approach consists in using the memory of the other processors to store the checkpoint, pairing nodes as "buddies", thus allowing to take advantage of the high bandwidth capability of the high speed network to design a scalable checkpoint storage mechanism [START_REF] Zheng | FTC-Charm++: an in-memory checkpoint-based fault tolerant runtime for Charm++ and MPI[END_REF][START_REF] Ni | Hiding checkpoint overhead in HPC applications with a semi-blocking algorithm[END_REF][START_REF] Dongarra | Revisiting the double checkpointing algorithm[END_REF][START_REF] Rajachandrasekar | A 1 PB/s file system to checkpoint three million MPI tasks[END_REF].

The MTBF for 10 6 nodes is set to 2 hours, and this value scales linearly with the number of components. Given these parameters, Figures 3a and3b shows (i) the execution time ratio of AlgoE over AlgoT, and (ii) the energy consumption ratio of AlgoT over AlgoE, both as a function of the number of nodes. Figures 3a and3b confirm the important gain in energy that can be achieved, namely up to 30% for a time overhead of only 12%. When the number of nodes gets very high (up to 10 8 ), then we observe that both energy and time ratios converge to 1. Indeed, when C becomes of the order of magnitude of the MTBF, then both periods T opt Time and T opt Energy become close to C to account for the higher failure rate.

Conclusion

In this short paper, we have provided a detailed analysis to compute the optimal checkpointing period, when the checkpointing activity can be partially overlapped with computations. We have considered two distinct objectives: either the goal is to minimize the total execution time, or it is to minimize the total energy consumption. Because of the different power consumption overheads due to computations and I/Os, we obtain different optimal periods.

We have instantiated the formulas with values derived from current and future Exascale platforms, and we have studied the impact of the power overhead due to I/O activity on the gains in time and energy. With current values, we can save more than 20% of energy with an MTBF of 300 min, at the price of an increase of 10% in the execution time. The maximum gains are expected for a platform with between 10 6 and 10 7 processors (up to 30% energy savings).

Our analytical model is quite flexible and can easily be instantiated to investigate scenarios that involve a variety of resilience and power consumption parameters.

Figure 1 -

 1 Figure 1 -Time and energy ratios as a function of ρ, with C = R = 10 min, D = 1 min, γ = 0, ω = 1/2, and various values for µ.

Figure 2 -

 2 Figure 2 -Ratios of the different strategies with C = R = 10 min, D = 1 min, γ = 0, ω = 1/2 as a function of µ and ρ.

Figure 3 -

 3 Figure 3 -Ratios of total energy and time for the two period strategies, as a function of the number of nodes, with µ = 120 min for 10 6 nodes, C = R = 1 min, D = 0.1 min, γ = 0, ω = 1/2.

RR n°8387

Acknowledgements

This work was supported in part by the ANR RESCUE project. A. Benoit and Y. Robert are with the Institut Universitaire de France.


Wu Ke

Garnier Carole

Shu Huazhong

Prostate segmentation on T2 MRI using Optimal Surface Detection

Keywords: Prostate, surface segmentation, T2 MRI

This paper deals with a T2 MRI prostate segmentation method. We assume to have an initial surface mesh obtained interactively or after a first rough segmentation. The surface of the prostate is then searched within the initial mesh neighborhood using the Optimal Surface Detection algorithm (OSD). This algorithm is based on the construction of a directed graph from the information obtained around the initial mesh. The optimal surface is then obtained by a graph cut. Three different cost functions for the graph have been explored, one based on the local gradient, another on a statistical model of shape and a third on a model of gradient profile. The parameters of this method have been tuned on 33 different T2 MRI volumes.

Introduction

The segmentation of the prostate surface in T2 MRI is a difficult task because several prostatic zones are not homogeneous and its surface has got little information (hyposignal of the capsule) with gray level values close to those of the neighboring organs surface (bladder and rectum). In this paper we will adapt a method previously developed for the semiautomated prostate segmentation on transrectal ultrasound images to the T2 MRI case [START_REF] Garnier | Prostate segmentation in HIFU therapy[END_REF]. This method uses the Optimal Surface Detection (OSD) [START_REF] Li | Efficient optimal surface detection: theory, implementation, and experimental validation[END_REF] to refine the prostate surface segmentation in the neighborhood of a first approximate surface obtained interactively or after a first rough segmentation. This paper will deal with the adaptation and the tuning of the OSD in the T2 MRI context. Two attempts to introduce statistical shape information will also be presented.

Material et methods

Optimal Surface Detection -OSD

As an initial stage we assume to have a first approximative prostate surface mesh. This initial mesh can be defined semi-automatically (e.g. an ellipsoid adjusted on 6 initial points on the boundary of the prostate volume -apex, base and 4 points within the central axial plane-defined interactively by a radiologist [START_REF] Garnier | Prostate segmentation in HIFU therapy[END_REF]) or after a first more automated rough segmentation.

The principle of the OSD algorithm is to consider the segmentation problem as the computation of a minimum s-t cut in an oriented weighted graph [START_REF] Li | Efficient optimal surface detection: theory, implementation, and experimental validation[END_REF]. The graph is built as follows. Each vertex i of the initial mesh serves as a graph node. From each vertex i a column of nodes is created along the surface normal N i and distributed on both side of the surface (Fig. 1). The search distances inside and outside the surface are set respectively to InD i and OutD i . Two auxiliary nodes s and t (source and sink respectively in Boykov's terminology [START_REF] Boykov | An experimental comparison of mincut/max-flow algorithms for energy minimization in vision[END_REF]) will allow the volume segmentation into two classes "inside" and "outside" of the surface.

These nodes are connected by links. The weights of the links will be used to encode different information (local image-based surface energy, shape constraints, gradient profile, ...). Several kind of links will be defined:

1) Intra-column links of infinite weights connect together two adjacent nodes of a same column from outside to inside. This construction insures that the surface intersects each column on only one of its nodes.

2) Inter-columns links of infinite weights will be used as a regularity constraint. One node N i,c of column i and position c on the column is connected to the node N j,max(0,c-∆ i,j ) of column j, with ∆ i,j a smoothness constraint (Fig. 1). The purpose of this constraint is to ensure that a surface can only move from distance less than ∆ from one column to another. ∆ will be adapted according to the distance between two neighbor vertexes p i and p j by ∆ i,j = ⌊δ p i -p j ⌋ with δ a constant.

3) The surface local energy is encoded as follows. In a first step a cost function C based on local image features and/or shape probability and/or some gradient profile model is defined. This function is then used to assign a weight w i,c to each node (with c = 0 the lowest node in a column):

w i,c = C i,c , if c = 0 C i,c -C i,c-1 , otherwise
The nodes with negative weights will be linked to the source s with a weight |w i,c | and the nodes with positive (or null) weights to the sink t with a weight w i,c .

4) The optimal surface is finally obtained with a min-cut/max-flow graph cut algorithm [START_REF] Boykov | An experimental comparison of mincut/max-flow algorithms for energy minimization in vision[END_REF].

Cost function

The definition of the cost function is the key point of the OSD algorithm. Three different classes of cost functions have been tested based on: 1) local image information (gradients), 2) a shape probability function and 3) a gradient profile model.

Image-based feature: gradients. The image-based surface energy is based on the smoothed gradient of the volume

V : E(x, y, z) = ∇(G σ * V (x, y, z)) = G(x, y, z)
with G σ a Gaussian filter with standard deviation σ. We tested 3 different cost functions: POS GDT promoting positive gradients (low to high value transitions which means the same direction as the surface normal):

C i,c = -G(x i,c , y i,c , z i,c ) s i,c +1 2 with s i,c = G(x i,c ,y i,c ,z i,c )•N i G(x i,c ,y i,c ,z i,c
) . EG GDT promoting negative gradients:

C i,c = -G(x i,c , y i,c , z i,c ) |s i,c -1| 2
. ALL GDT promoting any directions:

C i,c = -G(x i,c , y i,c , z i,c ) .

Shape probability function.

A shape model is first constructed from a training set of prostate contours manually extracted by an urologist on a set of prostate volumes (see § 3.1 for the data): 1) For each volume, a surface is built from the contours. For this, a local prostate-based reference system is first defined (the first axis is aligned from the apex to the base, the reference origin is set as the middle of the apex-base segment, the 2 other axes are adjusted to the prostate surface). Then the expert surface is regularly resampled within this local reference system in order to form a regular mesh.

2) The meshes of all the surfaces of the training set are aligned using the generalized Procustes algorithm [START_REF] Gower | Generalized Procrustes analysis[END_REF] implemented in VTK. This algorithm gives on one hand an average shape and on the other all the forms aligned on the average shape. A probability map is then estimated from the aligned shapes: a) a 3D distance map is build on each shape; b) the distance map is converted into a Gaussian distance map: exp -d 2 2σ 2 m with σ m a constant; and c) the probability map is deduced by averaging the Gaussian maps of all the aligned shapes. For the segmentation of a new prostate the probability map must be adjusted on the data. The model average shape is aligned onto the initial mesh of the data to segment. The estimated rigid transformation matrix (translation, rotation scaling) is finally applied on the probability map. The node cost value is the gradient-based cost multiplied by the probability at the node voxel position. Cost functions that add the probability to the gradient-based cost were also tested. As they lead to worse performance, results were not reported.

Gradient profile model. The idea is to build a statistical model of the gradient profiles variation from a learning set and then use this model to assign some costs [START_REF] Cootes | Statistical models of appearance for medical image analysis and computer vision[END_REF]. Surfaces of the training set are first aligned as described previously. At this stage, we have K shapes each sampled on N vertexes. N gradient profile models are then built as described in [START_REF] Cootes | Statistical models of appearance for medical image analysis and computer vision[END_REF]. A gradient profile g ij is estimated at each sample i (i = 1, . . . , N ) in each training surface j (j = 1, . . . , K). Each profile consists in the energy of 10 voxels on both sides of the surface. For each i we have a profile vector g i which can be normalized as:

g i → 1 j |g ij | g i . If
we assume that all profiles are distributed as a multivariate Gaussian, the appearance on one sample can be modeled by the mean profile ĝi and the covariance matrix S gi . The quality of the fit of a new profile g si to the model can be estimated by the Mahalanobis distance: f (g si ) = (g siĝi ) T S gi -1 (g si -ĝi ). For the segmentation of a new prostate, the initial mesh must be aligned on its local prostate based reference system and resampled in order to ensure correspondence to the model vertices. Gradients profiles g i,c are computed for each voxel that corresponds to a node N (i, c) in the graph. The node cost value is then defined as: C i,c = f (g i,c ) -max(f (g s )) where max(f (g s )) is the maximum Mahalanobis distance over all the computed profiles in the current image.

Results

Data and evaluation protocol

The algorithms were tested on 33 T2 axial MRI volumes from 5 different devices. We had images of different qualities and image size ranged from 256 × 256 × 24 to 720 × 720 × 20 voxels with pixel size ranging from 0.391 to 0.781mm and axial slice thicknesses from 2.63 to 4.50mm. On each of the 33 volumes, an urologist delineates manually the contours of the prostate. These contours will serve as ground truth or to build the shape and the gradient profile model. The MRI 3D images were interpolated along z before segmenting the prostate.

The performance of a specific segmentation scheme is evaluated with the Volume Overlap metric: VO = V M ∩V A V M ∪V A with V M and V A respectively the expert (manual) and the OSD (automatic) segmented volume. The mean in % and the standard deviation (std ) of all the 33 VOs serve as global performance metrics. However, it has to be noticed that this evaluation scheme is just one indicator of a performance. Only one manual segmentation per volume from one expert has been performed. Other studies with several segmentations per volume from one or several experts showed relative huge intra-and interobserver variabilities [START_REF] Fiorino | Intraand inter-observer variability in contouring prostate and seminal vesicles: implications for conformal treatment planning[END_REF]. So the results of our VO metric should be used only with caution and should only be considered as a relative performance measurement which gives high performance increases in some cases (e.g. cases 8, 11, 25) but also really important performance decreases (e.g. case 10).

Graph parameters setting

The graph parameters setting is performed on the classical OSD algorithm with the gradient-based cost function. Five parameters have to be tuned: InD and OutD the search distances inside and outside the surface, δ from the hard shape constraint, σ the Gaussian filter standard deviation and the promoted gradient direction (POS GDT, NEG GDT or ALL GDT). We tuned sequentially δ, σ, InD and OutD and the gradient direction. For a specific parameter, the value which gives the best performance is kept for the next parameter tuning. The optimal setting after our exhaustive search for the several parameters variations was: δ = 1/2, σ = 1mm, InD = -5mm, OutD = 15mm and the ALL GDT gradient direction. This setting gives a mean VO of 77.2% with a standard deviation of 5.5%.

The VO of 77.2% could be considered as relatively low. However, it should be noticed that due to the high interslice spacing of T2 MRI, the prostate is generally described on only 6 or 7 slices and that it is relatively difficult even for the radiologist or urologist to define the prostate shape in the planes tangent to the organ near the base or the apex. We also showed in a previous study [START_REF] Garnier | Prostate segmentation in HIFU therapy[END_REF] that the VO obtained between the prostate extracted on 3D US images (which had a better interslice spacing as T2 MRI) by two different experts or twice by the same expert presented a mean VO value of around 83%. Even if the score obtained on T2 MRI could not be compared to this obtained on the US volume, an analog behavior should be considered.

Shape probability and gradient profile model

In order to evaluate these algorithms, the training set used for the segmentation of a specific prostate volume was built from the expert manual segmentations of all the other 32 volumes. We first tuned the shape probability cost function based method on 7 volumes in order to find the optimal σ m parameter (σ m = 8).

Fig. 2 compares the VO obtained by the gradient based, shape probability and gradient profile model for all the 33 cases. The shape probability increases slightly the overlap in some cases but globally gives less performing results. A same conclusion can be drawn with the gradient profile modelbased cost function which gives some high performance increasing in some cases (e.g. cases 8, 11, 25) but also some really performance decreases (e.g. case 10). Some of the assumptions like a Gaussian distribution of the profiles are not always pertinent. Also because T2 MRIs of diseased prostates have a high variability in texture the gradient search (or model) has difficulty to separate the prostate boundary from its surroundings or its inside part.

Conclusion

The shape probability and gradient profile model can in some cases increase the performance of the segmentation but also can lead to very poor results. Such algorithms can be interesting in a segmentation scheme where the end-user (the urologist) can choose the best results in his sense between the several intermediary results obtained using the different cost functions. This work is part of the French MULTIP project supported by an ANR Grant (ANR-09-TECS-011-06).

Figure 1 :

 1 Figure 1: Graph construction.

Figure 2 :

 2 Figure 2: Comparison of local gradient-, shape probability-and gradient profile modelbased cost functions performance. Volume overlap metric for all the 33 volumes.

List of Figures


Sparse Sampling in Phase Space: Theory and Applications Ayush Bhandari * (ayush@MIT.edu) * Massachusetts Institute of Technology, 75 Amherst St, Cambridge, MA 02139, USA Abstract-The question that we seek to answer is: Given a continuous time sparse signal (in primal domain), is it possible to sample and reconstruct this sparse signal based on its characterization in a dual domainthe one which generalizes a number of already known transformations?

Phase Space transformations like the Linear Canonical Transform parametrically generalize a number of interesting transformations including the Fourier, Fresnel and fractional Fourier Transform among others. In this work, we develop a sampling and reconstruction methodology for sampling a sparse signal by analyzing their properties in Phase Space which is backward compatibe with several known transformations. This problem, as will be seen, can be recast as a parameter estimation problem that has its roots in Prony's method.

Given the general scope of analysis of sparse signals in transform domain, this work has applications in areas of spread spectrum methods for imaging, phase retrieval and computational photography.

Phase Space transformations are theoretically appealing for they generalize a number of existing transformations like the Fourier transform among others. One such construction includes the case of Linear Canonical Transform (LCT). The LCT (parametrized by

Λ Λ Λ ≡ a b c d with ad -bc = 1) of a signal, say x(t), is given by x Λ Λ Λ (ω) = LCT Λ Λ Λ {x} or, x Λ Λ Λ (ω) def = x (•) , φ Λ Λ Λ (•, ω) for b = 0, √ de j 1 2 (cdω 2 ) x(dω) for b = 0, (2)
where φ Λ Λ Λ (•, ω) is the phase space kernel defined in (1) of Table I and g, h denotes the inner-product in L2-sense. For Λ Λ Λ ∈ R 2 , the LCT is a unitary transformation and generalizes many transformations (Table I). Results applicable to the LCT domain can easily be extended to any of its special cases. The Fourier transform and the FrFT are the most notable byproducts. The LCT is equipped with an interesting composition property:

LCT Λ Λ Λ 2 {LCT Λ Λ Λ 1 {x}} = LCT Λ Λ Λ 3 {x} ⇒ Λ Λ Λ2 • Λ Λ Λ1 = Λ Λ Λ3
that leads to the definition of inverse-LCT of function, which is equivalent to

x (t) = LCT Λ Λ Λ -1 { x Λ Λ Λ } or, x(t) = x Λ Λ Λ (•), φ Λ Λ Λ -1 (t, •) for b = 0, √ ae -j 1 2 (cat 2 ) x(at) for b = 0, (3)
where Λ Λ Λ -1 is the inverse of Λ Λ Λ. If x (t) is the approximation of x(t), then x (t)x(t) 2 L 2 = 0 whenever (ωmb) ωs 2 (Nyquist rate for LCT domain), where ωs is the sampling frequency. When Λ Λ ΛFT = 0 1 -1 0 , all the aforementioned results take form of Shannon's sampling theorem. In this paper, we are interested in sampling a stream of K Dirac impulses-a signal with structure:

x (t) def = K-1 k=0 c k δ (t -t k ), (4)
with weights and arbitrary shifts, {c k , t k }. By constructing an aperiodic version of Fourier Series expansion (cf. [START_REF] Bhandari | Shift-Invariant and Sampling Spaces Associated with the Fractional Fourier Domain[END_REF]) for LCT domain, it turns out that (4) has a representation of form:

x (t) = 1 τ e -j a 2b t 2 k n=+∞ n=-∞    K-1 k=0 c k • e j a 2b t 2 k a k e -jnω 0 t k u n k    p[n]-Sum of K-complex exponentials e jnω 0 t . ( 5
) TABLE I LCT AS A GENERALIZATION OF OTHER TRANSFORMATIONS φ Λ Λ Λ (t, ω) = 1 √ -j2πb exp -j 2b at 2 + dω 2 -2ωt (1)
LCT Parameters (Λ Λ Λ) Corresponding Transform

cos θ sin θ -sin θ cos θ = Λ Λ Λ θ Fractional Fourier Transform (FrFT) 0 1 -1 0 = Λ Λ Λ FT Fourier Transform (FT) 1 b 0 1
Fresnel Transform

1 b  1
Bilateral Laplace Transform

1 -b 0 1 , b ≥ 0 Gauss-Weierstrass Transform 1 √ 2 0 e -π/2 -e -π/2 1
Bargmann Transform Discretizing (5) in m-points results in the following system:

x (m) e j a 2b m 2 x = 1 τ n=+∞ n=-∞ K-1 k=0 a k u n k p e jnω 0 m ⇔ x = D m×n DFT p
which is reminiscent of the problem tackled by Prony. For a full rank system, D -1 DFT x yields p and then the parameters are estimated. In what remains, we will recast this as a sampling problem [START_REF] Bhandari | Shift-Invariant and Sampling Spaces Associated with the Fractional Fourier Domain[END_REF] where we use a pre-filter to sample (4) and reconstruct it using [START_REF] Nayar | Fast separation of direct and global components of a scene using high frequency illumination[END_REF].

The finite-rate-of-innovation theory [START_REF]Special issue on compressive sensing[END_REF] is a special case of our result (set Λ = Λ θ ). Similarly, our work has close connections with Phase Retrieval problem of P. Jaming [START_REF] Jaming | Uniqueness results for the phase retrieval problem of fractional Fourier transforms of variable order[END_REF]. These ideas have applications in computational imaging [START_REF] Barsi | TOF Range Imaging[END_REF]. Given that the sparse signal spectrum can be uniquely characterized in [START_REF] Nayar | Fast separation of direct and global components of a scene using high frequency illumination[END_REF], the problem of source separation of of bandlimited + impulsive signal (4) can be efficiently tackled. This has applications in computational imaging [START_REF] Nayar | Fast separation of direct and global components of a scene using high frequency illumination[END_REF] as well as in mitigating impulsive noise [START_REF] Wolf | Redundancy, the discrete fourier transform, and impulse noise cancellation[END_REF]. The chirp-kernel of LCT might explain the success of spread spectrum imaging method of Wiaux et al. in [START_REF] Wiaux | Spread spectrum for imaging techniques in radio interferometry[END_REF]. In context of the recent results in line of Mathematical Theory of Super-Resolution of Candes et al. [START_REF] Candes | Super-Resolution from Noisy Data[END_REF], our problem can be reformulated as an optimization problem (noisy case).


Guillaume Aupy

Mathieu Faverge

Yves Robert

Jakub Kurzak

Piotr Luszczek

Jack Dongarra

Implementing a Systolic Algorithm for QR Factorization on Multicore Clusters with PaRSEC

Keywords: algorithm, QR factorization, PaRSEC, exascale

This article introduces a new systolic algorithm for QR factorization, and its implementation on a supercomputing cluster of multicore nodes. The algorithm targets a virtual 3D-array and requires only local communications. The implementation of the algorithm uses threads at the node level, and MPI for inter-node communications. The complexity of the implementation is addressed with the PaRSEC software, which takes as input a parametrized dependence graph, which is derived from the algorithm, and only requires the user to decide, at the high-level, the allocation of tasks to nodes. We show that the new algorithm exhibits competitive performance with state-of-the-art QR routines on a supercomputer called Kraken, which shows that high-level programming environments, such as PaRSEC, provide a viable alternative to enhance the production of quality software on complex and hierarchical architectures.

Implémentation d'un algorithme systolique pour la factorisation QR sur des clusters multicores à l'aide de PaRSEC Résumé : Cet article présente un nouvel algorithme systolique pour la factorisation QR, ainsi que son implémentation sur un cluster de noeuds multicoeurs. L'algorithme a été conc ¸u pour un tore-3D virtuel et ne demande que des communications locales. L'implémentation de cet algorithme utilise des threads au niveau des noeuds, ainsi que MPI pour les communication inter-noeuds. La complexité de l'implémentation a été maîtrisée grâce à l'utilisation du logiciel PaRSEC, qui prend en entrée un graphe de dépendances paramétrisé, dérivé de l'algorithme, et ne laisse à l'utilisateur que le choix de l'allocation haut-niveau des tâches aux noeuds. Le nouvel algorithme s'avère aussi efficace que des routines QR à la pointe de l'art sur le super-ordinateur Kraken, montrant ainsi que l'environnement PaRSEC est une excellente alternative pour accroître la production de logiciels de qualité sur des architectures complexes et hiérarchiques.

Mots-clés : algorithmes, factorisation QR, PaRSEC, exascale

Introduction

Future exascale machines are projected be massively parallel architectures, with 10 5 to 10 6 nodes, each node itself being equipped with 10 3 to 10 4 cores. At the node level, the architecture is a shared-memory machine, running many parallel threads on the cores. At the machine level, the architecture is a distributed-memory machine. This additional level of hierarchy dramatically complicates the design of new versions of the standard factorization algorithms, that are central to many scientific applications. In particular, the performance of numerical linear algebra kernels are at the heart of many grand challenge applications, and it is of key importance to provide highly efficient implementations of these kernels to leverage the potential of exascale platforms. This article introduces a new systolic algorithm for QR factorization on clusters of multicore nodes. The main motivation is to enhance the state-of-the-art algorithms, that use tile kernels and several elimination domains per panel, that enforce the inter-node communication between neighbors only. The systolic algorithm targets a 3D torus, which is the underlying interconnection topology of the contemporary and up-coming HPC systems. For instance, Blue Gene/L is a 3D torus of size 64 × 32 × 32 [START_REF] Adiga | An overview of the BlueGene/L supercomputer[END_REF], Kraken, a Cray XT 5, is a 3D torus of size 25 × 16 × 24 [START_REF]Kraken machine size[END_REF]. In addition, the Cray XT3 and XT4 also are architectures based on a 3D torus [START_REF] Bhatele | Application-specific topology-aware mapping for three dimensional topologies[END_REF]. Our systolic algorithm uses a two-level allocation of tile rows to the faces of the torus cube, in order to restrict the reduction tree for each panel to only local communication.

Implementing such a complex algorithm with low-level primitives would require non-trivial and error-prone programming effort. However, using the PaRSEC software [START_REF] Bosilca | DAGuE: A generic distributed DAG engine for high performance computing[END_REF] has enabled us to implement, validate, and evaluate the algorithm on Kraken, within a few weeks of development. Although we use a high-level environment, we report competitive performance results with state-of-the-art QR routines, thereby showing that PaRSEC provides a viable alternative to enhance the production of quality software prototypes on complex hierarchical architectures.

The rest of the paper is organized as follows. Section 2 provides background information on QR factorization algorithms, and surveys state-of-the-art algorithms in the literature. Section 4 presents the new systolic algorithm, while Section 5 provides additional details of its implementation using the PaRSEC software. Section 6 presents experimental results obtained on the Kraken supercomputer. Finally, we close with concluding remarks in Section 7.

Background

Tiled-QR Factorization

Algorithm 1: Generic QR algorithm.

begin for k = 0 to min(m, n) -1 do for i = k + 1 to m -1 do elim(i, CurPiv(i, k), k) RR n°8390
Algorithm 2: Elimination elim(i, CurPiv(i, k), k).

begin (a)

With TS (Triangle on top of square) kernels GEQRT(CurPiv(i, k), k) TSQRT(i, CurPiv(i, k), k)

for j = k + 1 to n -1 do UNMQR(CurPiv(i, k), k, j) TSMQR(i, CurPiv(i, k), k, j) (b) With TT (Triangle on top of triangle) kernels GEQRT(CurPiv(i, k), k) GEQRT(i, k) for j = k + 1 to n -1 do UNMQR(CurPiv(i, k), k, j) UNMQR(i, k, j) TTQRT(i, CurPiv(i, k), k) for j = k + 1 to n -1 do TTMQR(i, CurPiv(i, k), k, j)
The general shape of a tiled QR algorithm for a tiled matrix of m × n tiles, whose rows and columns are indexed from 0, is given in Algorithm 1. Here i and k are tile indices, and operate on square b × b tiles, where b is the block size. Thus, the actual size of the matrix is M × N, with M ≡ m × b and N ≡ n × b. The first loop index, k, is the panel index, and elim(i, CurPiv(i, k), k) is an orthogonal transformation that combines rows i and CurPiv(i, k) to zero out the tile in position (i, k). Each elimination elim(i, CurPiv(i, k), k) consists of two sub-steps: (i) first in column k, tile (i, k) is zeroed out (or eliminated) by tile (CurPiv(i, k), k), which is called the pivot; and (ii) in each subsequent column j > k, tiles (i, j) and (CurPiv(i, k), j) are updated; all these updates are independent and can be triggered as soon as the elimination is completed. The algorithm is entirely characterized by its elimination list, which is the ordered list of all the eliminations elim(i, CurPiv(i, k), k) that are executed.

To implement an orthogonal transformation elim(i, CurPiv(i, k), k), we can use either TT kernels or TS kernels, as shown in Algorithm 2. TT kernels are needed to allow for several eliminator tiles in a given column, but are less efficient than TS kernels. More detailed information on the various kernels is provided elsewhere [START_REF] Dongarra | Hierarchical QR factorization algorithms for multi-core clusters[END_REF]. In a nutshell, a tile can have three states: square, triangle, and zero. Transitions between these states occur through the following kernels:

• GEQRT is the transformation of one square tile to a triangle tile,

• TSQRT(i, CurPiv(i, k), k) is the transformation of a square tile (tile i) into a zero tile, using a triangle tile (tile CurPiv(i, k)) at step k, • TTQRT(i, CurPiv(i, k), k) is the transformation of a triangle tile (tile i) into a zero tile using a triangle tile (tile CurPiv(i, k)) at step k.

Related Work

While the advent of multi-core machines is somewhat recent, there is a long line of papers related to tiled QR factorization. Tiled QR algorithms have first been introduced in Buttari et al. [START_REF] Buttari | Parallel tiled QR factorization for multicore architectures[END_REF][START_REF] Buttari | A class of parallel tiled linear algebra algorithms for multicore architectures[END_REF] and Quintana-Ortí et al. [START_REF] Quintana-Ortí | Programming matrix algorithms-by-blocks for thread-level parallelism[END_REF] for shared-memory (multi-core) en-RR n°8390

vironments, with an initial focus on square matrices. The sequence of eliminations presented in these papers is analogous to the prior work [START_REF] Sameh | On stable parallel linear systems solvers[END_REF], and corresponds to reducing each matrix panel with a flat tree: in each column, there is a unique eliminator, namely the diagonal tile.

The introduction of several eliminators in a given column has a long history [START_REF] Sameh | On stable parallel linear systems solvers[END_REF][START_REF] Modi | An alternative Givens ordering[END_REF][START_REF] Pothen | Distributed orthogonal factorization: Givens and Householder algorithms[END_REF][START_REF] Da Cunha | New parallel (rank-revealing) QR factorization algorithms[END_REF][START_REF] Demmel | Communication-avoiding parallel and sequential QR and LU factorizations: theory and practice[END_REF][START_REF] Langou | Computing the R of the QR factorization of tall and skinny matrices using MPI Reduce[END_REF]. For shared-memory (multi-core) environments, recent work advocates the use of domain trees [START_REF] Hadri | Tile QR factorization with parallel panel processing for multicore architectures[END_REF] to expose more parallelism with several eliminators while enforcing some locality within domains. A recent paper [START_REF] Bouwmeester | Tiled QR factorization algorithms[END_REF] introduces tiled versions of the Greedy algorithm [START_REF] Cosnard | Complexity of parallel QR factorization[END_REF] and the Fibonacci scheme [START_REF] Modi | An alternative Givens ordering[END_REF], it shows that these algorithms are asymptotically optimal.

There are recent efforts for distributed-memory environments. The algorithm of [START_REF] Agullo | QR factorization of tall and skinny matrices in a grid computing environment[END_REF] uses a hierarchical approach: for each matrix panel, it combines two levels of reduction trees. First, several local binary trees are applied in parallel, one within each node, and then a global binary tree is applied for the final reduction across nodes. Yet another implementation [START_REF] Song | Scalable tile communication-avoiding QR factorization on multicore cluster systems[END_REF] also uses a hierarchical approach, and it also uses a 1D block distribution. The main difference is that the first level of reduction is performed with a flat tree within each node. Note that the hierarchical algorithm (HQR) used previously [START_REF] Dongarra | Hierarchical QR factorization algorithms for multi-core clusters[END_REF] can be parametrized to implement this original algorithm [START_REF] Song | Scalable tile communication-avoiding QR factorization on multicore cluster systems[END_REF] as well as a more efficient variant with cyclic layout. The HQR algorithm [START_REF] Dongarra | Hierarchical QR factorization algorithms for multi-core clusters[END_REF] is the reference algorithm for multi-level clusters: it provides a flexible approach, and allows one to use various elimination trees (Flat, Binary, Fibonacci or Greedy) at each level. [START_REF] Bosilca | DAGuE: A generic distributed DAG engine for high performance computing[END_REF] The SYSTOLIC-3D algorithm Platform and data layout -We first detail the 3D torus architecture. Within a p × q × r 3D torus, processor P a,b,c has a direct communication link with processors P a-1 mod p,b,c , P a+1 mod p,b,c , with processors P a,b-1 mod q,c , P a,b+1 mod q,c , and with processors P a,b,c-1 mod r , P a,b,c+1 mod r . We have a m × n tile matrix. Tiles are mapped as follows: we use a two-level cyclic distribution for the rows (directions a and b in the torus) and a cyclic distribution for the columns (direction c in the torus). The mapping is defined formally as follows: proc P abc is assigned all tiles T t,s such that t ≡ b mod q, t-b q ≡ a mod p and s ≡ c mod r. We give an example of the two-level cyclic distribution for the rows in Figure 1a, for a matrix with 27 rows mapped onto a 3 × 3 × r torus.

General description -As stated in Section 2.1, a tiled-QR algorithm is entirely defined by the ordered list of eliminations. The algorithm eliminates the tiles using a hierarchical approach, using the 3D torus to minimize inter-processor communication contention. In order to do this, pivots should be propagated across physical links in the torus, and only to neighbor nodes, before each elimination. Figure 1 describes the elimination of the first column of the matrix.

Consider a given step k of the factorization. The k-th tile column is distributed across a face of the cube, i.e. a square of p × q processors (those whose third index is c 0 ≡ k mod r). Let dimension a be "horizontal" and dimension b be "vertical". There are three levels of elimination in the algorithm:

1. The first level of elimination corresponds to local tiles and uses TS kernels. There are pq pivots in this step, one for each processor in the square, and they correspond to rows numbered k, k + 1, • • • , k + pq -1. These pivots are used to eliminate all local tiles within each processor, hence they do not require any commu-nication across the square. We use a flat tree reduction for this step, but other elimination trees could be chosen freely. This first elimination level is illustrated in Figure 1a when k = 0. 2. The second level of elimination consists of concurrent flat trees along the vertical dimension, and uses TT kernels (see Figure 1b). There are p pivots for this level, namely the k th elements of rows k, k + q, • • • , k + q(q -1). Each of these pivots will sequentially eliminate the q -1 subsequent tiles, which are located in the corresponding grid column. 3. The third level of the elimination consists of a single flat tree along the horizontal dimension (see Figure 1c). There remains a single pivot, in row k, that sequentially eliminates with TT kernels the q -1 remaining tiles. At the end of step k, row number k will have been routed through at most p + q -2 physical communication links. The communication pattern is the same for the other faces of the cube. The whole algorithm is summarized in Algorithm 3. // Note that from now on, we do not need to use GEQRT anymore, all the remaining tiles are triangles. / * Vertical FlatTree * / for

l 2 = i 2 + 1 to d -1 do TTQRT V(k + (l 2 -i 2 ), k, k); for j = k + 1 to n -1 do TTMQR V(k + (l 2 -i 2 ), k, k, j); for l 2 = 0 to i 2 -1 do TTQRT V(k + d 2 + (l 2 -i 2 ), k, k); for j = k + 1 to n -1 do TTMQR V(k + d 2 + (l 2 -i 2 ), k, k, j); for l = k + d to k + d 2 -1 by d do for x = l + 1 to (l -i 2 ) + d -1 do TTQRT V(x, l, k); for j = k + 1 to n -1 do TTMQR V(x, l, k, j); for x = l -i 2 to l -1 do TTQRT V(x, l, k); for j = k + 1 to n -1 do TTMQR V(x, l, k, j); / * Horizontal FlatTree * / for x = k + d to k + d 2 -1 by d do TTQRT H(x, k, k); for j = k + 1 to n -1 do TTMQR H(x, k, k, j);

Implementation with PaRSEC

This section details the implementation of the SYSTOLIC-3D algorithm using PaRSEC.

With an infinite number of resources, the scheduling could follow a greedy heuristic: RR n°8390 the execution would progress as fast as possible.

The elimination list of the algorithm is the composition of the reduction trees at all of the different levels. All eliminators are known before the execution. Each component of an elimination is triggered as soon as possible, i.e., as soon as all dependencies are satisfied: first we have the elimination of the tile, and then the updates in the trailing panels. Note that the overall elimination scheme is complex, and mixes the elimination of tiles at all levels. However, with a fixed number of resources, it is necessary to decide an order of execution of the tasks, hence to schedule them: this is achieved through the PaRSEC environment.

Introduction to PaRSEC

PaRSEC is a high-performance fully-distributed scheduling environment for systems of micro-tasks. It takes as input a problem-size-independent, symbolic representation of a Direct Acyclic Graph (DAG) of tasks, and schedules them at runtime on a distributed parallel machine of multi-cores. Data movements are expressed implicitly by the data flow between the tasks in the DAG representation. The runtime engine is then responsible for actually moving the data from one machine (node) to another, using an underlying communication mechanism such as MPI. A full description of PaRSEC, and the implementation of classical linear algebra factorizations in this environment, is provided elsewhere [START_REF] Bosilca | DAGuE: A generic distributed DAG engine for high performance computing[END_REF][START_REF] Bosilca | Flexible development of dense linear algebra algorithms on massively parallel architectures with DPLASMA[END_REF].

To implement any QR algorithm in PaRSEC, it is sufficient to give an abstract representation of all the tasks (eliminations and updates) that constitute the QR factorization, and how data flows from one task to the other. Since a tiled QR algorithm is fully determined by its elimination list, it suffices to provide a function, that the runtime engine is capable of evaluating, and that computes this elimination list. The PaRSEC object obtained in this way is generic: when instantiating a PaRSEC QR factorization, the user simply gives the size of the platform (p × q × r), defining a new DAG at each instantiation. Note that this DAG is not fully generated: PaRSEC keeps only a parametric representation of the DAG in memory, and interprets symbolic expressions at runtime to explicitly represent only the ready tasks at any given time. This technique is similar to the Parametrized Tasks Graphs [START_REF] Bosilca | DAGuE: A generic distributed DAG engine for high performance computing[END_REF], and to HQR [START_REF] Dongarra | Hierarchical QR factorization algorithms for multi-core clusters[END_REF].

At runtime, task executions trigger data movements, and create new ready tasks, following the dependencies defined by the elimination list. Tasks that are ready to compute are scheduled according to a data-reuse heuristic: each core will try to execute close successors of the last task it ran under the assumption, that these tasks require data that was recently touched by the completed task. This policy is tuned by the user through a priority function: among the ready tasks for a given core, the choice is done by following a hint from this function. To balance the load between cores, tasks on the same node in the algorithm (residing on the same shared memory machine) are shared between the computing cores, and a NUMA-aware job stealing policy is implemented. The user is responsible for defining the affinity between data and tasks, and to distribute the data between the computing nodes. Thus, he defines which tasks should execute on which node, and he is responsible for this level of load balancing. In our case, the data distribution is the data layout given in Section 4. Since all kernel operations modify a single tile (or a tile and its reflectors, which are distributed in the same way), we chose the strategy "owner computes" for the tasks: tasks' affinity is set to the node that owns the data that is going to be modified, and the input data might need to be transferred from other nodes. Because HQR uses a 2D-processor grid, we use T nodes configured as a (pq)×r 2D grid.

• SYSTOLIC-2D is a variant of SYSTOLIC-3D where q is set to 1 and then runs on a 2D grid of size (pq) × r. We introduced it for the sake of comparison with the HQR variants -SYSTOLIC-2D can be viewed as yet another HQR variant with a new high-level reduction tree.

We compare all the previous algorithms that were implemented with PaRSEC with the following algorithms from the literature [START_REF] Kurzak | Virtual systolic array for QR decomposition[END_REF] on the very same hardware:

• SYSTOLIC-1D is the virtual systolic array decomposition [START_REF] Kurzak | Virtual systolic array for QR decomposition[END_REF]. As its name indicates, it targets a 1D-linear array of processors. Note that SYSTOLIC-1D has been implemented using a hand-written communication engine over MPI -not PaRSEC.

• HPL 4/3 N 3 is the virtual performance of the High Performance Linpack LU factorization using the flops count of QR: O( 4 3 N 3 ). • LIBSCI QR is the QR factorization from ScaLAPACK used in the Cray Scientific Library.

• HPL 2/3 N 3 is the High Performance Linpack LU factorization with the actual flops count of LU: O( 2 3 N 3 ). For each set of results, we ran the different algorithms ten times, and we take the average performance over all these executions.

M = N = 10, 368 M = N = 20, 736 M = N = 41, 472 T p × q × r T p × q × r T p × q × r T p × q × r T p × q × r T p × q × r 4
Our decision to include performance numbers for HPL's LU factorization might seem controversial due to the fundamental differences between the LU and QR factorization algorithms including their numerical properties, operation-count, and the runtime behavior. However, from the end-user perspective, both LU and QR solve a system of linear equations, both are backward stable, and only an explicitly stated rule [START_REF] Dongarra | The LINPACK benchmark: Past, present, and future[END_REF] prohibits QR from scoring the entrants to the TOP500 ranking. With this in mind, we include results for the LU factorization, and include the case when we pretend that LU performs as many Flops as QR: O( 43 N 3 ) (this may be simply treated as time-to-run comparison) as well as the case where we report the actual performance rate based on the actual amount of floating point operations LU: O( 23 N 3 ).

Performance Results

The first observation is that PaRSEC-based algorithms (SYSTOLIC-3D, SYSTOLIC-2D and all HQR variants) always perform better than LIBSCI QR and HPL 4/3 N 3 , the QR factorization algorithms from Kraken's standard software stack. This observation holds for all matrix sizes, and this makes our main point: owing to the PaRSEC system, we have been able to experiment with a variety of complex, hierarchical algorithms, without paying the price of lengthy development effort. Note, how SYSTOLIC-3D , HQR variants, and SYSTOLIC-1D compare with each other. SYSTOLIC-3D has approximatively the same efficiency as HQR-BINARY and HQR-GREEDY on all matrix sizes. For matrices of size M = N = 10, 368, HQR-FLAT is 54% more efficient on 1536 cores (≈ 1700 GFlops compared to ≈ 1100 GFlops), and SYSTOLIC-1D (≈ 1900) GFlops is 73% more efficient. This difference increases with the size of the matrix: for M = N = 41, 472, HQR-FLAT reaches ≈ 16000 GFlops, and SYSTOLIC-1D reaches ≈ 21000 GFlops on 23868 cores, where SYSTOLIC-3D is bound by ≈ 10600 GFlops (half the performance of SYSTOLIC-1D). However, for M = N = 41, 472, and with a small number of cores, SYSTOLIC-3D performs better than SYSTOLIC-1D.

As mentioned earlier, it is infeasible to guarantee, that the assignment of Kraken nodes from our batch queue submissions form a true 3D torus. We expected the constraints to be less stringent when using a 2D torus, and this turns out quite true: SYSTOLIC-2D, the implementation of SYSTOLIC-3D on a 2D torus, performs very well for all matrix sizes, and is the best algorithm for larger matrices of size M = N = 41, 472.

In this article, we have presented a systolic QR factorization algorithm, SYSTOLIC-3D, which aims to minimize the amount of communication in the reduction trees. We have shown that mapping this systolic algorithm onto a 3D torus leads to a competitive factorization kernel with strong scaling capabilities. As of today, the main limitation to fully validate the experiments is the lack of possibility to reserve an actual 3D torus architecture on the Kraken supercomputer. Still, the performance of the new algorithm, together with its 2D counterpart are very encouraging. Both versions dramatically outperform LIBSCI QR and HPL 4/3 N 3 , the vendor QR factorization implementations on Kraken, and also HPL 2/3 N 3 , the widely-used LU factorization routine (despite its favorable flop count). This last observation fully demonstrates the usefulness of the PaRSEC system, which has enabled us to experiment with complex, hierarchical QR algorithms, without paying the price of lengthy and complex development effort of distributed memory software engineering.

Figure 1 -Algorithm 3 :

 13 Figure 1 -Elimination in the first panel (panel 0) of all tiles below the diagonal (rows 1 to 27) on a 3 × 3 processor square (face 0 of the 3D-torus).

1 .

 1 HQR-FLAT uses the FLATTREE reduction; 2. HQR-FIBO uses the FIBONACCI reduction; 3. HQR-BINARY uses the BINARYTREE reduction; 4. HQR-GREEDY uses the GREEDY reduction.

Figure 2 -

 2 Figure 2 -Performance of the various algorithms for different problem sizes.

Table 1 -

 1 Partition of the nodes into a 3D torus for each matrix size and each total number of nodes T .

	2 × 2 × 1 52 6 × 3 × 3	16 4 × 2 × 2 210 6 × 5 × 7	64	4 × 4 × 4	840 10 × 6 × 14
	12 3 × 2 × 2 80 5 × 4 × 4	48 4 × 4 × 3 320 8 × 5 × 8	192 8 × 4 × 6 1232 11 × 7 × 16
	18 3 × 3 × 2 96 6 × 4 × 4	80 5 × 4 × 4 405 9 × 5 × 9	336 7 × 6 × 8 1632 12 × 8 × 17
	28 5 × 2 × 3 128 8 × 4 × 4	112 6 × 4 × 5 486 9 × 6 × 9	480 8 × 6 × 10 1989 13 × 9 × 17
	42 7 × 2 × 3	168 7 × 4 × 6	648 9 × 6 × 12

RR n°8390

Acknowledgements Yves Robert is with Institut Universitaire de France. This work was supported in part by the ANR RESCUE project. It was also supported in part by grant #SHF-1117062: "Parallel Unified Linear algebra with Systolic ARrays (PUL-SAR)" from the National Science Foundation (NSF). The authors would like to thank the National Institute for Computational Sciences (NICS) for a generous time allocation on the Kraken supercomputer. Finally, the authors would also like to thank Thomas Hérault for his help with PaRSEC.

Implementation Details

The implementation of SYSTOLIC-3D in PaRSEC involves limited effort compared with other software strategies, that we are aware of. We only implemented a few functions that are used by PaRSEC to generate the dependency graph. They depend on the current elimination step k as follows:

1. CurPiv(i, k), returns the pivot to use for the row i at step k; 2. NextPiv(pivot, k,start), returns the next row which will use the row "pivot" as a pivot in step k after it has been used by row "start"; 3. PrevPiv(pivot, k,start), returns the previous row which used the row "pivot" as a pivot in step k before it has been used by row "start"; We have decomposed each one of these functions in two sub-functions: (i) a lowlevel function, which takes all the TS operations into account, and which calls the local FlatTree because operations are local to each node; and (ii) a high-level function, which takes all the TT operations into account, and where the pivot will "move" across the architecture. Using these functions, PaRSEC is able to construct a dependency graph between the different tiles in order to run the algorithm as efficiently as possible.

Experimental Evaluation

In this section, we report experimental results obtained on Kraken. We compare the SYSTOLIC-3D algorithm with a number of competing implementations such as vendor library routines and recent algorithms from literature.

Experimental Setup

All runs were done on the Kraken supercomputer at the National Institute for Computational Science [START_REF]Kraken machine size[END_REF]. The Kraken machine is a Cray XT5 system operated by the University of Tennessee and located in Oak Ridge, Tennessee, U.S.A. The entire system consists of 9048 computing nodes. The experiments presented here used up to 1989 nodes, which is about one fifth of the machine. Each node contains two 2.6 GHz sixcore AMD Opteron (Istanbul) processors, 16 GB of memory and the Cray SeaStar2+ interconnect.

We have compared SYSTOLIC-3D with several state-of-the-art algorithms, using three matrix sizes: (i) small matrices, of size M = N = 10, 368; (ii) medium matrices, of size M = N = 20, 736; and large matrices, of size M = N = 41, 472. Here is the list of the algorithms used for comparison:

• SYSTOLIC-3D is the algorithm described in this paper. Table 1 shows the 3D grid configuration (p, q, r) used for each matrix size (M = N) and for each total number of nodes T , where T = p × q × r. Note, that there is no guarantee, that the nodes assigned to the experiment will indeed form the desired 3D torus. They can be scattered across the machine. To the best of our knowledge, the only way to guarantee that assigned nodes indeed form a 3D torus would be to reserve the entire Kraken machine: something beyond our capabilities.

• HQR is the hierarchical QR factorization algorithm [START_REF] Dongarra | Hierarchical QR factorization algorithms for multi-core clusters[END_REF], which was also implemented using the PaRSEC software. We compare several variants of HQR, which all use the same FLAT-TREE low-level reduction tree, but which use different high level (or distributed) reduction trees [START_REF] Dongarra | Hierarchical QR factorization algorithms for multi-core clusters[END_REF]: