/*Translated by FOR_C, v3.4.2 (-), on 07/09/115 at 08:30:05 */
/*FOR_C Options SET: ftn=u io=c no=p op=aimnv s=dbov str=l x=f - prototypes */
#include
#include "fcrt.h"
#include "djacg.h"
#include
#include
/* PARAMETER translations */
#define LNTBUF 11
#define LRBUF 18
/* end of PARAMETER translations */
void /*FUNCTION*/ djacg(
long *mode,
long m,
long n,
double y[],
double f[],
double *fjac,
long ldfjac,
double yscale[],
double fac[],
long iopt[],
double wk[],
long lwk,
long iwk[],
long liwk)
{
#define FJAC(I_,J_) (*(fjac+(I_)*(ldfjac)+(J_)))
long int _l0, i, intbuf[LNTBUF], ircmp, irdel, irow, irowmx, itry,
j, jcol, nifac, ns, statea, stated, states;
double adiff, ay, cfact, del, delm, dfmj, diff, dmax,
facmax, facmin, fjacl, fmj, pert, rbuf[LRBUF], rdel, rmnfdf,
rmxfdf, savdel, savy, sdf, sf, sgn, t1, t2, u, u3qrt, u7egt,
uegt, umegt, uqrt, usqt;
/* OFFSET Vectors w/subscript range: 1 to dimension */
double *const F = &f[0] - 1;
double *const Fac = &fac[0] - 1;
long *const Intbuf = &intbuf[0] - 1;
long *const Iopt = &iopt[0] - 1;
long *const Iwk = &iwk[0] - 1;
double *const Rbuf = &rbuf[0] - 1;
double *const Wk = &wk[0] - 1;
double *const Y = &y[0] - 1;
double *const Yscale = &yscale[0] - 1;
/* end of OFFSET VECTORS */
/* Copyright (c) 2006, Math a la Carte, Inc.
*>> 2006-04-11 DJACG R. J. Hanson removed unneeded values saved.
*>> 2003-07-08 DJACG R. J. Hanson fixed bug checking FAC(*) < 0.
*>> 2003-07-08 DJACG R. J. Hanson fixed bug on MODE < 0 or > N.
*>> 2002-06-21 DJACG R. J. Hanson modified Salane Code
*--D replaces "?": ?JACG */
/* Main routine DJACG, computing numerical derivatives.
* -----------------------------------------------------------------
*
****BEGIN PROLOGUE DJACG
****DATE WRITTEN 860410
****Date Modified 060411
****CATEGORY NO. D4
****KEYWORDS NUMERICAL DIFFERENCING, Numerical Jacobians, Numerical
* Derivatives
****AUTHOR SALANE, DOUGLAS E., SANDIA NATIONAL LABORATORIES
* NUMERICAL MATHEMATICS DIVISION,1642
* ALBUQUERQUE, NM 87185
* Interface reworked by Richard J. Hanson, Rice University.
* Last changes made by Hanson after leaving Rice University.
*
****PURPOSE SUBROUTINE DJACG
*
* Subroutine DJACG uses finite differences to compute the Jacobian for
* a system of (M) equations in (N) variables. DJACG is designed
* for use in numerical methods for solving nonlinear problems where a
* Jacobian is evaluated repeatedly at neighboring arguments.
* For example in a Gauss-Newton method for solving non-linear least squares
* problems. DJACG is suitable for applications in which the Jacobian
* is a dense rectangular matrix, over- exactly or under-determined.
* The interface to the values of equations is with reverse
* communication. The code DJACG saves and restores values between
* calls. In particular DJACG is thread-safe. */
/* The design allows for computation of derivatives that are a sum of
* analytic and functional values that are not easily differentiated.
****DESCRIPTION
*
* SUBROUTINE PARAMETERS
*
* MODE.......An integer flag that directs user action with reverse
* communication. Enter the routine with MODE=0.
* DJACG returns with MODE=J when computing column J of
* the Jacobian. When MODE=0 on return, the Jacobian
* has been computed.
* N..........[in] THE NUMBER OF VARIABLES.
* M..........[in] THE NUMBER OF EQUATIONS.
* Y(*).......[inout] AN ARRAY OF DIMENSION N. THE POINT AT WHICH THE
* Jacobian IS TO BE EVALUATED.
* F(*).......[in] AN ARRAY OF DIMENSION M. THE EQUATIONS EVALUATED AT
* THE POINT Y. This must be defined when MODE .eq. 0 or
* just after the first perturbed column is computed,
* MODE .eq. 1.
* FJAC(*,*)..[inout] AN ARRAY OF SIZE LDFJAC BY N. FJAC
* CONTAINS THE Jacobian WHICH IS AN M BY N MATRIX.
* Columns that are accumulated must have the analytic
* part defined on entry or else be set to zero.
* Columns that are skipped can be defined either before
* or after the routine exits.
* LDFJAC.....[in] THE LEADING DIMENSION OF THE ARRAY FJAC AS DIMENSIONED
* IN THE CALLING ROUTINE. LDFJAC MUST BE GREATER OR
* EQUAL TO M.
* YSCALE(*)..[in] AN ARRAY OF LENGTH N. THE USER CAN PROVIDE
* REPRESENTATIVE MAGNITUDES FOR Y VALUES IN THE ARRAY
* YSCALE. THE USER CAN ALSO USE YSCALE TO PROVIDE
* APPROPRIATE SIGNS FOR THE INCREMENTS. */
/* FAC(*).....[inout] AN ARRAY OF LENGTH N. FAC(I) CONTAINS THE PERCENTAGE
* FACTOR FOR DIFFERENCING.
*
* If FAC(1)=0., the array FAC is set to the square root of
* machine unit roundoff on the first call to DJACG.
* Unless the user wishes to initialize FAC,
* the FAC array should not be altered between subsequent calls
* to DJACG.
*
* The user may provide FAC values if desired.
* If the user provides FAC values, FAC(i) should be
* set to a value between 0 and 1. DJACG will not permit FAC(i) to
* be set to a value that results in too small an increment. DJACG
* ensures that
* FACMIN <= FAC(i) <= FACMAX.
* For further details on how the code chooses FACMIN and FACMAX see
* the report [2].
*
*
* IOPT(*) [in] An integer array defining the methods used to
* compute the derivatives. The default is to use one sided differences.
* Entries in this array are interpreted as follows. */
/* [0] Use the current settings for all remaining variables. The starting
* setting is to simply use one sided differences.
* [k > 0] Use the current settings for all variables from the last
* specified up to and including variable $k$.
* [-1] Set to use one sided differences. (Not needed at the
* beginning since this is the default state.)
* [-2] Set to use central differences.
* [-3] Set to accumulate the result from whatever type of differences
* have been specified above into initial values provided in FJAC. This */
/* must be followed by a number $\geq 0$ as any number $< 0$ will turn
* this off.
* [-4] Skip variables. This must be followed by a value $\geq 0$
* as any number $< 0$ will turn this off. */
/* Re. state -2: For variable p, compute numerical derivatives using
* central divided differences. This will typically yield more
* accuracy than the one-sided differences, but with the expense
* of an additional function evaluation per variable.
* The increment used for central differencing
* is the default T=$macheps^{-1/6}$
* times the increment used in one-sided differencing.
* To change this factor for succeeding
* variables, assign a new value between calls with MODE .gt. 0
* in the storage location FAC(MODE). */
/* The default value T=$macheps^{-1/6}$ is based on the
* approximate relation T*FAC(J) =$macheps^{2/3}$.
* This value is near an optimal choice, under certain conditions
* on higher derivatives, provided FAC(J) = $macheps^{1/2}$.
* Sometimes larger or smaller values of T will give more accuracy
* than the default. */
/* For example if column 3 is not to be computed and column 4 is
* to be accumulated, then IOPT(1)=2, followed by:
* IOPT(2)=-4, IOPT(3)=3, IOPT(4)=-3, IOPT(5)=4,
* IOPT(6)=-1, IOPT(7)=0 */
/* wk(*)...[inout]a work array whose dimension is at least LWK.
* LWK.....[in] the required length of the work array. LWK=3*M+18.
* iwk(*)..[inout] an integer array whose dimension is at least LIWK.
* The first 10 positions of IWK contain diagnostics
* information.
*
* IWK(1) gives the number of times a function evaluation was computed.
*
* IWK(2) gives the number of columns in which three attempts were
* made to increase a percentage factor for differencing (i.e. a
* component in the FAC array) but the computed DEL remained
* unacceptably small relative to Y(JCOL) or YSCALE(JCOL). In such
* cases the percentage factor is set to the square root of the unit
* roundoff of the machine.
*
* IWK(3) gives the number of columns in which the computed DEL was
* zero to machine precision because Y(JCOL) or YSCALE(JCOL) was
* zero. In such cases DEL is set to the square root of the unit
* roundoff.
*
* IWK(4) gives the number of Jacobian columns which had to be
* recomputed because the largest difference formed in the column
* was close to zero relative to scale, where
*
* scale = max(|f (y)|,|f (y+del)|)
* i i
*
* and i denotes the row index of the largest difference in the
* column currently being processed. IWK(10) gives the last column
* where this occurred.
*
* IWK(5) gives the number of columns whose largest difference is
* close to zero relative to scale after the column has been
* recomputed.
*
* IWK(6) gives the number of times scale information was not
* available for use in the roundoff and truncation error tests.
* this occurs when
*
* min (|f (y)|,|f (y+del)|) = 0.
* i i
*
* where i is the index of the largest difference for the column
* currently being processed.
*
* IWK(7) gives the number of times the increment for differencing
* (DEL) was computed and had to be increased because (Y(JCOL)+DEL )
* -Y(JCOL)) was too small relative to Y(JCOL) or YSCALE(JCOL).
*
* IWK(8) gives the number of times a component of the FAC array was
* reduced because changes in function values were large and excess
* truncation error was suspected. IWK(9) gives the last column in
* which this occurred.
*
* IWK(9) gives the index of the last column where the corresponding
* component of the FAC array had to be reduced because excessive
* truncation error was suspected.
*
* IWK(10) gives the index of the last column where the differeence
* was small and the column had to be recomputed with an adjusted
* increment (see iwk(4)). the largest derivative in this column
* may be inaccurate due to excessive roundoff error.
*
* LIWK.[in] the length of the array iwk. LIWK >= 10+LNTBUF=21. */
/* Using reverse communication, DJACG will return control to
* the user each time a function value is required. The user must
* provide the values of the function f evaluated at y where y is the */
/* array return by DJACG. The user must set
*
* wk(i) = i-th component of f evaluated at y
*
* for i = 1,2,...,M. After wk is assigned, the user must recall DJACG.
* When DJACG returns to the user with
*
* MODE = 0,
*
* the computation of the Jacobian is complete.
*
* Roundoff and truncation errors.
*
* Subroutine DJACG takes advantage of the way in which the Jacobian
* is evaluated to adjust increments for differencing to control
* roundoff and truncation errors. The routine usually (but not always)
* requires one additional function evaluation to compute a column of the
* Jacobian. Also, the routine returns a variety of error diagnostics
* to warn users when computed derivatives may not be accurate.
*
* WARNING: DJACG does not guarantee the accuracy of the computed
* derivatives. In order to save on function evaluations, heuristic
* tecniques for increment adjustment and safeguarding increments are
* used. These usually work well, but are not guaranteed.
*
* In order to determine the accuracy of computed derivates, users should
* pay attention to the diagnostic array positions IWK(4), IWK(8), IWK(9)
* and IWK(10). Non-zeros in these position indicate
* possible large roundoff or truncation errors. Here is a summary.
*
* IWK(4) NONZERO => ROUNDOFF ERROR. IWK(10) LAST NOTED COLUMN.
* IWK(8) NONZERO => TRUNCATION ERROR. IWK(9) LAST NOTED COLUMN.
*
* WARNING: Some of the diagnostics returned can only be interpreted
* with a detailed knowledge of the routine. Nevertheless, they are
* provided to give users full access to the information produced by
* the subroutine.
*
****references (1)D.E. Salane and L. F. Shampine
* "An economical and efficient routine for computing
* sparse Jacobians", Report no. SAND85-0977,
* Sandia National Laboratories, Albuquerque,NM,87185.
* (2)D.E. Salane
* "Adaptive routines for forming Jacobians
* numerically", Report no. SAND86-1319, Sandia
* National Laboratories, Albuquerque, NM, 87185
*
****ROUTINES CALLED D1MACH
*
* REQUIRED FORTRAN INTRINSIC FUNCTIONS:
* ABS,MAX,MIN,SIGN,SQRT
* REQUIRED Math a la Carte FUNCTIONS:
* D1MACH
*
****END PROLOGUE DJACG
* IMPLICIT NONE */
/* SET CONSTANTS AND ALGORITHM PARAMETERS.
* PERT........IS USED IN EXTRAPOLATION TEST.
* FACMAX......FAC(I) <= FACMAX.
* EXPFMN......FAC(I) => FACMIN WHERE FACMIN = U ** EXPFMN.
* NIFAC.......ITERATION BOUND USED IN COMPUTING INCREMENT. */
/****FIRST EXECUTABLE STATEMENT DJACG */
if (*mode > 0 && *mode <= n)
goto L_90;
/* Do error checking on some parameters. A returned value of MODE < 0
* denotes that argument number -MODE has an error condition associated
* with it.
* SUBROUTINE DJACG(MODE,M,N,Y,F,
* . FJAC,LDFJAC,YSCALE,FAC,IOPT,
* . WK,LWK,IWK,LIWK)
* Variables 1-5, 6-10, 11-14
* Must have MODE=0 the first time. */
if (*mode < 0 || *mode > n)
{
/* This says that MODE < 0 or MODE > N is error. */
*mode = 1;
goto L_1000;
}
*mode = 2;
/* Must have M > 0 */
if (m <= 0)
goto L_1000;
*mode = 3;
/* Must have N > 0 */
if (n <= 0)
goto L_1000;
*mode = 7;
/* Must have LDFJAC .ge. M if N > 1. */
if (ldfjac < m && n > 1)
goto L_1000;
/* Must have FAC() >= 0. The value FAC(1)=0.
* skips this test on the rest of FAC(*). */
if (Fac[1] != 0.)
{
*mode = 9;
for (j = 1; j <= n; j++)
{
if (Fac[j] < 0.)
goto L_1000;
}
}
*mode = 10;
/* Run throught the IOPT() array making sure that the
* state changes and indices make sense. */
ns = 1;
L_1070:
;
/* A positive value .gt. N is an error. */
if (Iopt[ns] > n)
goto L_1000;
/* A value of 0 (or N) says to use the present setting for
* the rest of the variables */
if (Iopt[ns] == 0 || Iopt[ns] == n)
goto L_1080;
/* A postitive value < N must be followed by a zero, N, or a
* state change. */
if (Iopt[ns] > 0)
{
ns += 1;
goto L_1070;
}
if (Iopt[ns] < -4)
goto L_1000;
if (Iopt[ns] >= -2)
{
ns += 1;
goto L_1070;
}
/* An accumulation or skip must be followed by a non-negative index. */
if (Iopt[ns + 1] < 0)
goto L_1000;
ns += 1;
goto L_1070;
L_1080:
;
*mode = 12;
/* Must have LWK >= 3*M+LRBUF */
if (lwk < 3*m + LRBUF)
goto L_1000;
*mode = 14;
/* Must have LIWK >= 21 */
if (liwk < 21)
goto L_1000;
/* Since reverse communication, initialize JCOL only on first call. */
jcol = 0;
/* COMPUTE APPROPRIATE MACHINE CONSTANTS.
* */
u = DBL_EPSILON;
usqt = sqrt( u );
/* UQRT = U**P25
* U3QRT = U**P75
* UEGT = U**P125
* U7EGT = U**P875
* FACMIN = U**EXPFMN */
uqrt = sqrt( usqt );
t1 = 1.e0/usqt;
u3qrt = u/uqrt;
uegt = sqrt( uqrt );
umegt = 1.e0/uegt;
u7egt = u*umegt;
facmin = u3qrt;
pert = 2.0e0;
facmax = 1.e-1;
nifac = 3;
for (j = 1; j <= liwk; j++)
{
Iwk[j] = 0;
}
/* The value FAC=0. means that this code wants to initialize. */
if (Fac[1] == 0.e0)
{
for (j = 1; j <= n; j++)
{
Fac[j] = usqt;
}
}
/* Default to not accumulate: */
statea = 1;
/* Default to use one-sided differences. */
stated = 1;
cfact = 0.e0;
/* Default to not skip computing columns. */
states = 1;
/* Index for processing IOPT(*) states and changes. */
ns = 1;
L_30:
;
jcol += 1;
/* IF IOPT(NS) .eq. 0, the current states hold for the rest
* of the columns. */
L_1020:
;
if (Iopt[ns] != 0)
{
/* A state change will occur: */
if (jcol == Iopt[ns] + 1)
{
/* This signals a state change. Namely the last variable
* in the group has been processed and a change is required. */
ns += 1;
}
if (Iopt[ns] < 0)
{
/* Set for one-sided or central differences. */
if (Iopt[ns] >= -2)
{
stated = -Iopt[ns];
statea = 1;
}
else
{
/* Set for accumulation or skipping variables. */
if (Iopt[ns] == -3)
statea = 2;
if (Iopt[ns] == -4)
{
states = 2;
statea = 1;
}
}
ns += 1;
}
else
{
goto L_1030;
}
/* Still in range for a state, without change to that state. */
goto L_1020;
}
L_1030:
;
if (stated > 1)
{
if (cfact == 0e0)
{
/* CFACT=U**(-1.D0/6.D0)
* Compute U**(-1/6) using Newton's method; compute the cube
* root of U**(-1/2). This is to avoid the ** elmentary function
* which may not be thread-safe or re-entrant. */
cfact = umegt*umegt;
for (j = 1; j <= 14; j++)
{
cfact = (2.e0*cfact + t1/(cfact*cfact))/3.e0;
/* Converged, to over 7/8 precision. Escapes from loop. */
if (fabs( cfact*cfact*cfact*usqt - 1.e0 ) < u7egt)
goto L_1060;
}
L_1060:
;
}
}
/* Check if this column is to be skipped. */
if (states > 1)
goto L_170;
ircmp = 0;
itry = 0;
/* Check if this column is to be accumulated. */
if (statea > 1)
{
/* Save the initial contents of the array column for final accumulation. */
for (irow = 1; irow <= m; irow++)
{
Wk[irow + 2*m] = FJAC(jcol - 1,irow - 1);
}
}
L_50:
;
savy = Y[jcol];
/* Compute DEL. If DEL is too small increase FAC and recompute
* DEL. NIFAC attempts are made to increase FAC and find an
* appropriate DEL. If DEL can't be found in this manner, DEL is computed
* with FAC set to the square root of the machine precision(=USQT).
* If DEL is zero to machine precision because T is zero or
* YSCALE is zero, DEL is set to USQT.
* */
irdel = 0;
if (Yscale[1] != 0.e0)
{
ay = fabs( Yscale[jcol] );
}
else
{
ay = fabs( Y[jcol] );
}
sgn = 1.e0;
if (Yscale[1] != 0.e0)
sgn = sign( 1.e0, Yscale[jcol] );
delm = u7egt*ay;
for (j = 1; j <= nifac; j++)
{
del = Fac[jcol]*ay*sgn;
if (del == 0.e0)
{
del = usqt*sgn;
if (itry == 0)
Iwk[3] += 1;
}
t1 = Y[jcol] + del;
del = t1 - Y[jcol];
if (fabs( del ) < delm)
{
if (j >= nifac)
goto L_60;
if (irdel == 0)
{
irdel = 1;
Iwk[7] += 1;
}
t1 = Fac[jcol]*umegt;
Fac[jcol] = fmin( t1, facmax );
}
else
{
goto L_70;
}
L_60:
;
}
Fac[jcol] = usqt;
del = usqt*ay*sgn;
Iwk[2] += 1;
L_70:
;
/* See if central divided differences are intended. */
if (stated > 1)
{
del *= cfact;
}
savdel = del;
Y[jcol] += del;
Iwk[1] += 1;
/* With reverse communication, DJACG returns to calling routine
* for function values. Afterwards, DJACG transfers control
* to the statement after the comment - ENTRY POINT FOR REVERSE
* COMMUNICATION. */
/* Pack up scalars that need saving. Place in store within
* WK(*) and IWK(*). */
goto L_190;
L_80:
;
*mode = jcol;
return;
/* ENTRY POINT FOR REVERSE COMMUNICATON. */
L_90:
;
/* Unpack scalars that needed saving.
* Take from WK(*) and IWK(*) where they are stored. */
goto L_220;
L_100:
;
/* See if this is the second evaluation for a central difference.
* This is true STATED .gt. 1 .and. ITRY .eq. 1.
* and ITRY = 1. */
rdel = 1.e0/savdel;
if (stated == 1)
goto L_140;
/* IF ITRY = 1, proceed and compute the central divided difference. */
if (itry == 1)
goto L_120;
/* Otherwise ITRY = 0, so an alternate evaluation is needed. */
Y[jcol] = savy - del;
for (i = 1; i <= m; i++)
{
Wk[i + m] = Wk[i];
}
itry = 1;
Iwk[1] += 1;
goto L_190;
L_120:
;
/* Compute the central divided difference. */
rdel *= .5;
for (i = 1; i <= m; i++)
{
FJAC(jcol - 1,i - 1) = rdel*(Wk[i + m] - Wk[i]);
}
/* Restore the component of Y() and move on. */
Y[jcol] = savy;
itry = 0;
goto L_170;
/* Restore the component of Y() and move on. */
L_140:
;
Y[jcol] = savy;
/* Compute the Jacobian entries.
* Use largest elements in a column to determine scaling
* information for roundoff and truncation error tests.
* */
dmax = 0.e0;
irowmx = 1;
for (irow = 1; irow <= m; irow++)
{
diff = Wk[irow] - F[irow];
adiff = fabs( diff );
if (adiff >= dmax)
{
irowmx = irow;
dmax = adiff;
sf = F[irow];
sdf = Wk[irow];
}
fjacl = diff*rdel;
if (itry == 1)
{
Wk[irow] = fjacl;
Wk[irow + m] = FJAC(jcol - 1,irow - 1);
}
FJAC(jcol - 1,irow - 1) = fjacl;
}
/* If a column is being recomputed (ITRY=1), this section of the
* code performs an extrapolation test to enable the code to
* compute small derivatives more accurately.
* */
if (itry == 1)
{
t1 = Wk[irowmx + m];
t2 = Wk[irowmx]*Fac[jcol];
if (fabs( t2 ) < fabs( t1 )*pert)
{
t1 = Fac[jcol]*Fac[jcol];
Fac[jcol] = fmax( t1, facmin );
for (irow = 1; irow <= m; irow++)
{
fjacl = Wk[irow + m];
FJAC(jcol - 1,irow - 1) = fjacl;
}
}
}
fmj = fabs( sf );
dfmj = fabs( sdf );
rmxfdf = fmax( fmj, dfmj );
rmnfdf = fmin( fmj, dfmj );
/* If scale information is not available, perform no roundoff
* or truncation error tests.
* */
if (rmnfdf != 0.e0)
{
/* Test for possible roundoff error (first test)
* and also for possible serious roundoff error (second test).
* */
if (dmax <= (u3qrt*rmxfdf))
{
if (dmax <= (u7egt*rmxfdf))
{
if (itry == 0)
{
t1 = sqrt( Fac[jcol] );
Fac[jcol] = fmin( t1, facmax );
ircmp = 1;
Iwk[4] += 1;
Iwk[10] = jcol;
}
else
{
Iwk[5] += 1;
}
}
else
{
t1 = umegt*Fac[jcol];
Fac[jcol] = fmin( t1, facmax );
}
}
/* Test for possible truncation error.
* */
if (dmax > uqrt*rmxfdf)
{
t1 = Fac[jcol]*uegt;
Fac[jcol] = fmax( t1, facmin );
Iwk[8] += 1;
Iwk[9] = jcol;
}
}
else
{
Iwk[6] += 1;
}
/* If serious roundoff error is suspected, recompute the
* column.
*
* */
if (ircmp == 1)
{
ircmp = 0;
itry = 1;
goto L_50;
}
itry = 0;
/* Branch to this place if this column is skipped. */
L_170:
;
/* May have asked for this column to be accumulated. */
if (statea > 1 && states == 1)
{
for (irow = 1; irow <= m; irow++)
{
FJAC(jcol - 1,irow - 1) += Wk[irow + 2*m];
}
}
if (jcol < n)
goto L_30;
/* This is the signal that all Jacobian columns are computed. */
*mode = 0;
return;
/* Internal procedure to save local data.
* This step is needed to make the code re-entrant and threadsafe. */
L_190:
;
Rbuf[1] = facmin;
Rbuf[2] = facmax;
Rbuf[3] = cfact;
Rbuf[4] = usqt;
Rbuf[5] = uegt;
Rbuf[6] = umegt;
Rbuf[7] = uqrt;
Rbuf[8] = u3qrt;
Rbuf[9] = u7egt;
Rbuf[10] = u;
Rbuf[11] = cfact;
Rbuf[12] = savdel;
Rbuf[13] = pert;
Rbuf[14] = savy;
Rbuf[15] = ay;
Rbuf[16] = sgn;
Rbuf[17] = del;
Rbuf[18] = delm;
/* Save data values in calling program unit. When the code is
* entered again these are extracted from storage the store
* in the calling program. This eliminates the need to save
* any variables. */
for (i = 1; i <= LRBUF; i++)
{
Wk[i + 3*m] = Rbuf[i];
}
Intbuf[1] = nifac;
Intbuf[2] = i;
Intbuf[3] = j;
Intbuf[4] = itry;
Intbuf[5] = irdel;
Intbuf[6] = ircmp;
Intbuf[7] = jcol;
Intbuf[8] = ns;
Intbuf[9] = statea;
Intbuf[10] = stated;
Intbuf[11] = states;
/* Save data values in calling program unit. When the code is
* called again these are extracted. */
for (i = 1; i <= LNTBUF; i++)
{
Iwk[i + 10] = Intbuf[i];
}
goto L_80;
/* Internal procedure to unpack saved data. */
L_220:
;
/* Extract saved data from store in the calling program unit. */
for (i = 1; i <= LRBUF; i++)
{
Rbuf[i] = Wk[i + 3*m];
}
facmin = Rbuf[1];
facmax = Rbuf[2];
cfact = Rbuf[3];
usqt = Rbuf[4];
uegt = Rbuf[5];
umegt = Rbuf[6];
uqrt = Rbuf[7];
u3qrt = Rbuf[8];
u7egt = Rbuf[9];
u = Rbuf[10];
cfact = Rbuf[11];
savdel = Rbuf[12];
pert = Rbuf[13];
savy = Rbuf[14];
ay = Rbuf[15];
sgn = Rbuf[16];
del = Rbuf[17];
delm = Rbuf[18];
/* Extract saved data from store in the calling program unit. */
for (i = 1; i <= LNTBUF; i++)
{
Intbuf[i] = Iwk[i + 10];
}
nifac = Intbuf[1];
i = Intbuf[2];
j = Intbuf[3];
itry = Intbuf[4];
irdel = Intbuf[5];
ircmp = Intbuf[6];
jcol = Intbuf[7];
ns = Intbuf[8];
statea = Intbuf[9];
stated = Intbuf[10];
states = Intbuf[11];
goto L_100;
L_1000:
;
/* Error return when arguments are invalid. */
*mode = -*mode;
return;
#undef FJAC
} /* end of function */