User guide for JAKEF
Kenneth E. Hillstrom
Argonne National Laboratory
ABSTRACT
JAKEF is a language processor that accepts as input a sin-
gle precision or double precision FORTRAN subroutine defining
an objective function f(x) or a vector function F(x) and pro-
duces as output a single precision or double precision FORTRAN
subroutine defining the gradient of f(x) or the Jacobian of F(x).
JAKEF accepts as input ANSI 1977 standard FORTRAN subrou-
tines and produces as output ANSI 1977 FORTRAN subroutines.
JAKEF itself is written in ANSI 1977 standard FORTRAN, except
for the use of an extended character set. This report serves as
a user guide for the processor.
1. Introduction
A frequent plea received from people using nonlinear optimization
algorithms is for methods that are "derivative free", i.e., that do not
require them to provide derivatives. This is understandable since taking
derivatives is a tedious and error prone task. Our response to date has
been to include "derivative free" as well as derivative routines in the
first edition of the MINPACK[5] package of nonlinear optimization routines
developed in the Mathematics and Computer Science Division of Argonne
National Laboratory.
These so called "derivative free" routines have not dispensed with
taking derivatives but approximate them instead using differencing schemes
that sample the function at neighboring points. However, the method of
selecting these points may result in a poor approximation slowing the pro-
gress of the optimization routine toward a solution or even causing it to
fail. Also, differencing schemes are profligate in their function evalua-
tions. For example, computing the gradient of a function of n variables
requires n+1 function evaluations, whereas a hand-coded version, taking
---------------------------------------------------------------------------
Mathematics and Computer Science Division, Argonne Na-
tional Laboratory, Argonne, Illinois 60439. Work sup-
ported in part by the Applied Mathematical Sciences
Research Program (KC-04-02) of the Office of Energy
Research of the U.S. Department of Energy under Con-
tract W-31-109-Eng-38.
advantage of the considerable redundancy between computations of different
components of the gradient, usually requires far less computational effort.
Finally, numerical differencing regards the function as a black box, a
monolithic entity blind to the structure of the function.
Thus we would like to provide users with a method that would
relieve them of the need to code gradient and Jacobian routines and yet
take advantage of the more robust derivative optimization routines.
A technique that may fulfill the above requirements, called symbolic
or automatic differentiation, is implemented by JAKEF, a language processor
that accepts as input a FORTRAN subroutine defining a function and produces
as output a FORTRAN subroutine that computes the first derivatives of the
function.
A description of the JAKEF package is given in sections 2 through 10.
This is followed by a description of automatic differentiation, the algo-
rithm differentiator JAKE and the portable FORTRAN-77 translation of JAKE
into JAKEF.
2. Description of the JAKEF package
The JAKEF package of programs consists of interface subroutines JAKED1
and JAKEF1, the core JAKEF language processor package and run-time sin-
gle precision and double precision support packages. The user must supply
a driver in addition to the FORTRAN subroutine to be processed. In the
following sections the user is given
The rules and restrictions governing the preparation of the input
subroutine, illustrated by examples,
Instructions on how to use JAKEF,
The role of the support packages in running JAKEF-produced routines,
The JAKEF compilation and run-time error messages,
The storage requirements of JAKEF.
3. The JAKEF Input: How to Prepare Your Program
The JAKEF input consists of a single ANSI 1977 FORTRAN subroutine to
which a CONSTRUCT statement has been added and which is subject to certain
restrictions detailed in the next section.
In order to spare JAKEF and yourself considerable grief, the subrou-
tine should be successfully compiled by your FORTRAN compiler before sub-
mission to JAKEF for processing.
The CONSTRUCT statement informs JAKEF whether a one or two dimensional
array of partial derivatives, i.e., gradient or Jacobian, is requested and
where and in what pattern it is to be stored. It is the dimension of the
dependent variable vector that determines the number of rows in the array
(1 for a gradient); it is the dimension of the independent variable vector
that determines the number of columns in the array. Thus the shape of the
array is determined by the dependent and independent variables; the
declared dimensions of the array receiving the result must be compatible
with this. Note that JAKEF does not permit subscript arithmetic in dimen-
sioning the array. Finally, the array must be typed appropriately.
The CONSTRUCT statement, like the FORTRAN specification statements,
must be located after the SUBROUTINE statement and before the first execut-
able statement. The keyword 'CONSTRUCT' must be positioned in columns 1
through 9. The user is cautioned that any comment statement with the char-
acter string 'CONSTRUCT' in columns 1 through 9 will be incorrectly inter-
preted by JAKEF as a CONSTRUCT statement, triggering an error return.
Finally, the dependent and independent variables referred to in the CON-
STRUCT statement must be enclosed in parentheses.
The CONSTRUCT statement should be checked to make sure that the depen-
dent and independent variables referred to are the same, respectively, as
those in the subroutine parameter list. Any discrepancy will result in a
false translation but without any error indication.
The following examples illustrate different types of CONSTRUCT state-
ments.
SUBROUTINE ALPHA(N,X,Y,COEF)
INTEGER N
REAL Y,COEF
REAL X(N)
CONSTRUCT D(Y)/D(X) IN GRAD(N)
REAL GRAD
INTEGER I
Y = COEF
DO 10 I = 1, N
Y = Y*X(I)
10 CONTINUE
RETURN
END
Here Y is computed as a function of X. The CONSTRUCT statement informs
JAKEF that the gradient of Y as a function of X is desired, and that it is
to be stored in the array GRAD as follows:
GRAD(i) = DY/DX(i), i = 1,2,...,N.
By themselves, the variable names N, X, Y, COEF have no special meaning.
Also, as is shown in the following example, neither the name of the depen-
dent variable, the number or order or precision of the parameters in the
parameter list or even the form of the independent variables is presupposed
by JAKEF. However, the gradient must be named GRAD.
SUBROUTINE BETA(Q,P,R,PI)
DOUBLE PRECISION Q,P,R,PI
CONSTRUCT D(Q)/D(P,R) IN GRAD(2)
DOUBLE PRECISION GRAD
DOUBLE PRECISION T
T = DSIN(P*PI/4.0D0)
Q = DCOS(R)/T*R
RETURN
END
This CONSTRUCT statement will cause JAKEF to define
GRAD(1) = DQ/DP, GRAD(2) = DQ/DR,
with PI regarded as a constant. The next example introduces Jacobians with
varying dimensions.
SUBROUTINE GAMMA(M,N,X,F)
INTEGER M,N
DOUBLE PRECISION X(N),F(M)
CONSTRUCT D(F)/D(X) IN JACOB(M,N)
DOUBLE PRECISION JACOB
INTEGER I,J
DO 20 I = 1, M
F(I) = 0.0D0
DO 10 J = 1, N
F(I) = F(I) + X(J)**I
10 CONTINUE
20 CONTINUE
RETURN
END
This CONSTRUCT statement asks for the partial derivatives to be computed
and stored in an array, which must be named JACOB, as follows:
JACOB(i,j) = DF(i)/DX(j), i = 1,2,...,M; j = 1,2,...,N.
The next example introduces Jacobians with fixed dimensions.
SUBROUTINE DELTA(X,F)
REAL X(3),F(4)
CONSTRUCT D(F)/D(X) IN JACOB(4,3)
REAL JACOB
F(1) = SIN(X(1)+X(2)) + 10.0E0*X(2)
F(2) = (X(3) - X(2))**2
F(3) = (X(2) - 2.0E0*X(3))**2
F(4) = EXP(X(1) - X(3))
RETURN
END
This CONSTRUCT statement asks for the partial derivatives to be computed
and stored in an array, which must be named JACOB, as follows:
JACOB(i,j) = DF(i)/DX(j), i = 1,2,3,4; j = 1,2,3.
The following example introduces multiple Jacobians.
SUBROUTINE THETA(N,U,V,W)
INTEGER N
REAL U(N),V(N),W(N)
CONSTRUCT D(W)/D(U,V) IN JACOB(N,100)
REAL JACOB
INTEGER I
DO 10 I = 1, N
W(I) = U(I)*V(N-I+1)
10 CONTINUE
RETURN
END
This CONSTRUCT statement asks for the Jacobians DW/DU and DW/DV to be
computed and stored in JACOB as follows:
JACOB(i,j) = DW(i)/DU(j), i = 1,...,N; j = 1,...,N;
JACOB(i,N+j) = DW(i)/DV(j), i = 1,...,N; j = 1,...,N.
Since JAKEF does not permit subscript arithmetic in the derivative arrays,
the dimensions of the array JACOB were defined as (N,100) instead of
(N,2*N), imposing a limit of 50 on N.
In contrast, the following example produces a gradient.
SUBROUTINE SIGMA(M,N,V,W)
INTEGER M,N
REAL W
REAL V(M,N)
CONSTRUCT D(W)/D(V) IN GRAD(1000)
REAL GRAD
INTEGER I,J
W = 1.0
DO 20 J = 1, N
DO 10 I = 1, M
W = W*V(I,J)
10 CONTINUE
20 CONTINUE
RETURN
END
The storage pattern is
GRAD(I+(J-1)*M) = DW/DV(i,j), i = 1,...,M; j = 1,...,N.
So a matrix such as V(M,N) is really treated as the one dimensional array
of M*N contiguous storage locations it represents. Again, note that the
subscript arithmetic restriction does not permit a dimension specification
of M*N.
Next consider the subroutine OMEGA which introduces mixed precisions.
SUBROUTINE OMEGA(N,U,V,W)
INTEGER N
REAL W
REAL U(N),V(N)
CONSTRUCT D(W)/D(U) IN GRAD(N)
DOUBLE PRECISION GRAD
INTEGER I
W = 1.0
DO 10 I = 1, N
W = W*U(I) + V(I)
10 CONTINUE
RETURN
END
By indicating that the result vector GRAD is double precision, one asks
JAKEF to emit all derivative factors in double precision, to perform multi-
plication of factors in double precision and to extract the result in dou-
ble precision. At the same time, all variables that were single precision
in the original program remain single precision. Some intermediate values,
such as W*U(I), are computed in single precision in the original program
but in double precision when processed by JAKEF. This is a compromise among
accuracy, speed and ease of handling by JAKEF. Conversely, it is possible
to produce a single precision gradient from an otherwise double precision
computation.
4. Restrictions on the Input Subroutine.
The purpose of this section is to list the limitations of JAKEF with
regard to the input program. The restrictions are:
JAKEF does not correctly handle statement functions, and no warning is
issued. Thus they must not be used.
JAKEF does not recognize variables, constants or functions of type
complex. The type designation COMPLEX is treated as an illegal key-
word, triggering an error return.
The input program must not contain any of the following subroutine
names, as they are reserved for the members of the single precision
and double precision run-time support packages associated with JAKEF:
SPINIT EMIT0 EMIT1 EMIT2 SPGRAD SPCOPY SERRM
DPINIT DMIT0 DMIT1 DMIT2 DPGRAD DPCOPY DERRM
JAKEF does not correctly handle EQUIVALENCEs that involve the indepen-
dent variables. More precisely, EQUIVALENCEs that violate this res-
triction are misinterpreted without warning by the JAKEF differentia-
tor.
Dimensioning of the dependent and independent variables must occur in
a type statement and must conform to FORTRAN-66 standards. All other
arrays, dimensioned in a DIMENSION statement, follow FORTRAN-77
standards.
Variables with dimensions determined by a parameter statement must be
dimensioned in a DIMENSION statement. Also, the PARAMETER statement
must preceed the DIMENSION statement.
Statement numbers 70000 - 100000 are reserved for usage by JAKEF.
The length specifier "len" in a CHARACTER statement must be an
unsigned nonzero integer constant or an asterisk enclosed in
parenthesis. The entity specifier "nam" in a CHARACTER statement must
be a variable, function or array name.
Character substring expressions are not permitted.
Alternate returns are not permitted.
JAKEF does not correctly handle calls and function references that
involve the independent variables unless the function is one of the
following: SIN/DSIN, COS/DCOS, ATAN/DATAN, EXP/DEXP, ALOG/DLOG,
SQRT/DSQRT, TAN/DTAN, ASIN/DASIN, ABS/DABS, SINH/DSINH and COSH/DCOSH.
Calls and references to other functions are misinterpreted without
warning by the JAKEF differentiator. Note also that except for those
in the above list JAKEF declares with default typing all functions not
explicitly declared in the input program.
5. Calling JAKEF
This section contains the instructions for writing and compiling the
FORTRAN subroutine that defines the function and using the JAKEF-produced
derivative subroutine.
The process of differentiating a subroutine with JAKEF requires that
the subroutine which serves as data for JAKEF be programmed as described in
Section 3, and that a main program be written to communicate to JAKEF.
The easiest way to communicate to JAKEF is via the interface routine
JAKED1. The following program illustrates the use of JAKED1.
C
C SAMPLE PROGRAM.
C
INTEGER IWA(7600)
CHARACTER*1 CWA(4500)
INTEGER LCWA,LIWA
LIWA = 7600
LCWA = 4500
CALL JAKED1(IWA,LIWA,CWA,LCWA)
STOP
END
The JAKED1 parameters are defined as follows:
IWA is an integer input array of length LIWA.
LIWA is an integer input variable.
CWA is a character*1 input array of length LCWA.
LCWA is an integer input variable not smaller than 4500.
It is usually not possible to compute a priori precise dimensions for the
arrays IWA and CWA. The above example sets LIWA to 7600 because this set-
ting has been found to be adequate for most problems. LCWA must always be
set to at least 4500. If the dimensions of IWA or CWA are inadequate an
error message will be sent to the output. The user should then consult
the compile-time error messages listed in Section 8 to see which of the
dimensions of IWA or CWA should be increased.
Most of the storage required by JAKEF is for the differentiation pro-
cess. If the subroutine is of the size and complexity of subroutines ALPHA
through OMEGA listed in Section 3 then JAKEF will not require much storage.
In particular, these subroutines can be successfully processed by JAKED1 if
IWA has dimension 6100 and CWA has the (minimum) dimension 4500. More com-
plicated subroutines require additional space. For example, the WATSN sub-
routine listed in Appendix C is successfully processed by JAKED1 if IWA has
dimension 7300 and CWA has dimension 4500.
If the user desires complete control of storage specification in
JAKEF, he must write a main program to call JAKEF1, a subroutine which par-
titions work storage arrays according to its input parameters. In order to
do this he should consult the listing of the interface routine JAKED1 and
the explanation of its JAKEF1 argument settings in Appendix A. The JAKEF1
listing itself is located in Appendix B.
6. Using the JAKEF-produced derivative subroutine
The JAKEF output is a FORTRAN subroutine adhering to the ANSI 1966
standard as much as the input subroutine does; non-standard constructs in
DATA and FORMAT statements that are accepted by the FORTRAN compiler pass
through JAKEF processing without change. The output FORTRAN derivative
subroutine cannot be readily compared with a handwritten counterpart since
it consists largely of calls to a run-time support package of routines that
implements the differentiation process.
Certain errors (for example, an improper CONSTRUCT statement) will
produce incorrect FORTRAN syntax in the JAKEF output. Thus it is advisable
to check the syntax by compiling the JAKEF output.
It should be kept in mind that the JAKEF process of adjoining dif-
ferentiation operations at appropriate points of the user-coded input sub-
routine results in an output subroutine that simultaneously evaluates both
the objective function and its derivatives. When, as is often the case,
both are required in an optimization problem, a single call to the JAKEF-
produced subroutine suffices. On the other hand, an untimely function call
would have to be redirected by the user, through the calling sequence, to
dispose of the evaluated function innocuously.
Using the JAKEF output requires a driver with a calling sequence
matched to the parameter list of the JAKEF output subroutine. In particular
certain auxiliary arrays must be provided.
Consider gradient codes. The input subroutine ALPHA listed in Section
3 would be translated into the gradient subroutine ALPHAJ listed, in part,
below.
SUBROUTINE ALPHAJ(N,X,Y,COEF,GRAD,
* YGRAD,LYGRAD,RFS,IFS,LFS)
INTEGER LFS,IFS(LFS)
REAL RFS(LFS),TGRA(10)
INTEGER N,I,LYGRAD,IGRAD,RGRAD,IX
REAL X(N),Y,COEF,GRAD(N),YGRAD(LYGRAD)
IX = 11
CALL SPINIT(IX+N,LYGRAD)
.
.
END
Note that the name assigned by JAKEF to the gradient code is the input sub-
routine name with a J appended to it. If this name has more than five char-
acters JAKEF will truncate the name after the fifth character and append a
J.
In addition to N, X, Y, COEF and the array GRAD, subroutine ALPHAJ
requires the calling program to supply the work array YGRAD with dimension
LYGRAD, and the arrays RFS and IFS each with dimension LFS.
A Jacobian subroutine requires similar auxiliary arrays. For example,
the GAMMA subroutine listed in Section 3 is translated into the subroutine
GAMMAJ listed, in part, below.
SUBROUTINE GAMMAJ(M,N,X,F,JACOB,LJAC,
* YJACOB,LYJACO,RFS,IFS,LFS)
INTEGER LJAC
INTEGER LFS,IFS(LFS)
DOUBLE PRECISION RFS(LFS),TJACO(10)
INTEGER M,N,LYJACO,IF,IX
DOUBLE PRECISION X(N),F(M),JACOB(LJAC,N),YJACOB(LYJACO)
IF = 10
IX = IF+M
CALL DPINIT(IX+N,LYJACO)
.
.
END
In addition to M, N, X, F and the two-dimensional array JACOB, subroutine
GAMMAJ requires LJAC to be set to the leading dimension of the array JACOB,
the work array YJACOB to be supplied with dimension LYJACO and the arrays
RFS and IFS each with dimension LFS.
SPINIT, called from ALPHAJ, and DPINIT, called from GAMMAJ, check that
their second argument is at least equal to the first. Thus for ALPHAJ,
LYGRAD must be set at least equal to 11 + N since IX = 11. For GAMMAJ,
LYJACO must be set at least equal to 10 + M + N since IX = 10 + M. If
LYGRAD or LYJACO is too small, the run-time error message "VECTOR YGRAD OR
YJACOB IS TOO SMALL" is triggered in SPINIT or DPINIT. The setting of the
dimension parameter and the dimension of the corresponding array must then
be increased before another run is attempted.
RFS and IFS are work arrays required by other members of the run-time
support package invoked by JAKEF-produced routines. If these arrays are too
small the error message "OVERFLOW OF BUFFER" is written, execution is ter-
minated and the user must increase the setting of LFS and the corresponding
dimensions of RFS and IFS. See Section 7 for a description of the run-time
support packages and see Section 8 for a listing of the run-time error
messages.
Tests were conducted in order to determine reasonable estimates of LFS
settings. In the results we have noted the smallest value of LFS which per-
mits error-free execution of the JAKEF-generated code. For subroutine
ALPHAJ the results are as follows:
N LFS
10 31
20 61
30 91
Note that LFS depends linearly on N, and that LFS is of modest size for
reasonable values of N. Much larger settings are necessary, however, when
running programs containing do-loops, because the amount of space required
is directly proportional to the number of statements executed in a trace of
the run. This is illustrated by the JAKEF-produced Watson function gradient
code listed in Appendix C. This code contains two inner N-pass do-loops
within a 29 pass outer loop. The results are as follows:
N LFS
10 3179
20 6079
40 11879
Note that LFS again depends linearly on N but that the additional complex-
ity of the Watson function leads to a substantial increase in its size.
7. The Run-time Support Packages
The run-time support package consisting of the seven routines EMIT0..SERRM
or DMIT0..DERRM listed above is is required during the execution of
JAKEF-produced routines to implement the differentiation process. Both packages
are similar, so only the double precision package will be described. The members
of this package and their functions are as follows:
DPINIT - Initializes integer and double precision factor storage
counters and checks to see that the dimension LYGRAD or LYJACO of the
work vector YGRAD or YJACOB required by the JAKEF-generated code is at
least IX + N. If the dimension is too small the message "VECTOR YGRAD
OR YJACOB IS TOO SMALL" is printed and execution is terminated.
DMIT0, DMIT1 and DMIT2 record the derivatives of functions of a con-
stant, of one, and of two independent variables, respectively.
DPGRAD - Performs the factor or matrix multiplications that generate
the vector of partial derivatives in YGRAD.
DPCOPY - Copies the gradient requested from N locations in YGRAD,
beginning at IX + 1, to the single row of GRAD or to one of the rows
of the two-dimensional array JACOB.
DERRM - Writes an error message and aborts.
A listing of the double precision support package is given in Appendix D.
8. JAKEF Error Messages
Error messages from JAKEF can occur at compile time or at run time.
Any of these messages announces a fatal error which terminates processing.
A listing of the messages and suggested corrective actions can be found in
Sub-sections 8.1 and 8.2 below.
8.1. Compile-time Error Messages
Missing CONSTRUCT statement.
CONSTRUCT statement error.
Unrecognized keyword - FORTRAN input program error. Check FORTRAN source.
FORTRAN statement exceeds character limit. Increase LCWA.
FORTRAN statement array insufficiency - LCWA has been set less than its
minimum value 4500. Increase LCWA to at least 4500.
Lexical preprocessing - No continuation card expected here. FORTRAN input
program error. Check FORTRAN source.
Overflow of name table index space - Number of names exceeds name table
index space. Increase LIWA.
Overflow of name table space - Number of characters exceeds name table
space. Increase LCWA.
Syntax error in parsing phase - FORTRAN syntax error in input program.
Check FORTRAN source.
Parsing stack overflow - Increase LIWA.
Overflow of carry space - Number of characters exceeds carry space.
Increase LCWA.
Overflow of parse buffer - Increase LIWA.
Unknown token in maketree - FORTRAN input program error. Check FORTRAN
source.
Nested do-loop limit exceeded - Increase LIWA.
Parse tree has wrong structure - FORTRAN input program error. Check FORTRAN
source.
Wrong successor structure - FORTRAN input program error. Check FORTRAN
source.
Wrong go-to structure - FORTRAN input program error. Check FORTRAN source.
Wrong floating point constant - Digit outside of 0-9 range. Check FORTRAN
source.
Overflow of expression space - Too many expressions in expression tree.
Increase LIWA.
Save space overflow in alpha - Increase LIWA.
Variable limit exceeded - Increase LIWA.
Label limit exceeded - Increase LIWA.
Right-hand-side space overflow - Increase LIWA.
Overflow of statement space - Increase LIWA.
Overflow of flag space - Increase LCWA.
Subscripted scalar encountered - FORTRAN input program error. Check FORTRAN
source.
CONSTRUCT statement variable should be typed real or double precision.
TLIM is not large enough - More temporary variables are required. Increase
LCWA.
More than 999 temporaries required - FORTRAN input program is too large to
be processed by JAKEF.
Unrecognized unary operation - FORTRAN input program error. Check FORTRAN
source.
Wrong binary operation - FORTRAN input program error. Check FORTRAN source.
8.2. Run-time Error Messages
Overflow of buffer - Increase setting of LFS and the corresponding dimen-
sions of RFS and IFS.
Vector YGRAD or YJACOB is too small - Increase setting of LYGRAD or LYJACO
and the corresponding dimension of YGRAD or YJACOB.
9. Accuracy and Efficiency of JAKEF-Generated Codes
Test results show that the accuracy of JAKEF-generated codes is essen-
tially the same as that obtained using hand-coded counterparts, and much
better than the accuracy of difference approximations.
Test results also show that JAKEF-generated codes are markedly more
efficient than difference approximations when computing gradients and much
less efficient when computing Jacobians. When compared with hand-coded
versions, JAKEF-produced codes are always less efficient; this is primarily
due to the large number of bookkeeping operations required by the run-time
support routines.
10. JAKEF Storage Requirements
There are 5997 lines of comments and code in JAKEF and an additional
144 in each of the run-time support packages. The array storage usage of
the support packages derives from the LFS value set by the user in his
driver.
11. The Automatic Differentiation Technique
The basic idea in automatic differentiation is rather simple: given a
function code, the rules for differentiation are applied line for line to
give a list of instructions for evaluating the derivative. Details about
and examples illustrating automatic differentiation are given in the follow-
ing section.
Sophisticated systems for symbolic algebraic manipulation such as
MACSYMA and FORMAC offer facilities for differentiation of formulas. How-
ever, these systems have not been viewed as adequate solutions to the prob-
lem of taking derivatives for the following reasons:
Firstly, algebraic manipulation systems deal with formulas in closed
form; they are not designed to handle functions necessarily defined by
an algorithm. In particular, these systems cannot handle do-loops and
flow-control statements.
Secondly, a hand compilation and translation is required to get the
output from these systems into computer-executable form.
Finally, these systems by their very nature cannot optimize entire
gradient computations; they can only simplify individual derivatives,
Evaluating the gradient still takes O(n) function evaluations.
What is desired is a system that will accept the text of an algorithm A
embodying an objective function f(x) or a vector function F(x) and con-
struct from it the text of an algorithm A' that computes the gradient of
f(x) or the Jacobian of F(x). Two algorithm differentiators, AUGMENT[1] and
JAKE[2], were investigated. JAKE was selected because of its versatility
and the fact that it uses the computationally more efficient reverse mode
of automatic differentiation. For a recent survey of various techniques for
obtaining derivative values see [xx].
12. The JAKE Algorithm Differentiator
JAKE is a compiler that accepts as data a single precision or double
precision FORTRAN subroutine defining an objective function f(x) or a vec-
tor function F(x) and produces as output a single precision or double pre-
cision FORTRAN subroutine defining the gradient of f(x) or the Jacobian of
F(x).
12.1 The Algorithm Implementation JAKE
JAKE is a four-pass compiler consisting of
A lexical preprocessor,
A parser,
A tree building and flow analysis pass,
A differentiator and output construction pass.
The lexical preprocessor of JAKE reworks the input FORTRAN program to give
it a recognizable lexical structure.
The parser transforms the preprocessed input into a string of tokens in a
postfix representation of the program tree.
The tree building and flow analysis pass constructs a tree out of the post-
fix token string.
The differentiator identifies relevant assignment statements, that is,
statements involving the independent variables. Then, if necessary, it
analyzes them into component statements governed by a single differentia-
tion rule and augments each of these statements with a call to a member of
the run-time support package which implements the differentiation rule.
After completing the construction of the main body of the routine, JAKE
inserts calls to support package routines that complete the
differentiation. This results in a modified program tree which is then
printed in a form compatible with FORTRAN rules.
13. JAKEF - FORTRAN-77 Version of JAKE
As stated previously, the input and output language of JAKE is FORTRAN.
However, JAKE itself was written in C, a language well suited for writing
compilers.
Since our interest in JAKE was largely motivated by the desire to pro-
vide MINPACK users with an automatic method of obtaining derivatives, we
required a portable version of JAKE. Thus, for practical reasons, we needed
a version of JAKE written in a FORTRAN dialect, and the extensive character
manipulation features embodied in JAKE dictated a FORTRAN-77 version. Thus
JAKEF, a FORTRAN-77 version of JAKE, was written. JAKEF is portable in the
sense that it can be used without alteration at any installation that
recognizes an extended character set that includes lower case alphabetics
and the following characters:
? @ " # % & ~ ^ | _ [ ; < >.
The C to FORTRAN-77 transformation of JAKE was not completely
straight-forward; some innovations were required. In addition, there are
some features of JAKEF that may be of separate interest to the user and
thus both are listed below as follows:
The reworked FORTRAN input program described in Section 12 is put on
scratch file RFORTD created by JAKEF and is still available after
JAKEF processing.
The recursive calls to subroutines used in JAKE are not supported by
ANSI FORTRAN-77 and thus are simulated in JAKEF.
The C language parser of JAKE was automatically coded by YACC(Yet
Another Compiler Compiler) in an ASCII environment. Thus the parsing
scheme is ASCII-character-representation dependent, requiring a trans-
lation table in JAKEF to achieve portability. This table is automati-
cally constructed by interrogating the particular environment to
obtain the number representation of each character and associating it
with its ASCII counterpart.
Various C language bit manipulation features are simulated in JAKEF.
The benefit/cost tradeoff of the JAKEF translation is a portable compiler,
free of the confines of a specialized language, but slower in operation.
14. Summary
JAKEF produces derivative codes that are reliable and accurate - as
accurate as hand-coded counterparts and more accurate than difference
approximations. Timing tests indicate that the gradient codes are more
efficient and the Jacobian codes less efficient than differencing schemes.
Translators such as JAKEF are certainly labor saving devices, espe-
cially useful at the laboratory or investigative level where limited runs
rather than production runs are the rule. In addition, such translators are
indispensable tools in solving problems where analytic derivative accuracy
is required but hand coding of derivatives is difficult.
15. Acknowledgments
The author wishes to thank Burton S. Garbow, Jorge J. More', Gary A.
Roediger and Brian T. Smith for the many helpful suggestions in designing
and documenting JAKEF that have contributed to its robustness and ease of
use.
16. References
1. Kedem, G., Automatic Differentiation of Computer Programs, ACM
Transactions on Mathematical Software, Vol. 6, No. 2, June 1980,
Pages 150-165.
2. Speelpenning, B., Compiling Fast Partial Derivatives of Functions
given by Algorithms, UILU-ENG 80 1702, University of Illinois-Urbana,
January 1980.
3. Rall, Louis B., Automatic Differentiation: Techniques and Applica-
tions, Lecture Notes in Computer Science, 120, Springer-Verlag,
Berlin, Heidelberg, New York, 1981.
4. Warner, D. D., A Partial Derivative Generator, C.S. Tech Report,
Bell Laboratories, April 1975.
5. More', Jorge J., Garbow, Burton S., and Hillstrom, Kenneth E.,
User Guide for MINPACK-1, ANL-80-74, Argonne National Laboratory,
6. Griewank, A., 'On Automatic Differentiation',
Preprint MCS-P10-1088, Mathematics and Computer Science Division
Argonne National Laboratory.
Appendix A.
SUBROUTINE JAKED1(IWA,LIWA,CWA,LCWA)
INTEGER LIWA,LCWA
INTEGER IWA(LIWA)
CHARACTER*1 CWA(LCWA)
INTEGER TLIM,LSAVE,LSTM,LVAR,LVARD,LCRYSP,LOPNDS,
* LLBL,LEXPR,LRHS,LFLAG,LNMTBL,LPFILE,LINDX,
* LSTACK,NREAD,NWRITE,RFORT
REAL DCWA,DIWA,PCWA,PIWA
DATA DCWA,DIWA /4.5E3,7.6E3/
C
C JAKEF INTERFACE ROUTINE JAKED1.
C
C THIS ROUTINE CONNECTS EXTERNAL UNITS 5 AND 6 TO THE
C STANDARD INPUT AND OUTPUT. A SCRATCH FILE REQUIRED BY
C JAKEF IS ARBITRARILY CONNECTED TO UNIT 1. OTHER
C INSTALLATIONS MAY REQUIRE THE JAKEF INSTALLER TO CHANGE
C THESE DESIGNATIONS.
C
NREAD = 5
NWRITE = 6
RFORT = 1
C
C THE FOLLOWING DIMENSION PARAMETERS ARE SET TO FRACTIONS
C OF WORK ARRAY DIMENSIONS LIWA AND LCWA.
C
PIWA = FLOAT(LIWA)/DIWA
PCWA = FLOAT(LCWA)/DCWA
TLIM = 10*PCWA
TLIM = MIN(TLIM,999)
LSAVE = 50*PIWA
LSTM = 50*PIWA
LVAR = 200*PCWA
LVARD = 200*PIWA
LCRYSP = 100*PCWA
LOPNDS = 5*PIWA
LLBL = 10*PIWA
LEXPR = 1000*PIWA
LRHS = 50*PIWA
LFLAG = 1000*PCWA
LNMTBL = 500*PCWA
LPFILE = 500*PIWA
LINDX = 50*PIWA
LSTACK = 150*PIWA
CALL JAKEF1(NREAD,NWRITE,RFORT,TLIM,LSAVE,LSTM,LVAR,LVARD,
* LCRYSP,LOPNDS,LLBL,LEXPR,LRHS,LFLAG,LNMTBL,
* LPFILE,LINDX,LSTACK,IWA,LIWA,CWA,LCWA)
RETURN
C
C LAST CARD OF INTERFACE ROUTINE JAKED1.
C
END
The parameter settings in JAKED1 must satisfy the JAKEF language require-
ments described below.
NREAD,NWRITE,RFORT - NREAD and NWRITE are set to 5 and 6, values that
are often the standard I/O settings. RFORT is arbitrarily set to 1.
The JAKEF installer may be required to change these settings according
to the dictates of his installation.
TLIM - JAKEF-produced routines are programmed to reserve TLIM tem-
poraries for intermediate storage. TLIM should be set to at least the
number of temporaries anticipated, but must be less than 1000. Since
few if any routines produced by JAKEF will require 1000 temporaries,
this is not a significant restriction.
LSAVE - Several of the JAKEF subprograms require recursion, a language
feature not included in the ANSI 1977 standard. Thus these subprograms
are coded to simulate recursive calls, requiring information to be
saved before re-entry into the program and restored upon exit. This in
turn requires storage arrays for each item of information preserved
that are of dimension at least LSAVE times as large as the depth of
recursion. The depth of recursion varies directly with the length of,
and the order of, do-loop nesting in the routine being processed.
LSTM - Should be set at least equal to the number of statements in the
routine processed.
LVAR - Should be set at least equal to the number of characters in the
longest statement in the routine processed.
LVARD - Should be set at least equal to the number of variables in the
routine generated.
LCRYSP - JAKEF transmits, without change, any DATA and FORMAT state-
ments in the routine processed. LCRYSP must be set at least as large
as the total number of characters in these statements.
LOPNDS - Should be set at least equal to the depth of do-loop nesting
in the routine processed.
LLBL - Should be set at least equal to the number of labels in the
routine generated.
The following parameters provide the space required by the parsing and tree
building functions of JAKEF; their settings are difficult for the user to
estimate. Brief descriptions of the functions associated with the storage
being dimensioned and recommended settings are given below.
LEXPR - Should be set at least equal to the number of expressions in
the expression tree constructed by the JAKEF parser. A setting of 1000
is recommended.
LRHS - Should be set at least equal to the number of right hand side
variables in the statements represented in the expression tree. A set-
ting of 50 is recommended.
LFLAG - Should be set at least equal to the number of derivative fac-
tors emitted by the JAKEF differentiator. A setting of 1000 is recom-
mended.
LNMTBL - Should be set at least equal to the number of characters con-
tained in the variable, function and subroutine names stored in the
array NMTBL. These names were obtained from the routine processed and
from the JAKEF data base. A setting of 500 is recommended.
LPFILE - Should be set at least equal to the number of tokens in the
parse buffer. A setting of 500 is recommended.
LINDX - Should be set at least equal to the number of variable, func-
tion and subroutine names in the array NMTBL. A setting of 50 is
recommended.
LSTACK - Should be set at least equal to the number of tokens in the
parsing stack. A setting of 150 is recommended.
IWA,LIWA,CWA,LCWA - All of the arrays in JAKEF governed by the above
dimensions exist as partitions of the integer array IWA or the charac-
ter array CWA.
Appendix B.
SUBROUTINE JAKEF1(NREAD,NWRITE,RFORT,TLIM,LSAVE,LSTM,LVAR,
* LVARD,LCRYSP,LOPNDS,LLBL,LEXPR,LRHS,LFLAG,
* LNMTBL,LPFILE,LINDX,LSTACK,IWA,LIWA,CWA,
* LCWA)
INTEGER NREAD,NWRITE,RFORT,TLIM,LSAVE,LSTM,LVAR,LVARD,
* LCRYSP,LOPNDS,LLBL,LEXPR,LRHS,LFLAG,LNMTBL,LPFILE,
* LINDX,LSTACK,LIWA,LCWA
INTEGER IWA(LIWA)
CHARACTER*1 CWA(LCWA)
C
C THIS ROUTINE PARTITIONS THE INTEGER AND CHARACTER WORK
C ARRAYS IWA AND CWA USING THE INPUT DIMENSION PARAMETERS.
C
INTEGER INDX1,INDX2,INDX3,INDX4,INDX5,INDX6,INDX7,INDX8,
* INDX9,INDX10,INDX11,INDX12
INDX1 = 5*LSAVE + 4
INDX2 = INDX1 + 4*LEXPR + 2
INDX3 = INDX2 + LRHS
INDX4 = INDX3 + LPFILE
INDX5 = INDX4 + 16*LSTM + 6
INDX6 = INDX5 + 8*LVARD
INDX7 = INDX6 + 2*LSTACK + 2
INDX8 = INDX7 + LOPNDS
INDX9 = INDX8 + LLBL
INDX10 = 2*LFLAG + 1
INDX11 = INDX10 + LNMTBL + 1
INDX12 = INDX11 + LCRYSP
CALL JAKEF(IWA(1),LSAVE,IWA(INDX1+1),LEXPR,IWA(INDX2+1),LRHS,
* CWA(1),LFLAG,CWA(INDX10+1),LNMTBL,CWA(INDX11+1),
* LCRYSP,IWA(INDX3+1),LPFILE,IWA(INDX4+1),LSTM,
* CWA(INDX12+1),LVAR,IWA(INDX5+1),LVARD,
* IWA(INDX6+1),LSTACK,IWA(INDX7+1),LOPNDS,
* IWA(INDx8+1),LLBL,IWA(INDX9+1),LINDX,TLIM,RFORT,
* NREAD,NWRITE)
RETURN
C
C LAST CARD OF SUBROUTINE JAKEF1.
C
END
Appendix C.
C
C THE WATSON OBJECTIVE FUNCTION
C
SUBROUTINE DWATS(N,X,F)
CONSTRUCT D(F)/D(X) IN GRAD(N)
INTEGER N
DOUBLE PRECISION X(N)
DOUBLE PRECISION F
INTEGER I,J
DOUBLE PRECISION DFLOAT,D1,D2,GRAD,S1,S2,T,T1
DOUBLE PRECISION ZERO,ONE,C29
DATA ZERO,ONE,C29 /0.0D0,1.0D0,2.9D1/
F = ZERO
DO 30 I = 1, 29
D1 = DFLOAT(I)/C29
S1 = ZERO
D2 = ONE
DO 10 J = 2, N
S1 = S1 + DFLOAT(J-1)*D2*X(J)
D2 = D1*D2
10 CONTINUE
S2 = ZERO
D2 = ONE
DO 20 J = 1, N
S2 = S2 + D2*X(J)
D2 = D1*D2
20 CONTINUE
T = S1 - S2**2 - ONE
F = F + T**2
30 CONTINUE
T1 = X(2) - X(1)**2 - ONE
F = F + X(1)**2 + T1**2
RETURN
END
C
C THE WATSON OBJECTIVE FUNCTION.
C
SUBROUTINE DWATSJ(N,X,F,GRAD,YGRAD,LYGRAD,RFS,IFS,LFS)
INTEGER LFS,IFS(LFS)
DOUBLE PRECISION RFS(LFS),TGRA(66)
INTEGER N,I,J,LQ00,LQ01,LQ02,LQ03,LQ04,LQ05,LYGRAD,IGRAD,
*RGRAD,IX
DOUBLE PRECISION X(N),F,DFLOAT,D1,D2,GRAD(N),S1,S2,T,T1,ZERO,
*ONE,C29,YGRAD(LYGRAD)
DATA ZERO,ONE,C29/0.0D0,1.0D0,2.9D1/
IX = 71
CALL DPINIT(IX+N,LYGRAD)
CALL DMIT0(1,RFS,IFS,LFS)
F = ZERO
LQ00 = 1
LQ01 = 29
DO 90001 I = LQ00,LQ01
D1 = (DFLOAT(I))/C29
CALL DMIT0(2,RFS,IFS,LFS)
S1 = ZERO
D2 = ONE
LQ02 = 2
LQ03 = N
DO 90002 J = LQ02,LQ03
TGRA(2) = (DFLOAT(J-1))*D2
CALL DMIT1(IX+J,TGRA(2),6,RFS,IFS,LFS)
TGRA(1) = TGRA(2)*X(J)
CALL DMIT2(2,1.D0,6,1.D0,2,RFS,IFS,LFS)
S1 = S1+TGRA(1)
D2 = D1*D2
90002 CONTINUE
CALL DMIT0(3,RFS,IFS,LFS)
S2 = ZERO
D2 = ONE
LQ04 = 1
LQ05 = N
DO 90003 J = LQ04,LQ05
CALL DMIT1(IX+J,D2,6,RFS,IFS,LFS)
TGRA(1) = D2*X(J)
CALL DMIT2(3,1.D0,6,1.D0,3,RFS,IFS,LFS)
S2 = S2+TGRA(1)
D2 = D1*D2
90003 CONTINUE
CALL DMIT1(3,S2+S2,7,RFS,IFS,LFS)
TGRA(2) = S2**2
CALL DMIT2(2,1.D0,7,-(1.D0),6,RFS,IFS,LFS)
TGRA(1) = S1-TGRA(2)
CALL DMIT1(6,1.D0,4,RFS,IFS,LFS)
T = TGRA(1)-ONE
CALL DMIT1(4,T+T,6,RFS,IFS,LFS)
TGRA(1) = T**2
CALL DMIT2(1,1.D0,6,1.D0,1,RFS,IFS,LFS)
F = F+TGRA(1)
90001 CONTINUE
CALL DMIT1(IX+1,X(1)+X(1),7,RFS,IFS,LFS)
TGRA(2) = X(1)**2
CALL DMIT2(IX+2,1.D0,7,-(1.D0),6,RFS,IFS,LFS)
TGRA(1) = X(2)-TGRA(2)
CALL DMIT1(6,1.D0,5,RFS,IFS,LFS)
T1 = TGRA(1)-ONE
CALL DMIT1(IX+1,X(1)+X(1),8,RFS,IFS,LFS)
TGRA(3) = X(1)**2
CALL DMIT2(1,1.D0,8,1.D0,6,RFS,IFS,LFS)
TGRA(1) = F+TGRA(3)
CALL DMIT1(5,T1+T1,7,RFS,IFS,LFS)
TGRA(2) = T1**2
CALL DMIT2(6,1.D0,7,1.D0,1,RFS,IFS,LFS)
F = TGRA(1)+TGRA(2)
90000 CONTINUE
RGRAD = 0
CALL DPGRAD(YGRAD,LYGRAD,1,RGRAD,IGRAD,RFS,IFS,LFS)
CALL DPCOPY(GRAD,IGRAD,1,YGRAD(IX+1),N)
RETURN
END
Appendix D.
C******************************************************************
C
C DOUBLE PRECISION RUN-TIME SUPPORT PACKAGE.
C
C******************************************************************
SUBROUTINE DMIT0(DEPVAR,RFS,IFS,LFS)
INTEGER DEPVAR,LFS
INTEGER IFS(LFS)
DOUBLE PRECISION RFS(LFS)
INTEGER FSP,NFS
COMMON /FACTOR/ FSP,NFS
IF (FSP .GT. LFS) CALL DERRM('OVERFLOW OF BUFFER@')
FSP = FSP + 1
IFS(FSP) = FSP - 1
RFS(FSP) = DEPVAR
RETURN
C
C LAST CARD OF SUBROUTINE DMIT0.
C
END
C******************************************************************
SUBROUTINE DMIT1(IND1,DER1,DEPVAR,RFS,IFS,LFS)
INTEGER IND1,DEPVAR,LFS
INTEGER IFS(LFS)
DOUBLE PRECISION DER1
DOUBLE PRECISION RFS(LFS)
INTEGER FSP,NEWFSP,NFS
COMMON /FACTOR/ FSP,NFS
NEWFSP = FSP + 2
IF (NEWFSP .GT. LFS) CALL DERRM('OVERFLOW OF BUFFER@')
IFS(FSP+1) = IND1
RFS(FSP+1) = DER1
IFS(NEWFSP) = FSP
RFS(NEWFSP) = DEPVAR
FSP = NEWFSP
RETURN
C
C LAST CARD OF SUBROUTINE DMIT1.
C
END
C******************************************************************
SUBROUTINE DMIT2(IND1,DER1,IND2,DER2,DEPVAR,RFS,IFS,LFS)
INTEGER IND1,IND2,DEPVAR,LFS
INTEGER IFS(LFS)
DOUBLE PRECISION DER1,DER2
DOUBLE PRECISION RFS(LFS)
INTEGER FSP,NEWFSP,NFS
COMMON /FACTOR/ FSP,NFS
NEWFSP = FSP + 3
IF (NEWFSP .GT. LFS) CALL DERRM('OVERFLOW OF BUFFER@')
IFS(FSP+1) = IND1
RFS(FSP+1) = DER1
IFS(FSP+2) = IND2
RFS(FSP+2) = DER2
IFS(NEWFSP) = FSP
RFS(NEWFSP) = DEPVAR
FSP = NEWFSP
RETURN
C
C LAST CARD OF SUBROUTINE DMIT2.
C
END
C******************************************************************
SUBROUTINE DPCOPY(GRAD,IGRAD,STRIDE,BUFFER,N)
INTEGER IGRAD,STRIDE,N
DOUBLE PRECISION GRAD(STRIDE,N),BUFFER(N)
INTEGER I
DO 10 I = 1, N
GRAD(IGRAD,I) = BUFFER(I)
10 CONTINUE
IGRAD = IGRAD + N*N
RETURN
C
C LAST CARD OF SUBROUTINE DPCOPY.
C
END
C******************************************************************
SUBROUTINE DPGRAD(YGRAD,LYGRAD,ID,RGRAD,IGRAD,RFS,IFS,LFS)
INTEGER LYGRAD,ID,RGRAD,IGRAD,LFS
INTEGER IFS(LFS)
DOUBLE PRECISION YGRAD(LYGRAD),RFS(LFS)
INTEGER DEPVAR,FSP,FSPEND,I,IND,NFS,OLDFSP
DOUBLE PRECISION ONE,T,ZERO
COMMON /FACTOR/ FSP,NFS
DATA ZERO,ONE /0.0D0,1.0D0/
DO 10 I = 1, LYGRAD
YGRAD(I) = ZERO
10 CONTINUE
YGRAD(ID) = ONE
RGRAD = RGRAD + 1
IGRAD = RGRAD
FSPEND = FSP
20 CONTINUE
IF (FSPEND .EQ. 0) THEN
IF (NFS .EQ. 0) GO TO 40
CALL DERRM('NOTHING TO BE READ. THIS IS AN ERROR@')
END IF
DEPVAR = RFS(FSPEND)
OLDFSP = FSPEND - 1
FSPEND = IFS(FSPEND)
T = YGRAD(DEPVAR)
IF (T .EQ. ZERO) GO TO 20
YGRAD(DEPVAR) = ZERO
30 CONTINUE
IF (OLDFSP .LE. FSPEND) GO TO 20
IND = IFS(OLDFSP)
YGRAD(IND) = YGRAD(IND) + T*RFS(OLDFSP)
OLDFSP = OLDFSP - 1
GO TO 30
40 CONTINUE
RETURN
C
C LAST CARD OF SUBROUTINE DPGRAD.
C
END
C******************************************************************
SUBROUTINE DPINIT(NEED,HAVE)
INTEGER NEED,HAVE
INTEGER FSP,NFS
COMMON /FACTOR/ FSP,NFS
FSP = 0
NFS = 0
IF (NEED .GT. HAVE)
* CALL DERRM('VECTOR YGRAD OR YJACOB IS TOO SMALL@')
RETURN
C
C LAST CARD OF SUBROUTINE DPINIT.
C
END
C******************************************************************
SUBROUTINE DERRM(STR)
C
C THIS ROUTINE WRITES AN ERROR MESSAGE AND ABORTS.
C
CHARACTER*1 STR(72)
INTEGER I,L,NWRITE
DATA NWRITE /6/
WRITE (NWRITE,'FATAL ERROR IN DPJAKELIB')
DO 10 L = 1, 72
IF (STR(L) .EQ. '@') GO TO 20
10 CONTINUE
20 CONTINUE
L = L - 1
WRITE (NWRITE,'(72A1)') (STR(I),I = 1, L)
STOP
C
C LAST CARD OF SUBROUTINE DERRM.
C
END