Laurent Hascoet

INRIA Val Érie Pascual

A The Tapenade Automatic Differentiation tool: principles, model, and specification

Keywords: G.1.4 [Quadrature and Numerical Differentiation]: Automatic differentiation, D.1.2 [Automatic Programming]: Program transformation, D.3.4 [Processors]: Compilers, Preprocessors, F.3.2 [Semantics of Programming Languages]: Operational semantics, Program analysis, G.1.6 [Optimization]: Gradient methods Algorithms, Performance Source transformation, Adjoint compiler

Tapenade is an Automatic Differentiation tool which, given a Fortran or C code that computes a function, creates a new code that computes its tangent or adjoint derivatives. Tapenade puts particular emphasis on adjoint differentiation, which computes gradients at a remarkably low cost. This paper describes the principles of Tapenade, a subset of the general principles of AD. We motivate and illustrate on examples the AD model of Tapenade, i.e. the structure of differentiated codes and the strategies used to make them more efficient. Along with this informal description, we formally specify this model by means of Data-Flow equations and rules of Operational Semantics, making this the reference specification of the tangent and adjoint modes of Tapenade. One benefit we expect from this formal specification is the capacity to study formally the AD model itself, especially for the adjoint mode and its sophisticated strategies. This paper also describes the architectural choices of the implementation of Tapenade. We describe the current performances of Tapenade on a set of codes that include industrial-size applications. We present the extensions of the tool that are planned in a foreseeable future, deriving from our ongoing research on AD.

INTRODUCTION

Computing accurate derivatives of a numerical model f : IR n → IR m is a crucial task in many domains of Scientific Computing. This is a computation-intensive task for which research and development of software tools are most wanted.

Automatic Differentiation (AD) is a collection of techniques to obtain analytical derivatives of differentiable functions, in the case where these functions are provided in the form of a computer program. Instead of using a mathematical representation, explicit or implicit, of these functions, AD uses the computer program itself as the basis to compute the derivatives. In general, AD techniques are made available to the end-users through specific software libraries or through software tools. For most applications, e.g. for optimization, AD competes with hand-coded derivatives, with inexact derivatives such as Divided Differences, or with derivative-free methods such as Evolutionnary Algorithms or Stochastic Methods. Research on AD thus focuses on:

-designing techniques and models to compute efficient analytical derivatives of functions given as programs, -developing libraries and tools that make these techniques directly available to endusers, and Tapenade A:3 implementation. Section 6 gives short practical notes on Tapenade and installation directions, followed by a synthetic view of the performances of derivative code produced on a selection of applications of all sizes. Section 7 concludes with developments planned.

ELEMENTS OF AUTOMATIC DIFFERENTIATION

Given a computer algorithm P (identified with a piece of program) that implements a function F : X ∈ IR n → Y ∈ IR m , AD builds a new algorithm (a program piece) P that computes derivatives of F by computing the analytical derivative of each elementary mathematical operation in P. We first observe that any run-time trace of P {I 1 ; I 2 ; . . . I p ; } computes the composition of elementary functions, one per instruction I k ,

f p • f p-1 • • • • • f 1 ,
which we can identify to F . This is of course assuming that P is a correct implementation of F , i.e. the discretization and approximation employed in P are sufficiently accurate and do not introduce non-differentiability.

Let us clarify the correspondence between the mathematical variables (X, Y . . . ) and the program variables found in P. As imperative programs classically overwrite their variables to save memory space, let us call V the collection of all the program variables of P and consider that each instruction I k (partly) overwrites V. With these conventions (this run-time trace of) P is indeed the program:

original program P Initialize V with X (I1) V := f1(V) ... (I k ) V := f k (V) ... (Ip) V := fp(V) Retrieve Y from V
At any given location in P, the program variables V correspond to one particular set, or vector, of mathematical variables. We will call this vector X k for the location between instructions I k and I k+1 . The set V is actually large enough to accomodate X, Y , or each successive X k . At each location, V may thus "contain" more than the X k but only the X k play a role in the semantics of the program. The program instruction V := f k (V) actually means taking from V the mathematical variables X k-1 before the instruction and applying f k (X k-1 ) to obtain X k . After I k , V corresponds to X k . The Initialize with and Retrieve from instructions in the program sketch define X 0 as X and Y as X p .

Since we identify F with a composition of functions, the chain rule of calculus gives the first-order full derivative, i.e. the Jacobian:

F (X) = f p (X p-1 ) × f p-1 (X p-2 ) × . . . × f 1 (X 0 ) .
It is thus possible in theory to adapt algorithm P so that it computes F (X) in addition to F (X). This can be done simply by extending locally each instruction I k with a piece of code that computes J k = f k (X k-1 ) × J k-1 . This transformation is not limited to straight-line code and can be applied to any program P, so that the extended algorithm P takes the same control decisions as P. Of course, derivatives are valid only if the control does not change in an open neighborhood around X. Otherwise, the risk is that AD may return a derivative in cases where F is actually non-differentiable. Keeping this caveat in mind, the adapted algorithm can return J p , the complete Jacobian F (X).

However, the J k are matrices whose height and width are both of the order of the number of variables in the original P, and may require too much memory space.

To work around this difficulty, one observes that the derivative object that is needed for the target application is seldom the full Jacobian matrix, but rather one of

F (X) × Ẋ or Y × F (X)
where Ẋ is some vector in IR n and Y is some row-vector in IR m . Moreover when F (X) is needed explicitly, it is very often sparse and can therefore be retrieved from a relatively small number of the above projections. This motivates the so-called tangent and adjoint modes of AD:

-Tangent mode: evaluate Ẏ = F (X) × Ẋ, the directional derivative of F along direction Ẋ. It expands as

Ẏ = f p (X p-1 ) × f p-1 (X p-2 ) × . . . × f 1 (X 0 ) × Ẋ . ( 1
)
Since Ẋ is a vector, this formula is most efficiently evaluated from right to left i.e., using mathematical variables:

X 0 = X Ẋ0 = Ẋ X 1 = f 1 (X 0 ) Ẋ1 = f 1 (X 0 ) × Ẋ0 . . . X k = f k (X k-1 ) Ẋk = f k (X k-1 ) × Ẋk-1 . . . X p = f p (X p-1 ) Ẋp = f p (X p-1 ) × Ẋp-1 Y = X p Ẏ = Ẋp
An algorithm Ṗ for this evaluation is relatively easy to construct, as the derivative instructions follow the order of the original instructions. Keeping the original program variables V to hold the successive X k , and introducing a set of new program variables V of the same size as V to hold the successive Ẋk , Ṗ writes:

tangent differentiated program Ṗ Initialize V with X and V with Ẋ ( İ1) V := f 1 (V) × V (I1) V := f1(V) ... ( İk ) V := f k (V) × V (I k ) V := f k (V) ... ( İp) V := f p (V) × V (Ip) V := fp(V)
Retrieve Y from V and Ẏ from V

Notice that each derivative statement İk now precedes I k , because I k overwrites V.

Tapenade A:5

-Adjoint mode: evaluate X = Y ×F (X), the gradient of the scalar function Y ×F (X) derived from F and weights Y . It expands as

X = Y × f p (X p-1 ) × f p-1 (X p-2 ) × . . . × f 1 (X 0 ) . ( 2
)
Since Y is a (row) vector, this formula is most efficiently evaluated from left to right i.e., with mathematical variables:

X 0 = X X 1 = f 1 (X 0 ) . . . X k = f k (X k-1 ) . . . X p = f p (X p-1 ) Y = X p X p = Y X p-1 = X p × f p (X p-1 ) . . . X k-1 = X k × f k (X k-1 ) . . . X 0 = X 1 × f 1 (X 0 ) X = X 0
However, an algorithm that evaluates these formula is not immediate to construct, as the derivative instructions will follow the inverse order of the original instructions.

Similarly to the tangent mode, we want the adjoint program to use only the original program's variables V plus a corresponding set of new program variables V, of the same size as V, to hold the successive X k . In that case, we see that e.g. X k-1 contained in V will be overwritten by X k and thus lost, before it is needed to evaluate X k × f k (X k-1 ). We will see later how this problem is solved, but let us keep in mind that there is a fundamental penalty attached to the adjoint mode that comes from the need of a data-flow (and control-flow) reversal.

Let us compare the run-time costs of the tangent and adjoint modes. The tangent differentiated algorithm Ṗ costs only a small multiple of the run-time of the original P.

The ratio R t varies with the given P, typically between 1 and 3. Using a simplified cost model that counts the number of memory accesses (fetch or store), additions, multiplications, and nonlinear operations, R t is always less than 2.5 [Griewank and Walther 2008, Section 4.5]. Similarly, the adjoint differentiated algorithm P costs only a small multiple of the run-time of P. The ratio R a also varies with P, typically between 5 and 10. With the same simplified cost model, R a is bound by 5 instead of 2.5, due to the different memory access pattern. If only multiplication and nonlinear operations are counted, R a = R t , but in practice P suffers from the extra penalty coming from the data-flow reversal. Figure 1 will help us compare the costs of computing the complete Jacobian F (X) by employing the adjoint mode, the tangent mode, or even divided differences. We use no sparsity property.

-With the adjoint mode, we obtain F (X) row by row by setting Y successively to each element of the Cartesian basis of the output space IR m . The run time for the full Jacobian is thus m × R a × runtime(P).

( ) [

] n inputs m outputs Adjoint Tangent Fig. 1. Elements of the Jacobian computable by tangent AD and adjoint AD -With the tangent mode, we obtain F (X) column by column by setting Ẋ successively to each element of the Cartesian basis of the input space IR n . The run time for the full Jacobian is thus n × R t × runtime(P). -With divided differences (F (X +h Ẋ)-F (X))/h, we also get F (X) column by column, although with a much poorer accuracy. The full Jacobian costs (n + 1) × runtime(P).

When n is much larger than m, the adjoint mode is recommended. In particular, this is the case when gradients are needed, e.g. in optimization or in inverse problems.

There are typically very few optimization criteria (often m = 1), and on the other hand n is often large as the optimization parameters may be functions, shapes, or other complex objects. In that case, no matter R a being higher than R t , the adjoint mode of AD is the only reasonable option. This is the most flagrant situation where adjoint AD outperforms tangent AD or divided differences.

Considering the design and implementation of AD tools, there are two main ways to code the algorithms Ṗ and P namely, operator overloading and program transformation. The implementation choice for Tapenade is program transformation, and to motivate this requires an elementary description of these two ways.

-Operator Overloading: if the language of P permits, one can replace the types of the floating-point variables with a new type that contains additional derivative information, and overload the arithmetic operations for this new type so as to propagate this derivative information along. See e.g. [START_REF] Forth | An efficient overloaded implementation of forward mode Automatic Differentiation in MAT-LAB[END_REF][START_REF] Walther | Getting started with ADOL-C[END_REF]. Schematically, the AD tool boils down to a library that defines the overloaded type and arithmetic operations. This approach is both elegant and powerful. The overloaded library can be quickly redefined to compute higher-order derivatives, Taylor expansions, intervals. . . By nature, evaluation of the overloaded operations will follow the original order of P. This is fine for the tangent mode, but implies some acrobacy for the adjoint mode, bearing severe consequences on performance and/or loosing a part of the elegance of the approach. -Program Transformation: One can instead decide to explicitly build a new source code that computes the derivatives. This means parse the original P, build an internal representation, then build the differentiated Ṗ or P. See e.g. [START_REF] Carle | ADIFOR 3.0 overview[END_REF][START_REF] Utke | OpenAD/F: A modular, open-source tool for Automatic Differentiation of Fortran codes[END_REF]. This allows the tool to analyze P globally, e.g. its data-flow, to optimize the differentiated code. This is similar to a compiler, but it produces source code. This approach is more development-intensive than Operator Overloading. Program Transformation AD tools are perhaps slightly more fragile and need more effort to follow the constant evolution of programming constructs and styles. On the other hand Program Transformation is the choice approach for the adjoint mode of AD which requires control-flow and data-flow reversal and where global data-flow analysis is essential for efficiency.

Tapenade A:7

These elements of AD are about all the background we need to describe a tangent mode AD tool. For the adjoint mode however, we need to address the question of dataflow reversal. We saw that the equations of the adjoint mode do not transpose immediately into a program because the values X k are overwritten before they are needed by the derivatives. Basically, there are two ways to solve this problem, sketched in figure 2, and a variety of combinations between them. Fig. 2. Data-flow reversal with Recompute-All approach (left) and Store-All approach (right). Big black dot represents storage of X 0 , big white dots are retrievals. Small black dots are values being recorded before overwriting. Small white dots are corresponding restorations.

-Recompute-All: For each derivative instruction computing X

k-1 = X k × f k (X k-1 ),
denoted as ← -I k , we recompute the X k-1 that it requires by a repeated execution of the original code, from the stored initial state X 0 to instruction I k-1 that computes X k-1 . The memory cost is only the storage of X 0 . On the other hand, the run time cost is quadratic in p.

-Store-All: Each time instruction I k overwrites a part of V, we push this part of V onto a stack just before overwriting. Later, we restore these values by popping them just before executing ← -I k . We get the following sketch of program P:

adjoint differentiated program P (Store-All) Initialize V with X and V with Y push(out(I1)) (I1) V := f1(V) ... push(out(I k )) (I k ) V := f k (V) ... push(out(Ip-1)) (Ip-1) V := fp-1(V) ( ← - Ip ) V := V × f p (V) pop(out(Ip-1)) ... pop(out(I k )) ( ← - I k ) V := V × f k (V) ... pop(out(I1)) ( ← - I1 ) V := V × f 1 (V) Retrieve X from V
We see two successive sweeps in P. The forward sweep -→ P is essentially a copy of P augmented with storage of overwritten values. The backward sweep ← -P is the computation of the derivatives, in reverse order, augmented with retrieval of recorded values. We define out(I k ) to be the subset of the variables V that are effectively overwritten by I k . The memory cost is proportional to p, whereas the extra run time cost comes from stack manipulation, usually minimal and proportional to p. Due to retrievals, the exit V does not contain the original result Y . These data-flow reversal strategies can be adapted to programs P that are no longer straight-line, but that use control. This implies that the control decisions of the original code (branch taken, number of iterations, . . . ) or equivalently of the forward sweep, must be made available in reverse for the backward sweep. In a Store-All context, this can be done just like for data: control-flow decisions must be recorded at the exit of the control structure (conditional, loop,. . . ) and then retrieved in the backward sweep to control the backward execution. This can use the same stack as data-flow reversal. The conjunction of the recorded control and data values is called the trajectory.

In practice, neither pure Recompute-All nor pure Store-All can be applied to large programs, because of their respective cost in run-time or memory space. Trade-offs are needed, and the classical trade-off is called checkpointing.

-In the Recompute-All approach, checkpointing means choosing a part C of P and storing the state upon exit from this part. Recomputing can then start from this state instead of the initial state. This is sketched on the top-left part of figure 3. At the cost of storing one extra state, the run-time penalty has been divided roughly by two.

Checkpoints can be nested to further reduce the run-time penatly, as shown on the bottom-left of figure 3. -In the Store-All approach, checkpointing means choosing a part C of P and not recording the overwritten values during C. Before the backward sweep reaches ← -C , C is run again from a stored state this time with recording. This is sketched on the top-right part of figure 3. At the cost of storing one extra state and of running C twice, the peak memory used to record overwritten data is divided roughly by two. Checkpoints can be nested to further reduce the peak memory consumption, as shown on the bottomright of figure 3.

Notice on the bottom of figure 3 that the execution scheme for nested checkpoints become very similar. Recompute-All and Store-All are the two ends of a spectrum, with optimal checkpointing scheme(s) lying somewhere in between. A good placement of (nested) checkpoints is crucial for efficient adjoint differentiation of large codes. Natural candidates to become a checkpointed part are procedure calls and loop bodies, but any piece of code with a single entry point and a single exit point can be chosen. There is no formula nor efficient algorithm to find this optimal placement of checkpoints, except in the case [START_REF] Griewank | Achieving logarithmic growth of temporal and spatial complexity in reverse Automatic Differentiation[END_REF]] of a loop with a known number of iterations all of the same cost. A good enough default strategy is to apply checkpointing at the level of each procedure call. In practice, it is important to give the end-user the freedom to place checkpoints by hand. Good placements of checkpoints perform similarly: memory and run-time costs grow with the logarithm of the run-time of P or more precisely -the peak memory size during execution, to store states and record data, grows like the logarithm of the run-time of P.

-the maximum number of times a checkpointed piece of the program is re-executed, which approximates the slowdown factor of P compared to P, also grows like the logarithm of the run-time of P. This explains why the slowdown ratio R a of P compared to P can be larger by a few units than the ratio R t of Ṗ compared to P.

Actually, one need not store the entire memory state to checkpoint a program piece C. What must be stored is called the snapshot. In a Store-All context, we can see that a variable need not be in the snapshot if it is not used by C. Likewise, a variable need not be in the snapshot if it is not overwritten between the initial execution of C and the execution of its adjoint C. This can be refined further, cf section 5.4.

DESIGN CHOICES RATIONALE

This section presents the design choices made when developing Tapenade. When appropriate, we give references to the underlying techniques and contrast with the choices made by other AD tools.

Source Transformation: as the main focus or our research is on the adjoint mode and how to keep it efficient on very large applications, we decided to stick with the program transformation approach. This approach was the choice made by the ancestor of Tapenade, Odyssée, developed at INRIA in the years 1990 [START_REF] Rostaing | Automatic Differentiation in Odyssée[END_REF]]. An efficient adjoint mode must completely remodel the control-flow structure of the code to achieve control-flow reversal. While Operator Overloading AD tools such as ADOL-C [START_REF] Walther | Getting started with ADOL-C[END_REF] need to store a trace of all operations made by the code to run it backwards, Program Transformation AD tools can store as little as a history of branches taken, of size in practice small compared to the tape of intermediate values computed [START_REF] Naumann | Control flow reversal for adjoint code generation[END_REF]. Of equal importance is the fact that global data-flow analysis such as TBR [Hascoët et al. 2005] allow the AD tool to reduce significantly the quantity of intermediate values that need to be stored on that tape. Let us mention two other interesting AD approaches that we have left aside. One way to build an AD tool that is truly independent from the source language is to apply AD to the assembly code. The problem is that there are also many different assembly languages. Moreover, compilation has wiped out structural information that can be used to optimize the AD code. Another way is to embed the AD process directly inside a compiler [START_REF] Naumann | Computing adjoints with the NAGWare Fortran 95 compiler[END_REF], in its "middle end" which is very often independent both from the source language and from the target assembly language. This can take advantage from the data-flow analyses already performed by the compiler, but these are not global due to the principle of separate compilation. Also implementation of the adjoint mode appears hard, due to a complex internal representation specialized for compilation.

Complete support of imperative languages: to experiment and validate refinements of adjoint AD, we must apply it to real, large application codes such as the ones found in industry. Also, a widespread use of the tool in an industrial context would be definitely a mark of the quality of its algorithms. For both reasons, we wanted the tool to be used on real-size applications, not asking for heavy preprocessing nor hand preparation. In other words, we didn't want to indulge restrictions of the tool with respect to the size of the application, nor to the language effectively accepted (ini-tially Fortran 77, then Fortran 90, then C). Accepting large applications essentially means choosing a compact internal representation and a careful design of data-flow analysis, such that the complexity of data-flow analyses remain close to linear with respect to code size and, depending on the analysis, linear or quadratic with respect to the number of variables. Accepting the full application language means treating correctly memory constructs such as global variables, Fortran COMMON's, EQUIVALENCE's, and control constructs such as GOTO's, that are sometimes neglected when only studying principles. This supposes a precise model of the internal memory, and the use of control-flow graphs (flow graphs for short) rather than mere syntax trees.

Readability: the output of the tool is a new source code, in which the end-user should be able to recognize the structure of the original code. This is useful not only for debugging purposes. This requires that the internal representation of codes keeps information that a compiler normally discards, and which is used during the source reconstruction phase. As a consequence, we don't want the tool to begin with a set of "canonicalization" steps, unless they can be reversed in the end.

Source language independence: obviously, the application language should not be only Fortran, like it was for Odyssée. The differentiation algorithms must not be expressed in terms of one given application language, but rather on an abstract language, called IL for Imperative Language, that is the sole input and output of the tool. IL must be rich enough to represent all imperative programming constructs, including Object Oriented constructs for future extension of the tool to these languages. Whenever possible, a unique IL operator must be used to represent similar constructs of different languages. This design choice, also made by the AD tool OpenAD [START_REF] Utke | OpenAD/F: A modular, open-source tool for Automatic Differentiation of Fortran codes[END_REF][START_REF] Utke | Separating language dependent and independent tasks for the semantic transformation of numerical programs[END_REF], allows a refinement in the AD model to benefit at the same time to applications written in Fortran or C.

Scoping: the internal representation, based on flow-graphs, must be rich enough to represent variable scoping. Each basic block of these graphs points to a symbol table, and these symbol tables are nested to implement scoping. The hierarchy of nested symbol tables must provide separation between public/formal parameters symbol table and private/local variables symbol table, and a mechanism to import the public symbol table of an imported module.

Arrays as atoms: every static data-flow analysis must face the problem of arrays and array indices. To distinguish data-flow information of different cells of arrays, one must be able to evaluate and compare each array indices in the code. This is out of reach of a static analysis in general, one reason being nondecidability. This can be achieved only in special cases such as regular array indices, but at a high development cost. While array index analysis is certainly vital for some static analysis tools such as loop-parallelizing or loop-vectorizing compilers, it is questionable whether it is profitable for an AD tool. Even if we discover that only a fraction of an array needs a derivative counterpart, leading to a smaller derivative array, the index computations needed to implement the correspondence between an array cell and its derivative cell may be very costly and will certainly hamper readability of the differentiated code. As a result, the design choice is to consider arrays as atoms and to accept the incurred information "blurring". If an array access A(i) modifies data-flow information, this will conservatively affect information on the complete array A. As only one cell of A is affected in reality, a "must" information on A(i) will result in a weaker "may" information on A. There are two cases where we will afford a little more accuracy on array data-flow analysis:

-when an array is split in sequential fragments by means of a Fortran EQUIVALENCE or COMMON, we will consider each fragment as a separate memory zone for data-flow analysis. The data-flow information may then be more accurate on each fragment. -when every access to an array inside a loop has the same index expression, and this index is a simple function of the loop indices, e.g. with a constant offset, data-flow analysis will gain accuracy by temporarily individualizing array elements.

Distinguished structure components: In contrast to array elements, components of structures must be separated. Whereas elements of an array often share a similar behavior, we believe that different components of a structure probably behave differently. Therefore, data-flow information on one structure component must not affect information on the other components of the same structure.

Association by name: A memory reference, e.g. a[i].x, must be associated in some way to its differentiated counterpart. We choose association by name, i.e. to create new variables with new names. On our example, this means creating a new ad and referring to ad[i].x. In contrast, we don't choose association by address, i.e. to attach the derivative to the address of each reference by extending the type of atomic components. On our example, this would mean extending the type of the x component and referring to a[i].x.derivative. AD tools that rely on overloading naturally use association by address. Our choice strongly impacts the differentiated code [START_REF] Fagan | Data representation alternatives in semantically augmented numerical models[END_REF]]. In some cases such as pointers, it may appear awkward. Future versions of Tapenade will probably offer both options.

Context-sensitive and Flow-sensitive Data-Flow analysis: Accurate data-flow analysis must be flow-sensitive and context-sensitive [START_REF] Shin | Design and implementation of a Context-Sensitive, Flowsensitive activity analysis algorithm for Automatic Differentiation[END_REF], which is easier to achieve on a graph-based representation i.e. a call graph of flow graphs. The analyses perform a number of fixed-point sweeps on these graphs. However, different calling contexts may lead to different differentiated procedures, and there is a risk of combinatorial explosion. To avoid that, our choice is generalization instead of specialization: at the time of differentiation, only one differentiated subroutine is built for each original subroutine. In other words, the analyses implement context-sensitivity, but its effects at differentiation time are restricted.

Store-All in adjoint mode, Save-on-kill: For the adjoint mode, we choose the Store-All approach to recover intermediate values. Intermediate values are saved just before a statement may overwrite them (Save-on-kill). We don't consider the alternative choice of storing the partial derivatives themselves (preaccumulation). Ope-nAD [START_REF] Utke | OpenAD/F: A modular, open-source tool for Automatic Differentiation of Fortran codes[END_REF]] is also basically Store-All, unlike TAF [START_REF] Giering | Tangent linear and Adjoint Model Compiler[END_REF]] which is basically Recompute-All. This contrast is reduced by the use of checkpointing as shown in figure 3. Also, Tapenade applies when possible a limited form of local recomputing instead of storing.

THE TAPENADE AUTOMATIC DIFFERENTIATION MODEL

This section intends to give the reader familiarity with the tangent and adjoint codes produced by Tapenade, going from simple creation of differentiated variable names to the more sophisticated improvements of adjoint AD. Unless otherwise stated, the illustration codes use the Fortran syntax but they easily transpose to C.

Symbol names

If a variable v has a non-trivial derivative, this derivative must be stored. As Tapenade uses association by name, we need a new variable name which is built by adding a user-defined prefix and/or suffix to v. By default in tangent mode the new variable will be vd ("v dot"), and in adjoint mode vb ("v bar"). Similarly, Tapenade AD model introduces derivatives of COMMONS, procedures, interfaces, user-defined types, and modules, and must build new names for them. By default this will be done by appending " d" in tangent mode and " b" in adjoint mode. Figure 4 shows the transformation of the declaration part of a procedure. If a conflict is found with some existing symbol, then a "0" is appended, then "1" until finding a free symbol.

Types of differentiated variables

When a variable of type T has a non-trivial derivative, Tapenade builds a differentiated variable of a "differentiated" type T . For elementary types T of discrete nature such as booleans, integers, characters, strings. . . , there is no T . For the other elementary types such as REAL, double, COMPLEX. . . (possibly with modifiers), T = T . For arrays of element type T with non-trivial derivatives, the differentiated array has element type T and the same dimensions number and size. Likewise, pointers to variables of type T with non-trivial derivatives are differentiated into pointers to variables of type T . When T is a structured type1 , the differentiated type T is a new structured type deduced from T by erasing all components that never have a non-trivial derivative. Notice the choice of generalization vs specialization: if somewhere in the code, component x of a variable of type T has a non-trivial derivative, then type T will have a component x. In contrast, specialization would create several differentiated types for T , for a benefit that we judge negligible. We will repeat this choice of generalization for subroutines and functions (cf section 4.5). For instance in the tangent differentiation of figure 5, structured type VECTOR is differentiated into VECTOR D, and type components y and name do not appear in VECTOR D because there are never any derivatives to store there. Component x appears because ad needs it, z appears because ud(3) needs it.

Please notice the way references to structure components are differentiated. With our choice of association by name, the derivative of component u(3)%z is stored in ud(3)%z, i.e. the z component of the derivative name associated to u(3). If we chose association by address, the VECTOR D type would contain components x, y, z, and name, where x and z would be of type, say, ACTIVE REAL instead of REAL. Access to the primitive value would become u(3)%z%value and to the derivative u(3)%z%derivative.

Simple assignments

Consider an assignment I k . In tangent mode (cf equation ( 1)), derivative instruction

İk implements Ẋk = f k (X k-1 ) × Ẋk-1 , with initial Ẋ0 = Ẋ. In adjoint mode (cf equa- tion (2)), derivative instruction(s) ← - I k implements X k-1 = X k × f k (X k-1 ), with initial X p = Y .
In practice, the original code may overwrite variables, i.e. implement X k = f k (X k-1 ) through an assignment I k that puts the resulting X k into the memory locations that used to contain X k-1 . Tapenade will do the same for the derivative assignments: the old Ẋk-1 will be overwritten by the new Ẋk , and the old X k will be overwritten by the new X k-1 . Suppose for instance I k is a(i)=x*b(j) + COS(a(i)), with a, b, and x having nontrivial derivatives (cf section 4.6). It is straightforward to write the Jacobian of the corresponding function

f k : IR 3 → IR 3 a(i) b(j) x → x * b(j) + COS(a(i)) b(j) x which is the matrix -SIN(a(i)) x b(j) 0 1 0 0 0 1
from which we write the operations that İk and ← -I k must implement:

İk implements   ȧ(i) ḃ(j) ẋ   := -SIN(a(i)) x b(j) 0 1 0 0 0 1 ×   ȧ(i) ḃ(j) ẋ   , ← - I k implements a(i) b(j) x := a(i) b(j) x × -SIN(a(i)) x b(j) 0 1 0 0 0 1
, Tapenade thus produces the derivative instructions shown on figure 6.

Tapenade tangent Tapenade adjoint ad

(i) = xd*b(j) + x*bd(j) -ad(i)*SIN(a(i)) xb = xb + b(j)*ab(i) bb(j) = bb(j) + x*ab(i) ab(i) = -SIN(a(i))*ab(i)
Fig. 6. Differentiation of a simple assignment

Sequences of instructions and control

The expressions from the Jacobian, e.g. SIN(a(i)) in section 4.3 refer to values found in X k-1 . The differentiation model must ensure that these values are available when the derivative instructions are executed.

In tangent mode, this is achieved by putting the derivative instructions of I k just before a copy of I k . Thus the tangent differentiated code for a procedure P is basically a copy of P with derivative instructions inserted just before each assignment statement. Procedure calls are treated differently, cf section 4.5.

In Adjoint mode, our design choice is fundamentally the store-all strategy (cf section 3). Therefore control-flow reversal uses a two-sweeps architecture. The first sweep (the forward sweep) is basically a copy of P, augmented with recording of the control. This recorded control is used by the second sweep (the backward sweep) to orchestrate control-flow reversal. The forward sweep also records intermediate data values (the X k ), that will be used by the backward sweep to evaluate elements of the Jacobian.

As both control-and data-flow are reversed, the natural way to record is to use a stack that grows during the forward sweep and shrinks during the backward sweep. Because of this stack architecture, data values must be recorded after their last use forward, so as to be restored before their first use in the backward sweep. The simplest choice to achieve this (Save-on-kill), saves a data value immediately before an instruction overwrites it. Concerning control, the analogue of Save-on-kill consists in saving control values just before control-flow merges, rather than when it splits just after the forward control tests. One advantage is that this choice applies naturally to any flow graph, even not cleanly structured. At the same time, the reversed control-flow retains most of the structure of the original code. The storage strategies for control-flow and for data-flow are compatible, which means one stack is enough to store both. Figure 7 shows the Tapenade tangent and adjoint AD models on a representative example.

Procedures

We saw that differentiation operates at the level of an atomic instruction, e.g. an assignment, by adding to it one or a few new, differentiated instruction(s). Differentiation at the level of a procedure is different. In tangent mode, the differentiated procedure contains at the same time the original computation and its derivative instructions, and calls to the procedure are replaced by calls to the differentiated procedure. In adjoint mode, the original computation and the derivative computation are split between the forward and backward sweeps, thus creating two differentiated procedures, one for each sweep. Therefore in Tapenade's basic adjoint model, each original call to a procedure FOO gives birth to a call to some FOO FWD in the forward sweep, and to FOO BWD in the backward sweep. Notice however that this basic adjoint model is not Tapenade's default behavior. The default behavior applies checkpointing (cf section 4.10) to each call, yielding a different structure. The basic model is effectively what one gets when checkpointing is turned off.

The arguments of the differentiated calls are deduced naturally from the arguments of the original call. Arguments that have a nontrivial derivative get followed by an additional argument for this derivative. Our choice of generalization vs specialization ... bb = 0.0 CALL FOO_BWD(c,z,zb,t,tb,bb) CALL FOO_BWD(a,x,xb,2.5,dummb,yb) Fig. 8. Tangent and adjoint differentiation of procedure calls. Assume only x, y, z, and t have a nontrivial derivative has a consequence here again: we want only one differentiated version of each procedure, which must therefore accomodate the "envelope" of all possible call situations regarding the existence of a nontrivial derivative. In other words the differentiated procedure may expect a derivative for some formal argument while one particular call has no nontrivial derivative for this actual argument. In this case, this particular call must provide a dummy derivative argument which can be an immediate 0.0 or sometimes a placeholder argument that may have to be initialized with 0.0's.

Figure 8 illustrates this on the representative case of a function call. Assume only x, y, z, and t have a nontrivial derivative, so that only arguments 2 and 3 and result of FOO have nontrivial derivatives at some call site. This determines the arguments expected by the derivative FOO D and FOO BWD. The original result y becomes an extra output argument of FOO D. The adjoint mode FOO FWD has no extra argument as the forward sweep doesn't manipulate derivatives. The original result y is needless in the backward sweep and disappears from FOO BWD, which need not be a function any more. The adjoint data-flow reversal interacts with the parameter-passing discipline of the application language. This is transparent on a language that uses call by reference, or call by value-return like Fortran, but becomes apparent with call by value as used in C. An input argument of the original code can be passed by value, but the corresponding differentiated argument becomes an output of the adjoint differentiated procedure and must be passed by reference.

Activity

The past sections compose a naïve AD model that effectively computes the expected derivatives. It is the simplest AD model that Tapenade can actually obey, using a combination of command-line options that deactivate all further refinements. However, there are a number of efficiency improvements that can and should be done. They are the subject of the coming sections.

In the previous sections we used the notion of a nontrivial derivative. At a given location in the source, a variable has a "trivial" derivative if we know statically that this derivative is always null, or is useless for the desired final derivatives. Otherwise the variable has a nontrivial derivative and is said active. One could consider that all variables of a "differentiable" type (e.g. float, complex. . . ) are active everywhere. Obviously one can do much better with the help of a static data-flow analysis called activity analysis. In case of doubt, one must assume conservatively that a variable is active. Thus, active means "possibly active". Tapenade model does not distinguish array elements: if one array element becomes active, so does the complete array.

One important ingredient of activity analysis is independent and dependent arguments. Differentiation operates on an algorithm coded as a part (P in section 2) of a bigger program. The AD user should designate differentiation candidate arguments of P, according to the call and usage context of the differentiated code:

independent variables are the subset of the inputs that must be differentiated. By designating them, the end-user guarantees the context will provide their differentiated variables to P D in tangent mode, whereas in adjoint mode the context has the right to use the values found in their differentiated variables upon return from P B. -dependent variables are the subset of the outputs that must be differentiated. By designating them, the end-user guarantees the context will provide their differentiated variables to P B in adjoint mode, whereas in tangent mode the context has the right to use the values found in their differentiated variables upon return from P D.

Anticipating on section 5.2, we see that active variables are indeed variables that are at the same time varied, i.e. maybe influenced in a differentiable way by some independent, and useful, i.e. maybe influencing some dependent in a differentiable way.

-In tangent mode when v is not varied, or in adjoint mode when v is not useful, we know the current derivative of v is certainly null. The differentiated code is therefore simplified by forward partial evaluation. -In tangent mode when v is not useful, or in adjoint mode when v is not varied, we know the current derivative of v is not needed for the rest of the derivative computation. The differentiated code is therefore simplified by backward slicing.

Figure 9 illustrates the benefits of activity on the tangent AD model.

To further reduce the cost of the differentiated instructions, we decided that derivatives that become null need not be immediately reset to 0.0. We say that they are implicit-zero. This happened once in figure 9. This spares the reset statement and can be compensated by the next instruction that assigns to the derivative variable. However, the model must pay special attention when the control flow merges a branch in which the variable is implicit-zero with another branch in which the variable is active. The implicit-zero differentiated variable must be explicitly reset to 0.0 before the flow merges. All in all, this implicit-zero tactique yields a limited gain, at the cost of extra complexity in the code produced. It can be turned off by a command-line option.

Vector mode

It happens that derivatives are required for the same values of the input parameters, i.e. at the same point X in the input space, but for several sets of values for the derivatives Ẋ or Y . For example in section 2, we saw that the full Jacobian Matrix at a given point X can be obtained by repeatedly running the tangent code for each Ẋ in the Cartesian basis of the input space, or by repeatedly running the adjoint code for each Y in the Cartesian basis of the output space. At the same time, we observe that there are two sorts of instructions in a differentiated code, regardless tangent or adjoint:

-some instructions are in fact copies of the original code, and do not operate on derivative variables. -some instructions actually compute the derivatives, using values both from the original code and from other derivatives computation

The copies of the original code operate on the same values and return the same results regardless of the derivatives. Therefore they can be factored out. To this end, the socalled vector extension of Tapenade produces a code where all instructions that operate on derivatives are encapsulated into an extra loop on all different values of Ẋ.

Figure 10 illustrates vector tangent mode on a small procedure. All derivative variables receive an additional dimension, that represents the different directions Ẋ in the input space. In the case of a variable of a structured type, or of array type, the additional dimension is placed at the deepest location in the data structure, so that taking a part of the variable can easily translate into taking the corresponding part of the derivative.

Tapenade makes a data-dependency analysis on the differentiated code that allows differentiated instructions to safely move about the code. In particular, this allows the tool to group loops on differentiation directions when possible, and then merge them into a single loop. This reduces loop overhead, see figure 10. More importantly since these loops are essentially parallel, this builds larger parallel sections e.g. for a multithreaded execution context like OpenMP.

At present, there is no equivalent vector extension for the adjoint mode. However this extension could be built provided some sizeable application has a need for it.

Splitting and merging the differentiated instructions

From this section on, the features described concern only the adjoint mode of AD. These features are probably more complex and demanding regarding implementation. On the other hand, Tapenade clearly focuses on an efficient adjoint mode, and the benefits of this mode are certainly worth the effort.

The adjoint differentiation model applies elementary transformations to the differentiated code, such as splitting and merging instructions, to eliminate redundant operations and improve performance. Instructions merging relies on Tapenade running data-dependency analysis on the basic blocks of the differentiated code. These transformations remain localized to each basic block of the control-flow graph. Although they are local to each basic block, these transformations are not redundant with compiler optimizations because basic blocks of numerical programs are often larger than the peep-hole optimization window of the compiler. Also, the compiler doesn't know the special structure of an adjoint code. Let us consider as an example the following basic block:

res = (tau-w(i,j))*g(i,j)*(z(j)-2.0)/v(j) g(i,j) = res*v(j)
Differentiation of complex expressions often introduces common subexpressions. Tapenade does not look for common sub-expressions in the original code, but at least it should not introduce new ones. In the context of AD, the problem of common subexpression detection is simpler than in the general case because we know how derivative statements are structured and built, and therefore when new common subexpressions are introduced. On our current example, assuming that res, and only res, is nonactive outside the basic block, the adjoint statements of the naïve backward sweep could be as shown on the left of figure 11. However, from the structure of the deriva-Tapenade A:19 naïve adjoint split and merge resb = v(j)*gb(i, j) vb(j) = vb(j) + res*gb(i, j) gb(i, j) = 0.0 taub = taub +(z(j)-2.0)*g(i, j)*resb/v(j) wb(i, j) = wb(i, j) -g(i, j)*(z(j)-2.0)*resb/v(j) gb(i, j) = gb(i, j) +(z(j)-2.0)*(tau-w(i, j))*resb/v(j) zb(j) = zb(j) +(tau-w(i, j))*g(i, j)*resb/v(j) vb(j) = vb(j) -(tau-w(i, j))*g(i, j)*(z(j)-2.0)*resb/v(j)**2 resb = v(j)*gb(i, j) temp = (z(j)-2.0)/v(j) tempb0 = temp*g(i, j)*resb tempb = (tau-w(i, j)) *g(i, j)*resb/v(j) vb(j) = vb(j) +res*gb(i, j) -temp*tempb gb(i, j) = temp *(tau-w(i, j))*resb taub = taub + tempb0 wb(i, j) = wb(i, j) -tempb0 zb(j) = zb(j) + tempb Fig. 11. Splitting and merging adjoint differentiated instructions tive statements, Tapenade can find that three subexpressions are duplicated and are worth being precomputed namely, (z(j)-2.0)/v(j), g(i, j)*(z(j)-2.0)*resb/v(j), and (tau-w(i, j))*g(i, j)*resb/v(j). It is also apparent that the two statements incrementing vb(j) can be merged, and that one statement resets gb(i, j) to zero before it is later incremented. The right of figure 11 shows the differentiated code after all this splitting and merging.

Improving trajectory computation

The adjoint differentiation model of Tapenade applies several improvements to eliminate unnecessary instructions from the trajectory computation done by the forward sweep. This section describes three such improvements, that must be applied in order. We will illustrate this on the following piece of code, considered as the complete body of the differentiated procedure, so that its adjoint forward sweep is immediately followed by its backward sweep:

n = ind1(i) b(n) = SIN(a(n))-b(n) a(n) = a(n)+x c = a(n)*b(n) a(n) = a(n)*a(n+1) n = ind2(i+2) z(n) = z(n) + c
Figure 12 shows the successive adjoint codes after each improvement is performed. Diff-liveness is about removing original instructions of the forward sweep that were used for the original result but are not needed for any derivative. Since Tapenade considers that the only important result is now the derivatives, the original result is unimportant and we can perform slicing, backwards from the end of the forward sweep, to take away original computation not needed for derivatives.

At first sight, this slicing is of limited impact as it applies only "near" the end of the forward sweep. Moreover, some end-users require that the adjoint code also returns the original result, and therefore this slicing may need to be deactivated by a command-line option. However, the original result is definitely useless for every part of the code where checkpointing is applied (cf section 4.10). In this case this slicing is done systematically.

In the example of figure 12, we can see that the values of z(n) and c are not used in the derivatives computation. Therefore the corresponding assignments are not diff- 12. Diff-liveness, TBR, and recomputation improvements applied successively to a small example live, i.e. they are dead code for the adjoint. This also spares the PUSH/POP statements that save into the trajectory the values that are overwritten by these assignments.

naïve Diff-liveness TBR recomputation CALL PUSHINTEGER4(n) n = ind1(i) CALL PUSHREAL4(b(n)) b(n)=SIN(a(n))-b(n) CALL PUSHREAL4(a(n)) a(n) = a(n) + x CALL PUSHREAL4(c) c = a(n)*b(n) CALL PUSHREAL4(a(n)) a(n) = a(n)*a(n+1) CALL PUSHINTEGER4(n) n = ind2(i+2) CALL PUSHREAL4(z(n)) z(n) = z(n) + c CALL POPREAL4(z(n)) cb = zb(n) CALL POPINTEGER4(n) CALL POPREAL4(a(n)) ab(n+1) = ab(n+1) +a(n)*ab(n) ab(n) = b(n)*cb +a(n+1)*ab(n) CALL POPREAL4(c) bb(n) = bb(n) +a(n)*cb CALL POPREAL4(a(n)) xb = xb + ab(n) CALL POPREAL4(b(n)) ab(n) = ab(n) +COS(a(n))*bb(n) bb(n) = -bb(n) CALL POPINTEGER4(n) CALL PUSHINTEGER4(n) n = ind1(i) CALL PUSHREAL4(b(n)) b(n)=SIN(a(n))-b(n) CALL PUSHREAL4(a(n)) a(n) = a(n) + x CALL PUSHINTEGER4(n) n = ind2(i+2) cb = zb(n) CALL POPINTEGER4(n) ab(n+1) = ab(n+1) +a(n)*ab(n) ab(n) = b(n)*cb +a(n+1)*ab(n) bb(n) = bb(n) +a(n)*cb CALL POPREAL4(a(n)) xb = xb + ab(n) CALL POPREAL4(b(n)) ab(n) = ab(n) +COS(a(n))*bb(n) bb(n) = -bb(n) CALL POPINTEGER4(n) n = ind1(i) b(n)=SIN(a(n))-b(n) CALL PUSHREAL4(a(n)) a(n) = a(n) + x CALL PUSHINTEGER4(n) n = ind2(i+2) cb = zb(n) CALL POPINTEGER4(n) ab(n+1) = ab(n+1) +a(n)*ab(n) ab(n) = b(n)*cb +a(n+1)*ab(n) bb(n) = bb(n) +a(n)*cb CALL POPREAL4(a(n)) xb = xb + ab(n) ab(n) = ab(n) +COS(a(n))*bb(n) bb(n) = -bb(n) n = ind1(i) b(n)=SIN(a(n))-b(n) CALL PUSHREAL4(a(n)) a(n) = a(n) + x n = ind2(i+2) cb = zb(n) n = ind1(i) ab(n+1) = ab(n+1) +a(n)*ab(n) ab(n) = b(n)*cb +a(n+1)*ab(n) bb(n) = bb(n) +a(n)*cb CALL POPREAL4(a(n)) xb = xb + ab(n) ab(n) = ab(n) +COS(a(n))*bb(n) bb(n) = -bb(n) Fig.
To-Be-Recorded analysis (TBR) [Hascoët et al. 2005] deals with original values which, although they are diff-live, do not appear in any derivative. For instance if a variable is always used in linear expressions, it never appears in the derivative computations. In the example, this is the case for the initial n and b(n).

Finally, recomputation analysis looks for variables that are recorded in the trajectory, and that could be recomputed instead because their value is the result of a relatively cheap assignment and all inputs to this assignment are still available at the time of the POP. In this case, it is preferable to recompute the value by duplicating the assignment, rather than using precious memory space to record it. In the example, this is the case for the assignments to index n.

Checkpointing

Checkpointing (cf section 2) is the standard way to organize a tradeoff between run-time and memory consumption in adjoint AD. In our context that is based on the store-all strategy, checkpointing applied to a fragment C of an original code P=U;C;D amounts to replacing the normal two-sweep adjoint of P Tapenade A:21

P = - → P ; ← - P = - → U ; - → C ; - → D ; ← - D ; ← - C ; ← - U ; by P = - → U ; •; C; - → D ; ← - D ; •; - → C ; ← - C ; ← - U ,
where • and • respectively store and restore a snapshot i.e. "enough" values to guarantee that the execution context is the same for runs of C and of -→ C ; ← -C The gain is that the trajectory for C and for D are never stored at the same time in memory. The cost is the relatively small memory size of the snapshot, and the time for the extra run of C.

By default, the Tapenade adjoint AD model applies checkpointing to every procedure call. Therefore checkpoints are nested just like calls are nested. Checkpointing changes radically the differentiated code for a procedure call, as shown in figure 13. When the When checkpointing is applied to C, the adjoint code contains the sequence -→ C ; ← -C , in addition to the sequence -→ D ; ← -D already present before checkpointing is applied. Recalling that the code improvements such as diff-liveness (cf section 4.9) bear more fruit when the forward sweep is immediately followed by the backward sweep, we see that checkpointing opens more opportunity to simplify the adjoint code. In any case, when checkpointing is applied, an AD tool should restart the static data-flow analyses for the checkpointed part, in order to discover these additional simplifications.

The default checkpointing strategy of Tapenade can be amended in several ways. The directive $AD NOCHECKPOINT attached to a procedure call deactivates checkpointing on this particular call. To deactivate checkpointing on all calls, attach the directive to the procedure header itself. Conversely, the pair of directives $AD CHECKPOINT-START and $AD CHECKPOINT-END define an arbitrary portion of code (well-structured) on which checkpointing must be applied.

Concerning loops, directive $AD BINOMIAL-CKP attached to the loop header tells Tapenade to apply binomial checkpointing as defined in [START_REF] Griewank | Achieving logarithmic growth of temporal and spatial complexity in reverse Automatic Differentiation[END_REF]], which is the optimial strategy to apply to iterative loops such as time-stepping loops. Finally, directive $AD II-LOOP attached to the loop header tells Tapenade that the loop has Independent Iterations. Checkpointing can be therefore applied to each iteration of the loop, thus reducing the peak memory used for trajectory storage [Hascoët et al. 2002].

Adjoint of array notation

The adjoint model applies nicely to Fortran array notation. It exhibits a natural duality between the sum reductions and the spread operations that turn a scalar into an array. Masks and where guards are exchanged or preserved by the transformation. Figure 14 presents a few examples, which can be checked manually by writing the instruction's Jacobian like in section 4.3.

SPECIFICATION AND NOTES ON ACTUAL IMPLEMENTATION

We will give the specification of Tapenade, from its overall architecture down to the individual analyses and differentiation transformation. The architecture is fundamental, and is therefore described in every necessary detail. It corresponds exactly to the actual implementation. The analyses and differentiation are specified in a more abstract way, namely with data-flow equations and inference rules of Operational Semantics. They can be used for formal reasoning on the AD model. Their link with the actual implementation is discussed in the end of the section.

Internal Representation

Figure 15 summarizes the architecture of Tapenade, motivated by our design choices.

On its left, this architecture resembles that of a compiler, building an internal rep-

C parser

Fortran95 parser

Fortran77 parser

IL Tree

Internal representation builder

Internal Representation

Data-Flow analysis

Data-Flow info Differentiation

Internal Representation

IL Tree rebuilder

IL Tree

C printer

Fortran95 printer

Fortran77 printer

Fig. 15. Overall architecture of Tapenade resentation and running data-flow analysis on it. A big difference, though, is that ef-AD analyses must be global and therefore there is no such thing as separate compilation: all the source code must be parsed, then analyzed jointly. On the right, figure 15 differs from a compiler in that Tapenade produces its result back in the original programming language instead of as machine code. This imposes additional constraints to keep the result readable, and to maintain some degree of resemblance between source and differentiated codes. For instance unlike a compiler, the internal representation will store the order and arrangement of declarations inside a procedure. Although this order is useless for the data-flow analysis and differentiation phases, as declarations are compiled into symbol tables, it is used to regenerate declarations for the differentiated code in a matching order. The same applies to comments.

Language independence is enforced by delegating parsing and unparsing to separate language-specific modules, each of which communicate with the central part through a generic Imperative Language (IL). IL is an abstract language, with no concrete syntax, i.e. no textual representation. Parsing any source code produces an IL tree, whose nodes are IL operators. IL provides enough operators for all constructs of the accepted source languages, such as procedure definition, variable declaration, assignment, procedure call, variable reference. . . Overlapping is promoted: at present, 71 IL operators are common to Fortran and C, 55 are Fortran-specific and 26 are C-specific. In addition, 19 operators exist in provision for Object languages. Promotion of overlapping sometimes implies (reversible) tree transformation, such as using the more explicit C-style pointer access methods or using a generalized form of memory allocation.

However, the central part of Tapenade i.e. all data-flow analysis and actual differentiation, operate on an internal representation of programs which is more elaborate than plain IL trees. This internal representation is sketched in figure 16, and its motivation follows. As we chose to make data-flow analysis context-sensitive, the internal representation must provide easy access to the call graph. The call graph is on the left part of figure 16. The example illustrates the 3 kinds of arrows: Top contains M3 which contains P5, M5 imports (or uses) M3, P7 calls P5. In fact, the top level of each analysis consists of a sweep on the call graph, visiting each procedure possibly repeatedly until reaching a fixed-point. This also proves useful for differentiation itself, as the differentiated code is built procedure per procedure.

At the level of an individual procedure, the internal representation must essentially capture the flow of control. The natural representation is therefore a classical flow graph, whose nodes are Basic Blocks (sequential lists of instructions) and whose arrows are jumps from a Basic Block to another. In figure 16, all the flow graphs are gathered in the middle part. Notice that the flow graph is reduced to a single Basic Block for declarative units such as modules. When several arrows flow from one Basic Block, the instruction at the end of this Basic Block is the control instruction that decides the arrow taken at run-time. This uniform representation of control flow makes it easier for all data-flow analysis to sweep through the procedure code in a flow-sensitive way. It also makes it easier for the differentiation tool to generate a differentiated flow graph, leaving the details of how to create a code that implements this differentiated flow graph to the IL Tree rebuilder. In theory the declaration statements could be taken out of the flow graph, since they are not "executed". However we will keep them because their order is important. Moreover, keeping the original form and order of declarations enhances readability of the differentiated code.

Efficient access from a symbol name (variable, procedure, type . . . ) to its relevant information is ensured through Symbol Tables, sketched on the right part of figure 16. In particular data-flow analysis requires that each procedure call provides efficient access to the Call Graph node(s) of the called procedure, which can be several when overloading comes into play. This is done through the Symbol Table of the calling location. Symbol Tables use a hash code algorithm. Symbol Tables are organized as a tree to implement scoping. A Symbol Table can be copied elsewhere in the tree, with possible renaming, when another program element "imports" or "uses" it. Each Symbol contains in addition the corresponding node of the Call Graph, itself related to each of its callers and callees. This node in turn has a link to its flow graph.

The complete network of all these links as illustrated by figure 16 is not so intricate after all. This is about all we need to implement the data-flow analysis and program transformations of Tapenade.

At the level of an individual statement, the internal representation can remain the IL Tree itself, as the very simple control in this statement is adequately represented by this tree. However, when an instruction of the original code incorporates some degree of control such as an if-expression or a function call nested in an expression, then the internal representation must split this instruction into simpler instructions. In practice, this is done selectively and creates far less splitting than a systematic normalization step.

When a symbol is associated with storage (a variable, but also in the future an object or a class), it is given a (bunch of) unique integer number that represent its memory location(s). This allows all data-flow analysis to manipulate bitsets, based on these Tapenade A:25 integer numbers, called zones. A scalar or an array of scalars receive only one zone. On the other hand, a structured object receives one zone per "component" or "field". We believe this greatly contributes to efficiency of the data-flow analysis.

Data-Flow equations of program static analysis

Static data-flow analysis [START_REF] Aho | Compilers: Principles, Techniques and Tools[END_REF]] is an essential step to achieve efficient differentiation. In other words, almost every improvement to our Automatic Differentiation model required a new specific data-flow analysis in Tapenade. Figure 17 shows them.

Useful inputs Analysis

Result

LEGEND :

... used in all analyses below:

Source Internal Representation

Pointer analysis

Pointer targets

In-Out analysis

In-Out sets

Detection of Uninitialized uses

Warning messages

Detection of Aliasing

Warning messages

Diff-dependency analysis

Diff-dep relations Dependency analysis

Dep relations Activity analysis (Varied and Useful)

Active variables

Diff-liveness analysis

Diff-live variables

TBR analysis

TBR variables

Adjoint-out analysis

Adj-out variables

Tangent AD

Tangent Internal Representation

Adjoint AD Regarding notation, when a data-flow information Info is defined for a given location in the code, i.e. in our case before or after an arbitrary instruction I, we will write Info -(I) (resp. Info + (I)) the value of Info immediately before (resp. after) I.

Adjoint Internal Representation

To begin with, we need to run two classical, AD-independent data-flow analyses that are needed almost everywhere:

-Pointer analysis finds pointer destinations, which tell for each instruction and for each available pointer, the smallest possible set of possible targets of this pointer. This pointer analysis is not specific to AD and is outside the scope of this paper. -In-out analysis finds in-out sets, which contain for each procedure the smallest possible sets of variables that may or must be read or overwritten by this procedure. This composite piece of information is best represented by four booleans per variable. Each of them is true when, during some execution of the procedure, a part of the input value of the variable may be -read and then overwritten, -read and not overwritten, -not used but overwritten, -neither used nor overwritten.

The in-out sets also include, for each instruction I in each procedure, the sets Cst(I) (resp. Dead(I)) of variables that certainly remain constant (resp. unused) from this location till the exit from the procedure. As a side effect, in-out analysis makes it possible to detect potential aliasing and potential use of uninitialized variables, yielding warning messages especially valuable in the context of AD.

There is an issue of relative dependence between the two analyses. In some cases, pointer destinations analysis could benefit from the in-out sets, but the influence of pointer destinations on in-out analysis is far more important and dictates the order in which the analyses are run. Pointer destinations and in-out sets are necessary to almost all subsequent analyses in Tapenade.

Then come the AD-specific analyses. They are needed to implement the differentiation model described in sec. 4. Their relative dependences are shown in figure 17. We start with activity analysis, that allows for the improvements discussed in sec. 4.6. Then come diff-liveness and TBR analyses, allowing for the improvements discussed in sec. 4.9, and finally comes adjoint-out analysis, which is one ingredient of the snapshots used when checkpointing, cf section 4.10. For efficiency reasons we introduce two preliminary analyses which, for each procedure, compute useful summarized information used in activity, diff-liveness, and adjoint-out analyses: dependency analysis and diff-dependency analysis. TBR and adjoint-out analyses are needed only by the adjoint mode of AD. We will now give the data-flow equations of these AD-specific analyses.

-Activity analysis is the combination of a forward and a backward analysis. It propagates forward the Varied set of the variables that depend in a differentiable way on some independent input. Similarly, it propagates backwards the Useful set of the variables that influence some dependent output in a differentiable way. Since the relation "depends in a differentiable way of " is transitive on code sequences, the essential equations of the propagation are: Varied + (J)

Useful + (I) = J successor of I Useful -(J)
The final activity is just the intersection:

Active -(I) = Varied -(I) ∩ Useful -(I) Active + (I) = Varied + (I) ∩ Useful + (I)
-Diff-dependency analysis is the preliminary computation, for activity analysis, of the diff-dependency Diff-dep(P ) for each procedure P . If we define the composition of Diff-dep objects as: -Diff-liveness analysis finds instructions of the original program that, although perfectly useful to compute its results, are not needed to compute its derivatives. If only the derivatives are of interest in the end, then these instructions can be removed from the differentiated code. This is therefore a special case of dead code detection. Diff-liveness analysis propagates, backwards on the flow graph, the Diff-live set of variables that are needed to ultimately compute a derivative. At the tail of the differentiated code, it is initialized to ∅. A variable v is diff-live before an instruction I either if v is used in the derivative instruction I , or if some variable that depends on v is diff-live after I. Hence the data-flow equations:

Diff-live -(I ∞ ) = ∅ Diff-live -(I) = use(I ) ∪ (Dep(I) ⊗ Diff-live + (I))
that use Dep(I), the dependency relation between variables before I and after I. Like Diff-dep(I), Dep(I) is precomputed for efficiency when I is a call to a procedure P . We just need to add the closure data-flow equation for backward propagation:

Diff-live + (I) = J successor of I Diff-live -(J)
Diff-liveness must be extended to control: a statement that has a derivative or that computes a diff-live variable makes its surrounding control diff-live. A diff-live control makes all the variables it uses diff-live. -Dependency analysis is the preliminary computation of the dependency between inputs and outputs of each procedure P . We need the same composition operation ⊗ as was needed for Diff-dependency analysis. Only the elementary dependencies through a simple instruction, e.g. an assignment, is different: the left-hand side of the assignment depends on all variables occurring in the right-hand-side, including indices, and also on the indices occurring in the left-hand side itself. We can thus compute Dep(P ) by a forward sweep on P , on the model of Diff-dependency analysis:

Dep + (I 0 ) = Id Dep + (I) = Dep -(I) ⊗ Dep(I) Dep -(I) = J predecessor of I Dep + (J)
and the result

Dep(P ) is found in Dep -(I ∞ ).
-TBR analysis is specific to the adjoint mode of AD when using the Store-All strategy. It finds variables whose present value is used in a derivative instruction. This property is dynamic: it is reset to false when the variable is completely overwritten with a new value. The goal of TBR analysis is to find locations where a variable is overwritten, and we know that the overwritten value is needed for a derivative instruction. At such a location, our adjoint AD model (store-all) will insert storage of the value before it gets overwritten in the forward sweep, and will restore it at the corresponding location in the backward sweep. TBR analysis propagates, forward on the flow graph, the TBR set of variables whose present value is (actually may be) used in a derivative instruction. Initially, this set is empty at the beginning of the code to differentiate. A variable v enters the TBR set after instruction I if v is used in the derivative instruction I , and is removed from this set if I is present in the differentiated program and v is certainly and completely overwritten by I (v ∈ kill(I)). This refinement about I being present (or live) in the differentiated code comes directly from diff-liveness analysis:

isLiveForDiff(I) ⇐⇒ out(I) ∩ Diff-live + (I) = ∅
Hence the data-flow equations:

TBR + (I) = (TBR -(I) ∪ use(I )) \ kill(I) if isLiveForDiff(I), TBR -(I) ∪ use(I ) otherwise TBR -(I) = J predecessor of I TBR + (J)
-Adjoint-out analysis is specific to the adjoint mode of AD. It finds the variables whose value may be modified by the fragment of the adjoint code that consists of the sequel of the code to differentiate (i.e. the sequel of the forward sweep) followed by its adjoint (i.e. its backward sweep). In other words, seeing the code to differentiate P as the sequence P=U;I;D around instruction I, and recalling that its standard two-sweep Tapenade A:29

adjoint is:

P = - → U ; - → I ; - → D ; ← - D ; ← - I ; ← - U ,
then Adj-out + (I) contains the variables maybe modified by -→ D ; ← -D , and Adj-out -(I) those modified by

- → I ; - → D ; ← - D ; ← - I .
Adjoint-out analysis propagates the Adj-out sets backwards, with an initial value ∅ at the tail of P. A variable v enters the Adjout before instruction I if it is overwritten by the I that is contained in -→ I (when isLiveForDiff(I)).

← -I only overwrites derivative variables, and modifies no original variables. However, the possible store/restore mechanism controlled by the TBR analysis plays a role: if the value of a variable v is in TBR -(I), or will enter the TBR set because I uses it, then the Store-All strategy will make sure that whatever overwriting is done to v in -→ I ; -→ D it will be undone in ← -D ; ← -I . Therefore v can be removed from Adj-out -(I). Hence the data-flow equations:

Adj-out -(I ∞ ) = ∅ Adj-out -(I) =    (Adj-out + (I) ∪ out (I)) \ (TBR -(I) ∪ use(I )) if isLiveForDiff(I), Adj-out + (I) otherwise Adj-out + (I) = J successor of I Adj-out -(J)
We encounter a difficulty here: it turns out that the recomputation mechanism of section 4.9 uses the Adj-out sets, and modifies the TBR sets, thus causing a dependence cycle between the analyses. We will not detail the way we break this cycle, but we hope the present formalization can help us find a more accurate way to do so.

Inference Rules for Tangent differentiation

To specify actual differentiation of a procedure, we will focus on the case of wellstructured code, i.e. code where the control structure is apparent and embodied in the abstract syntax tree. This case is general enough, especially for industry codes that respect strict coding rules. For well-structured code, the differentiation algorithms can be specified in the formalism of Natural Semantics [START_REF] Kahn | Natural Semantics[END_REF]], a branch of Operational Semantics [START_REF] Plotkin | A structural approach to Operational Semantics[END_REF]] that uses inference rules to describe program transformation. The Natural Semantics rules that follow capture almost all the refinements presently implemented in Tapenade, except for a few optimizations that require a global vision of the code, discussed in section 5.5. In particular they describe exactly for both tangent and adjoint mode, how activity analysis benefits to actual differentiation and how the slicing resulting from diff-liveness analysis is done.

The algorithm described by these rules can be readily implemented in any programming language. Implementation is almost immediate in logical languages such as PROLOG. This yields a complete and faithful view of the actual implementation of tangent and adjoint AD in Tapenade, missing only a few global code optimization.

We will use a classical syntax for the inference rules: a possibly empty list of hypothesis predicates, on top of a fraction bar, on top of a conclusion predicate. The first hypothesis predicates can be boolean conditions that restrict applicability of the rule. The other hypotheses and the conclusion are rewrite predicates that turn a source tree pattern into a rewritten one. The source pattern is connected to the rewritten pattern by a → sign that can bear the name of the predicate. The predicate with no name can be read as "is differentiated as". Tree patterns may have tree variables at their leaves, displayed in upper-case ITALICS, that get instantiated by the corresponding sub-tree by matching, also called unification. Apart from tree variables, tree patterns are IL trees. We display them here with a concrete syntax, akin to a classical imperative language, and with a typewriter font. The inference rules we give are complete and mutually exclusive: one and only one rule is applicable to any given abstract syntax tree, so that no backtracking is needed during the tree transformation process.

The boolean conditions that appear in the rules only depend on the result of the static data-flow analysis. By default, they are checked on the source tree to be matched, and otherwise the boolean condition specifies as its argument the tree to be tested.

-isActiveExpr(EXPR) is true if EXPR is (or uses in a differentiable way) a reference to a variable that belongs to Active -. -isActiveCall(CALL) is true if one argument of this CALL is (or uses in a differentiable way) a reference to a variable that belongs to Active -or to Active + . -isLiveForDiff(I) is true if the instruction I produces a value that is diff-live, i.e. useful for the derivatives, i.e. out(I) ∩ Diff-live

+ (I) = ∅. -isRef(EXPR) is true if EXPR is a memory reference, e.g. x or A[i] but not x*2.0.
-isDiffProc(P) is true if isActiveCall is true for at least one call to P.

-isDiffFormalArg(ARG) is true if at least one call to the current procedure feeds into this formal argument an actual argument that belongs to Active -or to Active + . -isDiff(REF) is true if memory reference REF refers to a variable that is Active at some location in the program.

The rules also use a few primitives to create new symbol names, with the guarantee that the returned name is not already used in the original program.

-varName builds a name for the differentiated counterpart of a variable, generally appending "d" (tangent/dot) or "b" (adjoint/bar) to the original variable name. -procName builds a name for the differentiated counterpart of a procedure, generally appending " D" or " B" to the original procedure name. -newFreeVar builds a new name for a new variable of given type.

The differentiated program must contain declarations for these new symbols. The inference rules deal with that through the rewrite predicates (not detailled here):

DECLS → • DECLS and DECLS → DECLS
In tangent mode, a procedure is differentiated if at least one call to it is active. The envelope of all active calls determine which formal arguments are differentiated. This results in a single differentiated procedure. The corresponding inference rules are:

isDiffProc(P) P procName ------→ • P ARGS → • ARGS INSTRS → • INSTRS DECLS → • DECLS procedure P ( ARGS ) {DECLS; INSTRS} → procedure • P( • ARGS) { • DECLS; • INSTRS} isDiffFormalArg(ARG) ARG varName ------→ • ARG ARGS → • ARGS (ARG . ARGS) → (ARG , • ARG . • ARGS) !isDiffFormalArg(ARG) ARGS → • ARGS (ARG . ARGS) → (ARG . • ARGS)
() → () The rules for the sequence preserve the order of the original instructions:

Tapenade A:31 INSTR → • INSTR INSTRS → • INSTRS {INSTR ; INSTRS} → { • INSTR ; • INSTRS} {} → {}
Tangent differentiation is pushed deep down the control structures. There are special optimized cases when the control structure is not diff-live. Note that when the control structure is active, i.e. contains at least one active instruction, then it is labelled diff-live too. For instance for conditionals and loops:

isLiveForDiff() PART1 → • PART1 PART2 → • PART2 if (TEST) PART1 else PART2 → if (TEST) • PART1 else • PART2 !isLiveForDiff() if (TEST) PART1 else PART2 → {} isLiveForDiff() INSTRS → • INSTRS loop(CONTROL) INSTRS → loop(CONTROL) • INSTRS !isLiveForDiff() loop(CONTROL) INSTRS → {}
The typical atomic instruction is the assignment. Four cases must be considered, depending on the assignment being diff-live or not, and on the left-hand side reference having a derivative or not: REF:=EXPR → {} A procedure call is in general replaced by a call to the differentiated procedure, which contains both the original statements of the procedure and their differentiated counterparts. There are specialized rules when the particular call is active or not, or diff-live or not. Because of our generalization design choice, the activity of the actual argument may differ from that of the formal argument. This requires dummy differentiated arguments to be initialized and passed:

{} isRef() isDiff() EXPR ref --→ • EXPR EXPR actualArg ------→{}, • EXPR, {} !(isRef() isDiff()) typeOf(EXPR) newFreeVar --------→V V varName ------→ • V EXPR expr ---→ • EXPR EXPR actualArg ------→{ • V:= • EXPR}, • V, {}
Finally, tangent differentiation of a single expression applies the usual rules of calculus. For simplicity, only the rule for the product is presented here. Activity analysis allows us to define an improved rule that replaces the derivative of a non-active expression with zero. This opens the way for further simplifications of the differentiated expressions. We will not detail this (optional) simplification done by predicate simplify -----→. This completes the inference rules for the tangent mode of AD.

Inference Rules for Adjoint differentiation

The rules that we are going to list express all the adjoint differentiation model, except for a few global manipulations that require a global vision of the code (cf section 5.5). They express in particular -how TBR analysis is used to reduce the memory storage, -the exact implementation of user control on coarse-grain checkpointing, and of the other specific improvements described in section 4.10.

The rules can also capture a few other technical points that we decided to skip for clarity. These are the treatment of pointers and memory allocation, call-by-value arguments, and the recomputation of cheap intermediate variables described in section 4.9. Also, the general mechanism used by the inference rules on assignements to cope with possible aliasing is too systematic and is actually applied only when needed.

In addition to the boolean conditions and utility primitives already defined for the tangent mode, the inference rules for the adjoint mode use new conditions related to checkpointing:

-isJoint(CALL) is true if the choice for this CALL is joint differentiation i.e., to apply checkpointing to this CALL -isCalledJoint(P) is true if at least one call to this procedure P was selected for joint differentiation. -isCalledSplit(P) is true if at least one call to this procedure P was not selected for joint differentiation.

We also use two new primitives to create procedure names for the split mode:

-procNameFwd builds a name for the split forward sweep of the adjoint of a procedure, generally appending " FWD" to the original procedure name. -procNameBwd builds a name for the split backward sweep of the adjoint of a procedure, generally appending " BWD" to the original procedure name.

Checkpointing a piece of code INSTRS opens the way for simplification (cf section 4.9) because the forward sweep is immediately followed by the backward sweep. This requires that the data-flow analyses influenced by checkpointing (diff-liveness and TBR) be run again on the checkpointed program piece before it is differentiated. After differentiation of the checkpointed piece, the analyses results are reset to their previous values. This mechanism is expressed in the rules by "With CKPdataFlow(INSTRS) :"

In adjoint mode, a procedure can be differentiated in one of two ways (or both) depending on the checkpointing choice made by the end-user on each call. By default ("joint mode") checkpointing is applied to the call, which implies building one differentiated procedure which calls the forward sweep followed by the backward sweep of adjoint AD. Alternatively, the end-user can ("splitt mode") choose not to apply checkpointing on some calls, which implies building two differentiated procedures, one for the forward sweep and one for the backward sweep. It may happen that a procedure must be differentiated both in joint and split modes, thus creating three differentiated procedures. The inference rules are: The rules for an assignment distinguish the same cases as for the tangent mode. Actual computation of derivatives is placed into the backward sweep. Saving and restoring of the intermediate values is done here, through a push and pop of the variables RVARS overwritten by this assignment I that are needed in derivative expressions. Predicate restore ----→ computes RVARS = out(I) ∩ (TBR -(I) ∪ use(I )). Notice also the systematic introduction of the temporary variable V, which is necessary to cope with aliasing in the adjoint mode. Consider the assignment a(i) := 2*a(j) and assume there is no static way to decide whether i equals j or not. Introducing V is the only way to create correct adjoint instructions. However we can spare introduction of V when the aliasing question between left-and right-hand sides of the assignment is decidable statically. The inference rules here do not show this improvement.

!isLiveForDiff() isDiff(REF) REF ref --→ REF typeOf(REF) newFreeVar --------→V V varName ------→V EXPR, V expr ---→ INSTRS REF:=EXPR → {} V:=REF; REF:=0.0; INSTRS !isLiveForDiff() !isDiff(REF) REF:=EXPR → {} {}
The rules for differentiating a procedure call are different whether a joint or split differentiation is chosen. As joint mode actually means checkpointing the procedure call, we must compute the snapshots, i.e. the set SNP of variables that must be restored in order to repeat execution of the call, and possibly another set SBK of variables that must be restored after execution of the adjoint call P. There are many possible choices for SNP and SBK and no systematically better choice. The choice made by Tapenade and returned by predicate snapshot ------→ is for any checkpointed piece of code C:

SNP = out(C) ∩ TBR -(C) ∪ Diff-live -(C) SBK = ∅
The generated code uses procedure look, that restores values from the stack without actually popping them. In Joint mode, we can also make an extra simplification when the call is not diff-live: The next rule implements the possibility offered to the end-user to require application of checkpointing to any well-structured piece of source code: The rules for adjoint differentiation of the actual arguments of a procedure call are very similar to the tangent case. However, the initializations INITS before the call and the trailing POSTS instructions are different:

isActiveCall() isJoint() isLiveForDiff() call P(ARGS) snapshot ------→ SNP,SBK P procName ------→P ARGS actualArgs -------→INITS
INSTRS snapshot ------→ SNP,SBK With CKPdataFlow(INSTRS) : INSTRS → -----→ INSTRS ← ----- INSTRS checkpoint(INSTRS) → push (
() isDiff() EXPR ref --→EXPR EXPR actualArg ------→{}, EXPR, {} !(isRef() isDiff()) typeOf(EXPR) newFreeVar --------→V V varName ------→V EXPR, V expr ---→POSTS EXPR actualArg ------→{V:=0.0}, V, POSTS isRef
Finally, adjoint differentiation of a single expression catches the specificity of adjoint differentiation: an adjoint expression VAL is propagated through the expression EXPR down to its leaves (i.e. references to an active variable), updating the adjoint expression each time it goes down into an arithmetic operator (only the rule for the product is presented here). For each active leaf, a new assignment statement is produced that increments the adjoint variable of this leaf with the adjoint expression.

Implementation notes for analyses and transformations

Tapenade runs data-flow analysis to the code's internal representation, then applies differentiation to it. This results into the internal representation of a new code.

As each data-flow analysis runs on a Call Graph whose nodes are actually Flow Graphs, the data-flow equations that we provided are implemented with a number of sweeps on these graphs. Depending on the equations, these sweeps can be top-down or bottom-up (on the Call Graph), and forward or backward (on the Flow Graphs). For instance activity analysis requires a top-down sweep on the Call Graph, which performs for each Flow Graph a forward sweep for the Varied information and a backward sweep for the Useful information. Because of cycles, these sweeps must iterate until reaching a fixed point. For efficiency, these iterations visit the graph nodes in a well-chosen order (e.g. depth first spanning tree for a forward sweep) and maintain a wait list of nodes that require a new iteration. Implementation-wise, this very uniform structure of data-flow analyses allows each of them to inherit from a parent Java class that defines these sweeps.

Differentiation itself is implemented at the upper level (Call Graph) as a quite simple graph transformation: activity analysis has labelled the procedures that must be differentiated, and differentiation amounts just to creating the new, differentiated, procedure(s). When the application language supports a notion of modules, then differentiated modules are created accordingly to contain the differentiated procedures and declarations, together with the original ones to cope with visibility restrictions. There are other minor subtleties, irrelevant here.

At the intermediate level, each procedure is indeed stored as a graph (Flow Graph) and differentiation is therefore a graph transformation. In particular, this enables Tapenade to equally differentiate structured and unstructured code. This implementation is thus slightly different from its inference rules specification of sections 5.3 and 5.4, but for well-structured code the result of the transformation is identical to the code built with the inference rules. The implemented algorithm essentially consists of two phases:

(1) Build new, differentiated Basic Block(s) for each Basic Block of the original flow graph. At this deepest level of our program representation, which uses abstract syntax trees, the inference rules apply directly. (2) Connect the new Basic Blocks with Control Flow arrows according to the Control Flow arrows of the original flow graph. This is easy in tangent mode and in the forward sweep of the adjoint mode. In the backward sweep, the basic idea is to remember the control at each place where the flow merges in the forward sweep, and to use this memory to take the corresponding direction backwards in the backward sweep. This graph approach yields a better dead code elimination at the junction of forward and backward sweeps than can be expressed by the inference rules.

What Tapenade does that the inference rules do not express is the following non-local manipulations on the differentiated code:

-the reordering and the common subexpression elimination that take place in each differentiated expression. -the additional slicing and simplifications that occur in adjoint mode at the junction of forward and backward sweeps. -the implicit-zero mechanism that spares initialization of some derivatives, or at least postpones it till it is absolutely necessary. -in the adjoint mode, the reordering and fusion of increments to the same differentiated variable, sometimes known as incremental vs non-incremental adjoint mode.

Tapenade A:39

-in the backward sweep of the adjoint mode, the choice to recompute an intermediate variable instead of restoring it from the stack, as shown on the rightmost column of figure 12.

DEVELOPMENT STATUS AND PERFORMANCES

Development of Tapenade started in 1999 as a major redesign of the AD tool Odyssee [START_REF] Rostaing | Automatic Differentiation in Odyssée[END_REF]. In its present distributed version (v3.6, revision 4344), Tapenade source is mainly 138,000 lines of Java. The development environment also features 1,120 automated non-regression tests and 15 complete validation codes.

From the end-user point of view, Tapenade can be used as a web server at url http://www-tapenade.inria.fr:8080. Tapenade can also be installed locally on most architectures, which gives access to the full set of command-line options. All Tapenade documentation is available online with a user's guide, a tutorial, a FAQ, and downloading instructions. Users may subscribe to a dedicated mailing list. Downloading Tapenade is free for nonprofit academic research or education. Although still incomplete, the documentation describes the parameters of a differentiation request, the command-line options, and the AD directives that can be inserted in the source program. In addition to the differentiated code, Tapenade may return error and warning messages, discussed in the documentation.

Tapenade provides a special mode to validate its own output. This validation mode also provides support to locate possible bugs in the differentiated code. Validation of the tangent code is done by comparing the tangent derivatives with an approximation computed by divided differences. Validation of the adjoint code is done by computing a signature of the function's Jacobian matrix, using the tangent mode on one hand, and the adjoint mode on the other hand. Recall the notations of a program P that computes a function F : X ∈ IR n → Y ∈ IR m , call F (X) its Jacobian matrix at point X, Ẋ a column vector of IR n and Y a row vector of IR m . Examples of possible signatures for comparison are:

-The "checksum", that computes Y × F (X) × Ẋ for any arbitrary Ẋ and Y . By associativity, this can be computed as Y × (F (X) × Ẋ) using the tangent code, or as (Y × F (X)) × Ẋ using the adjoint code. -The "dot-product", that computes Ẋ * × F (X) * × F (X) × Ẋ, where the * denotes transposition, for an arbitrary Ẋ. This can be computed as the square norm of Ẏ = F (X) × Ẋ, which is obtained using the tangent code, or by setting Y = Ẏ * , then computing X = Y × F (X) using the adjoint mode, and finally computing the dot product Ẋ.X.

Derivatives are useful in most fields of Scientific Computing, and AD tools only provide a way to obtain them efficiently. In particular, gradients computed through the adjoint mode are one proeminent achievement of AD. For some classes of applications, there have been fruitful joint research between the Scientific Computing and AD specialists, to devise better computing algorithms that rely on specially tailored kinds of derivatives [START_REF] Dervieux | Optimization loops for shape and error control[END_REF]]. This meant joint progress for both Scientific Computing and AD. Still, derivatives are standalone mathematical objects, and the way they are used by scientific applications is disconnected in principle from the way an AD tool provides them. Therefore, from the AD tool point of view, it is not relevant here to go into the detail of particular Scientific Computing applications. To illustrate the performances of Tapenade, we have selected some of its applications that we are aware of. For each, after a brief description of the application, we will concentrate on the strictly AD-related performance figures i.e., the accuracy of the derivatives and the time and memory costs of the differentiated program.

-uns2D is a 2D steady-state Navier-Stokes CFD simulation, on an unstructured mesh with 3000 triangles, and which on this application converges after 566 iterations then outputs a handful of aerodynamic characteristics (lift, drag. . . ). -nsc2ke is a small unsteady Euler CFD simulation, adapted for the needs of the car industry. On this particular simulation, the size of the unstructured mesh is 1516 triangles and the number of time steps is fixed to 80. -lidar simulates the propagation of light across the atmosphere. As an inverse problem, it is used to estimate parameters of the atmosphere from measurements of the light received by a lidar. -nemo is a large "configuration" of the OPA [START_REF] Madec | OPA8.1 ocean general circulation model reference manual[END_REF]] ocean circulation model, here solving an inverse problem around the Antarctic. Mesh elements are 2 o wide, and this simulation runs for 400 time steps, i.e. about one month. -gyre is a reduced configuration of the OPA model on a rectangular basin meshed with 21824 cells. The test computes the derivatives of a short simulation of 10 time steps. -winnie is a simplified 2D version of the grisli model, that simulates the movement of polar ice caps. This particular unsteady simulation runs on a mesh with 250 elements, for 1000 time steps. -stics [START_REF] Lauvernet | Using Automatic Differentiation to study the sensitivity of a crop model[END_REF] simulates the growth of a crop during one year (350 time steps) depending on climate and agricultural practices. The goal of this particular study was to evaluate the sensitivity of the biomass produced on the agricultural parameters. -smac-sail [START_REF] Lauvernet | Multitemporal-patch ensemble inversion of coupled surface-atmosphere radiative transfer models for land surface characterization[END_REF]] are two coupled atmosphere-surface radiative transfer models used here to estimate biophysical variables, from canopy reflectance data observed from a satellite. When used to minimize a discrepancy function, the code can be seen as having a scalar output. For other uses, it can have many outputs, like 7801 here. -traces simulates the transport of pollutants through layers of soil and rock, carried by water. It was used here in an inverse problem fashion to estimate the permeability parameters of rock. Unfortunately the quantity of available measurements was very small in this application, so only 8 parameters could be estimated. For each application, table 18 gives the application language, the size in lines of the source part that was submitted for differentiation, and its average run-time T . It gives also the dimensions n and m of the independent input space and dependent output space and the differentiation time T AD to build the adjoint. The time to build the tangent code is a little shorter. For the tangent mode, the table gives the accuracy A t of tangent derivatives tt with respect to divided differences dd , which basically measures the inaccuracy of divided differences. Actually A t = -log 10 ( tt-dd dd ) so that better accuracy means higher A t . The table then shows the ratio R t of the time to obtain one tangent derivative Ẏ divided by the time T of the original code. For the adjoint mode, the table gives (the opposite log 10 of) the accuracy A a with respect to the tangent derivatives. This accuracy A a should only reflect the different computation order of the derivatives and should therefore be high. On nsc2ke and nemo, A a is not as good as expected (and so is A t ). We found this comes from the iterative linear solvers inside. Forcing extra iterations of the adjoint solvers fixes the problem, but further studies are needed before we automate this fix. The table then shows the ratio R a of the time to obtain one gradient X divided by T . These time ratios along with n and m allow one to check the superiority of the adjoint mode for computing gradients. Specifically for the adjoint mode, we finally give the peak size of the trajectory stack, and the total traffic of data pushed to and popped from this stack. Most of these applications make use of Tapenade directives on checkpointing to reduce these memory costs.

OUTLOOK

We presented the specification of the Automatic Differentiation tool Tapenade. This specification is naturally split in two parts, one about the necessary static analysis of the provided source program, and one about its actual transformation into the derivative source program. Tapenade follows a model of AD that is only one in a panorama of possible models. We motivated why we chose this model, and illustrate it on examples.

Tapenade has been applied with success to several large codes, both academic and industrial, and its adjoint mode has often been considered valuable by end users, especially by those having written an adjoint code by hand in the past. Much work remains to make it even more widely applicable. Our long list of future developments contains specific differentiation of well-known patterns such as solution of linear systems or fixed-point iterations: the standard AD strategy applied by Tapenade works, but research has come up with adapted strategies that are more efficient and accurate. Also, we are considering a tool interface that could give freedom to differentiate an intermediate variable wrt another. We also want to give the option of association by address instead of by name, and evaluate the performance benefit brought by a better memory locality. Support for repeated differentiation must be improved, in order to compute second derivatives by Tangent-over-Adjoint differentiation. This last point is gaining importance with the growing concern for Uncertainty Quantification.

We identify three more ambitious developments, each of them meaning actual research on some aspects of the model that are still not fixed. Most of the difficulty lies with the adjoint mode. These are (1) Message-passing communication, (2) Dynamic memory management, and (3) Object-Oriented languages.

Message-passing communication can be theoretically transposed in the backward sweep of the adjoint. This exhibits very elegant duality between "sends" and "receives", which are adjoint of one another. Non-blocking (or asynchronous) communication such as an mpi isend/mpi wait pair, causes the backward sweep to contain an mpi irecv/mpi wait pair, in reverse order. There are a number of open issues on the best way to detect matching "sends" and "receives", as well as matching asynchronous communication pairs. These questions are not safely implemented in Tapenade yet.

Dynamic memory management causes the adjoint backward sweep to use memory that was deallocated by the forward sweep. One answer is to allocate some memory again, but this can return a different chunk of memory. Therefore, saved addresses are inconsistent with the new chunk of memory. We investigate two complementary answers. First, saved addresses can also save their offset with respect to the allocated base address, and one can recompute the new addresses from the saved offsets. This amounts to dynamically monitoring allocation and deallocation, which is quite general but costly. In favorable cases, we plan to mitigate with a special efficient strategy that recomputes memory addresses on the fly when possible.

Object-Oriented languages imply a deep change in the internal representation of programs. Because of the Object point of view, association by address seems more adapted for these languages. It also seems reasonable to first address Java rather than C++, although implicit deallocation by Garbage Collection makes it harder to locate the adjoint allocation during the backward sweep.

Fig. 3 .

 3 Fig.3. Checkpointing with the Recompute-All (left) and Store-All (right) approaches. The shaded areas reproduce the basic pattern of the chosen approach. top: single checkpointing, bottom nested checkpointing.

Fig. 4 .

 4 Fig. 4. Names of differentiated symbols

  TYPE (VECTOR) u(10),v,w TYPE (VECTOR_D) ud(10) ... FUNCTION TEST_D(a,ad,b,bd,test) REAL TEST, TEST_D TYPE (VECTOR) a,b TYPE (VECTOR_D) ad,bd TEST_D = a%x*bd%x + ad%x*b%x + ud(3)%z TEST = a%x*b%x + u(3)%z ...

Fig. 5 .

 5 Fig. 5. Differentiated Data Types

  Fig. 7. Tapenade tangent and adjoint on some code with control original code tangent adjoint (no checkpointing) y = FOO(a,x,2.5) b = FOO(c,z,t) yd = FOO_D(a,x,xd,2.5,0.0,y) dumm = FOO_D(c,z,zd,t,td,b) y = FOO_FWD(a,x,2.5) b = FOO_FWD(c,z,t)... bb = 0.0 CALL FOO_BWD(c,z,zb,t,tb,bb) CALL FOO_BWD(a,x,xb,2.5,dummb,yb)

Fig. 9 .

 9 Fig. 9. Effect of activity analysis on Tangent AD model. Independents are a, b, c. Dependents are d, e. Partial evaluation follows from a becoming not varied after receiving constant 5.0. Slicing follows from e being not useful before being truncated by the call to floor.

  Fig.10. Tangent G D and vector tangent G DV derivatives of a procedure G. nbdirsmax is a static constant to be defined in module DIFFSIZES, that defines the size of the extra dimension. Argument nbdirs (≤nbdirsmax) is the actual number of differentiation directions for one particular run.

Fig. 13 .

 13 Fig. 13. The effect of checkpointing on the adjoint of a procedure call call to FOO is checkpointed, the forward sweep stores a snapshot which is specific for each call site, then calls plain FOO. The backward sweep restores the snapshot, then calls FOO B which contains both forward and backward sweeps for FOO.When checkpointing is applied to C, the adjoint code contains the sequence-→ C ; ← -C , in addition to the sequence -→ D ;← -D already present before checkpointing is applied. Recalling that the code improvements such as diff-liveness (cf section 4.9) bear more fruit when the forward sweep is immediately followed by the backward sweep, we see that checkpointing opens more opportunity to simplify the adjoint code. In any case, when checkpointing is applied, an AD tool should restart the static data-flow analyses for the checkpointed part, in order to discover these additional simplifications.The default checkpointing strategy of Tapenade can be amended in several ways. The directive $AD NOCHECKPOINT attached to a procedure call deactivates checkpointing on this particular call. To deactivate checkpointing on all calls, attach the directive to the procedure header itself. Conversely, the pair of directives $AD CHECKPOINT-START and $AD CHECKPOINT-END define an arbitrary portion of code (well-structured) on which checkpointing must be applied.Concerning loops, directive $AD BINOMIAL-CKP attached to the loop header tells Tapenade to apply binomial checkpointing as defined in[START_REF] Griewank | Achieving logarithmic growth of temporal and spatial complexity in reverse Automatic Differentiation[END_REF]], which is the optimial strategy to apply to iterative loops such as time-stepping loops. Finally, directive $AD II-LOOP attached to the loop header tells Tapenade that the loop has Independent Iterations. Checkpointing can be therefore applied to each iteration of the loop, thus reducing the peak memory used for trajectory storage[Hascoët et al. 2002].

Fig. 14 .

 14 Fig. 14. Adjoints of a few example Fortran array operations

Fig. 16 .

 16 Fig. 16. Internal representation of a program inside Tapenade

Fig. 17

 17 Fig. 17. The chain of data-flow analyses implemented in Tapenade

.

  Fig. 17. The chain of data-flow analyses implemented in Tapenade

Varied

  + (I) = Varied -(I) ⊗ Diff-dep(I) Useful -(I) = Diff-dep(I) ⊗ Useful + (I) in which the diff-dependency relation Diff-dep(I) between variables is defined by (v1.v2) ∈ Diff-dep(I) iff the variable v2 after I depends in a differentiable way on Tapenade A:27 the variable v1 before I. Composition ⊗ is defined naturally as: v2 ∈ S ⊗ Diff-dep(I) ⇐⇒ ∃v1 ∈ S|(v1.v2) ∈ Diff-dep(I)and likewise in the reverse direction Diff-dep(I) ⊗ S. When I is a simple statement such as an assignment, Diff-dep(I) can be computed on the fly: the left-hand-side variable depends differentiably on all right-hand-side variables except those that occur in a non-differentiable way such as indices or other address computation. All other variables are unchanged and therefore depend on themself only. When the lefthand-side is only a part of its variable v, e.g. an array element, we take the other elements of v into account, conservatively adding (v.v) into Diff-dep(I). When on the other hand I is a call to a procedure P , then we rely on the preliminary diffdependency analysis to provide us with Diff-dep(P ). The data-flow equations are closed with the equations at control-flow branches:Varied -(I) = J predecessor of I

(

  v1.v3) ∈ Diff-dep(A) ⊗ Diff-dep(B) ⇐⇒ ∃v2|(v1.v2) ∈ Diff-dep(A)&(v2.v3) ∈ Diff-dep(B) then we can compute Diff-dep(P ) by a forward sweep on P , exploiting again the natural transitivity of dependency. Calling I 0 the entry instruction of P , I ∞ its exit instruction, and Id the identity dependence relation, the data-flow equations are: Diff-dep + (I 0 ) = Id Diff-dep + (I) = Diff-dep -(I) ⊗ Diff-dep(I) Diff-dep -(I) = J predecessor of I Diff-dep + (J) and the result Diff-dep(P ) is found in Diff-dep -(I ∞ ).

  ----→{}, (),

ACM

  Transactions on Mathematical Software, Vol. V, No. N, Article A, Publication date: 0.

  ; VF:=F; VS:=S; loop(I=VF,T,VS) -----→ INSTRS; push(I-VS,VF,-VS) pop(VF,VT,VS); loop(I=VF,VT,VS)

  This completes the inference rules for the adjoint mode of AD.

  Table sees the symbols that are defined in its parent Symbol Table. Each Basic Block has a link to the deepest Symbol Table it can see. Searching for the definition of a symbol from a Basic Block starts from this Symbol Table, then if not found continues to the parent Symbol Table possibly until the "root" Symbol Table.Most symbols in the Symbol Table are in fact associated with a kind (variable, procedure. . . ) a type, and miscellaneous information. When the symbol is a procedure name,

	Unit: P7

root Symbol Table: M6: module M3: module Symbol Table of M3: P5: procedure Symbol Table of P5: A: variable, real(:) b: variable, real i: variable, integer Symbol Table of M3 (copy): P5: procedure Symbol Table of M6: P7: procedure v: variable, real Symbol Table of P7: i: variable, integer Local Symbol Table: tmp: variable, real

  Performances of Tapenade derivative code on representative applications

			original code	tangent		adjoint
		lines	T	n → m	TAD At	Rt	Aa	Ra peak traffic
			(s)		(s)			(Mb)	(Mb)
	uns2d (F77)	2000 1.3 14000 → 3	6 3.4 2.4 15.1 5.9	241	1243
	nsc2ke (F77)	3500 0.4 1602 → 5607	11 1.9 2.4	4.5 16.2	168	2806
	lidar (F90)	330 3.3	37 → 37	2 6.7 1.1 14.4 2.0	11	11
	nemo (F90)	55000 208 9100 → 1	95 3.0 2.0	8.1 6.5 1591 85203
	gyre (F90)	21000	30 21824 → 1	26 4.5 1.9 13.3 7.9	481 48602
	winnie (F90)	3700 0.6	3 → 1	8 1.4 1.7 13.7 5.9	421	614
	stics (F77)	17000 0.2	739 → 1467 206 8.6 2.4 15.3 3.9	155	186
	smac-sail (F77)	1700 2.2 1321 → 7801	6 5.9 1.0 10.5 3.1	2	21
	traces (F90)	19800 13.2	8 → 1	53 4.0 1.3 12.9 3.8	159	4390
	Fig. 18.

ACM Transactions on Mathematical Software, Vol. V, No. N, Article A, Publication date: 0.

Fortran 95 standard uses the term derived type. This terminology leads to confusion with "derivative", of course frequent in AD. We'll rather speak of a structured type ACM Transactions on Mathematical Software, Vol. V, No. N, Article A, Publication date: 0.

A:34 L. Hasco ët and V. Pascual isDiffProc(P) isCalledSplit(P) P procNameFwd ---------→ -→

One specificity of the adjoint AD model is the two-sweeps structure, forward then backward, of the adjoint code. The rules for the sequence and for the control structures reflect that by producing the pair of the forward and backward codes: