Philippe Virouleau

François Broquedis

Thierry Gautier
email: thierry.gautier@inrialpes.fr

Fabrice Rastello

Using data dependencies to improve task-based scheduling strategies on NUMA architectures

Keywords: OpenMP, task dependencies, benchmark, runtime systems, NUMA, XKAAPI, scheduling, work-stealing

come

Introduction

While non-uniform memory access (NUMA) architectures stand today as one of the most popular design to build large-scale shared memory machines, exploiting them at their full potential remains challenging. On such architectures, the memory is split into several NUMA nodes and both bandwidth and latency depend on which processor accesses specific data : accessing memory allocated locally is most of the time faster than accessing data allocated to remotely-located NUMA nodes. Controlling data locality over the application lifetime is one of the key steps to achieving both good performance and scalability on these architectures.

Task-based parallel programming environments like OpenMP have become very popular when it comes to program shared memory machines with hundreds of cores. Indeed, they offer ways of expressing massive fine-grain parallelism with a relatively low overhead. Most of them also come with facilities to dynamically perform load balancing of tasks over the processors. Even if such characteristics fill the need of generating more and more parallelism out of parallel applications, standard parallel programming environments still not explicitly address the problem of data locality on NUMA systems.

The runtime system plays a central role in the execution of a task-based parallel application. For example, it is responsible for assigning ready tasks to the target platforms' processors. It is also in charge of performing load balancing when a processor idles. Both these decisions should take the architecture topology into account in order to avoid NUMA-related performance penalties on the overall application performance.

The recent addition of data dependencies to the OpenMP tasking model provides the runtime system with very precise information about which part of an application accesses which variables. Thanks to these dependencies, the runtime system knows which memory areas are read or written by which task. As shown in this paper, the task scheduler can rely on this information when assigning tasks to processors to implement NUMA-aware strategies.

This paper describes several of these strategies we implemented inside the XKAAPI [START_REF] Gautier | Kaapi: a thread scheduling runtime system for data flow computations on cluster of multi-processors[END_REF] runtime system. We identified three major steps in the task scheduler workflow that may have an impact on parallel applications on NUMA systems : the data distribution, the assignment of ready tasks to the processors and the way the task scheduler browses the architecture topology to perform load balancing. This paper describes and evaluates them, showing how they impact the application performance on a 192-core NUMA machine. We also compare them to state-of-the-art task scheduling strategies taken from related works and implemented within XKAAPI.

The layout of this paper falls into six sections as follows. In Section 2, we first give some background on NUMA architectures and the task programming model with data dependencies. We then describe in Section 3 the ideas, strategies and implementation details that we used to improve the runtime performances for these applications. Section 4 is devoted to the presentation performances evaluation. We eventually present some related works in Section 5 before concluding.

NUMA architectures design and exploitation

Hardware background

Most of nowadays parallel shared memory architectures are built according to a NUMA design where the memory is physically split into several banks attached to processors. Many vendors assemble these banks in a hierarchical way, thus building shared memory machines embedding several hundreds of cores. Exploiting such architectures at their full potential requires a fine control of the execution of a parallel application, as accessing local memory is most of the time faster than accessing memory stored in a memory bank attached to a remote processor.

The machine we experimented on is an SGI UV2000 platform made of 24 NUMA nodes. Each NUMA node holds an 8-core Intel Xeon E5-4640 CPU for a total of 192 cores. We refer to this machine as Intel192 in the paper. The memory topology is organized by pairs of NUMA nodes connected together through Intel QuickPath Interconnect. These pairs can communicate together through a proprietary fabric called NU-MALink6 with up to two hops.

Table 1 shows the distances advertised by the hwloc library [START_REF] Broquedis | hwloc: A generic framework for managing hardware affinities in HPC applications[END_REF] that represents the communication time for different distances normalized to the time of a local communication. Distances named local and peer form a pair of NUMA nodes (through Intel QPI), other nodes are either one hop away or two hops away (through NUMALink6).

Software background

To exploit large-scale shared memory architectures, the application programmer needs:

1. to express massive fine grain parallelism to get the most out of the numerous processing units of the platform ; 2. to control the execution of the application, especially the way computations and data are distributed over the platform, to prevent the NUMA design to have a negative impact on the overall application performance.

Task-based parallel programming environments provide ways of expressing fine grain parallelism that can be dynamically assigned to processors at runtime. OpenMP [START_REF]OpenMP application program interface version 4[END_REF], the de-facto standard for shared-memory parallel programming, supports task parallelism with dependencies since revision 4.0.

A glimpse at OpenMP tasking An OpenMP task can be seen as an independent unit of work an OpenMP thread can execute. Tasks can be created by an OpenMP thread and executed by any thread of the same parallel region. As managing tasks at runtime is way cheaper than creating and synchronizing threads, the application programmer can take the parallelization of its application further, as he can now consider portions of code that were too fine grain to be parallelized using only threads. The synchronization of OpenMP 3.0 tasks is performed thanks to the taskwait keyword that waits for the completion of all the tasks generated from the current OpenMP parallel region. On one hand, the application programmer is responsible for creating and synchronizing OpenMP tasks explicitly. On the other hand, the runtime system is in charge of correctly assigning tasks to threads during the application execution.

OpenMP 4.0 pushes the concept of task further introducing the depend keyword to specify the access mode of each shared variable a task will access during its execution. Access modes can be set to either in, out or inout whether the corresponding variable is respectively read as input, written as output or both read and written by the

[m][k],A[n][k])\ 21 depend(inout: A[m][n]) 22 dgemm(NB, 23 &A[m][k],&A[n][k],&A[m][n]); 24 } 25 }
Fig. 1: Cholesky factorization with OpenMP-4.0 task dependencies considered task. This information is then processed by the underlying runtime system to decide whether a task is ready for execution or should first wait for the completion of other ones.

KASTORS Benchmark suite

Listing 1 shows the implementation of a Cholesky factorization implemented with OpenMP task dependencies. This factorization algorithm comes from the PLASMA library and is very similar to the one implemented in the KASTORS benchmark suite [START_REF] Virouleau | Evaluation of OpenMP Dependent Tasks with the KASTORS Benchmark Suite[END_REF]. Task dependencies support comes with several benefits. First, task dependencies involve decentralized, selective synchronization operations that should scale better than the broad-range taskwait-based approaches. In some situations, this way of programming unlocks more valid execution scenarios than explicitly synchronized tasks, which provides the runtime system with many more valid task schedules to choose from. For example, in the Cholesky factorization, many instances of the dtrsm, dsyrk and dgem BLAS computations can legally run concurrently when executing the version with task dependencies. Secondly, information about task dependencies also enables the runtime system to optimize further, such as improving tasks and data placement.

The way we execute task-based applications Most task-based programming environments rely on a work-stealing execution model, originally introduced in Cilk [START_REF] Frigo | The implementation of the cilk-5 multithreaded language[END_REF]. Work-stealing is indeed often considered when it comes to dynamically balance the workload among processing units. The work-stealing principle can be summarized as follows. An idle thread, called a thief, initiates a steal request to a randomly selected victim. On reply, the thief receives a copy of one ready task, leaving the original task marked as stolen. Coherency between a thief and its victim is ensured by a variant of Cilk's T.H.E protocol, also described in [START_REF] Frigo | The implementation of the cilk-5 multithreaded language[END_REF].

The runtime system we develop, called XKAAPI, also implements the work-stealing execution model to execute OpenMP task-based applications. The runtime creates a system thread, called a kproc, for each processing unit to be used. On a NUMA multicore machine, a processing unit is a core. A kproc creates tasks and pushes them on its own work queue, which is implemented as a stack. The enqueue operation is very fast : it takes around ten cycles on modern x86/64 processors [START_REF] Broquedis | Libkomp, an efficient openmp runtime system for both fork-join and data flow paradigms[END_REF]. As in Cilk, a running XKAAPI task can create children tasks. Depending on the number of tasks per thread, XKAAPI implements two strategies to find a ready task. If the number of tasks is lower than a threshold, XKAAPI follows the Cilk's work first principle. In that case, the thief iterates through the victim's stack queue from the least recently pushed task to the most recently one and it computes true data-flow dependencies for each task until a ready task is found. If the number of tasks is greater than the threshold, the data flow graph is built and the thief picks a task from the victim's list of ready tasks [START_REF] Bleuse | Data Flow Program in XKaapi: A New Affinity Based Algorithm for Heterogeneous Architectures[END_REF]. By the nature of our benchmarks, this latter strategy is de facto selected and we developed new NUMA aware strategies.

3 Using OpenMP tasks dependencies to improve tasks and data placement on NUMA machines

In this section, we describe how the runtime system can have a positive impact on the application execution using the information provided by data dependencies. The following sections describe the way we adapted the behavior of the runtime system to control data placement during the initialization phase, when data will be allocated and accessed for the first time, and how we modified the way tasks that perform the actual computations are dynamically assigned to processors while maximizing data locality.

Inside the XKAAPI task-based runtime system

This section describes some of the key internal structures and mechanisms of the XKAAPI runtime system.

The way XKAAPI models the architecture. XKAAPI sees the architecture topology as a hierarchy of places. A place is a list of tasks associated with a subset of the machine processing units. XKAAPI's places are very similar to the notion of shepherd introduced in [START_REF] Olivier | Openmp task scheduling strategies for multicore NUMA systems[END_REF], or ForestGOMP's runqueues [START_REF] Broquedis | ForestGOMP: an efficient OpenMP environment for NUMA architectures[END_REF]. XKAAPI most of the time only considers two levels of places : node-level places, which are bound to the set of processors contained in a NUMA node, and processor-level places, which are bound to a single processor of the platform. This way, at the processor level one place is associated to each of the physical cores, and at the NUMA node level one place is associated to each of the NUMA nodes.

The way XKAAPI enables ready tasks and steals them. The scheduling framework in XKAAPI [START_REF] Bleuse | Data Flow Program in XKaapi: A New Affinity Based Algorithm for Heterogeneous Architectures[END_REF] relies on virtual functions for selecting a victim, selecting a place to push a ready task and pushing a set of initial ready tasks. When a processor becomes idle the runtime system has to select a victim and calls a function, called WSselect for work-stealing select, to browse the topology to find a place from which stealing a task from the place task queue.

The completion of a task may unlock the execution of some of its children in the dependency graph. This means marking them as ready for execution and pushing at least one of them to a place. Once again, there are many ways of selecting the place where to push ready tasks, implemented in strategies we refer to as WSpush, for work-stealing push.

Before parallel computation begins, the runtime system can distribute (push) the set of ready tasks to multiple places, according to the strategy defined by WSpush_init.

These three functions are the main entry points to specify a scheduling algorithm in XKAAPI. Sections 3.3 to 3.5 describe strategies for these three points, all them were designed to explore the possibilities of the target NUMA architectures, to be able to evaluate which one are worth taking into account.

Controlling data distribution on a NUMA system

Controlling the way data are allocated on a NUMA system requires a good understanding of the underlying memory architecture. Application programmers can achieve this using dedicated tools or libraries, like libNUMA's numactl [START_REF] Wittmann | Optimizing ccnuma locality for task-parallel execution under openmp and TBB on multicore-based systems[END_REF], which can be used to set a default memory allocation policy for the whole application. For example, the --interleave=all memory policy spreads out all the memory pages of dynamically allocated variables, over all the NUMA nodes of the machine. This policy is widely used on NUMA systems in conjunction with dynamic parallelism, like taskbased programs, as it distributes the memory traffic over all the memory controllers, making processors "all equally bad" when it comes to memory access. To better control data placement, parallel application programmers are used to relying on the first-touch allocation policy, which is the default behavior for memory allocation on most Linux systems. This allows allocating memory pages when they are accessed for the first time.

To better control data distribution on NUMA systems, we propose two different approaches :

either the application programmer explicitly allocates data on specific NUMA nodes of the machine through a dedicated API we provide [START_REF] Durand | Proceedings of the 9th International Conference on OpenMP in the Era of Low Power Devices and Accelerators[END_REF] (omp_locality_domain_allocate_XXX) where XXX may be a bloc cyclic data distribution for one or two-dimensional arrays over MAMI [START_REF] Broquedis | ForestGOMP: an efficient OpenMP environment for NUMA architectures[END_REF]; or the application programmer only marks some regions of code that initialize data to give the runtime system the opportunity to map the corresponding tasks to make the first-touch allocation policy indirectly apply the data distribution we target. Indeed, Olivier et al. [START_REF] Olivier | Characterizing and mitigating work time inflation in task parallel programs[END_REF] have shown that specifying affinity for initialization tasks can lead to huge improvement over locality oblivious techniques. To avoid remote memory accesses, the threads must access the data during the computation phase the exact same way it was accessed during the initialization phase, which is very difficult to guarantee with dynamic task-based parallelism. We extend the OpenMP runtime in two ways. First, by adding functions to provide a dedicated API: omp_set_affinity to make the runtime map the next task to a specific NUMA node. Secondly, by extending scheduling heuristics to take into account task's dependencies to better map ready tasks.

During the application's execution, the runtime relies on system's get_mempolicy to determine on which physical node data are allocated. This information is then used to guide the way we perform task creation and load balancing.

Distribution of initial ready tasks : WSpush_init strategies

We refer to initial tasks when considering the sources of a task dependency graph, usually declared at the beginning of an OpenMP parallel region. These tasks are basically the first ones to be marked as ready and to be distributed over the platforms' places. We have implemented two initial tasks distribution strategies : cyclicnuma which distributes the tasks in a round-robin fashion over the NUMA nodes, and randnuma which randomly distributes the tasks over the NUMA nodes. Note that unlike numactl, the strategies we implemented consider the whole data appearing in the OpenMP task depend clause instead of working at the page level. In other words, while the two memory pages holding an 8K-wide array would be distributed on different nodes by numactl --interleave=all, they are always assigned to the same NUMA node when using one of our data distribution strategies.

Distribution of ready tasks : WSpush strategies

This section describes four different ways of pushing ready tasks to a NUMA system places. Two of them are data-oblivious while the other two rely on the dependencies expressed using the depend keyword on OpenMP tasks. The pLoc strategy makes a processor push ready tasks to its own place, while the pLocNum strategy makes a processor push ready tasks to the place of its NUMA node (local NUMA node). The pNumaW strategy pushes tasks on the node-level place corresponding to the NUMA node where most of their output data are allocated to (W stands for Write). The last WSpush strategy, called pNumaWLoc, behaves almost the same than pNumaW except that if the data are allocated to the NUMA node of the processor pushing the task, we directly push the task to this processor's place instead of pushing it to the node-level place (Loc stands for Local).

It's important to note that pLoc and pLocNum does not take initial data placement into account, while pNumaW and pNumaWLoc are both aware of where a task's data are physically allocated and which of them are written, thanks to the OpenMP depend keyword.

Dynamic load balancing using work-stealing : WSselect strategies

Another important step when implementing work-stealing is the selection of the victim processor we want to steal from. This section describes the selection strategies we implemented, that take the architecture memory hierarchy into account. The first two strategies, sRand and sRandNuma are similar to those studied in [START_REF] Olivier | Openmp task scheduling strategies for multicore NUMA systems[END_REF] and distinguish two levels of hierarchy : the processor level and the NUMA node level. sRand selects a random processor's place while sRandNuma selects a random NUMA node's place. We additionally implemented several strategies mixing both levels of hierarchy, described below.

-sProcNuma : First, we browse the processor's place. Upon failure, we browse the topology in the following order : we first browse one of the neighbor processors ; when all the neighbors have been visited, we browse the local NUMA place ; we continue by browsing all the processors' places from a random remote node and we eventually consider the place of its NUMA node. -sNumaProc : This strategy is similar, except we always look at the NUMA place before looking at the processors' place. -sProc : In this strategy the stealer will visit only the processors' places and its own NUMA place. -sNuma : In this strategy the stealer will visit only NUMA places and its neighbors.

Like proposed in [START_REF] Olivier | Characterizing and mitigating work time inflation in task parallel programs[END_REF], all these strategies come in two versions : a strict version in which we prevent processors from stealing from other NUMA nodes to improve data locality and a loose version where these restrictions do not apply.

Evaluation

We ran all our experiments on the Intel192 machine described in section 2.1. We evaluated our strategies using the KASTORS [START_REF] Virouleau | Evaluation of OpenMP Dependent Tasks with the KASTORS Benchmark Suite[END_REF] 1 benchmark suite. More specifically, we used the dependent tasks version of the blocked QR factorization (dgeqrf_taskdep), and of the blocked Cholesky factorization (dpotrf_taskdep). These applications rely on kernels from BLAS and LAPACK libraries provided by OpenBLAS 2.15. We used the OpenMP GCC compliant runtime libKOMP [START_REF] Broquedis | Libkomp, an efficient openmp runtime system for both fork-join and data flow paradigms[END_REF] based on XKAAPI runtime system. We tagged the version we used on XKAAPI's git repository2 in the branch public/europar2016. For all of the above applications, we used the GCC 5.2.0 compiler. We also made our execution log files public 3 , as well as all the scripts we used, so that anyone may reproduce our data analysis, and look at the other results we did not put in the figures.

Impact of the data distribution

We first evaluated the impact the initial data distribution has on the application performance. We did an evaluation for multiple matrix sizes and block sizes, as well as multiple combinations of WSpush and WSselect strategies. Figure 2 reports the results we obtained for the Cholesky application, on 32K-wide matrices divided into blocks of 512×512 elements. We observed similar behavior running Cholesky on different matrix sizes (16K to 64K) and block sizes (256 to 1024). The lower double dashed horizontal line is the GCC performance baseline using sequential initialization. The middle dashed line is the same experiment using numactl. The upper solid line is the GCC baseline using parallel initialization.

Using numactl provides an important performance gain compared to the sequential initialization. However using a parallel initialization, either controlled (cyclicnuma, randnuma) or not (GCC init-para), is necessary to significantly improve the performances, regardless of the strategies used.

The cyclicnuma distribution is the one that works best regardless of the strategies, and we will use it as the default strategy for the next experiments. Cholesky's performances using 32K matrices depending on the data distribution Fig. 2: Evaluating data distributions for multiple strategies (WSselect + WSpush)

GCC init-seq + Numactl GCC init-seq + Numactl GCC init-seq + Numactl GCC init-seq + Numactl GCC init-seq + Numactl GCC init-seq + Numactl GCC init-seq + Numactl GCC init-seq + Numactl GCC init-para GCC init-para GCC init-para GCC init-para GCC init-para GCC init-para GCC init-para GCC init-para GCC init-seq GCC init-seq GCC init-seq GCC init-seq GCC init-seq GCC init-seq GCC init-seq GCC init-seq

Impact of the stealing restriction

Given a data distribution, previous works [START_REF] Olivier | Characterizing and mitigating work time inflation in task parallel programs[END_REF] have shown that restricting the task execution to the node where the data are written leads to better data locality which may improve the application performance. However, this is heavily dependent on the algorithm the application implements. For instance, in the case of a Cholesky factorization, many tasks write to the diagonal tiles of the matrix comparatively to other tiles of the matrix. Therefore applying a steal restriction on these tasks will potentially lead to an important number of inactive processors. We evaluated both strict and loose versions of our work-stealing strategies and found out that preventing processors from stealing from other NUMA nodes can lead to a loss of performance by around 25% to more than 75% with respect to the same setup without the strict restriction. For the sake of brevity we did not include a figure for this, but the results of these experiments are included in the logs publicly published.

Overview of the strategies performances

We took a given data distribution, cyclicnuma, and compared the different strategies, without any steal restriction. The performance obtained running the Cholesky application executed by the libGOMP 4 runtime system (without modification) is considered as a baseline for these experiments. Once again, even if the performances we obtained are obviously not the same, the behavior of the different strategies comparatively to each others are similar for the different applications we ran. Figure 3 shows the results of the experiments for the Cholesky application on 32K-wide matrices divided into blocks of 512×512 elements (best configuration for this matrix size). The dashed horizontal line is the GCC performance baseline using parallel initialization.

It is first interesting to note that even very basic WSselect and WSpush strategies, like sRand+pLoc, obtained decent performances thanks to the data distribution. Also, given a selection strategy (e.g. sRandNuma), placing the task on the NUMA GCC GCC GCC GCC GCC GCC GCC Cholesky and QR performances for multiple sizes Fig. 4: Evaluating specific strategies on multiple sizes, using cyclicnuma WSpush_init node where the written data are allocated (pNumaW) behaves better than simply pushing the data to its NUMA node (pLocNum). However, assuming the tasks are being pushed using the pNumaWLoc strategy, focusing the place selection on only one level of the hierarchy (sProc or sNuma) fails to reach the same level of performance we obtained with naive strategies. On the contrary, taking into account both levels of the hierarchy (sProcNuma, sNumaProc) achieve similar performance that outperforms other strategies.

Strategies performance scaling

We eventually selected the three strategy combinations that outperformed the GCC baseline to evaluate their scalability depending on the size of the input matrix. These strategy combinations are :

-sRand + pLoc, which is a basic strategy that does not take the architecture topology into account ;

-sNumaProc + pNumaWLoc, which was the best strategy in our previous evaluation and is also equivalent to using sProcNuma ; -sRandNuma + pNumaW that performs random selection of node-level places.

Figure 4 reports their performances using a cyclicnuma distribution without steal restriction. The figure shows the performances using the best block size for each matrix size (which is, for our setup, 256 for a matrix size of 16384, and 512 for the others). As expected, combinations of strategies taking both the architecture hierarchy and data locality into account (sNumaProc + pNumaWLoc) achieve the best performances. The only exceptions are for small matrix sizes (16384), where there may just be not enough work to be able to take advantage of this strategy. We must note that simply distributing the data over the nodes enables the basic sRand + pLoc combination to achieve satisfying performances.

Finding the appropriate combination is highly application-dependent, therefore it is hard to give a solution for every type of application. However some general guidelines can be followed:

-The most critical part is the initial data distribution. Figure 2 showed it is absolutely necessary to use one. -Hierarchical strategies have a cost, so the problem size has to generate enough work and data transfer to see a benefit.

Related work

Numerous works focus on data locality and/or topology-aware task scheduling strategies for NUMA architectures. Clet-Ortega et al. [START_REF] Clet-Ortega | Evaluation of openmp task scheduling algorithms for large NUMA architectures[END_REF] studied different ways of decorating the architecture topology with task lists and how it impacts the performance of task-based applications on NUMA systems, promoting private per-threads lists of tasks browsed in a hierarchical way by work-stealing strategies. We somehow extended this work considering also node-level task lists. We showed considering these lists for pushing ready tasks and selecting work-stealing victims can help improving performance on NUMA systems. Olivier et al. [START_REF] Olivier | Openmp task scheduling strategies for multicore NUMA systems[END_REF] evaluated hierarchical task scheduling with respect to traditional centralized or distributed task schedulers. Creating a thread list, called shepherd, per NUMA node allowed their hierarchical scheduler to outperforms other approaches on several task-based applications. Tahan et al. [START_REF] Tahan | Towards efficient openmp strategies for non-uniform architectures[END_REF] also studied the behavior of task-based OpenMP applications on NUMA systems, extending the NANOS runtime system with two NUMA-aware task schedulers called DFWSPT and DFWS-RPT, taking into account the notion of task priority when pushing tasks to core-level queues. They also try to minimize the number of memory hops when performing load balancing. Drebes et al. [START_REF] Drebes | Topology-aware and dependence-aware scheduling and memory allocation for task-parallel languages[END_REF] proposed similar ideas in another dataflow programming model named OpenStream. This model has a focus on data streaming and has a lot of flexibility on their placement, but does not provide flexibility to the user. While the same kind of studies have been conducted in other contexts [START_REF] Terboven | Task-parallel programming on NUMA architectures[END_REF][START_REF] Wittmann | Optimizing ccnuma locality for task-parallel execution under openmp and TBB on multicore-based systems[END_REF][START_REF] Weng | Implementing openmp using dataflow execution model for data locality and efficient parallel execution[END_REF], none of them takes advantage of the OpenMP depend clause, which precisely indicates which data are read and written by a given task. As advertised by the results obtained by our sNumaProc+pNumaWLoc combined strategy, this information is worth taking into account when choosing a place to push ready tasks to.

Task-based programming environments like OpenMP have become a standard way to program large-scale NUMA systems. Indeed, they give the programmer ways of expressing massive fine-grain parallelism that can be dynamically mapped to the architecture topology at runtime. OpenMP recently evolved to deal with tasks dependencies describing the data a task reads as input and writes as output.

This paper presented several runtime-level strategies to efficiently assign tasks to processors on any NUMA architecture. We presented strategies assigning ready tasks to lists of tasks, called places, attached to processors and NUMA nodes. These strategies define the way a task-based runtime system pushes ready tasks to their initial place and the way idle processors browse the architecture topology to select a place to steal from. We considered several initial data distributions and evaluated different combinations of "push" and "select" strategies on a 192 core NUMA system, on linear algebra applications. We achieved the best performance with strategies taking into account both the architecture topology and the initial data placement obtained through OpenMP tasks dependencies.

A short-term future work will be to extend an OpenMP compiler to be able to identify the initialization tasks in a more OpenMP-friendly manner, like extending the task construct with a init clause. We also intend to experiment with more OpenMP 4 applications. In a longer term, we intend to move our focus to compile-time techniques able to infer and to attach valuable information on tasks, like an estimation of a task operational intensity, that could guide some of the runtime system's decisions regarding task scheduling and load balancing. We strongly believe a tight cooperation between the compiler and the runtime system is a key step to enhance the performance and scalability of task-based programs on large-scale platforms.

  for (size_t k=0; k < NB; ++k) { #pragma omp task shared(A) \ depend(inout: A[k][k]) dpotrf(NB,&A[k][k]); for (int m=k; m < NB; ++m) #pragma omp task shared(A)\ depend(in: A[k][k]) \ depend(inout: A[m][k]) dtrsm(NB,&A[k][k],&A[m][k]); 12 for (int m=k; m < NB; ++m) { 13 #pragma omp task shared(A)\ 14 depend(in: A[m][k]) \ 15 depend(inout: A[m][m]) 16 dsyrk(NB,&A[m][k], &A[m][m]); 17 18 for (int n=k; n < m; ++n) 19 #pragma omp task shared(A)\ 20 depend(in: A

Fig. 3 :

 3 Fig.3: Evaluating all strategies (WSselect + WSpush), using cyclicnuma WSpush_init

Table 1 :

 1 NUMA distances from node 0 advertised by the hwloc library on Intel192.

	NUMA nodes location	local	peer	one hop away	two hops away
	hwloc distances	1.0	5.0	6.5	7.9

git available at https://scm.gforge.inria.fr/anonscm/git/kastors/kastors.git, tag "tag-europar16"

https://scm.gforge.inria.fr/anonscm/git/kaapi/xkaapi.git

https://github.com/viroulep/europar-2016-public

Acknowledgments

This work is integrated and supported by the ELCI project, a French FSN ("Fond pour la Société Numérique") project that associates academic and industrial partners to design and provide software environment for very high performance computing.


Kristelle Bougot-Robin

Anne Talneau

Henri Benisty
email: henri.benisty@institutoptique.fr

Low-index nanopatterned barrier for hybrid oxide-free III-V silicon conductive bonding

Oxide-free bonding of a III-V active stack emitting at 1300-1600 nm to a silicon-on-insulator wafer offers the capability to electrically inject lasers from the silicon side. However, a typical 500-nm-thick silicon layer notably attracts the fundamental guided mode of the silicon + III-V stack, a detrimental feature compared to established III-V Separate-Confinement Heterostructure (SCH) stacks. We experimentally probe with photoluminescence as an internal light source the guiding behavior for oxide-free bonding to a nanopatterned silicon wafer that acts as a low-index barrier. We use a sub-wavelength square array of small holes as an effective "low-index silicon" medium. It is weakly modulated along one dimension (superperiodic array) to outcouple the resulting guided modes to free space, where we use an angle-resolved spectroscopy study. Analysis of experimental branches confirms the capability to operate with a fundamental mode well localized in the III-V heterostructures.

Introduction

The interest of hybrid silicon photonics, with III-V active material bonded onto a platform such as Silicon-On-Insulator (SOI), is undisputed [START_REF] Roelkens | III-V/silicon photonics for on-chip and inter-chip optical interconnects[END_REF][START_REF] Doylend | The evolution of silicon photonics as an enabling technology for optical interconnection[END_REF][START_REF] Keyvaninia | Heterogeneously integrated III-V/Si single mode lasers based on a MMI-ring configuration and triplet-ring reflectors[END_REF]. One first point was to get guidance from well mastered low-loss silicon structures, and to get active emission from the III-V epitaxial heterostructure stack (denoted EHS below), either in combination or for separate action, with a transfer in between. For these optical functions, bonding can make use of oxide, or polymer (e.g., benzocyclobutene BCB). Attempts have been pursued for bonding layer down to a few nm or of molecular nature [START_REF] Bakish | Self-assembled monolayer assisted bonding of Si and InP[END_REF]. GaAs-Si bonding without oxide is also sought actively [START_REF] Tanabe | III-V/Si hybrid photonic devices by direct fusion bonding[END_REF][START_REF] Tanabe | 1.3 μm InAs/GaAs quantum dot lasers on Si rib structures with current injection across direct-bonded GaAs/Si heterointerfaces[END_REF] for the useful 1300 nm emission. Oxide-free bonding [START_REF] Talneau | Atomic-plane-thick reconstruction across the interface during heteroepitaxial bonding of InP-clad quantum wells on Silicon[END_REF], is the only known technique to date that ensures not only a mechanical and thermal bond, but also an electrical bond, i.e. a low impedance non rectifying contact between silicon and EHS, as shown recently [START_REF] Pantzas | Electrical transport across the heterointerface of InP membranes bonded oxide-free on Si[END_REF] for InP-Si bonding (Refs [START_REF] Tanabe | III-V/Si hybrid photonic devices by direct fusion bonding[END_REF][START_REF] Tanabe | 1.3 μm InAs/GaAs quantum dot lasers on Si rib structures with current injection across direct-bonded GaAs/Si heterointerfaces[END_REF] show a non-rectifying Si-GaAs contact). Echoing the numerous novel uses of subwavelength structures in silicon photonics [START_REF] Schmid | Refractive index engineering with subwavelength gratings in silicon microphotonic waveguides[END_REF], oxide-free bonding of nanostructured silicon to III-V EHS was also demonstrated successfully [START_REF] Itawi | Void-free direct bonding of InP to Si: Advantages of low H-content and ozone activation[END_REF], with positive photonics tests [START_REF] Pang | Oxide-free InP-on-Silicon-on-Insulator nanopatterned waveguides: propagation losses assessment through end-fire and internal probe measurements[END_REF]. Besides a possible ohmic use, nanopatterns offer the capability of lateral guidance combined with additionally lateral thermal sinking in the upper silicon slab [START_REF] Pang | Nanostructured silicon geometries for directly bonded hybrid III-V-silicon active devices[END_REF].

In this work, we argue that when an electrical contact across the bonded Si/EHS interface is desired, the natural vertical modal engineering practiced with the oxide (SiO 2 ) or polymer barrier is absent. Then, the stack's fundamental mode (FM) may widely overlap the Si slab, hence a penalty for a laser structure, as the confinement factor Γ FM of the active quantum wells (QWs) is reduced. To circumvent this penalty, we show that a nanopatterned silicon consisting of an array of holes with a modest air fraction f acts as a low-index barrier between Si and EHS, and thus greatly assists laser operation in the established regime of III-V's Separate Confinement Heterostructures (SCH). To evidence this, we analyze the dispersion of the Si-EHS stack: we excite guided modes by the localized QWs photoluminescence and outcouple them thanks to a superperiodic modulation of the basic sub-wavelength nanopattern. We first detail the sample, the experiment, and illustrate the TE polarized mode distribution undergoing confinement by the low-index barrier. We finally present the main experimental results.

Sample and experimental set-up

We show a schematic cross-section of our stack in Fig. 1(a) [START_REF] Talneau | Atomic-plane-thick reconstruction across the interface during heteroepitaxial bonding of InP-clad quantum wells on Silicon[END_REF]. The silicon slab of the SOI is 500-nm-thick, the EHS has a 380-nm-thick SCH with QWs emitting at 1505 nm, sandwiched between 100 nm InP on the Si side, and 820 nm InP on the top (air) side to avoid that the FM visits the top surface. A representative cross-section of an oxide-free bonded SOI-EHS with nanopatterned silicon (here about 200nm) and QW revelation is shown in Fig. 1(b). In our photonic study, the etch-depth is 270 nm, leaving 230 nm un-etched Si.

In Fig. 1(c), we sketch an injected laser diode exploiting silicon to spread the current below the QWs, thus reducing Joule heating near the active junction. High doping is then needed in silicon. Without nanopattern, this doping would cause residual absorption in Si as the FM would overlap this region (usually a TE mode, fitted to exploit e-hh transitions in QWs). It is therefore welcome to repel the FM profile up into the EHS, a feature that, in addition, increases Γ FM for this mode. The nanopattern lowers the Si index ~3.48 to an homogenized in-plane effective index n hom,xy ~2.8-3.3, well suited for this purpose. We will denote n eff,z the effective indices of guided modes (pertaining to index profile in the zdirection). The chosen nanopattern is a square array of holes. For the sole use as a low-index layer, any sub-wavelength period a < λ/2n eff,z ~215 nm suffices, but no easy diagnostic of its photonic role at device scale (say > 50 µm) can be made [START_REF] Pang | Oxide-free InP-on-Silicon-on-Insulator nanopatterned waveguides: propagation losses assessment through end-fire and internal probe measurements[END_REF]. We therefore add a weak onedimensional periodicity along x, by enlarging every third hole. As shown in Fig. 1(d), there are three hole diameters used in three pair combinations, lithography and etching technologies being optimized to deliver vertical walls. The result is a set of three designs, labelled f 1 (82,66), f 2 (82,54) and f 3 (66,54), (d i ,D i ) being the two diameters in nm, with decreasing air fractions f i (thus, increasing n eff,xy ). These modulated patterns are essentially equivalent to a superposition of (i) arrays of uniform diameters 72 nm, 65 nm, 58 nm for f 1 , f 2 , f 3 respectively (based on the same average hole area, e.g. 72 2 ≈[82 2 + 66 2 + 66 2 ]/3) with (ii) a 3a-period modulation. The diffraction efficiency of this modulation is given by the relevant Fourier component A G of the xy dielectric map at wavevector G = G x = 2π/3a along x (actually not needed here).

Typically, we operate to extract modes of indices 3.0-3.2 at λ = 1500 nm around normal incidence, hence 3a ~480 nm. The air-filling factors are then f 1 = 0.158, f 2 = 0.128, f 3 = 0.104. Applied to a uniform medium of index 3.48 (Si) with infinite air cylinders in the microscopic limit (d << λ), they have homogenized indices n hom,xy of about 3.04, 3.13 and 3.19, calculated from a standard two-dimensional plane-wave expansion. We note that ellipsometry could seem a suitable technique for establishing the index sequence, but for transparent substrates, thick top layers, small patterns (~50 µm) and small InP bonded pieces (large sizes are not routine yet), it seems too awkward. In Fig. 1(e,f) the setup used for measurement is explained. It resembles that of [START_REF] Benisty | Investigation of Extracting Photonic Crystal Lattices for Guided Modes of GaAs-Based Heterostructures[END_REF] used for guided mode extraction studies. A red laser excites a photoluminescent spot inducing radially expanding guided waves. Those waves launched along x meet the effective medium (70 × 30 µm patches) in which the nanopattern modulation causes diffraction and outcoupling. Behind the objective and the beam splitter, a slit is translated to select the beam waves with a given in-plane direction k || ≡ k x according to the Ewald construction of Fig. 1(e), k || = k g -G, k g = n eff,z k o being a guided mode wavevector, and denoting k o = 2π/λ. A second splitter provides imaging with an infrared camera as well as a fibered path to a spectrometer equipped with an InGaAs cooled array.

Dispersion analysis

We first calculated the several TE and TM guided modes of the SOI + EHS stack using tabulated Si, SiO 2 and InP dispersion laws n(λ) without absorption. In Fig. 2(a), we vary n hom,xy ("effective medium", Fig. 2(b)) and plot all branches at a fixed wavelength λ = 1505 nm. We replaced the whole SCH by a single 380 nm-thick purely dispersive layer with typical Sellmeier parameters [START_REF] Mondry | Refractive Indexes of (Al, Ga, I n)As Epilayers on InP for Optoelectronic Applications[END_REF], Fig. 2(c), chosen to fit experimental data for all modes of the stack.

The most salient feature around our range, n hom,xy ~3.0, is the TE1/TE2 anticrossing. The presence of such an anticrossing is typical of a mode residing or not on one side of the low index barrier upon interaction with the mode on the other side. As for the fundamental mode, there is a broad smeared anticrossing immediately at n hom,xy ~3.4, which is indicative of the Si attraction for the mode when it is pure. We hope to mitigate this feature thanks to the low index effective medium (it would be less smeared with a first InP layer of >100 nm). Nevertheless, assessing the small shifts of n eff,z for TE0 with our set-up is not very reliable. Instead the TE1/TE2 anticrossing offers a remarkable signature of the role of n hom,xy .

To understand more explicitly the role of the low-index barrier, we plot on Fig. 2(d-g) the modal profiles of the three first TE modes for illustrative situations close to the four actual cases: no holes, f 3 , f 2 , f 1 . In this order the index n hom,xy decreases. We push the lower bound to 2.93, instead of the actual one (~3.00) to better visualize the fate of higher modes.

We first see that for pure silicon, the FM mostly lies within the silicon layer. This pertains also to the choice of a rather thick silicon slab (500 nm). It is a good typical value, and could be useful, either to carry intense current, or to sink substantial heat sideways off a structure's active junction. But we also see that as soon as air is introduced, for the f 3 sample and n hom,xy ~3.23, the FM is clearly centered in the SCH, and its overlap with the nanopattern is quite low. As more air is introduced, going from f 3 to f 2 and to f 1 , the overlap of the SCH guided mode with the pattern still diminishes, and all parts of the mode profile are pushed toward the SCH.

The origin of the anticrossing of modes TE1 and TE2 can also be inferred from the profiles: for f 3 , TE1 occupies the silicon and TE2 lies rather in the III-V stack (the EHS), whereas for f 1 , we observe almost the contrary: TE2 now manages to be in the silicon region, while TE1 has a substantial lobe in the EHS. These considerations will be useful to explain experimental data that will be presented in the following section.

Dispersion characterization and discussion

The experimental angular resolved spectra shown below are obtained with a ~350 µm slit scanned on ~8 mm across the beam behind the objective whose limits are k || /k o ≈ ± 0.40 (NA). These angular spectra are represented as color maps in Fig. 3(a-c) for a = 160 nm and Fig. 3(d-f) for a = 170 nm. There is a rich content, with lines of quite modulated intensity and some quite visible gaps: in Brillouin zone center k || /k o = 0 for forward/backward coupling of the same mode, or at other positions if distinct modes are interacting. This is as expected from a broad periodic multimode guide [START_REF] Olivier | Mini stopbands of a one-dimensional system: the channel waveguide in a two-dimensional photonic crystal[END_REF][START_REF] Benisty | Crossing of manifolds leads to flat dispersion: Blazed Littrow waveguides[END_REF][START_REF] Khayam | General recipe for flatbands in photonic crystal waveguides[END_REF][START_REF] Rosenblatt | Resonant Grating Waveguide Structures[END_REF]. However, our concern is simply the modal sequence of underlying unperturbed modes (smooth branches without their gaps) to evidence the role of n hom,xy according to the above analysis. We are of course aware that gaps and signal strength reveal the degree of interaction of the mode with the lattice, further mitigated by various interference effects [START_REF] David | Photonic crystal light-emitting sources[END_REF]. Hence, to fit our modeling, we first adjusted thicknesses and angle calibration from all measured samples (dots on Fig. 3(a-f), TE in dark blue, TM in green), which includes some other periods not shown here, taking notably high order modes (n eff,z <2) into account. This adjustment consolidated the above mentioned thicknesses and SCH dispersion function of Fig. 2 (c). We adjusted in this process the value of n hom,xy to reproduce in particular the first three modes of the sequence, which, with our resolution and diffraction geometry, amount to TE0, TE1, TE2 (marked with bluish half-transparent ribbons in Fig. 3).

The indices that we found correspond to the grey bands of Fig. 2(a), whose width accounts for the slight data scatter (Δn ~0.03). Relatively accurate values of effective indices n eff,z were deduced from the spectra. Since they were found to fit with the expectations (see the dotted fits on all spectra maps), it also means that the absolute values for all branches in Fig. 2 It reflects the weak profile amplitude at the grating locus, showing that this mode is rather centered on the SCH, as noted in the above analysis. But the diffracted intensity cannot be naively used as a probe of that interaction because it involves also the quite variable reflection of bottom-diffracted waves onto the SOI stack [START_REF] Itawi | Void-free direct bonding of InP to Si: Advantages of low H-content and ozone activation[END_REF].

Conclusion

We have evidenced the positive role of a nanopatterned silicon acting as an effective lowindex layer for the good operation of active optoelectronics of a III-V stack in the new silicon photonics configuration permitted by oxide-free void-free bonding: Electrical contact seems effective in this geometry [START_REF] Pantzas | Electrical transport across the heterointerface of InP membranes bonded oxide-free on Si[END_REF], but the possible penalty of excessively attracting the useful (lasing) mode away from the active SCH into the silicon has to be mitigated. We have evidenced that a nanopatterned silicon layer of typical thickness 250-300 nm, air fraction f = 0.1-0.2, with a square lattice period 150-180 nm could work well in this respect. The hole diameters involved, about 55-85 nm, are well suited to state-of-the-art technology, making this approach helpful for forthcoming silicon photonic devices.

Fig. 1 .

 1 Fig. 1. (a) Effective medium made of nanostructured silicon in a hybrid silicon-III-V heterostructure elaborated by oxide-free bonding; (b) Micrograph of the slab section with quantum-wells revealed; (c) scheme of electrical injection and fundamental mode profile allowed by the lowered index of nanostructured silicon; (d) Modulated square lattices with every third hole slightly bigger with the three pairs of diameters used; (e) Ewald diagram for diffraction in air; (f) set-up to get an angular spectrum by translating a slit in the imaging path.

Fig. 2 .

 2 Fig. 2. (a) Expected effective indices of TE and TM modes for the simplified heterostructure as a function of an isotropic effective-material low index slab of 270 nm replacing the silicon top part. Note the anticrossing behavior of the TE1 and TE2 modes around an effective-material index n~3; (b) Structure stack; (c) Average index of a typical SCH layer; (d) Mode profiles of the first three TE modes in the non-etched bonded stack at λ = 1505 nm; (e) Same with a nanostructured layer of ~isotropic index n hom,xy = n Si -0.25~3.23; (f) Same with n hom,xy = n Si -0.40 ~3.08; (g) Same with n hom,xy = n Si -0.55 ~2.93. The fundamental mode is repelled in the III-V stack for lower n hom,xy . Its overlap with the nanostructured layer decreases. The anticrossing behavior of modes TE1 and TE2 clearly appears as Si guidance is driven to lower effective indices by the diminishing n hom,xy .

  (a) along the vertical grey bands are the actual ones for this oxide-free bonded sample, helping to clarify other future investigations. Values of n hom,xy are somewhat more spread than predicted by the hole-diameter-based spanning [3.00 → 3.25] instead of [3.04 → 3.19]. We see that the three first modes' unmistakable evolution along the sequence f 1 →f 2 →f 3 [darker shades in grey bands, Fig. 2(a)], is well reproduced experimentally in the observed branch spacing. Note the weakness of the TE0 line through all data.

  Fig. 3. (a-c) Experimental angular spectra (abscissa: k || /k o ≡ numerical aperture) for the three filling factors (f 1 , f 2 , f 3 ) and for the period 480nm = 3 × 160 nm; (d-f) Same for the period 510 nm = 3 × 170 nm; calculated TE branches are the blue dots, TM ones are the green dots. Bluish semi-transparent rectangles at top and bottom of maps underline modes TE0, TE1, TE2.

Acknowledgments

We acknowledge the support of French ANR through the COHEDIO project.


Priscilla Brandao Silva

Jean-Mathieu Mencik

Roberto Arruda

José R F Arruda

On the use of the wave finite element method for passive vibration control of periodic structures

Keywords: wave finite element method, periodic structures, band gaps, passive vibration control, aircraft fuselage

In this work, a strategy for passive vibration control of periodic structures is proposed which involves adding a periodic array of simple resonant devices for creating band gaps. It is shown that such band gaps can be generated at low frequencies as opposed to the well known Bragg scattering effects when the wavelengths have to meet the length of the elementary cell of a periodic structure. For computational purposes, the wave finite element (WFE) method is investigated, which provides a straightforward and fast numerical means for identifying band gaps through the analysis of dispersion curves. Also, the WFE method constitutes an efficient and fast numerical means for analyzing the impact of band gaps in the attenuation of the frequency response functions of periodic structures. In order to highlight the relevance of the proposed approach, numerical experiments are carried out on a 1D academic rod and a 3D aircraft fuselage-like structure.

Introduction

The development of passive vibration control strategies of industrial periodic structures, which use band gaps, is a current research topic which aims at proposing efficient lightweight solutions for global reduction of the vibration levels. Actually, there exist two band gap mechanisms which are either of Bragg scattering type [START_REF] Sigalas | Elastic and acoustic wave band structure[END_REF][START_REF] Kushwaha | Acoustic band structure of periodic elastic composites[END_REF][START_REF] Kushwaha | Theory of acoustic band structure of periodic elastic composites[END_REF], or locally resonant type [START_REF] Liu | Locally resonant sonic materials[END_REF][START_REF] Goffaux | Evidence of Fano-like interference phenomena in locally resonant materials[END_REF]. Regarding Bragg scattering, band gaps represent zones of destructive interference between incident and reflected waves which occur when the wavelengths are of the same order of magnitude as the dimensions of the periodic cells. This means that, at low frequencies, such band gaps only occur in large-sized structures. In elastodynamics, this mechanism has been highlighted in numerous engineering *Corresponding author, Ph.D., E-mail: pbrandaos@fem.unicamp.br a Ph.D., E-mail: jean-mathieu.mencik@insa-cvl.fr b Ph.D., E-mail: arruda@fem.unicamp.br systems, which are usually named phononic crystals, such as rods, plates and cylindrical shells [START_REF] Bennett | Free wave propagation in periodically ring stiffened cylindrical shell[END_REF][START_REF] Lee | Analysis of wave propagation in a thin composite cylinder with periodic axial and ring stiffeners using periodic structure theory[END_REF][START_REF] Sorokin | Plane wave propagation and frequency band gaps in periodic plates and cylindrical shells with and without heavy fluid loading[END_REF], Goldstein et al. 2011). More recently, researchers in the field have moved their attention to locally resonant devices which are motivated by the possibility of designing band gaps at low frequencies, i.e., in case when the wavelengths exceed the dimensions of the periodic cells [START_REF] Liu | Locally resonant sonic materials[END_REF]. One can mention for instance locally resonant phononic crystals, also known as metamaterials, which exhibit effective negative elastic constants in certain frequency bands. The key idea here is to consider periodic arrays of locally resonant devices which are attached to periodic structures. Although the concept has been applied to academic structures such as strings, rods, beams, plates [START_REF] Xiao | Formation and coupling of band gaps in a locally resonant elastic system comprising a string with attached resonators[END_REF][START_REF] Xiao | Longitudinal wave band gaps in metamaterial-based elastic rods containing multi-degree-of-freedom resonators[END_REF][START_REF] Xiao | Flexural wave propagation in beams with periodically attached vibration absorbers: Band-gap behavior and band formation mechanisms[END_REF][START_REF] Wang | Locally resonant band gaps in flexural vibrations of a Timoshenko beam with periodically attached multioscillators[END_REF], its application to complex structures with true periodicities, such as aircraft fuselages, seems to constitute an open research topic.

Analyzing the feasibility of using periodic arrays of locally resonant devices for band gaps generation and passive vibration control of periodic structures is the motivation of the present work. For this purpose, the wave finite element (WFE) method is used, which provides an efficient means for computing the dispersion curves of the waves traveling along periodic structures [START_REF] Zhong | On the direct solution of wave propagation for repetitive structures[END_REF]Williams 1995, Mencik and[START_REF] Mencik | Multi-mode propagation and diffusion in structures through finite elements[END_REF]. Also, the WFE method is used for computing the forced response of complex periodic structures, at a low computational cost [START_REF] Mencik | New advances in the forced response computation of periodic structures using the wave finite element (WFE) method[END_REF][START_REF] Silva | Wave finite element based super-elements for forced response analysis of coupled systems via dynamic substructuring[END_REF].

The rest of the paper is organized as follows. In Section 2, the basics of the WFE method are recalled regarding the calculation of wave modes in periodic structures. Also, the strategy for computing the dispersion curves and the frequency response functions (FRFs) of periodic structures are presented. In Section 3, a comprehensive numerical analysis is carried out regarding a 1D homogeneous rod which is connected to periodic arrays of harmonic oscillators, or a single vibration absorbers. Additional numerical experiments are finally brought in Section 4 regarding a 3D fuselage-like structure equipped with a periodic array of locally resonant devices.

Modeling of periodic structures

WFE method

The WFE method focuses on calculating wave modes traveling along periodic structures, and further, the forced response of those structures by means of wave-based matrix systems of small size, which can be computed in a very fast way. Those structures are periodic in the sense that their FE mesh is composed of identical substructures/cells (with a length of Δ) which are assembled to each other along a certain straight direction (say, axis x), see Fig. 1. The computation of the wave modes involves considering the FE model of one single substructure [START_REF] Mencik | Multi-mode propagation and diffusion in structures through finite elements[END_REF]. One of the advantages of this numerical approach is that, as opposed to analytical approaches, periodic structures composed of complex and heterogeneous substructures can be taken into account.

Within the framework of band gap analysis, a periodic structure is designed so that each substructure contains a locally resonant device, as shown in Fig. 1. By considering the FE model of one substructure, a dynamic equilibrium equation is expressed in matrix form, as follows

[ ] , ) i 1 ( 2 F q K M = + + - h w (1)
where ω is the angular frequency; M, K and η are the mass matrix, stiffness matrix and loss factor is the condensed dynamic stiffness matrix of the substructure (B denoting the DOFs on the left and right boundaries). It is worth pointing out that the left and right boundaries of the substructures are meshed in the same way, i.e., with the same number n of DOFs. This particularly means that the vectors q L , q R , F L and F R have the same size, i.e, n×1 By considering two consecutive substructures k-1 and k, a state vector representation can be proposed as follows , ,

) 1 ( R ) ( R ) 1 ( L ) ( L - - = = k k k k Su u Su u (4)
where

[ ] T T T L L L F q u - = and [ ] T T T R R R F q u =
are 2n×1 state vectors. Also, S is a 2n×2n symplectic transfer matrix expressed by where β j has the meaning of wavenumbers, while the eigenvectors ϕ j are to be understood as wave shapes. The eigenproblem which needs to be solved for expressing μ j and ϕ j is written as follows ,

j j j f m f = S (6)
Notice that the wave shapes ϕ j are to be partitioned into 1 ń vectors of displacements/rotations ϕ qj and forces/moments ϕ Fj , which provide spatial descriptions for the wave motion over the substructure boundaries. The eigensolutions (μ j , ϕ j ) are referred to as the wave modes of the structure. In fact, there are twice as many wave modes as the number of DOFs used to discretize each substructure boundary. Due to the symplectic property of the matrix S, the eigensolutions appears in pairs as (μ j ,1/μ j ), which correspond to n right-going and n left-going wave modes. As some damping is always present in real systems, the magnitudes of the eigenvalues μ j are different from one. Hence the wave modes can be classified into right-going modes for which the magnitudes of μ j are less than one, and left-going modes for which the magnitudes of μ j are greater than one. By convention, those left-going and right-going wave modes are denoted as

{ } n j j j , , 1 ) , ( K = f m and { } n j j j , , 1 ) , ( K = ¢ ¢ f m , where { }. , , 1 / 1 n j j j K Î " = ¢ m m (7)
Also, there exist analytical relations which link the wave shapes ϕ′ j and ϕ j as follows

{ }, , , 1 , F F q q n j j j j j K Î " - = ¢ = ¢ f f f f R R (8)
where R relates a symmetry transformation matrix [START_REF] Mencik | On the low-and mid-frequency forced response of elastic structures using wave finite elements with one-dimensional propagation[END_REF]. Notice that Eq. ( 8) is true provided that the periodic structure is composed of symmetric substructures.

Dispersion curves

Within the WFE framework, the dispersion curve of each wave mode is obtained by analyzing the frequency evolution of the wavenumbers β j which are defined so that being the eigenvalues of the matrix S (Eq. ( 5)) and the substructure length, respectively. Tracking each wave mode over the frequency domain appears to be a crucial task aiming at providing a deep physical insight of the dynamic behavior of periodic structures. Indeed, it will be shown later in this paper that the dispersion curves can be used to predict the far-field dynamic behavior of periodic structures with a large number of substructures. The issue when assessing the frequency evolution of β j (or μ j ) consists in tracking each wave mode over the frequency domain, it being understood that many wave modes are to be computed at several discrete frequencies and there exists a priori no direct connection to link these modes between two consecutive frequencies. In other words, a given mode m defined at a frequency ω i may not match with the mode m defined at the previous frequency ω i-1 .

The wave mode tracking procedure that is considered here makes use of the symplectic orthogonality properties of the wave modes [START_REF] Zhong | On the direct solution of wave propagation for repetitive structures[END_REF] between two consecutive discrete frequencies ω i-1 and ω i spaced by a small step Δω. The procedure can be stated as follows [START_REF] Mencik | On the low-and mid-frequency forced response of elastic structures using wave finite elements with one-dimensional propagation[END_REF]. Given two wave modes l and m verifying μ l =1/ μ m at the frequency ω i-1 , the wave mode m at the frequency ω i is chosen so that

( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ï þ ï ý ü ï î ï í ì = - - - - i k i k i l T i l k i m i m i l T i l w f w f w f w f w f w f w f w f J J 1 1 1 1 max (9) where . ú û ù ê ë é - = n n n n 0 I I 0 J
(10)

Forced response computation

Predicting the frequency behavior of a bounded periodic structure, composed of N substructures and subjected to arbitrary boundary conditions on its left and right ends, consists in expanding the vectors of displacements/rotations and forces/moments, on the left or right boundary of a substructure k, using wave bases

n j j n j j , , 1 q , , 1 q } { } { L L U = = ¢ f f and n j j n j j , , 1 F , , 1 F } { } { L L U = = ¢ f f (Mencik 2014) , , , 1 , q q ) ( R 1 q 1 q ) ( L N k k N k k k N k k K = ¢ ¢ + = ¢ ¢ + = - + - - Q μ Φ Q μ Φ q Q μ Φ Q μ Φ q (11) , , , 1 , F F ) ( R 1 F 1 F ) ( L N k k N k k k N k k K = ¢ ¢ + = ¢ ¢ + = - - + - - Q μ Φ Q μ Φ F Q μ Φ Q μ Φ F ( 12
)
where Φ q , Φ′ q , Φ F and Φ′ F are n×n matrices expressed as

[ ] n q 1 q q f f L = Φ , [ ] n q 1 q q f f ¢ ¢ = ¢ L Φ , [ ] n F 1 F F f f L = Φ and [ ] n F 1 F F f f ¢ ¢ = ¢ L Φ . Also, [ ] T n Q Q L 1 = Q and [ ] T n Q Q ¢ ¢ = ¢ L 1 Q
denote the vectors of wave amplitudes at the left and right ends of the whole periodic structure, respectively. Finally, μ is the n×n diagonal matrix defined so that

{ } n j j , , 1 diag K = = m μ
where it is understood that the magnitudes j m are less than one (see above Eq. ( 7)).

By considering Eqs. ( 11) and ( 12), the condensed dynamic stiffness matrix of the whole periodic structure can be expressed in a straightforward way, as follows. Denote as ) (s D the dynamic stiffness matrix of the periodic structure which is condensed on its left and the right ends. Thus, the dynamic equilibrium equation of the periodic structure can be written as follows .

) ( R ) 1 ( L ) ( R ) 1 ( L ) ( ú û ù ê ë é = ú û ù ê ë é N N s F F q q D ( 13
)
In order to express the dynamic stiffness matrix

) (s D , the vectors of nodal displacements/rotations ) , ( ) ( R ) 1 ( L N q q
and forces/moments ) , (

) ( R ) 1 ( L N F F
on the left and right ends of the periodic structure are expressed by means of Eqs. ( 11) and (12) as follows

, q q q q ) ( R ) 1 ( L ú û ù ê ë é ¢ ú û ù ê ë é ¢ ¢ = ú û ù ê ë é Q Q Φ μ Φ μ Φ Φ q q N N N (14) and . F F F F ) ( R ) 1 ( L ú û ù ê ë é ¢ ú û ù ê ë é ¢ ¢ - - = ú û ù ê ë é Q Q Φ μ Φ μ Φ Φ F F N N N (15)
To derive the condensed dynamic stiffness matrix D (s) of the periodic structure, Eq. ( 14) is left multiplied by the following matrix ,

1 q 1 q ú û ù ê ë é ¢ - - Φ 0 0 Φ (16) which yields . q 1 q q 1 q ) ( R ) 1 ( L 1 q 1 q ú û ù ê ë é ú û ù ê ë é ¢ ¢ = ú û ù ê ë é ú û ù ê ë é ¢ * - - - - Q Q I μ Φ Φ μ Φ Φ I q q Φ 0 0 Φ N N N (17)
The motivation behind the use of the matrix ( 16) is to make the matrix occurring on the right hand side of Eq. ( 17) well-conditioned, as explained in [START_REF] Mencik | On the low-and mid-frequency forced response of elastic structures using wave finite elements with one-dimensional propagation[END_REF]. As it turns out, the vectors of wave amplitudes Q and Q′ can be expressed from Eq. ( 17) as .

) ( R ) 1 ( L 1 q 1 q 1 q 1 q q 1 q ú û ù ê ë é ú û ù ê ë é ¢ ú û ù ê ë é ¢ = ú û ù ê ë é - - - - * - * N N N q q Φ 0 0 Φ I μ Φ Φ μ Φ Φ I Q Q (18)
The condensed dynamic stiffness matrix D (s) of the periodic structure, Eq. ( 13), is readily derived from Eqs. ( 15) and ( 18), as follows .

1 q 1 q 1 q 1 q q 1 q F F F F ) ( ú û ù ê ë é ¢ ú û ù ê ë é ¢ ¢ ú û ù ê ë é ¢ ¢ - - = - - - - - Φ 0 0 Φ I μ Φ Φ μ Φ Φ I Φ μ Φ μ Φ Φ D N N N N s (19)
The determination of the vector of displacements/rotations [ ]

T T N T ) ( R ) 1 ( L q q
can be achieved in the same way as in the conventional FE procedure, i.e., by considering the inverse of the matrix D (s) , see Eq. ( 13). It is worth emphasizing that the consideration of assemblies composed of several periodic structures and other FE components, that may be connected to each other in arbitrary ways, can be addressed without any additional issues through FE procedures [START_REF] Silva | On the forced harmonic response of coupled systems via a WFE-based super-element approach[END_REF][START_REF] Silva | Wave finite element based super-elements for forced response analysis of coupled systems via dynamic substructuring[END_REF].

Comparisons between local and global attenuation effects

Vibration absorbers are mass-spring or mass-spring-damper systems which are widely used in the field of vibration control engineering. Since the early development of the theory of dynamic vibration absorbers [START_REF] Ormondroyd | The theory of the dynamic vibration absorber[END_REF], they have been applied in a wide range of applications. In many situations, a single vibration absorber punctually attached to one structure is considered so as to avoid the vibrations induced by some vibrational modes. The use of periodic arrays of these locally resonant devices, along structures, is very recent [START_REF] Thompson | A continuous damped vibration absorber to reduce broad-band wave propagation in beams[END_REF]. It has been first investigated by [START_REF] Kashina | Waveguide vibration reduction of longitudinal and flexural modes by means of a multielement structure of resonators[END_REF] who have proposed the use of a set of undamped resonators to reduce the longitudinal and flexural motions in beams and plates. This idea has received new attention with the publication of the work of [START_REF] Liu | Locally resonant sonic materials[END_REF], in which emphasis is placed on exploring the impact of using periodic arrays of resonant devices on wave propagation in artificial periodic composites known as phononic crystals. In such systems, band gaps can be optimized by tuning the resonator's natural frequency to the control procedure's targeted one.

In this section, the decrease of the vibration levels produced by adding a periodic array of resonant devices to a 1D homogeneous rod is analyzed by means of the WFE method. The efficiency in using such devices is highlighted through comparisons with the vibration levels produced by a single vibration absorber. While a single vibration absorber produces local vibration attenuation, a periodic array of small absorbers yields global vibration control. These phenomena are illustrated here regarding a simple 1D homogeneous rod structure having longitudinal vibrations (see Fig. 2), whose material and geometric properties are listed in Table 1. The left and right ends of the rod are, respectively, excited by a unitary axial harmonic force, and free. Here, the structure is either connected to a periodic array of nine spring-mass systems (Fig. 2(b)), or connected to one single spring-mass system (Fig. 2(c)). The dispersion curves and the FRFs of the uncoupled and coupled structure are computed using the WFE method. Within this framework, the substructure being considered (Section 2.1) represents either a single rod cell, or a rod cell with a spring-mass resonant system (see Fig. 2(b)). Notice that each substructure is modelled by means of the spectral element method (SEM) which provides an exact analytical expression of the condensed dynamic stiffness matrix D * (Eq. ( 3)) [START_REF] Doyle | Wave Propagation in Structures, 2nd Edition[END_REF][START_REF] Xiao | Longitudinal wave band gaps in metamaterial-based elastic rods containing multi-degree-of-freedom resonators[END_REF]. Within the analytical framework, the displacement field of a rod cell is described in a standard way by means of two longitudinal waves, i.e., two propagating ones traveling in opposite directions along a onedimensional system. The consideration of mass-spring-damper systems is simply achieved through basic kinematic compatibility and force equilibrium equation. To begin with, a periodic array of nine resonant devices (Fig. 2(b)) is considered so that their resonance frequency is 0.5 f ref = 790.6 Hz. Two test cases are considered which involve resonant devices whose total mass represents, respectively, 44% and 14% of that of the rod with the resonant devices. The related dispersion curves of the longitudinal wave traveling in the coupled rod are displayed in Fig. 3, along with the dispersion curve that concerns the uncoupled rod. Here the non-dimensional wavenumber βΔ is analyzed as a function of the non-dimensional frequency f/f ref , where f=2πω. As it can be seen in Fig. 3(a), the longitudinal wave in the rod without resonant devices is always propagative, as expected. Notice that the imaginary part of βΔ, although small, is not equal to zero due to damping effects. Regarding the resonant devices, three band gaps (i.e., attenuation zones) occur in Figs. 3(b)-(c interesting feature when using a periodic array of resonant devices is that it generates band gaps with global attenuation of the vibration levels (i.e., in frequency and space) along the structure. By comparing resonant devices of different masses to each other, one may notice that higher mass fractions produce higher levels of attenuation and larger attenuation zones (see Figs. Consider now the 2D maps of the FRFs of the rod with one single vibration absorber. The results are shown in Fig. 5 for different mass fractions (8%, 44% and 14%) and positions of the vibration absorber (x=L/2 and x=Δ/2). Those results are issued from a conventional analytical approach that consists in modeling each rod component (before and after the resonator) by means of a spectral element (Doyle 1997) -i.e., by modeling its dynamic stiffness matrix by means of two longitudinal waves propagating in right and left directions, respectively-and modeling the resonator by means of a spring-mass system. As it can be seen, band gaps formed by Bragg scattering mechanisms-i.e., for f / f ref = 1 and f / f ref = 2 -are not observed in any case, as opposed to the periodic array of resonant devices even though in case when the added masses involved are the same. Also, the attenuation zone remains restricted to a short frequency bandwidth around f / f ref = 0.5. Notice that the attenuation zone is local, i.e., concentrated at anti-resonance points. Moreover, it depends on the locations of the vibration absorber and excitation point on the structure. In contrast, due to their global attenuation effect, periodic arrays of resonant devices have the potentiality to be used for excitations located at arbitrary locations, i.e., which are not necessarily known. This feature is clearly highlighted in Fig. 6 in which the spatial variations in the reduction of the vibration levels produced, respectively, by the vibration absorber and a periodic array of resonant devices of same masses, are plotted at the target resonance frequency f = 0.5 f ref .

Band gap generation in a 3D aircraft fuselage-like structure

Problem description

In this section, one aims at showing that important reduction of the vibration levels of complex periodic structures can be achieved through the use of a lightweight periodic array of resonant devices whose parameters are adequately tuned. The periodic structure that is considered here represents a coarse approximation of a 3D aircraft fuselage, consisting in a thin cylindrical shell with axial stiffeners (stringers) and circumferential stiffeners (frames), see Fig. 7. Those stiffeners are periodically distributed over the cylinder, two consecutive circumferential (resp. axial) stiffeners being spaced of 0.4 m (resp. of 11.25°). The material and geometric characteristics of the stiffened cylindrical shell are listed in Table 2.

Within the framework of the WFE method (Section 2.1), a substructure is considered as shown in Fig. 7 It is excited by unitary harmonic radial and axial forces, as well as moments, which are applied on the left end at θ=45°, and it is free from excitations on its right end. Periodic arrays of resonant devices are considered so as to attenuate the vibration levels throughout the whole structure. The motivation behind the present study is to design lightweight resonant devices which may be effectively attached to a real aircraft fuselage. In the present study, each resonant device is modeled as a flexural Euler beam with a concentrated mass on one end and whose other end is connected to the stringers and frames of the periodic structure at the other end, as shown in Fig. 7(a). The first flexural resonances of those resonators are assigned to be 197 Hzwhich corresponds to one resonance frequency of the whole periodic structure -and 275 Hz, respectively. The characteristics of the resonators are listed in Table 3. Notice that the added mass induced by the periodic array of these resonant devices represents about 12.5% of the total mass to the original periodic structure, i.e., it remains small. Within the WFE framework, a substructure is thus considered which takes into account the FE models of the resonant beams along with those of the cylinder, frames and stringers. This yields 3,114 DOFs for modeling the substructure, with 246 DOFs over the left/right boundary, and 2,622 internal DOFs. By considering the symmetry boundary conditions, this yields n=240 right-going and left-going wave modes which are to be computed by means of the WFE method.

Numerical analysis and discussion

The dispersion curves of the periodic structure are assessed by means of the WFE method as shown in Fig. 8. Regarding Fig. 8(a) (without resonant devices), it should be noticed that one particular wave mode which propagates along the 3D aircraft fuselage-like structure becomes evanescent between 230 Hz and 240 Hz, which characterizes a band gap phenomenon. This might be explained by the fact that the circumferential stiffeners (frames) exhibit some flexural vibration modes in this frequency band. Regarding Fig. 8(b) (with resonant devices), it is seen that the added periodic array of resonant devices yields large evanescent parts around the target frequencies 197 Hz and 275 Hz, which means that the propagation of waves is strongly attenuated.

Also, 2D maps of the FRFs of the periodic structure are shown in Fig. 9. Here, the root mean Vibration decrease is also observed around the other target frequency 275 Hz, even though it does not correspond to a resonance frequency of the periodic structure. The attenuation between the FRFs of the controlled and uncontrolled structure are plotted in Fig. 10 in dB. Such a decrease in the vibration levels appears to be widely spread over the length of the structure, even near the excited cross-section. This phenomenon is clearly highlighted in Figs. 11-12, which highlight the spatial distributions of the displacements at 197 Hz and 275 Hz, respectively. This provides a clear evidence that the magnitude of the whole displacement field of the structure can be strongly decreased.

It is worth pointing out that the proposed WFE approach makes use of full wave mode bases for modeling periodic structures, i.e., it does not invoke any basis truncation or reduction of the interface DOFs. In other words, the WFE method provides, in theory, the same level of accuracy as the FE method for modeling periodic structures. However, the WFE method can greatly improve the computational efficiency of the FE method because it involves modeling a single substructure. Further evidence of the efficiency of the WFE method is shown by comparing the sizes of the FE and WFE matrix systems. Considering the present case, the corresponding FE model of the structure with resonant devices would have involved 114,966 DOFs, compared to 3,144 DOFs with the WFE method. This fully gives credit to the proposed approach.

Conclusions

In this work, a passive vibration control of periodic structures has been analyzed which involves considering periodic arrays of simple resonant devices. The WFE method has been used for this purpose which constitutes a fast and efficient means for assessing the dispersion curves and FRFs of periodic structures. It has been shown that the consideration of lightweight periodic arrays of resonant devices yields large and global decreases of the vibration levels, even at low frequencies. Such resonant devices have proved to be relevant for treating large-sized complex 3D structures such as aircraft fuselages. Future work may involve those resonant devices for passive control of internal acoustics in fluid-filled periodic structures.

Fig. 1

 1 Fig. 1 Schematics of a one-dimensional periodic structure composed of N substructures (a) and of a particular substructure (b)

  's theorem, the eigensolutions of the matrix S i.e., eigenvalues μ j and right eigenvectors ϕ j are referred to as waves traveling along the periodic structure. More precisely, what states Bloch's theorem is that the eigenvalues are expressed as

  Fig. 2 1D homogeneous rod: (a) without resonant devices; (b) with a periodic array of resonant devices; (c) with a single vibration absorber

  Fig. 4 Continued

  4(b)-(c)).

  (b), which consists in a quarter of the stiffened cylindrical shell with a length Δ=0.4 m. Symmetry boundary conditions are considered so as to model the structure along the whole circumference. Regarding the substructure, the cylinder and stiffeners are meshed by means of model of the 3D aircraft fuselage-like structure: (a) unit cell; (b) substructure; (c) whole structure Table 2 Characteristics of the stiffened cylindrical shell. Mindlin shell elements with six DOFs per node, i.e., translations in the x, y and z directions and rotations about the same axes. Here, the total number of DOFs used to discretize the substructure is 1,578, with 246 DOFs over the left/right boundary and 1,086 internal DOFs. The whole periodic structure is modeled by means of N=40 identical substructures connected along the x-direction, as shown in Fig. 7(c).

Fig. 8

 8 Fig. 8 Dispersion curves of the waves traveling along the 3D aircraft fuselage-like structure: (a) without resonant devices; (b) with the periodic array of resonant devices. (blue color) real part of β j Δ; (red color) imaginary part of β j Δ

Acknowledgments

The authors are grateful to the government funding agency Fundação de Amparo à Pesquisa de São Paulo, FAPESP, for the financial support provided for the present research work through processes 2010/17317-9 and 2013/23542-3. The authors are also thankful to the Program "Cátedras Francesas na UNICAMP" and the National Research Council, CNPq, through process 304482/2011-0.