/*Translated by FOR_C, v3.4.2 (-), on 07/09/115 at 08:30:10 */ /*FOR_C Options SET: ftn=u io=c no=p op=aimnv s=dbov str=l x=f - prototypes */ #include #include "fcrt.h" #include "sjacg.h" #include #include /* PARAMETER translations */ #define LNTBUF 11 #define LRBUF 18 /* end of PARAMETER translations */ void /*FUNCTION*/ sjacg( long *mode, long m, long n, float y[], float f[], float *fjac, long ldfjac, float yscale[], float fac[], long iopt[], float 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; float 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 */ float *const F = &f[0] - 1; float *const Fac = &fac[0] - 1; long *const Intbuf = &intbuf[0] - 1; long *const Iopt = &iopt[0] - 1; long *const Iwk = &iwk[0] - 1; float *const Rbuf = &rbuf[0] - 1; float *const Wk = &wk[0] - 1; float *const Y = &y[0] - 1; float *const Yscale = &yscale[0] - 1; /* end of OFFSET VECTORS */ /* Copyright (c) 2006, Math a la Carte, Inc. *>> 2006-04-11 SJACG R. J. Hanson removed unneeded values saved. *>> 2003-07-08 SJACG R. J. Hanson fixed bug checking FAC(*) < 0. *>> 2003-07-08 SJACG R. J. Hanson fixed bug on MODE < 0 or > N. *>> 2002-06-21 SJACG R. J. Hanson modified Salane Code *--S replaces "?": ?JACG */ /* Main routine SJACG, computing numerical derivatives. * ----------------------------------------------------------------- * ****BEGIN PROLOGUE SJACG ****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 SJACG * * Subroutine SJACG uses finite differences to compute the Jacobian for * a system of (M) equations in (N) variables. SJACG 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 squ * problems. SJACG 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 SJACG saves and restores values between * calls. In particular SJACG 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. * SJACG 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 DIMENSIONE * 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 PERCENTA * 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 SJACG. * Unless the user wishes to initialize FAC, * the FAC array should not be altered between subsequent calls * to SJACG. * * 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. SJACG will not permit FAC(i) to * be set to a value that results in too small an increment. DJAC * 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 difference * Entries in this array are interpreted as follows. */ /* [0] Use the current settings for all remaining variables. The startin * 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. T */ /* must be followed by a number $\geq 0$ as any number $< 0$ will tur * 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 * 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 condition * on higher derivatives, provided FAC(J) = $macheps^{1/2}$. * Sometimes larger or smaller values of T will give more accura * 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 comput * * 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 colu * 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, SJACG 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 SJACG. 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 SJACG. * When SJACG returns to the user with * * MODE = 0, * * the computation of the Jacobian is complete. * * Roundoff and truncation errors. * * Subroutine SJACG 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 * Jacobian. Also, the routine returns a variety of error diagnostics * to warn users when computed derivatives may not be accurate. * * WARNING: SJACG 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 R1MACH * * REQUIRED FORTRAN INTRINSIC FUNCTIONS: * ABS,MAX,MIN,SIGN,SQRT * REQUIRED Math a la Carte FUNCTIONS: * R1MACH * ****END PROLOGUE SJACG * 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 SJACG */ 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 SJACG(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 = FLT_EPSILON; usqt = sqrtf( u ); /* UQRT = U**P25 * U3QRT = U**P75 * UEGT = U**P125 * U7EGT = U**P875 * FACMIN = U**EXPFMN */ uqrt = sqrtf( usqt ); t1 = 1.e0/usqt; u3qrt = u/uqrt; uegt = sqrtf( 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.E0/6.E0) * 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 (fabsf( 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 compute * 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 = fabsf( Yscale[jcol] ); } else { ay = fabsf( Y[jcol] ); } sgn = 1.e0; if (Yscale[1] != 0.e0) sgn = signf( 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 (fabsf( del ) < delm) { if (j >= nifac) goto L_60; if (irdel == 0) { irdel = 1; Iwk[7] += 1; } t1 = Fac[jcol]*umegt; Fac[jcol] = fminf( 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, SJACG returns to calling routine * for function values. Afterwards, SJACG 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 = fabsf( 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 (fabsf( t2 ) < fabsf( t1 )*pert) { t1 = Fac[jcol]*Fac[jcol]; Fac[jcol] = fmaxf( t1, facmin ); for (irow = 1; irow <= m; irow++) { fjacl = Wk[irow + m]; FJAC(jcol - 1,irow - 1) = fjacl; } } } fmj = fabsf( sf ); dfmj = fabsf( sdf ); rmxfdf = fmaxf( fmj, dfmj ); rmnfdf = fminf( 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 = sqrtf( Fac[jcol] ); Fac[jcol] = fminf( t1, facmax ); ircmp = 1; Iwk[4] += 1; Iwk[10] = jcol; } else { Iwk[5] += 1; } } else { t1 = umegt*Fac[jcol]; Fac[jcol] = fminf( t1, facmax ); } } /* Test for possible truncation error. * */ if (dmax > uqrt*rmxfdf) { t1 = Fac[jcol]*uegt; Fac[jcol] = fmaxf( 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 */