Anthony Scemama

Michel Caffarel

Emmanuel Oseret

William Jalby

Quantum Monte Carlo for large chemical systems: Implementing efficient strategies for petascale platforms and beyond

Various strategies to implement efficiently QMC simulations for large chemical systems are presented. These include: i.) the introduction of an efficient algorithm to calculate the computationally expensive Slater matrices. This novel scheme is based on the use of the highly localized character of atomic Gaussian basis functions (not the molecular orbitals as usually done), ii.) the possibility of keeping the memory footprint minimal, iii.) the important enhancement of single-core performance when efficient optimization tools are employed, and iv.) the definition of a universal, dynamic, fault-tolerant, and load-balanced framework adapted to all kinds of computational platforms (massively parallel machines, clusters, or distributed grids). These strategies have been implemented in the QMC=Chem code developed at Toulouse and illustrated with numerical applications on small peptides of increasing sizes (158, 434, 1056 and 1731 electrons). Using 10k-80k computing cores of the Curie machine (GENCI-TGCC-CEA, France) QMC=Chem has been shown to be capable of running at the petascale level, thus demonstrating that for this machine a large part of the peak performance can be achieved. Implementation of large-scale QMC simulations for future exascale platforms with a comparable level of efficiency is expected to be feasible.

I. INTRODUCTION

Quantum Monte Carlo (QMC) is a generic name for a large class of stochastic approaches solving the Schrödinger equation by using random walks. In the last forty years they have been extensively used in several fields of physics including nuclear physics, [1] condensedmatter physics, [2] spin systems, [START_REF] Suzuki | Quantum Monte Carlo Methods in Condensed Matter Physics[END_REF], quantum liquids, [START_REF] Ceperley | [END_REF] infrared spectroscopy, [5,6] etc. In these domains QMC methods are usually considered as routine methods and even in most cases as state-of-the-art approaches. In sharp contrast, this is not yet the case for the electronic structure problem of quantum chemistry where QMC [START_REF] Hammond | Monte Carlo Methods in Ab Initio Quantum Chemistry[END_REF][START_REF] Caffarel | Encyclopedia of Applied and Computational Mathematics[END_REF] is still of confidential use when compared to the two well-established methods of the domain (Density Functional Theory (DFT) and post-Hartree-Fock methods). Without entering into the details of the forces and weaknesses of each approach, a major limiting aspect of QMC hindering its diffusion is the high computational cost of the simulations for realistic systems.

However -and this is the major concern of this work -a unique and fundamental property of QMC methods is their remarkable adaptation to High Performance Computing (HPC) and, particularly, to massively parallel computations. In short, the algorithms are simple and repetitive, central memory requirements may be kept limited whatever the system size, and I/O flows are negligible. As most Monte Carlo algorithms, the computational effort is almost exclusively concentrated on pure CPU ("number crunching method") and the execution time is directly proportional to the number of Monte Carlo steps performed. In addition, and this is a central point for massive parallelism, calculations of averages can be decomposed at will: n Monte Carlo steps over a sin-gle processor being equivalent to n/p Monte Carlo steps over p processors with no communication between the processors (apart from the initial/final data transfers). Once the QMC algorithm is suitably implemented the maximum gain of parallelism (ideal scalability) should be expected.

A most important point is that mainstream high-level quantum chemistry methods do not enjoy such a remarkable property. Essentially based on iterative schemes defined with the framework of linear algebra they involve the manipulation and storage of extremely large matrices and their adaptation to extreme parallelism and lowmemory-footprint implementation is intrinsically problematic. Now, in view of the formidable development of computational platforms, particularly in terms of number of computing cores (presently up to a few hundreds of thousands and many more to come) the practical bottleneck associated with the high computational cost of QMC is expected to become much less critical and thus QMC may become in the coming years a method of practical use for treating chemical problems out of the reach of presentday approaches. Following this line of thought a number of QMC groups are presently working on implementing strategies allowing their QMC codes to run efficiently on very large-scale parallel computers. [START_REF] Esler | [END_REF][START_REF] Kim | Proceedings of the 2010 Scientific Discovery through Ad-vanced Computing (SciDac) Conference[END_REF][START_REF] Gillan | Psi-k Highlight of the Month[END_REF] Essentially, most strategies rely on massive parallelism and on some efficient treatment ("linear-scaling"-type algorithms) for dealing with the matrix computations and manipulations that represent the most CPU-expensive part of the algorithm.

Here, we present several strategies implemented in the QMC=Chem code developed in our group at the University of Toulouse. [START_REF]Quantum Monte Carlo for Chemistry@Toulouse[END_REF] A number of actual simulations realized on the Curie machine at the French GENCI-TGCC-CEA computing center with almost ideal parallel efficiency in the range 10 000-80 000 cores and reaching the petascale level have been realized.

The contents of this paper is as follows. In The first section, a brief account of the QMC method employed is presented. Only those aspects essential to the understanding of the computational aspects discussed in this work are given. In second section, the problem of computing efficiently the Slater matrices at the heart of the QMC algorithm (computational hot spot) is addressed. A novel scheme taking advantage of the highly-localized character of the atomic Gaussian basis functions (not the molecular orbitals as usually done) is proposed. A crucial point is that the approach is valid for an arbitrary molecular shape (e.g. compact molecules), there is no need of considering extended or quasi-one-dimensional molecular systems as in linear-scaling approaches. The third section discusses the overall performance of the code and illustrates how much optimizing the single-core performance of the specific processor at hand can be advantageous. The fourth section is devoted to the way our massively parallel simulations are deployed on a general computational platform and, particularly, how fault-tolerance is implemented, a crucial property for any large-scale simulation. Finally, a summary of the various strategies proposed in this work is presented in the last section.

II. THE QUANTUM MONTE CARLO METHOD

In this work we shall consider a variant of the Fixed-Node Diffusion Monte Carlo (FN-DMC) approach, the standard quantum Monte Carlo method used in computational chemistry. Here, we shall insist only on the aspects needed for understanding the rest of the work. For a complete presentation of the FN-DMC method the reader is referred, e.g to [2], [START_REF] Hammond | Monte Carlo Methods in Ab Initio Quantum Chemistry[END_REF], or [START_REF] Caffarel | Encyclopedia of Applied and Computational Mathematics[END_REF] and references therein.

A. Fixed-node diffusion Monte Carlo (FN-DMC) Diffusion Monte Carlo. In a diffusion Monte Carlo scheme, a finite population of "configurations" or "walkers" moving in the 3N -dimensional space (N , number of electrons) is introduced. A walker is described by a 3Ndimensional vector R ≡ (r 1 , r 2 , ..., r N ) giving the positions of the N electrons. At each Monte Carlo step each walker of the population is diffused and drifted according to

R ′ = R + τ b(R) + √ τ η (1)
where τ is a small time-step, η a Gaussian vector (3N independent normally distributed components simulating a free Brownian diffusion), and b(R) the drift vector given by

b(R) ≡ ∇ψ T (R) ψ T (R) , (2)
where ψ T , the trial wave function, is a known computable approximation of the exact wavefunction. At the end of this drift/diffusion step each walker is killed, kept unchanged, or duplicated a certain number of times proportionally to the branching weight w given by

w = e -τ 2 [(EL(R ′ )-ET )+(EL(R)-ET )] (3)
where E T is some reference energy and E L the local energy defined as

E L (R) ≡ Hψ T (R) ψ T (R) . ( 4
)
The population is propagated and after some equilibrium time it enters a stationary regime where averages are evaluated. As an important example, the exact energy may be obtained as the average of the local energy.

The Fixed-Node approximation.

Apart from the statistical and the short-time (finite time-step) errors which can be made arbitrary small, the only systematic error left in a DMC simulation is the so-called Fixed-Node (FN) error. This error results from the fact that the nodes of the trial wavefunction [defined as the (3N -1)-dimensional hyper-surface where Ψ T (R) = 0] act as infinitely repulsive barriers for the walkers [divergence of the drift vector, Eq.( 2)]. Each walker is thus trapped forever within the nodal pocket delimited by the nodes of Ψ T where it starts from. When the nodes of ψ T coincide with the exact nodes, the algorithm is exact. If not, a variational fixed-node error is introduced. However, with the standard trial wavefunctions used, this error is in general small,[13] a few percent of the correlation energy for total energies.

B. Parallelizing FN-DMC

Each Monte Carlo step is carried out independently for each walker of the population. The algorithm can thus be easily parallelized over an arbitrary number of processors by distributing the walkers among the processors, but doing this implies synchronizations of the CPUs since the branching step requires that all the walkers have first finished their drifted-diffusion step.

To avoid this aspect, we have chosen to let each CPU core manage its own population of walkers without any communication between the populations. On each computing unit a population of walkers is propagated and the various averages of interest are evaluated. At the end of the simulation, the averages obtained on each processor are collected and summed up to give the final answers. Regarding parallelism the situation is thus ideal since, apart from the negligible initial/final data transfers, there are no communications among processors.

The only practical problem left with FN-DMC is that the branching process causes fluctuations in the population size and thus may lead to load-balancing problem among processors. More precisely, nothing prevents the population size from decreasing or increasing indefinitely during the Monte Carlo iterations. To escape from this, a common solution consists in forcing the number of walkers not to deviate too much from some target value for the population size by introducing a population control step. It is usually realized by monitoring in time the value of the reference energy E T via a feedback mechanism, see e.g. [START_REF] Umrigar | [END_REF]. The price to pay is the introduction of some transient load imbalances and inter-processor communications/synchronization to redistribute walkers among computing cores, inevitably degrading the parallel speed-up. This solution has been adapted by several groups and some tricks have been proposed to keep this problem under control. [START_REF] Esler | [END_REF][START_REF] Kim | Proceedings of the 2010 Scientific Discovery through Ad-vanced Computing (SciDac) Conference[END_REF][START_REF] Gillan | Psi-k Highlight of the Month[END_REF][START_REF] Krogel | Advances in Quantum Monte Carlo[END_REF] Here, we propose to avoid this problem directly from the beginning by employing a variant of the FN-DMC working with a constant number of walkers. Several proposals can be found in the literature, e.g. [16,17]. Here, we shall employ the method described in Ref. [17]. In this approach the branching step of standard DMC is replaced by a so-called reconfiguration step. Defining the normalized branching weights as follows

p k = w k M i=1 w i (5)
the population of walkers is "reconfigured" by drawing at each step M walkers among the M walkers according to the probabilities p k . At infinite population, this step reduces to the standard branching step where walkers are deleted or duplicated proportionally to the weight w. At finite M the normalization factor M i=1 w i is no longer constant and a finite-population bias is introduced. To remove this error and recover the exact averages, a global weight given as the product of the total population weights for all preceding generations must be included into the averages. Note that this algorithm enables the possibility to use small walker populations on each core since there is no finite-population bias (typically, we use 10 to 100 walkers per core). For all details the reader is referred to [17].

C. Critical CPU part

At each Monte Carlo step the CPU effort is almost completely dominated by the evaluation of the wavefunction Ψ T and its first and second derivatives (computational hot spot). More precisely, for each walker the values of the trial wavefunction, Ψ T , its first derivatives with respect to all 3N -coordinates [drift vector, Eq.( 2)], and its Laplacian ∇ 2 Ψ T [kinetic part of the local energy, Eq.( 4)] are to be calculated. It is essential that such calculations be as efficient as possible since in realistic applications their number may be very large (typically of the order of 10 9 -10 12 ).

A common form for the trial wavefunction is

Ψ T (R) = e J(R) K=(K ↑ ,K ↓ ) c K Det K ↑ (r 1 , ..., r N ↑ )Det K ↓ (r N ↑ +1 , ..., r N ). ( 6
)
where the electron coordinates of the N ↑ (resp. N ↓ ) electrons of spin ↑ (resp. ↓) have been distinguished, N = N ↑ + N ↓ . In this formula e J(R) is the Jastrow factor describing explicitly the electron-electron interactions at different levels of approximations. A quite general form may be written as

J(R) = α U (e-n) (r iα ) + i,j
U (e-e) (r ij )

+ αi,j

U (e-e-n) (r ij , r iα , r jα ) + ... where r ij = |r ir j | is the inter-electronic distance and r iα = |r i -Q α |, the distance between electron i and nucleus α located at Q α . Here U 's are simple functions and various expressions have been employed in the literature. The Jastrow factor being essentially local, short-ranged expressions can be employed and the calculation of this term is usually a small contribution to the total computational cost. As a consequence, we shall not discuss further the computational aspect of this term here.

The second part of the wavefunction describes the shell-structure in terms of single-electron molecular orbitals and is written as a linear combination of products of two Slater determinants, one for the ↑ electrons and the other for the ↓ electrons. Each Slater matrix is built from a set of molecular orbitals φ i (r) usually obtained from a preliminary DFT or SCF calculations. The N orb molecular orbitals (MOs) are expressed as a sum over a finite set of N basis basis functions (atomic orbitals, AOs)

φ i (r) = N basis j=1 a ij χ j (r) (8)
where the basis functions χ j (r) are usually expressed as a product of a polynomial and a linear combination of Gaussian functions. In the present work the following standard form is employed

χ(r) = (x -Q x ) nx (y -Q y ) ny (z -Q z ) nz g(r) (9)
with

g(r) = k c k e -γ k (r-Q) 2 . ( 10
)
Here Q = (Q x , Q y , Q z ) is the vector position of the nucleus-center of the basis function, n = (n x , n y , n z ) a triplet of positive integers, g(r) is the spherical Gaussian component of the AO, and γ k its exponents. The determinants corresponding to spin ↑-electrons are expressed as

Det K ↑ (r 1 , ..., r N ↑ ) = Det    φ i1 (r 1 ) . . . φ i1 (r N ↑ ) . . . . . . . . . φ iN ↑ (r 1 ) . . . φ iN ↑ (r N ↑ )    (11
) where K ↑ is a compact notation for denoting the set of indices {i 1 , ..., i N ↑ } specifying the subset of the molecular orbitals used for this particular Slater matrix. A similar expression is written for spin ↓-electrons.

In contrast to the calculation of the Jastrow factor, the evaluation of the determinantal part of the wavefunction and its derivatives is critical. To perform such calculations we employ a standard approach [START_REF] Hammond | Monte Carlo Methods in Ab Initio Quantum Chemistry[END_REF] consisting in calculating the matrices of the first and second (diagonal) derivatives of each molecular orbital φ i with respect to the three space variables l = x, y, z evaluated for each electron position r j , namely

D (1) l,ij ≡ ∂φ i (r j ) ∂x j l (12) D (2) l,ij ≡ ∂ 2 φ i (r j ) ∂x j l 2 (13)
and then computing the inverse D -1 of the Slater matrix defined as D ij = φ i (r j ). The drift components and the Laplacian corresponding to the determinantal part of the trial wavefunction are thus evaluated as simple vectorproducts

1 Det(R) ∂Det(R) ∂x i l = j=1,N D (1)
l,ij D -1 ji (14) 1 Det(R) ∂ 2 Det(R) ∂x i l 2 = j=1,N D (2) l,ij D -1 ji (15)
From a numerical point of view, the computational time T needed to evaluate such quantities as a function of the number of electrons N scales as O(N 3 )

T = αN 3 + βN 3 . ( 16
)
The first N 3 -term results from the fact that the N 2 matrix elements of the Slater matrices are to be computed, each element being expressed in terms of the N basis ∼ N basis functions needed to reproduce an arbitrary delocalized molecular orbital. The second N 3 -term is associated with the generic cubic scaling of any linear algebra method for inverting a general matrix.

III. EXPLOITING THE HIGHLY LOCALIZED CHARACTER OF ATOMIC BASIS FUNCTIONS

As seen in the previous section, one of the two computational hot spots of QMC is the calculation of the derivatives of the determinantal part of the trial wave function for each electronic configuration (r 1 , ..., r N ) at each Monte Carlo step. To be more precise, the N orb molecular orbitals (MO) used in the determinantal expansion (6) are to be computed (here, their values will be denoted as C 1 ) together with their first derivatives with respect to x, y, and z (denoted C 2 , C 3 , C 4 ) and their Laplacians (denoted C 5 ). Calculations are made in single precision using an efficient matrix product routine we describe now. The matrix products involve the matrix of the MO coefficients a ij , Eq.( 8) (here denoted as A) the matrix of the atomic Gaussian basis functions evaluated at all electronic positions, χ j (r i ) (denoted B 1 ), their first derivatives (denoted B 2 , B 3 , B 4 ) and Laplacians (denoted B 5 ). The five matrix products are written under the convenient form

C i = AB i i = 1, 5 (17)
Note that matrix A remains constant during the simulation while matrices B i and C i depend on electronic configurations. The matrix sizes are as follows: N orb × N for the C i 's, N orb × N basis for A, and N basis × N for B. In practical applications N orb is of the order of N while N basis is greater than N by a factor 2 or 3 for standard calculations and much more when using highquality larger basis sets. The expensive part is essentially dominated by the N basis multiplications. The total computational effort is thus of order N orb × N × N basis , i.e. ∼ O(N 3 ).

The standard approach proposed in the literature for reducing the N 3 -price is to resort to the so-called linearscaling or O(N )-techniques. [18][19][20][21][22][23] The basic idea consists in introducing spatially localized molecular orbitals instead of the standard delocalized (canonical) ones obtained from diagonalization of reference Hamiltonians (usually, Hartree-Fock or Kohn-Sham). Since localized orbitals take their value in a finite region of space -usually in the vicinity of a fragment of the molecule -the number of basis set functions N basis needed to represent them with sufficient accuracy becomes essentially independent of the system size (not scaling with N as in the case of canonical ones). In addition to this, each electron contributes only to a small subset of the localized orbitals (those non-vanishing in the region where the electron is located). As a consequence, the number of non-vanishing matrix elements of the C i matrices no longer scales as N orb × N ∼ N 2 but linearly with N . Furthermore, each matrix element whose computation was proportional to the number of basis set used, N basis ∼ N , is now calculated in a finite time independent of the system size. Putting together these two results, we are led to a linear dependence of the computation of the C i matrices upon the number of electrons.

Here, we choose to follow a different path. Instead of localizing the canonical molecular orbitals we propose to take advantage of the localized character of the underlying atomic Gaussian basis set functions. The advantages are essentially three-fold:

1. the atomic basis set functions are naturally localized independently of the shape of the molecule. This is the most important point since the localization procedures are known to be effective for chemical systems having a molecular shape made of well-separated sub-units (for example, linear systems) but much less for general compact molecular systems that are ubiquitous in chemistry.

the degree of localization of the standard atomic

Gaussian functions is much larger than that obtained for molecular orbitals after localization (see results below)

3. by using the product form, Eq.( 17), the localized nature of the atomic Gaussian functions can be exploited very efficiently (see next section).

In practice, when the value of the spherical Gaussian part g(r) of an atomic orbital function χ(r) is smaller than a given threshold ǫ = 10 -8 , the value of the AO, its gradients and Laplacian are considered null. This property is used to consider the matrices B 1 , . . . , B 5 as sparse. However, in contrast with linear-scaling approaches, the MO matrix A is not considered here as sparse. We shall come back to this point later. To accelerate the calculations, an atomic radius is computed as the distance beyond which all the Gaussian components g(r) of the atomic orbitals χ(r) centered on the nucleus are less than ǫ. If an electron is farther than the atomic radius, all the AO values, gradients and Laplacians centered on the nucleus are set to zero.

The practical implementation to perform the matrix products is as follows. For each electron, the list of indices (array indices in what follows) where g(r) > 0 is calculated. Then, the practical algorithm can be written as (where x += y denotes x = x + y).

C1 = 0. C2 = 0. C3 = 0. C4 = 0. C5 = 0. do i=1, Number of electrons do k=1, Number of non-zero AOs for electron i do j=1, Number of molecular orbitals C1(j,i) += A(j,indices(k,i))*B1(k,i) C2(j,i) += A(j,indices(k,i))*B2(k,i) C3(j,i) += A(j,indices(k,i))*B3(k,i) C4(j,i) += A(j,indices(k,i))*B4(k,i) C5(j,i) += A(j,indices(k,i))*B5(k,i) end do end do end do
This implementation allows to take account of the sparsity of the B matrices, while keeping the efficiency due to a possible vectorization of the inner loop. The load/store ratio is 6/5 (6 load-from-memory instructions, 5 store-to-memory instructions) in the inner loop : the elements of B n are constant in the inner loop (in registers), and the same element of A is used at each line of the inner loop (loaded once per loop cycle). As store operations are more expensive than load operations, increasing the load/store ratio improves performance as will be shown in the next section. Using this algorithm, the scaling of the matrix products is expected to drop from O(N 3 ) to a scaling roughly equal to O(N 2 ) (in a regime where N is large enough, see discussion in the next section). Let us now illustrate such a property in the applications to follow.

The different systems used here as benchmarks are represented in Figure 1. The trial wavefunctions used for describing each system are standard Hartree-Fock wavefunctions (no Jastrow factor) with molecular orbitals expressed using various Gaussian basis sets. System 1 is a copper complex with four ligands having 158 electrons and described with a cc-pVDZ basis set. System 2 is a polypeptide taken from reference [24] (434 electrons and 6-31G * basis set). System 3 is identical to System 2 but using a larger basis set, namely the cc-pVTZ basis set. System 4 is the 1ZE7 molecule from the Protein Data Bank (1056 electrons, 6-31G * ), and System 5 is the 1AMB molecule from the Protein Data Bank (1731 electrons, 6-31G * ).

Table IV shows the level of sparsity of the matrices A (A ij ≡ a ij ) and B 1 (B 1ij ≡ χ i (r j )) for the five systems (matrices B n with n > 1 behave as B 1 with respect to sparsity). As seen the number of basis set functions employed is proportional to the number of electrons with a factor ranging from about 2.2 to 6.8.

Regarding the matrix A of MO coefficients the results are given both for standard canonical (delocalized) MOs and for localized orbitals. To get the latter ones, different localization schemes have been applied. [25][26][27] However, they essentially lead to similar results. Here, the results presented are those obtained by using the Cholesky decomposition of the density matrix expressed in the AO basis set. [27] As seen the level of sparsity of the matrix A is low. Although it increases here with the system size it remains modest for the largest size (there are still about one third of non-zero elements). Of course, such a result strongly depends on the type of molecular system considered (compact or not compact) and on the diffuse character of the atomic basis set. Here, we have considered typical systems of biochemistry.

Next, the level of sparsity of the B matrices is illustrated. The percentage of non-zero values of χ i (r j ) has been obtained as an average over a Variational Monte Carlo run. In sharp contrast with MOs the atomic orbitals are much more localized, thus leading to a high level of sparsity. For the largest system, only 3.9% of the basis function values are non-negligible.

In the last column of the table the maximum number of non-zero elements obtained for all columns of the matrix (i.e, a value of electron position) during the entire Monte Carlo simulation is given. A first remark is that this number is roughly constant for all system sizes. A second remark is that the percentage of non-zero values is only slightly greater than the average, thus showing that the B matrices can be considered sparse during the whole simulation, not only in average. As an important consequence, the loop over the number of non-zero AOs for each electron in the practical algorithm presented above (loop over k index) is expected to be roughly constant as a function of the size at each Monte Carlo step. This latter remark implies for this part an expected behavior of order O(N 2 ) for large N . Let us now have a closer look at the actual performance of the code.

IV. OVERALL PERFORMANCE OF QMC=CHEM

When discussing performance several aspects must be considered. A first one, which is traditionally discussed, is the formal scaling of the code as a function of the system size N (N ∼ number of electrons). As already noted, due to the innermost calculation, products, and inversion of matrices, such a scaling is expected to be cubic, O(N 3 ). However, there is a second important aspect, generally not discussed, which is related to the way the expensive innermost floating-point operations are implemented and on how far and how efficiently the potential performance of the processor at hand is exploited. In what follows we shall refer to this aspect as "single-core optimization". It is important to emphasize that such an aspect is by no way minor and independent on the previous "mathematical" one.

To explicit this point, let us first recall that the compu-tational time T results essentially from two independent parts, the first one resulting from the computation of the matrix elements, T 1 ∼ αN 3 and the second one from the inversion of the Slater matrix, T 2 ∼ βN 3 . Now, let us imagine that we have been capable of devising a highly efficient linear-scaling algorithm for the first contribution such that T 1 ∼ ǫN << T 2 within the whole range of system sizes N considered. We would naturally conclude that the overall computational cost T ∼ T 2 is cubic. In the opposite case where a very inefficient linear-scaling algorithm is employed for the first part, T ∼ T 1 ≫ T 2 , we would conclude to a linear-scaling type behavior. Of course, mathematically speaking such a way of reasoning is not correct since scaling laws are only meaningful in the asymptotic regime where N goes to infinity. However, in practice only a finite range of sizes is considered (here, between 2 and about 2000 active electrons) and it is important to be very cautious with the notion of scaling laws. A more correct point of view consists in looking at the global performance of the code in terms of total CPU time for a given range of system sizes, a given compiler, and a given type of CPU core. Finally, a last aspect concerns the memory footprint of the code whose minimization turns out to be very advantageous. Indeed, the current trend in supercomputer design is to increase the number of cores more rapidly than the available total memory. As the amount of memory per core will continue to decrease, it is very likely that programs will need to have a low memory footprint to take advantage of exascale computers. Another point is that when less memory is used less electrical power is needed to perform the calculation: data movement from the memory modules to the cores needs more electrical power than performing floating point operations. Although at present time the power consumption is not yet a concern to software developers, it is a key aspect in present design of the exascale machines to come.

In this section, the results discussed will be systematically presented by using two different generations of Intel Xeon processors. The first processor, referred to as Core2, is an Intel Xeon 5140, Core2 2.33 GHz, Dual core, 4 MiB shared L2 cache. The second one, referred to as Sandy Bridge, is an Intel Xeon E3-1240 at 3.30 GHz, Quad core, 256 KiB L2 cache/core, 8 MiB shared L3 cache (3.4 GHz with turbo). Note also that the parallel scaling of QMC being close to ideal (see next section), single-core optimization is very interesting: the gain in execution time obtained on the single-core executable is directly transferred to the whole parallel simulation.

A. Improving the innermost expensive floating-point operations

For the Core2 architecture, the practical algorithm presented above may be further improved by first using the unroll and jam technique, [START_REF] Zingirian | High-Performance Computing and Networking[END_REF] which consists in unrolling the outer loop and merging multiple outer-loop iterations in the inner loop : For the Sandy Bridge architecture, the external body is unrolled four times instead of two, and the most internal loop is split in three loops: one loop computing C 1 , C 2 , a second loop computing C 3 , C 4 , and a third loop computing C 5 . The load/store ratio is improved from 6/5 to 6/2 and 5/1.

Then, all arrays were 256-bit aligned using compiler directives and the first dimensions of all arrays were set to a multiple of 8 elements (if necessary, padded with zeros at the end of each column) in order to force a 256bit alignment of every column of the matrices. These modifications allowed the compiler to use only vector instructions to perform the matrix products, both with the Streaming SIMD Extension (SSE) or the Advanced Vector Extension (AVX) instruction sets. The x86 64 version of the MAQAO framework [START_REF] Djoudi | Workshop on EPIC Architectures and Compiler Technology[END_REF] indicates that, as the compiler unrolled twice the third loop (C 5 ), these three loops perform 16 floating point operations per cycle, which is the peak performance on architecture.

Finally, to improve the cache hit probability, blocking was used on the first dimension of B n (loop over k). In each block, the electrons (columns of B) are sorted by ascending first element of the indices array in the block. This increases the probability that columns of A will be in the cache for the computation of the values associated with the next electron.

The results obtained using the Intel Fortran Compiler XE 2011 are presented in table I for both the Core2 and the Sandy Bridge architectures. The single-core doubleprecision Linpack benchmark is also mentioned for comparison. The results show that the full performance of the matrix products is already reached for the smallest system. However, as opposed to dense matrix product routines, we could not approach further the peak performance of the processor since the number of memory accesses scales as the number of floating point operations (both O(N 2 )): the limiting factor is inevitably the data access. Nevertheless, the DECAN tool [START_REF] Koliai | [END_REF] revealed that data access only adds a 30% penalty on the pure arithmetic time, indicating an excellent use of the hierarchical memory and the prefetchers. In table II the memory required together with the CPU time obtained for one Monte Carlo (VMC) step for the five systems are presented using both processors. The two expensive computational parts (matrix products and inversion) are distinguished. A first remark is that the trends for both processors are very similar so we do not need to make a distinction at this moment. A second remark is that the memory footprint of QMC=Chem is particularly low. For the biggest size considered (1731 electrons) the amount of RAM needed is only 313 MiB. Finally, another important remark is that at small number of electrons the multiplicative part is dominant while this is not the case at larger sizes. Here, the change of regime is observed somewhere between 400 and 1000 electrons but its precise location depends strongly on the number of basis functions used. For example, for systems 3 and 4 corresponding to the same molecule with a different number of basis functions, the multiplicative part is still dominant for the larger basis set (β-strand with cc-pVTZ) while it is no longer true for the smaller basis set (β-strand with 6-31G * ). In figure 2 a plot of the total computational time for the Sandy Bridge core as a function of the number of electrons is presented. A standard fit of the curve with a polynomial form N γ leads to a γ-value of about 2.5. However, as discussed above such a power is not really meaningful. From the data of table II it is easy to extract the pure contribution related to the inversion and a factor very close to 3 is obtained, thus illustrating that for this linear algebra part we are in the asymptotic regime. For the multiplicative part, the pure N 2 behavior is not yet recovered and we are in an intermediate regime. Putting together these two situations leads to some intermediate scaling around 2.5.

Sparsity

In our practical algorithm, for the matrix products we have chosen to consider the B matrices as sparse as opposed to the A matrix which is considered dense. The reason for that is that considering the matrix A sparse would not allow us to write a stride-one inner loop. In single precision, SSE instructions executed on Intel processors can perform up to 8 instructions per CPU cycle (one 4-element vector ADD instruction and one 4element vector MUL instruction in parallel). Using the latest AVX instruction set available on the Sandy Bridge architecture, the width of the SIMD vector registers have been doubled and the CPU can now perform up to 16 floating point operations per cycle. A necessary condition for enabling vectorization is a stride-one access to the data. This implies that using a sparse representation of A would disable vectorization, and reduce the maximum number of floating operations per cycle by a factor of 4 using SSE (respectively 8 using AVX). If matrix A has more than 25% (respectively 12.5%) non-zero elements, using a sparse representation is clearly not the best choice. This last result is a nice illustration of the idea that the efficiency of the formal mathematical algorithm depends on the core architecture.

Inversion step

Now, let us consider the inversion step which is the dominant CPU-part for the big enough systems (here, for about a thousand electrons and more). In Table I the performance in GFlops/s of the inversion step is presented for both processors. For comparisons the theoretical single-core peak and single-core Linpack performance are given. The first column gives the overall performance of the code while the second column is specific to the inversion part. As seen the performance of both parts increases with the number of electrons. For largest systems the performance represents more than 50% of the peak performance of each processor. For the largest system the whole code has a performance of about 54% of the peak performance for the Core2, and about 61% for the Sandy Bridge. The performance is still better for the inversion part: 60% for the Core2 and 67% for the Sandy Bridge.

Determinant calculation compared to spline interpolation

Most authors use three-dimensional spline representations of the molecular orbitals in order to compute in constant time the values, first derivatives and Laplacians of one electron in one molecular orbital, independently of the size of the atomic basis set. This approach seems efficient at first sight, but the major draw-back is that the memory required for a single processor can become rapidly prohibitive since each molecular orbital has to be pre-computed on a three-dimensional grid. To overcome the large-memory problem, these authors use shared memory approaches on the computing nodes, which implies coupling between the different CPU cores. In this paragraph, we compare the wall time needed for spline interpolation or computation of the values, first derivatives and Laplacians of the wave function at all electron positions. Version 0.9.2 of the Einspline package [31] was used as a reference to compute the interpolated values, gradients and Laplacians of 128 MOs represented on 23 × 21 × 29 single precision arrays. The "multiple uniform splines" set of routines were used. To evaluate the value, gradient and Laplacian of one molecular orbital at one electron coordinate, an average of 1200 CPU cycles was measured using LIKWID [START_REF] Treibig | 39th International Conference on Parallel Processing Workshops (ICPPW)[END_REF] on the Core2 processor versus 850 CPU cycles on the Sandy Bridge processor. Even if the interpolation is done using a very small amount of data and of floating point operations, it is bound by the memory latency. Indeed, the needed data is very unlikely to be in the CPU cache and this explains why the number of cycles per matrix element is quite large. As our code uses a very small amount of memory, and as the computationally intensive routines are very well vectorized by the compiler, the computation of the matrix elements is bound by the floating point throughput of the processor.

The number of cycles needed to build the C 1 . . . C 5 matrices is the number of cycles needed for one matrix element scaled by the number of matrix elements N 2 α + N 2 β . Table III shows the number of CPU cycles needed to build the full C 1 . . . C 5 matrices for a new set of electron positions using spline interpolation or using computation. The computation includes the computation of the values, gradients and Laplacians of the atomic orbitals (matrices B 1 . . . B 5 ) followed by the matrix products.

Using a rather small basis set (6-31G * ), the computation of the matrices in the 158-electron system is only 10% slower than the interpolation on the Core2 architecture. Using a larger basis set (cc-pVTZ), the computation is only 57% slower.

As the frequency is higher in our Sandy Bridge processor than in our Core2 processor, we would have expected the number of cycles of one memory latency to increase, and therefore we would have expected the Einspline package to be less efficient on that specific processor. One can remark that the memory latencies have been dramatically improved from the Core2 to the Sandy Bridge architectures and the number of cycles for the interpolation decreases.

The full computation of the matrix elements benefits from the improvement in the memory accesses, but also from the enlargement of the vector registers from 128 bits to 256 bits. This higher vectorization considerably reduces the number of cycles needed to perform the calculation such that in the worst case (the largest basis set), the full computation of the matrix elements takes as much time as the interpolation. In all other cases, the computation is faster than the spline interpolation. Finally, let us mention that as the memory controller is directly attached to the CPU, on multi-socket computing nodes the memory latencies are higher when accessing a memory module attached to another CPU (NUMA architecture).

V. PARALLELISM: IMPLEMENTING A UNIVERSAL, DYNAMIC, AND FAULT-TOLERANT SCHEME

Our objective was to design a program that could take maximum advantage of heterogeneous clusters, grid environments, the petaflops platforms available now and those to come soon (exascale).

To achieve the best possible parallel speed-up on any hardware, all the parallel tasks have to be completely decoupled. Feldman et al have shown that a naive implementation of parallelism does not scale well on commodity hardware. [START_REF] Feldman | [END_REF] Such bad scalings are also expected to be observed on very large scale simulations. Therefore, we chose an implementation where each CPU core realizes a QMC run with its own population of walkers independently of all the other CPU cores. The run is divided in blocks over which the averages of the quantities of interest are computed. The only mandatory communications are the one-to-all communication of the input data and the all-to-one communications of the results, each result being the Monte Carlo average computed with a single-core executable. If a single-core executable is able to start as soon as the input data is available and stop at any time sending an average over all the computed Monte Carlo steps, the best possible parallel speed-up on the machine can always be obtained. This aspect is detailed in this section.

A. Fault-tolerance

Fault-tolerance is a critical aspect since the mean time before failure (MTBF) increases with the number of hardware components : using N identical computing nodes for a singe run multiplies by N the probability of failure of the run. If one computing node is expected to fail once a year, a run using 365 computing nodes is not expected to last more than a day. As our goal is the use both of massive resources and commodity clusters found in laboratories, hardware failure is at the center of our software design.

The traditional choice for the implementation of parallelism is the use of the Message Passing Interface (MPI). [34] Efficient libraries are proposed on every parallel machine, and it is probably the best choice in most situations. However, all the complex features of MPI are not needed for our QMC program, and it does not really fit our needs: in the usual MPI implementations, the whole run is killed when one parallel task is known not be able to reach the MPI Finalize statement. This situation occurs when a parallel task is killed, often due to a system failure (I/O error, frozen computing node, hardware failure, etc). For deterministic calculations where the result of every parallel task is required, this mechanism prevents from unexpected dead locks by immediately stopping a calculation that will never end. In our implementation, as the result of the calculation of a block is a Gaussian distributed random variable, removing the result of a block from the simulation is not a problem since doing that does not introduce any bias in the final result. Therefore, if one computing node fails, the rest of the simulation should survive.

We wrote a simple Python TCP client/server application to handle parallelism. To artificially improve the bandwidth, all network transfers are compressed using the Zlib library,[35] and the results are transferred asynchronously in large packets containing a collection of small messages. Similarly, the storage of the results is executed using a non-blocking mechanism. The computationally intensive parts were written using the IRPF90 code generator, [36] to produce efficient Fortran code that is also easy to maintain. The architecture of the whole program is displayed in figure 3.

B. Program interface

Our choice concerning the interaction of the user with the program was not to use the usual "input file and output file" structure. Instead, we chose to use a database containing all the input data and control parameters of the simulation, and also the results computed by different runs. A few simple scripts allow the interaction of the user with the database. This choice has several advantages:

• The input and output data are tightly linked together. It is always possible to find to which input corresponds output data.

• If an output file is needed, it can be generated on demand using different levels of verbosity

• Graphical and web interfaces can be trivially connected to the program

• Simple scripts can be written by the users to manipulate the computed data in a way suiting their needs.

Instead of storing the running average as the output of a run, we store all the independent block-averages in the database, and the running averages are post-processed on demand by database queries. There are multiple benefits from this choice:

• Checkpoint/restart is always available

• It is possible to compute correlations, combine different random variables, etc, even when the QMC run is finished

• Combining results computed on different clusters consists in simply merging the two databases, which allows automatically the use of the program on computing grids [START_REF] Monari | Remote Instrumentation for eScience and Related Aspects[END_REF] • Multiple independent jobs running on the same cluster can read/write in the same database to communicate via the file system. This allows to gather more and more resources as they become available on a cluster, or to run a massive number of tasks in a best effort mode. [START_REF]When cluster resources are unused, a QMC job starts[END_REF] C. Error-checking

We define the critical data of a simulation as the input data that characterizes uniquely a given simulation. For instance, the molecular coordinates, the molecular orbitals, the Jastrow factor parameters are critical data since they are fixed parameters of the wave function during a QMC run. In contrast, the number of walkers of a simulation is not critical data for a VMC run since the results of two VMC simulations with a different number of walkers can be combined together. A 32-bit cyclic redundancy code (CRC-32 key) is associated with the critical data to characterize a simulation. This key will be used to guarantee that the results obtained in one simulation will never be mixed with the results coming from another simulation and corrupt the database. It will also be used to check that the input data has been well transferred on every computing node.

D. Program execution

When the program starts its execution, the manager process runs on the master node and spawns two other processes: a data server and a main worker process.

At any time, new clients can connect to the data server to add dynamically more computational resources to a running calculation, and some running clients can be terminated without stopping the whole calculation. The manager periodically queries the database and computes the running averages using all the blocks stored in the database. It controls the running/stopping state of the workers by checking if the stopping condition is reached (based for example on the wall-clock time, on the error bar of the average energy, a Unix signal, etc).

When running on super-computers, the main worker process spawns one single instance of a forwarder on each computing node given by the batch scheduler system using an MPI launcher. As soon as the forwarders are started the MPI launcher terminates, and each forwarder connects to the data server to retrieve the needed input data. The forwarder then starts multiple workers on the node with different initial walker positions.

Each worker is an instance of the single-core Fortran executable, connected to the forwarder by Unix pipes. Its behavior is the following : while ( .True. ) { compute_a_block_of_data(); send_the_results_to_the_forwarder(); } Unix signals SIGTERM and SIGUSR2 are trapped to trigger the send_the_results_to_the_forwarder procedure followed by the termination of the process. Using this mechanism, any single-core executable can be stopped immediately without losing a single Monte Carlo step. This aspect is essential to obtain the best possible speed-up on massively parallel machines. Indeed, using the matrix product presented in the previous section makes the CPU time of a block non-constant. Without this mechanism the run would finish when the last CPU finishes, and the parallel efficiency would be reduced when using a very large number of CPU cores.

While the workers are computing the next block, the forwarder sends the current results to the data-server using a path going through other forwarders. The forwarders are organized in a binary tree as displayed in figure 4: every node of the tree can send data to all its ancestors, in order to deal with possible failures of computing nodes. This tree-organization reduces the number of connections to the data server, and also enlarges the size of the messages by combining in a single message the results of many forwarders.

At the end of each block, the last walker positions are sent from the worker to the forwarder. The forwarder keeps a fixed-sized list of N kept walkers enforcing the distribution of local energies : when a forwarder receives a set of N walkers, it appends the list of new walkers to its N kept list, and sorts the N kept + N list by increasing local energies. A random number η is drawn to keep all list entries at indices ⌊η + i(N kept + N )/N kept ⌋, i = {1, . . . , N kept } After a random timeout, if the forwarder is idle, it sends its list of walkers to its parent in the binary tree which repeats the list merging process. Finally, the data server receives a list of walkers, merges it with its own list and writes it to disk when idle. This mechanism ensures that the walkers saved to disk will represent homogeneously the whole run and avoids sending all the walkers to the data server. These walkers will be used as new starting points for the next QMC run.

Using such a design the program is robust to system failures. Any computing node can fail with a minimal impact on the simulation:

• If a worker process fails, only the block being computed by this worker is lost. It does not affect the forwarder to which it is linked.

• If a forwarder fails, then only one computing node is lost thanks to the redundancy introduced in the binary tree of forwarders.

• The program execution survives short network disruption (a fixed timeout parameter). The data will arrive to the data server when the network becomes operational again.

• The disks can crash on the computing nodes : the temporary directory used on the computing nodes is a RAM-disks (/dev/shm).

• The shared file system can fail as the single-core static executable, the python scripts and input files are broadcast to the RAM-disks of the compute nodes with the MPI launcher when the run starts

• Redundancy can be introduced on the data server by running multiple jobs using the same database. Upon a failure of a data server, only the forwarders connected to it will be lost.

• In the case of a general power failure, all the calculations can be restarted without losing what has already been stored in the database.

Finally, we have left the possibility of using different executables connected to the same forwarder. This will allow a combined use of pure CPU executables with hybrid CPU/GPU and CPU/MIC executables, to use efficiently all the available hardware. The extension to hybrid architectures will be the object of a future work.

E. Parallel speed-up

The benchmarks presented in this section were performed on the Curie machine (GENCI-TGCC-CEA, France). Each computing node is a dual socket Intel Xeon E5-2680: 2×(8 cores, 20 MiB shared L3-cache, 2.7 GHz) with 64 GiB of RAM. The benchmark is a DMC calculation of the β-strand system with the cc-PVTZ basis set (Table IV) using 100 walkers per core performing 300 steps in each block. Note that these blocks are very short compared to realistic simulations where the typical number of steps would be larger than 1000 to avoid the correlation between the block averages.

Intra-node

The CPU consumption of the forwarder is negligible (typically 1% of the CPU time spent in the single-core executables). The speed-up with respect to the number of sockets is ideal. Indeed the single-core binaries do not communicate between each other, and as the memory consumption per core is very low, each socket never uses memory modules attached to another socket. When multiple cores on the same socket are used we observed a slow-down for each core due to the sharing of the L3cache and memory modules. Running simultaneously 16 instances of the single-core binaries on our benchmark machine yields an increase of 10.7% of the wall-clock time compared to running only one instance. For a 16-core run, we obtained a 14.4× speed-up (the Turbo feature of the processors was de-activated for this benchmark).

Inter-node

In this section the wall-clock time is measured from the very beginning to the very end of the program execution using the standard GNU time tool. Hence, the wall-clock time includes the initialization and finalization steps.

The initialization step includes

• Input file consistency checking

• Creating a gzipped tar file containing the input files (wave function parameters, simulation parameters, a pool of initial walkers), the Python scripts and static single-core executable needed for the program execution on the slave nodes

• MPI Initialization

• Broadcasting the gzipped tar file via MPI to all the slave nodes

• Extracting the tar file to the RAM-disk of the slave nodes

• Starting the forwarders

• Starting the single-core instances Note that as no synchronization is needed between the nodes, the computation starts as soon as possible on each node. The finalization step occurs as follows. When the data server receives a termination signal, it sends a termination signal to all the forwarders that are leaves in the tree of forwarders. When a forwarder receives such a signal, it sends a SIGTERM signal to all the single-core binary instances of the computing node which terminate after sending to the forwarder the averages computed over the truncated block. Then, the forwarder sends this data to its parent in the binary tree with a termination signal, and sends a message to the data server to inform it that it is terminated. This termination step walks recursively through the tree. When all forwarders are done, the data server exits. Note that if a failure happened on a node during the run, the data server never receives the message corresponding to a termination of the corresponding forwarder. Therefore when the data server receives the termination signal coming from the forwarders tree, if the data server is still running after a given timeout it exits.

We prepared a 10-minutes run for this section to compute the parallel speed-up curve as a function of the number of 16-core nodes given in Figure 5. The data corresponding to this curve are given in Table V. The reference for the speed-up is the one-node run. The speed-up for N nodes is computed as:

t CPU (N )/t Wall (N ) t CPU (1)/t Wall (1) (18)
The initialization time was 9 seconds for the single node run, and 22 seconds for the 1000 nodes run. The finalization time was 13 seconds for the single node run, and 100 seconds for the 1000 nodes run. Apart from the initialization and finalization steps (which obviously do not depend on the total execution 5). CPU time is the cumulated CPU time spent only in the Fortran executables, Wall time is the measured wall-clock time, including initialization and finalization steps (serial). The 10-minutes run were measured, and the longer runs are estimated from the 10-minutes run data. Two checks were measured for the 60-minutes runs with 1 and 100 nodes (in italics). time), the parallel speed-up is ideal. This allowed us to estimate the speed-ups we would have obtained for a 1hour run and for a 3-hours run. For instance, to estimate the one-hour run we added 50 minutes to the wall-clock time and 50 minutes×16×number of nodes ×0.99 to the CPU time. The 99% factor takes account of the CPU consumption of the forwarder for communications. Our simple model was checked by performing a one-hour run on one node, and a one-hour run on 100 nodes. An excellent agreement with the prediction was found: a 99.5× speed-up was predicted for 100 nodes and a 99.3× speedup was measured.

Finally, a production run was made using 76,800 cores of Curie (4,800 nodes) on the β-strand molecule with a cc-pVTZ basis set via 12 runs of 40 nodes, and a sustained performance of 960 TFlops/s was measured. All the details and scientific results of this application will be presented elsewhere. [START_REF] Caffarel | [END_REF] VI. SUMMARY Let us summarize the main results of this work. First, to enhance the computational efficiency of the expensive innermost floating-point operations (calculation and multiplication of matrices) we propose to take advantage of the highly localized character of the atomic Gaussian basis functions, in contrast with the standard approaches using localized molecular orbitals. The advantages of relying on atomic localization have been illustrated on a series of molecules of increasing sizes (number of electrons ranging from 158 to 1731). In this work, it is emphasized that the notion of scaling of the computational cost as a function of the system size has to be considered with caution. Here, although the algorithm proposed is formally quadratic it displays a small enough prefactor to become very efficient in the range of number of electrons considered. Furthermore, our implementation of the linear-algebra computational part has allowed to enlighten a fundamental issue rarely discussed, namely the importance of taking into consideration the close links between algorithmic structure and CPU core architecture. Using efficient techniques and optimization tools for enhancing single-core performance, this point has been illustrated in various situations. Remark that this aspect is particularly important: as the parallel speed-up is very good, the gain in execution time obtained for the singlecore executable will also be effective in the total parallel simulation.

In our implementation we have chosen to minimize the memory footprint. This choice is justified first by the fact that today the amount of memory per CPU core tends to decrease and second by the fact that small memory footprints allow in general a more efficient usage of caches. In this spirit, we propose not to use 3D-spline representation of the molecular orbitals as usually done. We have shown that this can be realized without increasing the CPU cost. For our largest system with 1731 electrons, only 313 MiB of memory per core was required. As a consequence the key limiting factor of our code is only the available CPU time and neither the memory nor disk space requirements, nor the network performance. Let us re-emphasize that this feature is well aligned with the current trends in computer architecture for large HPC systems.

Finally, let us conclude by the fact that there is no fundamental reason why the implementation of such a QMC simulation environment which has been validated at petaflops level could not be extended to exascale.

FIG. 1 :

 1 FIG. 1: Molecular systems used as benchmarks.

  do i=1, Number of electrons do k=1, Number of non-zero AOs for electron i, 2 do j=1, Number of molecular orbitals C1(j,i) += A(j,indices(k ,i))*B1(k ,i) + & A(j,indices(k+1,i))*B1(k+1,i) C2(j,i) += A(j,indices(k ,i))*B2(k ,i) + & A(j,indices(k+1,i))*B2(k+1,i) ... end do end do end do To avoid register spilling, the inner loop is split in two loops : one loop computing C 1 , C 2 , C 3 and a second loop computing C 4 , C 5 . The load/store ratio is improved from 6/5 to 5/3 and 4/2

FIG. 2 :

 2 FIG.2: Single-core scaling with system size.

FIG. 5 :

 5 FIG. 5: Parallel speed-up of QMC=Chem with respect to 16-core compute nodes (Reference is one 16-core node).

TABLE II :

 II

			Core2		Sandy Bridge
		Products Inversion Overall Products Inversion Overall
	Linpack (DP)		7.9 (84.9%)			24.3 (92.0%)
	Peak	18.6	9.3		52.8	26.4
	Smallest system 9.8 (52.7%) 2.6 (28.0%)	3.3	26.6 (50.3%) 8.8 (33.3%)	6.3
	β-Strand	9.7 (52.2%) 4.3 (46.2%)	3.7	33.1 (62.7%) 13.7 (51.2%) 13.0
	β-Strand TZ	9.9 (53.2%) 4.3 (46.2%)	4.5	33.6 (63.6%) 13.7 (51.2%) 14.0
	1ZE7	9.3 (50.0%) 5.2 (55.9%)	4.6	30.6 (57.9%) 15.2 (57.6%) 17.9
	1AMB	9.2 (49.5%) 5.6 (60.2%)	5.0	28.2 (53.4%) 16.2 (61.4%) 17.8
	TABLE I: Single core performance (GFlops/s) of the matrix products (single precision), inversion (double precision) and overall
	performance of QMC=Chem (mixed single/double precision). The percentage of the peak performance is given in parentheses.
	Core2 : Intel Xeon 5140, Core2 2.33GHz, Dual core, 4MiB shared L2 cache. Sandy Bridge : Intel Xeon E3-1240, Sandy Bridge
	3.30GHz, Quad core, 256KiB L2 cache/core, 8MiB shared L3 cache (3.4GHz with turbo).
		Smallest system β-Strand β-Strand TZ 1ZE7 1AMB
	RAM (MiB)		9.8	31	65	133	313
	Core2
	QMC step(s)	0.0062	0.0391	0.0524	0.2723 0.9703
	Inversion		15%	31%	21%	47%	58%
	Products		25%	23%	35%	21%	18%
	Sandy Bridge
	QMC step(s)	0.0026	0.0119	0.0187	0.0860 0.3042
	Inversion		12%	26%	17%	42%	52%
	Products		24%	22%	32%	21%	20%

Single-core memory consumption and elapsed time for one VMC step. Values in % represent the percentage of the total CPU time. Core2 : Intel Xeon 5140, Core2 2.33GHz, Dual core, 4MiB shared L2 cache. Sandy Bridge : Intel Xeon E3-1240, Sandy Bridge 3.30GHz, Quad core, 256KiB L2 cache/core, 8MiB shared L3 cache (3.4GHz with turbo).

TABLE III :

 III Number of million CPU cycles needed for the computation of the values, gradients and Laplacians of the molecular orbitals using the Einspline package and using our implementation for the Core2 and the Sandy Bridge micro-architectures. The ratio QMC=Chem / Einspline is also given.

							Core2	Sandy Bridge
						QMC=Chem Einspline ratio QMC=Chem Einspline ratio
			Smallest system	16.7		15.0	1.11	9.2	10.6	0.87
			β-Strand TZ	177.3	113.0	1.57	81.7	80.1	1.02
			1ZE7		783.5	669.1	1.17	352.0	473.9	0.74
			1AMB		2603.0	1797.8	1.45	1183.9	1273.5	0.93
		Master compute node
			Data Server		Slave Compute node
	Manager	Main worker thread	Network Thread		Forwarder
			I/O Thread		Forwarder	Worker	Worker	Worker
		Database	Worker	Worker	Worker
	FIG. 3: Overview of the QMC=Chem architecture.

  Connections of the forwarders with the data server.

	Forwarder	Forwarder	Forwarder	Forwarder	Forwarder	Forwarder	Forwarder	Forwarder
	Forwarder		Forwarder			Forwarder	Forwarder
		Forwarder				Forwarder
					Forwarder
					Data Server
	FIG. 4:

TABLE IV :

 IV System sizes, percentage of non-zero molecular orbital coefficients and average percentage of non-zero Atomic Orbital values. a Zero MO coefficients are those below 10 -5 . These are given for localized orbitals, and for canonical orbitals in parentheses. b Zero AO matrix elements are those for which the radial component of the basis function has a value below 10 -8 for given electron positions.

	Number of		10 minutes			60 minutes		180 minutes
	16-core				(estimated)	(estimated)
	Nodes	CPU (s) Wall (s) Speed-up CPU(s) Wall(s) Speed-up Speed-up
	1	9627	625	1.0	57147	3625	1.00	1.00
					57332	3629	1.00
	10	95721	627	9.9	570921	3627	9.98	9.99
	25	239628	629	24.7	1427628	3629	24.95	24.98
	50	477295	631	49.1	2853295	3631	49.85	49.95
	100	952388	636	97.2	5704388	3636	99.52	99.84
					5708422	3638	99.32
	200	1869182	637	190.5	11373182 3637	198.36	199.45
	400	3725538	648	373.3	22733538 3648	395.30	398.42
	500	4479367	641	453.7	28239367 3641	491.98	497.31
	1000	8233981	713	749.7	55753981 3713	952.50	983.86

TABLE V :

 V Data relative to the scaling curve (Figure

Acknowledgments

AS and MC would like to thank ANR for support under Grant No ANR 2011 BS08 004 01. This work was possible thanks to the generous computational support from CALMIP (Université de Toulouse), under the allocation 2011-0510, GENCI under the allocation GEN1738, CCRT (CEA), and PRACE under the allocation RA0824. The authors would also like to thank Bull, GENCI and CEA for their help in this project.


Abbas Chamseddine

Youmin Zhang

Camille-Alain Alain Rabbath

Didier Theilliol

Trajectory planning and replanning strategies applied to a quadrotor unmanned aerial vehicle

published or not. The documents may come   L'archive ouverte pluridisciplinaire

I. Introduction

Path or trajectory planning is a complex problem. It involves meeting the physical constraints of the unmanned vehicles, constraints from the operating environment and other operational requirements [START_REF] Tsourdos | Cooperative Path Planning of Unmanned Aerial Vehicles[END_REF]. The quadrotor helicopter is one of the unmanned aerial vehicles (UAVs) that has been considered for the trajectory planning problem. However, most of the existing works are optimization-based methods. A time-optimal motion planning approach is proposed in [START_REF] Lai | Time-optimal Control of a Hovering Quad-rotor Helicopter[END_REF] for a non-linear model of a quadrotor helicopter. In [START_REF] Meister | Adaptive Path Planning for a VTOL-UAV[END_REF], adaptive path planning algorithms are developed for an optimized trajectory. Several trajectory optimization algorithms are presented in [START_REF] Kuwata | Trajectory Planning for Unmanned Vehicles Using Robust Receding Horizon Control[END_REF] for a team of cooperating unmanned vehicles. Flatness has been also employed for the trajectory planning of the quadrotor helicopter: a method is presented in [START_REF] Bouktir | Trajectory Planning for a Quadrotor Helicopter[END_REF] to generate time-optimal trajectories for the quadrotor system whereas Cowling et al. [START_REF] Cowling | Optimal Trajectory Planning and LQR Control for a Quadrotor UAV[END_REF] use the differential flatness to pose the trajectory planning as a constrained optimization problem in the output space.

Posing the trajectory planning as an optimization problem has one major limitation. The on- 1 Post-doctoral Fellow, Department of Mechanical and Industrial Engineering, 1515 St. Catherine West. 2 Associate Professor, Department of Mechanical and Industrial Engineering, 1515 St. Catherine West. AIAA Senior Member. Corresponding author. 3 Defense Scientist, Department of National Defense, 2459 Boul. Pie XI North. 4 Professor, Centre de Recherche en Automatique de Nancy (CRAN) -CNRS UMR 7039, BP 70239.

1 board implementation of the trajectory planning approach requires a microcontroller that is capable of solving a nonlinear optimization problem in real time. Moreover, most of the exciting works on trajectory planning are only illustrated in the simulation framework. The main contribution of this work with respect to existing works is in: a) proposing a simple and effective trajectory planning/re-planning method that can be solved in a straightforward manner without the need to solve a complex optimization problem in real time; and b) successfully flight-testing the approach through an experimental application to a cutting-edge quadrotor helicopter UAV.

The Note is organized as follows. Section II presents the mathematical model of the quadrotor helicopter. Section III presents the trajectory planning/re-planning method where two types of constraints are investigated: flight envelope and actuator limitations. The flight envelope is represented by the pitch and roll angles and the actuator limitations are the maximal thrusts that can be generated by the rotors. To initialize trajectory re-planning in the presence of actuator faults, fault detection and diagnosis problem is also considered. This is achieved by implementing on-board an unscented Kalman filter-based fault detection and diagnosis scheme. Experimental results are given in Section IV followed by conclusions in Section V.

II. Dynamics of the Quadrotor UAV System

The quadrotor UAV available at the Networked Autonomous Vehicles (NAV) Lab in the Department of Mechanical and Industrial Engineering of Concordia University is the Qball-X4 testbed well as its nonlinear model can be found in [START_REF] Amoozgar | Fault-Tolerant Fuzzy Gain Scheduled PID for a Quadrotor Helicopter Testbed in the Presence of Actuator Faults[END_REF]. In this paper, the following simplified model will be employed for trajectory planning as well as for the design of the baseline LQR controller:

ẍ = θg; J 1 θ = u θ ÿ = -φg; J 2 φ = u φ (1) z = u z /m -g; J 3 ψ = u ψ
This simplified model is obtained by assuming hovering conditions (u z ≈ mg in the x and y directions) with no yawing (ψ = 0) and small roll and pitch angles. x, y and z are the coordinates of the quadrotor UAV center of mass in the earth-fixed frame. θ, φ and ψ are the pitch, roll and yaw Euler angles respectively. m is the mass, g is the gravitational acceleration and J i (i = 1, 2, 3) are the moments of inertia along y, x and z directions respectively. u z is the total lift generated by the four propellers and applied to the quadrotor UAV in the z-direction (body-fixed frame). u θ , u φ and u ψ are respectively the applied torques in θ, φ and ψ directions. The Qball-X4 is equipped with four outrunner brushless motors driven by pulse-width modulated (PWM) inputs. By assuming that the motor dynamics is a simple gain K, it is possible to express the relation between the lift/torques and the PWM inputs u i as follows:

u z = K(u 1 + u 2 + u 3 + u 4 ) u θ = KL(u 1 -u 2 ) u φ = KL(u 3 -u 4 ) u ψ = KK ψ (u 1 + u 2 -u 3 -u 4 ) (2)
where K and K ψ are positive constants and L is the distance from the center of mass to each motor.

III. Flatness-based Trajectory Planning/Re-planning

The objective of the trajectory planning/re-planning is to find the profile of the path that a quadrotor UAV will follow when moving from an initial to a final position while respecting system constraints in nominal and fault conditions. Differential flatness is the first step towards this objective as explained hereafter.

A. Differential Flatness of the Quadrotor Helicopter

Flatness can be defined as follows. A general nonlinear dynamic system ẋ = f (x, u), y = h(x) with x ∈ n and u ∈ m , is flat if and only if there exist variables F ∈ m called the flat outputs such that: x = Ξ 1 (F, Ḟ , ..., F (n-1) ), y = Ξ 2 (F, Ḟ , ..., F (n-1) ) and u = Ξ 3 (F, Ḟ , ..., F (n) ) [START_REF] Fliess | Flatness and Defect of Nonlinear Systems: Introductory Theory and Example[END_REF], [START_REF] Faiz | Trajectory Planning of Differentially Flat Systems with Dynamics and Inequalities[END_REF].

Ξ 1 , Ξ 2 and Ξ 3 are three smooth mappings and F (i) is the i th derivative of F . The parameterization of the control inputs u in function of the flat outputs F plays a key role in the trajectory planning problem: the nominal control inputs to be applied during a mission can be expressed in function of the desired trajectories. This allows to tune the profile of the trajectories to keep the applied control inputs below the actuator limits. The same idea applies for the roll and pitch angles that the system will achieve during the mission. For the quadrotor UAV simplified model given in (1), the system is flat with flat outputs F 1 = z, F 2 = x, F 3 = y and F 4 = ψ. In addition to x, y, z and ψ, the parameterization of θ and φ in function of the flat outputs is:

θ = F2 g ; φ = - F3 g (3)
The parameterization of the control inputs in function of the flat outputs is:

u z = m( F1 + g); u θ = J 1 F (4) 2 g ; u φ = -J 2 F (4) 3 g ; u ψ = J 3 F4 (4)
Let F * i be the reference trajectory for the flat output F i with i = 1, ..., 4. An important consequence of the differential parameterization ( 4) is that if the system is forced to follow the reference trajectories, then the nominal control inputs to be applied along the trajectories are:

u * z = m( F * 1 + g); u * θ = J 1 F (4) * 2 g ; u * φ = -J 2 F (4) * 3 g ; u * ψ = J 3 F * 4 (5)

B. Trajectory Parametrization

Before considering the trajectory planning problem, the reference trajectory F * i will be first time-parameterized. The total time of the mission is then tuned so that the system constraints are respected. In [START_REF] Ren | Virtual Structure Based Spacecraft Formation Control with Formation Feedback[END_REF], the authors use second-order systems for the reference trajectories in the framework of spacecrafts in formation flight. They also approximate the natural frequency as function of the rise time to define the transient specifications for the formation. However, they state that it is hard to choose a good rise time beforehand to achieve a good performance in the system. On one hand, if the rise time is small then the formation moves too fast and the formation members could fall far behind their desired positions. On the other hand, if the rise time is big then the formation moves too slowly and the mission cannot be achieved within short time. Inspired by that work, this paper proposes to solve the above-mentioned problem by using differential flatness. Let us first define the reference trajectories F * i (i = 1, ..., 4) as second-order systems:

F * i (t) = -2ξω n Ḟ * i (t) -ω 2 n F * i (t) + ω 2 n r i (t) (6)
where ξ is the damping ratio and ω n is the natural frequency. r i (t) (i = 1, ..., 4) are the references along z, x, y and ψ directions respectively and are defined as step functions of amplitude R i .

Consider a critically damped system by setting ξ = 1 and assuming zero initial conditions, it turns out that the solutions of the reference models given in (6) are:

F * i (t) = 1 -(1 + ω n t) e -ωnt R i ; i = 1, ..., 4 (7)
This work considers the actuator faults that may take place during the mission. When faults occur it is necessary to re-plan the nominal trajectories. This is because forcing the system to follow the same pre-fault trajectory may lead to saturation or even instability if the new limitations of the actuators are not taken into consideration. Eq. ( 7) is only valid for zero initial conditions. However, when faults occur and re-planning is needed, the initial conditions at the re-planning instant depend on the quadrotor status and they are not necessarily zero. With a critically damped system (ξ = 1)

and non-zero initial conditions, the solutions of the reference models given in (6) are:

F * i (t) = R i -(1 + ω n t) e -ωnt ∆ i + Ḟ * i (t rep ) te -ωnt ; i = 1, ..., 4 (8)
where

∆ i = R i -F * i (t rep )
is the remaining distance to the desired position and

F * i (t rep ) and Ḟ * i (t rep )
are the position and the velocity of the system at the instant of re-planning t rep . It can be shown that the presence of the term Ḟ * i (t rep ) te -ωnt in the reference trajectory (due to the non-zero initial velocity) makes it impossible to derive a solution for the trajectory re-planning in a straightforward manner. Nevertheless, it is always possible to find the solution to the trajectory re-planning problem by posing it as an optimization problem where the function to minimize is the mission time and the constraints are those of the flight envelope or the actuator limits. However, the objective of this work is to minimize as much as possible the calculation requirements for the quadrotor on-board microcomputer. To this end, it is proposed to proceed as following. Once a fault is detected by a Fault Detection and Diagnosis (FDD) scheme, the trajectory is re-planned with a sufficiently large mission time. This big mission time will result in a very slow reference trajectory serving in two directions. First, it will give sufficient time to diagnose (isolate and identify) the fault before taking any action since the system will almost stay in hovering position. Second, the system velocity will drop to zero and therefore, the last term of ( 8) can be set to zero. This allows to easily re-plan the trajectory as will be shown in the sequel.

C. Fault Detection and Diagnosis: Unscented Kalman Filter for Parameter Estimation

Actuator faults can be modeled as u f i = (1 -w i )u i for i = 1, ..., 4. w i represents the loss of effectiveness in the i th rotor. w i = 0 denotes a healthy rotor, w i = 1 denotes a complete loss of the i th rotor and 0 < w i < 1 represents a partial loss of control effectiveness. As will be shown in the subsequent sections, trajectory re-planning in fault cases requires the knowledge of the fault amplitude w i . Thus, an FDD module is needed to detect, isolate and identify the fault. This is achieved by employing a parameter estimation-based FDD where the unscented Kalman filter [START_REF] Julier | A New Extension of the Kalman Filter to Nonlinear Systems[END_REF] is used to estimate the parameters w i . It should be noted that for FDD purpose, the nonlinear model [START_REF] Amoozgar | Fault-Tolerant Fuzzy Gain Scheduled PID for a Quadrotor Helicopter Testbed in the Presence of Actuator Faults[END_REF] is employed. In the fault-free case, the estimates of w i (i = 1, ..., 4) are close to 0 and a deviation from 0 will be an indication for a fault occurrence. In the subsequent sections, the trajectory planning/re-planning approach will be explained for two cases: flight envelope constraints and actuator constraints.

D. Trajectory Planning/Re-planning Under Flight Envelope Constraints

A trajectory with a small travel time involves aggressive maneuvers with big roll and pitch angles. A controller based on the nonlinear model would be able to guarantee the system stability.

However, a controller based on a linearized model (such as (1)) may not be able to keep system's stability when the quadrotor goes outside the linear operating zone. To avoid this latter situation, it is possible to impose bounds on the maximal pitch and roll angles that are attained during the system's travel from the initial to the final position. According to (3), the pitch angle of the system when it is forced to follow the reference trajectory is given by:

θ * = F * 2 g (9)
Using ( 9) and the double time derivative of (7), one can find the nominal pitch angle θ * in function of the time t, the natural frequency ω n , and the amplitude R 2 as following:

θ * = - ω 2 n g (tω n -1) e -ωnt R 2 (10)
For a critically damped second-order system (ξ = 1) and for a settling time specification of 2% of the steady state value, it is possible to approximate the natural frequency as:

ω n ≈ 5.83 t s (11)
where t s is the settling time of the reference model. Eq. ( 10) reads then:

θ * = - 5.83 2 gt 2 s 5.83t t s -1 R 2 e - 5.83t t s (12)
It is necessary at this stage to find a relation between the settling time t s and the total time of the mission t t . For a critically damped system (ξ = 1), it is reasonable to assume that the settling time approximately equals the total time of the mission (t s ≈ t t ). The objective turns out to determine t s such that -ρ θ θ max ≤ θ * ≤ ρ θ θ max where θ max is the maximal pitch angle that the system is allowed to make during its travel. Since ( 9) is derived using the simplified model ( 1), the parameter 0 < ρ θ < 1 is introduced to create a safety margin to deal with model uncertainties. One way to calculate t s so that -ρ θ θ max ≤ θ * ≤ ρ θ θ max is to determine where the maximal pitch is taking place and then tune t s such that the angle constraints are respected. To determine where the maximal pitch is taking place, let us calculate first the time derivative of the pitch angle with respect to time:

dθ * dt = 5.83 3 gt 3 s 5.83t t s -2 R 2 e - 5.83t t s (13)
By setting [START_REF] Julier | A New Extension of the Kalman Filter to Nonlinear Systems[END_REF] to zero, one finds that the maximal pitch angle is taking place at t = 2t s /5.83.

Plugging this value of t in [START_REF] Faiz | Trajectory Planning of Differentially Flat Systems with Dynamics and Inequalities[END_REF] gives the extrema:

θ * Ext = - 5.83 2 gt 2 s e 2 R 2 (14)
where the extrema collectively denote the minima and the maxima of a function. It is also necessary to check the value of the pitch angle at the beginning and at the end of the mission, i.e. for t = 0 and t = t s . These are given by: which respectively corresponds to t = 0, t = 2t s /5.83, and t = t s . The solution to consider is the maximal one among the three possible solutions. It is easy to see however that the solution corresponding to t = 0 is the maximal one and hence it is the one to be considered. Similarly to the study carried out for the pitch angle θ, it is possible to derive the solution for the roll angle φ which is associated to the reference r 3 (t) along the y-direction as following:

θ * (0) = 5.
t s ≥ 5.83 2 |R 3 | gρ φ φ max (17)
In the fault case, the trajectory re-planning can be achieved by replacing R 2 and R 3 in ( 16) and ( 17) with the remaining distances to travel

∆ 2 = R 2 -F 2 (t rep ) and ∆ 3 = R 3 -F 3 (t rep ) respectively.
Remark 1 The solution obtained under flight envelope constraints does not explicitly consider the actuator limits. Therefore, it is necessary to calculate t s by considering actuator constraints as explained in the subsequent section. If the solution obtained under actuator constraints is bigger than the one obtained under flight envelope constraints then the former must be considered since the actuator constraints are more restrictive than those of the flight envelope.

E. Trajectory Planning/Re-planning Under Actuator Constraints

If aggressive maneuvers are allowed, the constraints on the pitch and roll angles can be ignored.

However, it is necessary to consider actuator constraints. When actuators hit their limits and cannot deliver the actuation inputs desired by the controller this results in degraded performance or can even lead the closed-loop system to an unstable behavior [START_REF] Henrion | Control of Linear Systems Subject to Input Constraints: a Polynomial Approach[END_REF]. The control inputs to be applied along the reference trajectories are given in [START_REF] Bouktir | Trajectory Planning for a Quadrotor Helicopter[END_REF]. In this work, the system is assumed to not changing altitude during the mission and thus F * 1 is constant and F * 1 = 0. It is also assumed that the system is not yawing, thus F * 4 = 0 and F * 4 = 0. By using (5), the control inputs are then:

u * z = mg; u * θ = J 1 F (4) * 2 g ; u * φ = -J 2 F (4) * 3 g ; u * ψ = 0 (18)
The relation between the lift/moments and the PWM inputs is given by [START_REF] Lai | Time-optimal Control of a Hovering Quad-rotor Helicopter[END_REF]. Taking the inverse of (2) and using (18) gives the PWM inputs to be applied along the reference trajectories. For instance, the first PWM input u * 1 (t) is:

u * 1 (t) = 1 4K mg + J 1 2gKL F (4) * 2 (t) (19)
Since F * 2 (t) is given in (7) it is possible to calculate the fourth derivative of F * 2 (t) with respect to time. Then, approximating the natural frequency as ω n ≈ 5.83/t s gives the nominal PWM inputs in function of t and t s as:

u * 1 (t) = 1 4K mg - J 1 2gKL 5.83 4 t 4 s 5.83t t s -3 R 2 e - 5.83t t s (20)
The objective is to determine the settling time t s so that the nominal PWM inputs are less than the maximal allowable input u max , i.e. u * i (t) ≤ u max ∀t ∈ [0, t s ] for i = 1, ..., 4. Similarly as before, the extrema of u * i (t) must be determined as well as its terminal values at t = 0 and t = t s .

For u * 1 (t), the extrema is taking place at t = 4t s /5.83 and its value is:

u * 1Ext = 1 4K mg - J 1 2gKL 5.83 4 e 4 t 4 s R 2 (21)
The terminal values at t = 0 and t = t s are:

u * 1 (0) = 1 4K mg + 3J 1 2gKL 5.83 4 t 4 s R 2 and u * 1 (t s ) = 1 4K mg - J 1 2gKL 2.83 × 5.83 4 t 4 s e 5.83 R 2 (22)
From ( 21) and ( 22), one can determine t s for which u * 1 (t) ≤ ρ u u max . That is:

t s ≥ max      - 5.83 4 J 1 R 2 2gKLe 4 ūmax 1 4 ; 3 × 5.83 4 J 1 R 2 2gKLū max 1 4 ; - 2.83 × 5.83 4 J 1 R 2 2gKLe 5.83 ūmax 1 4      (23)
where ūmax = ρ u u max -mg/4K and ρ u works similarly to ρ θ . Similar development can be carried out for the three remaining PWM inputs u * 2 (t), u * 3 (t) and u * 4 (t). After ignoring the complex values of t s , four solutions are obtained and the maximal one is the solution to be considered.

Trajectory re-planning consists in solving (23) by replacing R 2 with the remaining distance to travel ∆ 2 = R 2 -F 2 (t rep ) and ūmax with the new actuator constraint ū1 max = (1 -ŵ1 )ρ u u maxmg/4K. ŵ1 is the estimate of effectiveness loss in the first actuator provided by the unscented Kalman filter of Section III C. Based on the obtained solution t s , the damaged system must decide whether to continue the mission, return to the base or land safely while taking into consideration the remaining fuel or battery power and/or the attainable speed.

IV. Experimental Results

The trajectory planning/re-planning approach is applied to the quadrotor UAV testbed. The baseline controller for the Qball-X4 is an LQR controller which is designed based on the simplified model [START_REF] Tsourdos | Cooperative Path Planning of Unmanned Aerial Vehicles[END_REF]. The fault-free trajectory planning is illustrated in Figures 1(a) and 1(b) with R 2 = -1.5 m and R 3 = 0. The figures show the pitch angle and position tracking as well as the system's height which is set to 0.3 m and the PWM inputs to the four rotors. The reference trajectory along the x-direction is given by (7) whereas the desired pitch angle is generated by the LQR controller and its evolution is predicted in (9) thanks to flatness. The trajectory re-planning is illustrated in Figure 2(a) where a 40% loss of control effectiveness is taking place in the second motor (w 2 = 0.4) at time instant t = 31.5 s. The fault makes the system to loose its height and to fall behind the reference trajectory. This sudden jump in the tracking error along x and z directions indicates an eventual occurrence of a fault. Once the fault is detected, the trajectory is temporarily re-planned at about t = 32 s with a t s = 900 s. shows that the fault is successfully identified where 1 -w 2 = 0.6 and 1 -w 1 = 1 -w 3 = 1 -w 4 = 1.

The temporary trajectory lasts (for illustration purpose) for about 30 s where the system is held in place before re-planning the trajectory with t s = 7 s. The re-planned trajectory is slower than the initial pre-fault trajectory to take into consideration the new actuator limits after fault occurrence. This paper proposes a trajectory planning/re-planning approach where reference trajectories are defined as second-order systems. The main advantage of the proposed approach is that solutions are obtained in a straightforward manner by solving simple equations in the fault-free and fault conditions. Therefore, it has very few calculation requirements and it is very suitable for systems with limited calculation capabilities such as small-scale UAVs with single-board microcomputers.

  [7]. It is developed by Quanser Inc. under the financial support of an NSERC (Natural Sciences and Engineering Research Council of Canada) Strategic Project Grant lead by Concordia University with Quanser Inc. as one of three industrial partners. A detailed description of the Qball-X4 as

Figure 1 (Fig. 1

 11 Fig. 1 Trajectory planning under flight envelope and actuator constraints

  Figure 2(b)

Fig. 2

 2 Fig. 2 Trajectory re-planning and fault identification

Acknowledgments

This work is supported by the Natural Sciences and Engineering Research Council of Canada (NSERC) Postdoctoral Fellowship (PDF) program and the NSERC Strategic Project Grant (STPGP 350889-07). Support from Quanser Inc. and colleagues from Quanser Inc. for the development of the Qball-X4 UAV testbed is also highly appreciated.