T Leroy

J Chauvin

N Petit

G Corde

MOTION PLANNING CONTROL OF THE AIRPATH OF AN S.I. ENGINE WITH VALVE TIMING ACTUATORS

Keywords: Airpath control, SI Engine, VVT, Motion planning

We address the control of the airpath of a turbocharged S.I. engine equipped with Variable Valve Timing actuators (VVT). Compared to standard configurations, the engine does not possess any external EGR (Exhaust Gas Recirculation) loop. Rather, VVT are used to produce internal EGR, providing similar beneficial effects in terms of emissions reduction. The airpath dynamics takes the form of a single mono-dimensional air balance in the intake manifold. In this equation, the VVT act as a disturbance by impacting on the air mass flow through the inlet valves. This impact can be estimated from real-time measurements. We use this information in a motion planning based control strategy by, successively, turning the driver's torque demand into a trajectory generation problem for the air mass contained in the intake manifold, and then deriving an intake manifold pressure trajectory. Supportive simulation results show the relevance of this approach and suggest ways of further improvements.

INTRODUCTION

Lately, downsizing (reduction of the engine size) has appeared as a major way of reducing fuel consumption of SI engines. Downsizing lowers pumping losses under partial load and, thus, saves energy (see [START_REF] Lecointe | Downsizing a gasoline engine using turbocharging with direct injection[END_REF]). To compensate the induced power losses, turbochargers are used. For such engines, Variable Valve Timing (VVT) devices are used to improve combustion efficiency and reduce NO x emissions.

In this paper, we focus on two control objectives. What the driver primarily expects from the engine control system is to respond to its torque demand. Since SI engine dynamics are mostly represented by the airpath, controlling the torque is equivalent to controlling the air mass aspirated in the cylinder. Besides, it is also desired to reduce engine pollutant emissions. Three-way catalysts (TWC) are used to reduce CO, HC and NO x created by the combustion. However, to efficiently convert these species to water and carbon dioxide, TWC need that the air/fuel ratio resides in a very narrow region around the stoichiometric value 1 (see Figure 1) 1 . In practice, the air mass in the cylinder, which is not measured, can be estimated to adjust the injected fuel quantity so that the air/fuel ratio equals the stoichiometric ideal value. 1 The air/fuel ratio under consideration here is the ratio of trapped air mass and fuel mass in the cylinder when the valves are closed, it is defined by λ On classic SI engines, a variable named the "volumetric efficiency" is used to estimate the air mass in the cylinder through balance equations.

Usually, its values are calibrated under steady state running engine conditions and ambient temperature and pressure conditions (see [START_REF] Chevalier | On the validity of mean value engine models during transient operation[END_REF], [START_REF] Andersson | Air Charge Estimation in Turbocharged Spark Ignition Engines[END_REF], (Chevalier et al., 2000a), [START_REF] Eriksson | Modeling of a turbocharged SI engine[END_REF] for example).

The volumetric efficiency is modelled as a nonlinear function of the engine speed, which is usually easily available from measurements, and the intake pressure.

This calibration task is much more involved in the case of our downsized engines. Indeed, the position of the VVT actuator affects the exhaust manifold pressure and the residual gas mass trapped inside the cylinder. This has an influence on the volumetric efficiency which, as has been discussed, is the core variable when computing the fuel injection feedforward term. Figure 2 reports experiments stressing the strong influence of VVT on the volumetric efficiency. Experimentally, it is difficult to obtain reliable look-up tables of the volumetric efficiency for every engine speed, intake pressure and VVT po-sition. Yet, observers can be used to solve an input estimation problem (see [START_REF] Stotsky | Application of input estimation techniques to charge estimation and control in automotive engines[END_REF] for a global overview) and to estimate the trapped air mass in the cylinder. Volumetric efficiency modelling and determination by the measurement of the intake manifold air flow, the manifold air pressure and temperature has been adressed previously in the literature (see [START_REF] Andersson | Air-tocylinder observer on a turbocharged SI-engine with wastegate[END_REF], (Chevalier et al., 2000a), [START_REF] Chauvin | Air path estimation on Diesel HCCI engine[END_REF], [START_REF] Grizzle | Improved cylinder air charge estimation for transient air fuel ratio control[END_REF]). These approaches lead to an observer (usually an extended Kalman filter) relying on a mean-value model of the intake flow (at the TDC time-scale). This estimation of the trapped air mass is not suffisant to get the expected amount of fuel. Actually, because of computational and injection duration, fuel control signal is only applied about 2 TDC later (see (Chevalier et al., 2000a)). To compensate the delay, port air mass flow must be predicted. The prediction is then given as a feedforward to compute a fuel mass set point (M F F f uel = 1 14.6 M pred air ). This strategy is the most used for air/fuel ratio control in SI engine without VVT devices. This is more difficult in VVT SI engines. Indeed, there is no uniqueness of the relationship between intake pressure and air mass trapped (see Figure 2), so a simple prediction of the intake pressure from a prediction of the throttle position is no more possible. That is why we propose an alternative approach. We simply use a planned trajectory of the air mass as a feedforward term in the fuel injection (M F F f uel = 1 14.6 M mp air ). A simple nonlinear observer reconstructs the volumetric efficiency. It is only used to convert the air mass control problem into an intake pressure control objective. Convergence of the observer guarantees that, asymptotically, the air mass set point is reached. The main advantage of substituting the estimated value of the air mass with a (varying) set point is that it allows an easy synchronization of the air and fuel controllers. We send the future reference signal (two TDC ahead, see Figure 4) to the fuel injection system and the current reference to the air path controller. This synchronization improves the matching of the air mass and the fuel mass, yielding a better tracking of the air/fuel ratio reference. This strategy appears efficient as is demonstrated by simulation results.

Pman (P a) ηvol (-) V V Tin = 40, V V Texh = 0 V V Tin = 20, V V Texh = 0 V V Tin = 0, V V Texh = 0 V V Tin = 0, V V Texh = 20 V V Tin = 0, V V Texh = 40 0 0.5 1 1.5 2
The paper is organised as follows. In Section 2, we present the reference model of the intake manifold. We consider a mass balance equation and a model of the aspirated flow. In Section 3, we present our control structure. In Section 4, we propose an observer that estimates a correcting term in the volumetric efficiency appearing in the aspirated flow. It allows to take into account variables that are not considered in the data fit (VVT, exhaust pressure, exhaust temperature, ...). In Section 5, we detail the motion planning technique and the feedforward control strategy. Section 6 briefly reports the model used for simulation and presents obtained results. We finally conclude and give future directions in Section 7.

AIR PATH MODELLING

Intake

Mass balance in the intake manifold

Notations are given in Table 1. We consider the intake manifold as a fixed volume for which the thermodynamic states (pressures, temperatures, composition) are assumed homogeneous. We neglect variations of temperature in the volume (as in [START_REF] Guzzella | Introduction to Modeling and Control of Internal Combustion Engine Systems[END_REF], [START_REF] Andersson | Air-tocylinder observer on a turbocharged SI-engine with wastegate[END_REF] for example), i.e. Ṫman = 0. The mass balance in the intake manifold gives

Ṗman = RT man V man (D air -D asp ) (1)
where P man and T man are measured by sensors located in the intake manifold, and D air can be obtained by an air flow sensor in the intake of the engine as represented in Figure 3. D asp is the air flow aspirated into the cylinders. Subsequently we model the air mass flow through inlet valves [START_REF] Heywood | Internal Combustion Engine Fundamentals[END_REF],

D asp = η vol (P man , Φ) V cyl RT man N e 30 P man (2)
where η vol (P man , Φ) is the volumetric efficiency. It depends on the intake manifold pressure, P man , and on other parameters such as the engine speed, N e , and the VVT positions, i.e. Φ (N e , V V T in , V V T exh ).

Reference model

Let x P man and y P man be the state and the measurement of the system, respectively. We note Pressure in the intake manifold

P a R Ideal gas constant J.kg -1 .K -1 S th
Openning throttle section m 2 Tman Temperature in the intake manifold

K V cyl Cylinder volume m 3 Vman Intake manifold volume m 3 V V T in Intake VVT position • CA V V T exh Exhaust VVT position • CA η vol Volumetric efficiency - λ Air/fuel ratio - Φ Φ = (Ne, V V T in , V V T exh ) -
the throttle opening section u S th . The mass air flow through the throttle now writes D air uf (x) where f (x) is the mass flow rate [START_REF] Heywood | Internal Combustion Engine Fundamentals[END_REF])

f (x) = κP int                    x P int 1 γ 2γ γ -1 1 - x P int γ-1 γ if x Pint > 0.528 γ 2 γ + 1 γ+1 γ-1 else (3) where κ = 1 √
RTman , P int is the upstream pressure of the throttle, and γ is the specific heats ratio in the intake manifold.

The air mass trapped in the cylinder is given by z M air . A simple integration gives z = D asp 30 Ne . With the notations reported in Table 2, the reference system is Air mass trajectory generation (block A) The torque demand is directly specified by the driver through the acceleration pedal. This information is converted into a trapped air mass set point, z sp , thanks to a two dimensional look-up table (experimentally determined on a test bench) having torque and engine speed as inputs. This information is then filtered to create a feasible trajectory of trapped air mass. We do not detail this filter here. Its structure is discussed in (Chauvin et al., 2006a). We note z mp the obtained air mass trajectory.

     ẋ = α int (uf (x) -η vol (x, Φ)β int x) y = x z = η vol (x, Φ)β int 30 N e x (4)
Air observer (block B) This block contains an observer that estimates the air mass trapped into the cylinders ẑ from engine speed, N e , intake flow, D air , intake pressure, y, and intake temperature, T man , measurements. It uses a model of the intake manifold and takes VVT variations into account. This is presented in Section 4. This observer also estimates the fluctuation of the volumetric efficiency.

Pressure trajectory calculation (block C) Here, the desired air mass trajectory z mp is turned into a pressure trajectory x mp from the dynamics equation of the intake manifold. This block uses the estimation ẑ and the reconstructed volumetric efficiency variation δη coming from block B. In this paper, we do not tackle turbocharged conditions (intake pressure higher than 1 bar) which will be discussed in a next paper. However, we take into account the upstream throttle pressure in simulation, P int , which is varying because of compressor rotation. This block is detailed in Section 5.

Dynamics model inversion (block D) This block

generates an open-loop control law u F F to keep the aspirated air mass equal to the desired air mass, z mp . This feedforward term is determined by inverting the intake manifold dynamics. This computation is presented in Section 5.

Stoichiometric air/fuel ratio (block E) This block computes a fuel mass set point to keep the air/fuel ratio to the stoichiometric value.It is simply determined by applying the following stoichiometric relationship between air and fuel

M F F f uel = 1 14.6 z mp 2 T DC ahead

AIR MASS ESTIMATION

We now detail the design of the air observer (block B). Previously η vol (x, Φ) was defined as the cylinder filling volumetric efficiency. For sake of simplicity, it is usually assumed that η vol depends only on engine speed and intake pressure. To model the influence of the VVT, we consider a correcting term δη. Then, equation (2) becomes

D asp = (η vol (x, N e ) + δη)β int x (5)
This approach follows several ideas proposed in [START_REF] Andersson | Air-tocylinder observer on a turbocharged SI-engine with wastegate[END_REF], [START_REF] Stotsky | Application of input estimation techniques to charge estimation and control in automotive engines[END_REF], [START_REF] Chauvin | Air path estimation on Diesel HCCI engine[END_REF]. The correcting term δη is used to estimate the trapped air mass z thanks to the relationship z = D asp 30

Ne . We propose to estimate δη with an observer using intake pressure measurement.

Let x 1

x and x 2 δη be the states of the observer. The measurement is noted y. Considering equations ( 1) and ( 5), the reference model of the intake manifold writes

   ẋ1 = α int (D air -(η vol (x 1 , N e ) + x 2 )β int x 1 ) ẋ2 = 0 y = x 1 (6) We propose the following observer    ẋ1 = α int (D air -(η vol (y, N e ) + x2 )β int x1 ) -L 1 (x 1 -y) ẋ2 = L 2 (x 1 -y) (7) where (L 1 , L 2 ) ∈ (R + \ {0}) 2 . The state-error is x = x -x. It satisfies    ẋ1 = α int β int (-(η vol (y, N e ) + x 2 )x 1 -x1 x2 ) -L 1 x1 ẋ2 = L 2 x1
(8) There, we see the importance of using the air mass flow measurement D air instead of the throttle flow uf (x 1 ) in the observer structure. Thanks to this arbitrary choice in the observer design, the D air term is cancelled out in the error system (8) and the control variable u does not appear anymore. Observer design and controller design can be treated separately.

Gains are chosen as follows

L 1 = α int β int l 1 and L 2 = α int β int x1 l 2
where l 1 and l 2 are positive constants. With this choice, the error system writes

   ẋ1 = α int β int (-(η vol (y, N e ) + x 2 + l 1 )x 1 -x 1 x2 ) ẋ2 = α int β int x1 l 2 x1
Theoretical convergence proof is discussed in [START_REF] Leroy | Air path estimation for a turbocharged SI engine with variable valve timing[END_REF].

AIR MASS CONTROL

The air path control strategy we propose consists of managing the trapped air mass in the cylinders. The idea is not to improve the existing control in terms of speed or control effort but to master the air mass trajectory with the intention of improving air/fuel ratio control during torque transients. Trapped air mass depends on intake pressure and other parameters (see discussion in Section 2). In this paper, we do not consider the control of the VVT actuators but we which to counteract their influences on trapped air mass. The only actuator left is thus the throttle.

Motion planning

It is desired to accurately control the trapped air mass in the cylinders. A trajectory is generated from the air mass set point z sp . The filter we use is presented in (Chauvin et al., 2006a). This yields an air mass motion planning trajectory z mp . Yet, the balance equation ( 4) is expressed in terms of the intake pressure x. So, we need to cast the air mass motion planning trajectory into an intake pressure motion planning trajectory x mp . This is not quite simple because the relationship between these two variables is uncertain. Let us write

z = ϕ(x) = (η vol (x, N e ) + δη)β int 30 N e x ( 9
)
where δη is given by the air observer. Equation ( 9) can be inverted with respect to x, yielding

x mp = ϕ -1 (z mp ) (10)
For that purpose, a Newton algorithm with constant step can be used.

The generated pressure trajectory highly depends on the air observer results. We will see in Section 6 that transient performance losses may occur if the observer converges too slowly.

Dynamics model inversion

System (4) is invertible. Thus, an analytic expression of the input can be derived from the state variable x and its first derivative history. In fact,

u = 1 f (x) ẋ α int + ϕ(x) N e 30
In this last expression, α int and N e are given by sensor measurements. The unique open-loop control law u F F corresponding to any desired z mp trajectory is

u F F = 1 f (x mp ) ẋmp α int + z mp N e 30 ( 11
)
where x mp = ϕ -1 (z mp ).

SIMULATION RESULTS

Tests setup

The proposed motion planning control strategy has been tested in simulation. The simulator under consideration is a high frequency engine model developed on AMESim platform (Le [START_REF] Berr | Modelling of a turbocharged SI engine with variable camshaft timing for engine control purposes[END_REF].

Results

Figures 5 presents experimental open-loop results.

We test the control strategy on a torque trajectory at constant engine speed (1500 rpm). VVT actuators are not plotted in the figure, but they move according to the engine speed and the air mass trapped. We now detail these experiments.

z sp z mp ẑ z ref λ sp λ meas , M F F f uel = 1 14.6 z mp 2 T DC ahead λ meas , M F F f uel = 1 14.6 ẑ u F F u meas
x mp x meas A first order dynamics is given to the throttle actuator.

1) From t=12s to t=14s: we have here a classic tip-in. Torque starts at 20 Nm and reaches 100 Nm. Figure 5 (a) shows the trapped air mass set point z sp and the filtered z mp trajectory. The trajectory generation filter is tuned to be smooth enough to avoid generating a non feasible trajectory, and not too slow, to keep a certain reactivity in torque response. Air estimation ẑ given by the observer and AMESim reference trapped air mass z ref are also plotted in the figure. We can see that convergence of the observer through the reference is pretty slow. This is the real bottleneck of our strategy here. Relatively poor convergence has an impact that can be seen in Figure 5 (c). Actually, an overshoot is created on the x mp trajectory. This can be explained by looking at equation ( 9), x mp depending on δη. We come across this large overshoot in throttle command u F F in Figure 5 (d) and of course on reference air mass trapped in the cylinder z ref in Figure 5 (a). Figure 5 (b) gives the fuel/air ratio measurements in the case we give the estimated air mass ẑ and in the case we give the 2 TDC ahead generated trajectory x mp for the fuel set point computation. There we see the relevance of the presented method looking at the very low excursions of the air/fuel ratio.

2) From t=14s to t=16s: we have now a classic tip-out. Torque starts at 100 Nm and decreases to 20 Nm. Observation errors are almost the same as for tip-in transient except that the reference trapped air mass z ref is much better. Actually, after the transient, u F F solely depends on z mp and f (x mp ). Now considering chocked flow conditions (see equation ( 3)), u F F does not more depends on x mp . The overshoot is therefore avoided.

CONCLUSION

This paper presents a motion-planning based method of controlling the air path of a downsized SI engine equipped with VVT actuators.

Simulation results stress the efficiency of the approach. However, the relative slowness of the air observer generates over and undershoots on reference trapped air mass during torque transients. Improving convergence of the observer seems a natural direction for further investigations. When this problem is resolved, we plan to add a feedback term to the throttle control to cancel further modelling errors. Another PI controller could also be considered for the fuelpath control.

Fig. 1 .

 1 Fig. 1. Three-way catalytic converter (TWC) efficiency as a function of air/fuel ratio.

Fig. 2 .

 2 Fig. 2. Experimental volumetric efficiency, for different VVT positions, as a function of the intake manifold pressure at constant engine speed (2000 rpm).

Fig. 3 .

 3 Fig. 3. Engine scheme, VVT are used to produce internal EGR.

Fig. 4 .

 4 Fig. 4. Feedforward control strategy for the airpath system. The unique closed loop is through the observer. Its estimates are used in the derivation of the open-loop control law u F F . Block B uses equation (7), block C uses equation (10) and block D uses equation (11).

Fig. 5 .

 5 Fig. 5. Simulation results: torque transient from 20 to 100 Nm and then to 20 Nm at 1500 rpm. (a) Solid: air mass set point z sp , Dashed: generated trajectory z mp , Dotted: observed trapped air mass ẑ, Dashdot: reference AMESim trapped air mass z ref . (b) Solid: stoichiometric fuel/air ratio, Dashed: measured fuel/air ratio when z mp is the feedforward of the fuel injection set point computation, Dotted: measured fuel/air ratio when ẑ is the feedforward of the fuel injection set point computation. (c) Solid : Generated pressure trajectory x mp , Dashed : Measured intake pressure x meas . (d) Solid: opening throttle area feedforward command u F F , Dashed: measured openning throttle area u meas .A first order dynamics is given to the throttle actuator.

Table 1 .

 1 Nomenclature

	α int	RTman Vman and β int	V cyl RTman	Ne 30 . The control is

Table 2 .

 2 Variables description

	Variable	Symbol	Variable	Expression
	u x	S th Pman	α int β int	RTman Vman V cyl Ne RTman 30
	y	Pman
	z	M air
	3. CONTROL STRATEGY
	Figure 4 details the control strategy we propose.
	We now present its five constitutive blocks.

TABLE LOOK

 LOOK

	torque Driver's	-UP
	request
		z mp 2 T DC ahead


Jean-Stéphane Varré

Bertil Schmidt

Stéphane Janot

Mathieu Giraud

Manycore high-performance computing in bioinformatics

Mining the increasing amount of genomic data requires having very efficient tools. Increasing the efficiency can be obtained with better algorithms, but one could also take advantage of the hardware itself to reduce the application runtimes.

Since a few years, issues with heat dissipation prevent the processors from having higher frequencies. One of the answers to maintain Moore's Law is parallel processing. Grid environments provide tools for effective implementation of coarse grain parallelization. Recently, another kind of hardware has attracted interest: multicore processors.

Graphic processing units (GPUs) are a first step towards massively multicore processors. They allow everyone to have some teraflops of cheap computing power in its personal computer.

The CUDA library (released in 2007) and the new standard OpenCL (specified in 2008) make programming of such devices very convenient. OpenCL is likely to gain a wide industrial support and to become a standard of choice for parallel programming. In all cases, the best speedups are obtained when combining precise algorithmic studies with a knowledge of the computing architectures. This is especially true with the memory hierarchy: the algorithms have to find a good balance between using large (and slow) global memories and some fast (but small) local memories.

In this chapter, we will show how those manycore devices enable more efficient bioinformatics applications. We will first give some insights into architectures and parallelism. Then we will describe recent implementations specifically designed for manycore architectures, including algorithms on sequence alignment and RNA structure prediction. We will conclude with some thoughts about the dissemination of those algorithms and implementations: are they today available on the bookshelf for everyone?

Introduction

Most of the algorithms presented in this book require powerful computers to run. These problems are intrinsically complex to solve, and, above all, the amount of input data follows an exponential curve. With next-generation sequencers (NGS), the data produced every day is growing faster than before. The exponential grows faster, much faster than the growth of computational power.

The challenge for computer scientists and bioinformaticians is to think about new methods and technologies able to deal with this incredible amount of data. A first solution is to design better algorithms with more elaborate data structures. However, even if a good algorithm is known, the growth of data still impose to have better supports of execution. To increase the computational power, one may use several computers or processors in a grid. Another possibily is to use multicore processors. Now manufacturers manage to put several processors in one: this is what is called dual-core, quad-core, or, more generally, multicore processors.

In this chapter, we will see how these recent advances in massively multicore processors, also called manycore processors, help to treat the growing flow of bioinformatics data. Manycore processors like Graphics Processing Units (GPUs) offer a huge density of computational units, for an extremely cheap price. They are more and more often directly included in personal computers, making their power available without extra cost. Thus, developments of methods using them have an impact stronger than ever. But the drawback of this technology is that the design of algorithms must take into account the constraints of these architectures. In fact, completely new distributed algorithmics are required, giving new challenges for bioinformaticians. But it is worth the effort: results published in 2008 or 2009 achieve speedups up to 100 x compared to serialized one-core algorithms. . . on commodity GPUs costing less that $500 ! In the following, we briefly describe the evolution of processors, that leads to manycore processors. The evolution of power consumption and heat dissipation led to the so-called power wall, that limits the rise of frequencies and commands us to use more and more cores in parallel. In Methods, we present manycore processors and their programming. In Results, we detail algorithms in different fields of bioinformatics, explaining where gains can or cannot be obtained. Finally, we discuss whether these solutions can now be used by everyone, on real applications.

A small history of processors

A processor, or Central Processing Unit (CPU), is roughly made of three units: the memory unit that stores input data, intermediate results and output data, the instruction unit that decodes instructions and the processing unit, also called arithmetic logic unit (ALU), that actually realizes the computations (Figure 1a) at regular clock cycles. The following paragraphs explore two dimensions contributing to the "computational power" of processors: their complexity, driven by the number of transitors, and the frequency of their clock cycles.

Moore's Law. The number of transistors used to build a chip has dramatically increased, from 2300 for the Intel 4004, released in 1971, to several billions today. Moore's law, formulated in [Moore, 1965] and precised in [Moore, 1975], states that this number doubles every two years. This law has been almost exactly verified since 1960 (Figure 2). More than a natural observation, this is a self-fulfilling prophecy that drives the semiconductor industry.

With more transistors, it is possible to increase the data width (from 4 bits in 1970 to 128 bits or more today), to build more complex operators with more complex instruction sets, and to enhance memory and caches management. This rise of the number of transistors also made possible several improvements in processor design, including instructions pipelines (dividing operators into shorter sections, processing more than one instruction at a time and thus increasing the frequency), superscalar architectures (processing several instructions at a time), out-of-order execution (permuting instructions to better use processor pipelines).

Figure 1: (a) Von Neumann architecture: at each clock cycle, one instruction is fetched from the memory, decoded, and then executed on the ALU. Results are stored back in the registers (very fast, but small, local memory) or in the main memory (slower access). Real processors also contain one or several caches (not shown here) between the registers and the main memory. (b) A processor with two cores. Note that each core has its own instruction decoding, and so they can operate independently as long as they do not access the same places in the memory. (c) A SIMD processor. A single instruction decoding unit controls several processing units. Only one flow of instructions is executed, but most of the surface area of the chip is devoted to actual computation. Frequencies and the "power wall". A common misformulation of the Moore's law is that the "computational power" of the processors doubles every 18 or 24 months. Of course, all techniques cited in the previous paragraph aim to obtain more computational power. Most notably, between 1970 and 2000, the frequencies of the processor steadily increased, resulting in major speed gains. Between 1990 and 2000, Intel estimates that the 75 x gain in computational power was composed by a 13 x gain due to the frequencies and a 6 x gain due to processor design improvements [Gelsinger, 2001].

Since a few years, issues with heat dissipation prevent the processors from having higher frequencies. Figure 2 displays some power consumptions of common CPUs. The last Intel mono-core processor, the Pentium 4, achieved a power density of about 100 W/cm 2 . If this increase in the frequencies had continued between 2000 and 2010, the thermal density of some processors could today approach the one of a rocket nozzle [Gelsinger, 2001].

Multicore processors. One of the answers to maintain the growth of computational power despite frequencies limits is parallel processing. Parallelism is not new: research in distributed algorithms has more than 40 years, with many specialized architectures for high-performance computing. Today, this field of research reaches commodity hardware through grids (several processors or computers) or multicore processors (several cores in a processor).

The idea of multicore processors is to have two, four, eight, or even more processors in one (Figure 1b). The high majority of personal computers are provided with this kind of processor at the end of the 2000's. Basically, this can be seen as a "grid on a chip", each core having its own instruction flow. This allows to execute either the same piece of code or different instructions on the different cores, on the same data or on different data.

Data-parallelism and SIMD. To compute the same piece of code on different data, another idea can be implemented. Instead of several processors, we can use a single instruction unit and several processing units (Figure 1c). In this case, the same piece of code is executed at the same time by several processing units. This SIMD paradigm (Single Instruction, Multiple Data) dates from the supercomputers of the 1980's. It becomes really interesting if the same task has to be executed on different data. No surface area is wasted in unnecessary control logic: the processor can be filled with efficient computing units.

Towards manycore processors

Multicore processing and SIMD processing can be combined to obtain a higher level of parallelism (Figure 3). This is what is implemented in current manycore processors, including Graphical Processing Units (GPUs) and Cell/BE processor. Depending on the design of the processor, emphasis is put either on the multicore or on the SIMD. GPU processors. Today, GPUs allow everyone to have some teraflops of cheap computing power in its personal computer. High-end GPUs, for no more than $500, embed far more arithmetic units than a CPU of the same price. For example, the NVIDIA GTX 285 has 30 cores, each one being a SIMD unit on 8 x 32 bits at about 600 MHz, whereas the ATI Radeon 4890 has 10 cores, each one being a SIMD unit on 80 x 32 bits, at about 850 MHz.

The GPUs also embed memories with different sizes and capabilities. For example, the NVIDIA GTX 285 has 1 or 2 GB of global memory, additional global memories (constant and texture memories), and, for each core, 16 KB of local memories. A detailed presentation on some GPUs is available in the review of [START_REF] Fatahalian | A Closer Look at GPUs[END_REF].

The CPU/GPU convergence. Recent trends blur the line between GPUs and CPUs: CPUs have more and more cores, and cores in GPUs have more and more functions. Some processors are now designed both for graphics and high performance computing.

The Cell Broadband Engine Architecture [START_REF] Gschwind | Synergistic processing in cell's multicore architecture[END_REF], or Cell/BE, is the chip included in the Playstation 3. It was designed by IBM and Sony and released in 2006. The Cell includes one CPU-like core, the Power Processor Element (PPE), and 8 GPU-like cores, the Synergistic Processing Elements (SPE), each one being a SIMD unit on 8 x 16 bits.

Larrabee [Seiler et al., 2008] is a Intel graphic processor project. It can be considered as a hybrid between a CPU and a GPU. Basically, each Larrabee core is a simple Pentium core running an extended version of the x86 instruction set, but it also includes a 16 x 32 bits SIMD vector processing unit. The Larrabee was supposed to be released in 2010, but Intel recently announced (in December 2009) that the release is delayed, without giving any new launch date.

Intel argues that the Larrabee processor is more flexible than GPUs. However, next-generation GPUs from NVIDIA (Fermi architecture) and ATI will also have greater flexibility. In 2006, AMD (CPU manufacturer) bought ATI (GPU manufacturer), and their next project, Fusion, is also expected to mix CPU-like and GPU-like cores.

General Purpose computation on GPU. Since the 1980's, graphics hardwares have been used for scientific computations (see for example [START_REF] Trendall | General calculations using graphics hardware with applications to interactive caustics[END_REF]). The revolution concerning today's "General-Purpose Computation on Graphical Processor Units" (GPGPU) is mainly due to the increasing number of cores, making those chips potentially very efficient, and the availability of programming libraries for developers non specialized in graphics computing (see in Methods). Some specialized GPUs do not have any video output: such cards are completely dedicated to general computation !

Methods

All manycore processors presented in the previous paragraph have a large number of cores, each of them being made of a lot of processing sub-units. An algorithm for such hardware thus requires a high level of parallelism, that is computations to dispatch to several cores and sub-units. This section details how to program such parallelism, and the next section explains in what kind of bioinformatics applications such parallelism can be used to obtain interesting speedups.

2.1 From GPU tweaks to OpenCL.

General-Purpose Computing on GPU, or GPGPU, was firstly done by tweaking graphics primitives. Operations on color components and pixels have been used to perform, for example, linear algebra computations. Such techniques require expertise in graphics. One of the first popular abstractions proposed for GPU programming is BrookGPU [START_REF] Buck | Brook for gpus: Stream computing on graphics hardware[END_REF], developed from 2003 at the Stanford University. It allows to compile and run code in the Brook Stream processing language on NVIDIA and ATI cards.

The CUDA libraries, released by NVIDIA in 2007 1 , deeply popularized GPGPU. CUDA is an extension of plain C/C++ (Figure 4), and does not require any knowledge in graphics. Hundreds of applications in various computer science domains are now using CUDA. In CUDA/OpenCL programs, there are some special functions, called kernels, that are executed on the GPU (Figure 4). Due to the architecture of the GPU, such functions are limited: for example, as there is no stack, there are no recursive functions. Ideally, the most intensive computing tasks of an application should be mapped into a kernel. The same kernel is then executed many times by different work-items, working on different data sets. The host program calls the kernel, but also ensures that data is transferred to and back from the GPU.

Programming SIMD work-items.

The CUDA/OpenCL kernels on Figure 4 contain an addition instruction,

c[i] = a[i] + b[i].
As all work-items execute this same instruction, this model is exactly the SIMD paradigm (Figure 1c). This is implicit SIMD: the compiler generates code that simultaneously controls all sub-units.

Similar SIMD instructions are available in the common CPUs, but on a much smaller data width. For example, the original SSE instruction set, published in 1999, allows to consider 4 "packed" 32bits floats in one 128-bit machine word. The SSE assembly instruction ADDPS adds two float vectors with one single operation. This instruction can also be compiled from a C "intrinsic" instruction mm add sd. Using SSE assembly or intrinsic operations is explicit SIMD. Even if it provides more control to the developer, it is not as useful as CUDA or OpenCL for the programmer.

for (i=0; i<n; i++) { if (d[i] > 53) c[i]=a[i]+b[i]; else c[i]=d[i]; } for (i=0; i<n; i++) { unsigned int mask=(d[i] > 53) ? 0xffffffff:0; c[i] = (a[i]+b[i]) & mask) // then | (d[i] & ~mask); // else }

Branches and divergence.

In the SIMD model, at one clock cycle, every arithmetic unit executes the very same instruction, allowing efficient operations on large vectors. But this is not always wanted. For example, depending on a condition in a loop, one could want to execute different instructions (Figure 5, on the left).

Such branches can be software simulated in SIMD extensions, by actually running both branches, and finally selecting with a mask the correct value (Figure 5, on the right). Although half of the computations are "wasted" in those branches, a good overall speedup can still be obtained. On some hardware, it is possible to temporarily disable some processing sub-units into a SIMD group. This leads to hardware simulated branches: at a given time, only some sub-units can be working in the SIMD processor [START_REF] Lorie | Method for conditional branch execution in SIMD vector processors[END_REF].

In CUDA/OpenCL, the developer directly writes a code similar to the left of Figure 5, and the compiler produces code suitable for the GPU, including software or hardware simulation. When, within a same work-group, all work-items do not pass the same branch, a divergence occurs and some work-items are stalled (Figure 6) [START_REF] Coon | [END_REF]Lindholm, 2008, Collange et al., 2009]. Of course, a code with a lot of divergences between work-items results in very poor performance. But, knowing this limitation, it is more simple to write such instructions in CUDA/OpenCL than to explicitly design SIMD masks.

Work-groups.

Current GPUs are not purely SIMD, because work-items are gathered into independent work-groups. There is some global memory available for their communication. This logical organization implied by CUDA/OpenCL has different physical implementations in actual hardware. For example, the exact relation between the work-groups depends on the efficiency of the global memory and its caching mechanisms. On the current GPU architectures, it is advised to design truly independent work-groups, whereas the Larrabee, with global cache coherence, will allow a more traditional programming with independent threads.

Results

Recent parallelizations on GPUs for sequence analysis problems achieve speedups up to 100 x compared to a serialized one-core version. This section reviews recent results in different bioinformatics applications.

The first part of this section deals with Smith-Waterman sequence alignments. Manycore processors deliver speeds enabling to run exact dynamic programming (DP) computations rather than incomplete heuristics. The second part explains other algorithms manipulating sequence data, including RNA algorithms and algorithms on weight matrices. Those applications use string algorithmics and DP computations similar to those used by sequence alignments, but on other objects or with different dependencies. Finally, the last part presents a few other applications where kernels deal with non-sequence data.

Smith-Waterman sequence alignments

In this section we describe how the Smith-Waterman (SW) algorithm for scanning of protein sequence databases can be efficiently mapped onto some manycore architectures. Mapping onto any SIMD architecture requires choosing a fine-grained SIMD vectorization. Mapping to the Cell/BE or a GPU architecture further requires coarse-grained distribution on available cores.

Smith-Waterman Algorithm. The Smith-Waterman (SW) algorithm computes the optimal local pairwise alignment [START_REF] Smith | Identification of common molecular subsequences[END_REF] of two given sequences S 1 and S 2 of length l 1 and l 2 using DP with the following recurrence relations.

     E(i, j) = max{E(i, j -1) + g e , H(i, j -1) + g o } F (i, j) = max{F (i -1, j) + g e , H(i -1, j) + g o } H(i, j) = max{0, E(i, j), F (i, j), H(i -1, j -1) + sbt(S 1 [i], S 2 [j])}
where sbt() is a substitution matrix such as BLOSUM62 [START_REF] Henikoff | Amino Acid Substitution Matrices form Protein Blocks[END_REF], g o is the gap opening penalty, and g e is the gap extension penalty. The above recurrences are computed for 1 ≤ i ≤ l 1 and 1 ≤ j ≤ l 2 and are initialized as

H(i, 0) = H(0, j) = E(i, 0) = F (0, j) = 0 for 0 ≤ i ≤ l 1 and 0 ≤ j ≤ l 2 .
The score of the optimal local pairwise alignment is the maximal score in matrix H (maxScore). The actual alignment can be found by a traceback procedure. However, for SW-based protein sequence database scanning, we just need to compute maxScore for each query/database sequence pair. Database sequences are then ranked according to their maxScore value and the top hits are displayed to the user. Note that the score-only computation can be done in linear space and does not require storing the full DP matrix. The data dependency in the SW DP matrix (Figure 7) implies that all cells in the same minor diagonal can be computed in parallel. Mapping onto SIMD registers. Between 1997 and 2007, three approaches have been proposed to vectorize the SW algorithm on CPUs with SIMD instructions sets (Figure 8): (i) minor diagonal approach [Wozniak, 1997], (ii) column-based approach with sequential memory layout [START_REF] Rognes | Six-fold speed-up of smithwaterman sequence database searches using parallel processing on common microprocessors[END_REF]; and finally (iii) column-based approach with striped memory layout [Farrar, 2007].

In order to calculate H(i, j), the value sbt(S 1 [i], S 2 [j]) needs to be added to H(i -1, j -1). The main challenge for any vectorized SW implementation is to avoid performing this table lookup for each element in a SIMD register. Therefore, all three approaches shown in Figure 8 calculate a query profile parallel to the query sequence beforehand. In the minor diagonal approach, the query profile is then used to arrange a SIMD register with all required sbt()-values. The advantage of this Figure 9: Inter-task vs intra-task parallelization, from [Liu et al., 2009a].

approach is that it does not require any conditional branches in the inner loop. However, it has the disadvantages that this SIMD register has to be updated in each iteration step and that registers are not fully utilized for diagonals at the beginning and at the end. Both column-based approaches have shown higher efficiency due to the simplified dependency relationship and parallel loading of the vector scores from memory. However, a disadvantage introduced by processing in column-based order with sequential memory layout is that a conditional branch is introduced in the inner loop for computing matrix F .

The advantage of the striped layout compared to the sequential layout is that data dependencies between vector registers are moved outside the inner loop. For instance, when calculating vectors for the DP matrices H or F with the sequential layout, the last element in the previous vector has to be moved to the first element in the current vector. When using the striped query layout, this needs to be done just once in the outer loop when processing the next database sequence character.

Implementation on Cell/BE. To fully exploit the capability of the Cell/BE, a parallel SW implementation has to take advantage of the SIMD registers of each SPE. Overall, the column-based approach with striped memory layout is the most efficient of the three approaches and is consequently utilized in publicly available Cell/BE SW implementation such as SWPS33 [START_REF] Szalkowski | SWPS3 -fast multi-threaded vectorized smith-waterman for IBM Cell/B.E. and x86/SSE2[END_REF] and CBESW4 [START_REF] Wirawan | CBESW: Sequence alignment on the Playstation 3[END_REF]. However, it should be mentioned that actual performance depends on the utilized scoring scheme (i.e. substitution matrix and gap penalties).

For coarse-grained parallelization across SPEs, the database can be partitioned into non-overlapping workloads of similar size which are then distributed from the PPE to the different SPEs using multi-threading. Furthermore, due to the limited SPE local store (256KB) a partitioning of the DP matrix for long query sequences is required in SWPS3 [START_REF] Szalkowski | SWPS3 -fast multi-threaded vectorized smith-waterman for IBM Cell/B.E. and x86/SSE2[END_REF].

Intra-task or inter-task parallelization on GPUs. There are two basic approaches to map SW-based protein sequence databases scanning onto manycore GPUs: (i) inter-task and (ii) intratask. Considering the SW alignment of a query sequence/database sequence pair as a basic task, the inter-task approach assigns each task to one work-item while the intra-task approach assigns each task to one work-group (Figure 9). Implementation results have shown that inter-task parallelization generally achieves higher performance at the cost of higher memory consumption. Therefore, CUDASW++5 [Liu et al., 2009a] uses this method for most alignments with the exception of very long database sequences, which use the intra-task approach. To achieve good load balancing within a work-group, all work-items within the same work-group can be assigned database sequences of roughly equal length.

Memory optimization on GPUs. In the CUDA programming model, global memory has to be accessed in a coalesced fashion to achieve high efficiency. Inter-task parallelization uses a coalesced global memory access pattern. A memory slot is allocated to a work-item in a work-group and is indexed top-to-bottom. Therefore, the access to the memory slot uses the same index for all work-items in a work-group. The coalesced global memory arrangement is used for both database sequences and intermediate results of the SW computation.

Fast local memory and registers can be exploited in the SW computation as follows. Instead of computing only one SW cell by each work-item for each global memory access, a whole cell block of size n × n can be computed. In this case, the computation of n cells in a column (or row) of a cell block only requires one load and one store operation to the global memory instead of n load and n store operations. Since one global access takes hundreds of clock cycles, this method leads to a significant performance improvement. However, the size of the cell block is limited by the amount of local memory and registers available; e.g. CUDASW++2.0 uses a size of 8 x8 on a GTX280. Furthermore, the gap penalties and query sequence are stored in a constant global memory, while the scoring matrix is loaded to local memory.

Peak performance SSEARCH [START_REF] Pearson | Improved tools for biological sequence comparison[END_REF] CPU (2.0 GHz Xeon) ∼ 0.1 GCUPS non-SIMD SW implementation [Farrar, 2007] CPU (2.0 GHz Xeon) 3.0 GCUPS striped layout, SSE2 SIMD instructions CBESW [START_REF] Wirawan | CBESW: Sequence alignment on the Playstation 3[END_REF] Cell/BE 3.6 GCUPS striped layout, long queries (800) SWPS3 [START_REF] Szalkowski | SWPS3 -fast multi-threaded vectorized smith-waterman for IBM Cell/B.E. and x86/SSE2[END_REF] Cell/BE 9.3 GCUPS striped layout, long queries [Farrar, 2009] Cell/BE 16 GCUPS striped layout, long queries (32K) SWCUDA [START_REF] Manavski | CUDA compatible GPU cards as efficient hardware accelerators for Smith-Waterman sequence alignment[END_REF] GeForce 8800 GTX ∼ 1.5 GCUPS inter-task [START_REF] Ligowski | An efficient implementation of Smith-Waterman algorithm on GPU using CUDA, for massively parallel scanning of sequence databases[END_REF] GeForce 9800 GX2 7.5 / 14.5 GCUPS CUDASW++2.0 [Liu et al., 2009a] GTX 280/295 16.8 / 28.8 GCUPS inter-task + intra-task, short/long queries [START_REF] Pearson | Improved tools for biological sequence comparison[END_REF], the best implementations now achieve more than a 100 x speedup.

However, peak performances do not measure every aspect of the implementations. For example, in the study of [Liu et al., 2009a], the authors compare the performance of their CUDASW++ on a single-GPU GTX 280 and a dual-GPU GTX 295 versus the performance of SWPS3. SWPS3 achieves a peak performance of around 9.3 GCUPS for longer queries, but is very inefficient for short queries. Due to the inter-task parallelization, the CUDASW++2.0 single-GPU version has a relatively constant performance with an average of 16.3 GCUPS. The dual-GPU version performance is relatively low for short query sequences, but increases up to 28.8 GCUPS for longer ones.

Algorithms on sequence data

The algorithms discussed in this section alse use kernels working on raw sequence data, but address other bioinformatics problems. The first two algorithms concern DP recurrences, whereas the last one is about weight matrices computations.

RNA folding. Determining the folding of RNAs is important for the study of noncoding RNAs. The secondary structure is a succession of base pairings, most often A/T, C/G and G/U. This secondary structure is a basis of the tertiary structure of the RNA. Moreover, comparative genomics has shown that, through evolution, this structure is more conserved than the nucleic sequence.

The Nussinovs algorithm [START_REF] Nussinov | Algorithms for loop matchings[END_REF] finds the RNA structure with the maximum number of base pairs. However, complete RNA folding algorithms are typically based on energy minimization, and include energies of stacking regions (or helices), bulge loops, internal loops, hairpin loops and multiple loops, for example as in the full Turner model [START_REF] Matthews | Expanded sequence dependence of thermodynamic parameters improves prediction of rna secondary structure[END_REF], included in the mfold/unafold packages [Zuker, 2003].

In [START_REF] Rizk | GPU accelerated RNA folding algorithm[END_REF], the authors develop an optimized GPU implementation of mfold/ unafold, by computing in parallel all cells from the diagonal of the dynamic programming matrix (Figure 10). At each cell (i, j), if bases i and j can pair, huge computations are run to search for additional RNA structures. However, in the mfold/unafold implementation, only 6 over the 16 possible combinations of bases launch these computations. That means that only 6/16 cells of the DP matrix lead to further computations. This breaks the SIMD model of the GPU, as neighbor work-items would probably diverge. To prevent this, the authors pack together these 6/16 cells, and compute all of them in the same time. Finally, they obtain a 17 x speedup on a GTX 280 over a one core version on a 2.66 GHz Xeon.

Figure 10: Dependencies in the RNA folding DP matrix, as implemented in [START_REF] Rizk | GPU accelerated RNA folding algorithm[END_REF]. Each cell depends on its underlying triangle. All cells in the same major diagonal can be computed in parallel. Generic dynamic programming. The Algebraic Dynamic Programming (ADP) framework [START_REF] Steffen | Versatile and declarative dynamic programming using pair algebras[END_REF] is able to encode different DP problems. Starting from an abstract grammar (Figure 11) and an evaluation algebra, the ADP compiler generates the DP dependencies and recurrences, and finally translates them into C (Figure 12). This abstraction is especially useful when the optimization problem includes lots of subcases, resulting in large DP recurrences.

In [START_REF] Steffen | GPU parallelization of algebraic dynamic programming[END_REF], a CUDA backend of the ADP compiler is proposed6 . In all ADP programs, all results are combined from results of shorter subsequences. Therefore, the calculation of a table element (i, j) depends only on results that lie in the "underlying triangle" under (i, j), and, [START_REF] Steffen | GPU parallelization of algebraic dynamic programming[END_REF]. The same abstract grammar and algebra, in RNAfold.lhs, is used to produce plain C or CUDA code. as in the case of RNA folding, such elements can be computed in parallel (Figure 10).

The authors apply this technique on several RNA tools. On pknotsRG, a RNA pseudo-knot detection application [START_REF] Reeder | pknotsRG: RNA pseudoknot folding including near-optimal structures and sliding windows[END_REF], they obtain speedups up to 7 x. On RNA folding, they achieve a speedup up to 9 x, far below the implementation of [START_REF] Rizk | GPU accelerated RNA folding algorithm[END_REF]. This solution is thus less efficient than a manually crafted implementation but the advantage is definitely the genericity of the approach, even if for now only a few DP problems are encoded like that.

Position Weight Matrices algorithms. Position Weight Matrices (PWMs) model approximates patterns from a set of sequences, as for instance in transcription factor binding sites, splicing sites or protein domains. Given a finite alphabet Σ and a positive integer m, a PWM M is a matrix with |Σ| rows and m columns (Figure 13). The PWM associates to each word u = u 1 u 2 . . . u p the score m p=1 M (p, u p ). Given a score threshold α, it is possible to compute the probability to achieve a score greater than α: this is the P-value.

There are three usual problems on PWMs: finding positions where a pattern occurs (the scan), assessing a P-value for each occurrence, and comparing patterns. In [Giraud and Varré, 2009], a CUDA parallelization achieves speedups of 21 x for the brute-force parallelization of the scan and 77 x for the P-value computation (on a NVIDIA GTX 280)7 .

The parallelization of the scan is straightforward, distributing positions of the text to different work-items. To reduce memory transfers, work-items in a same work-group process interleaved positions (left of Figure 14).

Computing the P-value is NP-hard [START_REF] Touzet | Efficient and accurate p-value computation for position weight matrices[END_REF]. The usual algorithm proceeds by recursion [Staden, 1989], computing the score distribution Q M , where Q M (s) is the probability to achieve exactly the score s. Then, for a given score s ′ , the P-value is obtained with the equation: P-value M (s ′ ) = s≥s ′ Q M (s). The distributed algorithm of [Giraud and Varré, 2009] computes an approximation of the score distribution Q M by splitting the matrix into submatrices. Q is computed for each submatrix and the resulting score distributions are merged (right of Figure 14). Here the brute-force computation of Q for each submatrix takes advantage of the GPU architecture: each work-item evaluates a large set of 4 ℓ words (where ℓ is the length of the submatrix), with no divergence, and no communication nor memory accesses during the main computation. The best speedups are obtained for larger matrices. Unfortunately, for smaller matrices, which are the ones stored into databases, the execution time remains lower using the CPU.

Other applications

We now briefly review other bioinformatics algorithms that have been parallelized on GPUs. Even if some problems still input sequence data, here the kernels do not directly deal with sequence data.

Indexing structures. As some string algorithms are intrinsically complex, several heuristics have been proposed on CPUs, including the popular BLAST [START_REF] Altschul | Basic local alignment search tool[END_REF] for SW sequence alignments. Heuristics often build indexes when preprocessing a part of the input. Some researches try to put this kind of heuristics on manycore processors.

MUMmer [START_REF] Kurtz | Versatile and open software for comparing large genomes[END_REF] is an heuristic that uses maximal exact matches (MEM) to seed SW sequence alignments. The MEM detection uses a large suffix tree [Gusfield, 1997]. The paper [START_REF] Schatz | Highthroughput sequence alignment using graphics processing units[END_REF] proposes a CUDA parallelization of MUMmer, splitting the text in 8 Mb pages. [START_REF] Schatz | Optimizing data intensive GPGPU computations for DNA sequence alignment[END_REF] further optimizes this implementation. This latter article contains an exhaustive study of 2 7 = 128 combinations of design choice for their CUDA implementation8 . [START_REF] Shi | Accelerating error correction in high-throughput short-read dna sequencing data with cuda[END_REF] elaborates an error correction algorithm for read mapping, based on spectral alignment. Their algorithm uses a specific data structure, a Bloom filter, that is a set of hashing tables with multiple keys. The parallelization is done on the querying of this Bloom filter. Authors report speedups up to 19 x compared to an equivalent CPU algorithm.

Phylogeny. [START_REF] Charalambous | Initial experiences porting a bioinformatics application to a graphics processor[END_REF] was the first bioinformatics application on GPU. Authors study the RAxML phylogenetics program, and, after profiling, decide to parallelize a particular loop. The speedup (1.2 x) was obtained with a Brook kernel on a NVIDIA GeForce FX 5700. [START_REF] Suchard | Many-core algorithms for statistical phylogenetics[END_REF]] studies statistical phylogenetics. The goal of such methods is to evaluate the likelihood of a set of sequences given a phylogenetic tree and a model of evolution. The difficulty is in the computation of the likelihood itself. Their implementation uses very precise memory accesses, obtaining a speedup up to 60 x on a NVIDIA GTX 280 GPU9 .

Multiple sequence alignment. ClustalW [START_REF] Larkin | Clustal W and Clustal X version 2.0[END_REF] is a multiple sequence alignment program. This progressive alignment technique aims to align the sequences following a guided tree, computed from distance matrices. [START_REF] Liu | GPU-ClustalW: Using graphics hardware to accelerate multiple sequence alignment[END_REF] proposes a parallelization of the first phase that builds pairwise alignments, with a method similar to the one presented in the first part of this section.

The study of [START_REF] Liu | Parallel reconstruction of neighbor-joining trees for large multiple sequence alignments using cuda[END_REF] elaborates a parallelization of the second phase, the Neighbor-Joining algorithm [START_REF] Saitou | The neighbor-joining method: a new method for reconstructing phylogenetic trees[END_REF]. Here the parallelization is on the distance matrix computation, where each cell can be computed independently. The study includes a discussion on the optimal number of work-items and work-groups, as well as on a good memory management taking advantage of the symmetry of the distance matrix. A speedup up to 26 x is finally obtained on a NVIDIA GTX 280 card.

Motif finding. CUDA-MEME [START_REF] Liu | Cuda-meme: Accelerating motif discovery in biological sequences using cuda-enabled graphics processing units[END_REF] is a CUDA implementation of the popular MEME tool for finding regulatory regions in DNA sequences (so-called motifs). It achieves a speedup of up to 20.5 over the sequential MEME code on a GTX280.

Hidden Markov Models profiles. In [START_REF] Walters | Evaluating the use of gpus in liver image segmentation and hmmer database searches[END_REF], the authors propose a GPU version of HMMER [Eddy, 1998]. HMMER is a popular software that aims to compute a Hidden Markov Model (HMM) from a set of aligned protein sequences allowing to search for occurrences of similar proteins in databases. The paralellization focused on the heart of the hmmmsearch tool that uses the Viterbi algorithm. A speedup up to 38 x is obtained depending on the size of the HMM.

Cell molecules simulation. [START_REF] Roberts | Long time-scale simulations of in vivo diffusion using GPU hardware[END_REF] simulates diffusion of molecules in the cell10 . The cell topology is modeled as a lattice, and particles follow a random walk. The speedup is not so high (2.4 x on a NVIDIA GTX 280), but this study has interesting remarks on the limits of the GPU for this problem.

Discussion

The previous section shows that the development of manycore algorithms for bioinformatics has started, taking into account architecture details of such processors. Here we point out what difficulties might be encountered, and discuss how to exploit the processing power of manycore processors in bioinformatics analysis. We also mention some challenges in designing parallel algorithms.

Challenges in parallel algorithmics

Most of the applications presented in this chapter exhibit data-parallelism on sequence data: the same instruction flow is applied to every chunk of data. This is especially the case for DP algorithms (Smith-Waterman sequence comparison, RNA folding, Algebraic Dynamic Programming), but also for other algorithms on sequence data. In all those algorithms, different work-items work on different small sections of the input data. These algorithms are the best ones to parallelize, with almost no divergence between work-items in a SIMD sub-unit. As mentioned in the previous section, researches try to further optimize those algorithms by better memory access patterns or better work-item/work-group balancing.

Another active field of research is to design elaborated data structures intrinsically parallel, as, for example, parallel structures for trees allowing large independent parallel executions. Works on suffix trees [START_REF] Schatz | Highthroughput sequence alignment using graphics processing units[END_REF]Trapnell, 2009] are a first step in this direction, but the speedups are not so high for the moment (around 5 x). Techniques can be borrowed from graphics community, where researches on ray-tracing also try to conceive good parallel algorithms for trees [START_REF] Popov | Stackless kd-tree traversal for high performance GPU ray tracing[END_REF].

Of course, new chips with more independent cores could allow more functionalities , but at the price of less dense computational units. Anyway, there is place for research on improved parallel algorithms, using better the resources of manycore architectures.

Challenges for bioinformatics analysis

Final users of algorithms and methods developed for manycore processors are supposed to be biologists or bioinformaticians who analyze data. Nevertheless, most of these programs are research prototypes rather than end-user applications. In some cases, they need better integration to be, for example, compliant with existing programs, using the same data formats for input and output.

Biomanycores. Biomanycores11 is intended to be a collection of GPU and manycore bioinformatics tools [Varré et al., 2009] interfaces to Bio * frameworks like BioJava [START_REF] Holland | BioJava: an open-source framework for bioinformatics[END_REF], BioPerl [START_REF] Stajich | The Bioperl toolkit: Perl modules for the life sciences[END_REF], and Biopython [START_REF] Cock | Biopython: freely available Python tools for computational molecular biology and bioinformatics[END_REF] (Figure 15). We mentioned in the previous section several parallel bioinformatics applications already available, especially for sequence alignment, that can provide important speedups with common GPUs. However, those programs are almost not used, or even known, by biologists or bioinformaticians. Biomanycores is based on Bio* frameworks, that are widely used. These frameworks are a collection of tools that allow the user to create specific analysis pipelines for his data. With Biomanycores, the user can keep his pipeline and increase his capacity of analysis by replacing the sequential tool by a parallel one. For example, suppose that a biologist has a large number of sequences to analyse, and wants to format the results. By simply replacing the call to the alignment function by a call to the parallel one, he can have a speedup of around 10x. Concretely, Biomanycores is a repository of open-source parallel bioinformatics code (in Opencl or Cuda) and interfaces to use these programs from Bio* frameworks. It tries to bridge the gap between researches in highperformance-computing and platforms used by bioinformaticians and biologists, by giving access to high-performance prototypes through Bio* frameworks. The project started with three different applications: Smith-Waterman [START_REF] Manavski | CUDA compatible GPU cards as efficient hardware accelerators for Smith-Waterman sequence alignment[END_REF], pKnotsRG and weight matrices scan (see in section Results). Contributors are welcome to upload their own applications. Of course, the speedup using those APIs is lower than the one obtained by using directly the applications or a manually crafted pipeline, but the advantage is the ease of use inside standard frameworks.

Conclusion

Initiated by multicore and graphics processors, the trend of having more cores in a processor will surely keep on during the next years. As explained in the introduction, in the future, the increase of computational power will come from a higher number of processors (either in CPUs and in GPUs) rather than higher frequencies. On the other hand, the amount of biological data produced is growing more and more rapidly, and increases the need for high performance bioinformatics applications. Parallel algorithms could then become the standard, since parallelism is necessary to fully exploit the power of such processors. New tools like CUDA or OpenCL, as detailed in the second section, makes the use of GPUs for general programmning easier, allowing the user to benefit from the computational power of GPUs.

In the third section, we have studied some parallel bioinformatics applications. Parallelism is particularly suitable for bioinformatics algorithms based on dynamic programming, like sequence comparisons. Those applications can achieve good speedups, even on common (and cheap) GPUs. But the use of parallelism is not always straightforward, and more work is needed, for example to design parallel data structures or to find new parallel methods. Bringing parallelism to biologists and bioinformaticians is another challenge, the Biomanycores project is a first step toward this goal.

We believe that the developpement of manycore processors (on commodity hardware) opens a new era in parallel processing, and is an opportunity for high performance computing in bioinformatics.

Figure 2 :

 2 Figure 2: Number of transistors and power consumption of some Intel chips. Scales on Y-axis are logarithmic. Both values have grown at an exponential rate between 1970 and 2005. Data from http://www.intel.com/.

Figure 3 :

 3 Figure3: A manycore processor, like current GPUs, combining several cores, each one being a high-density SIMD unit. Again, real processors also contain local memories or caches at different levels between the registers and the main memory.

Figure 4 :

 4 Figure 4: Vector addition in CUDA (left) and OpenCL (right). The addv function is the kernel, that is the function executed on each ALU. Each work-item executes this addition on different data, as showed by the instruction modifying the value of i. Full codes include initialization of the card, transfer of arrays a and b from the host to the GPU, and transfer of the resulting array c back to the GPU.

Figure 5 :

 5 Figure 5: Simulating branching in SIMD. The branching on the left can be simulated by the code on the right. Both branches are executed, and the mask finally chooses the good value for each item.

Figure 6 :

 6 Figure 6: Hardware simulation on branches in a SIMD unit with 8 sub-units. Here only workitems 2 and 6 execute the conditional code: all other work-items are stalled, reducing the overall performance.

Figure 7 :

 7 Figure 7: Dependencies in the Smith-Waterman DP matrix. Each cell depends on its left, upper, and upper-left neighbors. All the cells in the same minor diagonal can be computed in parallel.

Figure 8 :

 8 Figure 8: SIMD vectorization approaches: (i) minor diagonal; (ii) column-based with sequential layout; (iii) column-based with striped layout. Each blank box depicts a parallel profile stored into a SIMD register, here with 4 different values.

Figure 11 :

 11 Figure 11: Excerpt from the ADP grammar for RNA folding. Starting from the axiom (struct), several productions detail how to generate different RNA structure components.

Figure 12 :

 12 Figure12: ADP workflow, from[START_REF] Steffen | GPU parallelization of algebraic dynamic programming[END_REF]. The same abstract grammar and algebra, in RNAfold.lhs, is used to produce plain C or CUDA code.

Figure 13 :Figure 14 :

 1314 Figure 13: A Position Weight Matrix (PWM) modeling the CREB1 transcription factor binding site (from the JASPAR database). The coefficients indicate affinities between letters and positions at the binding site: positive if the letter is over-represented in the binding site, else negative.From[Giraud and Varré, 2009].

Table 1 :

 1 Peak performance of several implementations of Smith-Waterman in GCUPS (Billion Cell Updates Per Second). All Cell/BE and GPU implementations were published in 2008 or 2009. Note that the average performance greatly depends on the input size, as well on the score matrix and gap penalties. The 14.5 and 16.1 GCUPS of Ligowski and CUDASW++ are on dual-GPU cards.Performance comparison. Table1lists some peak performance achieved by recent SW implementations, measured in GCUPS (Billion Cell Updates Per Second). In two years, this active field of research has brought a 10 x improvement on almost equivalent hardware. Compared to the non-SIMD plain SW implementation of

  . The goal is both to gather manycore programs and to propose

	agatcttaCGTAGTGACGTCtgccatgg agatctgGCGGGTGACGTGttgccatgg agatcttgGCGGGTGACGTTtctccatgg	M (i, x) = log 2	frequency of letter x at position i background frequency of letter x
	... agatctcggggTATGTTGACGCCatgg agatCTTCGTGACGTTctttgccatgg aGATCTTGACGTCcgcaggtaccatgg	A [ -4.85826 -0.06247 -4.85826 -0.46381 • • • ] C [ 0.37462 0.03957 -0.46793 -0.46793 • • • ] G [ 0.03957 -0.18232 0.22107 0.82531 • • • ]
		T [ 0.22314	0.22314	0.78009 -4.85826 • • • ]

The next big step was the specification of a new standard, OpenCL, by a large consortium in December 2008 2 . OpenCL is very close to CUDA, and is likely to gain a wide industrial support. OpenCL offers a unique programming interface to deal with different manycore processors, allowing a same code to be compiled and to run on different chips. But a compiler is needed for each architecture. In 2009, three different implementations, by NVIDIA, AMD/ATI and Apple, were already available, and OpenCL could become in the following years a standard for parallel programming.1 http://www.nvidia.com/cuda 2 http://www.khronos.org/opencl 6

http://www.inf.ethz.ch/personal/sadam/swps3/

http://sourceforge.net/projects/cbesw/

http://cudasw.sourceforge.net/

http://bibiserv.techfak.uni-bielefeld.de/adp/cuda.html

http://bioinfo.lifl.fr/cudapwm/

http://www.mummergpu.sourceforge.net

http://beagle-lib.googlecode.com/

http://www.ks.uiuc.edu/Research/gpu/

http://www.biomanycores.org