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