### Running Example

The Farfallen release is available here.

This is an elaboration of the running example from
Radul, A., Pearlmutter, B., and Siskind, J.M., `AD in Fortran, Part 1: Design,' submitted to 6th International Conference on Automatic Differentiation, Advances in Automatic Differentiation, Lecture Notes in Computational Science and Engineering, Springer:Berlin, 2012.
which finds an equilibrium of a two-player game with continuous scalar strategies.

This page shows in detail how the Farfel program from Listing 1 in the paper transforms to Fortran77, using Tapenade 3.6 and adifor to perform the AD transformations.

The intermediate result of each stage of the transformation is given by a link to a complete Farfel program. The final result is a runnable Fortran77 program.

1. We start with the code from Listing 1

Initial Farfel program: html text

2. First, we extract nested subprogram definitions to the top level, starting with `fprime` from `argmax`
1. We transform
```      function argmax(f, x0, n)
function fprime(x)
fprime = deriv1(f, x)
end
argmax = root(fprime, x0, n)
end
```
into
```      function argmax_fprime(x)
argmax_fprime = deriv1(f, x)
end

function argmax(f, x0, n)
argmax = root(argmax_fprime, x0, n)
end
```
by extracting the nested subprogram `fprime` to the top level and renaming it argmax_fprime.

Intermediate state: html text

2. Then, since argmax_fprime has a free reference to f, we add it as a closure argument:
```      function argmax_fprime(x, f)
argmax_fprime = deriv1(f, x)
end
```
Intermediate state: html text

3. Then, since argmax_fprime is passed as an external to root, we add the closure argument at the call site:
```      argmax = root(argmax_fprime, f, x0, n)
```
Intermediate state: html text

4. Then we specialize root to accept the needed closure argument:
```      function root_1(f, f1, x0, n)
x = x0
do 1669 i = 1, n
call deriv2(f, f1, x, y, yprime)
1669 x = x - y / yprime
root_1 = x
end
```
and adjust the call site to use the correct version of `root`:
```      argmax = root_1(argmax_fprime, f, x0, n)
```
Intermediate state: html text

5. Then we specialize deriv2 for the same reason as `root`:
```      subroutine deriv2_1(f, f1, x, y, yprime)
external f
y = f(x, f1)
end adf(yprime = tangent(y))
end
```
and adjust the call site accordingly:
```      call deriv2_1(f, f1, x, y, yprime)
```
Valid Farfel program: html text

3. Then we constant propagate argmax_fprime from the first argument to the call to root_1 all the way to call site in deriv2_1, eliminating the arguments and parameters along the way as well as the external declaration:
```      function root_1(f1, x0, n)
x = x0
do 1669 i = 1, n
call deriv2_1(f1, x, y, yprime)
1669 x = x - y / yprime
root_1 = x
end

subroutine deriv2_1(f1, x, y, yprime)
y = argmax_fprime(x, f1)
end adf(yprime = tangent(y))
end

argmax = root_1(f, x0, n)
```
Valid Farfel program: html text

4. Now, we apply this same process to `f` from eqlbrm.
1. The subroutine
```      subroutine eqlbrm(biga, bigb, astar, bstar, n)
external biga, bigb
function f(astar)
function g(a)
function h(b)
h = bigb(astar, b)
end
bstar = argmax(h, bstar, n)
g = biga(a, bstar)
end
f = argmax(g, astar, n) - astar
end
astar = root(f, astar, n)
end
```
becomes
```      function eqlbrm_f(astar)
external biga, bigb
function g(a)
function h(b)
h = bigb(astar, b)
end
bstar = argmax(h, bstar, n)
g = biga(a, bstar)
end
eqlbrm_f = argmax(g, astar, n) - astar
end

subroutine eqlbrm(biga, bigb, astar, bstar, n)
astar = root(eqlbrm_f, astar, n)
end
```
by extracting the nested subprogram f to the top level and renaming it to `eqlbrm_f`.

Intermediate state: html text

2. Then, since f has free references to biga, bigb, bstar, and n, we add them as closure arguments:
```      function eqlbrm_f(astar, biga, bigb, bstar, n)
function g(a)
function h(b)
h = bigb(astar, b)
end
bstar = argmax(h, bstar, n)
g = biga(a, bstar)
end
eqlbrm_f = argmax(g, astar, n) - astar
end
```
Intermediate state: html text

3. Then, since eqlbrm_f is passed as an external to root, we add the closure arguments at the call site:
```      astar = root(eqlbrm_f, biga, bigb, bstar, n, astar, n)
```
Intermediate state: html text

4. Then we specialize root again to accept closure arguments for its parameter`f`:
```      function root_2(f, biga, bigb, bstar, x0, n)
x = x0
do 1669 i = 1, n
call deriv2(f, x, y, yprime)
1669 x = x - y / yprime
root_2 = x
end
```
and adjust the call site to call the specialization:
```      astar = root_2(f, biga, bigb, bstar, astar, n)
```
Note that we eliminate the extra copy of n in both the subprogram definition and the call site.

Intermediate state: html text

5. Then we specialize deriv2 again:
```      subroutine deriv2_2(f, biga, bigb, bstar, n, x, y, yprime)
external f
y = f(x, biga, bigb, bstar, n)
end adf(yprime = tangent(y))
end
```
and redirect its call site:
```      call deriv2_2(f, biga, bigb, bstar, n, x, y, yprime)
```
Valid Farfel program: html text

5. Then we constant propagate eqlbrm_f from the first argument to the call to root_2 all the way to call site in deriv2_2, eliminating the arguments and parameters along the way as well as the external declaration:
```      astar = root_2(biga, bigb, bstar, astar, n)

function root_2(biga, bigb, bstar, x0, n)
x = x0
do 1669 i = 1, n
call deriv2_2(biga, bigb, bstar, n, x, y, yprime)
1669 x = x - y / yprime
root_2 = x
end

subroutine deriv2_2(biga, bigb, bstar, n, x, y, yprime)
y = eqlbrm_f(x, biga, bigb, bstar, n)
end adf(yprime = tangent(y))
end
```
Valid Farfel program: html text

6. Now, we apply this same process to `g` from eqlbrm_f:
1. The function
```      function eqlbrm_f(astar, biga, bigb, bstar, n)
external biga, bigb
function g(a)
function h(b)
h = bigb(astar, b)
end
bstar = argmax(h, bstar, n)
g = biga(a, bstar)
end
eqlbrm_f = argmax(g, astar, n) - astar
end
```
becomes
```      function eqlbrm_f_g(a)
external biga, bigb
function h(b)
h = bigb(astar, b)
end
bstar = argmax(h, bstar, n)
eqlbrm_f_g = biga(a, bstar)
end

function eqlbrm_f(astar, biga, bigb, bstar, n)
eqlbrm_f = argmax(eqlbrm_f_g, astar, n) - astar
end
```
by extracting the nested subprogram g to the top level.

Intermediate state: html text

2. Then, since g has free references to astar, biga, bigb, bstar, and n, we add them as closure arguments:
```      function eqlbrm_f_g(a, astar, biga, bigb, bstar, n)
external biga, bigb
function h(b)
h = bigb(astar, b)
end
bstar = argmax(h, bstar, n)
eqlbrm_f_g = biga(a, bstar)
end
```
Intermediate state: html text

3. Then, since eqlbrm_f_g is passed as an external to argmax, we add the closure arguments at the call site:
```      eqlbrm_f = argmax(eqlbrm_f_g, astar, biga, bigb, bstar, n, astar, n) - astar
```
Intermediate state: html text

4. Then we specialize argmax to accept the closure arguments for `eqlbrm_f_g`:
```      function argmax_1(f, astar, biga, bigb, bstar, x0, n)
argmax_1 = root_1(f, astar, biga, bigb, bstar, n, x0, n)
end
```
and change the call site to call the specialization:
```      eqlbrm_f = argmax_1(eqlbrm_f_g, astar, biga, bigb, bstar, astar, n) - astar
```
Note that we eliminate the extra copy of n in both the subprogram definition and the call site.

Intermediate state: html text

5. Then we specialize root_1 to accept those same closure arguments:
```      function root_1_1(f1, astar, biga, bigb, bstar, x0, n)
x = x0
do 1669 i = 1, n
call deriv2_1(f1, astar, biga, bigb, bstar, n, x, y, yprime)
1669 x = x - y / yprime
root_1_1 = x
end
```
and adjust its call site:
```      argmax_1 = root_1_1(f, astar, biga, bigb, bstar, x0, n)
```
Note that we eliminate the extra copy of n in both the subprogram definition and the call site.

Intermediate state: html text

6. Then we specialize deriv2_1 to accept the needed closure arguments:
```      subroutine deriv2_1_1(f1, astar, biga, bigb, bstar, n, x, y, yprime)
y = argmax_fprime(x, f1, astar, biga, bigb, bstar, n)
end adf(yprime = tangent(y))
end
```
and change the call site to call the specialization:
```      call deriv2_1_1(f1, astar, biga, bigb, bstar, n, x, y, yprime)
```
Intermediate state: html text

7. Then we specialize argmax_fprime for these closure arguments:
```      function argmax_fprime_1(x, f, astar, biga, bigb, bstar, n)
argmax_fprime_1 = deriv1(f, astar, biga, bigb, bstar, n, x)
end
```
and retarget its call site:
```      y = argmax_fprime_1(x, f1, astar, biga, bigb, bstar, n)
```
Intermediate state: html text

8. Then we specialize deriv1 for these closure arguments:
```      function deriv1_1(f, astar, biga, bigb, bstar, n, x)
external f
y = f(x, astar, biga, bigb, bstar, n)
end adf(deriv1_1 = tangent(y))
end
```
and retarget its call site:
```      argmax_fprime_1 = deriv1_1(f, astar, biga, bigb, bstar, n, x)
```
Valid Farfel program: html text

7. Then we constant propagate eqlbrm_f_g from the first argument to the call to argmax_1 all the way to call site in deriv1_1, eliminating the arguments and parameters along the way as well as the external declaration:
```      eqlbrm_f = argmax_1(astar, biga, bigb, bstar, astar, n) - astar

function root_1_1(astar, biga, bigb, bstar, x0, n)
x = x0
do 1669 i = 1, n
call deriv2_1_1(astar, biga, bigb, bstar, n, x, y, yprime)
1669 x = x - y / yprime
root_1_1 = x
end

subroutine deriv2_1_1(astar, biga, bigb, bstar, n, x, y, yprime)
y = argmax_fprime_1(x, astar, biga, bigb, bstar, n)
end adf(yprime = tangent(y))
end

function argmax_fprime_1(x, astar, biga, bigb, bstar, n)
argmax_fprime_1 = deriv1_1(astar, biga, bigb, bstar, n, x)
end

function argmax_1(astar, biga, bigb, bstar, x0, n)
argmax_1 = root_1_1(astar, biga, bigb, bstar, x0, n)
end

function deriv1_1(astar, biga, bigb, bstar, n, x)
y = eqlbrm_f_g(x, astar, biga, bigb, bstar, n)
end adf(deriv1_1 = tangent(y))
end
```
Valid Farfel program: html text

8. Now, we apply this same process to `h` in eqlbrm_f_g:
1. The function
```      function eqlbrm_f_g(a, astar, biga, bigb, bstar, n)
external biga, bigb
function h(b)
h = bigb(astar, b)
end
bstar = argmax(h, bstar, n)
eqlbrm_f_g = biga(a, bstar)
end
```
becomes
```      function eqlbrm_f_g_h(b)
external bigb
eqlbrm_f_g_h = bigb(astar, b)
end

function eqlbrm_f_g(a, astar, biga, bigb, bstar, n)
external biga
bstar = argmax(eqlbrm_f_g_h, bstar, n)
eqlbrm_f_g = biga(a, bstar)
end
```
by extracting the nested subprogram h to the top level.

Intermediate state: html text

2. Then, since h has free references to astar and bigb, we add them as closure arguments:
```      function eqlbrm_f_g_h(b, astar, bigb)
external bigb
eqlbrm_f_g_h = bigb(astar, b)
end
```
Intermediate state: html text

3. Then, since eqlbrm_f_g_h is passed as an external to argmax, we add the closure arguments at the call site:
```      bstar = argmax(eqlbrm_f_g_h, astar, bigb, bstar, n)
```
Intermediate state: html text

4. Then we specialize argmax again, to accept closure arguments for `eqlbrm_f_g_h`:
```      function argmax_2(f, astar, bigb, x0, n)
argmax_2 = root_1(f, astar, bigb, x0, n)
end
```
and change its call site to refer to the specialization:
```      bstar = argmax_2(eqlbrm_f_g_h, astar, bigb, bstar, n)
```
Intermediate state: html text

5. Then we specialize root_1 again, for those same closure arguments:
```      function root_1_2(f1, astar, bigb, x0, n)
x = x0
do 1669 i = 1, n
call deriv2_1(f1, astar, bigb, x, y, yprime)
1669 x = x - y / yprime
root_1_2 = x
end
```
and change its call site:
```      argmax_2 = root_1_2(f, astar, bigb, x0, n)
```
Intermediate state: html text

6. Then we specialize deriv2_1 to accept the self-same closure arguments:
```      subroutine deriv2_1_2(f1, astar, bigb, x, y, yprime)
y = argmax_fprime(x, f1, astar, bigb)
end adf(yprime = tangent(y))
end
```
and adjust its call site:
```      call deriv2_1_2(f1, astar, bigb, x, y, yprime)
```
Intermediate state: html text

7. Then we specialize argmax_fprime again, for these closure arguments:
```      function argmax_fprime_2(x, f, astar, bigb)
argmax_fprime_2 = deriv1(f, astar, bigb, x)
end
```
and adjust its call site:
```      y = argmax_fprime_2(x, f1, astar, bigb)
```
Intermediate state: html text

8. Then we specialize deriv1 again, for these closure arguments:
```      function deriv1_2(f, astar, bigb, x)
external f
y = f(x, astar, bigb)
end adf(deriv1_2 = tangent(y))
end
```
and adjust its call site:
```      argmax_fprime_2 = deriv1_2(f, astar, bigb, x)
```
Valid Farfel program: html text

9. Then we constant propagate eqlbrm_f_g_h from the first argument to the call to argmax_2 all the way to the call site in deriv1_2, eliminating the arguments and parameters along the way as well as the external declaration:
```      bstar = argmax_2(astar, bigb, bstar, n)

function root_1_2(astar, bigb, x0, n)
x = x0
do 1669 i = 1, n
call deriv2_1_2(astar, bigb, x, y, yprime)
1669 x = x - y / yprime
root_1_2 = x
end

subroutine deriv2_1_2(astar, bigb, x, y, yprime)
y = argmax_fprime_2(x, astar, bigb)
end adf(yprime = tangent(y))
end

function argmax_fprime_2(x, astar, bigb)
argmax_fprime_2 = deriv1_2(astar, bigb, x)
end

function deriv1_2(astar, bigb, x)
y = eqlbrm_f_g_h(x, astar, bigb)
end adf(deriv1_2 = tangent(y))
end

function argmax_2(astar, bigb, x0, n)
argmax_2 = root_1_2(astar, bigb, x0, n)
end
```
Valid Farfel program: html text

10. Now we notice that the unspecialized root, deriv2, argmax, root_1, deriv2_1, argmax_fprime, and deriv1 are unused and eliminate them.

Valid Farfel program: html text

Now we are done closure-converting nested subprogram definitions.
11. The next step is to canonicalize `adf` blocks to be single subroutine calls. We transform
```      function deriv1_1(astar, biga, bigb, bstar, n, x)
y = eqlbrm_f_g(x, astar, biga, bigb, bstar, n)
end adf(deriv1_1 = tangent(y))
end
```
into
```      function deriv1_1(astar, biga, bigb, bstar, n, x)
y = eqlbrm_f_g(x, astar, biga, bigb, bstar, n)
end
end adf(deriv1_1 = tangent(y))
end
```
and then closure-convert (by re-applying the above process) to get
```      subroutine deriv1_1_adf1(x, astar, biga, bigb, bstar, n, y)
y = eqlbrm_f_g(x, astar, biga, bigb, bstar, n)
end

function deriv1_1(astar, biga, bigb, bstar, n, x)
call deriv1_1_adf1(x, astar, biga, bigb, bstar, n, y)
end adf(deriv1_1 = tangent(y))
end
```
Valid Farfel program: html text

12. By the same mechanism, we transform
```      function deriv1_2(astar, bigb, x)
y = eqlbrm_f_g_h(x, astar, bigb)
end adf(deriv1_2 = tangent(y))
end
```
into
```      subroutine deriv1_2_adf2(x, astar, bigb, y)
y = eqlbrm_f_g_h(x, astar, bigb)
end

function deriv1_2(astar, bigb, x)
call deriv1_2_adf2(x, astar, bigb, y)
end adf(deriv1_2 = tangent(y))
end
```
Valid Farfel program: html text

13. and
```      subroutine deriv2_1_1(astar, biga, bigb, bstar, n, x, y, yprime)
y = argmax_fprime_1(x, astar, biga, bigb, bstar, n)
end adf(yprime = tangent(y))
end
```
into
```      subroutine deriv2_1_1_adf3(x, astar, biga, bigb, bstar, n, y)
y = argmax_fprime_1(x, astar, biga, bigb, bstar, n)
end

subroutine deriv2_1_1(astar, biga, bigb, bstar, n, x, y, yprime)
call deriv2_1_1_adf3(x, astar, biga, bigb, bstar, n, y)
end adf(yprime = tangent(y))
end
```
Valid Farfel program: html text

14. and
```      subroutine deriv2_1_2(astar, bigb, x, y, yprime)
y = argmax_fprime_2(x, astar, bigb)
end adf(yprime = tangent(y))
end
```
into
```      subroutine deriv2_1_2_adf4(x, astar, bigb, y)
y = argmax_fprime_2(x, astar, bigb)
end

subroutine deriv2_1_2(astar, bigb, x, y, yprime)
call deriv2_1_2_adf4(x, astar, bigb, y)
end adf(yprime = tangent(y))
end
```
Valid Farfel program: html text

15. and finally
```      subroutine deriv2_2(biga, bigb, bstar, n, x, y, yprime)
y = f(x, biga, bigb, bstar, n)
end adf(yprime = tangent(y))
end
```
into
```      subroutine deriv2_2_adf5(x, biga, bigb, bstar, n, y)
y = f(x, biga, bigb, bstar, n)
end

subroutine deriv2_2(biga, bigb, bstar, n, x, y, yprime)
call deriv2_2_adf5(x, biga, bigb, bstar, n, y)
end adf(yprime = tangent(y))
end
```
Valid Farfel program: html text

16. Now we transform canonicalized `adf` blocks into calls to subroutines that will be generated by AD. For Tapenade, that looks like the following:
1. Transform
```      function deriv1_1(astar, biga, bigb, bstar, n, x)
call deriv1_1_adf1(x, astar, biga, bigb, bstar, n, y)
end adf(deriv1_1 = tangent(y))
end
```
into
```      function deriv1_1(astar, biga, bigb, bstar, n, x)
x_g1 = 1.0
astar_g1 = 0.0
bstar_g1 = 0.0
call deriv1_1_adf1_g1(x, x_g1, astar, astar_g1, biga, bigb, bstar,
+bstar_g1, n, y, deriv1_1)
end
```
Intermediate state: html text

2. Transform
```      function deriv1_2(astar, bigb, x)
call deriv1_2_adf2(x, astar, bigb, y)
end adf(deriv1_2 = tangent(y))
end
```
into
```      function deriv1_2(astar, bigb, x)
x_g2 = 1.0
astar_g2 = 0.0
call deriv1_2_adf2_g2(x, x_g2, astar, astar_g2, bigb, y, deriv1_2)
end
```
Intermediate state: html text

3. Transform
```      subroutine deriv2_1_1(astar, biga, bigb, bstar, n, x, y, yprime)
call deriv2_1_1_adf3(x, astar, biga, bigb, bstar, n, y)
end adf(yprime = tangent(y))
end
```
into
```      subroutine deriv2_1_1(astar, biga, bigb, bstar, n, x, y, yprime)
x_g3 = 1.0
astar_g3 = 0.0
bstar_g3 = 0.0
call deriv2_1_1_adf3_g3(x, x_g3, astar, astar_g3, biga, bigb, bstar,
+bstar_g3, n, y, yprime)
end
```
Intermediate state: html text

4. Transform
```      subroutine deriv2_1_2(astar, bigb, x, y, yprime)
call deriv2_1_2_adf4(x, astar, bigb, y)
end adf(yprime = tangent(y))
end
```
into
```      subroutine deriv2_1_2(astar, bigb, x, y, yprime)
x_g4 = 1.0
astar_g4 = 0.0
call deriv2_1_2_adf4_g4(x, x_g4, astar, astar_g4, bigb, y, yprime)
end
```
Intermediate state: html text

5. And transform
```      subroutine deriv2_2(biga, bigb, bstar, n, x, y, yprime)
call deriv2_2_adf5(x, biga, bigb, bstar, n, y)
end adf(yprime = tangent(y))
end
```
into
```      subroutine deriv2_2(biga, bigb, bstar, n, x, y, yprime)
x_g5 = 1.0
bstar_g5 = 1.0
call deriv2_2_adf5_g5(x, x_g5, biga, bigb, bstar, bstar_g5, n, y, yprime)
end
```
Intermediate state: html text

17. Now we propagate the constant functions `gmbiga` and `gmbigb` from the call to eqlbrm in the main program to everything it (indirectly) calls: deriv1_1_adf1, deriv1_1, deriv1_2_adf2, deriv1_2, deriv2_1_1_adf3, deriv2_1_1, deriv2_1_2_adf4, deriv2_1_2, deriv2_2_adf5, deriv2_2, root_1_1, root_1_2, root_2, argmax_fprime_1, argmax_fprime_2, argmax_1, argmax_2, eqlbrm_f_g_h, eqlbrm_f_g, eqlbrm_f, and eqlbrm, thereby freeing all programs to be differentiated from external declarations.

Valid Fortran77 program ready for AD with Tapenade: html text

Finally, we invoke Tapenade to perform the actual AD:
```tapenade -root deriv1_2_adf2   -d -o eqlbrm42 -diffvarname _g2 -difffuncname _g2 eqlbrm42.f
tapenade -root deriv2_1_2_adf4 -d -o eqlbrm42 -diffvarname _g4 -difffuncname _g4 eqlbrm42{,_g2}.f
tapenade -root deriv1_1_adf1   -d -o eqlbrm42 -diffvarname _g1 -difffuncname _g1 eqlbrm42{,_g2,_g4}.f
tapenade -root deriv2_1_1_adf3 -d -o eqlbrm42 -diffvarname _g3 -difffuncname _g3 eqlbrm42{,_g2,_g4,_g1}.f
tapenade -root deriv2_2_adf5   -d -o eqlbrm42 -diffvarname _g5 -difffuncname _g5 eqlbrm42{,_g2,_g4,_g1,_g3}.f
```
Note that the AD transforms must be invoked in the reverse of the call graph order, so that appropriate derivatives are available to be further transformed.

This work was supported, in part, by Science Foundation Ireland grant 09/IN.1/I2637, National Science Foundation grant CCF-0438806, the Naval Research Laboratory under Contract Number N00173-10-1-G023, and the Army Research Laboratory accomplished under Cooperative Agreement Number W911NF-10-2-0060. Any views, opinions, findings, conclusions, or recommendations contained or expressed in this document or material are those of the author(s) and do not necessarily reflect or represent the views or official policies, either expressed or implied, of SFI, NSF, NRL, the Office of Naval Research, ARL, or the Irish or U.S. Governments. The U.S. Government is authorized to reproduce and distribute reprints for Government purposes, notwithstanding any copyright notation herein.