Raphael Bleuse

Safia Kedad-Sidhoum

Florence Monna

G R Égory

Mouni É2

Denis Trystram

SCHEDULING INDEPENDENT TASKS ON MULTI-CORES WITH GPU ACCELERATORS

Keywords: Scheduling, Approximation algorithms, Parallel heterogeneous systems

More and more computers use hybrid architectures combining multi-core processors and hardware accelerators like GPUs (Graphics Processing Units). We present in this paper a new method for scheduling efficiently parallel applications with m CPUs and k GPUs, where each task of the application can be processed either on a core (CPU) or on a GPU. The objective is to minimize the maximum completion time (makespan). The corresponding scheduling problem is NP-hard, we propose an efficient approximation algorithm which achieves an approximation ratio of 4 3 + 1 3k . We first detail and analyze the method, based on a dual approximation scheme, that uses dynamic programming to balance evenly the load between the heterogeneous resources. Then, we present a faster approximation algorithm for a special case of the previous problem, where all the tasks are accelerated when affected to GPU, with a performance guarantee of 3 2 for any number of GPUs. We run some simulations based on realistic benchmarks and compare the solutions obtained by a relaxed version of the generic method to the one provided by a classical scheduling algorithm (HEFT). Finally, we present an implementation of the 4/3-approximation and its relaxed version on a classical linear algebra kernel into the scheduler of the xKaapi runtime system.

Introduction

Most of the computing systems available today include parallel multi-core chips sharing a large memory with additional hardware accelerators [START_REF] Lee | Debunking the 100x gpu vs. cpu myth: an evaluation of throughput computing on cpu and gpu[END_REF]. There is an increasing complexity within the internal nodes of such parallel systems, mainly due to the heterogeneity of the computational resources. In order to take advantage of the benefits offered by heterogeneity in terms of performance, effective and automatic management of the hybrid resources will be more and more important for running any applications. These new hybrid architectures have given rise to new scheduling problems consisting in allocating and sequencing the computations on the different resources such that a given objective is optimized.

The main challenge is to create adequate generic methods and software tools that fulfill the requirements for optimizing the performances. In the field of High Performance Computing (HPC), one of the most studied optimization problem is the minimization of the maximum completion time (makespan) of the schedule, which is the objective addressed in this paper. Some Polynomial Time Approximation Schemes (PTAS) exist for problems of minimizing the makespan on heterogeneous processors [START_REF] Hochbaum | A polynomial approximation scheme for scheduling on uniform processors: using the dual approximation approach[END_REF][START_REF] Bonifaci | Scheduling unrelated machines of few different types[END_REF], but their running times make them impractical for solving scheduling problems on actual computing platforms.

In the field of parallel processing, there exist a huge number of papers dealing with implementations of ad hoc algorithms using GPUs or hybrid architectures. They expand over several aspects of parallelism from operating system, runtime, application implementation or languages. However, only few of them focus on the intermediate problem of scheduling on hybrid platforms [START_REF] Pinel | Solving very large instances of the scheduling of independent tasks problem on the gpu[END_REF]. Most of the works in the literature consist in studying the gains and performances of parallel implementation of some specific numerical kernels [START_REF] Agullo | Qr factorization on a multicore node enhanced with multiple gpu accelerators[END_REF][START_REF] Song | Enabling and scaling matrix computations on heterogeneous multi-core and multi-gpu systems[END_REF], specific applications like multiple alignments of biological sequences [START_REF] Boukerche | A hardware accelerator for the fast retrieval of dialign biological sequence alignments in linear space[END_REF], or molecular dynamics [START_REF] Phillips | Adapting a message-driven parallel application to gpuaccelerated clusters[END_REF]. The existing scheduling algorithms and tools are usually not well-suited for general purpose applications since the internal hardware organization of a GPU highly differs from a CPU and thus, the GPU should be considered as a new type of resource in order to determine efficient approaches. Scheduling is usually done on a case by case basis and often offers good performances, however, it lacks high-level mechanisms that provide transparent and efficient schedules for any application. Some actual runtime systems include the basic mechanisms for developing scheduling algorithms like OMPSS [START_REF] Bueno | Productive programming of gpu clusters with ompss[END_REF], StarPU [START_REF] Augonnet | Starpu: A unified platform for task scheduling on heterogeneous multicore architectures[END_REF] or xKaapi [START_REF] Gautier | Xkaapi: A runtime system for dataflow task programming on heterogeneous architectures[END_REF]. Several scheduling algorithms have been implemented on top of these systems. An online algorithm with a performance guarantee [START_REF] Chen | Online scheduling on a cpu-gpu cluster[END_REF] has recently been developed for CPU-GPU platforms, but there is no performance guarantee for any offline problem on these systems.

Our objective within this work is to propose a new method for scheduling independent tasks on hybrid CPU-GPU architectures designed for HPC. The considered input is a set of independent sequential tasks whose execution times are known. This hypothesis is realistic, since some computing systems such as StarPU have a module which estimates the execution times at compile time. The method that we propose in this work determines the allocation and schedule of the tasks to the computing units, CPUs and GPUs. We present and analyze in detail this methodology for the case of scheduling n sequential tasks on m cores (CPUs) and k GPUs. This leads to an efficient approximation algorithm which achieves a ratio 4 3 + 1 3k + ǫ using dual approximation [START_REF] Hochbaum | Using dual approximation algorithms for scheduling problems theoretical and practical results[END_REF] with a dynamic programming scheme. The computational cost of this algorithm is in O n 2 k 3 m 2 per step of dual approximation. We also present an approximation algorithm for the more specific problem of scheduling a set of taskson the same hybrid platform that are all accelerated when processed on a GPU, although they can have different accelerations. This algorithm also relies on the dual approximation technique and achieves a 3 2 performance ratio for a time complexity in O (mn log n). As the method with a 4 3 approximation ratio is rather costly, in the perspective of an integration into actual runtime systems, we derive a relaxed algorithm and compare it experimentally with one of the most popular algorithms (HEFT [START_REF] Topcuoglu | Performance-effective and low-complexity task scheduling for heterogeneous computing[END_REF]), as well as with the approximation algorithm presented for the case where all the tasks are accelerated on GPUs.

The outline of the paper is as follows. A formal description of the scheduling problem with k GPUs is provided in Section 2, and some related works are presented. We propose in Section 3 a new approach for solving the problem. We propose in Section 4 an algorithm for the scheduling problem with k GPUs where all tasks are accelerated when processed on GPU. We report the results of experiments in Section 5 where a relaxed version of the method presented in Section 3 and the algorithm from Section 4 are compared to the classical HEFT algorithm on simulations built from realistic workloads. The experimental analysis shows that the proposed methods have a more stable behavior than HEFT for a similar performance. Then, an implementation of both the generic method and its relaxed version on a real run-time system have been realized and tested on a classical Linear Algebra kernel. Finally, a synthesis is provided in Section 6, which opens some perspectives. Some of this work was presented at the HeteroPar 2013 conference [START_REF] Kedad-Sidhoum | Scheduling independent tasks on multicores with gpu accelerators[END_REF], but the approximation algorithm for the specific problem of scheduling a set of tasks that are all accelerated when processed on a GPU is an entirely new contribution, and the experiments were extended in comparison to the original paper, where only partial simulations for the main algorithm were presented. Practically, the cost of the 4 3 algorithm is too expensive since the number of independent tasks is too low in the considered kernel. The relaxed 2-approximation is shown to be as good as HEFT, but provides a more stable behavior and a lower volume of communications.

Problem Definition and Related Works

We consider a multi-core parallel platform composed of m identical CPUs and k identical GPUs. The m CPUs are considered independent from the GPUs that are commanded by some extra driving CPUs, not mentioned here because they do not execute any task. An application is composed of n independent tasks denoted by T 1 , . . . , T n . The tasks are considered sequential, meaning that they only are processed on one processor of any type, CPU or GPU. Each of these tasks has two processing times depending on which type of processor it is allocated to. The processing time of task T j is denoted by p j if T j is processed on a CPU and p j if it is processed on a GPU. The acceleration factor of task T j will be given by the ratio pj pj . We assume that both processing times of a task are known in advance as it is commonly admitted. An accurate estimation can be obtained at compile time for regular numerical applications in HPC. The objective here is to minimize the makespan C max of the whole schedule. It is defined as the maximum completion time of the last finishing task on both CPUs and GPUs, C max = max(C CP U max , C GP U max ). The problem will be denoted as (P m, P k) || C max .

If both processing times are equal (p j = p j ) for all j = 1, . . . , n, the problem (P 1, P 1) || C max is equivalent to the classical P 2 || C max problem, which is NP-hard [START_REF] Garey | Bounds for multiprocessor scheduling with resource constraints[END_REF]. Thus, the problem of scheduling with GPUs is also NP-hard. Our objective is to propose efficient approximation algorithms with a performance guarantee. Recall that for a given problem, the approximation ratio ρ A of an algorithm A solving this problem is defined as the maximum over all the instances I of the ratio f (I) f * (I) where f is any minimization objective and f * is its optimal value [START_REF] Hochbaum | Using dual approximation algorithms for scheduling problems theoretical and practical results[END_REF].

The problem considered here is a new problem, harder than the classical scheduling problem on uniform machines Q || C max , but better results can be expected if the problem is considered on its own rather than as an unrelated scheduling problem R || C max [START_REF] Blazewicz | Handbook on Scheduling, From Theory to Applications[END_REF]. Lenstra et al. [START_REF] Lenstra | Approximation algorithms for scheduling unrelated parallel machines[END_REF] proposed a PTAS for the problem R || C max with running time bounded by the product of (n + 1) m/ǫ and a polynomial of the input size. Let us notice that if m is not fixed, then the algorithm is not fully polynomial. The authors also proved that unless P = N P , there is no polynomialtime approximation algorithm for R || C max with an approximation factor less than 3 2 and they presented a 2-approximation algorithm. This algorithm is based on rounding an optimal solution of the preemptive version of the problem obtained by an integer linear program. Shmoys and Tardos [START_REF] Shmoys | An approximation algorithm for the generalized assignment problem[END_REF] generalized this technique to obtain the same approximation factor for the generalized assignment problem. Furthermore, they generalized the rounding technique to hold for any fractional solution. Recently, Shchepin and Vakhania [START_REF] Shchepin | An optimal rounding gives a better approximation for scheduling unrelated machines[END_REF] introduced a new rounding technique which yields an improved approximation factor of 2 -1 m for R || C max for a similar time complexity as [START_REF] Lenstra | Approximation algorithms for scheduling unrelated parallel machines[END_REF]. According to our knowledge, this is so far the best approximation result for this problem. However, the prohibitive computational cost of these algorithms prevents their usage on actual computing platforms.

It is worth noticing that if all the tasks of the problem have the same acceleration on the GPUs, the problem reduces to the classical Q || C max problem, with two machine speeds. For Qm || C max , Friesen [START_REF] Friesen | Tighter bounds for lpt scheduling on uniform processors[END_REF] proved that the approximation ratio of the well-known Longest Processing Time (LPT) scheduling policy satisfies 1.52

C LP T max /C * max 5/3.
The first PTAS for Q || C max was given by Hochbaum and Shmoys [START_REF] Hochbaum | A polynomial approximation scheme for scheduling on uniform processors: using the dual approximation approach[END_REF]. The overall running time of the algorithm is O((log m + log(3/ǫ))(m/ǫ)(n/ǫ)1/ǫ). However, these solving methods would only work for specific instances of the problem of scheduling on hybrid platforms, where the acceleration factors pj pj would be equal to a constant, which is not the case in practice.

Another direction is to consider the problem of scheduling on unrelated machines of few different types. Indeed, the R || C max reference problem can be refined to better fit the constraints of the hybrid platforms. Bonifaci and Wiese [START_REF] Bonifaci | Scheduling unrelated machines of few different types[END_REF] presented a PTAS to solve a scheduling problem with unrelated machines of few different types. The tools used in their solving method are somewhat similar to the ones used for solving R || C max , and the rounding phases of the algorithm require a significant amount of time, raising the time complexity of the algorithm to an impractical level, even when only two types of machines are considered, as it is for (P m, P k) || C max . There is a need to consider other algorithms than these PTAS to design algorithms that could be implemented on actual platforms. A PTAS with a reasonable time complexity has been developed for the online version of the problem of the assignment of sporadic tasks on hybrid platforms [START_REF] Nélis | A ptas for assigning sporadic tasks on two-type heterogeneous multiprocessors[END_REF]. However, an offline version of the problem with non-periodic tasks has not been studied and the algorithm cannot be trivially extended to the problem (P m, P k) || C max . On another side, Imreh [START_REF] Imreh | Scheduling problems on two sets of identical machines[END_REF] presented different greedy algorithms for the problem of scheduling on two sets of identical machines, with varying approximation ratios including 2 + m-1 k and 4 -2 m (for m CPUs and k GPUs). These algorithms are fast enough for being implemented in modern platforms, nevertheless the approximation ratios of these algorithms are quite high since usually the number of GPUs is much lower than the number of CPUs.

From a practical perspective, the scheduling strategy is a key point for the performance of an application. Tuning scheduling algorithms for a specific case (problem and computer architecture) is common. Fast heuristics without performance guarantee are often used on computing platforms, time efficiency being the crucial factor. However, simple strategies are not sufficient to guarantee the performance for more general cases potentially far from the specific one. The performance portability is difficult to achieve when the number of CPUs and GPUs varies or when the speedup of the various parts of the application is evolving during the execution.

Our objective within this work is to build a bridge between purely theoretical algorithms with good performance guarantees and practical low cost heuristics. Thus, we propose a tradeoff solution with a provable performance guarantee and a reasonable time complexity.

A New Algorithm

3.1. Rationale of the solving method. The proposed algorithm is based on the dual approximation technique [START_REF] Hochbaum | Using dual approximation algorithms for scheduling problems theoretical and practical results[END_REF]. A g-dual approximation algorithm for a generic problem takes a real number λ (guess) as an input, assumes that there exists a schedule of length at most λ and either delivers a schedule of makespan at most gλ, or answers correctly that there exists no schedule of length at most λ. A binary search is used to try different guesses to approach the optimal makespan as follows: we first take an initial lower bound B min and an initial upper bound B max of the optimal makespan. We start by solving the problem with a λ equal to the average of these two bounds and then the bounds are updated as follows:

• If the algorithm returns "NO", then λ becomes the new lower bound.

• If the algorithm returns a schedule of makespan at most gλ, then there exists a schedule of makespan at most λ and λ becomes the new upper bound. The number of iterations of the binary search is bounded by log 2 (B max -B min ) . Hence, a g-dual approximation algorithm can be converted, by bisection search, in a g(1 + ǫ)-approximation algorithm with a similar running time.

We target g = 4 3 + 1 3k . Let λ be the current real number input (guess) for the dual approximation. The key point is to show how it is possible to build a schedule of length at most 4λ 3 + λ 3k , starting from the assumption that there exists a schedule of length lower than λ.

The idea is to partition the set of tasks on the CPUs into two sets, each consisting in two shelves: a first set with a shelf of length λ and the other of length λ 3 , and a second set with two shelves of length 2λ 3 . The partition ensures that the makespan on the CPUs is lower than 4λ 3 . The same partition can be applied to the set of tasks on the GPUs, with the smallest shelf of length λ 3 + λ 3k instead of λ 3 . Since the tasks are independent, the scheduling strategy is straightforward when the assignment of the tasks has been determined and yields directly a solution of length at most 4λ 3 . The main problem is to assign the tasks in each shelf on the CPUs or on the GPUs in order to obtain a feasible solution. This will be done using dynamic programming. The main steps are summarized in the following algorithmic scheme:

(1) Compute the guess λ = Bmin+Bmax where B min (resp. B max ) is a lower (resp. upper) bound of the optimal makespan. (2) Search for an allotment of the tasks such that:

• the total load (work) on CPUs is at most mλ,

• the total load (work) on GPUs is at most kλ,

• the tasks alloted to the CPUs (resp. GPUs) whose processing time is strictly greater than 2λ 3 occupy a maximum number of CPUs (resp. GPUs) denoted by µ (resp. κ).

• the tasks alloted to the CPUs (resp. GPUs) whose processing time is strictly greater than λ 3 and lower than 2λ 3 can be assigned two by two to a maximum number of CPUs (resp. GPUs) denoted by µ ′ /2 (resp.κ ′ /2). The total number of CPUs (resp. GPUs) must not exceed

m (resp. k), i.e. µ + µ ′ /2 ≤ m (resp. κ + κ ′ /2 ≤ k),
• the tasks assigned to the CPUs (resp. GPUs) with processing time lower that λ 3 can be scheduled such that the induced makespan on the CPUs (resp. GPUs) will be at most equal to 4λ 3 (resp. 4λ 3 + λ 3k ). (3) If such an allotment does not exist, adjust the bound B min to λ and restart the process (Step 1). (4) If such an allotment exists, build the corresponding schedule with sets of shelves such that the makespan is lower than 4 3 λ, adjust the bound B max to λ and restart the process.

3.2.

Structure of an Optimal Schedule of length at most λ. We introduce an allocation function π(j) of a task T j which corresponds to the processor where the task is processed. The set C (resp. G) is the set of all the CPUs (resp. GPUs). Therefore, if a task T j is assigned to a CPU, we can write π(j) ∈ C. We define W C as the computational area of the CPUs on the Gantt chart representation of a schedule, i.e. the sum of all the processing times of the tasks allocated to the CPUs: W C = j / π(j)∈C p j . This corresponds to the computational load of all the CPUs.

To take advantage of the dual approximation paradigm, we have to make explicit the consequences of the assumption that there exists a schedule of length at most λ. We state below some straightforward properties of such a schedule. They should give the insight for the construction of the solution.

Property 1. In an optimal solution, the execution time of each task is at most λ, and the computational area on the CPUs (resp. GPUs) is at most mλ (resp. kλ).

Property 2. In an optimal solution, if there exist two tasks executed on the same processor such that one of these tasks has an execution time greater than 2λ 3 , then the other one has an execution time lower than λ 3 . Property 3. Two tasks with processing times on CPU (resp. GPU) greater than λ 3 and lower than 2λ 3 can be executed successively on the same CPU (resp. GPU) within a time at most 4λ 3 . The basic idea of the solution that we propose comes from the analysis of the shape of an optimal schedule. From Property 2, the tasks whose execution times on CPU (resp. on GPU) are strictly greater than 2λ 3 do not use more than m CPUs (resp. k GPUs), and hence can be executed concurrently in the first set in a shelf denoted by S 1 (resp. S 5 ). We denote by µ (resp. κ) the number of CPUs (resp. GPUs) executing these tasks. These shelves are represented for the CPUs in Figure 1. The tasks whose execution times are lower than 2λ 3 and strictly greater than λ occupied by S 1 (resp. S 5 ) from Property 1. Moreover, from Property 3, 2(m -µ) (resp. 2(k -κ)) of these tasks on CPU (resp. GPU) can be executed in time at most 4λ 3 on the remaining (m -µ) CPUs (resp. (k -κ) GPUs) in the second set and fill two shelves S 3 and S 4 (resp. S 7 and S 8 ) of equal length 2λ 3 . The tasks remaining to be assigned to the CPUs (resp. GPUs) have a processing time lower than λ 3 . The µ (resp. κ) longest remaining tasks are assigned to the first set on the CPUs (resp. GPUs) in another shelf denoted by S 2 (resp. S 6 ). The length of S 2 is λ 3 and S 6 has a length of λ 3 + λ 3k . W L (resp. W R ) will denote the computational area on the CPUs (resp. GPUs) remaining idle after this allocation in the schedule of length 4λ 3 + λ 3k . W L corresponds to the stripped area in Figure 1. Regarding the question of how the remaining tasks fit in the constructed schedule, we state the following lemma: Lemma 1. The tasks remaining to be assigned on the CPUs after the construction of S 1 , S 2 , S 3 , S 4 fit in the remaining free computational space W L between these shelves.

Proof. The tasks remaining to be assigned after the construction of S 1 , . . . , S 4 all have a processing time lower than λ 3 by construction and they necessarily fit into the remaining computational space W L , otherwise the schedule would not satisfy Property 1. The following algorithm can be used to schedule these tasks:

• Consider the remaining tasks ordered by decreasing order of processing time on CPU T 1 , . . . , T f , f being the total number of tasks remaining to be allocated. • At each step i, i = 1, . . . , f , allocate task T i to the least loaded processor, at the latest possible date. Update its load. At each step, the least loaded processor has a load at most λ; otherwise it would contradict the fact that the total work area of the tasks is bounded by mλ (according to Property 1). Hence, the idle time interval on the least loaded CPU has a length at least equal to λ 3 and can contain the task T i , which proves the correctness of the scheduling algorithm.

The question of how the remaining tasks to be assigned on GPUs fit in the constructed schedule will be addressed later in the paper (see Lemma 2) .

3.3.

Partitioning the Tasks into Shelves. In this section, we detail how to fill the shelves (see Figure 1) on the CPUs (resp. GPUs) by specifying an initial assignment of the tasks to the processors.

In order to obtain a 2-sets and 4-shelves schedule on the CPUs (resp. GPUs), we look for an assignment satisfying the following constraints:

• (C 1 ) The total computational area W C on the CPUs is at most mλ.

• (C 2 ) The set T 1 of tasks on the CPUs with an execution time strictly greater than 2λ 3 in the allotment, to be scheduled in S 1 , uses a total of at most m processors. We still denote by µ the number of processors they use.

• (C 3 ) The set T 2 of tasks on the CPUs with an execution time lower than 2λ 3 and strictly greater than λ 3 in the allotment, to be scheduled in S 3 or S 4 , uses a total of at most 2(m -µ) processors.

• (C 4 ) The total computational area on the GPUs is at most kλ.

• (C 5 ) The set T 3 of tasks on the GPUs with an execution time strictly greater than 2λ 3 in the allotment, to be scheduled in S 5 , uses a total of at most k processors. We still denote by κ the number of processors they use.

• (C 6 ) The set T 4 of tasks on the GPUs with an execution time lower than 2λ 3 and strictly greater than λ 3 in the allotment, to be scheduled in S 7 or S 8 , uses a total of at most 2(k -κ) processors.

Let us notice that if constraints (C 3 ) and (C 6 ) are satisfied, then constraints (C 2 ) and (C 5 ) will also be satisfied. Hence, constraints (C 2 ) and (C 5 ) are relaxed.

We define for each task T j a binary variable x j such that x j = 1 if T j is assigned to a CPU or 0 if T j is assigned to a GPU. Determining if an allotment satisfying (C 1 ), (C 3 ), (C 4 ) and (C 6 ) exists reduces to solving a three-dimensional knapsack minimization problem that can be formulated as follows:

W * C = min n j=1 p j x j (1) s.t. 1 2
2λ/3 pj >λ/3

x j + pj >2λ/3

x j m (2) 1 2 2λ/3 pj >λ/3 (1 -x j ) + pj >2λ/3 (1 -x j ) k (3) n j=1 p j (1 -x j ) kλ (4) x j ∈ {0, 1} (5)
Equation ( 1) represents the minimal workload on all the CPUs. Constraint (2) (resp. Constraint ( 3)) imposes that no more than m (resp. k) tasks can be executed on the CPUs (resp. GPUs) with a processing time greater than 2λ 3 , we note µ = pj >2λ/3

x j (resp. κ = pj >2λ/3

(1 -x j )) their number; and that there cannot be more than 2(m -µ) (resp. 2(k -κ)) tasks on the CPUs (resp. GPUs) with a processing time lower than 2λ 3 and greater than λ 3 (cf. Constraint (C 3 ) (resp. (C 6 )))). Constraint (4) imposes an upper bound on the computational area on the GPUs which is kλ (cf. (C 4 )).

We propose a dynamic programming algorithm in O n 2 m 2 k 3 to solve the knapsack problem. For this purpose, we first discretize the processing times of the tasks on the GPUs. We introduce ν j = pj λ/(3n) to represent the number of integer time intervals of length λ 3n required for a task T j if it is executed on the GPUs, as shown in Figure 2. N = π(j)∈G ν j denotes the total integer number of these intervals on the GPUs. We thus define the error on the processing time of each task ǫ j = p j -ν j λ 3n induced by this time discretization.

This result allows us to consider only N states in the dynamic programming regarding the workload on the GPUs. The error ǫ j on each task is at most λ 3n , so if all the tasks were assigned to one of the GPU, we would have underestimated the processing time on this GPU by at most n λ 3n = λ 3 . Constraint (4) becomes:

(6) N = π(j)∈G ν j 3kn
The approximated computational area of the GPUs is at most kλ and thus, the full computational area on GPU remains lower than kλ + λ 3 . This allows us to answer the question of how the remaining tasks to be assigned on GPUs fit in the constructed schedule: Lemma 2. The tasks remaining to be assigned on the GPUs after the construction of S 5 , S 6 , S 7 , S 8 fit in the remaining free computational space W R between these shelves. 00 00 00 00 00 00 00 00 00 00 00 00 Proof. The proof is similar to the one of Lemma 1. If we modify the starting time of the tasks of S 6 , currently λ, so that all the working processors complete their tasks at 4λ 3 + λ 3k , creating an idle time interval between the end of S 5 and the starting time of S 6 , the load of a GPU is equal to 4λ 3 + λ 3k minus the length of the idle time interval.

With the same algorithm as for Lemma 1, the only problem that may occur is if a task T i remaining to be allocated cannot be completed before the starting time of the tasks of S 6 . But at each step, the least loaded processor has a load at most λ + λ 3k since the total work area of the tasks is bounded by k λ + λ 3k . Hence, the idle time interval on the least loaded GPU has a length at least λ 3 and can contain the task T i .

We define W C (j, µ, µ ′ , κ, κ ′ , N ) as the minimum sum of all the processing times of the tasks on the CPUs when the first j tasks are considered, with among the tasks on the CPUs (resp. GPUs), µ (resp. κ) of them having processing times greater than 2λ 3 and µ ′ (resp. κ ′ ) with λ 3 < p j (resp. p j ) 2λ 3 , and where N time intervals are occupied on the GPUs.

We use a dynamic programming algorithm to compute the value of W C (j, µ, µ ′ , κ, κ ′ , N ) with the values of W C with j -1 tasks considered that were previously computed. The optimal value of the computational area W C on the CPUs will be given by

W * C = min W C (n, µ, µ ′ , κ, κ ′ , N ) 0 µ m, 0 µ ′ 2(m-µ), 0 κ k, 0 κ ′ 2(k-κ), 0 N 3kn
If W * C is greater than mλ, then there exists no solution with a makespan at most λ, and the algorithm answers "NO" in the dual approximation framework. Otherwise, the guess λ is large enough, we construct a feasible solution with a makespan at most 4λ 3 + λ 3k , with the shelves and the corresponding µ, µ ′ , κ, κ ′ and N values. The dynamic programming algorithm represents one step of the dual-approximation algorithm, with a fixed guess λ. A binary search is then used to try different guesses to approach the optimal makespan as explained in Section 3.1.

, we conclude that the value of λ is too small and adjust the bounds of our binary search accordingly.

If W C mλ + λ 2 , we can construct a feasible schedule with a makespan lower than 3 2 λ with the previous algorithm. Indeed, the number of tasks in C 1 is lower than m so we can build a shelf S 1 as we did in Section 3, occupying |C 1 | CPUs, with a length at most λ. The same arguments given in the proof of Lemma 1 can be used for building a shelf S 2 of length λ 2 and all the tasks from C 2 can be fitted in the schedule as before. For the GPUs, the algorithm makes sure that the number of tasks in G 1 is lower than k and that W G does not go over the bound of kλ + λ 2 so shelves similar to S 1 and S 2 can be built easily. However, we did not have to make any discretization on the processing times of the tasks assigned to the GPUs here, so, contrary to Section 3, we get the same performance ratio of 3 2 for any number k of GPUs. The time complexity of an algorithm based on this principle is in O (mn log n).

Experiments

We run some experiments in order to show the efficiency of the proposed algorithmic scheme. They are two-fold, we first run a relaxed version of our algorithm by simulations on random instances and compared them to the classical reference algorithm HEFT [START_REF] Topcuoglu | Performance-effective and low-complexity task scheduling for heterogeneous computing[END_REF] used on several actual systems (HEFT stands for Heterogeneous Earliest Finishing Time). Then, we implement the algorithm on the top of the xKaapi runtime system. 5.1. Preliminary: Analysis of HEFT. HEFT proceeds in two phases, starting by a prioritization of the tasks that are sorted by decreasing average execution time and then the processor selection is obtained with the heterogeneous earliest finish time rule. Lemma 3. The worst case performance ratio of HEFT is larger than m/2.

Proof. We show on the following instance that the prioritizing phase can provide a schedule whose makespan is far from the optimum. Let us consider an instance with a list of the following tasks: m tasks such that p = 1 and p = ǫ and for i = 0, • • • , m -1: a single task of type A such that p = 1 -i/m and p = 1 -i/m and m -1 tasks of type B such that p = 1 -i/m and p = 1/m 2 . All these tasks are executed faster on the GPUs.

For this instance, HEFT fills first the m CPUs. Then, it fills alternatively the GPU with one task of type A and the m CPUs with m tasks of size B. HEFT ends with a makespan equal to m/2 + 3/2 -1/m. It is easy to check that the optimal makespan is equal to 1.

5.2.

Simulations. The dual approximation based algorithm for problem (P m, P k) || C max presented in Section 3 provides a performance ratio of 4 3 + 1 3k with a reasonable time complexity. However, this running cost is not comparable to the one of HEFT which basically only needs to sort the tasks (n log n). This method can be modified in order to obtain a slightly worse performance ratio of 2 in a much slower time O(n 2 k), which would make it comparable to HEFT in terms of running time and still provide a performance guarantee. The idea is based on leaving aside the constraints ordering the tasks into shelves. The only constraint that remains is the one on the computational area on the GPUs being lower than kλ, λ being the current guess of the dual approximation. With the optimal computational area on the CPUs under this constraint determined by dynamic programming, we build a schedule with a makespan lower than twice the optimal value. This 2-approximation algorithm, denoted in short by DP in what follows, was implemented and compared to HEFT by simulations based on various classes of instances. Moreover, for the special case where all the tasks are accelerated, we implemented the algorithm presented in Section 4, denoted by Accel, which provides a performance ratio of 3 2 with a time complexity of O (mn log n). All these algorithms are implemented in C++ programming language and run on a 3.4 GHz PC with 15.7 Gb RAM.

We report below a series of experiments run on random instances of various sizes: from 10 to 1000 tasks, with a step of 10 tasks, 2 a CPUs, a varying from 0 to 6, and 2 b GPUs, b varying from 0 to 3. For each combination of these sizes, 30 instances were considered, bringing us to a total of 10500 tested instances. The processing times on the CPUs are randomly generated using the uniform distribution U [START_REF] Augonnet | Starpu: A unified platform for task scheduling on heterogeneous multicore architectures[END_REF]100]. The distribution of the acceleration factors on the GPUs has been measured in [START_REF] Seifu | Scheduling on heterogeneous cluster environments[END_REF] using the classical numerical kernels of Magma [START_REF] Agullo | Numerical linear algebra on emerging architectures: The plasma and magma projects[END_REF] in a multi-core multi-GPU machine hosted by the Grid'5000 infrastructure experimental platform [START_REF] Bolze | Grid'5000: A large scale and highly reconfigurable experimental grid testbed[END_REF]. We extracted a distribution of the acceleration factors which reflects the qualitative speed-up on classical numerical kernels: we assign to each task an acceleration factor α j = pj pj of 1/15 or 1/35 with a probability of 1/2. The resulting processing times on the GPUs are thus p j = α j p j . Since in this generation scheme all the tasks of these instances are accelerated on GPU, DP, HEFT and Accel were all compared on these instances. The running time of the three algorithms is always under one second, even for the largest instances. We calculated the mean and maximal deviations of the makespans of the solutions returned by these algorithms from the lower bound of the makespan derived from the binary search of the approximation algorithm, over all the instances. As can seen in Table 1, the maximal deviations of DP are usually below the maximal deviations of HEFT and more importantly, these deviations respect the theoretical performance guarantee in the case of DP whereas the maximal deviations of HEFT sometimes go over the 100% barrier corresponding to a performance ratio of 2. The same can be said for Accel, with maximal deviations staying below the 50% barrier corresponding to a performance ratio of 3 2 . Figure 3 shows that in average, DP even outperforms HEFT for large instances. However, Accel remains slightly above HEFT, staying close to its 3 2 bound but never going above it, contrary to HEFT which does not have a performance guarantee, as seen in Table 1. This better performance ratio of 3 2 with its low cost make Accel preferable to DP for practical intensive use with still a better maximum performance guarantee.

5.3.

Experiments on a real run-time. We investigate in this section the practical use of the algorithms of the previous family. We target classical linear algebra kernels, since they are extensively used and they generate a loop of independent tasks. Here, the computations consist in a series of independent tasks. DP (with ratio 2), HEFT and 4 3 -approximation algorithms were implemented in the scheduler of the xKaapi run-time system [START_REF] Gautier | Xkaapi: A runtime system for dataflow task programming on heterogeneous architectures[END_REF].

Implementation of the 4 3 -approximation algorithm. In most linear algebra applications, the block size is the most important characteristic to maximize the performance. Indeed, the block size has a direct impact on memory transfers between the host and the accelerators, cache effects and task graph size .

In the following experiments of Figure 4, we study the variation of the computation time as a function of the block size for the same matrix size of a Cholesky factorisation extracted from the MAGMA library. We use a single GPU and the matrix is decomposed over simple square blocks. The experiments have been conducted on a quad-core Intel i7-3840QM with hyperthreading and a Nvidia Quadro K1000M GPU. As the block size decreases, the number of independent tasks increases. Thus, the computation time of the scheduling using the 4 3 -approximation algorithm increases quadratically with the number of tasks too. As a result, the scheduling time dominates the execution time saved for large block size (which corresponds to small number of tasks). The 4 3 -approximation algorithm is therefore usable mostly for cases where the computation time is larger than the scheduling time. It is probably not the best suited algorithm for linear algebra kernels.

Practical issues: DP versus HEFT. We compare now the relaxed version of the algorithm (DP) with HEFT on a machine with several GPUs. The experiments have been conducted on a heterogeneous, multi-GPU system composed of two six-core Intel Xeon X5650 CPUs running at 2.66 GHz with 72 GB of memory. This parallel system is enhanced with eight NVIDIA Tesla C2050 GPUs (Fermi architecture) of 448 GPU cores (scalar processors) running at 1.15 GHz each (2688 GPU cores total) with 3 GB GDDR5 per GPU (18 GB The implementation of DP was combined with an improved local mapping for minimizing the data transfers [START_REF] Bleuse | Scheduling data flow program in xkaapi: A new affinity-based algorithm for heterogeneous architectures[END_REF]. With 8 GPUs, DP outperforms HEFT both in the raw performance and memory transfers (see Table 2). The execution times are close to each other in all cases, but DP has the three major following advantages:

• the guarantee on the scheduling results,

• a decrease in the volume of communication,

• no need of accurate communication model.

Concluding Remarks

In this paper, we presented a new scheduling algorithm using a generic methodology (in the opposite of specific ad hoc algorithms) for hybrid architectures (multicore machine with GPUs). We proposed fast algorithms with a constant approximation ratio in the case of independent tasks and an faster version in the special case where all the tasks are accelerated when assigned to GPU. A ratio of 4 3 + 1 3k + ǫ is achieved for k GPUs in the general case, and a ratio of 3 2 + ǫ when the tasks are all accelerated on GPU. The main idea of the approach is to determine an adequate partition of the set of tasks on the CPUs and the GPUs using a dual approximation scheme. A simulation and experimental analysis on a real run-time systemm have been provided to assess the computational efficiency of the proposed methods. The main conclusion is that these algorithms are stable because of their approximation guaranties, however, the high running time is often dominated by the cost of the scheduling it-self, leading to inefficiency if the number of tasks is too small. According to our experimental setting, the relaxed version with approximation ratio equal to 2 was the best trade-off.

As further investigations of this work, we plan to extend the analysis to more generic problems where the tasks are dependent with larger number of tasks. We also plan on testing the robustness of the presented algorithms to perturbations in the execution times that are only estimated on real-life computing platforms.

2

 2

2 Figure 1 .

 21 Figure 1. Partitioning the set of tasks on the CPUs into two sets of two shelves, the first one occupying µ CPUs, the second m -µ CPUs.

1 Figure 2 .

 12 Figure 2. Rounded allocation of two tasks T 1 with p 1 = 6.5 and T 2 with p 2 = 4.7 on a GPU

Figure 3 .

 3 Figure 3. Mean deviations of DP and HEFT for various n

Figure 4 .

 4 Figure 4. Execution time of a Cholesky factorization scheduled by DP, 4/3 and HEFT for various block sizes, on 3 hyperthreaded core and a single GPU

Table 1 .

 1 Maximal deviations (%) for DP and HEFT

	n	120	160	220	260	360	380	660	700	760	780	920	940
	DP	76.88	72.73	70.37	69.14	70.00	70.00	67.42	50.82	42.77	54.47	91.77	63.07
	HEFT	123.53	98.44	92.55	91.90	110.37	91.78	113.48	98.10	98.77	103.15	116.46	96.31
	Accel	46.15	42.86	50.00	41.18	37.82	43.59	36.36	40.91	32.52	37.04	48.24	34.65

Table 2 .

 2 in total). It has 4 PCIe switches to support up to 8 GPUs. When 2 GPUs share a switch, their aggregated PCIe bandwidth is bounded by the one of a single PCIe 16x. Perfomance of DP and HEFT for Cholesky factorization with m=4 CPUs and k=8 GPUs

	Algorithm Gflops Memory transfer/GB
	HEFT	535	2.62
	DP	565	1.91

Univ. Grenoble-Alpes, LIG, F-38000 Grenoble, France., CNRS, LIG, F-38000 Grenoble, France, Inria

[START_REF] Hochbaum | A polynomial approximation scheme for scheduling on uniform processors: using the dual approximation approach[END_REF] Sorbonne Universités, UPMC Univ.Paris 06, UMR 7606, LIP6, F-75005, Paris, France

R. BLEUSE, S. KEDAD-SIDHOUM, F. MONNA, G. MOUNI É, AND D. TRYSTRAM

This work was partially supported by the French Ministry of Defense (DGA)..

Cost Analysis. Solving the dynamic programming algorithm for a fixed value of λ requires to consider O n 2 m 2 k 3 states, since 1 j n, 1 µ m, 1 µ ′ 2(m -µ), 1 κ k, 1 κ ′ 2(k -κ), and 0 N 3kn. Therefore, the time complexity of each step of the binary search is O n 2 m 2 k 3 .

Analysis of a Special Case

We consider, in this section, a version of the problem (P m, P k) || C max where all the tasks are accelerated when assigned to GPU, since this is the case in most applications, i.e p j p j for j = 1, . . . , n (all the tasks do not have the same acceleration factor on GPU). When considering this special case in the algorithm presented in the previous section, we note no amelioration in the time complexity of the algorithm, whereas the problem was simpler. Therefore we present an new algorithm for this case, based on a similar scheme, but with a much lower time complexity for an approximation ratio of 3 2 . The algorithm is also based on the dual approximation technique. At each step, we have a guess on the optimal makespan. Let us consider one step of the dual approximation scheme and let as before λ denote the current guess.

The idea is to divide the set of tasks T into four sets of tasks, two of them whose tasks will be assigned to a CPU, C 1 , C 2 , and the other two whose tasks will be assigned to a GPU, G 1 and G 2 . We denote the cardinality of set C i (resp. G i ) by

The first step of the algorithm consists in a preliminary assignment of the tasks of the whole set T to two of the sets: G 1 and C 2 . This assignment of each task T j is done by considering the value of the processing time p j of T j on CPU. The possible cases are described below:

• Otherwise, λ 2 < p j , task T j is assigned to G 1 . This assignment does not guarantee that the resulting schedule will have a makespan lower than 3 2 λ, even if there exists a schedule of makespan lower than λ. However, if |G 1 | > k + m, there are too many tasks with a processing time on CPU greater than λ 2 to fit in a schedule of makespan lower than λ. The dual approximation rejects this guess λ.

In the second step of the algorithm, in order to achieve the desired performance ratio, we have to reassign some of the tasks allocated to C 2 to G 2 and some of the tasks allocated to G 1 to C 1 . We order the tasks of sets C 2 and G 1 in decreasing order of the computational surface change induced when a task T j changes from a CPU to a GPU, p j -p j . One exception has to be made for set G 1 : some tasks can have a processing time on CPU larger than λ. These tasks are too big to fit on the CPUs with the current guess. They cannot be reassigned and are put at the end of G 2 , no matter the impact they can have on the computational surfaces. The tasks of the two sets will be examined in this new order. We can note that at most m tasks of set G 1 can be reassigned to C 1 . The following steps are therefore repeated at most m + 1 times, as long as we have

With this assignment, the computational area on the CPUs has been reduced to a minimum with the constraint of keeping the computational area on the GPUs lower than kλ + λ 2 . Therefore, the value of W C obtained by our algorithm is smaller than the value of the computational area on the CPUs of the optimal schedule, the most accelerated tasks having been assigned to the GPUs. Therefore, if W C > mλ + λ


Raphaël Bleuse
email: raphael.bleuse@imag.fr

Thierry Gautier
email: thierry.gautier@inrialpes.fr

João V F Lima
email: jvlima@inf.ufsm.br

Grégory Mounié
email: gregory.mounie@imag.fr

Denis Trystram
email: denis.trystram@imag.fr

Scheduling data flow program in XKaapi: A new affinity based Algorithm for Heterogeneous Architectures

Keywords: heterogeneous architectures, scheduling, cost models, dual approximation scheme, programming tools, affinity

Efficient implementations of parallel applications on heterogeneous hybrid architectures require a careful balance between computations and communications with accelerator devices. Even if most of the communication time can be overlapped by computations, it is essential to reduce the total volume of communicated data. The literature therefore abounds with ad hoc methods to reach that balance, but these are architecture and application dependent. We propose here a generic mechanism to automatically optimize the scheduling between CPUs and GPUs, and compare two strategies within this mechanism: the classical Heterogeneous Earliest Finish Time (HEFT) algorithm and our new, parametrized, Distributed Affinity Dual Approximation algorithm (DADA), which consists in grouping the tasks by affinity before running a fast dual approximation. We ran experiments on a heterogeneous parallel machine with twelve CPU cores and eight NVIDIA Fermi GPUs. Three standard dense linear algebra kernels from the PLASMA library have been ported on top of the XKaapi runtime system. We report their performances. It results that HEFT and DADA perform well for various experimental conditions, but that DADA performs better for larger systems and number of GPUs, and, in most cases, generates much lower data transfers than HEFT to achieve the same performance.

Introduction

With the recent evolution of processor design, the future generations of processors will contain hundreds of cores. To increase the performance per watt ratio, the cores will be non-symmetric with few highly powerful cores (CPU) and numerous, but simpler, cores (GPU). The success of such machines will rely on the ability to schedule the workload at runtime, even for small problem instances.

One of the main challenges is to define a scheduling strategy that may be able to exploit all potential parallelisms on a heterogeneous architecture composed of multiple CPUs and multiple GPUs. Previous works demonstrate the efficiency of strategies such as static distribution [START_REF] Song | A scalable framework for heterogeneous GPU-based clusters[END_REF][START_REF] Tomov | Towards dense linear algebra for hybrid GPU accelerated manycore systems[END_REF], centralized list scheduling with data locality [START_REF] Bueno | Productive Programming of GPU Clusters with OmpSs[END_REF], cost models [START_REF] Agullo | Lu factorization for accelerator-based systems[END_REF][START_REF] Agullo | QR Factorization on a Multicore Node Enhanced with Multiple GPU Accelerators[END_REF][START_REF] Augonnet | Automatic calibration of performance models on heterogeneous multicore architectures[END_REF][START_REF] Augonnet | StarPU: a unified platform for task scheduling on heterogeneous multicore architectures[END_REF] based on Heterogeneous-Earliest-Finish-Time scheduling (HEFT) [START_REF] Topcuoglu | Performance-effective and low-complexity task scheduling for heterogeneous computing[END_REF], and dynamic for some specific application domains [START_REF] Bosilca | DAGuE: A generic distributed DAG engine for High Performance Computing[END_REF][START_REF] Hermann | Multi-GPU and Multi-CPU Parallelization for Interactive Physics Simulations[END_REF]. Locality-aware work stealing [START_REF] Gautier | XKaapi: A Runtime System for Data-Flow Task Programming on Heterogeneous Architectures[END_REF], with a careful implementation to overlap communication by computation [START_REF] Lima | Exploiting Concurrent GPU Operations for Efficient Work Stealing on Multi-GPUs[END_REF], improves significantly the performance of compute-bound linear algebra problems such as matrix product and Cholesky factorization.

Nevertheless, none of the above cited works considers scheduling strategies from the viewpoint of a compromise between performance and locality. In this paper, we propose a scheduling algorithm based on dual approximation [START_REF] Kedad-Sidhoum | Scheduling independent tasks on multi-cores with gpu accelerators[END_REF] that uses a performance model to predict the execution time of tasks during scheduling decision. This algorithm, called Distributed Affinity Dual Approximation (DADA), is able to find a compromise between transfers and performance. It is parametrized by α for tuning this trade-off. The main advantage of dual approximation algorithms is their theoretical performance guarantee as they have a constant approximation ratio. On the contrary, the worst case of HEFT can be arbitrarily bad [START_REF] Kedad-Sidhoum | Scheduling independent tasks on multi-cores with gpu accelerators[END_REF].

We compare these two different scheduling strategies for data-flow task programming. These strategies are implemented on top of the XKaapi scheduling framework with performance models for task execution time and transfer prediction. The contributions of this paper are first the design and implementation of dual approximation scheduling algorithms (with and without affinity) and second their evaluation in comparison to the well-known HEFT algorithm on three dense linear algebra algorithms in double precision floating-point operations from PLASMA [START_REF] Buttari | A class of parallel tiled linear algebra algorithms for multicore architectures[END_REF]: namely Cholesky, LU, and QR. To our knowledge, this paper is the first report of experimental evaluations studying the impact of data transfer model and contention on a machine with up to 8 GPUs.

The main lesson of this work is that scheduling algorithms need extra information in order to take the right decisions. Such extra information could be obtained in a precise communication model to predict processing time of each task or in a more flexible information such as the affinity in DADA. Even if HEFT remains a good candidate for scheduling such linear algebra kernels, DADA is highly competitive against it for multi-GPU systems: the experimental results demonstrate that it achieves the same range of performances while reducing significantly the communication volume.

The remainder of this paper is organized as follows. Section 2 provides an overview of XKaapi runtime system, describes the XKaapi scheduling framework and the cost model applied for performance prediction. Section 3 details the two studied scheduling strategies. Section 4 presents our experimental results on a heterogeneous architecture composed of 12 CPUs and 8 GPUs. In Section 5 we briefly survey related works on runtime systems, scheduling strategies and performance prediction. Finally, Section 6 concludes the paper and suggests future directions.

Scheduling framework in XKaapi

The XKaapi5 data-flow model [START_REF] Gautier | KAAPI: A thread scheduling runtime system for data flow computations on cluster of multi-processors[END_REF] -as in Cilk, Intel TBB, OpenMP-3.0, or OmpSs [START_REF] Bueno | Productive Programming of GPU Clusters with OmpSs[END_REF] -enables non-blocking task creation: the caller creates the task and proceeds with the program execution. Parallelism is explicit while the detection of synchronizations is implicit [START_REF] Gautier | KAAPI: A thread scheduling runtime system for data flow computations on cluster of multi-processors[END_REF]: dependencies between tasks and memory transfers are automatically managed by the runtime system.

XKaapi runtime system is structured around the notion of worker : it is the internal representation of kernel threads. It executes the code of the tasks and takes local scheduling decisions. Each worker owns a local queue of ready tasks. Our interface is mainly inspired by work stealing scheduler and composed of three operations that act on workers' queues of tasks: pop, push and steal. In our previous work, we demonstrated the efficiency of XKaapi locality-aware work stealing as well as the corresponding multi-GPU runtime support [START_REF] Gautier | XKaapi: A Runtime System for Data-Flow Task Programming on Heterogeneous Architectures[END_REF] using specialized implementation of these operations. A new operation, called activate, has been defined to push ready task to a worker's queue.

Execution flow

The sketch of the execution mechanism is the following: at each step, either the own queue of worker is not empty and the worker uses it; or the worker emits a steal request to a randomly selected worker in order to get a task to execute. According to the dependencies between tasks, once a worker performs a task, it calls the activate operations in order to activate the successors of the task which become ready for execution.

The XKaapi runtime system gets information from each internal events (such as start-end of task execution, or start-end of communication toward GPU) to calibrate the performance model and corrects erroneous predictions due to unpredictable or unknown behavior (e.g. operating system state or I/O disturbance). StarPU [START_REF] Augonnet | StarPU: a unified platform for task scheduling on heterogeneous multicore architectures[END_REF] uses similar runtime measurements in order to correct the performance predictions in their HEFT implementation.

All of our scheduling strategies follow this sketch. Every worker terminates its execution when all the tasks of the application have been executed.

Pop, Push, Steal and Activate Operations

A framework interface for scheduling strategies is not a new concept in heterogeneous systems. Bueno et al. [START_REF] Bueno | Productive Programming of GPU Clusters with OmpSs[END_REF] and Augonnet et al. [START_REF] Augonnet | StarPU: a unified platform for task scheduling on heterogeneous multicore architectures[END_REF] described a minimal interface to design scheduling strategies with selection at runtime. However, there is little information available on the comparison of different strategies. Most of them reported performance on centralized list scheduling and performance models. Our framework is composed of three classical operations in the work stealing context, plus an action to activate tasks when their predecessors have completed.

-The push operation inserts a task into a queue. A worker can push a task into any other workers' queue. -A pop removes a task from the local queue owned by the caller worker.

-A steal removes a task from the queue of a remote worker. It is called by an idle thread -the thief -in order to pick tasks from a randomly selected worker -the victim. -The activate operation is called after the completion of a task. The role of this operation is to allocate the tasks that are ready to be executed. Hence, most of the scheduling decision are done during this operation.

Performance Model

Cost models depend on a certain knowledge of both application algorithm and the underlying architecture to predict performance at runtime. In order to predict performance, we designed a StarPU [START_REF] Augonnet | Automatic calibration of performance models on heterogeneous multicore architectures[END_REF] like performance model for task execution time and communication. Our task prediction relies on an historybased model, and transfer time estimation is based on asymptotic bandwidth. They are associated with scheduling strategies that are based on task completion time such as HEFT and DADA with and without affinity.

In order to balance efficiently the load, for each processor XKaapi maintains a shared time-stamp of the predicted time when it has completed its tasks. The completion date of the last executed task is also kept. The update and incrementation of the time-stamps are efficiently implemented with atomic operators.

Scheduling Strategies

This section introduces the scheduling strategies designed on top of the XKaapi scheduling framework. We consider a multi-core parallel architecture with m homogeneous CPUs and k homogeneous GPUs. First, we describe our implementation of HEFT [START_REF] Topcuoglu | Performance-effective and low-complexity task scheduling for heterogeneous computing[END_REF]. Then, we recall the principle of the dual approximation scheme [START_REF] Hochbaum | Using dual approximation algorithms for scheduling problems theoretical and practical results[END_REF]. We propose a new algorithm -Distributed Affinity Dual Approximation (DADA) -based on this paradigm which takes into account the affinity between tasks.

In the following we denote by p CP U i the processing time of task T i on a CPU and p GP U i on a GPU. We define the speedup S i of task T i as the ratio

S i = p CP U i / p GP U i .

HEFT within XKaapi

The Heterogeneous Earliest-Finish-Time algorithm (HEFT), proposed by [START_REF] Topcuoglu | Performance-effective and low-complexity task scheduling for heterogeneous computing[END_REF], is a scheduling algorithm for a bounded number of heterogeneous processors. Its Algorithm 1: HEFTactivate operation.

Input : A list of ready tasks Ti LR Output: Tasks Ti pushed to the worker's queues time complexity is in O(n 2 • (m + k)). It has two major phases: task prioritizing and a worker selection. Our XKaapi version of HEFT implements both phases during the activate operation. The task prioritizing phase computes for all ready tasks T i its speedup S i relative to an execution on GPU. Next, it sorts the list of ready tasks by decreasing speedups. Whereas the original HEFT rule sorts the tasks by decreasing upward rank (average path length to the end), our rule gives priority on minimizing the sum of the execution times. In the worker selection phase, the algorithm selects tasks in the order of their speedup S i and schedules each task on the worker which minimizes the completion time. Algorithm 1 describes the basic steps of HEFT over XKaapi.

1 foreach Ti ∈ LR do 2 Si ← p CP U i / p GP U

Dual Approximation and Affinity

Dual Approximation Let us recall first that a ρ-dual approximation scheduling algorithm considers a guess λ (which is an estimation of the optimal makespan) and either delivers a schedule of makespan at most ρλ or answers correctly that there exists no schedule of length at most λ [START_REF] Hochbaum | Using dual approximation algorithms for scheduling problems theoretical and practical results[END_REF]. The process is repeated by a classical binary search on λ up to a precision of ǫ. We target ρ = 2. The dual approximation part of Algorithm 2 consists in the following steps:

-Choice of the initial guess λ (lines 2 and 4); -Extract the tasks which fit only into GPUs (p CP U i > λ), and symmetrically those which are dedicated to CPUs (line 9); -Keep this schedule if the tasks fit into λ (line 12). Otherwise, reject it if there is a task larger than λ on both CPUs and GPUs (line 15); -Add to the tasks allocated to the GPU those which have the largest speedup S i up to overreaching the threshold λ (line 19) which guarantees the ratio ρ = 2; -Put all the remaining tasks in the m CPUs and execute them using an earliest-finish-time scheduling policy (line 19).

Algorithm 2: DADAactivate operation.

Input : A list of ready tasks Ti LR Output: Tasks Ti pushed to the worker's queues Affinity DADA builds a compromise taking into account both raw performance and transfers. The principle consists in two successive phases: a first local phase targeting the reduction of the communications through the abstraction described below and a second phase which counter-balances the induced serialization aiming at a global balance. Any algorithm optimizing the makespan could be used for the second phase. We use a basic dual-approximation. In order to gain a finer control, the length of the first phase is controlled by a parameter (denoted by α, 0 ≤ α ≤ 1). A value of 0 for α means that the affinity is not taken into account: DADA is then a basic dual-approximation. While at the opposite a value close to 1 allows a length up to λ for the first phase, thus giving a greater weight to affinity.

1 lower ← 0 2 upper ← i max(p CP U i , p GP U i ) 3 while (upper -lower) > ǫ do
Each pair (task, computation resource) is given an affinity score. Maximizing the score over the whole schedule enables to consider local impacts. The affinity scores are computed using extra information of the runtime system. In our implementation, they were computed using the amount of data updated by each task. For instance, a task that writes or modifies a data stored on a resource R has a high score and is prone to be scheduled on R.

Experimental setup: Platform & Benchmarks

Platform All experiments have been conducted on a heterogeneous, multi-GPU system. It is composed of two hexa-core Intel Xeon X5650 CPUs running at 2.66 GHz with 72 GB of memory. It is enhanced with eight NVIDIA Tesla C2050 GPUs (Fermi architecture) of 448 GPU cores (scalar processors) running at 1.15 GHz each (2688 GPU cores total) with 3 GB GDDR5 per GPU (18 GB total). The machine has 4 PCIe switches to support up to 8 GPUs. When 2 GPUs share a switch, their aggregated PCIe bandwidth is bounded by the one of a single PCIe 16x. Experiments using up to 4 GPUs avoid this bandwidth constraint by using at most 1 GPU per PCIe switch.

Benchmarks All benchmarks ran on top of a GNU/Linux Debian 6.0.2 squeeze with kernel 2.6.32-5-amd64. We compiled with GCC 4.4 and linked against CUDA 5.0 and the library ATLAS 3.9.39 (BLAS and LAPACK). All experiments use the tile algorithms of PLASMA [START_REF] Buttari | A class of parallel tiled linear algebra algorithms for multicore architectures[END_REF] for Cholesky (DPOTRF), LU (DGETRF), and QR (DGEQRF). The QUARK API [START_REF] Yarkhan | Quark users' guide: Queueing and runtime for kernels[END_REF] has been implemented and extended in XKaapi to support task multi-specialization: the XKaapi runtime system maintains the CPU and GPU versions for each PLASMA task. At the task execution, our QUARK version runs the appropriate task implementation in accordance with the worker architecture. The GPU kernels of QR and LU are based on previous works from [START_REF] Agullo | Lu factorization for accelerator-based systems[END_REF][START_REF] Agullo | QR Factorization on a Multicore Node Enhanced with Multiple GPU Accelerators[END_REF] and adapted from PLASMA CPU algorithm and MAGMA from [START_REF] Tomov | Towards dense linear algebra for hybrid GPU accelerated manycore systems[END_REF]. Each running GPU monopolizes a CPU core to manage its worker. The remaining CPU cores are involved in the application computations.

Methodology Each experiment has been executed at least 30 times for each set of parameters and we report on all the figures (Fig. 1, 2, 3 and 4) the mean and the 95% confidence interval. The factorizations have been done in double precision floating-point operations with a PLASMA internal block (IB ) of size 128 and tiles of size 512. For each of them, we plot the highest performance obtained on various matrix sizes with the discussed scheduling strategies.

In the following, DADA(α) represents DADA parametrized by α. We denote by DADA(α)+CP the algorithm using Communication Prediction as supplementary information. HEFT strategy always computes the earliest finish time of a task taking into account the time to transfer data before executing the task.

Impact of the affinity control parameter α

This section highlights the impact of the affinity control parameter α on the compromise between performance and data transfers. The measures have been done with the Cholesky decomposition on matrices of size 8192 × 8192 and 16384 × 16384. However, we present only results for the smallest size as we observe similar behaviors for both matrix sizes. Fig. 1 shows both performance (Fig. 1(a) and 1(b)) and total memory transfers (Fig. 1(c) and 1(d)) for several values of α with respect to the number of GPUs. Both metrics are shown without (Fig. 1(a) and 1(c)) and with (Fig. 1(b) and 1(d)) communication prediction taken into account. Once affinity is considered (i.e. α = 0), the higher the value of α, the better the policy scales. Using as little information as possible (i.e. DADA(0) and no communication prediction), the policy performance does not scale with more than two GPUs due to a too huge amount of transfers.

Comparison of scheduling strategies

We present in this section the results for the three kernels with matrix size 8192 × 8192. Other tested sizes have the same behavior. The idea is to evaluate the behavior of each strategy with different work loads. Both performance and data transfers of the policies introduced above: HEFT, DADA(0), DADA(α) and DADA(α)+CP are studied. are the policies with the lowest bandwidth footprint up to 6 GPU. Yet, as the number of GPU grows, the use of communication prediction allows to reduce the communication volume with sustained high performances.

Experimental evaluation

Fig. 3 reports the behavior of the LU factorization (DGETRF). It studies both performance results (Fig. 3(a)) and total memory transfers (Fig. 3(b)). Apart from the performance of DADA+CP for six CPUs and six GPUs (with a large confidence interval), all scheduling policies sustain the same performance. Data transfers seem to have a little impact on performance. However, DADA(α)+CP generates less memory movements than other strategies. DADA(0) is the most costly policy while DADA(α) and HEFT have similar impacts.

The total memory transfers have the same shape than for the Cholesky factorization. Still, the gap between the curves is widening: DADA(α)+CP is 3.5 less demanding in bandwidth than HEFT for only a slowdown of about 1.13 in performance for 8 GPU.

Finally, Fig. 4 reports the behavior of the QR factorization (DGEQRF) with respect to the number of GPUs used. Both performance results (Fig. 4(a)-) and total memory transfers (Fig. 4(b)) are studied. All dual approximations (DADA(0), DADA(α), DADA(α)+CP) behave the same and are outperformed by HEFT. Even the low transfer footprint of both DADA(α) is not able to sustain performance. It seems that the dependencies between tasks for QR factorization have a strong impact on the schedule computed by all dual approximation algorithms. We are still investigating this particular point.

Discussion

Communication prediction Affinity is a viable alternative to communication modeling. Indeed, DADA without communication prediction is comparable to HEFT in terms of performance. Moreover, affinity based policy combined with communication prediction allows to reduce further more memory transfers (up to a factor 3.5 when compared to HEFT).

Comparison with work stealing scheduling algorithm For the sake of completeness, we also tested the work stealing algorithm. However we did not plot the results in previous figures for the sake of readability. We briefly discuss them now. The naive work stealing algorithm is cache unfriendly, especially with small matrices as its random choices are heavily penalizing [START_REF] Gautier | XKaapi: A Runtime System for Data-Flow Task Programming on Heterogeneous Architectures[END_REF]. On the contrary, the affinity policies proposed here are suitable for this case. When scheduling for medium and large matrix sizes, the impact of modeling inaccuracies grows. Model oblivious algorithms such as work-stealing behave well by efficiently overlapping communications and computations while HEFT is induced in error by the imprecise communication prediction. Hence, our approach is much more robust than work stealing and HEFT since it does not need a too precise communication model and adapts well to various matrix sizes.

Related Works

StarPU [START_REF] Augonnet | StarPU: a unified platform for task scheduling on heterogeneous multicore architectures[END_REF], OmpSs [START_REF] Bueno | Productive Programming of GPU Clusters with OmpSs[END_REF] and QUARK [START_REF] Yarkhan | Quark users' guide: Queueing and runtime for kernels[END_REF] are programming environments or libraries that enables to automatically schedule tasks with data flow dependencies. OmpSs is based on OpenMP-like pragmas while StarPU and QUARK are C libraries. QUARK does not schedule tasks on multi-GPU architecture and implements a centralized greedy list scheduling algorithm. OmpSs locality-aware scheduling, similar to our data-aware heuristic from [START_REF] Gautier | XKaapi: A Runtime System for Data-Flow Task Programming on Heterogeneous Architectures[END_REF], computes an affinity score based on where the data is located and its size. Then, the task is placed on the highest affinity resource or in a global list, otherwise. The StarPU scheduler uses the HEFT [START_REF] Topcuoglu | Performance-effective and low-complexity task scheduling for heterogeneous computing[END_REF] algorithm to schedule all ready tasks in accordance with the cost models for data transfer and task execution time [START_REF] Augonnet | Automatic calibration of performance models on heterogeneous multicore architectures[END_REF]. Our data transfer model is based on the StarPU model with minor extension. In the context of dense linear algebra algorithms, PLASMA [START_REF] Buttari | A class of parallel tiled linear algebra algorithms for multicore architectures[END_REF] provides fine-grained parallel linear algebra routines with dynamic scheduling through QUARK, which was conceived specially for numerical algorithms on multi-CPU architecture. MAGMA [START_REF] Tomov | Towards dense linear algebra for hybrid GPU accelerated manycore systems[END_REF] implements static scheduling for linear algebra algorithms on heterogeneous systems composed of GPUs. Recently it has included some methods with dynamic scheduling in multi-CPU and multi-GPU sytems on top of StarPU, in addition to the static multi-GPU version. In [START_REF] Song | A scalable framework for heterogeneous GPU-based clusters[END_REF] the authors based their Cholesky factorization on 2D block cyclic distribution with an owner compute rule to map tasks to resources. DAGuE [START_REF] Bosilca | DAGuE: A generic distributed DAG engine for High Performance Computing[END_REF] is a parallel framework focused on multi-core clusters and supports single-GPU nodes. Other papers reported performance results of task-based algorithms with HEFT cost model scheduling on heterogeneous architectures for the Cholesky [START_REF] Augonnet | StarPU: a unified platform for task scheduling on heterogeneous multicore architectures[END_REF], LU [START_REF] Agullo | Lu factorization for accelerator-based systems[END_REF], and QR [START_REF] Agullo | QR Factorization on a Multicore Node Enhanced with Multiple GPU Accelerators[END_REF] factorizations. All the results report evaluation of single floating point arithmetics with up to 3 GPUs. Due to the small number of GPUs, such studies cannot observe contention and scalability.

Conclusion

We presented in this paper a new scheduling algorithm on top of the XKaapi runtime system. It is based on a dual approximation scheme with affinity and has been compared to the classical HEFT algorithm for three tile algorithms from PLASMA on an heterogeneous architecture composed of 8 GPUs and 12 CPUs. Both algorithms attained significant speed up on the three dense linear algebra kernel. Moreover, if HEFT achieves the best absolute performance with respect to DADA on QR, while DADA has similar or better performances than HEFT on Cholesky and LU for large numbers of GPU. Nevertheless, DADA allows to significantly reduce the data transfers with respect to HEFT. More interesting, thanks to its affinity criteria DADA can introduce communication in the scheduling without too precise communication cost model which are required in HEFT to predict the completion time of tasks.

We would like to extend the experimental evaluations on robustness of scheduling with respect to uncertainties in cost models, especially on the communication cost which is very sensitive to contentions that may appear at runtime. Another interesting issue would be to study other affinity functions.

i 3 end 4 6 Schedule

 46 Sort all ready tasks Ti by decreasing speedup Si 5 foreach Ti ∈ LR do Ti on the worker wj achieving the earliest finish time 7 push of Ti into queue of worker wj 8 Update processor load time-stamps on worker wj 9 end

4 λ 5 begin local affinity phase 6 8 begin global balance phase 9 Schedule

 45689 ← (upper + lower)/ 2 Schedule tasks of LR per affinity score on its affinity processor, loading each processor up to overreaching αλ 7 end LR to minimize finish time using λ as hint 10 if tasks do fit into (2 + α)λ then 11 19 Push each task Ti of LR on queue of worker wj based on the last fitting schedule and update processor load time-stamps

  Performance of DADA(α).

  Memory transfer of DADA(α).

  Memory transfer of DADA(α)+CP.

Fig. 1 .

 1 Fig. 1. Impact of parameter α on Cholesky (DPOTRF) with matrix of size 8192 × 8192.

Fig. 2

 2 reports the behavior of the Cholesky decomposition (DPOTRF) with respect to the number of GPUs used. It studies both performance results (Fig.2(a)) and total memory transfers (Fig.2(b)). All scheduling algorithms have similar performances. DADA(α)+CP scales slightly better with the number of GPU. As expected DADA(α)+CP and DADA(α) Performance (8192 × 8192).

  Memory Transfer (8192 × 8192).

Fig. 2 .

 2 Fig. 2. Benchmarks of Cholesky (DPOTRF).

  Performance (8192 × 8192).

  Memory Transfer (8192 × 8192).

Fig. 3 .

 3 Fig. 3. Benchmarks of LU (DGETRF).

  Performance (8192 × 8192).

  Memory Transfer (8192 × 8192).

Fig. 4 .

 4 Fig. 4. Benchmarks of QR (DGEQRF).

http://kaapi.gforge.inria.fr

Acknowledgments

This work has been partially supported by the French Ministry of Defense -DGA, the ANR 09-COSI-011-05 Project Repdyn and CAPES/Brazil.