# Porting Parallel Applications to Heterogeneous Supercomputers: Libraries and Tools Can Make it Transparent

Jean-Yves Vet, Patrick Carribault, Albert Cohen

► To cite this version:

Jean-Yves Vet, Patrick Carribault, Albert Cohen. Porting Parallel Applications to Heterogeneous Supercomputers: Libraries and Tools Can Make it Transparent. Compiler, Architecture and Tools Conference (CATC), Sep 2016, Haifa, Israel. hal-04955498

HAL Id: hal-04955498

https://hal.science/hal-04955498v1

Submitted on 19 Feb 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Porting Parallel Applications to Heterogeneous
Supercomputers: Libraries and Tools Can Make it
Transparent

Jean-Yves Vet¹, Patrick Carribault¹, Albert Cohen²

¹ CEA, DAM, DIF, F-91297 Arpajon, France

² INRIA and Ecole Normale Supérieure, Paris, France

Abstract

This paper presents an integrated approach to port existing
parallel applications to heterogeneous supercomputers. Our
techniques apply transparently to high-performance library
calls and parallel regions in applications. Following a hier-
archical bulk-synchronous task decomposition, one may in-
crementally optimize data locality and performance of realis-
tic workloads combining multiple scheduling policies across
dynamic computation phases. Modifying a small number of
lines of code, we propose an integrated solution to orchestrate
the hardware acceleration of an existing OpenMP or multi-
paradigm MPI+OpenMP application.

1. Introduction

Multi-paradigm parallel programming models, such as the
combination of MPI and OpenMP, are popular in supercom-
puting applications. Yet compute nodes continue to evolve, to
grow in complexity, embedding multiple memory spaces and
heterogeneous compute units. To reach high performance,
a regular multi-paradigm MPI+OpenMP is no longer suffi-
cient since spatial locality should be strengthened and hetero-
geneous compute resources harnessed simultaneously. But
many legacy codes are based on MPI and/or OpenMP and are
composed of hundreds of thousands lines of code. It may nei-
ther be worth it nor cost effective to rewrite those applications
to fully take advantage of new architectures. In this context,
we believe scientists will need a framework to incrementally
adapt their applications to better support several kinds of up-
coming supercomputer nodes.

We study a typical heterogeneous compute node from the Tera-100 and Curie supercomputers in the TOP500; Figure 1 presents an abstract view of one such node. It is composed of compute units (two CPU sockets and two NVIDIA Graphics Processing Units (GPUs) named accelerators in this context) and different memory spaces. This architecture exhibits two different hardware strategies to deal with the integration of those memory spaces:

1. Memory with hardware coherency. This leads to Non-Uniform Memory Access (NUMA) architectures: the cores packaged on the same CPU die experience a higher

<figure><img src="image_2.png" /><figcaption>Figure 1: Heterogeneous and hierarchical compute node with NUMA and NUIOA effects</figcaption></figure>

latency when accessing data residing in memory space
managed by another CPU socket. Figure 1 illustrates
this NUMA effect between the two multicore processors
(NUMA node).

2. Memory without hardware coherency. Heterogeneous architectures frequently involve non-coherent memories and distributed address spaces. Feeding the functional units of accelerators, e.g., the Intel Xeon Phi or high-end GPUs, requires distinct memory technology such as GDDR memory, offering several times the bandwidth of DDR. But only a limited amount of this fast memory is placed on the compute node, because of its cost. Data residing in the memory of an accelerator may be transferred to the memory spaces managed by the CPUs or belonging to other accelerators, but those operations should be triggered (manually or automatically) by software. Figure 1 illustrates this configuration between the two GPUs with the notion of Non-Uniform Input/Output Access (NUIOA) [14] nodes.

Complexity introduced by such architecture motivates the
design of runtime interoperability frameworks and heteroge-
neous execution models [14, 15, 23, 26, 44]. But, to the best
of our knowledge, those frameworks do not answer one key
question for legacy codes: how to adapt an existing parallel
application, with its high-performance libraries and kernels,
to take advantage of a hierarchical and heterogeneous dis-
*tributed system? And as a derivative question: how to make the transition as smooth and automated as possible for the scientists and maintainers of the application?*

Thus, this paper puts together and extends existing approaches to augment parallel applications in order to fully exploit heterogeneous compute nodes. In particular, we integrate mechanisms to express task-based parallelism, work-data association, to support automatic load balancing across architectures, and to distribute the computations while taking advantage of temporal and spatial locality, avoiding costly data transfers and replication. Compared to related work, we extend existing approaches to handle both NUMA/NUIOA effects and to deal with multiple phases inside an application. It leads to the design of a hierarchical and heterogeneous task scheduler called **H3LMS** —*Harnessing Hierarchy and Heterogeneity with Locality Management and Scheduling*— and a software memory manager called **COM-PAS** —*Coordinate and Organize Memory Placement and Allocation for the Scheduler*. We implemented H3LMS and COMPAS into the *MultiProcessor Computing* (MPC) framework [21]. MPC is based on the coupling of a thread-based MPI and an OpenMP implementation. This transparent shared-memory approach facilitates the interactions with a task-based programing model on a compute node. Studying a realistic supercomputing application, we show that our integrated solution is particularly well suited to the transparent acceleration of scientific library calls and well-tuned kernels developed as separate modules.

The paper is organized as follows: Section 2 discusses related work on heterogeneous scheduling. Sections 3 and 4 present our runtime mechanisms. Section 5 exposes the scheduling policies. Section 6 describes the implementation and the experimental results, before we conclude in Section 7.

## 2. Related Work

One of the first approaches to the dynamic management of heterogeneous compute resources is called XKAAPI — *Kernel for Adaptive, Asynchronous Parallel and Interactive programming interface*. Based on the *Athapascan* programming model [25], it provides a task-based system with low-overhead creation, optimized synchronization and fine-grain scheduling. The design is based on collaborative work stealing [17] using dedicated locks [28] to reduce local invalidations. The runtime allows nested tasks but no specific load balancing for tasks with different levels. This platform has first been extended to support accelerators [29]: the scheduling decision is driven by *weight*. Task stealing is performed according to locality inspired by A. Umut et al. [8] based on accesses to the same *partition*. But it does not take into account the complex NUMA aspect of such nodes. The second extension [26] brought software caches and LRU policy for line eviction. The scheduling strategy has been updated to take into account cache contents.

OMPSs [22] gathers the developments of the StarSs model—*Star SuperScalar* [43]: it relies on OpenMP-style directives to declare heterogeneous tasks. The source-to-source compiler *Mercurium* [16] lowers the directives while the user-thread library *Nanos++* [47] takes care of data transfers and scheduling. StarSs leads to several sub-projects (SMPSs [38], CellSs [39, 41], ClusterSs [44] and GPUSs [15, 23]). The latter is an extension of CellSs to support GPUs: a new directive allows to specify the target architecture when declaring a task. Like OpenACC 2 or OpenMP 4, there is no support for dynamic heterogeneous load balancing. Finally, software caches control the data retention in remote memories and follow the LRU eviction policy.

StarPU [13, 14] is a software stack integrating a task-scheduling platform for heterogeneous nodes. With the help of C extensions to define *codelets*, it allows to dynamically schedule tasks between different target architectures. Software caches with MSI coherence protocol and LRU eviction policy increase the temporal locality of remote memories. Good performance results have been demonstrated on dense linear algebra [11]. Scheduling policies can be selected among a set of strategies. Most of them are variants of the HEFT approach (*Heterogeneous Earliest Finish Time* [46]) distributing the tasks on the worker that will finish its execution early. For this purpose, statistics are collected to estimate the execution time of each task [12] and data transfers. Based on greedy work distribution, there is no strategy to trade the performance of one kernel in favor of the following one, and there is no work stealing beyond the first task distribution. Also, there is no mechanism to increase spatial locality by taking into account the NUMA aspects of the target node.

Charm++ is a C++ library providing optimized mechanisms for inter-nodes communications. An integrated scheduler execute parallel objects called *chares* [31]. Extensions have been developed to support CELL processors [33] and GPUs [32]. Minas is a library designed to automatically or explicitly place data in memory of a system with NUMA effects. Many policies are available to distribute pages on the NUMA nodes. Positive impacts have been studied when used in conjunction with Charm++ [40]. We propose to go one step further and design such a library with COMPAS so that meta-information may be retrieved by a task scheduler. Indeed, those information should be taken into account when tasks are distributed to the thread workers.

Gelado et al. present GMAC, an (Asymmetric Distributed Shared Memory) (ADSM) designed to reduce the coherence traffic and make data management easier to program with heterogeneous systems [27]. Page faults are used to eagerly transfer data. However, the overhead of page faults may be avoided with task level parallelism by associating data with in/out information (e.g., StarPU “lets” or StarSs pragmas).

Finally, our paper focuses on the integration of exiting mandatory techniques (e.g., task-based mechanisms like XKaapi and StarPU) to exploit heterogeneous resources in-
side parallel multi-paradigm applications. This differentiates our contributions from the above algorithm- and language-centric works. Since we consider full applications with dynamic and phased behavior, our findings are complementary to those of runtime systems for heterogeneous implementations of numerical libraries [3].

# 3. Multi-Granularity Spatial Locality

Compute nodes exposing heterogeneous resources usually have physically distinct memories and no hardware coherence. Even if software coherence is possible, taking care of data locality is crucial to reach high performance. The picture is complicated with the NUMA topology and cache level hierarchies. We propose runtime mechanisms to exploit spatial locality in such nodes.

## 3.1. Multi-Granularity Compute Tasks

We leverage task parallelism to balance the load between heterogeneous resources. Each task captures the following information: (i) the set of data it accesses, and (ii) one function pointer per target architecture. The latter can be programmed by the developer, extracted from an optimized library, or outlined from the body of a task directive by the compiler like in OpenACC or OpenMP 4.0.

Decomposable tasks allow to adapt the workload to the target architecture: number of cores, compute power of each core, multiple levels of parallelism available, etc. We introduce the notion of *super-task* illustrated in Figure 2. When a super-task is scheduled on a multicore CPU, it spawns small tasks to feed the cores with finer grain parallelism. On the other hand, when scheduled on a GPU, the dedicated programming model is in charge of mapping a specific super-task implementation on the accelerator.

<figure><img src="image_3.png" /><figcaption>Figure 2: Runtime decomposition of a super-task</figcaption></figure>

From our runtime point of view, the super-task is not further decomposed when running on a GPU. The choice of the right decomposition granularity is a tough issue that is not addressed in this document. It is a tradeoff between load balancing, data distribution, data transfers and raw performance of nested computations. For example, the PLASMA library relies on the QUARK scheduler [49] and includes the right granularity for linear-algebra kernels. This information is exploited by the MAGMA library to allow an heterogeneous scheduling with data decomposition. Machine learning techniques can be used to automatically find the right set of parameters to reach higher performance. Libraries like

ATLAS [1, 48], PHiPAC [6, 18], FFTW [24, 36] and Spar-sity [5, 30] rely on this approach.

<figure><img src="image_4.png" /><figcaption>Figure 3: Basic flat heterogeneous execution model</figcaption></figure>

Figure 3 illustrates the basic execution model of our heterogeneous scheduler. This example exposes *N* CPU cores and one accelerator (e.g., GPU). The runtime creates one worker thread per core: small units represent the worker that will execute some work on CPU cores while the accelerator unit uses one core to deal with the underlying accelerator. The generated super-tasks are stored inside a shared list common to each compute unit located on the node. Each worker fetches first one super-task from this list. It is then decomposed into smaller tasks and distributed to the other units of the same type. Work stealing is available between these units.

## 3.2. Hierarchical Affinity

The concept of super-task is useful to deal with multiple compute units controlling different architectures. But the execution model in Figure 3 is too simplisitic for complex nodes: a shared task list may lead to performance loss because of contention and the lack of spatial reuse. To increase spatial locality of super tasks, we present an abstract organization, based on the hardware topology to determine the number of super-task lists and their location on compute nodes.

<figure><img src="image_5.png" /><figcaption>Figure 4: Compute units, teams and poles of a compute node composed of two CPUs and one accelerator</figcaption></figure>

Figure 4 illustrates our abstract organization of logical units. These are driven by CPU cores and execution flows are depicted as threads with a *S*-shape (as already detailed in Figure 3). Small gray boxes inside each unit represent a list of local small tasks to be executed by the corresponding core.
Notice that the size of those boxes are proportional to the granularity of the tasks executed by the corresponding unit. This organization exposes multiple super-task lists to maximize node locality when creating and scheduling super-tasks. We define the following groups in this abstract representation:

* **Compute unit** (constant granularity). Two kinds of logical units come into play: *small* and *accelerator*. The first is associated to one specific CPU core; it requires a step of decomposition into smaller tasks before execution. The second type manages an accelerator like a GPU. It is linked to a CPU core as well to drive the corresponding device (transferring data, launching kernels).

* **Work team** (shared physical memory). Logical units sharing the same physical memory are grouped into work teams, made of units of the same type and granularity (either small or accelerator). Two GPU devices may not belong to the same team because their memories are isolated. Each CPU socket also forms one team, assuming each processor represents one NUMA node.

* **Work pole** (memory affinity). A pole gathers teams with the largest affinity to NUMA or NUIOA memory. The goal is to limit, within a pole, the penalty to access distributed memories. Each pole controls as many lists of super-tasks as there are NUIOA nodes inside the pole. Teams associated with units with a uniform access to the main memory can get work from any of those lists. On the other hand, teams related to a device memory have privileged access to one specific list. Section 4 will show it improves temporal locality.

* **Pole hierarchy** (load balancing). Finally, poles are linked together regarding the NUMA factor characterizing their communication. As a result, poles form a hierarchy, favoring local accesses. This hierarchy will be used to drive work stealing for load balancing.

To create this abstract organization, we need to discover the hardware topology of the compute node. We rely on the *HWLOC* library [4, 19] to model nested memories, cache layout, NUMA configuration and factors, affinity with accelerators, etc. This organization represents the building blocks of the heterogeneous runtime platform H3LMS.

### 3.3. Initial Page Placement

This abstract organization favors affinity of super-tasks regarding the NUIOA aspect (links between CPUs and GPUs) and the NUMA layout. Super-tasks are placed inside the most local list and work stealing takes into account the spatial locality by following this representation. The scheduler needs to take into account the NUMA placement of pages to determine the best memory affinity of a super-task, and to steer it to the right team, reducing the number of remote accesses. The *libnuma* [34] library can enforce a memory allocation policy onto NUMA nodes (e.g., round robin). But in our context it is better to keep the control of data allocation and

placement to fully take into account spatial locality in our heterogeneous runtime. For this purpose, we introduce our module named **COMPAS—Coordinate and Organize Memory Placement and Allocation for Scheduler**. Its role is to drive the main heterogeneous scheduler according to spatial locality. COMPAS allows to track initial data placement with high-level patterns dedicated to the task scheduler. Information about initial NUMA affinity are stored in meta-data associated with data blocks manipulated by super-tasks, which are then assigned to a pole that minimizes the number of remote NUMA accesses. The choice of the super-task list inside a pole depends on parameters described in the next section. For that purpose, it is necessary to consider the block size (amount of data manipulated by the super-task) and the kind of access (read, write or read/write).

## 4. Temporal Locality for Data Management

We present an automatic data management scheme to avoid spurious transfers and optimize temporal locality. We first discuss mechanisms needed to maintain data coherence with the extension of the *COMPAS* module. Then, we present the advantages of designing a global eviction policy for data stored in software caches, coupled to the main task scheduler.

### 4.1. Data Management and Layout

To enable automatic data management, the application developer must describe the data blocks associated with the coherence granularity and super-task access patterns. Each block should contain enough data to avoid issues of transferring small amounts on dedicated buses. Furthermore, the placement should be in sync with the initial decomposition described in Section 3.3.

Figure 5 presents the updated COMPAS module to handle automatic data management. COMPAS handles the initial placement of memory pages according to the developer's hints or on the programming model. This is represented by the *E* box in Figure 5. With this allocation, the *resolver* is in charge of managing the allocated memory in blocks and can retrieve the block corresponding to one virtual address (Figure 5-F). Then, each data block is associated with metadata representing the description of this block on each target memory as depicted in box A of the figure.

The coherence state of each block is managed by the *directory* part depicted in Figure 5-B. It keeps track of which block is up-to-date in which memory on the whole compute node. It relies on a simple coherence protocol handling multiple NUMA nodes and device memories with a *write-back* policy. It is inspired from the *MESI* protocol with slight differences. the invalid state (I) is replaced by the delete state (D) for remote memories to avoid wasting space.

This mechanism enables automatic data management with the task scheduler: when creating a task, the virtual address
<figure><img src="image_6.png" /><figcaption>Figure 5: COMPAS data management</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>Figure 6: Scheduling strategies on 2-block memory</figcaption></figure>

ranges are provided with a read and/or write status. The re-
solver is then in charge to determine which blocks are in-
volved, according to these addresses. Once the blocks are
identified, it is possible to know where each block is cur-
rently located and in which state with the help of the directory
storing the coherence protocol of each block. This informa-
tion is useful to determine in which pole and list a super task
should be placed when created (see Section 5). Finally, when
scheduling the super-task, a team is chosen, associated with
a target architecture: by checking the coherence state of this
architecture for each block manipulated by this super-task, it

is possible to know if those blocks are present or not and from
where is should be transferred. Coherence states are updates
accordingly after transferring the necessary data.

Such automatic data management has several advantages:
(i) the application developer does not have to handle the data
transfer, it is just necessary to express the task parallelism
with the parts of the memory that will be manipulated by this
task (like in OpenACC or OpenMP 4), (ii) load balancing is
possible because the runtime system is aware, at any time,
where each block is and in which state and (iii) when trans-
ferring a block to another memory, the data layout can be
adapted to increase the performance of the super-task (this is
the role of the sub-part D in Figure 5).

## 4.2. Software Cache for Device Memories

The COMPAS module now enables automatic data manage-
ment on a heterogeneous compute node. But the different
memories do not have the same bandwidth and capacity. It
is necessary to deal with the overcapacity issue when trying
to transfer some data on a device whose memory is already
saturated. For this purpose, a *software cache* can be used: we
attach one software cache per remote memory (usually related
to one accelerator) to maintain the data copies of blocks that
currently reside in this memory. Linked to the coherence pro-
tocol described in the previous section, it reduces the number
of data transfers and increases the performance of kernels that
exploit the temporal locality. All remote memories are there-
fore handled with software caches, disabling the possibility to
manually manage them. When the cache is full, it means that
the capacity of the target memory is reached. Blocks should
therefore be evicted from this memory to make room for new
blocks needed for a future super-task. Thus, one should select
a *victim block* to be removed from the cache: the *Least Re-
cently Used* (LRU) policy is generally used. We demonstrate
here the significance of the interaction between the software
caches and the execution order of super tasks, leading to a
new scheduling policy based on cache statistics.

When the software cache is full, the eviction policy has a
crucial role to keep the temporal data locality and therefore
maintain the performance. Wrong victim choices may lead
to additional data transfers in the future and therefore perfor-
mance loss because some data have been evicted before they
were reused. Software caches usually rely on eviction poli-
cies based on information gathered during the application ex-
ecution to determine the next victim. The LRU policy is one
example: every time a block is used by a kernel scheduled on
the corresponding target architecture, it is marked as used re-
cently. Therefore, this block will not be the next victim. The
*reuse distance* is thus a metric that is used to select a victim.
But this reuse distance depends on the super-tasks execution
order. Inside one bulk of tasks, between two synchronization
barriers, the main scheduler has the opportunity to choose
any of the created super-tasks to schedule it on the target ar-
chitecture. Figure 6 illustrates this issue by presenting two different strategies with one remote memory containing only two blocks. On the left-hand side, greedy scheduling with an LRU metric determines the next victim. On the right-hand side, a locality-driven scheduler is used: this strategy takes into account the LRU value of each block located inside the memory to decide which kernel to execute next. Both schedulers decide first to execute the super-tasks manipulating A, then B, C and finally D before reaching the first synchronization barrier. While the greedy scheduler performs identically after this synchronization, the locality-aware scheduler will first executes the super-task touching the block D to minimize its reuse distance. Then the task corresponding to C is scheduled to exploit the temporal locality with the previous bulk. Finally, the remaining tasks (using A then B) are launched but they both have been evicted and therefore need to be transferred back. The bottom part of the figure gives the statistics on the total number of transfers from main memory to the device and the other way around. The locality-aware greedy scheduler is able to reduce the amount of transferred data by exploiting the temporal locality.

The scheduling strategy needs to take into account the current state of software caches to decide which task should be executed on which units. We implement this method by sorting super-tasks in lists. As described in Section 3.2, lists of super-tasks are created between two teams depending on the NUIOA affinity. Furthermore, Figure 3 page 3 detailed the main basic scheduling algorithm that is used in our runtime to execute the super-tasks of a list. The scheduling algorithm we propose is a mix of work distribution and work stealing. Given the bulk granularity, it can be seen as performing a fork/join operation every time we execute a bulk (this can be compared to an OpenMP parallel region with a dynamically-scheduled loop). Thus, at the beginning of each bulk, the runtime system is aware of all the super-tasks that will have to be executed before reaching the barrier, exiting this programming model. Because of this global view, the runtime can take some time to expose the super-tasks to drive the scheduling process. This is the main purpose of sorting each task list: in each list, the super-tasks are sorted according to their affinity with the corresponding memory.

<figure><img src="image_8.png" /><figcaption>Figure 7: Super-task affinity</figcaption></figure>

Figure 7 illustrates this concept with one super-task list and two workers: a small unit (CPU core) on the left-hand side and an accelerator unit (CPU core for a GPU) on the right-hand side. To drive the scheduler, super-tasks are sorted according to their affinity with either the NUMA node related to the CPU core (on the left-hand side) or the memory corresponding to the GPU on the right-hand side.

Performing this sort before waking up the workers allows to skew the scheduling process by distributing a priori the super-tasks between the CPUs and the GPUs depending on the current data locality (temporal and spatial).

$$A(\text{supertask, list}) = \sum_{i=1}^{\text{blocs}} L(i, \text{cache}) \quad (1)$$

A: Affinity regarding a remote team.

L: 1 if the block is located inside the remote memory, 0 otherwise.

More formally, (1) is used to sort the super-task inside a list. The value A determines if this super-task has a high affinity with the corresponding remote memory modeled through the cache. L(i, cache) will be set to 1 if the block i is already held inside the cache, meaning that no transfer is required.

## 4.3. Multiple Library Calls

To port legacy codes faster our framework embeds compute-intensive libraries ready to be used on the top of H3LMS and COMPAS. They rely on vendor-optimized libraries such as Intel MKL or Nvidia cuBLAS. This section discusses the exploitation of temporal locality when invoking consecutively kernels defined in one or several libraries. Indeed, heterogeneous libraries offer optimized kernels that exploit all computational resources available on the target nodes but they do not offer the possibility to exploit software caches throughout many calls. In homogeneous mode, data are placed inside only one memory. Therefore, consecutive calls to a library do not imply data transfers. For example, the linear-algebra library Nvidia cuBLAS lets the programmer deal with data placement and transfers. But in case of dynamic heterogeneous execution, the data placement is not predetermined because the runtime system is in charge of moving data according to current scheduling for load-balancing. Therefore, the classical way of using such libraries is to guarantee the data coherence inside the main memory before and after invoking a kernel. The library runs the kernel with load-balancing, moving data around. After exiting the kernel, modified data are transferred back to the main memory to maintain the coherence for the rest of the program execution. This mode adds extra communications between memories and involves a large number data transfers.

MORSE [2] is an interface designed to link runtime systems to linear-algebra libraries. This kind of interface increases the collaboration between heterogeneous runtime designers and kernel developers. However, the interface do not allow to benefit from software caches between consecutive kernel calls. We propose to wrap libraries we want to use with an interface similar to MORSE. Our solution offer a temporal locality feature by using COMPAS to detect which data blocks should be transferred. Moreover, heterogeneous kernels may be defined as well as super-task and sub-task decomposition. It allows: (i) to rely on the same runtime system for
the library and for the application and (ii) to benefit from the software caches between consecutive kernels. Performances are evaluated in Sections 6.2 and 6.4.

# 5. Locality-Driven Scheduling

Until now, the H3LMS scheduler is always based on work stealing. But this can be counterproductive: while keeping the focus on load-balancing will increase the performance of the current bulk, moving data around to finish the work earlier, the next kernel may perform poorly because the data layout has changed, inducing additional movement. To tackle this issue, we propose two scheduling strategies: the first one minimizes the time to execute the current super-tasks while the second one maximizes the data locality.

## 5.1. Scheduling Strategy for Load Balancing

When a super-task is created during the bulk instantiation, H3LMS should decide in which list it should be added. The super-task has then to be inserted inside the list according to (1). This first strategy distributes work in a way that maximizes data locality when placing the super-task in the list and then the worker executes these tasks and steals work from another list. To pick the right list, we apply a metric inspired from the owner computes rule [37, 42]. Because super-tasks should contain a large amount of work (suitable for one accelerator or a full CPU processor), one approach is to choose the execution units which minimize the amount of data to transfer, even if this computation has a cost. Transferring data from one memory to another is certainly more costly than computing metrics to decide where to place a super-task.

When a super-task is generated, the resolver part of COM-PAS translates the address ranges to blocks. Then, there are two possibilities: (i) one of these blocks is located inside one software cache on the node or (ii) all the blocks are up-to-date only in the main memory. For the first case, we define the $R_1$ index associated with caches. For each software cache, this index is computed and the choice is made according to the accelerator which maximizes $R_1$. Blocks carrying a read/write dependency are be more important than read-only or write-only blocks because they might entail more traffic: one transfer before executing the super-task for the read access and one transfer back to update another memory if necessary. Therefore, the super-task list associated with the accelerator with the largest $R_1$ index for this super-task will include it. The following formula details the computation of $R_1$:

$$R_1(\text{supertask, cache}) = \sum_{i=1}^{B(\text{supertask})} S(i) \times D(i) \times L(i, \text{cache}) \qquad (2)$$

$R_1$: Distribution index for one super-task on a software cache.

B: Set of blocks manipulated by the super-task.

S: Size of the block (in bytes).

D: Dependence type (read only: 1, write only: 1, read/write: 2).

L: 1 if the block is present in the cache, 0 otherwise.

On the other hand, if all blocks are located in main memory (i.e., $L(i, cache) = 0$ for each block $i$ and for each cache), then we define the $R_2$ index, computed as to avoid remote accesses by identifying the right target pole.

$$R_2(\text{super-task, pole}) = \sum_{i=1}^{B(\text{super-task})} S(i) \times D(i) \times N(i, \text{pole}) \qquad (3)$$

$R_2$: Distribution index for one super-task on a pole.

B: Set of blocks manipulated by the super-task.

S: Size of the block (in bytes).

D: Dependence type (read only: 1, write only: 1, read/write: 2).

N: NUMA presence (1 for an up-to-date block in the NUMA node corresponding to pole, 0 otherwise).

Figure 8 illustrates the main super-task distribution strategy to fill the different lists according to the indices defined in this section: $R_1$ and $R_2$. The right-hand side depicts an abstract heterogeneous architecture. It embeds two 12-core AMD Opteron 6164HE processors and two Nvidia Geforce GTX 470 accelerators. Notice that these two GPUs are connected to the same NUIOA node. For the representation, there are two poles (one per socket) and each pole contains at least 2 teams representing the NUMA nodes inside the processors. The first pole dedicates two cores for accelerator workers, leading to the allocation of two super-task lists: one linked to each software cache corresponding to the accelerators. The left-hand side presents one super-task that has been generated manipulating four blocks: $B1$, $B2$, $B3$ and $B4$. The distribution algorithm is described below the super-task definition: the first check is related to the presence of at least one block inside one of the two software caches. If this is correct, the $R_1$ index is computed for each software cache (for lists $LS1$ and $LS2$ only because the $LS3$ list is not connected to an accelerator). It leads to the results $V11$ and $V12$ which determines the target list (or a random choice if $V11 = V12$). If none of the blocks are present inside one software cache, the $R_2$ index is evaluated for each pole. It results $V21$ and $V22$ determining the target pole. If this pole is the first one (Pole I), then a random list is chosen between $LS1$ and $LS2$.

## 5.2. Scheduling Strategy for Locality

The previous strategy computes indices to choose the right target list according to the memory-block placement for one super-task. But the performance of device processors are lowered by costly data transfers. Those movements are amplified by data going back and forth during different compute phases: some parts of the application are compute intensive while other parts are more suitable for CPUs. Still, scheduling this second category only on CPU cores would involve data transfers and may compromise the performance of the next compute-intensive kernel.
<figure><img src="image_9.png" /><figcaption>Figure 8: Super-task creation and scheduling</figcaption></figure>

Thus, we propose another task distribution that is suitable for non-compute intensive phases. To maximize data reuse, this distribution hampers load balancing and performs only the task distribution that we explained before by putting the tasks directly into local list depending on the availability of blocks inside a memory. This strategy will be highlighted in Section 6.4 on the PN benchmark.

# 6. Implementation and Experiments

We integrated our mechanisms to MPC to support existing applications parallelized with MPI and OpenMP. We present four experimental results: (i) an increased spatial locality on heterogeneous single-precision matrix multiply operations (SGEMM) (ii) an increased temporal locality with multiple SGEMM calls (iii) High-Performance Linpack to test porting an application to benefit from heterogeneous capabilities of MPC+H3LMS, and (iv) PN, an application from the french *Atomic Energy Commission* (CEA) showing the benefits of hybrid scheduling strategies, as described in Section 5.

## 6.1. Spatial Locality on SGEMM

Figure 9 depicts the execution time to perform single-precision matrix multiply operations (SGEMM kernel) on the heterogeneous platform illustrated in Figure 1.

We compare here a state-of-the-art library MAGMA [3, 9, 10, 45] with our approach: H3LMS (with super-tasks) with or without the COMPAS module to drive and keep track of memory pages. For each version, the kernels rely on optimized dedicated libraries: *Nvidia cuBLAS* for the Nvidia GPUs and *Intel MKL* for the Intel Xeon CPUs. Matrices are split in 1024 × 1024 blocks to parallelize and dispatch the work between the compute units. For the H3LMS curve, super-tasks with blocks of size 1024×1024 are created for the accelerator unit, decomposed in 256 × 256 blocks for small

units (CPU cores). The MAGMA curve is a single granularity scheduler which only considers 1024 × 1024 blocks using the MAGMA kernel. This library relies on the StarPU scheduler and the task definition is done with the MORSE interface [2] which is designed to link runtime systems to linear-algebra libraries like MAGMA.

Finally, the last curve (H3LMS + COMPAS) presents our multi-granularity scheduler with a memory page distribution on NUMA nodes. For each matrix, the initial placement follows the super-task granularity: each block corresponds to one memory page and those blocks are allocated in a round robin fashion on the NUMA nodes. Super-tasks are assigned to the pole where the read/write matrix lies (i.e., *C* matrix considering the operation *C* = *C* + *A* × *B*). Thus, half of the blocks are assigned to each NUMA node, and super-tasks are evenly assigned to each pole. It minimizes the transfers needed to perform the sub-operations.

For small matrices, MAGMA outperforms H3LMS because tasks are only executed on accelerators. This short-cut mechanism could be applied to our super-task scheduler when the amount of work is too small. For larger matrices (6144 × 6144 and above), super-tasks allows to reach higher performance (H3LMS curve). The multiple granularities allows to have a better load balancing between each computing resource with only a small amount of information. This approach gains up to 20% more performance compare to the state-of-the-art library (MAGMA). Using COMPAS to have a high-level control on memory placement and to drive the super-task assignment (on the abstract organization: to pole and lists) adds 5% of performance.

## 6.2. Temporal Locality with Multiple Library Calls

To study the possibility to exploit temporal locality between consecutive library calls, we explored the extension of the MORSE interface to retain block information. It allows: (i)
<figure><img src="image_10.png" /><figcaption>Figure 9: Heterogeneous SGEMM, 2 CPUs and 2 GPUs</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>Figure 10: 5 SGEMMs accumulating in the same matrix</figcaption></figure>

to rely on the same runtime system for the library and for the application and (ii) to benefit from the software caches between consecutive kernels.

Our experiment calls the SGEMM kernel 5 times, accumulating the results inside the same *C* matrix.

Figure 12 shows the results: the regular MORSE interface with StarPU without exploiting temporal locality between two SGEMM kernels reaches just above 600 GFLOPS. The second curve exhibits the modified MORSE interface to perform data retention between the SGEMM kernels and the MAGMA library relying on StarPU for compute tasks. This modification reaches 1 TFLOPS: data retention yields speedup of 66% with modified MAGMA library. The two other versions are with H3LMS and COMPAS which naturally keep the temporal locality as we detailed in this section. The difference between these two flavors is the initial data placement on one NUMA node or one the two NUMA nodes in a round robin fashion. For H3LMS, super-tasks are using 1024 × 1024 blocks (like MAGMA) and are decomposed into 256 × 256 blocks for small tasks inside the same team. Thus, the multi-granularity scheduler allows to gain more performance especially if the initial placement of memory pages is spread on the NUMA nodes to drive spatial locality. Thus we can reach up to 1 139 GFLOPS leading to an overall improvement

of 87% compared to regular original MAGMA library, without data retention.

## 6.3. High-Performance Linpack (HPL)

The performance test *HPL* (*High Performance Linpack* [7]) solves a dense linear system $Ax = B$ by inversing matrix *A* with an optimized factorization. This benchmark is used to compare supercomputers for the *TOP500*.

The purpose of this evaluation is to check the relevance of our approach inside an existing application with very few modifications. Starting with the optimized CPU version of Linpack (*HPL 2.1* [7]), we modify the calls to the BLAS library to perform heterogeneous matrix multiplies. Parts of the code outside the DGEMM kernels are statically executed on a CPU core outside H3LMS. Therefore, data coherence ensures that up-to-date versions of each block are available inside the main memory before finishing a BLAS 3 operation. This highlights spatial locality but not temporal locality, since the code is not managed by H3LMS+COMPAS between multiple consecutive calls to the library.

We modified only 3 lines of code for this purpose as shown in Figure 11. First of all, the allocation of matrices were updated to let COMPAS pins the memory pages to increase the transfer bandwidth. Then, the function calls for BLAS 3 kernels are replaced by synchronous functions linked to library based on H3LMS. This library does not exploit the full software-cache features because our platform is unaware of the rest of the code outside BLAS 3 operations. Coherence is thus done at the end of each call: every block located on a GPU is transferred back to the main memory. For better performance, our matrix-multiply library relies on optimized kernels in *Intel MKL 12* and *Nvidia CUBLAS 4.2*.

include/hpl.blas.h

172 #define HPL_dgemm (...) h3lms.blas.dgemm_ext.sync(_VA_ARGS..., 4096,-1,256,-1)

173 #define HPL_dtrsm h3lms.blas.dtrsm_sync

testing/ptest/HPL_pdtest.c

164 vptr = (void*)compas.malloc(NULL, COMPAS_PLOCKED, ((size_t)(ALGO->align) + (size_t)(mat ld+1) * (size_t)(mat nq)) * sizeof(double));

<caption>Figure 11: HPL 2.1 code updates</caption>

We ran the *Linpack* application on one heterogeneous node of Tera-100 with 2 *Intel Xeon* quad-core processors, 2 *Nvidia Tesla M2090* GPUs (see Figure 1), and 24 GB of main memory. We compared 4 different scenarios:

1. Sequential execution on one CPU core. BLAS3-based kernels rely on the sequential version of *Intel MKL 12* library

2. Parallel homogeneous execution on 8 *Intel Xeon* core. Here, the parallel version of *Intel MKL 12* library is used based on the Intel OpenMP runtime.

3. Parallel homogeneous execution with MPC+H3LMS on the 8 CPU cores, using the H3LMS-based multi-granularity
<figure><img src="image_12.png" /><figcaption>Figure 12: HPL 2.1 Tera-100 heterogeneous node</figcaption></figure>

DGEMM. Super-tasks contain 4096 × *m* data blocks, with *m* the amount of modified data in one column per iteration, and are decomposed into 256 × *m* tasks.

<figure><img src="image_13.png" /></figure>

4 Execution with MPC+H3LMS on 6 Xeon cores and 2 Nvidia GPUs (the 2 remaining cores drive their corresponding GPU)

Figure 12 depicts the results obtained on the previous scenarios. The main HPL parameters are *N* = 46,080, *NB* = 512 (block size), *P* = *Q* = 1, with “methods” WC10L2L2. Matrix elements are stored column-major. Figure 12-b presents the performance of each scenario. The multi-granularity approach with H3LMS leads to a small performance drop because of scheduling overhead. Indeed, the scenario 3 is 5% slower than the MKL version (scenario 2). But the speedup is 7.46× compared to sequential execution (scenario 1), leading to a parallel efficiency of 93.25%. Exploiting all heterogeneous resources (scenario 4) reduces the execution to 135 seconds (482.4 GFLOPS) representing a speedup of 7.88× compared homogeneous parallel execution (scenario 3) and a total speed up of 55.03 compared to sequential execution (scenario 1). By modifying only 3 lines of code, the *Linpack* benchmark may benefit from a good acceleration by using heterogeneous resources with H3LMS+COMPAS.

Notice that heterogeneous performance is about 30% slower than state-of-start static results by Fatica et al. [35] with a larger matrix (108 032 × 108 032 with 96 GB of memory). Indeed, HPL performance is shaped by the size of the input parameters, mainly the matrix size. With a larger matrix, the DGEMM kernels are more time-consuming, increasing the parallel efficiency of the overall execution. In the scenario 4, profiling shows that only 50% of execution time is spent in DGEMM. Increasing the matrix size would definitely gain in performance. Other improvement could be done: for example exploiting a better parallelism of the *TRSM* kernel or tightly coupling to MPI to allow spatial locality between

kernels (including software cache exploitation).

## 6.4. PN Application

PN is a multi-paradigm application developed at CEA solving the linear particle transport equation with deterministic resolution based on spherical harmonics approximation [20]. This application is parallelized with the MPI and OpenMP programming models and exploits the Intel MKL library for matrix computation. It has been selected as a representative benchmark of realistic complexity for the design of future supercomputing architectures and software stacks, and is thus similar in spirit to many applications in the CORAL benchmark suite.<fnref n="1" /> This section describes how we port PN to our runtime platform H3LMS+COMPAS, incrementally accelerating it on complex heterogeneous compute nodes. The main goal is to benefit from all heterogeneous compute resources by minimizing the modifications of the source code (incremental approach).

The algorithm works on a 2D (*x* and *z*) Cartesian mesh. After a profiling phase with input parameter *N* = 15, two functions consume about 86% of execution time: *numerical flux x* and *numerical flux z*. Those functions perform almost the same sequence of operations but on different dimensions: *x* and *z*. But these sequences have no dependencies and can therefore be executed in parallel. The first transformation consists in merging the two bodies into one new function called *numerical flux*. It is composed of 4 main parts:

- Part 1: Two independent matrix multiplies with large matrices.

- Part 2: Four independent matrix multiplies with small matrices.

- Part 3: Consecutive loops operating on each cell of the mesh.

- Part 4: Two independent matrix multiplies with large matrices.

We modified the code to adapt to the H3LMS+COMPAS platform. First, as a preliminary step (B0), COMPAS needs to deal with data allocation and register the blocks for granularity. Thus, 22 lines were modified for the 22 matrices. To drive the spatial locality, data exclusively related to *x* dimension and *z* dimension are placed on different NUMA nodes. The direct effect would be to bias the scheduling against a specific accelerator. Then, during the first step (B1), we modified the calls to the matrix multiplies: the first approach ports only large matrices operations (Part 1 and Part 4) to the heterogeneous library we designed for the HPL benchmark in Section 6.3. 4 lines were modified for the DGEMM calls. Kernels on small matrices (Part 2) are just parallelized for CPU cores since the amount of work is not worth the transfer to GPUs. With this approach, no temporal locality is possible since there is some part of the main algorithm (Part 3) that is not managed by H3LMS. Therefore, the next step (B2) was to port these loops with super-task to allow temporal locality through software caches. Finally, we modified the scheduling

<footnote n="1">¹https://asc.llnl.gov/CORAL-benchmarks.</footnote>
1 /* Part 1 */
2 DGEMM [read: A<sub>X</sub>(136×136), B<sub>X</sub>(x×z×136)] [write: C<sub>X</sub>(x×z×136)]
3 DGEMM [read: A<sub>Z</sub>(136×136), B<sub>Z</sub>(x×z×136)] [write: C<sub>Z</sub>(x×z×136)]
4
5 /* Part 2 */
6 DGEMM [read: A<sub>X</sub>(136×136), B<sub>L</sub>X(x×136)] [write: CL<sub>X</sub>(x×136)]
7 DGEMM [read: A<sub>Z</sub>(136×136), B<sub>L</sub>Z(z×136)] [write: CL<sub>Z</sub>(z×136)]
8 DGEMM [read: A<sub>X</sub>(136×136), B<sub>R</sub>X(x×136)] [write: CR<sub>X</sub>(x×136)]
9 DGEMM [read: A<sub>Z</sub>(136×136), B<sub>R</sub>Z(z×136)] [write: CR<sub>Z</sub>(z×136)]
10 BARRIER
11
12 /* Part 3 */
13 SUPER-TASK [read: D<sub>X</sub>(1024), C<sub>X</sub>(x×z×136)] [write: E<sub>X</sub>(x×z×136)]
14 SUPER-TASK [read: D<sub>Z</sub>(1024), C<sub>Z</sub>(x×z×136)] [write: E<sub>Z</sub>(x×z×136)]
15 BARRIER
16 SUPER-TASK [read: D<sub>X</sub>(1024), CL<sub>X</sub>(x×136)] [read-write: E<sub>X</sub>(x×z×136)]
17 SUPER-TASK [read: D<sub>Z</sub>(1024), CL<sub>Z</sub>(z×136)] [read-write: E<sub>Z</sub>(x×z×136)]
18 BARRIER
19 SUPER-TASK [read: D<sub>X</sub>(1024), CR<sub>X</sub>(x×136)] [read-write: E<sub>X</sub>(x×z×136)]
20 SUPER-TASK [read: D<sub>Z</sub>(1024), CR<sub>Z</sub>(z×136)] [read-write: E<sub>Z</sub>(x×z×136)]
21 BARRIER
22
23 /* Part 4 */
24 DGEMM [read: F<sub>X</sub>(136×136), E<sub>X</sub>(x×z×136)] [write: G<sub>X</sub>(x×z×136)]
25 DGEMM [read: F<sub>Z</sub>(136×136), E<sub>Z</sub>(x×z×136)] [write: G<sub>Z</sub>(x×z×136)]
26 BARRIER
27
28 FLUSH.SOFTWARE_CACHES

<figure><img src="image_14.png" /><figcaption>**Figure 14: Running numericaloil on 1 Tera-100 heterogeneous node with 1536 × 1536 mesh in double precision**</figcaption></figure>

policy to reinforce locality (B3). The result is shown with the pseudocode 13. To sum up we went into the following steps to port the function numericaloil:

A1 Reference - Sequential execution with one 1 CPU core.

A2 Reference - Parallel execution on all CPU cores (Intel OpenMP).

B1 Heterogeneous execution on large DGEMM (Part 1 and Part 4), parallel homogeneous execution for the rest (Part 2 and Part 3). After each large DGEMM, data are flushed out of software caches.

B2 Heterogeneous execution with H3LMS on DGEMM and compute tasks defined by hand (Part 1, Part 3 and Part 4). Temporal locality can be exploited during those steps.

B3 Heterogeneous execution with H3LMS on matrix multiplies and compute tasks defined by hand (Part 1, Part 3 and Part 4). Load-balancing scheduling is used for DGEMM (Part 1 and Part 4) while locality-only scheduling is used for hand-defined tasks (Part 3) to minimize data movements.

C Theoretical - Heterogeneous execution (no transfer).

**Evaluation.** Figures 14-a and 14-b expose the experimental results of these scenarios. The execution time includes the time to perform all the timesteps of the simulation. The A2 version performs only homogeneous parallel tasks bringing a speedup of 5.84× compared to the sequential execution (A1). By exploiting the heterogeneous linear-algebra library for Part 1 and Part 4, the two GPUs execute kernels in parallel with the CPU cores (scenario B1): it leads to an execution time of 84.29 seconds compared to 154.21 seconds using only CPUs. Results for version C are much better because they ignore data transfers, yielding a theoretical performance limit. But this experiment shows that performance is shaped by data movement: in B1, the part of transfer not covered by computation is 61%. Hence performing Part 3 in a heterogeneous way improves performance (about 19% with scenario B2). But allowing aggressive load balancing for the Part 3 may induce a large amount of data transfers. The gain in performance might be better, but it changes the data layout for the next operations (Part 4) which is more time consuming. We tested multiple scheduling strategy on this application leading to version B3: by running Part 3 in a heterogeneous way without load balancing (only driven by locality), the transfer ratio is only 42% of the execution time. Finally, speedup reaches 2.64× when comparing the homogeneous approach A2 with the best heterogeneous scenario B3.

7. Conclusion and Future Work

We presented an integrated solution to port homogeneous
parallel applications to heterogeneous supercomputers, op-
timizing their temporal and spatial locality. We designed
and implemented the H3LMS platform and the COMPAS
module integrated to the production MPC framework. This
framework allows to incrementally wrap high-performance
library calls and application-specific parallel regions, orches-
trating the computations and data transfers on complex com-
pute node topologies. Multiple scheduling policies inside the
same application are needed to deal with the locality and par-
allelism trade-offs, and these policies may depend on the dif-
ferent phases of an application. We demonstrated competitive
performance results through a small number of incremental
changes to realistic multi-paradigm applications.

As a future work, we plan to provide language exten-
sions (e.g., C #pragma) based on OpenACC/OpenMP 4 di-
rectives to describe data allocation patterns and create super-
tasks and inner-tasks with directives. We are also working on
adding some features of the Asymmetric Distributed Shared
Memory GMAC (introduced in Section 2) to our solution.
It would remove the necessity to explicitly flush software
caches when entering a region where data are not managed
by H3LMS+COMPAS.
References

[1] Automatically Tuned Linear Algebra Software (ATLAS). http://math-atlas.sourceforge.net.

[2] Matrices Over Runtime Systems @ Exascale (MORSE). http://icl.cs.utk.edu/projectsdev/morse.

[3] Matrix Algebra on GPU and Multicore Architectures (MAGMA). http://icl.cs.utk.edu/magma.

[4] Portable Hardware Locality (hwloc). http://www.open-mpi.org/projects/hwloc.

[5] SPARSITY. http://www.cs.berkeley.edu/~yelick/sparsity.

[6] The PHiPAC (Portable High Performance ANSI C). http://www1.icsi.berkeley.edu/~bilmes/Philipac.

[7] HPL - A Portable Implementation of the High-Performance Linpack Benchmark for Distributed-Memory Computers - Version 2.1. http://netlib.org/benchmark/hpl, 2012.

[8] Umut A Acar, Guy E Blelloch, and Robert D Blumofe. The data locality of work stealing. *Theory of Computing Systems*, 35(3):321–347, 2002.

[9] E. Agullo, C. Augonnet, J. Dongarra, H. Ltaief, R. Namyst, S. Thibault, and S. Tomov. A hybridization methodology for high-performance linear algebra software for GPUs. In *GPU Computing Gems*. Morgan Kaufmann, 2010.

[10] E. Agullo, J. Demmel, J. Dongarra, B. Hadri, J. Kurzak, J. Langou, H. Ltaief, P. Luszczek, and S. Tomov. Numerical linear algebra on emerging architectures: The plasma and magma projects. In *Journal of Physics: Conference Series*, volume Vol. 180, 2009.

[11] Emmanuel Agullo, Cédric Augonnet, Jack Dongarra, Mathieu Faverge, Julien Langou, Hatem Ltaief, and Stanimire Tomov. LU factorization for accelerator-based systems. In *9th ACS/IEEE International Conference on Computer Systems and Applications (AICCSA 11)*, Sharm El-Sheikh, Egypt, 2011.

[12] Cédric Augonnet, Samuel Thibault, and Raymond Namyst. Automatic Calibration of Performance Models on Heterogeneous Multicore Architectures. In *3rd Workshop on Highly Parallel Processing on a Chip (HPPC 2009)*, Delft, Pays-Bas, August 2009.

[13] Cédric Augonnet, Samuel Thibault, Raymond Namyst, and Pierre-André Wacrenier. StarPU: A unified platform for task scheduling on heterogeneous multicore architectures. In *Euro-Par’09*.

[14] Cédric Augonnet, Samuel Thibault, Raymond Namyst, and Pierre-André Wacrenier. StarPU: a unified platform for task scheduling on heterogeneous multicore architectures. *Concurrency and Computation: Practice and Experience*, 23(2):187–198, 2011.

[15] Eduard Ayguadé, Rosa M Badia, Francisco D Igual, Jesús Labarta, Rafael Mayo, and Enrique S Quintana-Ortí. An extension of the stars programming model for platforms with multiple gpus. In *Euro-Par 2009 Parallel Processing*, pages 851–862. Springer, 2009.

[16] J. Balart, A. Duran, M. Gonzàlez, X. Martorell, E. Ayguadé, and J. Labarta. Nanos mercurium: a research compiler for openmp. In *European Workshop on OpenMP (EWOMP’04)*. Pp, pages 103–109, 2004.

[17] Xavier Besseron, Christophe Laferrière, Daouda Traore, and Thierry Gautier. X-kaapi: Une nouvelle implémentation extrême du vol de travail pour des algorithmes í¿ grain fin. *19èmes Rencontres Francophones Du Parallélisme (RenPar 2009)*.

[18] J. Bilmes, K. Asanović, J. Demmel, D. Lam, and C.W. Chin. PHiPAC: A portable, high-performance, ANSI C coding methodology and its application to matrix multiply. LAPACK working note 111, University of Tennessee, 1996.

[19] Francois Broquedis, Jérôme Clet-Ortega, Stéphanie Moreaud, Nathalie Furmento, Brice Goglin, Guillaume Mercier, Samuel Thibault, and Raymond Namyst. hwloc: A generic framework for managing hardware affinities in hpc applications. In *Proceedings of the 2010 18th Euromicro Conference on Parallel, Distributed and Network-based Processing*, PDP ’10, pages 180–186, Washington, DC, USA, 2010. IEEE Computer Society.

[20] Thomas A. Brunner and James Paul Holloway. Two-dimensional time dependent riemann solvers for neutron transport. *J. Comput. Phys.*, 210(1):386–399, November 2005.

[21] Patrick Carribault, Marc Pérache, and Hervé Jourdren. Enabling low-overhead hybrid mpi/openmp parallelism with MPC. In *Beyond Loop Level Parallelism in OpenMP: Accelerators, Tasking and More, Proc. 6th International Workshop on OpenMP, IWOMP 2010*, Tsukuba, Japan, pages 1–14, 2010.

[22] Alejandro Duran, Eduard Ayguade, Rosa M Badia, Jesús Labarta, Luis Martinell, Xavier Martorell, and Judit
Planas. Ompss: a proposal for programming heterogeneous multi-core architectures. *Parallel Processing Letters*, 21(02):173–193, 2011.

[23] Roger Ferrer, Judit Planas, Pieter Bellens, Alejandro Duran, Marc Gonzalez, Xavier Martorell, Rosa M Badia, Eduard Ayguade, and Jesus Labarta. Optimizing the exploitation of multicore processors and GPUs with OpenMP and OpenCL. In *Languages and Compilers for Parallel Computing*, pages 215–229. Springer, 2011.

[24] Matteo Frigo and Steven G. Johnson. The design and implementation of FFTW3. *Proceedings of the IEEE*, 93(2):216–231, 2005. Special issue on “Program Generation, Optimization, and Platform Adaptation”.

[25] François Galilée, Gerson GH Cavalheiro, J-L Roch, and Mathias Doreille. Athapascan-1: On-line building data flow graph in a parallel language. In *Parallel Architectures and Compilation Techniques, 1998. Proceedings. 1998 International Conference on*, pages 88–95. IEEE, 1998.

[26] Thierry Gautier, Joao Vicente Ferreira Lima, Nicolas Maillard, Bruno Raffin, et al. Xkaapi: A runtime system for data-flow task programming on heterogeneous architectures. In *27th IEEE International Parallel & Distributed Processing Symposium (IPDPS)*, 2013.

[27] Isaac Gelado, John E. Stone, Javier Cabezas, Sanjay Patel, Nacho Navarro, and Wen-mei W. Hwu. An asymmetric distributed shared memory model for heterogeneous parallel systems. In *ASPLOS’10*.

[28] Danny Hendler, Itai Incze, Nir Shavit, and Moran Tzafrir. Flat combining and the synchronization-parallelism tradeoff. In *Proceedings of the 22nd ACM symposium on Parallelism in algorithms and architectures*, pages 355–364. ACM, 2010.

[29] Everton Hermann, Bruno Raffin, François Faure, Thierry Gautier, and Jérémie Allard. Multi-gpu and multi-cpu parallelization for interactive physics simulations. In *Euro-Par 2010-Parallel Processing*, pages 235–246. Springer, 2010.

[30] Eun jin Im and Katherine Yelick. Optimizing sparse matrix vector multiplication on smps. In *In Ninth SIAM Conference on Parallel Processing for Scientific Computing*, 1999.

[31] Laxmikant V. Kale and Sanjeev Krishnan. Charm++: a portable concurrent object oriented system based on c++. *SIGPLAN Not.*, 28(10):91–108, October 1993.

[32] David M Kunzman and Laxmikant V Kalé. Programming heterogeneous clusters with accelerators using object-based programming. *Scientific Programming*, 19(1):47–62, 2011.

[33] David M. Kunzman, Gengbin Zheng, Eric Bohm, James C. Phillips, and Laxmikant V. Kale. Charm++ simplifies coding for the cell processor. In *Proceedings of the 2006 ACM/IEEE conference on Supercomputing, SC '06*, New York, NY, USA, 2006. ACM.

[34] A. K. Leen. A NUMA API for Linux. Technical report, Novell, Technical Linux Whitepaper, 2005.

[35] E. Phillips M. Fatica. CUDA accelerated linpack on clusters. NVIDIA Corporation, 2010.

[36] Akira Nukada and Satoshi Matsuoka. Auto-tuning 3-d fft library for cuda gpus. In *Proceedings of the Conference on High Performance Computing Networking, Storage and Analysis*, SC '09, pages 30:1–30:10. ACM, 2009.

[37] Rice University. Dept. of Computer Science, D. Callahan, and K. Kennedy. *Compiling Programs for Distributed-memory Multiprocessors*. Technical report. Rice University, 1988.

[38] Josep M Perez, Rosa M Badia, and Jesus Labarta. A dependency-aware task-based programming environment for multi-core architectures. In *Cluster Computing, 2008 IEEE International Conference on*, pages 142–151. IEEE, 2008.

[39] Josep M Perez, Pieter Bellens, Rosa M Badia, and Jesus Labarta. CellSs: Making it easier to program the cell broadband engine processor. *IBM Journal of Research and Development*, 51(5):593–604, 2007.

[40] Laércio L Pilla, Christiane Pousa Ribeiro, Daniel Cordeiro, and Jean-François Méhaut. Charm++ on numa platforms: the impact of smp optimizations and a numa-aware load balancer. In *4th workshop of the INRIA-Illinois Joint Laboratory on Petascale Computing. Urbana, IL, USA*, 2010.

[41] Judit Planas, Rosa M Badia, Eduard Ayguadé, and Jesus Labarta. Hierarchical task-based programming with StarSs. *International Journal of High Performance Computing Applications*, 23(3):284–299, 2009.

[42] A. Rogers and K. Pingali. Process decomposition through locality of reference. PLDI '89.

[43] Vladimir Subotić, Steffen Brinkmann, Vladimir Marjanović, Rosa M. Badia, Jose Gracia, Christoph Niethammer, Eduard Ayguade, Jesus Labarta, and Mateo Valero. Programmability and portability for exascale: Top down programming methodology and tools with StarSs. *Journal of Computational Science*, (0):–, 2013.
[44] Enric Tejedor, Montse Farreras, David Grove, Rosa M Badia, Gheorghe Almasi, and Jesus Labarta. ClusterSs: a task-based programming model for clusters. In *Proceedings of the 20th international symposium on High performance distributed computing*, pages 267–268. ACM, 2011.

[45] Stanimire Tomov, Jack Dongarra, and Marc Baboulin. Towards dense linear algebra for hybrid gpu accelerated manycore systems. *Parallel Comput.*, pages 232–240, June 2010.

[46] Haluk Topcuoglu, Salim Hariri, and Min-you Wu. Performance-effective and low-complexity task scheduling for heterogeneous computing. *Parallel and Distributed Systems, IEEE Transactions on*, 13(3):260–274, 2002.

[47] Daouda Traore. Nanos++ User Manual. https://pm.bsc.es/projects/nanox.

[48] Clint Whaley, Antoine Petitet, and Jack J. Dongarra. Automated empirical optimization of software and the atlas project. *PARALLEL COMPUTING*, 27:2001, 2000.

[49] A. YarKhan, J. Kurzak, and J. Dongarra. Quark users' guide: Queueing and runtime for kernels. Technical report, Innovative Computing Laboratory, University of Tennessee, 2011.