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.

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

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

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 parameterf:
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

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

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

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

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

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

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

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

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

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

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

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

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

into
function deriv1_1(astar, biga, bigb, bstar, n, x)
y = eqlbrm_f_g(x, astar, biga, bigb, bstar, n)
end
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

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

into
y = eqlbrm_f_g_h(x, astar, bigb)
end

function deriv1_2(astar, bigb, x)
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

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

Valid Farfel program: html text

14. and
subroutine deriv2_1_2(astar, bigb, x, y, yprime)
y = argmax_fprime_2(x, astar, bigb)
end

into
y = argmax_fprime_2(x, astar, bigb)
end

subroutine deriv2_1_2(astar, bigb, x, y, yprime)
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

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

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

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)
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

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)
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

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.