Alfredo Buttari

Fine-grained multithreading for the multifrontal QR factorization of sparse matrices *

The advent of multicore processors represents a disruptive event in the history of computer science as conventional parallel programming paradigms are proving incapable of fully exploiting their potential for concurrent computations. The need for different or new programming models clearly arises from recent studies which identify fine-granularity and dynamic execution as the keys to achieve high efficiency on multicore systems. This work presents an approach to the parallelization of the multifrontal method for the QR factorization of sparse matrices specifically designed for multicore based systems. High efficiency is achieved through a fine-grained partitioning of data and a dynamic scheduling of computational tasks relying on a dataflow parallel programming model. Experimental results show that an implementation of the proposed approach achieves higher performance and better scalability than existing equivalent software.

Introduction

The QR factorization is the method of choice for the solution of least-squares problems arising from a vast field of applications including, for example, geodesy, photogrammetry and tomography [START_REF] Rice | PARVEC workshop on very large least squares problems and supercomputers[END_REF][START_REF] Björck | Numerical methods for Least Squares Problems[END_REF].

The cost of the QR factorization of a sparse matrix, as well as other factorizations such as Cholesky or LU, is strongly dependent on the fill-in generated, i.e., the number of nonzero coefficients introduced by the factorization. Although the QR factorization of a dense matrix can attain very high efficiency because of the use of Householder transformations (see Schreiber and Van Loan [START_REF] Schreiber | A storage-efficient WY representation for products of Householder transformations[END_REF]), early methods for the QR factorization of sparse matrices were based on Givens rotations with the objective of reducing the fill-in. One such method was proposed by Heath and George [START_REF] George | Solution of sparse linear least squares problems using Givens rotations[END_REF], where the fill-in is minimized by using Givens

The Multifrontal QR Factorization

The multifrontal method was first introduced by Duff and Reid [START_REF] Duff | The multifrontal solution of indefinite sparse symmetric linear systems[END_REF] as a method for the factorization of sparse, symmetric linear systems and, since then, has been the object of numerous studies and the method of choice for several, highperformance, software packages such as MUMPS [START_REF] Amestoy | MUMPS: a general purpose distributed memory sparse solver[END_REF] and UMFPACK [START_REF] Davis | Algorithm 832: UMFPACK V4.3 -an unsymmetricpattern multifrontal method[END_REF]. At the heart of this method is the concept of an elimination tree, introduced by Schreiber [START_REF] Schreiber | A new implementation of sparse Gaussian elimination[END_REF] and extensively studied and formalized later by Liu [START_REF] Liu | The role of elimination trees in sparse factorization[END_REF]. This tree graph describes the dependencies among computational tasks in the multifrontal factorization. The multifrontal method can be adapted to the QR factorization of a sparse matrix thanks to the fact that the R factor of a matrix A and the Cholesky factor of the normal equation matrix A T A share the same structure under the hypothesis that the matrix A is Strong Hall (for a definition of this property see, for example, Bjork's book [3, page 234]). Based on this equivalence, the elimination tree for the QR factorization of A is the same as that for the Cholesky factorization of A T A. In the case where the Strong Hall property does not hold, the elimination tree related to the Cholesky factorization of A T A can still be used although the resulting QR factorization will perform more com-putations and consume more memory than what is really needed; alternatively, the matrix A can be permuted to a Block Triangular Form (BTF) where all the diagonal blocks are Strong Hall.

In a basic multifrontal method, the elimination tree has n nodes, where n is the number of columns in the input matrix A, each node representing one pivotal step of the QR factorization of A. Every node of the tree is associated with a dense frontal matrix that contains all the coefficients affected by the elimination of the corresponding pivot. The whole QR factorization consists in a topological order (i.e., bottom-up) traversal of the tree where, at each node, two operations are performed:

• assembly: a set of rows from the original matrix is assembled together with data produced by the processing of child nodes to form the frontal matrix. This operation simply consists in stacking these data one on top of the other: the rows of these data sets can be stacked in any order but the order of elements within rows has to be the same. The set of column indices associated with a node is, in general, a superset of those associated with its children;

• factorization: one Householder reflector is computed and applied to the whole frontal matrix in order to annihilate all the subdiagonal elements in the first column. This step produces one row of the R factor of the original matrix and a complement which corresponds to the data that will be later assembled into the parent node (commonly referred to as a contribution block). The Q factor is defined implicitly by means of the Householder vectors computed on each front; the matrix that stores the coefficients of the computed Householder vectors, will be referred to as the H matrix from now on.

In practical implementations of the multifrontal QR factorization, nodes of the elimination tree are amalgamated to form supernodes. The amalgamated pivots correspond to rows of R that have the same structure and can be eliminated at once within the same frontal matrix without producing any additional fill-in in the R factor. The elimination of amalgamated pivots and the consequent update of the trailing frontal submatrix can thus be performed by means of efficient Level-3 BLAS routines through the WY representation [START_REF] Schreiber | A storage-efficient WY representation for products of Householder transformations[END_REF]. Moreover, amalgamation reduces the number of assembly operations increasing the computations-to-communications ratio which results in better performance. The amalgamated elimination tree is also commonly referred to as assembly tree.

Figure 1 shows some details of a sparse QR factorization. The factorized matrix 1 is shown on the left part of the figure where a letters are used for representing the original matrix coefficients and dots for representing the fill-in introduced by the factorization. On the top-right part of the figure, the structure of the resulting R factor is shown. The elimination/assembly tree is, instead Figure 1: Example of multifrontal QR factorization. On the left side the matrix with the original coefficient represented as a and the fill-in coefficients introduced by the factorization as dots. On the upper-right part, the structure of the resulting R factor. On the right-bottom part the elimination tree; the dashed boxes show how the nodes are amalgamated into supernodes. reported in the bottom-right part: dashed boxes show how the nodes can be amalgamated into supernodes with the corresponding indices denoted by bigger size numbers. The amalgamated nodes have the same row structure in R modulo a full, diagonal block. It has to be noted that in practical implementations the amalgamation procedure is based only on information related to the R factor and, as such, it does not take into account fill-in that may eventually appear in the H matrix; for example, the amalgamation of the pivots 3 and 4 generates the h 10,3 and h 11,3 fill-in coefficients although no extra fill-in appears in the R factor. The supernode indices are reported on the left part of the figure in order to show the pivotal columns eliminated within the corresponding supernode and on the top-right part to show the rows of R produced by the corresponding supernode factorization.

In order to reduce the operation count of the multifrontal QR factorization, two optimizations are commonly applied:

1. once a frontal matrix is assembled, its rows are sorted in order of increasing index of the leftmost nonzero. The number of operations can thus be reduced, as well as the fill-in in the H matrix, by ignoring the zeroes in the bottom-left part of the frontal matrix;

2. the frontal matrix is completely factorized. Despite the fact that more Householder vectors have to be computed for each frontal matrix, the overall number of floating point operations is lower since frontal matrices are smaller. This is due to the fact that contribution blocks resulting from the complete factorization of frontal matrices are smaller.

h r r r 9 h h c c 5 h h h c 10 h h h 11 h h h 2 h h h 6 h h h 3 h

Figure 2: Assembly and factorization of the frontal matrices associated with supernodes 1, 2 and 3 in Figure 1. Coefficients from the original matrix are represented as a, those in the resulting R and H matrices as r and h, respectively, the fill-in coefficients as dots and the coefficients in the contribution blocks as c with a subscript to specify the supernode they belong to. Pivotal columns are shaded in gray.

Figure 2 shows the assembly and factorization operations for the supernodes 1, 2 and 3 in Figure 1 when these optimization techniques (referred to as Strategy 3 in Amestoy et al. [START_REF] Amestoy | Multifrontal QR factorization in a multiprocessor environment[END_REF]) are applied. Note that, because supernodes 1 and 2 are leaves of the assembly tree, the corresponding assembled frontal matrices only include coefficients from the matrix A. The contribution blocks resulting from the factorization of supernodes 1 and 2 are appended to the rows of the input A matrix associated with supernode 3 in such a way that the resulting, assembled, frontal matrix has the staircase structure shown in Figure 2 (bottom-middle). Once the front is assembled, it is factorized as shown in Figure 2 (bottom-right).

A detailed presentation of the multifrontal QR method, including the optimization techniques described above, can be found in Amestoy et al. [START_REF] Amestoy | Multifrontal QR factorization in a multiprocessor environment[END_REF].

The multifrontal method can achieve very high efficiency on modern computing systems because all the computations are arranged as operations on dense matrices; this reduces the use of indirect addressing and allows the use of efficient Level-3 BLAS routines which can achieve a considerable fraction of the peak performance of modern computing systems.

The factorization of a sparse matrix is preceded by a preprocessing phase, commonly referred to as the analysis phase, where a number of (mostly symbolic) operations are performed on the matrix such as row and column permutations to reduce the amount of fill-in and the determination of the elimination tree or the symbolic factorization to estimate the amount of memory needed during the factorization phase. The cost of the analysis phase is O(|A| + |R|) and can be considered negligible with respect to the matrix factorization.

The rest of this paper is based on the assumption that the analysis phase is already performed, and thus it only focuses on the factorization; specifically, it is assumed that a fill-reducing permutation of the input matrix and the corresponding assembly tree have been computed.

Fine-grained, asynchronous multithreading

Sparse computations are well known for being hard to parallelize on sharedmemory, multicore systems. This is due to the fact that the efficiency of many sparse operations, such as the sparse matrix-vector product, is limited by the speed of the memory system. This is not the case for the multifrontal method; since computations are performed as operations on dense matrices, a favorable ratio between memory accesses and computations can be achieved which reduces the utilization of the memory system and opens opportunities for multithreaded, parallel execution. This property, resulting from the use of Level-3 BLAS operations (i.e., matrix-matrix operations), is commonly known as the surface-to-volume [START_REF] Dongarra | Numerical Linear Algebra for High-Performance Computers[END_REF] property because, for a problem of size n (the size of a frontal matrix, for the multifrontal algorithm), n 3 operations are performed on n 2 data.

In a multifrontal factorization, parallelism is exploited at two levels:

• tree-level parallelism: computations related to separate branches of the assembly tree are independent and can be executed in parallel (see Liu [START_REF] Liu | The role of elimination trees in sparse factorization[END_REF]Proposition 10.1]);

• node-level parallelism: if the size of a frontal matrix is big enough, its factorization can be performed in parallel by multiple threads.

The classical approach to shared-memory parallelization of QR multifrontal solvers [START_REF] Matstoms | Parallel sparse QR factorization on shared memory architectures[END_REF][START_REF] Amestoy | Multifrontal QR factorization in a multiprocessor environment[END_REF][START_REF] Davis | Algorithm 915, SuiteSparseQR: Multifrontal multithreaded rank-revealing sparse QR factorization[END_REF] is based on a complete separation of the two sources of concurrency described above. The node parallelism is delegated to multithreaded BLAS libraries and only the tree parallelism is handled at the level of the multifrontal factorization. This is commonly achieved by means of a task queue where a task corresponds to the assembly and factorization of a front. A new task is pushed into the queue as soon as it is ready to be executed, i.e., as soon as all the tasks associated with its children have been completed. Threads keep polling the queue for tasks to perform until all the nodes of the tree have been processed.

Although this approach works reasonably well for a limited number of cores or processors, it suffers scalability problems mostly due to two factors:

• separation of tree and node parallelism: the degree of concurrency in both types of parallelism changes during the bottom-up traversal of the tree; fronts are relatively small at leaf nodes of the assembly tree and grow bigger towards the root node. On the other hand, tree parallelism provides a high level of concurrency at the bottom of the tree and only a little at the top part where the tree shrinks towards the root node. Since the node parallelism is delegated to an external multithreaded BLAS library, the number of threads dedicated to node parallelism and to tree parallelism has to be fixed before the execution of the factorization. Thus, a thread configuration that may be optimal for the bottom part of the tree will result in a poor parallelization of the top part and vice versa. Although some recent parallel BLAS libraries (for example, Intel MKL) allow to change the numbers of threads dynamically at run-time, it would require an accurate performance modeling and a rigid thread-to-front mapping in order to keep all the cores working at any time. Relying on some specific BLAS library could, moreover, limit the portability of the code.

• synchronizations: the assembly of a front is an atomic operation. This inevitably introduces synchronizations that limit the concurrency level in the multifrontal factorization; most importantly, it is not possible to start working on a front until all of its children have been fully factorized.

The limitations of the classical approach discussed above can be overcome by employing a different parallelization technique based on fine granularity partitioning of data and operations combined with a dataflow model for the scheduling of tasks, as described in the following subsections. This approach was already applied to dense matrix factorizations by Buttari et al. [START_REF] Buttari | A class of parallel tiled linear algebra algorithms for multicore architectures[END_REF] and extended to the supernodal Cholesky factorization of sparse matrices by Hogg et al. [START_REF] Hogg | A DAG-based sparse Cholesky solver for multicore architectures[END_REF].

Fine-grained computational tasks definition

In order to handle both tree and node parallelism in the same framework, a block-column partitioning of the fronts is applied as shown in Figure 3 (left) and five elementary operations defined:

1. activate: the activation of a frontal matrix corresponds to computing its structure (row/column indices, staircase structure, etc.) and allocating the memory needed for it;

2. panel: this operation amounts to computing the QR factorization of a block-column; Figure 3 (middle) shows the data modified when the panel operation is executed on the first block-column;

3. update: updating a block-column with respect to a panel corresponds to applying to the block-column the Householder reflectors resulting from the panel reduction; The multifrontal factorization of a sparse matrix can thus be defined as a sequence of tasks, each task corresponding to the execution of an elementary operation of the type described above on a block-column or a front. The tasks are arranged in a Directed Acyclic Graph (DAG) such that the edges of the DAG define the dependencies among tasks and thus the relative order in which they have to be executed. Figure 4 shows the DAG associated with the subtree defined by supernodes one, two and three for the problem in Figure 1 for the case where the block-columns have size one2 ; the dashed boxes surround all the tasks that are related to a single front.

The dependencies in the DAG are defined according to the following rules (an example of each of these rules is presented in Figure 4 with labels on the edges):

• d1: no other elementary operation can be executed on a front or on one of its block-columns until the front is not activated;

• d2: a block column can be updated with respect to a panel only if the corresponding panel factorization is completed;

Figure 4: The DAG associated with supernodes 1, 2 and 3 of the problem in Figure 1; for the panel, update and assemble operations, the corresponding block-column index is specified. For this example, the block-column size is chosen to be one.

• d3: the panel operation can be executed on block-column i only if it is up-to-date with respect to panel i -1;

• d4: a block-column can be updated with respect to a panel i in its front only if it is up-to-date with respect to the previous panel i -1 in the same front;

• d5: a block-column can be assembled into the parent (if it exists) when it is up-to-date with respect to the last panel factorization to be performed on the front it belongs to (in this case it is assumed that block-column i is up-to-date with respect to panel i when the corresponding panel operation is executed);

• d6: no other elementary operation can be executed on a block-column until all the corresponding portions of the contribution blocks from the child nodes have been assembled into it, in which case the block-column is said to be assembled;

• d7: since the structure of a frontal matrix depends on the structure of its children, a front can be activated only if all of its children are already active;

This DAG globally retains the structure of the assembly tree but expresses a higher degree of concurrency because tasks are defined on a block-column basis instead of a front basis. Moreover, it implicitly represents both tree and node parallelism which allows to exploit both of them in a consistent way. Finally, it removes unnecessary dependencies making it possible, for example, to start working on the assembled block-columns of a front even if the rest of the front is not yet assembled and, most importantly, even if the children of the front have not yet been completely factorized.

Scheduling and execution of tasks

The execution of the tasks in the DAG is controlled by a dataflow model; a task is dynamically scheduled for execution as soon as all the input operands are available to it, i.e., when all the tasks on which it depends have finished. The scheduling of tasks can be guided by a set of rules that prioritize the execution of a task based on, for example,

• data locality: in order to maximize the reuse of data into the different level of the memory hierarchy, tasks may be assigned to threads based on a locality policy [START_REF] Hogg | A DAG-based sparse Cholesky solver for multicore architectures[END_REF];

• fan-out: the fan-out of a task in the DAG defines the number of other tasks that depend on it. Thus, tasks with a higher fan-out should acquire higher priority since their execution generates more concurrency. In the case of the QR method described above, panel factorizations are regarded as higher priority operations over the updates and assemblies.

A scheduling technique was developed aiming at optimizing the reuse of local data in a NUMA system while still prioritizing tasks that have a high fan-out. Although the multifrontal method is rich in Level-3 BLAS operations (mostly in the update tasks), there is still a considerable amount of Level-2 BLAS operations (within the panel tasks) and symbolic or memory ones (the activate, assembly and clean tasks) whose efficiency is limited by the speed of the memory system. As the number of threads participating to the factorization increases, the relative cost of these memory-bound operations grows too high and there are fewer opportunities to hide this cost by overlapping these slow operations with faster ones. In addition, some frontal matrices, especially at the bottom of the tree, may be too small to achieve the surface-to-volume effect in Level-3 BLAS operations. Therefore, in order to improve the scalability, it is important to perform these memory-bound operations as efficiently as possible and this can be achieved by executing each of them on the core which is closest to the data it manipulates. The proposed method is based on a concept of ownership of a front: the thread that performs the activate operation on a front becomes its owner and, therefore, becomes the privileged thread to perform all the subsequent tasks related to that front. By using methods like the "first touch rule" (memory is placed on the NUMA node which generates the first reference) or allocation routines which are specific for NUMA architectures [START_REF] Broquedis | hwloc: A generic framework for managing hardware affinities in hpc application[END_REF], the memory needed for a front can be allocated in the NUMA node which is closest to its owner thread. At the moment no front-to-thread mapping is performed, and thus the ownership of a front is dynamically set at the moment when the front is activated.

The scheduling and execution of tasks is implemented through a system of task queues: each thread is associated with a task queue containing all the executable tasks (i.e., tasks whose dependencies are already satisfied) related to the fronts it owns.

The pseudo-code in Figure 5(left) illustrates the main loop executed by all threads; at each iteration of this loop a thread:

1. checks whether the number of tasks globally available for execution has fallen below a certain value (which depends, e.g., on the number of threads) and, if it is the case, it calls the fill queues routine, described below, which searches for ready tasks and pushes them into the corresponding local queues;

2. picks a task. This operation consists in popping a task from the thread's local queue. In the case where no task is available on the local queue, an architecture aware work-stealing technique is employed, i.e., the thread will try to steal a task from queues associated with threads with which it shares some level of memory (caches or DRAM module on a NUMA machine) and if still no task is found it will attempt to steal a task from any other queue. The computer's architecture can be detected using tools such as hwloc [START_REF] Broquedis | hwloc: A generic framework for managing hardware affinities in hpc application[END_REF].

3. executes the selected task if the pick task routine has succeeded.

The tasks are pushed into the local queues by the fill queues routine whose pseudo-code is shown in Figure 5 (right). At every moment, during the factorization there exists a list of active fronts; the fill queues routine goes through this list looking for ready tasks on each front. Whenever one such task is found, it is pushed on the queue associated with the thread that owns the front. If no task is found related to any of the active fronts, a new ready front (if any) is scheduled for activation; the search for an activable front follows a postorder traversal of the assembly tree, which provides a good memory consumption and temporal locality of data. Because the ownership of such a front is not yet defined, the activate task is pushed on the queue attached to the thread that executes the fill queues routine. Simultaneous access to the same front in the fill queues routine is prevented through the use of locks. An efficient use of tree-level parallelism makes it hard, if not impossible, to follow a postorder traversal of the tree which results in an increased memory consumption with respect to the sequential case [START_REF] Liu | On the storage requirement in the out-of-core multifrontal method for sparse factorization[END_REF][START_REF] Guermouche | Impact of reordering on the memory of a multifrontal solver[END_REF]. It is important to note that the proposed scheduling method tries to exploit node-level parallelism as much as possible and dynamically resorts to tree-level parallelism by activating a new node only when no more tasks are found on already active fronts. This keeps the tree traversal as close as possible to the one followed in the sequential execution and avoids the memory consumption to grow out of control.

Although tasks are always popped from the head of each queue, they can be pushed either on the head or on the tail which allows to prioritize certain tasks. In the implementation described in Section 4, the panel operations are always pushed on the head because the corresponding nodes in the execution DAG have higher fan-out.

The size of the search space for the fill queues routine is, thus, proportional to the number of fronts active, at a given moment, during the factorization. The size of this search space may become excessively large and, consequently, the relative cost of the fill queues routine excessively high, in some cases like, for example, when the assembly tree is very large and/or when it has many nodes of small size. Two techniques are employed in order to keep the number of active nodes limited during the factorization:

• Logical pruning. As the target of this work are systems with only a limited number of cores, extremely large assembly trees provide much more tree-level parallelism than what is really needed. A logical pruning can thus be applied to simplify the tree: all the subtrees whose relative computational weight is smaller than a certain threshold are made invisible to the fill queues routine. When one of the remaining nodes is activated, all the small subtrees attached to it are processed sequentially by the same thread that performs the activation. As shown in Figure 6, this corresponds to identifying a layer in the assembly tree such that all the subtrees below it will be processed sequentially. This layer has to be as high as possible in order to reduce the number of potentially active nodes but low enough to provide a sufficient amount of tree-level parallelism on the top part of the tree.

Figure 6: A graphical representation of how the logical amalgamation and logical pruning may be applied to an assembly tree.

• Tree reordering. Assuming that the assembly tree is processed sequentially following a postorder traversal, the maximum number of fronts active at any time in the subtree rooted at node i, P i is defined as the maximum of two quantities:

1. nc i + 1, where nc i is the number of children of node i. This is the number of active nodes at the moment when i is activated;

2. max j=1,...,nci (j -1 + P j ). This quantity captures the maximum number of active fronts when the children of node i are being processed.

In fact, at the moment when the peak P j is reached in the subtree rooted at the j-th child, all the previous j -1 children of i are still active.

In order to minimize the maximum number of active nodes during the traversal of the assembly tree it is, thus, necessary to minimize, for each node i, the second quantity, which is achieved by sorting all of its children j in decreasing order of P j [START_REF] Liu | On the storage requirement in the out-of-core multifrontal method for sparse factorization[END_REF]. The effect of this reordering on an example tree is illustrated in Figure 7. If the tree is traversed in the order shown in the left part of the figure, P 19 = 10 (the nodes active at the moment when the peak is reached are highlighted with a thick border); instead, if the tree is reordered as in the right part of the figure following the method described above P 19 is equal to four. Although no guarantee is given that a postorder is followed in a parallel factorization, this tree reordering technique still provides excellent results on every problem that has been tested so far. Besides, by reducing the number of active nodes, this reordering also helps in reducing the consumed memory although it will not be optimal in this sense as the actual size of the frontal matrices is not taken into account (the reordering technique for memory consumption minimization is described in Liu and Guermouche et al. [START_REF] Liu | On the storage requirement in the out-of-core multifrontal method for sparse factorization[END_REF][START_REF] Guermouche | Impact of reordering on the memory of a multifrontal solver[END_REF]).

Figure 7: The effect of leaves reordering on the number of active nodes.

Both the tree pruning and reordering are executed during the analysis phase.

Blocking of dense matrix operations

It is obviously desirable to use blocked operations that rely on Level-3 BLAS routines in order to achieve a better use of the memory hierarchy and, thus, better performance. The of blocked operations, however, introduces additional fill-in in the Householder vectors due to the fact that the staircase structure of the frontal matrices cannot be fully exploited. It can be safely said that it is always worth paying the extra cost of this additional fill-in because the overall performance will be drastically improved by the high efficiency of Level-3 BLAS routines; nonetheless it is important to choose the blocking value that gives the best compromise between number of operations and efficiency of the BLAS. This blocking size, which defines the granularity of computations, has to be chosen with respect to the block-columns size used for partitioning the frontal matrices, which defines the granularity of parallel tasks.

Figure 8: The effect of internal blocking on the generated fill-in. The light gray dots show the frontal matrix structure if no blocking of operations is applied whereas the dark gray dots show the additional fill-in introduced by blocked operations.

Figure 8 shows as dark gray dots the extra fill-in introduced by the blocking of operations (denoted as ib for internal blocking) with respect to the partition- ing size nb on an example frontal matrix. Details on how to choose these bock sizes will be provided in Section 4.1.

Experimental results

The method discussed in Section 3 was implemented in a software package called qr mumps3 (or qrm for the sake of brevity). The code is written in Fortran2003 and OpenMP is the technology chosen to implement the multithreading. Although there are many other technologies for multithreaded programming (e.g., pThreads, Intel TBB, Cilk or SMPSS), OpenMP offers the best portability since it is available on any relatively recent system. The experiments were run on a set of fifty matrices from the UF Sparse Matrix Collection [START_REF] Davis | The university of Florida sparse matrix collection[END_REF]; they were chosen from the complete set of over and under-determined problems by excluding those that are rank-deficient (because qr mumps cannot currently handle them), those that are too small to evaluate the scalability or that are too big for being factorized on the computer platforms described below. To these, a large matrix from a meteorology application from the HIRLAM4 research program was added for a total of 51 test matrices. In the case of under-determined systems, the transposed matrix is factorized, as it is commonly done to find the minimum-norm solution of a problem. All of the results presented below were produced without storing the H matrix in order to extend the test set to very large matrices that couldn't otherwise be factorized on the available computers. Figure 9 shows, for all these matrices, the operation count obtained when a COLAMD [START_REF] Davis | A column approximate minimum degree ordering algorithm[END_REF] fill-reducing column permutation is applied and the corresponding qrm sequential factorization time on the AMD system on the left and the right side, respectively; data are sorted in increasing order in each subfigure. Out of these 51 problems, a subset of 12 matrices, described in Table 1, was selected and analyzed in details in order to highlight the effectiveness of the proposed algorithms and methods.

#

Experiments were run on two architectures whose features are listed in Table 2:

• AMD Istanbul: this system is equipped with four hexa-core AMD Istanbul processors clocked at 2.4 GHz. Each of these CPUs has six cores and is attached to a DRAM module through two DRAM controllers and the CPUs are connected to each other through HyperTransport links in a ring layout. If, say, core 0 on processor 1 needs some data which is located in the memory attached to processor 0, a request has to be submitted to the DRAM controllers of processor 0 and, when this request is processed, the corresponding data is sent through the HyperTransport link. Concurrent access to the same memory causes conflicts which slow-down the access to data.

•

Choosing the block size

As explained in Section 3.3, the execution time of the factorization operation depends on the nb and ib parameters which define the granularity of tasks (and thus the amount of concurrence) and the blocking of dense linear algebra operations, respectively. The optimal values for these two parameters depend mostly on the input matrix structure, the architectural features of the underlying computer and the number of working threads. Clearly, small ib values yield lower fill-in but inefficient BLAS operations and vice versa; on the other hand, small nb values deliver more parallelism but increase the complexity and the cost of the scheduling and set an upper bound to the ib parameter which may not be optimal.

Figure 10 shows how the execution time is affected by different choices of the two blocking parameters on the Power6 machine for two different matrices. It can be seen that, for both matrices, smaller nb values are required in order to achieve the shortest execution time on 32 cores compared to the 16 cores case. Although it is very expensive to determine experimentally the optimal values for nb and ib, Figure 10 as well as experiments conducted on the other test matrices and architecture confirm that any reasonable (according to common knowledge on BLAS libraries) choice for the two parameters falls very close to the optimum. Choosing nb to be a multiple of ib, generally delivers better results.

For all the 51 matrices in the test set and for the two architectures in Table 2, the blocking values were chosen among three different combinations, i.e., (nb,ib)=(120, 120), (120, 60) or (60, 60) (this last one is not shown in Figure 10) as those that delivered the shortest factorization time using all the cores available on the system. All the experimental results presented in the rest of this section are related to that choice.

Understanding the memory utilization

As described in Section 3.2, the scheduling of tasks in qrm is based on a method that aims at maximizing the locality of data in a NUMA environment. The purpose of this section is to provide an analysis of the effectiveness of this approach. This analysis was conducted on the AMD system, always using all of the 24 cores available, with the PAPI [START_REF] Mucci | PAPI: A portable interface to hardware performance counters[END_REF][START_REF] Terpstra | Collecting performance data with PAPI-C. Tools for High Performance Computing[END_REF] tool. Recalling the architectural characteristics described in Section 4, the efficiency of the scheduling technique has been evaluated based on three metrics:

• the completion time;

• the amount of data transferred on the HyperTransport links5 ;

• the number of conflicts on the DRAM controllers 6Experiments were run on all the matrices in the test set and with three different settings for the memory policy:

• no locality: this is a code variant where the multiple task queues are replaced with a single one shared by all threads. Tasks are, thus, pushed and popped from this queue regardless of their affinity with the data placement in the memory system;

• locality: this corresponds to the scheduling strategy described in Section 3.2, i.e., each task is pushed on the queue attached to the thread which owns the related front;

• round robin: this setting uses the same code variant of the "locality" one but, in this case, allocated memory pages are interleaved in a round robin fashion over all the DRAM modules. This is achieved using the numactl [START_REF] Kleen | An NUMA API for Linux[END_REF] tool with the "-i all" option. Since the code does not control the placement of data in the memory system, the ownership of fronts does not make sense anymore and, consequently, the data locality is completely destroyed. Table 3: Analysis of the memory usage for the matrix factorization on the AMD system with 24 threads.

Table 3 shows the results of these experiments for a subset of the matrices in Table 1. It can be seen that the locality aware scheduling strategy described in Section 3.2 provides better execution times with respect to a naive dynamic scheduling policy: the improvement may be as high as almost 20% (for the EternityII E matrix) and quite often around 10%. This can be explained by the reduced amount of data transferred over the HyperTransport links as shown in the middle of the table. However, the best execution times are achieved with the round robin distribution of memory allocations. In this case, the amount of data transfers is higher than the other two cases since the frontal matrices are completely scattered over the NUMA nodes; nonetheless, this more even distribution of data reduces the number of conflicts on the DRAM controllers (see Table 3 (bottom)) and provides a better use of the memory bandwidth. This behavior is coherent to what was observed on the HSL MA87 [START_REF] Hogg | A DAG-based sparse Cholesky solver for multicore architectures[END_REF] code which uses a similar approach for the Cholesky factorization of sparse linear systems.

The considerable performance improvement resulting from the interleaved memory allocation suggests that reducing the memory conflicts may be more important than minimizing the amount of data transferred over the Hyper-Transport links. A closer look to the behavior of matrices e18 and flower 7 4 (columns three and four in Table 3) shows, instead, that both these objectives are very important for the efficiency of the code. The factorization of these two matrices takes roughly the same time in the "locality" and "no locality" cases but the memory interleaving provides only a small improvement to the second matrix due to a bigger increase of the data traffic on the HyperTransport links.

Maximizing the data locality and minimizing the memory conflicts are not conflicting objectives although it may be rather complicated to achieve both of them at the same time on a very heterogeneous workload such as a sparse 4: The effect of tree pruning and tree reordering on the matrix factorization time on the AMD system with 24 threads.

factorization. It has to be noted that in the proposed locality aware scheduling policy, the placement of data and the ownership of the associated tasks is defined on a front basis. Because the number of fronts becomes smaller than the number of working threads when the factorization approaches the root front, a lot of work stealing and memory contention take place on the top part of the assembly tree where most of the work is done. Further improvements can be obtained by defining the data placement and the tasks affinity on a block-column basis; despite this would make some operations such as the front assembly much more complex, it is reasonable to expect that it will allow to more evenly and efficiently distribute the data and thus improve the locality of reference and reduce the memory contention at the same time. This is the object of ongoing work.

The effect of tree pruning and reordering

The tree reordering and pruning techniques have been evaluated on the test matrices. For the tree pruning, the initial threshold was set to 0.01 which means that all the subtrees whose weight is smaller than 1% of the total factorization workload are pruned off. If the remaining tree does not provide enough treelevel parallelism, the threshold is divided by 2 and a new pruning is done on the original assembly tree. More precisely this procedure is iterated until the number of leaves in the pruned tree is bigger than twice the number of working threads. Clearly, the optimal values for both the starting threshold and the stopping criterion depend on the structure of the tree and, therefore, on the specific input matrix; the values described above were determined experimentally and were found to work well on the large set of test matrices previously described. Table 4 shows the experimental results related to a subset of matrices for which these techniques have proved to be particularly effective. The experiments show that, when applied separately, these two methods provide considerable benefits in some specific cases. The tree reordering yields good improvements on very unbalanced and irregular trees: this is the case of the LargeRegFile and sls matrices. The tree pruning, proved to be very effective on all the problems but particularly on those with extremely large trees and extremely small frontal matrices such as the cont11 l and the sls matrices.

It can be observed that the pruning clearly reduces the need for sorting as well as its effectiveness. Nonetheless on some matrices (see, particularly, the first three columns in Table 4) the best execution time is achieved when both techniques are applied. It is important to note that the size of the pruned tree increases with the number of working threads and, therefore, the improvements provided by the sorting, when both techniques are applied, are likely to be more important for higher degrees of parallelism.

Absolute performance and Scaling

This section shows experimental results aiming at evaluating the efficiency and scaling of the proposed approach to the parallelization of the multifrontal QR method.

The qrm code was compared to the SuiteSparseQR [START_REF] Davis | Algorithm 915, SuiteSparseQR: Multifrontal multithreaded rank-revealing sparse QR factorization[END_REF] (referred to as spqr) released by Tim Davis in 2009; because the spqr is based on Intel TBB which is only available for x86 systems, this comparison is only made on the AMD Istanbul system. For both packages, the COLAMD matrix permutation was applied in the analysis phase to reduce the fill-in and equivalent nodes amalgamation methods were used so that the differences between the produced assembly trees can be considered negligible. Note that other ordering tools based on nested dissection (such as METIS [START_REF] Karypis | A fast and high quality multilevel scheme for partitioning irregular graphs[END_REF] or SCOTCH [START_REF] Pellegrini | Scotch: A Software Package for Static Mapping by Dual Recursive Bipartitioning of Process and Architecture Graphs[END_REF]) will produce more balanced and wider assembly trees which means better tree parallelism. Using such orderings will, thus, change the behavior of both packages as well as the outcome of the comparison; however, it is reasonable to expect that differences will not be substantial. Both packages are based on the same variant of the multifrontal method (that includes the two optimization techniques discussed in Section 2) and, thus, the number of floating point operations done in the factorization and the number of entries in the resulting factors are comparable. For both codes, runs were executed with numactl memory interleaving where available (i.e., only on the AMD Istanbul machine).

Figure 11 shows time and memory profiles [START_REF] Dolan | Benchmarking optimization software with performance profiles[END_REF] for both codes using 24 threads. For any given code, a data point (x, y) on the profile means that the code is no worse than x times the best of the two codes for y problems. This means that the (1, y) point gives the number y of problems on which the associated code was found to be the best. The time profile shows that qrm is faster than spqr on 43 problems out of 51; the remaining eight problems are very small, i.e., their factorization takes less than one second. On basically half the problems in the test set qrm was more than twice as fast as spqr. The memory profile is more difficult to interpret. Because qrm is based on a much more eager execution model it is reasonable to expect that it consumes more memory. Indeed, as explained in Section 3.2, the scheduling method exploits as much as possible the node parallelism and resorts on tree parallelism, by activating new nodes, only when needed; this keeps the tree traversal close to the one followed in the sequential case and limits the memory consumption growth in parallel. Moreover, the tree reordering technique, as explained in the same section, helps reducing the memory footprint. As a result, qrm achieves, in general, a smaller memory consumption than spqr on the 51 test matrices as shown in the memory Figure 11: Time and memory profiles comparing the qrm and spqr factorizations for the 51 test matrices on the AMD Istanbul system. profile of Figure 11. Experimental data show that most of the matrices where qrm has a better memory consumption are relatively small, which may probably be explained with some overhead in spqr that becomes negligible for bigger size problems. Note that, in both codes, the memory footprint is defined as the peak memory consumption reached during the factorization, including the input matrix and all the allocated memory areas of any type. On bigger problems, the two codes have a similar memory consumption and no clear winner can be identified. If the H matrix was kept in memory, which is not the case for the results reported in Figure 11, the difference between the two codes would be even smaller because the relative weight of the contribution blocks, responsible for the increased memory consumption in parallel, would be lower.

Tables 5 and6 show detailed performance results for a selected subset of problems. These where chosen to be the biggest ones not belonging to the same family in the complete test set and correspond to matrices 3-12 in Table 1. For each architecture, results were measured for number of threads equal to multiples of the number of cores per NUMA node, i.e., six and eight for the AMD system and IBM system, respectively.

The number of threads participating in the factorization in the spqr code is given by the product of the number of threads that exploit the tree parallelism times the number of threads in the BLAS routines. For each fixed number of threads, the best results among all the possible combinations are reported in Table 5. As discussed in Section 3, this rigid partitioning of threads may result in suboptimal performance; choosing a total number of threads that is higher than the number of cores available on the system may yield a better compromise. This obviously does not provide any benefit to qrm. The last line in Table 5 shows, for spqr, the factorization times for the best combination of tree 1, for both codes, on both systems for core numbers from one to the maximum available.

and node parallelism; for example, for the sls matrix (number 11) the shortest factorization time is achieved by allocating 22 threads to the tree parallelism and 7 to the BLAS parallelism for a total of 154 threads. Table 5 shows that the proposed approach clearly outperforms the classical approach to parallelization of the QR multifrontal method on the AMD system as qrm achieves speedups of more than three over spqr. Figure 12 shows the best, worst and average speedup achieved by the two codes on the AMD Istanbul system and by qrm only on the IBM Power6. The best cases are the TF17 and the Hirlam matrices on the AMD and IBM systems, respectively, whereas the worst are the sls and cont11 l ones, respectively; the best and worst case matrices are the same for qrm and spqr. Figure 13 shows the fraction of the peak performance achieved by qrm on the two systems in the best and worst cases which correspond to the TF17 and cont11 l (number 12 and 3 in Table 1) matrices, respectively. The absolute performance and scalability of qrm is compared to that of the matrix multiply operation (DGEMM) and the QR factorization (DGEQRF) of square dense matrices of size 40,000 which is big enough to achieve asymptotic performance using any number of cores on the target machines.

The experimental results reported above show that qrm achieves a good scalability up to 24 cores, the best speedup being more than 19 for the TF17 matrix on the AMD system. The scalability of the qrm code tends to level-off past 24 threads. A closer analysis of the execution profiles shows that this may be due to a lack of parallelism explained by the fact that panel operations, which lie on the DAG's critical path, are extremely inefficient. This problem becomes more evident when frontal matrices tend to be "tall and skinny" (i.e., having many more rows than columns) and thus the relative cost of panel operations is even higher; this is, for example, the case of the sls matrix on which qrm achieves a poor scalability and absolute performance, although still doing better than spqr. It has to be noted that DGEQRF performance in Figure 13 is related to the factorization of a square matrix; when run on strongly overdetermined matrices the DGEQRF routine achieves a very poor absolute performance and scaling as well.

Two techniques may be used to overcome this limitation and are the object of ongoing research work:

• 2D node parallelism: following the methods presented Buttari et al. [START_REF] Buttari | A class of parallel tiled linear algebra algorithms for multicore architectures[END_REF] and Hadri et al. [START_REF] Hadri | Tile QR factorization with parallel panel processing for multicore architectures[END_REF], a 2D partitioning could be applied to frontal matrices; this delivers finer granularity and thus better parallelism although it requires more logic to correctly handle the assembly operations and to keep track of the status of the global factorization.

• tree preprocessing: the method proposed by Hadri et al. [START_REF] Hadri | Tile QR factorization with parallel panel processing for multicore architectures[END_REF] for the QR factorization of tall and skinny dense matrices relies on the idea of parallelizing the computations according to a tree reduction pattern. This idea can be easily adapted to the multifrontal QR factorization of sparse matrices, where the concept of tree is already present.

For example, the assembly tree in Figure 14 (left) could be rearranged into the tree in Figure 14 (right) which provides more tree parallelism and makes the root front less overdetermined. It has to be noted that the two intermediate nodes that appear in the assembly tree on the right, do not receive any coefficient from the original matrix in the assembly operation and do not provide any coefficient to the global R factor. Because the staircase structure of frontal matrices is exploited, no additional floating point operations are introduced by this modification which is completely transparent to the factorization phase and only involves some symbolic operations during the problem analysis.

Figure 14: Fronts splitting in order to increase parallelism.

Summary

A novel approach to the parallelization of the multifrontal QR factorization of a sparse matrix has been presented. A fine grained partitioning is applied to the frontal matrices and computational tasks are defined as the sequential execution of an elementary operation on a data partition. The tasks are arranged into a DAG where the edges define the dependencies among them and, thus, the order in which they have to be executed. This DAG globally retains the shape of the classical assembly tree in the multifrontal method but exposes the concurrency in both the tree and node types of parallelism allowing for a consistent use of the two. Following a dataflow execution model, the tasks in the DAG can be scheduled dynamically depending on different policies. The scheduling of tasks is done according to a technique which aims at minimizing the transfer of data between NUMA nodes of large multicore systems and experimental results were presented which prove its effectiveness. Nonetheless, through a low-level performance profiling based on the PAPI tool, it has been shown that minimizing the conflicts on the memory system is more important that maximizing the data locality although both play an important role in the scalability of a parallel, multithreaded code; a possible way of achieving both these objectives based on a finer data placement was suggested (this is currently under investigation). Two techniques to reduce the size of the scheduling search space have been described. Finally, experimental results have been presented on two different architectures for 51 test matrices to prove the efficiency of the proposed approach. These results show that an actual implementation of the proposed methods achieves good scalability and good absolute performance up to 32 cores outperforming by a factor of up to three the best equivalent software currently available. The causes of poor efficiency on some of the test cases have been identified and possible solutions have been presented which are the object of ongoing research work.

Figure 3 (Figure 3 :

 33 Figure 3: Block-column partitioning of a frontal matrix (left) and panel and update operations pattern (middle and right, respectively); dark gray coefficients represent data read by an operation while black coefficients represent written data.

Figure 5 :

 5 Figure 5: Pseudo-code showing, on the left, the main execution loop and, on the right, the instructions used to fill the task queues.

Figure 9 :

 9 Figure9: Some details of the 51 matrices in the test set. On the left side, the operation count obtained with the COLAMD ordering and on the right side the qrm sequential factorization time (in seconds) on the AMD Istanbul system. In each graph, matrices are sorted in increasing order of the reported measure.

Figure 10 :

 10 Figure 10: The effect of the blocking sizes on the factorization time of the e18 and Rucci1 matrices on the Power6 computer.

Figure 13 :

 13 Figure 13: Fraction of the peak performance achieved by qrm compared to the QR factorization of a dense matrix.

Table 1 :

 1 A subset of the test matrices. The operation count is related to the matrix factorization with COLAMD column permutation.

		Mat. name	m	n	nz op. count
						(Gflops)
	1 image interp	240000	120000	711683	30.2
	2 LargeRegFile	2111154	801374	4944201	84.7
	3	cont11 l	1468599	1961395	5382999	184.5
	4	EternityII E	11077	262144	1572792	544.0
	5	degme	185,501	659415	8127528	591.9
	6	cat ears 4 4	19020	44448	132888	716.1
	7	Hirlam	1385270	452200	2713200	2339.9
	8	e18	24617	38602	156466	3399.5
	9	flower 7 4	27693	67,593	202218	4261.2
	10	Rucci1	1977885	109900	7791168	12768.5
	11	sls	1748122	62729	6804304	22716.2
	12	TF17	38132	48630	586218	38203.1

Table 2 :

 2 IBM Power6 p575: one node of the larger Vargas supercomputer installed at the IDRIS supercomputing institute is equipped with 16 dualcore Power6 processors clocked at 4.7 GHz. Processors are grouped in sets of four called MCMs (Multi Chip Module) and each node has four MCMs for a total of 32 cores. Processors in an MCM are fully connected as well as MCMs in a node. Test architectures.

	All the tests were run with real data in double precision.

Table 5 :

 5 Factorization times, in seconds, on the AMD Istanbul system for qrm (top) and spqr (bottom). The first row shows the matrix number.

						IBM Power6
	Matrix	3	4	5	6	7	8	9	10	11	12
		th.
		1	36.2 51.9	60.5	75.2 227.3 290.2 426.1 1156 2142 3121
		2	24.5 27.6	32.0	38.6 117.2 151.7 192.9	600 1259 1656
	qrm	4 8	10.6 14.0 5.9 7.3	16.0 8.2	19.4 9.9	57.1 28.6	76.5 39.3	98.0 49.5	304 155	672 348	850 458
		16	3.9	4.1	4.7	5.7	15.4	21.3	27.6	84	197	235
		24	4.5	3.1	3.7	4.3	11.3	15.5	22.8	62	183	174
		32	6.4	2.6	3.3	4.0	9.5	12.9	18.9	52	182	141

Table 6 :

 6 Factorization times, in seconds, on the Power6 system for qrm. The first row shows the matrix number. Figure 12: Best, worst and average speedup achieved on matrices 3-12 in Table

The rows of A are sorted in order of increasing index of the leftmost nonzero in order to show more clearly the computational pattern of the method on the input data.

Figure 4 actually shows the transitive reduction of the DAG, i.e., the direct dependency between two nodes is not shown in the case where it can be represented implicitly by a path of length greater than one connecting them.

http://buttari.perso.enseeiht.fr/qr_mumps

http://hirlam.org

This quantity was measured as the number of occurrences of the PAPI HYPERTRANSPORT LINKx:DATA DWORD SENT event (where x is 0, 1, 2, 3 since each processor has 4 HyperTransport links) which counts the number of double-words transferred over the Hy-perTransport links.

[START_REF] Buttari | A class of parallel tiled linear algebra algorithms for multicore architectures[END_REF] This quantity was measured as the number of occurrences of the PAPI DRAM ACCESSES PAGE:DCTx PAGE CONFLICT event (where x is 0, 1 since each processor has two DRAM controllers) which counts the number of conflicts occurring on the DRAM controllers.

Acknowledgments

I wish to express my gratitude to the members of the MUMPS team for the countless advices they gave me on the implementation of multifrontal methods and to Chiara Puglisi for sharing with me her deep knowledge of the sparse QR factorization techniques.

I also wish to thank Dan Terpstra for helping me with using the PAPI tool and interpreting the results it gave.

This work was performed using HPC resources from GENCI-IDRIS (Grant x2012065063).

AMD Istanbul


Maximilien Naveau

Justin Carpentier

Sébastien Barthelemy

Olivier Stasse

Philippe Souères

P Soueres

METAPOD Template META-programming applied to dynamics: CoP-CoM trajectories filtering

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

I. INTRODUCTION

Computing dynamic state estimate for a complex redundant robot such as humanoid robot is a crucial technique for controlling it. For instance, Inverse Dynamics control heavily relies on the computation of the inertia matrix [START_REF] Mansard | A dedicated solver for fast operational-space inverse dynamics[END_REF], [START_REF] Prete | Prioritized optimal control[END_REF]. When considering balance, this matrix is of primary interest. Indeed most of the real-time approaches relies on the Linear Inverted Pendulum Model to generate a Center-Of-Mass reference trajectory which complies with the balance criteria provided by the Center-Of-Pressure and the contacts on the ground. Inverse Dynamics is necessary to counteract the inertial effect involved by the limbs: either by modifying the center-of-mass (CoM) reference [START_REF] Nishiwaki | Online walking control system for humanoids with short cycle pattern generation[END_REF], either by following the momentum reference trajectories. More advanced approaches are using Model Predictive Control (MPC) [START_REF] Feng | 3d walking based on online optimization[END_REF], [START_REF] Tassa | Control-limited differential dynamic programming[END_REF] and need to compute the whole robot dynamic over a time horizon.

Since 2006, processor speed is limited to 3∼4 GHz, because transistors have reached their physical limits mostly due to heat dissipation [START_REF] Zhirnoz | Limits to binary logic switch scalinga gedanken model[END_REF]. The Moore Law is still valid because as the chipsets integrate more cores, the transistor number continues to grow not because of their size, but because of their density. They are overheating when CPU clock is increased, and the focus is now more on concurrent software programming to take advantage of multiple cores. In the case of MPC, it mainly consists in having multiple instances starting with different initial solutions. Boosting *This work was supported by the European project KOROIBOT FP7-ICT-2013-10/611909, the French National project OSEO ROMEO-2 and by AIRBUS/Future of the Aircraft Factory.

1 are with the CNRS, LAAS,7 av. du Colonel Roche, F-31400, Toulouse, France, Univ de Toulouse, LAAS, F-31400, Toulouse, France mnaveau,jcarpent,ostasse,soueres@laas.fr 2

Sébastien

Barthelemy is with Aldebaran Robotics sbarthelemy@aldebaran.com the code running concurrently will also increase the performances. It is now necessary to implement code to compute dynamics efficiently while taking into account the underlying computing architecture.

A. Problem formulation

The rigid multibody dynamic model of the robot is assumed to be known and called model. According to [START_REF] Featherstone | Rigid Body Dynamics Algorithms[END_REF], the dynamical equations can be written under the canonical form:

H(model , q) q + C(model , q, q) = τ

where q, q et q are respectively the generalized vectors for positions, velocities and accelerations. H is the generalized inertia matrix, C is the matrix representing the centrifugal, gravitational and Coriolis forces. We should then compute as fast as we can:

• τ by using the Recursive Newton-Euler Algorithm (RNEA) with the following inputs q, q et q. • H(model , q) by using the Composite-Rigid-Body Algorithm (CRBA). To represent position, velocities, acceleration, and forces, spatial algebra is nowadays one of the most used representation in robotics. Introduced by Roy Featherstone [START_REF] Featherstone | Rigid Body Dynamics Algorithms[END_REF], it allows a correct mathematically algebra notation to represent the relationships between quantities of a rigid multibody system.

B. Applications

1) Walking: In the context of humanoid robotics, and following the method initiated by Nishiwaki et al. [START_REF] Nishiwaki | Online walking control system for humanoids with short cycle pattern generation[END_REF] for reactive walking, RNEA can be used to filter dynamically a Center-of-Mass trajectory. Typically, for a humanoid robot, the goal is to compute 160 iterations for 30 DOFs robot in less than 200 µs. This corresponds to a 1.6 s window for an optimal controller as the one used classically for the HRP-2 robot, with a 200 Hz sampling.

To reach this objective, a RNEA call should last no more than 1.25 µs. The result achieved by our library is 11.32 µs on a i7-2820QM CPU of 2.3 GHz. This imply to use a subsampling mechanism such as the one proposed by Nishiwaki et al. [START_REF] Nishiwaki | Online walking control system for humanoids with short cycle pattern generation[END_REF].

One may argue that sub-sampling this trajectory can decrease the overall complexity for walking (Nishiwaki uses a 20 ms time step). However we anticipate a demand on such information with the recent increase in work considering extreme motion involving strong angular moments (such as the Parkour [START_REF] Dellin | A framework for extreme locomotion planning[END_REF]). In such situation linearization will be only correct at high frequency.

2) Whole-body control: Staying in the context of humanoid robotics, and by now considering a control architecture based on the inverse dynamics [START_REF] Sentis | A whole-body control framework for humanoids operating in human environments[END_REF], [START_REF] Mansard | A dedicated solver for fast operational-space inverse dynamics[END_REF], computing the inertia matrix at a very high speed may be necessary. Following the previous remark on sub-sampling, this is probably not required for motion with a slow dynamic.

C. State of the art

A classical approach is to implement RNEA and CRBA using instances of classes and structures, and to iterate over each element of the kinematic tree. RBDL (Rigid Body Dynamic Library) and KDL (Kinematic Dynamic Library) are using this approach. RBDL is following specifically the directives suggested by Featherstone in his book [START_REF] Featherstone | Rigid Body Dynamics Algorithms[END_REF] by implementing the spatial algebra through C++ operators overloading. A different approach is to use symbolic computation. One of the first software implementing this approach is Symoro+ by Khalil et al. [START_REF] Khalil | Symoro+: A system for the symbolic modelling of robots[END_REF]. The software is able to produce reduced models and to generate code where unnecessary computation are stripped out. SD/FAST is a commercially available toolbox used by Boston Dynamics and the MIT leg laboratory [START_REF] Inc | Sd/fast[END_REF]. HuMAnS is a toolbox written by P.-B. Wieber [START_REF] Wieber | The humans toolbox, a homogeneous framework for motion capture, analysis and simulation[END_REF] relying on Maple for manipulating the robot kinematic tree. ROBOTRAN is a JAVA application which performs symbolic manipulation and generates code integrated in MATLAB and the Simulink toolbox [START_REF]ROBOTRAN: Symbolic generator of multibody systems[END_REF]. In each case, the robot model is written in a specific language to be simplified and finally generated in C++.

In Chapter 10 of his book, Featherstone indicates that a potential disadvantage of symbolic simplification is that it does not use vector-based arithmetic of CPUs, and that if the code is too large then it does not fit the computer's instruction cache. In this work the first point is addressed by using Eigen, a widely known linear algebra library carefully written to exploit vector-based floating-point units. The second point is solved by using template programming which allows each instanciation of the template to invoke common code while avoiding indirect calls.

More generally, in this work, we propose to use metaprogramming to exploit directly the C++ capabilities to perform symbolic computation at compile time without any additional library or software other than the BOOST libraries. The current solution is completed with a C++ program which is parsing URDF files and generates appropriate structures to decouple model and instance of robots. The current limitation of the approach is the case where model modification is needed. A solution (not yet implemented in the current software solution) is an abstract layer switching from an optimized implementation to a generic implementation (RBDL or KDL). Another more complex solution is to compile the new model and link dynamically to this new library.

Our contribution is to have written, to our knowledge, one of the first meta-programming library aiming at computing inverse dynamic for the robots. The application targeted in this paper is humanoid robot walking. By putting together the power of functional programming with C++ speed, our library is able to achieve performance necessary for real-time control.

II. TEMPLATE METAPROGRAMMING FOR INVERSE DYNAMICS

This work takes its root from the fact that code generated by symbolic formal computation has in general better performances than generic code for complex robotic models. This comes for the inner structure of the problem itself. For instance, the robot inertia matrix is sparse because the subblocks matrices correspond to the kinematic branches. Then when no coupling exists between some branches, the related sub-blocks mainly is equal to zero. Formal computation detects such branches and suppresses the related computation. We are proposing to make the same kind of operations by exploiting C++ template meta-programming.

The result is a template library able to develop an optimized computational tree to evaluate the inverse dynamics of a robot class named Robot and to apply it onto the specific instance named robot, while considering the following complete state (q, q, q). During execution, the model specificities are used to update robot internal variables according to (q, q, q).

A. Simple example of functional programming using C++

A famous example of functional programming in C++ is the factorial computation: At compile time, when the compiler finds factorial<4> it directly computes 24 by successive substitutions. The same type of mechanism is used to develop the RNEA algorithm from the model Robot. In order to achieve a decoupling between the model and the robot instances Metapod is using a template based container proposed by boost::fusion which allows an array like syntax with templates 1 . The kinematic structure is then embedded through indexes that are found classically in an array, but here declared through boost::fusion.

B. Spatial algebra and meta-programming

The spatial algebra can be easily implemented in metaprogramming thanks to static polymorphism. This is allowing to keep a very generic code quite similar to Oriented Object programming but without the additional cost induced by virtual methods. We can also avoid errors by using strongly typed operators. It has also been possible to implement the optimization described in Appendix 2 of Featherstone's book [START_REF] Featherstone | Rigid Body Dynamics Algorithms[END_REF]. This is especially useful to simplify computation when parts of the robot are planar. One of the unitary test used in the library available through Github is the planar arm described in [START_REF] Spong | Robot Modeling and Control[END_REF].

C. Visitor design pattern to compute inverse dynamic

The depth first exploration of the robot structure is performed thanks to the visitor design pattern. This design pattern has been used in all algorithms where an iterative exploration of the tree structure is necessary. Finally, thanks to the design pattern coupled to boost::fusion which make possible to aggregate different types, it is possible to decouple the model from its instance.

D. Benchmarks

We compared the performances of Metapod with RBDL2 by computing the inertia matrix and the inverse dynamics on the CPUs listed in Table I. First we would like to emphasize that the numbers are only a snapshot and are subject to modifications, since the writing of this paper numerous improvements have been found. Some are currently tested for RBDL and for Metapod lead to the writing of a new library. The reader is invited to look at the link given in the footnote for further discussion. In each case, the library was recompiled using the GNU gcc compiler with the optimization options (-O3 -NDEBUG) trying to generate the code producing the best performances. The model used for benchmarking is the sample humanoid model provided by Tokyo University in OpenHRP. For this specific work we concentrated in two algorithms which are used in controlling a humanoid robot, the Composite Rigid Body Algorithm (CRBA) which computes the inertia matrix, and the Recursive Newton-Euler Algorithm (RNEA) which provides the torques necessary to counterbalance the gravity.

Table II provides a quick view on the performance difference between Metapod and RBDL on various CPUs. The next section will detail some aspects related to the CPUs which will try to explain that the values provided here are to be taken cautiously and are likely to highly vary according to the CPU, the context of its overload, the compiler and the operating system. For this reason, we provide in Table II a ratio according to the CPU between the time provided by Metapod and RBDL. Note that CPU-2 corresponds to the HRP-2 real chipset with its real-time operating system. The interested user should perform an identical comparison for its own robot to have accurate time measurement. Here we have put the most favorable time measurements for both algorithms. Despite this high variability, Metapod is almost always better in performances than RBDL.

We have started some comparisons with symbolic computation libraries, but they usually imply the use of a costly third-part software. A comparison with HuManS on a structure similar to the model used here gives 8.35µs for the CRBA algorithm on CPU-1. To achieve this result on Metapod with CPU-1 it was necessary to use the compiler options given in the remainder of the paper. They changed nothing for HuManS. As the later one does not use Eigen, this implies that, as mentionned by Featherstone, the use of vector-based arithmetic has a strong impact when computing the inertia matrix. Thanks to the authors of Symoro+, and if the paper is accepted, we hope to augment the comparison with this software.

E. CPUs

The Intel company divides its chipsets into 3 categories: Core, Xeon and Atom. With its low power consumption and low heat dissipation the Atom processor is aimed for embedded systems. The first and second families are used respectively for desktop applications and high power computation. The Xeon Phi category is targeted for many cores applications.

1) TurboBoost: To speed up the performance of the Core and Xeon families their coprocessors embed a technology that Intel names TurboBoost. It consists in increasing the CPU speed while respecting safety rules regarding heat dissipation and voltage. The maximum speed that can be reached is called the maximum turbo frequency. For instance regarding Table II CPU-1 and CPU-3 can reach respectively 3.4 GHz (from 2.3 GHz) and 3.8 Ghz (from 3.6 GHz). There is no guarantee that the maximum frequency is always reached as it depends on the current status on the chipsets. The Atom chipsets do not have the TurboBoost technology as well as the Core2 Duo E7500 (CPU-2) currently inside the HRP-2 humanoid robot. This is reflected by the numbers given in Table II, where the raw times are 10 times faster on Xeon chipsets than on Atom. Interestingly one would have expected that CPU-1 with the turbo-boost would be at best near 0.9 the speed of CPU-3, as it is the ratio of their respective Turbo-Boost frequency. This relationship is not verified between CPU-1 and CPU-3. In addition the description of the TurboBoost technology in the Sandy Bridge micro architecture [START_REF] Rotem | Power-management architecture of the intel microarchitecture code-named sandy bridge[END_REF] shows a complex pattern depending on the current overall status of the CPU. It is worth noting that the operating system can also request TurboBoost by sending a signal (P-state request) if it detects an overloading status.

2) Branch prediction: Performances depend on another key technology: branch prediction. The chipset maintains counters to predict the branches behavior. Depending on the statistics of the branch, the predictor tries to guess what will we be the most probable code to be executed by the CPU. The code and its data are then prefetched from the memory while the CPU is decoding the current operation. Branch predictors can also detect cycles and save in cache intermediate variables to avoid costly access memory. Branch predictors is also able to indirect branch where multiple address are possible targets. The interest of the symbolic approach is to avoid as much as possible to make branches in the code. It facilitates the work of the branch predictor as well as the memory pre-fetching mechanism. This may explain why in Table II INDIRECT BRANCHES MISPREDICTED. TESTS REALIZED ON CPU-1 scription of their branch predictors. For this reason, users have usually to reverse engineer the behavior of the chipset in order to create branch predictor aware compilers [START_REF] Wang | Program interferometry[END_REF], [START_REF] Uzelac | Experiment flows and microbenchmarks for reverse engineering of branch predictor structures[END_REF]. We used Valgrind framework to analyse cache access and branch prediction. The result is compiled in Table III and in Table IV. From Table III Metapod is doing more misprediction than RBDL for the Level 1 cache both for instruction and data, but less for the last level cache. As a last level misprediction implies a memory access, with the slower bus frequency it is the most costly. Therefore this result is in favor of Metapod. In addition Table IV shows clearly that branch misprediction is far lower for Metapod than for RBDL.

3) Compiler: The importance of the compiling options differs according to the library and the CPUs. It was particularly strong on the Atom chipset with a gain of 20µs for Metapod and RBDL. On RBDL, the impact of compiling option was null on CPU-1 and CPU-3 whereas it is improved by 10 % on Metapod. The compiler options used to exploit vector-based arithmetic and loop optimizations are3 : -msse -msse2 -msse3 -march=corei7 -mfpmath=sse -fivopts -ftreeloop-im -fipa-pta. In addition during the linking phase it is possible to strip out useless code and perform wholeprogram optimization with the following option: -fwholeprogram. Again the static structure of the code generated through Metapod simplifies the work of optimization strategies applied by the compiler (here gcc).

III. DYNAMIC FILTERING FOR WALKING

To illustrate the use of Metapod on a humanoid robot, we consider an application where the robot has to evolve in a factory. AIRBUS/Future of the Aircraft factory is currently evaluating the potential of humanoid robots in such context. The use case considered here is the HRP-2 humanoid robot bringing an electric screw driver to an assembly line. On the other hand, in the context of the Koroibot project aiming at studying human walking to improve humanoid robot motion, we have created a stair climbing motion primitive. The behavior consisting in bringing a tool while climbing the stairs is then a combination of the two motion primitives. We show next that, thanks to Metapod, the dynamical consistency can be realized in real time.

A. Dynamic filtering on a sub-sampled walking pattern generator

This work is based upon the real time walking control system described by Nishiwaki in [START_REF] Nishiwaki | Online walking control system for humanoids with short cycle pattern generation[END_REF]. One key ingredient is to generate real-time dynamically stable CoM and Center-of-Pressure (CoP) trajectories using the cart-table model [START_REF] Kajita | Biped walking pattern generation by using preview control of zero-moment point[END_REF]. Then to correct the inertial effects induced by the legs movement Nishiwaki uses a subsampled dynamical filter. Indeed, the results tested on HRP-2 have shown that the CoP trajectory followed by the robot are not the desired one, the difference can be as big as 15 mm. To avoid this problem, the solutions consist in designing a controller regulating the linear and angular momentum [START_REF]Resolved momentum control: Humanoid motion planning based on the linear and angular momentum[END_REF] in addition to a dynamical filter aiming at correcting the CoM trajectory in order to take into account the dynamics of the robot. The first kind of solution is mostly used to reject instantaneous perturbations. However, it uses additional degrees of freedom that could have been used for manipulation tasks for example. For this reason, in our specific case, we have only implemented the dynamical filter as the commercially available stabilizer integrate already a perturbation rejection strategy. In the following the dynamic filter depicted by Nishiwaki and implemented in this paper is described in details:

• The walking pattern generator (WPG) defines a trajectory for the CoM, the feet and the CoP. The desired CoP trajectory is noted CoP * . • The dynamic filter starts by computing the joint trajectories using the analytical inverse kinematics. We make the assumption that the CoM and the free flyer are rigidly connected.

(q, q, q) = IK(ponctual model, c, ċ, c, X f ) with c, ċ, c et X f being respectively the position, the velocity, the acceleration of the CoM and the feet position.

• With the joint trajectory, it is possible to compute the inverse dynamics and to find the CoP matching up to the real robot motion. We call it the multi-body CoP noted CoP M B .

(f, τ ) = ID(model complet, q, q, q)

CoP M B x = - τy mg CoP M B y = τx mg CoP M B z = 0 ∆CoP = CoP * -CoP M B
• This provides an error between CoP * and CoP M B . This error is computed over a time window and is injected in a preview control (PC) in the shape of an LQR described by Kajita and al. [START_REF] Kajita | Biped walking pattern generation by using preview control of zero-moment point[END_REF]. The result of this step is an error for the CoM in position, velocity and acceleration. • Finally, the new CoM trajectory is used to compute the joint trajectory using the Inverse Kinematics. This dynamic filter is based on the Newton-Raphson algorithm using the CoM trajectory as free variables. This method does not guarantee the convergence. However during tests on HRP-2 with other WPG, we observed that the CoP trajectory is corrected in one iteration. The goal is to implement this method on the WPG of Andrei Herdt [START_REF] Herdt | Online walking motion generation with automatic footstep placement[END_REF] and Morisawa's [START_REF] Morisawa | Experimentation of Humanoid Walking Allowing Immediate Modification of Foot Place Based on Analytical Solution[END_REF].

B. Application on HRP-2

We propose to use this filter on the WPG proposed by Morisawa and al. in [21]. This algorithm uses an analytical form of the CoM and CoP trajectories assuming that the later is a third order polynomial. The Morisawa's WPG can either define off line the whole trajectory at once, or change online the foot steps given as reference.

To use the filter with its full potential on the online Morisawa's WPG we need to compute the trajectory over 1.6 s in advance, i.e. two steps forward. The computation of the trajectory is done every 5 ms including the control, so the remaining time for the WPG is 3 ms. Moreover, because the robot HRP-2 is controlled at 200 Hz, we need to compute 320 RNEA to use the whole preview. Our implementation, as fast as it is, does not allow this kind of computation, hence we used a strategy proposed by Nishiwaki, i.e. we sub sampled the preview window at a period of 0.05 s which leads to 32 RNEA every 0.1 s. We rather compute the whole dynamic filter for the offline WPG of Morisawa. This study is based on the offline WPG. We computed 1.6 s more than the initial trajectory and we computed the multi-body CoP for each sampling time.

The results obtained are illustrated in the Fig. 1. The graphs here contains each :

• the trajectories of the reference CoP,

• the multi-body CoP computed with RNEA,

• and the corrected CoP obtained with the RNEA using the corrected CoM The upper graph depicts the evolution on the axis X in function of the time and the lower graph shows the evolution of the Y axis in function of the time. This data are extracted from a straight walking using the offline WPG of Morisawa.

We can here easily see the influence of the filter, i.e. the corrected multi-body CoP is almost fused with the reference one. In terms of distance, the maximum error between the reference CoP and the multi-body CoP is 15.3 mm while the maximum distance after correction is 1.7 mm. In average the error before correction is 5.7 mm and after it is 0.7 mm.

In term of computation time, we executed the algorithm that creates a straight flat walking on the HRP-2 robot. The results are 15.0 us for the average time used to compute the RNEA during the whole trajectory and 342.0 us for the maximum duration of one iteration of RNEA. This peak corresponds to the first iterations of the algorithm, corresponding to memory allocation.

We have done a series of experiments using the dynamic filter including the straight walking and climbing stairs. For the climbing we had four scenarios: one simple walk, a walk while swinging the arm and holding a lightweight tool in one hand, another one holding the tool with 2 hands. In the last one the tool is hold with the hands lifted over the head. HRP-2 went through all the pattern has depicted in the companion video. We also succeded to make the robot going down the stairs with the lightweight in two hands, however we did not manage to make HRP-2 perform the rest of the scenario while going down. One of the raisons that explain this behaviour is that the dynamic filter places the CoM in a way such that the CoP is inside the support foot. To compensate for the arms being in front, the CoM is send slightly backward. The leg is then stretched and create a singularity. This shows the limit of this approach which calls for a more advanced development such as optimal control, whole-body predictive control or planning.

In this section we have shown that one iteration of a Newton-Raphson algorithm can make the CoP converge towards its reference and it could be applied to WPG focusing on the under-actuated part of the robot.

IV. CONCLUSION

We have presented a C++ library for efficient humanoid dynamic computation without relying on any particular symbolic computational tool. The achieved performances shows some advantage over the state-of-the art dynamic library RBDL mostly on ATOM processor and for the inertia matrix computation, which are relevant for robotics application. On recent desktop computer, the ratio of the gain is not so obvious and in general the time achieved by both library is not significantly different for inverse dynamics. Some preliminary results show that we have a computational speed similar to symbolic computation. This is achieved through a mix of code generation and C++ meta-programming. Its use is established in the context of dynamic filtering for real-time walking gait generation, and more specifically applied to the humanoid robot HRP-2. The library is available under GPLlicense at the following url: https://github.com/ laas/metapod.

  template <int n> struct factorial { enum { value=n * factorial<n-1>::value};}; template <> struct factorial<0> { enum { value = 1 }; };

Fig. 1 .

 1 Fig. 1. Results obtained on the Morisawa offline WPG

TABLE I CPUS

 I Metapod provides a 25-30 % increase in efficiency on the CRBA algorithm. Unfortunately chipsets providers, in general, do not give a detailed de-USED FOR THE BENCHMARKING. CPU-1 IS A LAPTOP, CPU-2 IS ONE USED IN HRP-2, CPU-3 IS A DESKTOP, AND CPU-4 IS A NETPC. CPU-2 IS USING ONLY ONE CORE DUE TO THE REAL-TIME OPERATING SYSTEM RECOGNIZING ONLY ONE CORE.

	Nickname		CPU name			Frequency	Core number	Cache size	Distribution
	CPU-1	Intel(R) Core(TM) i7-2820QM CPU	2.3 GHz	8	8 Mb	Ubuntu 12.04 LTS
	CPU-2	Intel(R) Core2(TM) Duo E7500		2.8 GHz	1	3 Mb	Ubuntu 10.04 LTS
	CPU-3	Intel(R) Xeon(R) CPU E3-1240 v3		3.4 GHz	8	8 Mb	Ubuntu 12.04 LTS
	CPU-4	Intel(R) Atom(TM) CPU N720		1.6 Ghz	1	512 Kb	Ubuntu 12.04 LTS
	Library	Algorithms	CPU-1			CPU-2		CPU-3	CPU-4
			µs	ratio	µs		ratio	µs	ratio	µs	ratio
		CRBA	7.89	0.77	12.97	0.69	5.08	0.76	61.06	0.76
	Metapod	RNEA	9.28	0.97	13.85	0.75	4.93	0.82	60.98	0.90
		CRBA	10.16		18.59		6.60		79.68
	RBDL	RNEA	9.5		18.66		5.99		67.68

TABLE II PROCESSING

 II TIME FOR THREE BASIC ALGORITHMS USED IN CONTROL OF HUMANOID ROBOT: INERTIA MATRIX COMPUTATION (CRBA), INVERSE DYNAMICS (RNEA). THE COMPUTATION HAVE BEEN TESTED ON 4 CPUS.

	Library	Ir	I1mr	ILmr	Dr	D1mr	DLmr	Dw	D1mw	DLmw
	RBDL	8,295,618,744	39,643	3,995	3,447,245,971	2,936,723	7,322	1,617,919,310	607,848	3,447
	Metapod	4,529,574,180	69,802,287	2,722	2,154,169,130	9,709,886	5,040	1,034,216,699	27,501,390	1,136

TABLE III

 III Ir: NUMBER OF EXECUTED INSTRUCTIONS, I1m1r: INSTRUCTION READ MISSES FOR CACHE L1, ILmr: INSTRUCTION READ MISSES FOR LAST LEVEL CACHE, Dr: NUMBER OF CACHE READ, D1mr: DATA READ MISSES FOR CACHE L1, DLmr: DATA READ MISSES FOR LAST LEVEL CACHE, Dw: NUMBER OF CACHE WRITE, D1mw: DATA WRITE MISSES FOR CACHE L1, DLmw: DATA WRITE MISSES FOR LAST LEVEL CACHE. TESTS REALIZED Bi: INDIRECT BRANCHES EXECUTED Bim:

					ON CPU-1
	Library	Bc	Bcm	Bi	Bim
	RBDL	245,662,040	7,604,170	45,177,857	269,654
	Metapod	98,901,411	977,079	14,105,063	726
			TABLE IV
	Bc: CONDITIONAL BRANCHES EXECUTED, Bcm: CONDITIONAL
	BRANCHES MISPREDICTED,

A naive description of the approach is available at https://github.com/ laas/metapod/wiki/Naive-implementation---First-approach

The benchmark programs for rbdl evaluation are available at https:// github.com/laas/metapod/wiki/Tests-realized-in-comparison-with-RBDL.

We kindly invite the interested reader to follow the gcc manual for a more detailed explanation of the options

ACKNOWLEDGMENT

The Metapod library was made possible thanks to the work of Maxime REIS. Antonio EL-KHOURY took an active part in the improvement of Metapod. We warmly thank them.