# Scalable and portable LU factorization with partial pivoting on top of runtime systems

Alycia Lisito, Mathieu Faverge, Matthieu Kuhn, Florent Pruvost, Pierre Ramet

## ► To cite this version:

Alycia Lisito, Mathieu Faverge, Matthieu Kuhn, Florent Pruvost, Pierre Ramet. Scalable and portable LU factorization with partial pivoting on top of runtime systems. IPDPS25 - 39th IEEE International Parallel and Distributed Processing Symposium, Jun 2025, Milan, Italy. hal-04984070

HAL Id: hal-04984070

https://inria.hal.science/hal-04984070v1

Submitted on 10 Mar 2025

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published 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
# Scalable and portable LU factorization with partial pivoting on top of runtime systems

Alycia Lisito*†, Mathieu Faverge†, Matthieu Kuhn*, Florent Pruvost†, Pierre Ramet†

*Eviden

Pessac, France

†Univ. Bordeaux, CNRS, Bordeaux INP, INRIA, LaBRI, UMR 5800, F-33400 Talence, France

**Abstract**—Task-based runtime systems have demonstrated efficiency in leveraging the capabilities of large, heterogeneous architectures. Many linear algebra algorithms and applications have been implemented on top of runtime systems to increase their performance. However, the High Performance Linpack (HPL) benchmark, used by the TOP500 to rank supercomputers, has not yet been successfully implemented using task-based runtime systems. In this paper, we explore solutions to implement efficient LU factorization with partial pivoting using the sequential task-flow programming model. We show that, due to the pivoting strategy, this algorithm generates a large number of very small tasks, which usually overload the runtime system and make it inefficient. We propose two solutions to improve the efficiency and reduce the number of tasks. First, we apply well-known blocking strategies in the context of task-based algorithms. Secondly, we explore batching techniques to reduce the number of tasks submitted to the runtime system. Moreover, in distributed architectures, partial pivoting generates many reductions on the critical path throughout the factorization which needs to be carefully handled to reach high performance. Two task-based reduction algorithms are proposed to express these operations and improve the runtime reactivity on the critical path. These proposals have been implemented in the dense linear algebra library CHAMELEON on top of the STARPU runtime system. Experiments conducted on our cluster with these optimizations show that our LU with partial pivoting asymptotically reaches the performance of the non-pivoting algorithm.

***Index Terms***—LU Partial pivoting, runtime system, task-based application, HPL

## I. INTRODUCTION

Nowadays, modern machines embed several computing hardware solutions. Indeed, nine out of the ten fastest supercomputers in the world include nodes that are equipped with both CPUs and GPUs [1], [2], this heterogeneity results in a higher implementation cost [3]. Efficiently exploiting all the available computing resources requires implementing several levels of parallelism and each material comes with a large number of software solutions. This presents challenges for performance portability and code maintainability in hybrid parallel systems. Over the past two decades, task-based paradigms have frequently been shown as an efficient solution to address modern computing architectures as they provide a good level of performance portability without increasing the code complexity. Therefore, several linear algebra libraries implementing task parallelism have emerged [3]–[8] and more complex applications have been built on top of these [9]–[11].

Although numerous dense linear algebra algorithms have been implemented and well-studied libraries exist, the LU factorization with partial pivoting remains one of the few algorithms that does not have an efficient task-based implementation for heterogeneous systems. The distinctive nature of this algorithm, which integrates both large-grain efficient tasks and global small-grain reductions, renders it less compatible with task-based programming models. Nevertheless, the LU factorization remains the primary component of the High-Performance Linpack (HPL) benchmark [12], which is used to rank the fastest supercomputers in the world through the Top500 [1]. There is a lack of implementations of the LU factorization with partial pivoting available on distributed heterogeneous hardware, and even fewer efficient implementations. However, there are several highly optimized implementations of the HPL benchmark usually provided by the constructors to benchmark their hardware. These implementations are carefully tailored for each hardware. To the best of our knowledge, none of these solutions have been implemented on task-based runtime systems to provide performance portability on various architecture hardware. They usually follow a regular fork-join computation pattern that is carefully tuned with highly optimized kernels on the fastest available hardware, e.g. the GPUs. Additionally, these implementations are not heterogeneous, in the sense that they are not able to efficiently and fully exploit both CPUs' and GPUs' computational power at the same time.

In this paper, our goal is to use the task-based programming model to propose a portable and efficient implementation of the LU factorization with partial pivoting (GETRF) on distributed CPU systems. Our final objective is to be able to efficiently exploit heterogeneous nodes with this task-based GETRF in order to extend it to a task-based HPL benchmark. Based on existing task-based programming models and the various works that have been proposed using runtime systems, we study multiple ways to express the LU factorization with partial pivoting. The GETRF raises both problems of granularity and low arithmetic intensity because of the partial pivoting. We propose to increase the computational intensity and to modify the granularity of the panel operations by coupling an inner blocking strategy to a task batching mechanism as our first contribution. Thanks to these techniques, we can reduce the number of tasks submitted to the runtime and
adjust their size which minimizes the runtime overhead and
allows us to exploit the runtime efficiently. The other challenge
resides in the large number of reduction operations that
generate synchronization barriers throughout the algorithm.
The second contribution is various task-based implementations
of these reductions using the *Cyclic Shift* algorithm for the
communications ordering through two strategies. One uses the
runtime's tasks to do the MPI communications and the other
does the MPI communications in one task to perform the
**ALLREDUCE**. These strategies greatly improve the reactivity
of the reductions that are on the critical path and are essential
to increase the performances on distributed architectures. Our
last contribution is to propose a solution to control the task
flow of the algorithm to avoid overflowing the runtime and to
regulate the memory usage in distributed architecture. We used
a lookahead in order to control the number of tasks in the sys-
tem and temporary workspaces to reduce the number of ready
receptions at each instant of the factorization and thus reduce
the memory overhead. Handling those challenges is crucial
to being able to exploit the heterogeneity of the machines
efficiently with a runtime system, which is our long-term
goal. This work is made available in the **CHAMELEON** [3],
[13] dense linear algebra library on top of the **STARPU** [14]
runtime system. While all of these are known techniques, they
are applied here for the first time in the context of task-based
programming models to achieve a portable and efficient LU
factorization with partial pivoting.

Section II presents related work on dense linear algebra on
top of runtime systems. Section III introduces the LU factor-
ization without partial pivoting with task-based programming
models. Section IV describes a naive approach to implement
a task-based LU factorization with partial pivoting. Section V
details our contributions to improve on this starting point, and
Section VI presents the performance study of the propositions.

## II. RELATED WORK

Many dense linear algebra libraries have recently been
developed on top of runtime systems to address modern
architectures. These libraries rely on expressing algorithms as
a set of operations applied on tiles, or small blocks of the
matrix, and different ways to express the direct acyclic graph
of dependencies where edges are data dependencies, and nodes
are elementary kernels. This principle has been introduced
in [4] in the PLASMA [15] library. This library expresses
the algorithms thanks to the sequential task-flow paradigm
on top of the QUARK [16] runtime system. In this paradigm,
tasks are asynchronously submitted to a runtime system and
data used by the tasks are annotated with access keywords:
Input, Output, Inout. The runtime system automatically infers
from these annotations the Directed Acrylic Graph (DAG)
of the application and schedules the tasks accordingly. The
LU factorization has been widely studied in this library and
multiple variants, including the classic partial pivoting, have
been proposed [17]. All these implementations, however, are
targeting single homogeneous nodes.

Thanks to GPU kernels proposed by the MAGMA [15]
library, the LU factorization with partial pivoting has been
extended to heterogeneous systems in both the PLASMA [18]
and the CHAMELEON [19] libraries. Still, these solutions are
for non-distributed computations. The CHAMELEON library
supports distributed heterogeneous systems for Cholesky and
QR factorizations [13], as well as for LU factorization with-
out pivoting. Similarly, PLASMA has been converted to the
SLATE [8], [20] library to exploits the OPENMP and MPI
standards on distributed architectures. However, this library
does not provide a LU factorization with partial pivoting
despite a tournament pivoting implementation [21].

The DPLASMA [22] library implemented on top of the
PARSEC [23] runtime system is another alternative solution to
SLATE and CHAMELEON. The algorithms are implemented
using the parameterized task-graph (PTG) programming model
specific to PARSEC and introduced in [24]. This library is
one of the few that proposes an LU factorization with partial
pivoting (GETRF_PTGPANEL in the code) while its imple-
mentation has seemingly never been studied in a paper. This
implementation describes the algorithm with fine grain tasks
to perform the pivoting strategy and the swap of the submatrix
update. The GETRF implemented in DPLASMA only uses
inner blocking technique to improve the task granularity and
lookahead to regulate the task flow. The reductions are done
through tasks with PARSEC using the *Cyclic Shift* algorithm
for the communication sequence. This work is the closest to
what we propose in this paper, although the task batching
strategy and MPI communications within one task for the
reductions are not used in the DPLASMA implementation of
GETRF.

CHARM++ is another programming model using an object-
based paradigm and is message-driven [25]. In this model,
the programmer application is over-decomposed to express all
parallelism. The authors proposed a scalable implementation
of the LU factorization on top of this model that reached 67%
efficiency on a quarter of the former Jaguar supercomputer.
However, this implementation did not combine tasks to reduce
the runtime overhead and did not optimize reduction opera-
tions.

Focusing only on the LU factorization, many HPL imple-
mentations have been developed by various constructors to
validate the performance of their clusters such as :intel
HPL, AMD CPU HPL, AMD GPU rocHPL [26], nvidia HPL.
All of these implementations include optimizations such as
inner blocking, lookahead strategies, and eventually dedicated
broadcast algorithms, but none of them are runtime based
which makes them not easily portable to other architectures.
Note also that the HPL algorithm simultaneously performs the
LU factorization and the solve of one single right-hand side.
Therefore, the permutation does not need to be applied on the
full matrix A, and the first triangular solve can be combined
with the matrix factorization.

Batching techniques introduced with the Batch BLAS ker-
nels [27] API, or replication strategies can be employed to
optimize the algorithm such as the TRSM replication used
in [28].

In this paper, we explore solutions to express these optimizations in the task-based programming model in order to reach high efficiency on various CPU architectures.

### III. BACKGROUND

In order to solve the equation $Ax = b$, where the matrix $A$ is general, a common strategy is to factorize $A$ into the product of a lower and an upper triangular matrix $L$ and $U$, with an operation called `GETRF` (General TRIangular Factorization). There are different ways to factorize $A$, we will first talk about the blocked factorization without pivoting to simplify the algorithm (`GETRF_NOPIV`). This factorization without pivoting provides us a baseline performance as it generates as many flops as the factorization with partial pivoting. The tiled `GETRF_NOPIV` algorithm consists of three nested loops as shown by Algorithm 1. At each iteration $k$, algebraic operations are applied to the blocks, also called tiles ($k$ iterating on the diagonal blocks and $nt$ (resp. $mt$) being the number of blocks per row (resp. column) of the matrix). First, a `GETRF_NOPIV` factorization is performed on the diagonal block $A_{k,k}$. Then, the blocks in the $k^{th}$ row, respectively column, are updated by solving a triangular system (`TRSM`) using the resulting upper part (`TRSM_RU`, $XU_{k,k} = A_{i,k}$), respectively lower part (`TRSM_LL`, $L_{k,k}X = A_{k,j}$), of the diagonal block factorization. Finally, the trailing submatrix tiles are updated with general matrix-matrix multiplications (`GEMM`, $A_{i,j} = A_{i,j} - A_{i,k} * A_{k,j}$).

**Data:** $A$

**Result:** $A \leftarrow L \times U$

for $k = 0$ to $\min(mt, nt) - 1$ do

$A_{k,k} \leftarrow \text{GETRF}(A_{k,k})$

for $m = k + 1$ to $mt - 1$ do

$A_{m,k} \leftarrow \text{TRSM_RU}(A_{k,k}, A_{m,k})$

end

for $n = k + 1$ to $nt - 1$ do

$A_{k,n} \leftarrow \text{TRSM_LL}(A_{k,k}, A_{k,n})$

for $m = k + 1$ to $mt - 1$ do

$A_{m,n} \leftarrow \text{GEMM}(A_{m,k}, A_{k,n}, A_{m,n})$

end

end

end

**Algorithm 1:** Tiled algorithm of the `GETRF_NOPIV` factorization on a matrix $A$ composed of $mt \times nt$ tiles.

Once the algorithm is expressed through tiled operations such as seen in Algorithm 1, it can be exploited through various task-based programming models where each block operation is seen as a task. The tasks can be described through a specific language with all its incoming and outgoing dependencies as done in PARSEC [23]. Or, more commonly, additional information are added to the sequential algorithm with pragmas or extra parameters to express the data accesses of the tasks. The runtime (OPENMP [29], STARPU [14], Legion [30]...) automatically infers the dependencies and schedule the tasks accordingly. Figure 1 shows such a generated DAG for a matrix of $4 \times 4$ tiles.

<figure><img src="image_2.png" /><figcaption>Fig. 1: DAG generated by the `GETRF_NOPIV` algorithm applied to a matrix with $nt = mt = 4$.</figcaption></figure>

This representation helps improve the portability of the application because it separates many concerns that have to be addressed on large architectures. Indeed the algorithm becomes oblivious of the architectures it runs on as opposed to tuned implementation of HPL for example. Having such a representation may help test various mapping strategies without having to consider rewriting the communications. And as opposed to regular fork-join implementations, scheduling strategies can help to overlap computation and communication while deciding the best available architecture to run on. The complexity of handling the hardware is delegated to the runtime system and the application programmer focuses on expressing the algorithm in the most efficient way to generate parallelism.

Different strategies exist to schedule the execution of the tasks of this graph. The best scheduling strategy depends on the computer's architecture and the chosen optimization criteria (for instance, the memory or the energy consumption). In the graph in Figure 1, we can see that the `GEMM` (green) task (2,3) of the iteration $k = 1$ depends on the `TRSM` (blue) tasks (1,3) and (2,1) of the same iteration but also on the `GEMM` (green) task (2,3) of the previous iteration $k = 0$. This means that iteration 1 can start before iteration 0 finishes (depending on the chosen scheduling strategy), which results in automatic look-ahead. Nevertheless, handling these tasks dependencies is complicated on heterogeneous and / or distributed architectures. They might induce data transfer from a node to another, or from CPU to GPU. Having the runtime system managing them transparently is a relief for the application programmer.

For this reason, we focus our work on expressing the `GETRF` algorithm on top of runtime systems to explore the limits of such a programming model and its potential to offer efficient implementations of such algorithms mixing fine and
large grain tasks. The LU factorization with partial pivoting is implemented in the CHAMELEON library that covers already a large spectrum of dense linear algebra subroutines such as BLAS operations and one-sided factorizations. CHAMELEON supports multiple runtime systems, but we focus in this work on the STARPU runtime system. Indeed, it provides several scheduling strategies that are relevant to our algorithm, as well as support for heterogeneous and distributed architectures.

## IV. GETRF WITH PARTIAL PIVOTING

In the scalar `GETRF_NOPIV` operation, every diagonal element is used to normalize the corresponding columns, which can lead to numerical errors if one of these elements is close to zero. In order to preserve the numerical stability, pivoting strategies have been developed and the most well-known and generally used is the partial pivoting strategy. By searching in each column for the greatest element in absolute value as a pivot, this algorithm ensures a practical numerical stability of the algorithm and avoids divide-by-0 issues unless the matrix presents some pathological properties. Thus, the final solution can be written as $PA = LU$ where $P$ is a permutation matrix. We now focus the study on the panel factorization, before discussing the update of the trailing submatrix.

a) Panel factorization: Partial pivoting generates a reduction operation at each column of the matrix in order to find the best pivot in each of them before pursuing the factorization. The `GETRF_NOPIV` and `TRSM_RU` operations of the same iteration $k$ in Algorithm 1 are therefore replaced by Algorithm 2, this is referenced as the panel factorization. In this paper, $mt$ is the number of tiles in the column of the matrix $A$, and $nb$ is the size of these tiles (number of rows and columns). If we consider each operation as a task, for each column, we first have a set of at most $mt$ tasks to search for the maximum element of the column of the block which generates at most $mt$ potential pivots. These tasks are involved in the reduction operation `REDUX_PIVOT` to compute the pivot of the column. In shared memory, runtime systems such as STARPU offer the `REDUCE` keyword to specify that a data is involved in a reduction and thus optimize this operation potentially creating a few more tasks. In distributed memory, communication tasks may be automatically inserted into the graph to operate this reduction. We will detail this operation in subsection V-B.

Once the pivot is computed, the `SWAP_ROWS` operation, which permutes the row of the pivot with the diagonal row, has to be performed. However, the task submission thread is not aware of the selected pivot, thus a `SWAP_ROWS` task has to be submitted for each tile of the panel even if only two tiles are involved. This ensure that the data dependencies are respected for the tiles that exchange rows. These are the orange tasks in Figure 2 and Algorithm 2.

Finally, the main column can be divided by the selected pivot and the remaining part of each tile can be updated. These operations are represented in pink in Figure 2 and Algorithm 2. The `UPDATE` tasks are the only tasks of the panel factorization that contain computations through Level 2 BLAS operations.

**Data:** $k, A_k: nt, k$

**Result:** $A_k: nt, k \leftarrow P_k \times L_k: nt, k \times U_k, k$

$P_k \leftarrow \text{Id}$

**for** $c = 1$ **to** $nb-1$ **do**

**for** $m = k$ **to** $mt-1$ **do**

$p(k) \leftarrow \text{REDUX_PIVOT}$ `SEARCH_MAX}(A_{m,k}, c)$

**end**

**for** $m = k$ **to** $mt-1$ **do**

$[A_{k,k}, A_{m,k}] \leftarrow \text{SWAP_ROWS}(A_{k,k}, A_{m,k}, c, p(k))$

**end**

**for** $m = k$ **to** $mt-1$ **do**

$A_{m,k} \leftarrow \text{UPDATE}(A_{m,k}, A_{k,k}, c)$

**end**

**end**

$P_k \leftarrow \text{IPIV_TO_PERM}(p())$

**Algorithm 2:** Factorization of the panel $k$ of the `GETRF` algorithm.

<figure><img src="image_3.png" /><figcaption>Fig. 2: Panel factorization with and without partial pivoting.</figcaption></figure>

The task-based implementation of the panel to cover all the possible pivot selection generates a huge number of fine grain tasks with low arithmetic intensity that are difficult to schedule efficiently. In this first approach, we indeed generate $3 \times nb \times (mt - k)$ tasks for the $k^{th}$ panel factorization, while the non-pivoting algorithm only generates $mt - k$ tasks. This number can first be reduced to $(nb+1) \times (mt - k)$ by merging together the `SWAP_ROWS` and `UPDATE` tasks of the step $c$ with the `REDUX_PIVOT` tasks of step $c+1$ steps as shown on Figure 2. We consider this solution as our reference task-based implementation and discuss optimization in section V. As shown by Table I, this solution of the panel factorization with partial pivoting generates $nb+1$ times the number of tasks of the panel factorization without pivoting.

<table><caption>TABLE I: Number of tasks generated during the factorization without pivoting and with partial pivoting considering $mt = nt$</caption><thead><tr><th></th><th>Panel factorization</th><th>Submatrix update</th></tr></thead><tbody><tr><td>No Pivoting</td><td>$\frac{1}{2}mt(mt+1)$</td><td>$\frac{1}{3}(mt-1)mt(mt+1)$</td></tr><tr><td>Partial Pivoting</td><td>$\frac{1}{2}mt(mt+1)(nb+1)$</td><td>$(mt-1)mt(mt+1)$</td></tr></tbody></table>
b) Submatrix update: Similarly to the panel factorization, describing the submatrix update of the GETRF algorithm with a task-based approach requires describing the permutation of the rows induced by the pivots selection in the panel. Algorithm 3 lists the tasks submitted to perform both the permutation of the rows and the resulting update. To ease the row exchange, the pivot array is converted to a permutation and inverse permutation array in the final task of the panel factorization (See IPV_TO_PERM Algorithm 2). Thus, the inverse permutation array contains the indices of the final rows of $U_{k,k:nt}$ and the permutation array contains the destination indices of each row in $A_{k,k:nt}$. As for the panel, by default, the task submission process has to cover all eventualities by submitting a collecting task (SWAP_GET) and a setting task (SWAP_SET) for each tile of the column. Both operations are separated in order to give priorities to the collecting tasks which are critical to the computations of the $U_{k,n}$ while the setting tasks are less critical and could be merged with the GEMM updates.

The number of tasks similarly grows in each panel from $mt-k$ tasks to $3 \times (mt-k)$ (see Table I). Note that we create a copy of $A_{k,k:nt}$ in order to remove the write after read dependency existing between the reduction operation (SWAP_GET) and the SWAP_SET tasks.

**Data:** $k, A_{k:mt,k:nt}, P_k$
**Result:** $A_{k:mt,k:nt} \leftarrow L_{k:mt,k:nt} \times U_{k:mt,k:nt}$
**for** $n = k + 1$ **to** $nt - 1$ **do**
> Acopy$_{k,n} \leftarrow COPY(A_{k,n})$
> $A_{k,n} \leftarrow SWAP_GET(A_{k,n}, P_k)$
> **for** $m = k + 1$ **to** $mt - 1$ **do**
> > $A_{k,n} \leftarrow REDUX_PERM$ SWAP_GET($A_{m,n}, P_k$)
> > $A_{m,n} \leftarrow SWAP_SET(Acopy_{k,n}, P_k)$
> **end**
> $A_{k,n} \leftarrow TRSMSouth(A_{k,n}, A_{k,n})$
> **for** $m = k + 1$ **to** $mt - 1$ **do**
> > $A_{m,n} \leftarrow GEMM(A_{m,k}, A_{k,n}, A_{m,n})$
> **end**
**end**

**Algorithm 3:** Permutation and update of the trailing submatrix at iteration $k$ of the GETRF algorithm.

## V. CONTRIBUTIONS

The naive task-based implementation of the GETRF algorithm proposed in the previous section generates multiple issues that can be tricky for a runtime system to handle correctly. For example, the panel factorization in the no-pivoting variant block-wise becomes then block-wise × column-wise, which generates two main issues: the small granularity of the tasks and their large number. In GETRF_NOIV, for a panel of size $m \times n$, the operation count is $\frac{2}{3}mn^2$ for $t$ tasks (with $t$ the number of tiles in the column), while in GETRF, the number of operations remains the same for $(n+1) \times t$ tasks. The number of operations per tasks is thus divided by the tile size. This significantly increases the risk of performance loss because of the runtime overhead.

According to the STARPU recommendations, the execution time of a task should lie between 1ms and 10ms on the

targeted architecture in order to maximize its efficiency and minimize the runtime task overhead [31]. In the naive algorithm, for typical block sizes, the average execution time of a panel task is 0.1ms which is too small for STARPU and slows down the execution of our LU factorization.

Moreover, the REDUX_PIVOT slows down the algorithm at each column iteration as every worker has to wait for each other and the next column iteration cannot start before this operation is done. The scheduling and reactivity on these tasks are essential.

During the trailing submatrix update, the REDUX_PERM are not as critical and can be performed in parallel. However, their number is so important that it can overload the communication thread of the runtime system and slow down the computations.

We now present a set of optimizations to improve this original task-based implementation.

### A. Blocking and batching the panel's tasks

It has been known for a long time that Level 3 BLAS offers better arithmetic intensity than Level 2 BLAS and it is important to block operations within dense linear algebra kernels to get good performances. The first proposition to improve the panel's tasks granularity is to use inner blocking as in the LAPACK library [32]. Inner blocking introduces a new level of block division *ib*: each block of the panel is split into tasks of size $nb \times ib$. This allows us to limit the use of Level 2 BLAS operations on each sub-block of the panel and use Level 3 BLAS operations to update the rest of the panel, hence improving the arithmetical intensity of a few panel tasks. Note that the number of tasks remains the same but some operations are grouped together for better efficiency. This is illustrated by the surface of the Level 2 BLAS operations (UPDATE in pink) on Figure 3 which is larger on the left without inner blocking than on the right with inner blocking. Every *ib* column, updates are replaced by Level 3 BLAS illustrated in blue and green on the right of the figure (respectively TRSM South and GEMM). Since the tasks always depend on the full tile, this change only requires modifying the kernel that is submitted.

Once the kernel efficiency has been increased, the large number of tasks still remains an issue. To reduce it, we exploit a batching solution to group small tasks into larger ones. To implement this, we use the capability of STARPU to submit tasks that depend on an array of data. This allows us to group several block dependencies together, and to submit one task for several blocks, applying the same calls to the same kernels in a sequential loop instead of submitting one task per block in the panel. STARPU also proposes an automatic batching strategy [33]. But due to the criticality of the tasks on the panel, we prefer to manually decide which tasks should be grouped together. We will name the batch size *b* as the number of tasks grouped together for the rest of this paper. The batching technique divides the number of tasks generated by *b*, leading to *b* times larger tasks, which reduces the runtime overhead. Figure 4 compares the tasks generated during the first panel factorization without batching (Figure 4a) and with
<figure><img src="image_4.png" /><figcaption>Fig. 3: Example of a panel factorization on a 27 × 27 matrix with blocks of size $nb = 9$ and inner blocks of size $ib = 3$.</figcaption></figure>

a batching (Figure 4b) of size $b = 3$ on top of the inner blocking strategy. The algorithm without batching generates 77 tasks (black squares Figure 4a) and the algorithm with a batching of size 3 generates 28 tasks (black frames in Figure 4b).

<figure><img src="image_5.png" /><figcaption>Fig. 4: Example of a panel factorization on a matrix with 11 × 11 blocks of size 6 × 6.</figcaption></figure>

An analysis of our batching approach is conducted on a square matrix of size $N = 115200$ with $nb = 576$, $ib = 32$ and $b = 12$ and presented in Table II. The parameters ($nb$, $b$ and $ib$) are issued from the tuning done in Section VI. This analysis compares the number of tasks, their size and granularity generated on a full-matrix factorization w/o and w/ inner blocking, as well as w/o and w/ batching.

Focusing first on Figure 3, we observe that the inner blocking strategy allows for a few tasks to reach a higher arithmetic intensity without changing the number of tasks submitted to the runtime system. In Table II, this reflects in a higher overall performance with and without batching.

In addition, independently of the inner blocking strategy, Table II shows that the batching mechanism reduces the number of tasks as expected by a ratio slightly lower than the batch size (11.37). It also increases the average time spent

in each task to reach a level closer to the recommendations of STARPU.

Figure 5 confirms that both solutions correctly group the tasks to reach a better suited number of operations to comply with the recommendation. The 1ms bar shows the theoretical amount of flops required in a task to be computed in 1ms. Since the Level 2 BLAS tasks are memory bound, the time actually spent is larger than the one estimated through the number of flops. Despite this observation, combining the inner blocking and the batching mechanism offers the best possible performances, showing a speedup of more than 6 over the naive approach. Note that the diagonal tasks are not batched for the solution without inner blocking due to different data access between the diagonal and offdiagonal tasks as well as the TRSM operations for the inner blocking version. This explains the difference in the number of tasks between the two versions with the batching and without the batching in Table II.

<table><caption>TABLE II: Panel task statistics of a matrix of size 115200 with a block size of 576</caption><thead><tr><th>inner-blocking<br>batching b = 12</th><th>no<br>no</th><th>no<br>yes</th><th>yes<br>no</th><th>yes<br>yes</th></tr></thead><tbody><tr><th scope="row">Tasks nbr</th><td>11 597 700</td><td>1 125 727</td><td>11 601 100</td><td>1 023 536</td></tr><tr><th scope="row">Time (ms/task)</th><td>0.49971</td><td>3.52709</td><td>0.50559</td><td>3.90406</td></tr><tr><th scope="row">GFlops/s</th><td>1.065241</td><td>1.179923</td><td>2.430522</td><td>6.558846</td></tr></tbody></table>

<figure><img src="image_6.png" /><figcaption>Fig. 5: Task size repartition of the Ppanel factorization of a matrix of size 115200 with a block size of 576</figcaption></figure>

## B. Optimizing reactivity on AllReduce operations

The REDUX_PIVOT operations performed during the panel factorization lie on the critical path of the GETRF. Indeed, SWAP_ROWS and UPDATE operations in the panel can only be performed once the previous REDUX_PIVOT operations is over. This operation is performed at two levels: locally inside each MPI process, and then, globally between the MPI processes. Inside each MPI process, the reduction is efficiently managed by the runtime system. However, between the MPI processes, the runtime system does not necessarily handle the reduction efficiently and it can be more efficient to perform the reduction at the application level (in CHAMELEON) instead. Different strategies can be explored to perform the reduction between the MPI processes at the application level.
When STARPU handles the reduction at the distributed level, through the insertion of additional tasks, it can automatically start the distributed reduction process as soon as the locally reduced data is computed. However, CHAMELEON features an automatic pruning mechanism that conflicts with the reduction algorithm from STARPU which can no longer correctly detect all the MPI processes involved in the reduction. This mechanism ensures that a task is not submitted to the runtime if the current MPI process is not involved. Therefore, removing this mechanism would induce another task overhead. Secondly, the batching mechanism makes it difficult to correctly submit the required tasks to enable the fully automatic detection of MPI communications by STARPU. A memory allocation per batched task is required in order to submit the array of blocks batched to the runtime at the task insertion, this array is freed by the runtime at the end of the task execution. In order to enable the fully automatic detection of MPI communications by STARPU, we would have to submit batched tasks that will not be executed. The difficulty is to free the memory allocated for these arrays as it will not be done by the runtime. We have to wait until the task has been handled by the runtime which adds complexity to our algorithm and greatly increase the memory usage. Therefore, we investigate multiple distributed ALLREDUCE implementations (for both the REDUX_PIVOT and REDUX_PERM) made at a higher level, in CHAMELEON, relying on STARPU only for the preliminary local reductions.

There exists many ALLREDUCE algorithms in literature, among them we chose the Cyclic Shift algorithm [34], known to be one of the most efficient ALLREDUCE algorithms. When *p* processes are involved, the Cyclic Shift algorithm limits the number of MPI communications to one round of *log₂(*p*)* instead of two rounds of *log₂(*p*)* communications for the naive reduction and communications for the broadcast.

**Data:** $p_{panel}$, myrank
$me = where_i_am(p_{panel})$

$np = size_of(p_{panel})$

$np_iter = np$

$shift = 1$

while $np_iter > 1$ do

$p_{send} = p_{panel}[(me + shift)\%np]$

$p_{recv} = p_{panel}[(me - shift)\%np]$

MPI_Isend(data_send, myrank, p_send)

MPI_Irecv(data_recv, myrank, p_recv)

my_redux(data_send, data_recv, myrank, p_recv)

$shift = shift << 1$

$np_iter = ceil(np_iter/2)$

end

**Algorithm 4:** ALLREDUCE with the Cyclic Shift algorithm.

To implement the Cyclic Shift algorithm, we first store the IDs of the processes involved in the panel when submitting the tasks in an array called *ppanel*. As shown in Algorithm 4, the current MPI process *myrank* = *ppanel*[me] sends

(resp. receives) data to (resp. from) the MPI process $p_{send} = p_{panel}[me + shift]$ (resp. $p_{recv} = p_{panel}[me - shift]$). The reduction between the local data $data_{send}$ and the distant data $data_{recv}$ is done before the next iteration in order to send the result of this reduction to the next process.

Figure 6 depicts how the data *piv*₀ is propagated to the other processes during the ALLREDUCE form the point of view of the process 0. The data is propagated to $2^{i-1}$ processes at the iteration *i*. For instance before the iteration 4 all the processes from 0 to 7 have their data reduced with the data *piv*₀.

<figure><img src="image_7.png" /><figcaption>Fig. 6: Example of the ALLREDUCE iterations with the Cyclic Shift algorithm on 25 MPI processes from the point of view of the process rank 0 for the reduction of the pivot.</figcaption></figure>

We first propose a task-based implementation of ALLREDUCE operations letting STARPU infer and handle the MPI communications. Figure 7 presents this solution. The submission thread goes through the loop of the Cyclic Shift algorithm to submit: one task to send the local data (blue rectangles), and one task to receive the data (purple rectangles) which also performs the reduction operation (mauve diamond). The communication thread handles the MPI_Isend and MPI_Irecv communications and once the reception communication is completed, the reduction operation is executed on an execution thread.

<figure><img src="image_8.png" /><figcaption>Fig. 7: ALLREDUCE version with STARPU tasks.</figcaption></figure>

This solution remains very close to the task-based programming model: we let STARPU handle all data exchanges and schedule the communications and computations to perform the ALLREDUCE. In fact, the different steps of the ALLREDUCE can be performed asynchronously between the MPI processes involved without blocking any execution thread. However, this again generates many small tasks to perform a critical
operation, and the reactivity of these tasks may be impacted
by the fact that communications and computations are not
performed by the same worker.

A second implementation we propose is to reduce the
number of tasks during the ALLREDUCE to a single task.
This solution follows the idea proposed by the Task-Aware
MPI (TAMPI) library [35]. The authors proposed an interface
to improve the interoperability between the runtime and the
MPI engine to enable efficient MPI communications within
a task without conflicting with the internal communications
of the runtime system. However, to follow this new API
some changes are required at the runtime level. STARPU
is not TAMPI compatible yet, so we chose to use a dif-
ferent MPI communicator to submit the communications of
the ALLREDUCE. In order to perform more than one MPI
communication at the same time, we have to activate the
MPI_THREAD_MULTIPLE support of MPI. Figure 8 presents
this solution. A single task (in green) is submitted to the
runtime system, and when a worker starts working on it, it
will perform both the communication and the reduction. This
helps to greatly reduce the latency of the kernel as everything
is performed by the same worker.

<figure><img src="image_9.png" /><figcaption>Fig. 8: ALLREDUCE version with the MPI communication in
one STARPU task.</figcaption></figure>

This second solution fixes the task overhead problem but
it introduces scheduling problems. Indeed, if one execution
thread $t_0$ on the MPI process $p_0$ starts one ALLREDUCE,
it will have to wait for the execution threads of the MPI
processes $p_1$ to $p_{np}$ to start that same ALLREDUCE. As a
result, $t_0$ cannot start another task while waiting for the MPI
communications, which results in a waste of resources. This
problem is mitigated on the panel for the REDUX_PIVOT
operation because this operation is highly critical and its
priority almost guarantees that all MPI processes involved will
try to work on it as soon as possible. However, when the same
algorithm is applied on the trailing submatrix update with the
REDUX_PERM operation, the synchronization problem may
arise much more often as we have no guarantee that all MPI
processes will perform the updates at similar times and in the
same order despite similar priorities.

### C. Task flow regulation

In [13], the authors introduced a mechanism in STARPU
to limit the number of submitted tasks by doing automatic
pruning at the runtime level, and/or at the application level.

This aims to reduce the number of tasks on each MPI
processes by considering only the ones they are involved
with. We inherit from this feature in a distributed environment
by using the CHAMELEON library with STARPU. They also
propose to introduce task quotas by setting the minimum and
maximum number of tasks in the systems through STARPU
environment variables. This reduces the number of submitted
communications jointly with the memory overhead of the
reception buffers by preventing the runtime system from
allocating them too much in advance. In the GETRF algorithm,
setting such limits accurately is difficult since the number
of tasks is extremely large and varies significantly from one
iteration to the other (due to the *nb* factor in the panel tasks).

To solve this issue, we introduce temporary workspaces
which can hold up to *l* panels on each MPI processes,
where *l* is a look-ahead parameter. The broadcast of the
factorized panel is then replaced by copies into the temporary
workspaces, implemented by inserting STARPU distributed
copy tasks. These copies are made according to a ring topology
over the involved MPI processes.

Thus, the runtime system does not allocate reception buffers
in advance, and the write-after-read dependencies automati-
cally regulate the flow of tasks and the communications.

Note that the look-ahead parameter is used only for the
broadcast of the panel (*L* part). The broadcast of the row of
tiles (*U* part) that is required after the SWAP_GET reduction
does not require look-ahead as each iteration depends on the
previous one. However, we replaced the reduction algorithm
by an ALLREDUCE algorithm as in the panel in order to
generate a copy of the $A_{k,k:nt}$ tiles. Then, the TRSM LL
is replicated on all MPI processes that are involved in the
panel computation, similarly to what is proposed by [28].
This prevents the MPI processes from having to wait for the
TRSM LL to be performed on the MPI processes that own
the first row and for the results to be broadcast.

## VI. PERFORMANCE ANALYSIS

The experiments have been conducted on our local cluster
of the Plafrim<fnref n="1" /> platform. The Intel nodes are equipped with
two Intel Cascade Lake 18-cores running at 2.6 GHz and
192 GB of memory and the AMD nodes are equipped with
two AMD Zen4 Genoa 24-cores running at 2.5 GHz and 256
GB of memory. The environment (Intel MKL 2023 for the
BLAS kernels and gcc 11.4) is handled with guix [36],
STARPU 1.4.7, NEWMADELEINE 2024, MKL/SCALAPACK
2024 and DPLASMA (c8da9621) are compiled within this
environment. The scheduler *prio* of STARPU is used for the
experiments. The MPI library used is NEWMADELEINE [37].
The experiments on distributed architecture use one MPI
process per socket and one worker thread per core, thus 18
threads per MPI process on the Intel nodes and respectively
24 threads per MPI process on the AMD nodes. The data
mapping follows the classic 2D block-cyclic matrix distribu-
tion. The CHAMELEON implementation has been developed

<footnote n="1">¹https://www.plafrim.fr</footnote>
recently proposed in STARPU [38] to describe the factorization as a two levels algorithm and better control the amount of tasks submitted to the system. The second aspect that needs to be covered is the support of heterogeneous architectures with efficient reduction operations. Incidentally, the proposed solutions can be made generic and directly integrated into a runtime system.

## ACKNOWLEDGEMENT

Experiments presented in this paper were carried out using the PlaFRIM experimental testbed, supported by Inria, CNRS (LABRI and IMB), Université de Bordeaux, Bordeaux INP and Conseil Régional d'Aquitaine (see https://www.plafrim.fr).

## REFERENCES

[1] J. J. Dongarra, "Performance of various computers using standard linear equations software," Computer Technical Report Number CS - 89 - 85, vol. 20, 11 2002. [Online]. Available: http://www.netlib.org/benchmark/performance.ps

[2] "The 63rd edition of the top500," 6 24. [Online]. Available: https://top500.org/lists/top500/2024/06/

[3] E. Agullo, C. Augonnet, J. 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 Jade Edition, ser. Applications of GPU Computing Series, W.-m. W. Hwu, Ed. Boston: Morgan Kaufmann, 2012, pp. 473-484. [Online]. Available: https://www.sciencedirect.com/science/article/pii/B9780123859631000344

[4] A. Buttari, J. Langou, J. Kurzak, and J. J. Dongarra, "A class of parallel tiled linear algebra algorithms for multicore architectures," Parallel Comput., vol. 35, no. 1, pp. 38-53, Jan. 2009. [Online]. Available: http://dx.doi.org/10.1016/j.parco.2008.10.002

[5] E. Chan, F. G. Van Zee, P. Bientinesi, E. S. Quintana-Orti, G. Quintana-Orti, and R. Van de Geijn, "Supermatrix: a multithreaded runtime scheduling system for algorithms-by-blocks," in Proceedings of the 13th ACM SIGPLAN Symposium on Principles and practice of parallel programming. ACM, 2008, pp. 123-132.

[6] F. Song, A. YarKhan, and J. J. Dongarra, "Dynamic task scheduling for linear algebra algorithms on distributed-memory multicore systems," in Proceedings of the Conference on High Performance Computing Networking, Storage and Analysis, ser. SC '09. New York, NY, USA: ACM, 2009, pp. 19:1-19:11. [Online]. Available: http://doi.acm.org/10.1145/1654059.1654079

[7] S. Tomov, J. J. Dongarra, and M. Baboulin, "Towards dense linear algebra for hybrid gpu accelerated manycore systems," Parallel Computing, vol. 36, no. 5-6, pp. 232-240, Jun. 2010.

[8] M. Gates, J. Kurzak, A. Charara, A. YarKhan, and J. J. Dongarra, "Slate: design of a modern distributed and accelerated linear algebra library," in Proceedings of the International Conference for High Performance Computing, Networking, Storage and Analysis, 2019, pp. 1-18.

[9] J. M. C. Carpaye, J. Roman, and P. Brenner, "Design and analysis of a task-based parallelization over a runtime system of an explicit finite-volume cfd code with adaptive time stepping," Journal of Computational Science, vol. 28, pp. 439-454, 2018.

[10] S. Abdulah, H. Ltaief, Y. Sun, M. G. Genton, and D. E. Keyes, "Exageostat: a high performance unified software for geostatistics on manycore systems," IEEE Transactions on Parallel and Distributed Systems, vol. 29, no. 12, pp. 2771-2784, 2018.

[11] A. Charara, H. Ltaief, D. Gratadour, D. Keyes, A. Sevin, A. Abdelfatih, E. Gendron, C. Morel, and F. Vidal, "Pipelineing computational stages of the tomographic reconstructor for multi-object adaptive optics on a multi-gpu system," in SC'14: Proceedings of the International Conference for High Performance Computing, Networking, Storage and Analysis. IEEE, 2014, pp. 262-273.

[12] J. J. Dongarra, "Performance of various computers using standard linear equations software," ACM SIGARCH Computer Architecture News, vol. 20, no. 3, pp. 22-44, 1992.

[13] E. Agullo, O. Aumage, M. Faverge, N. Furmento, F. Pruvost, M. Sergent, and S. Thibault, "Achieving high performance on supercomputers with a sequential task-based programming model," IEEE Transactions on Parallel and Distributed Systems, 2017. [Online]. Available: https://inria.hal.science/hal-01618526

[14] C. Augonnet, S. Thibault, R. Namyst, and P.-A. Wacrenier, "StarPU: a unified platform for task scheduling on heterogeneous multicore architectures," Concurrency and Computation: Practice and Experience, vol. 23, no. 2, pp. 187-198, 2011. [Online]. Available: https://inria.hal.science/inria-00550877

[15] E. Agullo, J. Demmel, J. 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," Journal of Physics: Conference Series, vol. 180, no. 1, p. 012037, jul 2009. [Online]. Available: https://dx.doi.org/10.1088/1742-6596/180/1/012037

[16] A. YarKhan, J. Kurzak, and J. J. Dongarra, "QUARK users' guide: queuing and runtime for kernels, version 1.0," University of Tennessee Innovative Computing Laboratory, Knoxville, Tennessee 37996, technical report UT-ICL-11-02, April 2011.

[17] S. Donfack, J. J. Dongarra, M. Faverge, M. Gates, J. Kurzak, P. Luszczek, and I. Yamazaki, "A survey of recent developments in parallel implementations of Gaussian elimination," Concurrency and Computation: Practice and Experience, vol. 27, no. 5, pp. 1292-1309, May 2014. [Online]. Available: https://inria.hal.science/hal-00986948

[18] J. Kurzak, P. Luszczek, M. Faverge, and J. J. Dongarra, "LU factorization with partial pivoting for a multicore system with accelerators," IEEE Transactions on Parallel and Distributed Systems, vol. 99, no. PrePrints, p. 1, 2012. [Online]. Available: https://inria.hal.science/hal-00809657

[19] E. Agullo, C. Augonnet, J. J. Dongarra, M. Faverge, J. Langou, H. Ltaief, and S. Tomov, "LU factorization for accelerator-based systems," in 9th ACS/IEEE International Conference on Computer Systems and Applications (AICCSA 11), Sharm El-Sheikh, Egypt, Jun. 2011. [Online]. Available: https://inria.hal.science/hal-00654193

[20] A. YarKhan, J. Kurzak, P. Luszczek, and J. J. Dongarra, "Porting the plasma numerical library to the openmp standard," International Journal of Parallel Programming, vol. 45, pp. 612-633, 2017.

[21] R. Alomairy, M. Gates, S. Cayrols, D. Sukkari, K. Akbudak, A. YarKhan, P. Bagwell, and J. J. Dongarra, "Communication avoiding lu with tournament pivoting in slate," UTK/ICL, Tech. Rep. 18, ICL-UT-22-01, 2022-01 2022.

[22] G. Bosilca, A. Bouteiller, A. Danalis, M. Faverge, A. Haidar, T. Herault, J. Kurzak, J. Langou, P. Lemariner, H. Ltaief, P. Luszczek, A. YarKhan, and J. J. Dongarra, "Flexible development of dense linear algebra algorithms on massively parallel architectures with dplasma." Anchorage, Alaska, USA: IEEE, 2011, pp. 1432-1441.

[23] G. Bosilca, A. Bouteiller, A. Danalis, T. Herault, P. Lemariner, and J. J. Dongarra, "Dague: a generic distributed dag engine for high performance computing," Parallel Computing, vol. 38, no. 1-2, pp. 27-51, 00-2012 2012.

[24] M. Cosnard, E. Jeannot, and T. Yang, "Compact dag representation and its symbolic scheduling," Journal of Parallel and Distributed Computing, vol. 64, no. 8, pp. 921-935, 2004.

[25] L. Kale, A. Arya, A. Bhatele, A. Gupta, N. Jain, P. Jetley, J. Lifflander, P. Miller, Y. Sun, R. Venkataraman et al., "Charm++ for productivity and performance: a submission to the 2011 hpc class ii challenge," Parallel Programming Laboratory, Tech. Rep., pp. 11-49, 2011.

[26] N. Chalmers, J. Kurzak, D. McDougall, and P. T. Bauman, "Optimizing high-performance linpack for exascale accelerated architectures," 2023.

[27] A. Abdelfatih, T. Costa, J. J. Dongarra, M. Gates, A. Haidar, S. Hammarling, N. J. Higham, J. Kurzak, P. Luszczek, S. Tomov, and M. Zounon, "A set of batched basic linear algebra subprograms and lapack routines," ACM Trans. Math. Softw., vol. 47, no. 3, Jun. 2021. [Online]. Available: https://doi.org/10.1145/3431921

[28] N. Chalmers, J. Kurzak, D. Mcdougall, and P. Bauman, "Optimizing high-performance linpack for exascale accelerated architectures," in Proceedings of the International Conference for High Performance Computing, Networking, Storage and Analysis, ser. SC '23. New York, NY, USA: Association for Computing Machinery, 2023. [Online]. Available: https://doi.org/10.1145/3581784.3607066

[29] E. Ayguadé, N. Copty, A. Duran, J. Hoeflinger, Y. Lin, F. Massaioli, X. Teruel, P. Unnikrishnan, and G. Zhang, "The design of openmp tasks," IEEE Transactions on Parallel and Distributed Systems, vol. 20, no. 3, pp. 404-418, 2008.
[30] M. Bauer, S. Treichler, E. Slaughter, and A. Aiken, "Legion: expressing locality and independence with logical regions," in *SC Conference on High Performance Computing Networking, Storage and Analysis, SC '12, Salt Lake City, UT, USA - November 11 - 15, 2012*, 2012, p. 66.

[31] S. Thibault, “On runtime systems for task-based programming on heterogeneous platforms,” Habilitation à diriger des recherches, Université de Bordeaux, Dec. 2018. [Online]. Available: https://inria.hal.science/tel-01959127

[32] E. Anderson, Z. Bai, C. Bischof, L. S. Blackford, J. Demmel, J. J. Dongarra, J. Du Croz, A. Greenbaum, S. Hammarling, A. McKenney et al., LAPACK users' guide. SIAM, 1999.

[33] C. Rossignon, H. Pascal, O. Aumage, and S. Thibault, “A NUMA-aware fine grain parallelization framework for multi-core architecture,” in *PDSEC - 14th IEEE International Workshop on Parallel and Distributed Scientific and Engineering Computing - 2013*, Boston, United States, May 2013. [Online]. Available: https://inria.hal.science/hal-00858350

[34] M. Ruefenacht, M. Bull, and S. Booth, “Generalisation of recursive doubling for allreduce: now with simulation,” *Parallel Computing*, vol. 69, pp. 24–44, 2017. [Online]. Available: https://www.sciencedirect.com/science/article/pii/S0167819117301199

[35] K. Sala, X. Teruel, J. M. Perez, A. J. Peña, V. Beltran, and J. Labarta, “Integrating blocking and non-blocking mpi primitives with task-based programming models,” *Parallel Computing*, vol. 85, pp. 153–166, 2019.

[36] L. Courtès, “Building a secure software supply chain with gnu guix,” *arXiv preprint arXiv:2206.14606*, 2022.

[37] O. Aumage, E. Brunet, N. Furmento, and R. Namyst, “NewMadeleine: a Fast Communication Scheduling Engine for High Performance Networks,” in *Workshop on Communication Architecture for Clusters (CAC 2007), workshop held in conjunction with IPDPS 2007*, Long Beach, California, United States, Mar. 2007. [Online]. Available: https://inria.hal.science/inria-00127356

[38] M. Faverge, N. Furmento, A. Guermouche, G. Lucas, R. Namyst, S. Thibault, and P.-A. Wacrenier, “Programming heterogeneous architectures using hierarchical tasks,” *Concurrency and Computation: Practice and Experience*, vol. 35, no. 25, 2023. [Online]. Available: https://hal.science/hal-04088833