Mathieu Faverge

Nathalie Furmento

| Abdou Guermouche

Gwenolé Lucas
email: gwenole.lucas@inria.fr

Raymond Namyst

| Samuel Thibault

Pierre-André Wacrenier

Programming Heterogeneous Architectures Using Hierarchical Tasks

Keywords: Multicore, accelerator, GPU, heterogeneous computing, task graph, programming model, runtime system, dense linear algebra 1

Task-based systems have become popular due to their ability to utilize the computational power of complex heterogeneous systems. A typical programming model used is the Sequential Task Flow (STF) model, which unfortunately only supports static task graphs. This can result in submission overhead and a static task graph that is not well-suited for execution on heterogeneous systems. A common approach is to find a balance between the granularity needed for accelerator devices and the granularity required by CPU cores to achieve optimal performance. To address these issues, we have extended the STF model in the STARPU [1] runtime system by introducing the concept of hierarchical tasks. This allows for a more dynamic task graph and, when combined with an automatic data manager, it is possible to adjust granularity at runtime to best match the targeted computing resource. That data manager makes it possible to switch between various data layout without programmer input and allows us to enforce the correctness of the DAG as hierarchical tasks alter it during runtime. Additionally, submission overhead is reduced by using large-grain hierarchical tasks, as the submission process can now be done in parallel. We have shown that the hierarchical task model is correct and have conducted an early evaluation on shared memory heterogeneous systems using the CHAMELEON [2] dense linear algebra library.

However, the task-based paradigm presents several challenges when trying to efficiently exploit heterogeneous platforms. The computing resources of heterogeneous platforms have diverse characteristics and requirements. For example, GPU devices typically perform best with large data sets, while conventional CPU cores reach peak performance with fine-grain kernels working on a reduced memory footprint. Additionally, systems typically have a larger number of CPU units than GPUs. This means that having more small tasks may be necessary to increase performance.

Several efforts have been made to address this problem, such as finding the best trade-off between the optimal granularity of each device [START_REF] Agullo | A Hybridization Methodology for High-Performance Linear Algebra Software for GPUs[END_REF][START_REF] Buttari | A Class of Parallel Tiled Linear Algebra Algorithms for Multicore Architectures[END_REF][START_REF] Bosilca | Flexible Development of Dense Linear Algebra Algorithms on Massively Parallel Architectures with DPLASMA[END_REF][START_REF] Valero-Lara | sLASs: A fully automatic auto-tuned linear algebra library based on OpenMP extensions implemented in OmpSs[END_REF], or by aggregating CPU cores to process a task that was meant to be executed by an accelerator like a GPU [START_REF] Cojean | Resource aggregation for task-based Cholesky Factorization on top of modern architectures[END_REF][START_REF] Maroñas | Worksharing Tasks: An Efficient Way to Exploit Irregular and Fine-Grained Loop Parallelism[END_REF]. Alternatively, some preliminary work has considered splitting the tasks on CPU cores [START_REF] Wu | Hierarchical DAG scheduling for Hybrid Distributed Systems[END_REF]. While these approaches can be efficient in specific contexts such as dense linear algebra, they suffer from the fact that the task graph is static, making it impossible to select an alternative granularity for a given operation at runtime. For example, when designing linear algebra solvers based on low-rank approximation algorithms, it is almost impossible to predict the right DAG to ensure good numerical accuracy [START_REF] Carratala-Saez | Exploiting nested task-parallelism in the H-LU factorization[END_REF][START_REF] Akbudak | Tile Low Rank Cholesky Factorization for Climate/Weather Modeling Applications on Manycore Architectures[END_REF][START_REF] Carratalá-Sáez | Tiled Algorithms for Efficient Task-Parallel ℌ-Matrix Solvers[END_REF][START_REF] Augonnet | A Hierarchical Fast Direct Solver for Distributed Memory Machines with Manycore Nodes[END_REF][START_REF] Börm | Semi-Automatic Task Graph Construction for H-Matrix Arithmetic[END_REF].

In general, most modern task-based runtime systems suffer from a lack of dynamism in task-graph generation. Some programming models, such as [START_REF] Bauer | Legion: expressing locality and independence with logical regions[END_REF][START_REF] Wu | Hierarchical DAG scheduling for Hybrid Distributed Systems[END_REF][START_REF] Faverge | Designing LU-QR hybrid solvers for performance and stability[END_REF], support just-in-time DAG submission, but the generation either follows static rules or requires a significant amount of programming effort. Another major limitation of using task-based runtime systems is the overhead imposed by the runtime system to manage large applications. This is closely linked to the size of the internal DAG representation as tasks are submitted. Many runtime systems that rely on the STF model [START_REF] Allen | Optimizing Compilers for Modern Architectures: A Dependence-Based Approach[END_REF] (such as OPENMP, STARSS, STARPU) build this representation at submission time and may be negatively impacted by a large number of non-ready tasks. On the other hand, programming models where the task graph is discovered using a high-level description of dependencies never build this representation and are less affected, but tend to be more difficult to use. An example of such a programming model is the parameterized task graph (PTG) [START_REF] Cosnard | SLC: Symbolic scheduling for executing parameterized task graphs on multiprocessors[END_REF] used by PARSEC [START_REF] Bosilca | PaRSEC: Exploiting Heterogeneity to Enhance Scalability[END_REF]. Finally, regarding the sequential task flow model, task submission is generally done by a single thread, which may represent another bottleneck for the execution of large task graphs.

In this paper, we propose a new type of task, called hierarchical tasks, which can transform themselves into a new task graph dynamically at runtime. Programmers only need to provide hints indicating which tasks can be transformed into hierarchical tasks. The runtime system can then delay the submission of parts of the task graph to support dynamic implementation selection, parallelize the task insertion process, and greatly reduce the number of tasks in the runtime system. Our approach is similar to what is done in OPENMP for nested task-based parallelization schemes. However, we extend this approach to handle heterogeneous platforms while expressing fine-grain dependencies. In order to express more parallelism in sub-DAGs inserted with hierarchical tasks, we might want to partition and unpartition data. This is possible thanks to an advanced data manager, which can dynamically and asynchronously change the data layout. The system performs these operations automatically by inserting data management tasks that enforce data consistency. We also rely on this data manager to ensure that a DAG involving hierarchical tasks is consistent with a DAG created with the sequential task flow model, and therefore correct.

The proposed model associated with these hierarchical tasks address the abovementioned issues: 1) How to make the task graph more dynamic? 2) How to reduce the overhead of the runtime system? 3) How to overcome the intrinsic limitation of the sequential task flow submission process? While this model is generic and targets distributed heterogeneous architectures, this paper focuses on an initial implementation for shared memory heterogeneous architectures.

The contribution of this paper is two-fold: 1) We present an advanced data management engine that supports asynchronous data layout modification, and 2) We show how we extend the sequential task flow model to support hierarchical tasks and present our implementation within the STARPU runtime system.

The remainder of the paper is organized as follows. We present a general context with the sequential task flow model and the STARPU runtime system. Afterward, we describe the proposed advanced data manager, a crucial feature for implementing the hierarchical task model. We then introduce our hierarchical task paradigm and show the correctness of the model. Finally, we present a set of experiments showing the interest in the approach when performing dense linear algebra operations implemented within the CHAMELEON library [START_REF] Agullo | A Hybridization Methodology for High-Performance Linear Algebra Software for GPUs[END_REF].

RELATED WORK

In recent years, a significant amount of attention has been devoted to designing advanced runtime systems for modern heterogeneous architectures. Most of these efforts fall within the task-based paradigm, which is a natural way for scientific programming workflow (e.g. CILK [START_REF] Frigo | The implementation of the Cilk-5 multithreaded language[END_REF], OPENMP or INTEL TBB [START_REF] Reinders | Intel Threading Building Blocks: Outfitting C++ for Multi-Core Processor Parallelism[END_REF] for multicore machines, APC [START_REF] Hartley | Improving performance of adaptive component-based dataflow middleware[END_REF], CHARM++ [START_REF] Kunzman | Programming heterogeneous clusters with accelerators using object-based programming[END_REF][START_REF] Acun | Parallel Programming with Migratable Objects: Charm++ in Practice[END_REF], HPX [START_REF] Kaiser | HPX -The C++ Standard Library for Parallelism and Concurrency[END_REF], KAAPI/XKAAPI [START_REF] Hermann | Multi-GPU and Multi-CPU Parallelization for Interactive Physics Simulations[END_REF], LEGION [START_REF] Bauer | Legion: expressing locality and independence with logical regions[END_REF], PARSEC [START_REF] Bosilca | PaRSEC: Exploiting Heterogeneity to Enhance Scalability[END_REF], STARPU [START_REF] Augonnet | StarPU: A Unified Platform for Task Scheduling on Heterogeneous Multicore Architectures[END_REF] or STARSS [START_REF] Ayguadé | An Extension of the StarSs Programming Model for Platforms with Multiple GPUs[END_REF] for heterogeneous configurations).

The common point among these runtime systems is that they all use high-level descriptions of dependencies to build the task graph at runtime and then schedule the corresponding computations on available resources. Various approaches are used to construct the task graph. For instance, most of the previously mentioned runtime systems rely on the so-called Sequential Task-Flow model (e.g. OPENMP, STARSS, STARPU) to build the task graph. By using data access modes and a sequential submission order, dependencies between tasks can be inferred through data dependency analysis [START_REF] Allen | Optimizing Compilers for Modern Architectures: A Dependence-Based Approach[END_REF]. On the other hand, runtime systems such as PARSEC are based on the parameterized task graph programming model (PTG) [START_REF] Cosnard | SLC: Symbolic scheduling for executing parameterized task graphs on multiprocessors[END_REF], where the task graph is unrolled at runtime, using a high-level description of the dataflow corresponding to the computations. Other runtime systems use different paradigms to express computations, such as LEGION, which describes logical data regions representing the data flow and dependencies between tasks. The various programming models have different levels of ease of use and the amount of overhead they impose on the underlying runtime system. For example, the sequential task flow model is easy to use as programmers only need to provide the sequential implementation of their application and then add data access modes. However, this method comes with higher overhead in the runtime system. On the other hand, the parameterized task graph approach requires users to express their computations in a subtle high-level formalism where the dataflow is explicitly described, but has less overhead on the runtime system [START_REF] Hoque | Dynamic Task Discovery in PaRSEC: A Data-Flow Task-Based Runtime[END_REF].

Several efforts have been made to address the problem of reducing the overhead of task-based runtime systems, which mainly focus on those based on the sequential task flow model or increasing the amount of parallelism provided by these systems. In [START_REF] Álvarez | Advanced Synchronization Techniques for Task-Based Runtime Systems[END_REF], the authors analyzed the limiting factors in the scalability of a task-based runtime system and proposed individual solutions for each of the identified challenges, such as a wait-free dependency system and a scalable scheduler design based on delegation instead of work-stealing. Other approaches focus on advanced dependency management. For example, in [START_REF] Elshazly | Accelerated execution via eager-release of dependencies in task-based workflows[END_REF], the authors proposed an eager approach for releasing data dependencies, where tasks are launched for execution as soon as their data requirements are met instead of waiting for predecessor tasks to finish execution. In [START_REF] Maroñas | Worksharing Tasks: An Efficient Way to Exploit Irregular and Fine-Grained Loop Parallelism[END_REF], worksharing tasks were introduced, which internally leverage worksharing techniques to exploit fine-grained structured loop-based parallelism without requiring a barrier. The closest contribution to our proposition, with regard to task dependencies, was introduced in [START_REF] Perez | Improving the Integration of Task Nesting and Dependencies in OpenMP[END_REF] with the concept of weak dependencies, which is an extension of the OPENMP model that supports fine-grained dependencies not only between sibling tasks but also between tasks with any family relationship. Our contribution is a generalization of the weak dependency concept to the heterogeneous case where the underlying hardware does not ensure memory consistency, thus requiring an advanced data manager (see Section 4). Additionally, in ParSEC [START_REF] Wu | Hierarchical DAG scheduling for Hybrid Distributed Systems[END_REF] and XKaapi [START_REF] Gautier | XKaapi: A Runtime System for Data-Flow Task Programming on Heterogeneous Architectures[END_REF], some preliminary work targeting heterogeneous architectures has considered splitting tasks when assigned to CPU cores.

Hierarchical tasks support

Runtime

Fine-grain Automatic Data Heterogeneity Dependencies Management Several efforts have been made to allow task-based runtime systems to have more dynamic capabilities in terms of task management and generation. In TaskFlow [START_REF] Huang | Taskflow: A Lightweight Parallel and Heterogeneous Task Graph Computing System[END_REF], advanced tasking schemes are introduced, including dynamic, composable, and conditional tasking. Dynamic tasking, in particular, allows for the dynamic generation of a sub-DAG from a given task. However, synchronization is added at the end of each hierarchical task to ease dependency management, and programmers are responsible for data management and changing data layout as needed. In [START_REF] Kim | IRIS: A Portable Runtime System Exploiting Multiple Heterogeneous Programming Systems[END_REF], the IRIS runtime is introduced, which can perform dynamic task partitioning either through user input or automatically via a polyhedral compiler, but no details are provided on how dependencies are handled in this context. Finally, libtask, an advanced runtime system supporting hierarchical tasks in the context of low-rank linear algebra solvers is presented in [START_REF] Augonnet | A Hierarchical Fast Direct Solver for Distributed Memory Machines with Manycore Nodes[END_REF]. This work introduces hierarchical tasks, and dependencies are expressed at the finest level. The correctness of the produced DAG is considered through extra automatic dependencies, but data management is straightforward as data partitioning is performed statically at the beginning of the execution.

PaRSEC [8] ✗ ✗ ✓ OmpSs [30] ✓ ✗ ✗ TaskFlow [32] ✗ ✗ ✓ IRIS [33] ✗ ✓ ✓ libtask [12] ✗ ✓ ✓
Table 1 provides a summary of how the runtime systems mentionned in this section compare with respect to the functionality we want to achieve with hierarchical tasks. We focus on three main features. The first one is the presence of fine-grain dependencies between distinct sub-graphs. The second one is a data manager that can adjust the data layout around hierarchical tasks without user supervision. The third one is the support of heterogeneous architectures.

BACKGROUND

The subsection below presents the basics of the sequential task flow model and the corresponding sequential consistency. We then briefly present our target runtime system STARPU.

The Sequential Task Flow Model

The Sequential Task Flow (STF) programming model entirely relies on Sequential Consistency using only implicit dependencies. The STF model, therefore, simply consists of submitting a sequence of tasks through a non-blocking function call that delegates the execution of the task to the runtime system. Upon submission, the runtime system adds the task to the current DAG, and its dependencies are automatically computed through data dependency analysis [START_REF] Allen | Optimizing Compilers for Modern Architectures: A Dependence-Based Approach[END_REF]. The actual execution of the task is then postponed until its dependencies are satisfied. This paradigm is sometimes referred to as Superscalar since it mimics the operation of superscalar processors where instructions are issued sequentially from a single stream but can be executed in a different order and, possibly, in parallel depending on their mutual dependencies. Figure 1 shows a dummy sequential algorithm and its corresponding STF version. Instead of making three function calls (F, G, H), the equivalent STF submits the three corresponding tasks. The data onto which these functions operate as and their access mode (Read, Write or Read/Write) are also specified. Because task G accesses data x after task F has accessed it in Write mode, the runtime infers a dependency between tasks F and G. Similarly, a dependency is inferred between tasks F and H due to data y. Figure 1 (right) shows the DAG corresponding to this STF dummy code. In the STF model, one thread is responsible for submitting the tasks; we refer to this thread as the master thread. Worker threads instead achieve the execution of tasks. The function called at the end of the STF pseudo-code is simply a barrier that prevents the master thread from continuing until all the submitted tasks are executed.

We present in Sections 4 and 5 how hierarchical tasks can be introduced in task-based runtime systems relying on the STF model without breaking the sequential consistency (i.e. ensuring a task execution order coherent with a sequential execution) to ensure the correctness of the DAG.

The STARPU Runtime System

STARPU [START_REF] Augonnet | StarPU: A Unified Platform for Task Scheduling on Heterogeneous Multicore Architectures[END_REF] is a library that provides programmers with a portable interface for scheduling dynamic graphs of tasks onto a heterogeneous set of processing units (i.e. CPUs and GPUs). The two basic principles of STARPU are, first, that tasks can have several implementations for some or each of the various heterogeneous processing units available in the machine. Secondly the runtime system handles necessary data transfers to these processing units transparently. STARPU tasks are defined as multiversion kernels, gathering the different implementations available for CPUs and GPUs, associated to a set of input/output data. To avoid unnecessary data transfers, STARPU allows multiple copies of the same registered data to reside simultaneously in different memory locations as long as it is not modified. Asynchronous data prefetching and write-back is also used to hide memory latencies.

STARPU is a platform for developing, tuning and experimenting with various task scheduling policies in a portable way. It provides an abstract view of the machine by relying on the notion of worker to describe the computing resource (e.g. CPU, GPU) in charge of executing the tasks. Implementing a scheduler consists of creating a set of queues, associating them with the different processing units, and defining the code that is triggered each time a new task is ready to be executed, or when a processing unit is about to go idle. Various designs can be used to implement queues (e.g. FIFOs or stacks), and they can be organized according to different topologies. Several built-in schedulers are available, ranging from greedy and work-stealing based policies to more elaborate schedulers implementing variants of the Minimum Completion Time (MCT) policy [START_REF] Topcuouglu | Performance-effective and low-complexity task scheduling for heterogeneous computing[END_REF].

AUTOMATIC DATA MANAGEMENT

Data handling is at the heart of STARPU to automatically infer dependencies between tasks in the STF model and automatically manage data transfers between the different memory banks of a distributed/heterogeneous system. Applications must register the data that are handled by the tasks. To do so, STARPU provides an opaque data structure called handle, which is an abstract view of registered data. Handles along an access mode (read-only, read-write, ...) are used as task parameters. A task must access data through the associated handle.

To ease data manipulation, STARPU brings the notion of data filter, a tool to partition data associated with a handle into sub-data parts associated with new sub-handles. Indeed, instead of registering all data sub-sets independently, it is often more convenient to register a large piece of data and to recursively partition it.

Once a handle is partitioned, we can remark that the same data can be designated simultaneously by several different handles. This may be very convenient when programmers want to use data with several views. Moreover, data in read-only access mode can advantageously be accessed simultaneously at different partitioning levels by several tasks on other devices. However, when data is accessed in write access mode, this access must be exclusive for coherency purposes. This property is ensured by STARPU when a single partitioning is used for data, but may be violated when several handles point to the same data. To deal with this problem, STARPU provides functions to invalidate other handles ensuring they can not be used to access their underlying data, and unpartition sub-handles back into the main handle to gather the sub-data.

However, users' feedback showed that maintaining consistency between different partitions is difficult. This is why we propose a mechanism to automate the management of several simultaneous partitions. This is done by letting STARPU automatically insert partition or unpartition tasks as needed.

We will now present the modifications we made in STARPU to allow this automatic data management. First, programmers need to define the partitioning scheme through the plan operation, which aims to declare the partitioning to STARPU, and can be seen as the declaration of a new set of sub-handles. One can even partition recursively and use handles at different levels of the recursion. Once a plan is performed, it is possible to submit tasks using the initial handle or any of the sub-handles, even if the actual partitioning has not been done yet. Furthermore, several partitioning schemes can be planned simultaneously as illustrated in Figure 2 (see the lines 12 to 25). The data manager will then handle the actual partitioning tasks and data coherency. The given example shows a matrix on which two partition plans are defined, one with vertical stripes and one with horizontal stripes. The matrix is first initialized through its root handle, then modified using vertical partitioning, and checks are performed in horizontal and vertical stripes.

This mechanism is very natural to use for the programmer because it allows them to deal with the registration and partitioning of their data once and for all. They are then free to use any handle within the tasks they insert and the runtime system will adapt the data layout by itself.

At runtime, STARPU will introduce coherency synchronization. When a task is ready to be executed, STARPU must ensure that the partition associated with each handle it uses is valid. STARPU will allow different partitioning to coexist if data is accessed in read-only mode. As soon as data is accessed in read-write mode, STARPU will automatically (and recursively) unpartition sub-data and activate only the partitioning leading to the handle being written to. To illustrate the behavior of the code given in Figure 2, we provide in Figure 3 the evolution of the data layout and the task graph by emphasizing the submission of automatic partition/unpartition tasks by the runtime system. We can see that after executing line 27 (see Figure 3a), the plan operations have been performed (see the dotted matrix layouts), and the initialization task has been submitted. When the first task using a vertical stripe, which modifies the matrix, is submitted on line 31, the runtime system will automatically insert the corresponding partitioning task (see Figure 3b). The same scheme is then applied at the submission of the tasks working on the horizontal layout (at line 35) and vertical layout (at line 37) in read-mode. One should note that 𝐶 𝑣 1 and 𝐶 𝑣 2 share the same vertical layout as 𝑉 1 and 𝑉 2 , so no partition operation is needed for these tasks. On the contrary, tasks 𝐶 𝐻 1 and 𝐶 𝐻 2 do not share any handles with those using the vertical layout. However, the data manager knows these handles share a joint ancestor (the whole matrix). Thus, it will insert, the unpartition/partition tasks as needed to make the data available to the tasks using the horizontal layout. This is illustrated in Figure 3c where the 𝑈 𝑣 and 𝑃 ℎ tasks are inserted, making the tasks using the horizontal layout depend on them. Finally, when the partition needs to be cleaned, the final unpartition task is inserted (see Figure 3d).

The previous example illustrates the general behavior of the data manager. More precisely, during the submission of tasks, each handle in the partitioning hierarchy can be either inactive (one cannot access the piece of data), read-write-active (one can read/write to the piece of data or a sub-part of it), or read-only-active (one can only read from the piece of data or a subpart of it). The main handle at the root of the partitioning hierarchy is always read-write-active. Each handle in the hierarchy, when active, is additionally either unpartitioned (one can read/write the piece of data itself), read-write-partitioned (one can only write to the sub-pieces of data), or read-only-partitioned (one can read the piece of data or sub-pieces of data) ; when it is partitioned, its children sub-handles in the hierarchy are active. When submitting a task that accesses a handle within the hierarchy, STARPU will automatically ensure that the handle is active. This possibly requires recursively making its ancestors active by submitting partitioning tasks for them, possibly starting right from the root handle of the hierarchy. This also possibly requires recursively submitting unpartitioning tasks for some sub-handles which were previously written to. In the case of the transition from Figure 3b to Figure 3c, STARPU indeed had to submit the unpartition of the root handle, and repartition it. It should be noted that all this mechanism is done completely asynchronously: STARPU just inserts data coherency synchronization tasks within the task graph.

𝐼 𝑃

THE HIERARCHICAL TASK PARADIGM

To extend the sequential task flow model with hierarchical tasks, we first need to define what precisely a hierarchical task is. Formally, a hierarchical task is simply a regular task that can, at runtime, submit a sub-DAG instead of performing actual computations. Processing a hierarchical task consists of submitting its corresponding task sub-graph, its outgoing dependencies can be released at the end of that submission process. Furthermore, to ensure portability concerning heterogeneous platforms, coherency synchronization tasks are submitted along the sub-graph to ensure a correct execution by connecting the sub-DAG with the rest of the DAG. In the end, hierarchical tasks represent an elegant answer to 1) the problem of adapting the granularity of tasks to the device executing them, 2) the question of the reduction of the number of active tasks in the runtime system, 3) the problem of the dynamic selection of the implementation of a given operation in the application. Introducing hierarchical tasks in a task-based runtime system like STARPU needs to respect the following constraints, which aim at having a general implementation of such a paradigm.

1. The depth of the hierarchy is not limited (a hierarchical task can spawn a DAG that contains hierarchical tasks).

2. Programmers express their task graph at the highest level and only annotate some tasks as possibly hierarchical.

3. Data management needs to be transparent to programmers. [START_REF] Bosilca | Flexible Development of Dense Linear Algebra Algorithms on Massively Parallel Architectures with DPLASMA[END_REF]. Task dependencies always have to be inferred at the finest grain (i.e. the deepest level of the hierarchy).

Properties 1 and 2 ensure a general programming interface with no limitations from the programmer's point of view. Property 3 is related to data layout management. It is necessary for the heterogeneous case where memory is not shared between all the computing resources, meaning that explicit data movements must be performed. This property will be handed to the advanced data manager introduced in Section 4. Finally, Property 4 is related to parallelism management. It only says that we do not want a fork-join approach, as there is no barrier at the end of hierarchical tasks.

To use hierarchical tasks, the programmer inserts a task that works on the coarsest level of data and adds a few hints to make it hierarchical. These hints consist of a graph generating function and the sub-data on which the sub-graph will work. Figure 4 is an example of how the code from Figure 2 could be extended to use hierarchical tasks. Instead of using for loops to submit tasks working on the sub-data, we move them into two new functions modify_gen_dag and check_gen_dag and submit a single task, working on the main data, in their place. These functions, with the appropriate sub-data, are passed as arguments to each corresponding task insertion, making these tasks hierarchical for STARPU.

[ Figure 5a is an execution scenario for a given task graph, where some tasks have been tagged as possibly hierarchical (shown as blue nodes). The state of each task (i.e. node in the graph) is described by its border: 1) a ready task has a green border (all dependencies are met), 2) a not-ready task has a red border (some dependencies are unsatisfied), 3) an already executed task has a black border. Thus, we can see in Figure 5a that 𝑇 1 has completed its execution, making 𝑇 2 and 𝐻 1 ready, while the remaining tasks are not ready for execution. 𝑇 2 and 𝑇 3 execute as normal tasks, while 𝐻 1 is processed, i.e. its corresponding sub-DAG is submitted, resulting in Figure 5b. The dependency between 𝐻 1 and 𝐻 2 is then released, making 𝐻 2 ready for processing. Furthermore, we can see that after the processing of 𝐻 2 (see Figure 5c), the dependencies between the resulting submitted tasks are inferred by the runtime system at the deepest level of the hierarchy. This will allow for a good pipelining of the processing of the two task graphs resulting from 𝐻 1 and 𝐻 2 .

We now have to consider how the data coherency will be achieved between the DAG and the sub-DAGs. Introducing hierarchical tasks in a task-based runtime system requires changing the granularity of data dynamically at runtime each time a hierarchical task has to be executed. Indeed, to be able to express parallelism within the hierarchical task, the task's input data has to be either partitioned or unpartitioned (see Section 4 for more details) depending on the level of the task in the hierarchy. An example of this aspect is provided in Figure 6, where a DAG composed of a single hierarchical task (implementing a matrix product (GEMM)) is provided. We can see that when task 𝐻 1 is processed, it submits its corresponding DAG. However, to exhibit parallelism, data has to be partitioned (in this case, the original matrices are partitioned in a 4 by 4 layout). This illustrates the fact that introducing hierarchical tasks requires advanced partitioning/unpartitioning operations for each input and output data.

FIGURE 6

Example of a hierarchical task corresponding to a matrix product operation 𝐶 = 𝐶 + 𝐴 × 𝐵 and its corresponding DAG.

Thus, the general approach is to have a data management operation (either partitioning or unpartitioning) preceding each individual task (resp. hierarchical task) when needed. We propose to automatically insert a data management task ahead of a task requiring data not in the correct layout by relying on the data manager introduced in Section 4. Coming back to Figure 5, we can see, in Figure 5b, the insertion of the partitioning task 𝑃 ahead of the sub-graph produced by 𝐻 1 . This can be done either during the processing of 𝐻 1 or delayed until the first task of the task graph corresponding to 𝐻 1 is ready for processing. Conversely, the unpartitioning task 𝑈 is inserted before the final task 𝑇 4 . The insertion of such a task can be delayed as late as the moment when all the task dependencies of 𝑇 4 are fulfilled. It is important to emphasize that the insertion of the unpartition task needs to know whether the task 𝑇 4 is hierarchical or regular. These data management tasks are inserted automatically by the runtime ahead of the task requiring access to data whenever the current layout of the data is not the one needed by the task (see Section 4 for more details). We can also notice that there is no data management task between the sub-graphs produced by 𝐻 1 and 𝐻 2 since they share the same data layout. Finally, it is essential to emphasize that hierarchical tasks are executed when their dependencies are fulfilled. However, the actual computations tasks submitted by these hierarchical tasks are executed whenever they are ready. Thus, we need to ensure a correct order of the actual computations.

We show in the next section how we ensure that the DAG generated when relying on hierarchical tasks is coherent with the sequential task flow model, even though we are now able to make task submission parallel, since the processing of various hierarchical tasks can be performed concurrently.

Ensuring the Correctness of the DAG

We will now show why the hierarchical task model to extend the sequential task flow model produces a correct DAG whatever the depth of the hierarchy and the interactions between regular and hierarchical tasks. Firstly, as stated above, the sequential task flow model infers the dependencies from data access modes of individual tasks while relying on sequential consistency (the insertion order of tasks provides a natural ordering to infer task dependencies). Thus, when introducing hierarchical tasks, one major issue is ensuring that the DAG produced after processing all the hierarchical tasks is correct with respect to a sequential submission scheme (in the end, we want the DAG to be equivalent to the one using a single submission thread). Indeed, when relying on hierarchical tasks, the submission can now be made in parallel, while in the sequential task flow model, the submission is made by a single entity. In other words, how can the correctness of the DAG be ensured?

𝑃 𝑈 𝑇 1

(c) Insertion of the unpartition task when 𝑇 1 is ready.

FIGURE 7

Example of a scenario where a task follows a hierarchical task.

We must show that the dependencies respect the sequential task flow model to answer this question. We discuss four simple scenarios, which are building blocks for any general DAG to show its correctness.

§ ¦ ¤ ¥ 𝑇 𝑇
This is a straightforward scenario that corresponds to the sequential task flow model. The sequential consistency and the submission order of the tasks will ensure the correctness of the DAG.

§ ¦ ¤ ¥ 𝑇 𝐻
This represents the case, where a regular task in the original DAG is followed by a hierarchical one that depends on it. The sequential task flow model used for submitting the original DAG ensures that the hierarchical task cannot start its processing before the regular task has been completed, ensuring the dependencies' correctness.

§ ¦ ¤ ¥ 𝐻 𝑇
This represents the case where a hierarchical tasks precedes a regular task in the original DAG. An example of this scenario is provided in Figure 7. The main problem in this scenario is that the regular task is by construction submitted before the tasks resulting from the hierarchical task (𝐻 1 in Figure 7). This may violate the order required by sequential consistency. However, in this specific situation, the hierarchical task has changed the data layout before it starts its processing (see Figure 7b). Thus, the task following the hierarchical task (𝑇 1 in Figure 7) will request the data layout to be changed. The data manager will then automatically submit data management tasks to turn back data to their original layout. These data management tasks will be inserted ahead of the task in the DAG and will depend on the data produced by the DAG resulting from the processing of the hierarchical task. This corresponds to the configuration depicted in Figure 7c where the unpartition task 𝑈 is inserted so that 𝑇 1 depends on 𝑈 , which needs the data produced by the previous tasks. Therefore, the data management tasks will ensure that the regular task 𝑇 cannot start its execution before the completion of the DAG submitted by the hierarchical task. This illustrates the data manager's central role, which has been discussed in Section 4. This represents the case where in the original DAG, a hierarchical task is followed by another one that depends on it. In this case, since the dependency between the two hierarchical tasks is not released until the first one has completed its processing, the tasks resulting from the two hierarchical tasks are correctly ordered, making the dependencies between these tasks coherent with the sequential consistency. This is illustrated in Figure 8 where initially, two hierarchical tasks 𝐻 1 and 𝐻 2 are submitted (see Figure 8a). Then 𝐻 1 is processed (see Figure 8b). Note that in the example, we assume that the data was previously unpartitioned, and thus a data partitioning task 𝑃 1 is needed before the DAG corresponding to 𝐻 1 . Afterward, 𝐻 2 is processed (see Figure 8c) and it does not require any data layout modification. Note that, each individual task produced by a hierarchical task can itself be hierarchical, and the same rules can be applied recursively to ensure the correctness of the DAG. This is illustrated in Figure 8d where the first task submitted by 𝐻 1 , which will be referred to as 𝐻 11 , is decided to be hierarchical (at runtime) and is executed. We can also see the partitioning task 𝑃 2 , which was automatically inserted by the data manager to partition the data further. The resulting task graph is coherent with the sequential task flow paradigm. Finally, note that the unpartitioning 𝑈 2 operation corresponding to 𝑃 2 is submitted in a second step when the final tasks of 𝐻 1 require data (see Figure 8e). This automatic insertion of 𝑈 2 is vital to turn back data into their original layout before the final tasks of 𝐻 1 can start their execution.

EXPERIMENTAL EVALUATION

To illustrate the potential of hierarchical tasks, we apply them in a dense linear algebra context using the CHAMELEON library [START_REF] Agullo | A Hybridization Methodology for High-Performance Linear Algebra Software for GPUs[END_REF]. The objective is to validate the hierarchical task model by using it to create correct DAGs with various level of task granularity. We use various matrix partitioning schemes to decide where to place hierarchical tasks to adress certain limitations in the applications we chose and potentially improve their performances. To do so, we adapted the matrix descriptors to describe hierarchical data structures. Three partitioning schemes have been implemented for these experiments. They are described further down for the experiment they were used in. As explained in Section 4, all these partitions are only planned, and the decision to submit a hierarchical task is taken at runtime. That decision is based on the data structure of the tiles involved in the computations. A sub-graph will be submitted if there is a partitioning plan for all of them. Otherwise, the tasks will operate at the highest joint data level.

Concerning hierarchical variants, we will use the following notation: 𝑥∕𝑦∕𝑧∕.... This means that each initial tile is of size 𝑥 and is partitioned according to one of the three schemes mentioned above into tiles of size 𝑦, which are in turn split into tiles of size 𝑧 etc.

The following experiments were conducted on two architectures, presented in Figure 9 shows the performance of a single dense DGEMM kernel for different matrix sizes (the test involves three square matrices) on both platforms. We present with a solid line the actual performance for both a single CPU and a single GPU. We can see that on the CPU side, the processor's peak performance is achieved for matrix sizes not exceeding 500. Providing larger matrices does not improve the absolute performance on a single core. On the GPU side, we can observe that the larger the matrix size, the higher the performance. This experiment aims to find the best trade-off for the task granularity between the CPU cores and GPU devices. Having a granularity of tasks not large enough will limit GPU performance. On the other hand, using very coarse grain tasks limits parallelism and may affect CPU performance. Thus, in future experiments, we will mainly use a tile size of 2880 (resp. 5700) on the INTEL-V100 (resp. AMD-A100) platform.

FIGURE 9

Performance of a matrix product (DGEMM) kernel. In the ideal setting a single kernel runs on a single ressource (CPU core or GPU device). In the noisy setting every CPU core runs 10 independent DGEMM kernels.

The dotted line is the normalized performance of a batch of independent DGEMM kernels spread evenly across all CPU cores. This accounts for the performance degradation resulting from memory contention in shared memory hierarchy levels (L3 cache, for example), and for the CPU frequency reduction due to the higher number of concurrent resources exploited.

On both architectures, STARPU and CHAMELEON are compiled with GCC 10.3.0, Nvidia CUDA 11.2, and the Intel MKL 2020. STARPU has been configured to use a single stream per GPU and to pipeline four events per stream. Unless otherwise specified, we use all CPU cores and GPU devices available on the node.

Overhead

To evaluate the overhead induced by hierarchical tasks, we consider the graph of a matrix-matrix multiplication (GEMM) using a tile size of 960. The partitioning scheme used in this experiment recursively splits each matrix tile. An example is shown in Figure 11a with one recursion level: each matrix tile is divided into a finer tile size. Figure 10 compares the graph submission time per computational task in two configurations. The '960' curve represents the non-hierarchical case. The '960/960' curve shows the worst possible scenario: the DAG is composed only of hierarchical tasks and each submits exactly one task when processed. This doubles the number of tasks submitted and heavily increases the workload of the data manager, making the submission time per computational task roughly 3.5 times slower. Finally, the '2880/960' curve is a more realistic scenario, where the graph is first submitted at coarse grain (with a tile size of 2880) and then refined down to the same granularity as the previous configurations (960). In this case, each hierarchical task submits ⌈2880∕960⌉ 3 = 27 regular tasks when processed, thus amortizing the overhead induced by the management of hierarchical tasks.

Dense matrix product (GEMM)

For the matrix product kernel 𝐶 = 𝛼𝐴𝐵 + 𝛽𝐶, the main objective of the hierarchical scheme is to balance the workload between CPU cores and GPU devices based on their respective computing power. A custom partitioning scheme is used. We first define a computing power for each resource,  𝐶𝑃 𝑈 , and  𝐺𝑃 𝑈 respectively, and then we estimate the computing power ratio of the CPU resources: We finally compute the number of tiles for which the partition planning has to be done with respect to 𝑅: ⌊𝑅 × 𝑁𝑡 2 ⌋, where 𝑁𝑡 2 is the total number of tiles of the 𝐶 matrix. Table 3 summarizes the different values of 𝑅 we selected for our experiments. 𝑅1 is computed from the best performances reached for tiled matrix multiplication using all the CPU cores without the GPU devices and vice-versa. 𝑅2 is computed from the plateau values from Figure 9, using the noisy setting in the CPU case. We evaluate the behavior of the GEMM operation on those matrices, using one and two GPUs (Figures 12a and12b). On the INTEL-V100 platform, this strategy allows the hierarchical versions to outperform the standard ones on medium matrix sizes by providing the CPU cores of the node with tasks working on smaller tiles. When using large tiles (2880) on bigger matrix sizes, STARPU's scheduler can start affecting more work to the CPUs without penalizing execution time and catches up with the hierarchical versions. On the AMD-A100 platform, the performance ratio between the CPU cores and the GPU devices is so high that the contribution of the CPU cores is marginal. Once the matrices become large enough, the percentage of smaller tiles becomes a disadvantage and negatively impacts the performance. This could be compensated by partitioning fewer tasks as the matrix size increase. However, even in this situation, the hierarchical version using 𝑅1 can achieve similar performance to standard CHAMELEON. A common trend between both platforms is that the 𝑅1 ratio provides better results than 𝑅2, because it was computed using values more representative of the disparity in computing power between CPUs and GPUs for this type of operation. Generally speaking, the dense matrix product is a very regular operation that is rather unfavorable to the use of hierarchical tasks. Despite that, the hierarchical variants have a good behavior and can achieve high performance with a simple partitioning strategy.

𝑅 = # 𝐶𝑃 𝑈 ×  𝐶𝑃 𝑈 (# 𝐶𝑃 𝑈 ×  𝐶𝑃 𝑈 ) + (# 𝐺𝑃 𝑈 ×  𝐺𝑃 𝑈 ) .

Platform

FIGURE 12

Comparison of the hierarchical strategy to the original one on the matrix-matrix multiplication (DGEMM) kernel, with a fixed percentage of hierarchical tasks (see Table 3).

We present experimental results in the forthcoming section illustrating the hierarchical tasks' behavior on more advanced dense linear algebra operations to illustrate the flexibility of hierarchical tasks better. First, we consider operations relying on Cholesky decomposition (POTRF): POSV (linear system solving, in this case, of a single vector) and POINV (matrix inversion). These operations have complex task graphs, and in the case of POINV, validate the anti-dependency problem (WRITE after READ). Furthermore, we provide an experimental evaluation using the LU factorization without partial pivoting (GETRF_nopiv) which has a wider task graph than the Cholesky factorization and for which the criticality of some tasks (mainly those corresponding to the tiles on the diagonal) can strongly affect performance if delayed.

FIGURE 13

Performance evaluation of Cholesky type operations (DPOTRF, DPOSV, DPOINV) with diagonal distribution of the hierarchical tasks.

Dense Cholesky factorization-based operations

The partitioning scheme we implemented for the following experiments splits recursively the matrix along the diagonal, subdiagonal and superdiagonal, to mimic -Matrices algorithms and use finer granularities on the critical path. An example is shown in Figure 11b. Figure 13 shows the results of operations relying on the Cholesky decomposition: POTRF (actual Cholesky decomposition), POSV (linear system solving, with a single right-hand side) and POINV (matrix inversion). Contrary to the previous DGEMM experiment, this partitioning scheme doesn't simply attempt to supply all resources with the right amount of work. Instead, it tries to improve performance on the critical path of the operation, because CHAMELEON requires all POTRF kernels (which are on the critical path of the factorization and located on the diagonal tiles) on CPU cores. This is illustrated in Figure 13 where we can see that the granularity having the best behavior used for the standard version use smaller tiles than the one relying on the hierarchical tasks on both platforms. Thanks to hierarchical tasks, we can partition the tiles along the diagonal and split those large tasks into sub-graphs with a smaller granularity, allowing for better CPU utilization on the critical path. For all the kernels and platforms, enabling hierarchical tasks enhances performance for most matrix sizes. We can also observe that the regular version can catch up for huge matrix sizes. The sudden drop observed at the end of some non-hierarchical curves is explained by a conflict between the STARPU scheduler data prefetching and eviction in GPU memory, as explained in [START_REF] Stanisic | Modeling and Simulation of a Dynamic Task-Based Runtime System for Heterogeneous Multi-Core Architectures[END_REF]. The experimental results illustrate the interest in hierarchical tasks for tackling the granularity problem of heterogeneous architectures. This highlights the fact that hierarchical tasks can find a better trade-off to fully exploit the whole platform: each resource is used with a task granularity adapted to its characteristics.

Dense LU factorization without pivoting

We provide in Figure 14 an experimental evaluation of hierarchical tasks to enhance the dense LU factorization without pivoting (DGETRF_nopiv). This dense linear algebra kernel exhibits a high amount of parallelism. It requires the tiles on the diagonal to be processed as fast a possible to release many dependencies. Once again, we can observe that hierarchical tasks improve performance on both platforms. We can also observe that the non-hierarchical version of CHAMELEON achieves lower performance because of the critical tasks on the diagonal: since the matrix partitioning is uniform, the best tiling uses a large granularity to exploit the GPU devices fully. Finally, we can observe that the gap between the hierarchical tasks and standard Chameleon is more significant on the more heterogeneous AMD-A100 platform. This is mainly because, on this platform, large tiles (5760) are needed to exploit the potential of the A100 devices fully. Thus, not using hierarchical tasks to adapt the granularity for each device penalizes performance. All in all, the experimental results illustrate the flexibility of the hierarchical task model and its ability to find a better tradeoff to tackle the granularity problem for heterogeneous platforms.

FIGURE 15

Performance evaluation of Cholesky type operations (DPOTRF, DPOSV, DPOINV) with diagonal distribution of the hierarchical tasks. The non-hierarchical variant relies on multiple streams. Each point of the plot is the best performance obtained for each variant considering tile size, number of streams, and partitioning scheme (for the hierarchical variant).

Going further

In the previous sections, we fairly compared the standard STARPU task paradigm and the hierarchical one. The results showed the benefits of relying on hierarchical tasks from the performance point of view. Still, we should compare the hierarchical variant with the standard one when using multiple GPU streams to be conclusive. Indeed, GPU streams are, in a sense, a way to achieve task parallelism within the GPU device: kernels submitted on different streams can be executed concurrently on the GPU device. This may enhance GPU occupancy and thus improve performance. In the context of our experimental testbed, using GPU streams will allow the baseline variant to find a better trade-off between the granularities adapted to CPU cores and GPU devices. It will thus be possible to use smaller tile sizes while relying on multiple streams to increase the device efficiency while shortening the critical path (the larger the tile size, the longer the critical path). We report in Figures 15 and16 the experimental results obtained when relying on streams for the non-hierarchical variant on both platforms using all resources. For both figures, each dot in the plot corresponds to the best performance obtained for the variant while considering different tile sizes (320, 640, 960, 1920, 2880, 5760), number of streams (1, 2, 4, 8) and partitioning scheme for the hierarchical variant (1 and 2 levels of hierarchy). It is important to note that we did not use streams for the hierarchical variant. First, we can observe that the use of GPU streams enhanced the performance of the non-hierarchical variant. We can also see in Figure 15 that both variants have a very close behavior for all three considered kernels. Moreover, the baseline approach relying on streams is slightly better for small to medium sizes, and the hierarchical variant is the best performing one for medium to large matrix sizes.

From the LU factorization without pivoting point of view, we can see in Figure 16 the behavior of both variants. As for the Cholesky-based kernels, the non-hierarchical variant's behavior was improved by using streams. However, for this kernel, the hierarchical variant outperforms the non-hierarchical one. This is mainly because even when using streams, the granularity used by the non-hierarchical variant does not allow exploiting CPU cores fully, affecting the progression on the critical path. These final results illustrate the interest in the hierarchical task paradigm from the performance point of view. The flexibility for choosing task granularities depending on the targeted resources allows one to exploit the underlying hardware's potential fully.

CONCLUSION

The growing complexity of high performance computing platforms has led to the design of highly efficient runtime systems, which aim to abstract the underlying platform. These runtime systems use various programming models to express computations. One of the most widely used models is the so-called sequential task flow model, where the task graph corresponding to the computation is submitted sequentially by a single entity. The runtime system then infers dependencies based on the data accessmodes provided by programmers at task submission. While the sequential task flow model has been shown to be effective for highly regular applications, it has limitations when computations become irregular or the task graph cannot be fully known at submission time. Additionally, the entirely sequential submission process can result in significant overhead for complex applications. FIGURE [START_REF] Allen | Optimizing Compilers for Modern Architectures: A Dependence-Based Approach[END_REF] Performance evaluation of the LU decomposition (DGETRF_nopiv) with diagonal distribution of the hierarchical tasks. The non-hierarchical variant relies on multiple streams. Each point of the plot is the best performance obtained for each variant considering tile size, number of streams, and partitioning scheme (for the hierarchical variant).

In this paper, we propose an extension of the sequential task flow model and an upgrade of the underlying runtime system to overcome these limitations. Our approach introduces a new type of task, called the hierarchical task, which can submit a new sub-graph of tasks at runtime. This allows for the possibility of parallel submission and delays the submission of parts of the task graph, enabling dynamic implementation selection. Additionally, to ensure that the parallel submission process still produces a valid directed acyclic graph (DAG), we introduce a new automatic data manager that handles data layout dynamically by submitting data management tasks at the appropriate time. We also demonstrate the correctness of our model and present a set of results that illustrate the potential of the hierarchical task model.

We plan to expand upon this work in several ways. First, we need to examine the hierarchical tasks from a scheduling perspective. We aim to answer the question of when a hierarchical task should be processed. The solution to this problem must take into account the number of tasks in the system and the work assigned to each resource. We will also consider the problem of choosing which sub-graph to submit when a hierarchical task is processed. To select the most appropriate implementation, we require advanced performance models that must be designed. Finally, the task graph resulting from a hierarchical task must be scheduled efficiently. Overall, addressing all of these issues will improve the behavior of the runtime system when using hierarchical tasks. More broadly, we aim to investigate how this model can be used to implement advanced irregular algorithms, such as linear algebra solvers based on low-rank approximation, sparse solvers and iterative solvers.

In the long-term, we believe that extending the hierarchical task model to the distributed memory context will be an elegant solution to the scalability of task-based runtime systems. 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).

FIGURE 1

 1 FIGURE 1Pseudo-code for a dummy sequential algorithm (left), corresponding STF version (center) and subsequent DAG (right).

#FIGURE 2

 2 FIGURE 2 Automatic unpartitioning example.

𝐼 ( a )

 a Root handle initialization and partition plan.

FIGURE 3

 3 FIGURE[START_REF] Buttari | A Class of Parallel Tiled Linear Algebra Algorithms for Multicore Architectures[END_REF] DAG and data layout corresponding to the code provided in Figure2. Dotted border stands for inactive, solid border stands for active. Red border stands for read-write partitioned. Green border stands for read-only partitioned or unpartitioned.

FIGURE 4

 4 FIGURE 4 Excerpt from Figure 2, rewritten to use hierarchical tasks.

𝑇

  Processing of 𝐻 1 .

𝑃 1 𝐻 2 ( b )

 2b Processing of 𝐻 1 .

𝑃 1 ( c )

 1c Processing of 𝐻 2 .

𝑃 1 𝑃 2 ( d )

 2d First task of 𝐻 1 is hierarchical and is executed.

2 ( e )

 2e Automatic insertion of the 𝑈 2 task when data needs to be accessed.

FIGURE 8

 8 FIGURE 8Example of a chain of two hierarchical tasks.

FIGURE 10

 10 FIGURE 10Submission cost of computational tasks for DGEMM with all tiles partitioned on INTEL-V100.

( a )

 a Full partitioning used in the overhead experiment. (b) 3 levels diagonal partitioning.

FIGURE 11

 11 FIGURE 11Illustration of the matrices layout used in the experiments.

  Performance results on the AMD-A100 platform.

FIGURE 14

 14 FIGURE[START_REF] Bauer | Legion: expressing locality and independence with logical regions[END_REF] Performance evaluation of the LU decomposition (DGETRF_nopiv) with diagonal distribution of the hierarchical tasks.

TABLE 1

 1 Summary of of hierachical tasks support in task-based runtime systems.

  Processing of 𝐻 1 (i.e. submission of the corresponding task graph). 𝐻 2 is not processed yet. Processing of 𝐻 2 . Dependencies between 𝐻 1 and 𝐻 2 expressed at the deepest level. Example of a DAG with 2 hierarchical tasks and 4 regular tasks.

	𝑇 2	𝑇 3	𝑇 2	𝑇 3
	1	𝑇 4	𝑇 1		𝑇 4
	𝐻 1	𝐻 2	𝑃	𝐻 2
	(a) Initial DAG. 𝑇 2 𝑇 3 (b) 𝑇 1 𝑃 𝑇 4 (c) 𝑇 1 (d) Automatic insertion of the unpartition task 𝑈 . 𝑇 2 𝑇 3 𝑃 𝑈 𝑇 4
	𝐻 1
				𝐵 11 𝐵 12
		𝐵
				𝐵 21 𝐵 22
			𝐴 11 𝐴 12	𝐶 11	𝐶 12
	𝐴	𝐶
			𝐴 21 𝐴 22	𝐶 21	𝐶 22
	(a) Initial hierarchical task and its corresponding data layout.

FIGURE 5 (b) Data partitioning and the corresponding task graph.

Table 2 .

 2

	Name	Processor	GPU	Memory
	INTEL-V100 2 x INTEL XEON GOLD 6142, 16 cores, 2.6GHz	2 x NVIDIA V100 (16GB) 384GB
	AMD-A100 2 x AMD ZEN3 EPYC 7513, 32 cores, 2.6GHz	2 x NVIDIA A100 (40GB) 512GB

TABLE 2

 2 Characteristics of the considered platforms.

TABLE 3

 3 Performances in GFlop/s used to compute 𝑅1 and 𝑅2 for each individual computational unit, and resulting percentages of tiles partitioned for the CPUs in the hierarchical version in Figure12.

DATA AVAILABILITY STATEMENT

The data that support the findings of this study are available from the corresponding author upon request.

CONFLICT OF INTEREST STATEMENT

The authors declare no potential conflict of interests.

ACKNOWLEDGEMENT

This paper is an extended version of [START_REF] Faverge | Programming Heterogeneous Architectures Using Hierarchical Tasks[END_REF] presented at the HeteroPar'22 workshop. This work is supported by the french ANR through the Solharis project under the grant (ANR-19-CE46-0009). Experiments presented in this paper were carried out using


Hachim Azzahhafi

Moussa El Yahyaoui

Ali El Moussati

Ghaïs El Zein

Ana Garcia Armada

Enhanced 4x4 MIMO RoF Architecture for 5G mmWave Indoor Applications at 60 GHz Unlicensed Band Authors

Keywords: 5G NR, RoF, MIMO, mmWave, 60 GHz, TSV channel

Multiple Input Multiple Output (MIMO) optical architecture based on Radio over Fiber (RoF) is proposed for 5G millimeter Wave (mmWave) indoor applications in the unlicensed 60 GHz band. First-order and second-order sidebands, generated by Dual-Drive Mach-Zehnder Modulator (DD-MZM), are exploited as optical subcarriers in order to reduce the number of Local Oscillators (LOs) implemented. Four MIMO-OFDM signals are multiplexed using Dual-Parallel Mach-Zehnder Modulator (DP-MZM) configured to perform Optical Carrier Suppression (OCS) modulation. 60 GHz spacing is ensured between each data signal and its corresponding subcarrier which allows simple heterodyne detection at Radio Access Point (RAP) level. In addition, we have proposed a 4x4 MIMO extension of Triple-S-and-Valenzuela (TSV) radio channel model based on geometrical approach. We demonstrate that the proposed RoF system with adequate configuration can meet the 5G New Radio (NR) requirement in terms of Error Vector Magnitude (EVM). Moreover, Bit Error Rate (BER) performance is evaluated with and without 5G NR channel coding, for transmission over 25 km of optical fiber followed by 3 m of radio channel, by exploiting the overall allocated band at 60 GHz which allows to achieve a very high data rate up to 122.5 Gb/s.

Introduction

The 5th generation (5G) has changed the concepts of wireless mobile communication systems that were focused only on increasing data rate and providing services for personal use. The 5G technology is developed to respond to the huge demand for new applications and services for both personal and professional uses, such as the Internet of Things (IoT), Augmented Reality (AR), Cloud Computing and Gaming, Smart Building, Industry Automation and Mission Critical Applications [START_REF] Husenovic | ITU, Setting the scene for 5G: Opportunities & Challenges[END_REF]. To meet these diverse challenges, 5G systems are expected to support three typical use cases, namely enhanced Mobile Broadband (eMBB), massive Machine Type Communications (mMTC) and Ultra Reliable Low Latency Communications (URLLC) [START_REF] Husenovic | ITU, Setting the scene for 5G: Opportunities & Challenges[END_REF].

Key enabling radio access technologies such as millimeter Wave (mmWave) communication, Multiple Input Multiple Output (MIMO) technology, Orthogonal Frequency Division Multiplexing (OFDM) numerology and Low-Density Parity-Check (LDPC) codes have been specified by 3rd Generation Partnership Project (3GPP), the organization in charge of standardization and specification of the 5G, under the name of New Radio (NR) technologies [START_REF]3GPP, Technical Specification Group Services and System Aspects[END_REF][START_REF] Wong | Key Technologies for 5G Wireless Systems[END_REF]. The main aim to adopt frequencies in the mmWave bands (30-300 GHz), is not only to mitigate the burden on the sub-6 GHz frequency band, but also to exploit its large bandwidth which directly leads to increase the data rate and reduce latency with higher numerologies. Moreover, the propagation loss in the mmWave, which was seen previously as a drawback, is considered now the main feature to support short-range applications in ultra-dense and indoor environments [START_REF] Gao | MmWave massive-MIMO-based wireless backhaul for the 5G ultra-dense network[END_REF][START_REF] Kim | MIMO-Supporting Radio-Over-Fiber System and its Application in mmWave-Based Indoor 5G Mobile Network[END_REF].

Currently, the licensed 28 GHz band is the adopted mmWave band for deploying 5G NR technologies. However, to ensure more capacity and flexibility in terms of technologies and applications, the 3GPP has started scaling 5G NR design to support unlicensed bands, specifically the 60 GHz band. The specification of this band, known as NR Unlicensed (NR-U), has been recently included as a part of release 17 [START_REF]3GPP, Technical Specification Group Services and System Aspects[END_REF]. The 60 GHz band, exploited by well-known wireless standards such as IEEE 802.15.3c [START_REF] Baykas | IEEE 802.15.3c: the first IEEE wireless standard for data rates over 1 Gb/s[END_REF], IEEE 802.11ad [START_REF] Perahia | IEEE 802.11ad: Defining the Next Generation Multi-Gbps Wi-Fi[END_REF], and IEEE 802.11ay [START_REF] Ghasempour | IEEE 802.11ay: Next-Generation 60 GHz Communication for 100 Gb/s Wi-Fi[END_REF], offers four channels with high bandwidths up to 8.64 GHz [START_REF] Ghasempour | IEEE 802.11ay: Next-Generation 60 GHz Communication for 100 Gb/s Wi-Fi[END_REF]. Several works have dealt with the coexistence of 5G NR-U and IEEE 802.11 technologies [START_REF] Lagen | New Radio Beam-Based Access to Unlicensed Spectrum: Design Challenges and Solutions[END_REF][START_REF] Hirzallah | 5G New Radio Unlicensed: Challenges and Evaluation[END_REF][START_REF] Daraseliya | Coexistence analysis of 5g nr unlicensed and wigig in millimeter-wave spectrum[END_REF][START_REF] Lu | 5G-U: Conceptualizing Integrated Utilization of Licensed and Unlicensed Spectrum for Future IoT[END_REF]. The fact that the radio coverage of 60 GHz band is limited to a few meters and waves at this frequency cannot penetrate walls makes it more suitable for 5G In-Building applications. Moreover, to ensure 5G connectivity in the entire building, the radio coverage needs to be expanded. The most viable solution to reach this aim is based on using a network of Radio Access Points (RAPs) supported by Radio over Fiber (RoF) technology. By deploying this solution, a building can be covered by multiple small cells served by low-cost RAPs. These RAPs are linked, using RoF technology, to a centralized Base Station (BS) where high-complexity signal processing is performed.

RoF technology has attracted attention in recent years because of its many benefits such as the large bandwidth and low attenuation offered by the optical fiber. Moreover, the RAPs are simplified thanks to centralized processing at BS. Another advantage of RoF is the possibility to exploit the existing optical fiber networks [START_REF] Armada | Radio-over-Fiber Aided Base Station Coordination for OFDM[END_REF]. To distribute signals between RAPs and BS, there are three different categories of RoF, namely Base Band over Fiber (BBoF), Intermediate Frequency-over-Fiber (IFoF), and Radio Frequency-over-Fiber (RFoF) [START_REF] Beas | Millimeter-Wave Frequency Radio over Fiber Systems: A Survey[END_REF]. The latter is the most suitable in our case since the 60 GHz signals are processed at BS and transported via an optical fiber to RAPs where the 60 GHz signals are photo-detected and directly transmitted on the air interface using antennas. By adopting this configuration, the RAPs complexity and cost are reduced. In addition, the fact that all complex treatments are carried out at the centralized BS makes this system easy to upgrade, which is a fundamental feature of 5G systems [START_REF] Akshatha | A Centralized SDN Architecture for the 5G Cellular Network[END_REF].

To guarantee high capacity and reliable communications, RoF systems can be endowed with MIMO technology. The straightforward solution to transport independent MIMO signals is to use as many optical links as there are MIMO channels, which leads to an increase in deployment costs. Several works have dealt with this challenge by using different multiplexing schemes to transport mmWave MIMO signals over a single optical link, including Wavelength Division Multiplexing (WDM) [START_REF] Capmany | Multiwavelength single sideband modulation for WDM radio-over-fiber systems using a fiber grating array tandem device[END_REF], Polarization Division Multiplexing (PDM) [START_REF] Huang | Direct-Detection PDM-OFDM RoF System for 60-GHz 2×2 MIMO Wireless Transmission Without Polarization Tracking[END_REF] and Optical Subcarriers Multiplexing (OSCM) [START_REF] -T. Lin | 60-GHz optical/wireless MIMO system integrated with optical subcarrier multiplexing and 2x2 wireless communication[END_REF].

In the WDM scheme, each MIMO channel requires an optical source and optical modulator. Moreover, high frequency Local Oscillators (LOs) are required to transport mmWave signals using this scheme. The complexity of this scheme increases by increasing the number of antennas in the MIMO system. In [START_REF] Huang | Direct-Detection PDM-OFDM RoF System for 60-GHz 2×2 MIMO Wireless Transmission Without Polarization Tracking[END_REF], a PDM scheme has been proposed to transport a 2x2 MIMO signal at 60 GHz, where a dual tone laser source is used to feed two independent Dual-Drive Mach-Zehnder Modulators (DD-MZMs). In addition, two LOs are implemented to perform frequency translation of OFDM signals. The authors in [START_REF] -T. Lin | 60-GHz optical/wireless MIMO system integrated with optical subcarrier multiplexing and 2x2 wireless communication[END_REF] used OSCM to transport 2x2 MIMO over Single Mode Fiber (SMF) by using Dual-Parallel Mach-Zehnder Modulators (DP-MZM) and two LOs. However, these architectures are limited to transport only 2x2 MIMO signals. The 4x4 MIMO version of these systems requires using two optical links or a combination with WDM scheme. These solutions increase the complexity of the system and are not cost-effective. In our previous works [START_REF] Yahyaoui | On the capacity of MIMO-OFDM based diversity and spatial multiplexing in Radio-over-Fiber system[END_REF][START_REF] Yahyaoui | New millimeter wave generation scheme for MIMO-OFDM based Radio-over-Fiber system[END_REF], we have introduced two optical schemes to transport 2x2 MIMO signals at 60 GHz, where two MIMO-OFDM signals combined with two optical subcarriers are generated and transmitted over optical fiber. These proposed RoF systems achieve a data rate up to 70 Gb/s.

In this work, we further improve our system by proposing a new 4x4 MIMO scheme which allows to reach a higher data rate up to 122 Gb/s. This proposal uses one optical link fed by one laser source and only three LOs. Frequency-Division Multiplexing (FDM) is performed to combine four OFDM signals by using two LOs. These signals are fed to a DP-MZM biased to perform Optical Carrier Suppression (OCS) modulation which leads to generate four independent optical MIMO signals. Our architecture allows the use of simple heterodyne detection using Photodiodes (PDs) to generate 60 GHz signals to be transmitted over the radio channel. For this purpose, four optical subcarriers, corresponding to the first and the second order optical sidebands, are generated by a DD-MZM driven by the third LO. This optical modulator is configured to perform OCS modulation with the third-order sidebands suppression. The use of both the first and the second order sidebands as optical subcarriers in this work, is supported by the fact that the Signal to Noise Ratio (SNR) does not depend on the power level of the optical subcarriers in the heterodyne detection technique at receiver [START_REF] Witteman | Detection and Signal Processing[END_REF]. The frequencies delivered by the LOs are well chosen to guarantee 60 GHz between each OFDM signal/subcarrier pair. The optical subcarriers are combined with the optical OFDM-MIMO signals by using an optical coupler. This resulting optical signal is transmitted over fiber to the RAP.

Moreover, to evaluate the performance of our RoF system for indoor applications, we introduce a 4x4 MIMO extension of the Triple-S and Valenzuela (TSV) channel model based on geometrical approach. This model was first proposed as a realistic Single Input Single Output (SISO) channel to model indoor environments for the 802.15.3c standard [START_REF] Shoji | A Modified SV-Model Suitable for Line-of-Sight Desktop Usage of Millimeter-Wave WPAN Systems[END_REF]. For performance evaluation, channel coding schemes are implemented in this work, according to the 5G NR specifications for multiplexing and channel coding [START_REF] Gpp ; G ; Nr | Multiplexing and channel coding Release 16[END_REF].

The rest of this paper is organized as follows: Section 2 presents the proposed 4x4 MIMO optical architecture. Section 3 introduces the 4x4 MIMO extension of the TSV channel model. In Section 4, we present the obtained results and discussions. Finally, conclusions are given in Section 5.

Proposed 4x4 MIMO optical architecture

The proposed 4x4 MIMO RoF architecture transmitting signals from the BS to the RAP is presented in Error! Reference source not found., where four baseband complex signals, obtained by performing OFDM modulation are converted to electrical signals by using I/Q mixers. The first two obtained signals, corresponding to the transmit antennas Tx1 and Tx2, are up-converted in frequency using two local oscillators LO1 and LO2 of frequencies 𝑓 LO1 = 12 GHz and 𝑓 LO2 = 36 GHz, respectively. These signals are then frequency multiplexed as shown in inset (a). The same processing is performed on the second two signals, corresponding to the transmit antennas Tx3 and Tx4, the corresponding multiplexed signal is shown in inset (a'). The two resulting multiplexed signals are coupled through a 90° hybrid coupler. The phase-shifted outputs of the hybrid coupler, shown in insets (b) and (b'), drive the optical modulator based on DP-MZM. This component is described in detail in subsection 2.1.

The modulated optical MIMO signals at the output of DP-MZM, shown in inset (c), are coupled with four optical subcarriers generated using DD-MZM modulator by exploiting first-order and second-order sidebands. This modulator is driven by a signal delivered from a third local oscillator at 𝑓 LO3 = 24 GHz. The DD-MZM configuration is mathematically demonstrated in subsection 2.2. Subsequently, the resulting optical MIMO signal, corresponding to inset (e), is transmitted through the optical fiber to the RAP.

The received optical signal is interleaved to separate the four optical components of the 4x4 MIMO signal. Indeed, the frequencies of the local oscillators employed in this design were chosen to be multiples of 12 GHz, allowing us to separate all components using a demultiplexer with a frequency spacing of 12 GHz. The channels are then merged with their respective subcarriers to achieve 60 GHz spacing. Finally, four PDs are placed to perform heterodyne detection to generate four electrical signals at 𝑓 RF = 60 GHz. These signals are then filtered, amplified, and transmitted over the MIMO wireless channel, as shown in inset (f).

Data modulation using DP-MZM

Our proposal is based on OCS for data modulation using DP-MZM. To perform this kind of modulation, a continuous optical carrier generated at 𝑓 𝑐 , by the Laser Diode (LD), is fed to the DP-MZM, which is biased by using appropriate parameters to realize OCS modulation. The DP-MZM used in this work is based on [START_REF] Yang | A general analytical method for suppressing the third-order intermodulation in microwave photonic link based on dual-parallel Mach-Zehnder modulator[END_REF]. This modulator consists of three sub-MZMs which can be controlled separately to change their working points. A detailed schematic of the DP-MZM with OCS configuration is shown in Fig. 2.

The first output of the hybrid coupler, shown in inset (a) of Fig. 2 superposed in a manner to have constructive and destructive interference allowing the separation of the overlapped components and creating the signal shown in inset (d) of Fig. 2. The main parameters of the DP-MZM are shown in Table 1.

Optical subcarriers generation

Our proposal is based on optical heterodyne detection at the RAP to generate four electrical signals to be transmitted at 𝑓 RF = 60 GHz using MIMO antennas as depicted in Error! Reference source not found.. In addition to the optical data-signals generated by DP-MZM, the optical heterodyne technique needs to have four optical subcarriers positioned at well-chosen frequencies to be detected by the PDs. To achieve that, we have opted for an optical subcarrier generation using DD-MZM powered by the LD and driven by the electrical output of the local oscillator at 𝑓 𝐿𝑂3 = 24 GHz. Two optical subcarriers are positioned at 𝑓 𝑐 -𝑓 𝐿𝑂3 and 𝑓 𝑐 + 𝑓 𝐿𝑂3 which represent the first-order sidebands and the other subcarriers are positioned at 𝑓 𝑐 -2𝑓 𝐿𝑂3 and 𝑓 𝑐 + 2𝑓 𝐿𝑂3 which represent the second-order sidebands. This electrical signal is expressed as

𝑉 𝐿𝑂3 (𝑡) = 𝑉 𝑚 cos(2𝜋𝑓 𝐿𝑂3 𝑡) (1)
where 𝑉 𝑚 is the peak amplitude of the signal.

The optical field of the signal generated by the DD-MZM is expressed as follows [START_REF] Armada | Radio-over-Fiber Aided Base Station Coordination for OFDM[END_REF] 𝐸 ] 𝑒 𝑗2𝜋𝑓 𝑐 𝑡

(2)

where 𝐸 𝐿𝑎𝑠𝑒𝑟 (𝑡), 𝑃 𝐿𝑎𝑠𝑒𝑟 , and 𝑓 𝑐 = 𝜔 𝑐 2𝜋 are the optical field, the transmit power and the optical carrier frequency of the LD, respectively. 𝑡 𝑎𝑡𝑡𝑛 is the insertion loss of the modulator and 𝑉 𝜋 𝐷𝐷 is the switching voltage of the DD-MZM, i.e., the voltage which results in a phase change of 𝜋 radians in the signal propagating in the arm where it is applied. 𝑉 1 (𝑡) and 𝑉 2 (𝑡) are the driving signals applied in the first arm and in the second arm of the DD-MZM respectively, and are expressed as

𝑉 1 (𝑡) = 𝑉 𝑚 cos(2𝜋𝑓 𝐿𝑂3 𝑡) + 𝑉 𝑏𝑖𝑎𝑠 𝐷𝐷 (3) 𝑉 2 (𝑡) = 𝑉 𝑚 cos(2𝜋𝑓 𝐿𝑂3 𝑡 + 𝜃) (4)
where 𝑉 𝑏𝑖𝑎𝑠 𝐷𝐷 is the bias voltage added to 𝑉 𝐿𝑂3 (𝑡) in the first arm of the DD-MZM (the bias voltage of the second arm is fixed to 0 volt). 𝜃 is the phase of the electrical phase shifter placed at the input of the second arm as depicted in Error! Reference source not found.. ] 𝑒 𝑗2𝜋𝑓 𝑐 𝑡 (5)

By inserting

For simplicity, Eq. ( 5) can be rewritten as

𝐸 𝐷𝐷-𝑀𝑍𝑀 (𝑡) = 𝐸 𝑎𝑡𝑡𝑛 [𝑒 (𝑗𝑚cos(2𝜋𝑓 𝐿𝑂3 𝑡)+𝑗𝜋𝛾) + 𝑒 𝑗𝑚cos(2𝜋𝑓 𝐿𝑂3 𝑡+𝜃) ] 𝑒 𝑗2𝜋𝑓 𝑐 𝑡
= 𝐸 𝑎𝑡𝑡𝑛 [𝑒 𝑗𝑚cos(2𝜋𝑓 𝐿𝑂3 𝑡) 𝑒 𝑗𝜋𝛾 + 𝑒 𝑗𝑚cos(2𝜋𝑓 𝐿𝑂3 𝑡+𝜃) ] 𝑒 𝑗2𝜋𝑓 𝑐 𝑡 [START_REF]3GPP, Technical Specification Group Services and System Aspects[END_REF] where 𝐸 𝑎𝑡𝑡𝑛 = √ 𝑃 𝐿𝑎𝑠𝑒𝑟 𝑡 𝑎𝑡𝑡𝑛 2

, 𝑚 = 𝜋 𝑉 𝑚 𝑉 𝜋 𝐷𝐷 is the modulation index of the DD-MZM and 𝛾 = 𝑉 𝑏𝑖𝑎𝑠 𝐷𝐷 𝑉 𝜋 𝐷𝐷 is its normalized bias voltage. From Jacobi-Anger expansion [START_REF] Ujang | Irregular Shifting of RF Driving Signal Phase to Overcome Dispersion Power Fading[END_REF], Eq. ( 6) becomes

𝐸 𝐷𝐷-𝑀𝑍𝑀 (𝑡) = 𝐸 𝑎𝑡𝑡𝑛 ∑ +∞ 𝑛=-∞ 𝑗 𝑛 𝐽 𝑛 (𝑚)(𝑒 𝑗𝑛(2𝜋𝑓 𝐿𝑂3 𝑡) 𝑒 𝑗𝜋𝛾 + 𝑒 𝑗𝑛(2𝜋𝑓 𝐿𝑂3 𝑡) 𝑒 𝑗𝑛𝜃 ) 𝑒 𝑗2𝜋𝑓 𝑐 𝑡 = 𝐸 𝑎𝑡𝑡𝑛 ∑ +∞ 𝑛=-∞
𝑗 𝑛 𝐽 𝑛 (𝑚)(𝑒 𝑗𝜋𝛾 + 𝑒 𝑗𝑛𝜃 ) 𝑒 𝑗2𝜋(𝑓 𝑐 +𝑛𝑓 𝐿𝑂3 )𝑡 [START_REF] Baykas | IEEE 802.15.3c: the first IEEE wireless standard for data rates over 1 Gb/s[END_REF] where 𝐽 𝑛 (𝑚) is the 𝑛 𝑡ℎ Bessel function of the first kind. By observing Eq. ( 7), the optical field generated from the DD-MZM consists of an optical carrier and sidebands with infinite order. However, the order of most important sidebands, in terms of power, can be adjusted by changing the amplitude of the signal emitted by LO and consequently the modulation index m. For our case, the sidebands of 0 𝑡ℎ , 1 𝑠𝑡 , 2 𝑛𝑑 , and 3 𝑟𝑑 order are the most important, while the sidebands of order 𝑛 ≥ 4 are neglected. Then we list the optical fields of these sidebands 𝐸 -3 = 𝐸 𝑎𝑡𝑡𝑛 . 𝐽 -3 (𝑚). (𝑒 𝑗𝜋𝛾 + 𝑒 -𝑗3𝜃 ). 𝑒 𝑗2𝜋(𝑓 𝑐 -3𝑓 𝐿𝑂3 )𝑡 . 𝑒 𝑗 𝜋 2 𝐸 -2 = 𝐸 𝑎𝑡𝑡𝑛 . 𝐽 -2 (𝑚). (𝑒 𝑗𝜋𝛾 + 𝑒 -𝑗2𝜃 ). 𝑒 𝑗2𝜋(𝑓 𝑐 -2𝑓 𝐿𝑂3 )𝑡 . 𝑒 𝑗𝜋

J o u r n a l P r e -p r o o f

Journal Pre-proof In order to realize the OCS modulation combined with the third-order sidebands suppression, the conditions below should be verified

𝑒 𝑗𝜋𝛾 + 1 = 0 𝑒 𝑗𝜋𝛾 + 𝑒 -𝑗3𝜃 = 0 (9) 𝑒 𝑗𝜋𝛾 + 𝑒 𝑗3𝜃 = 0 𝑒 𝑗𝜋𝛾 + 𝑒 𝑗𝑛𝜃 ≠ 0 with 𝑛 = 1; 2.
By solving this system of equations, we have

𝛾 = 𝑉 𝑏𝑖𝑎𝑠 𝐷𝐷-𝑀𝑍𝑀 𝑉 𝜋 𝐷𝐷 = 2𝑘 + 1 with 𝑘 ∈ ℤ and 𝜃 = 2𝜋 3 (𝑚𝑜𝑑 2𝜋) (10)
Therefore, the DD-MZM must be biased at 𝑉 𝑏𝑖𝑎𝑠 𝐷𝐷 = 𝑉 𝜋 𝐷𝐷 and the phase shifter placed must be configured to have 𝜃 = 120°.

4x4 MIMO extension of the TSV channel model

TSV channel model

We consider the TSV channel model, proposed by the National Institute of Information and Communication Technology (NICT) [START_REF] Sawada | Proposal of novel statistic channel model for millimeter wave WPAN[END_REF]. This model was used, in the framework of IEEE 802.15.3c standardization, as a realistic indoor channel model for mmWave band at 60 GHz. Indeed, the TSV model is based on statistical models validated by several experimental results and analysis for different types of environments, such as Office, Library, Kiosk and Desktop, for both Line of Sight (LoS) and non-LoS (NLoS) links.

In this work, we are interested in the LoS desktop environment. For this environment, the transmission channel behavior can be modeled by merging two well-known channel models. For LoS path components, the two-path model is dominantly observed, as for the NLoS path components, the SV model is adopted [START_REF] Shoji | A Modified SV-Model Suitable for Line-of-Sight Desktop Usage of Millimeter-Wave WPAN Systems[END_REF]. Then, the Channel Impulse Response (CIR) of the TSV model is expressed in function of the CIR of the two-path model, ℎ two (𝑡), and the CIR of the SV model ℎ sv (𝑡) as follows ℎ tsv (𝑡, 𝜑) = ℎ two (𝑡) + ℎ sv (𝑡, 𝜑) [START_REF] Hirzallah | 5G New Radio Unlicensed: Challenges and Evaluation[END_REF] where ℎ two (𝑡) is expressed as follows 1) 𝐺 𝑅𝑥 (1) + √ 𝐺 𝑇𝑥 (2) 𝐺 𝑅𝑥 (2)

ℎ two (𝑡) = 𝛽𝛿(𝑡) = ( 𝜆 4𝜋𝐷 ) 2 | √ 𝐺 𝑇𝑥 (
𝑅 0 exp [𝑗 2𝜋 𝜆 2ℎ 𝑇𝑥 ℎ 𝑅𝑥 𝐷 ]| 𝛿(𝑡) (12)
where 𝛿(. ), 𝛽, 𝜆, ℎ 𝑇𝑥 , ℎ 𝑅𝑥 , 𝐷 and 𝑅 0 denote the Dirac delta function, the path gain for LOS component, the wavelength, the height of Tx antenna, that of Rx antenna, the Tx-Rx distance separation and the reflection coefficient of the desk, respectively. 𝐺 𝑇𝑥 (1) , 𝐺 𝑇𝑥 (2) , 𝐺 𝑅𝑥 (1) and 𝐺 𝑅𝑥 (2) are the Tx antenna gain for the direction of direct path direction (path 1), the Tx antenna gain for the J o u r n a l P r e -p r o o f reflected path direction (path 2), the Rx antenna gain for the direct path direction, and the Rx antenna gain for the reflected path direction, respectively.

Then, ℎ sv (𝑡) is expressed as a function of time-invariant random parameters as follows

ℎ sv (𝑡, 𝜑) = ∑ 𝐿 𝑙=1 ∑ 𝑀 𝑙 𝑚=1
𝛼 𝑙𝑚 𝛿(𝑡 -𝑇 𝑙 -𝜏 𝑙𝑚 )𝛿(𝜑 -𝛹 𝑙 -𝜓 𝑙𝑚 ) [START_REF] Lu | 5G-U: Conceptualizing Integrated Utilization of Licensed and Unlicensed Spectrum for Future IoT[END_REF] where 𝐿, 𝑀 𝑙 , 𝑇 𝑙 , 𝜏 𝑙𝑚 , 𝛹 𝑙 and 𝜓 𝑙𝑚 denote the number of clusters, the number of rays in the 𝑙 th cluster, the arrival time of the 𝑙 th cluster, the arrival time of the 𝑚 th ray in the 𝑙 th cluster, the angle of arrival of the 𝑙 th cluster and that of the 𝑚 th ray in the 𝑙 th cluster, respectively. 𝛼 𝑙𝑚 is a statistical complex-envelope gain of the m th ray in the l th cluster, where the amplitude follows a log-normal distribution, and the phase follows a uniform distribution [START_REF] Shoji | A Modified SV-Model Suitable for Line-of-Sight Desktop Usage of Millimeter-Wave WPAN Systems[END_REF].

MIMO Extension

To extend the TSV channel model to support MIMO scenarios, a geometrical approach is considered, as shown in Fig. 3. It can be noted that for a transmission in mmWave band, the effect of the distance between different Tx and Rx antennas is significant. For simplicity, we consider the uniform linear antenna arrays, where the Tx and Rx antennas are distributed vertically. For 𝑁 𝑅 × 𝑁 𝑇 MIMO systems, the received signal 𝑌 is described by

𝑌 = 𝑯 𝑋 + 𝑛 ( 14
)
where 𝑋, 𝑛, 𝑁 𝑇 and 𝑁 𝑅 are the transmitted signal, the additive channel noise, the number of transmit antennas and the number of receive antennas, respectively. 𝑯 is the 𝑁 𝑅 × 𝑁 𝑇 channel matrix represented by

𝑯 = ( ℎ 1,1 ℎ 1,2 ⋯ ℎ 1,𝑁 𝑇 ℎ 2,1 ℎ 2,2 ⋯ ℎ 2,𝑁 𝑇 ⋮ ⋮ ⋱ ⋮ ℎ 𝑁 𝑅 ,1 ℎ 𝑁 𝑅 ,2 ⋯ ℎ 𝑁 𝑅 ,𝑁 𝑇 ) (15)
Each element ℎ 𝑖,𝑘 (𝑖 ≤ 𝑁 𝑅 and k ≤ 𝑁 𝑇 ) of the matrix H can be expressed according to the TSV channel model. For LoS components, the gain 𝛽 𝑖,𝑘 between the 𝑘 𝑡ℎ transmit antenna and the 𝑖 𝑡ℎ receive antenna is given by

𝛽 𝑖,𝑘 = ( 𝜆 4𝜋𝐷 ) 2
| √ 𝐺 𝑇𝑥 (1) 𝐺 𝑅𝑥 (1) exp [-𝑗 2𝜋 (𝑑 𝑖,𝑘 (1) -𝑑 1,1

+ √ 𝐺 𝑇𝑥 (2) 𝐺 𝑅𝑥 (2)

𝑅 0 exp [𝑗 2𝜋 𝜆 2ℎ 𝑇𝑥 ℎ 𝑅𝑥 𝐷 ] exp [-𝑗 2𝜋 (𝑑 𝑖,𝑘 (2) -𝑑 1,1 (2) ) 𝜆 ]| (16)
where for p = 1,2,

[𝑑 𝑖,𝑘 (𝑝) ] 2 = [𝑑 1,1 (𝑝) -(𝑖 -1)Δ 𝑅𝑥 𝑠𝑖𝑛 𝜃 𝑅𝑥1 (𝑝) + (𝑘 -1)Δ 𝑇𝑥 𝑠𝑖𝑛 𝜃 𝑇𝑥1 (𝑝) ] 2 + [(𝑘 -1)𝛥 𝑇𝑥 𝑐𝑜𝑠 𝜃 𝑇𝑥1 (𝑝) + (𝑖 -1)𝛥 𝑅𝑥 𝑐𝑜𝑠 𝜃 𝑅𝑥1 (𝑝) ] 2 (17)
where 𝑑 𝑖,𝑘 (𝑝) is the distance between the 𝑘 𝑡ℎ transmit antenna and the 𝑖 𝑡ℎ receive antenna along path 𝑝. Δ 𝑇𝑥 , Δ 𝑅𝑥 , 𝜃 𝑇𝑥1 (𝑝) and 𝜃 𝑅𝑥1 (𝑝) are the Tx antenna separation, the Rx antenna separation, the angle of departure (AoD) of the path p and the angle of arrival (AoA) of the path p, respectively. For the NLoS path components, each element of matrix H can be calculated using SV channel model in the same way as in a SISO scenario.

Results and discussion

A global simulation of our proposed system is carried out using MATLAB/OptiSystem cosimulation platform proposed and used in our previous works [START_REF] Yahyaoui | On the capacity of MIMO-OFDM based diversity and spatial multiplexing in Radio-over-Fiber system[END_REF][START_REF] Yahyaoui | New millimeter wave generation scheme for MIMO-OFDM based Radio-over-Fiber system[END_REF]. The electrical upconversion and the optical link are implemented in OptiSystem, while all other processes are performed using MATLAB. Data bits are generated and encoded according to the 5G NR specifications for multiplexing and channel coding [START_REF] Gpp ; G ; Nr | Multiplexing and channel coding Release 16[END_REF]. The encoding processing includes cyclic redundancy check (CRC) attachment, code block segmentation, LDPC coding and rate matching. The encoded bits are mapped using QPSK, 16-QAM, and 64-QAM. The complex symbols are fed to the OFDM modulator to generate four baseband MIMO signals. Channel estimation pilots and Cyclic Prefix (CP) are added to the signals, then a Digital-to-Analogue Conversion (DAC) is performed. We have used the parameters defined in Table 2 to generate the OFDM MIMO signal, which occupies the 7,2 GHz bandwidth. The resulting signals fed our optical architecture specified above. The main parameters of this architecture used for simulation are shown in Table 3. The 60 GHz signals at the output of PDs are transmitted over the 4x4 MIMO channel. To achieve this, a channel matrix is generated according to the 4x4 MIMO extension of the TSV model. At the receiver, the inverse processes are done. Channel estimation at pilot positions with arithmetic averaging is implemented according to the Least Squares (LS) criterion. The obtained channel coefficients are then interpolated to obtain the estimated channel coefficients at data positions. To reduce MIMO channel effects, Zero-Forcing Successive Interference Cancellation (ZF-SIC) detector is implemented. Firstly, we have configured the optical architecture without Additive White Gaussian Noise (AWGN) and without the TSV MIMO channel. The choice of the transmit laser power 𝑃 𝐿𝑎𝑠𝑒𝑟 has been found by simulation as illustrated in Fig. 4, where the Error Vector Magnitude (EVM) of the received 16-QAM OFDM signals is calculated for 25 km of optical fiber length for different values of transmitting laser power injected in the optical modulators. We can see that a value of 𝑃 𝐿𝑎𝑠𝑒𝑟 = 7 dBm offers the best performance in terms of reducing the degradation caused by the nonlinearity of the optical component and optical fiber. We have demonstrated in Fig. 6, the effect of the optical fiber length on the EVM performance of the optical link, using only the AWGN channel for different values of energy per bit to noise power spectral density ratio (Eb/No). The obtained values are compared with the EVM threshold of 12.5% defined by 3GPP requirements for 16-QAM modulation [START_REF] Gpp ; G; Nr | User Equipment (UE) radio transmission and reception; Part 1: Range 1 Standalone[END_REF]. This required level is achieved for an optical fiber with a length ranging from 5 to 25 km and an Eb/No ratio of 5 dB.

In order to add the effect of the MIMO radio channel, we have generated the channel matrix for different Δ 𝑇𝑥 and Δ 𝑅𝑥 antenna separations. Then, we evaluated the performance of the 16-QAM 4x4 MIMO system, in terms of EVM versus antenna separation, without using the optical link, as shown in Fig. 7. As we can see, EVM performance changes quickly by changing antenna separation. This variation is due to the significant phase rotation created by the difference in distance between the two signal paths at 60 GHz. To guarantee sufficient field decorrelation in MIMO antennas, we have chosen the antennas separation of Δ 𝑇𝑥 = 8 mm and Δ 𝑅𝑥 = 8 mm, which correspond here to the minimum value of EVM. The chosen value of the antennas separation of 8 mm represents 1,6 times the wavelength of the RF carrier at 60 GHz. To show the effect of radio channel length on performance, we have calculated the EVM in function of the channel length, as shown in Fig. 8. A strong degradation is noticed from the length of 6 m. This degradation is due to high level of attenuation present in the desktop environment of TSV channel at 60 GHz. In the rest of this paper, we have chosen to use the radio channel with length of 3 m. This distance corresponds to the average size of a small office room. Also, it is considered in several works and validated by experimental results as indicated in reference [START_REF] Shoji | A Modified SV-Model Suitable for Line-of-Sight Desktop Usage of Millimeter-Wave WPAN Systems[END_REF]. J o u r n a l P r e -p r o o f

Next, we have evaluated the performance of our global RoF system that includes both the effect of the optical fiber of 25 km and the radio MIMO channel. Fig. 9 To improve the quality of transmission and to meet 5G NR specifications, we have implemented the channel coding blocks. Fig. 10 presents the BER performance of 64-QAM using CRC, LDPC and rate matching for different coding rates. As shown, a significant improvement in terms of BER is achieved comparing with the uncoded 64-QAM. A BER of 10 -6 can be reached with Eb/No of 10 dB using coding rate of 1/2. In addition, using these coding blocks can improve the data rate of the global system. Indeed, the data rate of transmitting 64-QAM signal with coding rate of 1/2 which corresponds to 70 Gb/s, exceeds the data rate of the uncoded QPSK signal with a BER performance enhancement. The same observation is noted, compared to uncoded 16-QAM, for 64-QAM with coding rates of 3/4 and 7/8, which correspond to data rates of 105 Gb/s and 122.5 Gb/s, respectively.

Conclusion

In this paper, we presented a new 4x4 MIMO RoF architecture for 5G indoor mmWave applications. The proposal allows the transmission of 60 GHz OFDM MIMO signal over the optical fiber and performs heterodyne detection at the RAP level to reduce its cost and complexity. The optical architecture is based on DP-MZM for data optical modulation and DD-14 MZM for optical subcarrier generation by exploiting its optical sidebands. To evaluate this scheme under a realistic indoor scenario, we have proposed a 4x4 MIMO extension of the wellknown TSV channel model. For performance evaluation purposes, we have configured the optical power injected in the optical modulators and the MIMO antennas separation based on simulation results. The global RoF system is investigated in terms of BER versus of Eb/No for uncoded QPSK, 16-QAM, and 64-QAM schemes. Good performance is observed for uncoded QPSK and 16-QAM. To meet 5G NR specifications, we have implemented the channel coding techniques proposed by 3GPP. Significant improvement in terms of BER of 64-QAM is achieved using different coding rates. Also, data rate enhancement up to 122.5 Gb/s can be noticed by comparing coded 64-QAM with uncoded QPSK and uncoded 16-QAM under enhanced BER performance.

  drives differentially the first sub-MZM biased at its switching voltage 𝑉 𝑏𝑖𝑎𝑠1 = 𝑉 𝜋 𝐷𝑃 which leads to perform double-side band (DBS) modulation and generate the first DBS signal shown in inset (b) of Fig. 2, whereas the second output, in inset (a') of Fig. 2, drives the second sub-MZM using the same configuration i.e., 𝑉 𝑏𝑖𝑎𝑠2 = 𝑉 𝜋 𝐷𝑃 and a second DBS signal shown in inset (b') of Fig. 2 is generated. The third sub-MZM, which acts as an optical phase shifter, is biased at 𝑉 𝑏𝑖𝑎𝑠3 = 𝑉 𝜋 𝐷𝑃 2 to have a 𝜋 2 phase difference between the two DSB Signals. By using this configuration, these DBS signals are

Fig. 1 .

 1 Fig. 1. The proposed optical architecture.

Fig. 2 .

 2 Fig. 2. Detailed schematic of the DP-MZM with OCS configuration.

Fig. 3 .

 3 Fig. 3. MIMO extension of the TSV channel Model.

Fig. 5 (

 5 Fig.5(a) shows the power spectrum of the optical MIMO signals modulated using the DP-MZM. As we can see, the four signals are well placed at frequencies 𝑓 𝑐 -𝑓 𝐿𝑂2 , 𝑓 𝑐 -𝑓 𝐿𝑂1 , 𝑓 𝑐 + 𝑓 𝐿𝑂1 , and 𝑓 𝑐 + 𝑓 𝐿𝑂2 . Fig.5(b) shows the optical subcarriers, that correspond to the 1 𝑠𝑡 and 2 𝑛𝑑 sidebands, generated by the DD-MZM. As shown, the optical carrier of frequency 𝑓 𝑐 = 193.41 THz is suppressed. The power spectrum of the combined signals to be transmitted over the optical fiber is shown in Fig.5 (c). It can be seen that the spacing of 60 GHz is guaranteed between each datasignal spectrum and its corresponding subcarrier.

Fig. 4 .

 4 Fig. 4. Effect of the transmit laser power on the EVM of 16-QAM OFDM signals for 25 km of optical fiber length.

Fig. 5 .

 5 Fig. 5. (a) The optical spectrum of data-signals generated by DP-MZM, (b) The optical subcarriers generated by the DD-MZM, (c) The combined signals to be transmitted over fiber.

Fig. 7 .

 7 Fig. 7. EVM performance for different MIMO antenna separations.

Fig. 8 .

 8 Fig. 8. EVM in function of radio channel length at Eb/No of 20 dB.

Fig. 6 .

 6 Fig.6. Effect of the optical length on the EVM of the system using only the AWGN channel.

Fig. 9 .

 9 Fig. 9. BER performance of the global RoF system using uncoded QPSK, 16-QAM and 64-QAM.

Fig. 10 .

 10 Fig. 10. BER performance of 64-QAM using 5G NR channel coding blocks for different rates.

Table 1 .

 1 Main parameters of the DP-MZM.

	Parameter	value
	Extinction ratio	30 dB
	Switching voltage 𝑉 𝜋 𝐷𝑃	4 volt
	Bias voltage 1	4 volt
	Bias voltage 2	4 volt
	Bias voltage 3	2 volt

  Eq. (3) and Eq. (4) into Eq. (2), it is obtained

	𝐸 𝐷𝐷-𝑀𝑍𝑀 (𝑡) = √ 𝑃 𝐿𝑎𝑠𝑒𝑟 𝑡 𝑎𝑡𝑡𝑛 2	[𝑒	𝑗	𝜋 𝑉 𝜋 𝐷𝐷 (𝑉 𝑚 cos(2𝜋𝑓 𝐿𝑂3 𝑡)+𝑉 𝑏𝑖𝑎𝑠 𝐷𝐷 )	+ 𝑒	𝑗	𝜋 𝑉 𝜋 𝐷𝐷 𝑉 𝑚 cos(2𝜋𝑓 𝐿𝑂3 𝑡+𝜃)

7

  𝐸 -1 = 𝐸 𝑎𝑡𝑡𝑛 . 𝐽 -1 (𝑚). (𝑒 𝑗𝜋𝛾 + 𝑒 -𝑗𝜃 ). 𝑒 𝑗2𝜋(𝑓 𝑐 -𝑓 𝐿𝑂3 )𝑡 . 𝑒 -𝑗 𝜋 2 𝐸 0 = 𝐸 𝑎𝑡𝑡𝑛 . 𝐽 0 (𝑚). (𝑒 𝑗𝜋𝛾 + 1). 𝑒 𝑗2𝜋𝑓 𝑐 𝑡 (8)𝐸 1 = 𝐸 𝑎𝑡𝑡𝑛 . 𝐽 1 (𝑚). (𝑒 𝑗𝜋𝛾 + 𝑒 𝑗𝜃 ). 𝑒 𝑗2𝜋(𝑓 𝑐 +𝑓 𝐿𝑂3 )𝑡 . 𝑒 𝑗 𝜋 2 𝐸 2 = 𝐸 𝑎𝑡𝑡𝑛 . 𝐽 2 (𝑚). (𝑒 𝑗𝜋𝛾 + 𝑒 𝑗2𝜃 ). 𝑒 𝑗2𝜋(𝑓 𝑐 +2𝑓 𝐿𝑂3 )𝑡 . 𝑒 𝑗𝜋 𝐸 3 = 𝐸 𝑎𝑡𝑡𝑛 . 𝐽 3 (𝑚). (𝑒 𝑗𝜋𝛾 + 𝑒 𝑗3𝜃 ). 𝑒 𝑗2𝜋(𝑓 𝑐 +3𝑓 𝐿𝑂3 )𝑡 . 𝑒 -𝑗 𝜋 2

Table 2 .

 2 Parameters used to generate OFDM MIMO signal.

	Parameter	Value
	Sampling rate	10 GS/s
	FFT/IFFT size	512
	Number of used data subcarriers	336
	Number of pilot subcarriers	32
	Cyclic Prefix	64

Table 3 .

 3 Main parameters of the optical architecture simulation.

	Parameter	Value

  shows the obtained BER as a function of Eb/No for QPSK with data rate of 46,64 Gb/s, 16-QAM with data rate of 93,28 Gb/s and 64-QAM with data rate of 140 Gb/s, without using 5G NR channel coding. As we can see, a BER of 10 -3 can be achieved by QPSK, and 16-QAM at Eb/No of 17 dB and 22 dB, respectively.

Declaration of interests

The authors declare that they have no known competng fnancial interests or personal relatonships ☒ that could have appeared to influence the work reported in this paper.

The authors declare the following fnancial interests/personal relatonships which may be considered ☐ as potental competng interests:


Rayan Ziane

Sylvain Kahane

Martine Vanhove

Bruno Guillaume

A morph based and a word based treebank for beja

ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.