/*Translated by FOR_C, v3.4.2 (-), on 07/09/115 at 08:31:18 */ /*FOR_C Options SET: ftn=u io=c no=p op=aimnv s=dbov str=l x=f - prototypes */ #include #include "fcrt.h" #include "dxrk8.h" #include #include #include /* PARAMETER translations */ #define HCAU 25 #define HMAX 11 #define INTSAV 41 #define KACT12 (10*LREVF + 1) #define KACT13 (10*LREVF + 1) #define KACT14 (10*LREVO + 1) #define KACT15 (10*(LHMAXU + 100) + 2) #define KACT16 (10*(LHMIN + 100) + 2) #define KACT17 (10*(LEPS + 100) + 2) #define KACT20 (10*KSTEPM + 2) #define KACT27 (10*(LERTRY + 100) + 3) #define KACT28 (10*(LINTX + 100) + 2) #define KFAIL 28 #define KRKAEV 5 #define KRKEOS 21 #define KRKFR2 13 #define KRKIFS 28 #define KRKNOE 7 #define KRKRE1 3 #define KRKREV 6 #define KRKTS 23 #define KRKXG 25 #define KSTEP 27 #define KSTEPM 34 #define KSTEPX 36 #define LAST0 NGHI #define LBF 13 #define LEPS 1 #define LERTRY 14 #define LHMAXU 2 #define LHMIN 3 #define LINTX 9 #define LOCEOS 18 #define LOCEPR 11 #define LOCERR 44 #define LOCINT 38 #define LOCOUT 14 #define LOCTY 30 #define LOCXP 12 #define LREVF 23 #define LREVO 24 #define LTFIN 42 #define LTXTAB 8 #define LTXTAC 27 #define LTXTAD 90 #define LTXTAE 156 #define LTXTAF 209 #define LTXTAG 250 #define LTXTAH 319 #define LTXTAI 349 #define LTXTAJ 418 #define LTXTAK 448 #define MECONT 50 #define MEEMES 52 #define MEFVEC 61 #define MEIVEC 57 #define MERET 51 #define METEXT 53 #define NEGPOS 6 #define NGHI 32 #define NTGS 17 #define NUMINT 39 #define NXGS 16 /* end of PARAMETER translations */ void /*FUNCTION*/ dxrk8( double ts[], double y[], double opt[], long idat[], double dat[], double work[]) { LOGICAL32 gotint; long int _l0, i, ibeg, ierr, ii, itest, j, k, kat, kstate[6-(2)+1], kuse[6-(2)+1], kwork, l, ldat, leq, leql[6-(2)+1], lerr1, lidat, llerr, lo, lwork, maxler, minler, n, neq, neqnoe, next; double diffts, err1, err2; static char mtxtaa[3][182]={"DXRK8$BUndefined option.$EEquation indexes for groups are out of order or out of range.$EOption required inside equation group is outside the last group.$EOption specification made ", "twice in the same group..$EError tolerance vectors overlap in DAT.$EOptions must be turned on with 1.D0 and off with 0.D0, you used $F.$EIDAT(5) must be at least $I.$ECould not find", " space in DAT for all options. IDAT(6) must be > $I.$EIDAT(7) must be at least $I.$EThe paramters used in option 27, must have the first at least$ 4 times bigger than the second.$E"}; static char mtxtab[1][71]={"OPT up to the point where we are seeing problems.$NOPT: $BIDAT(4:7):$B"}; static long mact[5]={MEEMES,47,0,0,MECONT}; static long mactv[7]={METEXT,MEFVEC,0,METEXT,MEIVEC,4,MERET}; static long lermsg[10]={LTXTAB,LTXTAC,LTXTAD,LTXTAE,LTXTAF,LTXTAG, LTXTAH,LTXTAI,LTXTAJ,LTXTAK}; static long ierror[10]={14,11,21,13,17,26,18,19,20,27}; static long kact[28-(0)+1]={1,12,22,23,23,22,22,21,32,81,42,81, KACT12,KACT13,KACT14,KACT15,KACT16,KACT17,81,81,KACT20,71,62, 62,63,72,72,KACT27,KACT28}; /* OFFSET Vectors w/subscript range: 1 to dimension */ double *const Dat = &dat[0] - 1; long *const Idat = &idat[0] - 1; long *const Ierror = &ierror[0] - 1; long *const Lermsg = &lermsg[0] - 1; long *const Mact = &mact[0] - 1; long *const Mactv = &mactv[0] - 1; double *const Opt = &opt[0] - 1; double *const Ts = &ts[0] - 1; double *const Work = &work[0] - 1; double *const Y = &y[0] - 1; /* end of OFFSET VECTORS */ /* Copyright (c) 1996 California Institute of Technology, Pasadena, CA. * ALL RIGHTS RESERVED. * Based on Government Sponsored Research NAS7-03001. * Time-stamp: <2009-10-19 12:24:07 root> *>> 2009-10-15 DXRK8 Krogh Latest version of new stepsize control. *>> 2008-12-22 DXRK8 Krogh Many changes to put in new stepsize control. *>> 2008-02-28 DXRK8 Krogh Lots of changes in usage, bug fixes. *>> 2007-09-11 DXRK8 Krogh Fixed flag on max step output prior to final. *>> 2007-09-11 DXRK8 Krogh Fixed DAT(MAXINC) getting too small. *>> 2007-09-10 DXRK8 Krogh Some options were not processed. *>> 2007-09-10 DXRK8 Krogh Stored 0 in IDAT at end Err. Tol. Info. *>> 1997-12-18 DXRK8 Krogh Modifications for MATH77 library. *>> 1997-02-24 DXRK8 Krogh Initial code. *--D replaces "?": ?XRK8,?XRK8A,?XRK8F,?XRK8I,?XRK8N,?XRK8O,?XRK8X,?MESS * * Solves a system of first order ordinary differential equations * y'=f(t,y), using an eXplicit Runge-Kutta method of order 8 due to * Dormand & Prince. * This code was obtained starting with the code DOP853 by Hairer & * Wanner, Version of April 25, 1996. That code is described in * E. Hairer, S.P. Norsett and G. Wanner, Solving ordinary differential * equations I. Nonstiff Problems. 2nd Edition, Springer-Verlag (1993). * * This version written by Fred T. Krogh, has a different user interface, * avoids the function evaluations required for interpolation on steps * that do no interpolations, provides a G-Stop feature, uses a slightly * different algorithm for stepsize selection, and provides a number of * other options. * * ************* User Subprograms & calls to DXRK8I or DXRKIG *********** * * If reverse communication is not specified, for computing derivatives, * the user must code a subroutine of the form: * subroutine dxrk8f(TS, Y, F, IDAT) * And if reverse communication is not used for output, the user must * code a subroutine of the form: * subroutine dxrk8o(T, Y, IDAT, DAT) * If G-stops are used, then the user will be told to call dxrk8g from * time to time. Call looks different depending on whether they are made * from dxrk8f or from the users code as a result of using reverse * communication. * From dxrk8f: call dxrk8g(T, Y, F, IDAT) * Reverse Comm: call dxrk8g(TS, Y, WORK, IDAT) * After a return from dxrk8, if one is continuing an integration, one * should call dxrk8a instead of dxrk8 as follows * call dxrk8a(TS, Y, F, IDAT, DAT, WORK) * where F is only used if one has specified the option that uses it for * storing derivatives. * * ************************ Calling Sequence Arguments ****************** * * With the very simplest usage one sets: TS(1) = initial t;, TS(2) = 0 * (gives automatic selection of initial stepsize); TS(3) = final t; y = * initial values for y; IDAT(4) = NEQ = number of equations; IDAT(1) = 0 * (indicates starting); IDAT(5) to IDAT(7) to values giving the space * available in various arrays; and OPT(1) = 0.D0 to indicate no options * are begin specified. This gives an integration from TS(1) to TS(3) * with a default error tolerance, which gives an absolute and relative * error tolerance = (floating point epsilon)**.75. * * TS [in/out] Array holding information connected with the * independent variable. * TS(1) The current time, t. Initially the starting time. * TS(2) Current stepsize, h. If h=0.D0 it is selected automatically. * TS(3) Final time, tfinal. Program won't integrate past this time. * TS(k>3) If G-Stops are used G_i is stored in TS(i+3). The value of * TS(1) associated with previous G's is store in the next * location, and the previous G's follow it. The user need not * touch these latter locations unless they are doing something * unusual such as changing the definiation of G. * Y [in/out] Vector containing the current values of the dependent * variables, y. Dimension of NEQ, see IDAT(4). * OPT [in] This vector is used for the specification of options. * Options are indicated with a floating point integer followed by 0 or * more arguments. The arguments are indicated in parentheses in the * order in which they must be given. If an argument is followed by an * "=...", the "..." gives the default value for this option. Following * the arguments, the next option is specified. * Options are indicated with integers below, even though they must * of course be supplied as floating point numbers. For clarity in * specifying these options one may want to use the names for these * parameters. The start of the declarations below gives declarations * for defining these parameters. Other parameters used for accessing * IDAT and DAT are defined there also. If one is accessing these * arrays we recommend that it be done using these parameters. * * XRKEND=0 End of options. * * ************* Options applying to groups of equations. ************** * An Option which refers to a group of equations applies to all groups * after the option is specified, until another option overrides it. * When another option overrides a given option, the prior option * applies to equations with indexes from (the last value of A1 when * the option was specified + 1) to (the last value of A1 when the * option was overridden). There is an implied A1=0 at the start of * the option processing. * * XRKEQG=1 (A1=NEQ) Index of last equation in a group of equations. If * the end of the options is reached without getting this option, with * A1=NEQ, then it is as if this option appeared with A1 = NEQ = */ /* IDAT(4) just before the end of the options. This option is only * necessary if one wants to treat groups of equations differently in * some way. If this option is specified more than once, the values * of A1 must be strictly increasing. * * Following the practice in Hairer, Norsett, and Wanner, referenced * above, when referring to orders, we are referring to global order. * The local order is 1 greater than the global order. To generate an * error estimate, this code obtains E3 and E5 of orders 3 and 5 * respectively. We want a estimator of order 8. This is obtained by * forming h * E5 * (E5 / E3), where h is the stepsize. (Actually we * generate the squares of (these estimates / accuracy requested).) It * is this order 8 estimator that is used in controlling the stepsize. * (Actually E3 is (.01 * the average of recent value to what the error * would be with a third order method). When E5 / E3 is > 1 the code * uses h*E5 as the error estimate, and if this occurs too frequently * may give a diagnostic that the system appears stiff. This diagnostic * will only be given once. When forming error estimates we have * vectors v3_i and v5_i, i = 1, 2, ..., NEQ. For each of these v's we * form a sum of the form ||v||_e = sum_{i=1}^NEQ (v_i / tol_i)^2, where * tol_i is the absolute error tolerance being requested on the i-th * equation of the system. The code attempts to keep the order 8 error * estimate divided by the rquested accuracy times DAT(LERTRY) (nominal * value 100) at 1, and repeats the step if the square toot of this * quatity is >DAT(LERTRY+1) (nominal value is 4). If using an absolute * error tolerance or a relative error tolerance when the solution may * get very close to zero, one should insure that the absolute error * part of the tolerance does not underflow or overflow, when it is * squared. No checks are made to insure that this is the case. * * XRKAE=2 (A2) For the group in which this error tolerance applies, *tol_i (see just above) is A2. XRKRE1=3 (A3=e**(.75D0), B3=e**(.75D0)) *(tol_i**2) is (A3**2 + (B3**2) * (sum in group of (y_i)**2)). Note *that this is the default error test with A3 = B3 = e, where e is the *smallest number > 0 such that the floating point addition 1 + e gives a *result different from 1. XRKRE2=4 (A4, B4) tol_i**2 is (A3**2 + *(B3*y_i)**2). XRKAEV=5 (A5) Similar to option 2, except "A2" is now a *vector, with values stored starting in DAT(A5), where A5 > LOCTY + *2*NEQ. (The parameter LOCTY has the value 30.) This has the same *effect as option 2 followed by option 1 with A1 increased by 1.D0 each *time, for the number of times there are equations in this set. Space *used in DAT for storing this kind of error control information should *be contiguous and specified with increasing values of A5. XRKREV=6 *(A6) as for option 5, except applied to option 4. Thus there are *values ("A3", "B3", "A3", ...) stored starting at DAT(A6). If both *options 5 and 6 are used, the space used by the vectors for option 6 *must not overlap with the space used by vectors for 5. XRKNOE=7 Don't *accumulate anything into the measure of the error for equations in this *set. XRKDIG=8 (A8=0.D0) A8 = 0.D0, turns off diagnostic output for an *equation group, =1.D0 turns on output giving information used in *computing the error measure for the groups. XRKIN=9 (A9=1.D0) A9 = *1.D0 turns on interpolation, = 0.D0 turns it off. If one is doing *interpolation to an output point, no valid results are returned for *equations which are in groups with A9 = 0.D0. If you have G-Stops, *make sure that all the y's used in computing them are interpolated for. *XRKXP=10 (A10=0.D0) A10 = 1.D0 causes extra precision to be used in *accumulating Y, =0.D0 uses standard floating point precision. (This *option implies that extra precision will also be used in accumulating *TS(1). If one only wants the extra precision in TS(1), follow this *option with option 1, and A1 = 0.D0, and then option 10 with A10=0.D0.) *XRK???=11 Reserved (???) ********** End of Options applying to groups *of equations. *********** * * Other options * XRKFR1=12 Use reverse communication for computing f(t,y), rather than * having a call made to DXRK8F. When a return is made to the calling * program requesting an evaluation of f(t,y), the user should store * f(t, y) in F, and call DXRK8A. * XRKFR2=13 Use reverse communication for computing f(t,y), rather than * having a call made to DXRK8F. Store f(t, y) in WORK(IDAT(2)). * This saves copying F to WORK(IDAT(2)) inside DXRK8A. * XRKOR=14 Use reverse communication for output rather than having a * call made to DXRK8O. * XRKMAX=15 (A15=|tfinal-current value in TS(1)|, tfinal = TS(3)) A15 * gives absolute value of the maximum stepsize allowed. * XRKMIN=16 (A16=0.D0) A16 gives the minimum absolute stepsize allowed. * XRKEPS=17 (A17=mach_eps) A17 gives the precision of the floating point * arithmetic. If one knows that there is a loss of precision in * computing the derivatives, one may want to set this to a value * larger than the default value. * XRK???=18 Reserved (Use arc length as independent variable?) * XRK???=19 Reserved (???) * * XRKMXS=20 (A20) A20 = maximum number of steps between output points. * If this option is not set, 100,000 steps are taken before giving * output with IDAT(1) = 20. * XRKEOS=21 Give output at end of all steps. * XRKTI=22 (A22) Interpolate the solution at t = A22. A new value may * be set for the output point when output at this point is indicated. * XRKTS=23 (A23) As for option 22, except return value at the end of * the first step at or past A23 point. * XRKTDT=24 (A24, B24) As for option 22, except after giving the output, * if one does not change the output point, the value of B24 is added * to the current output point to get a new one. * XRKXG=25 (A25). A25 gives the number of extrapolatory G-Stops. * XRKIG=26 (A26). A26 gives the number of interpolatory G-Stops. * XRKEC=27 (A27=4.6D0, B27=.16D0) Logic in the code attempts to select * the stepsize so that -log(||Error Estimate||^2) is A27, and will * reject a step if ||Error Estimate||^2 > B27. ||Error Estimate|| is * computed by dividing error estimates by requested tolerances. */ /* These numbers will be adjusted upwards if the requested error can * not be obtained. * * XRKIFS=28 (A28=.1D0) Fraction of a step by which an interpolation * point can be outside the current integration interval. * * IDAT [in/out] Integer vector containing all the integer information * needed to do the integration. * IDAT(1), IDAT(2) and IDAT(3) define the state of the integration * when DXRK8 (or DXRK8A) is called, and when it makes a return to * the user, or calls one of the user coded routines. Changing * IDAT in any way other than described here gives undefined * results. IDAT(3) is only used as noted, and should not be * changed by the user. * IDAT(1) = 0 Starting an integration, not set on return from DXRK8. * IDAT(2) and IDAT(3) are not used when starting. * IDAT(1) = 1 F = f(T,Y) * If reverse communication is used replace T with TS(1). * If the reverse communication used is one that does not require * an extra copy, replace F with WORK(IDAT(2)). * If extrapolating G-Stops are being used, all such G's should * be computed prior to computing f, and if dxrk8g sets IDAT(2) * < 0, one should immediately return to the integrator. * IDAT(1) = 2 Computing G. (IDAT(3) is 0, except where indicated) * IDAT(2) = 0, compute G's indicated below and call dxrk8g. (In * many cases your code will be cleaner if you simply compute * all of the G's in all cases and this is allowed.) * IDAT(3) is used as follows in this case. * = 0 First check, compute all interpolating G's. * > 0 A zero of G(IDAT(3)) has been found, compute all G's. * =-1 Check of extrapolating G-Stops on partial advanced * step, compute all G's. * IDAT(2) > 0, the index of the G whose 0 is currently being * iterated for. Compute at least this G, and call dxrk8g. * IDAT(3) is 0 if in the middle of an iteration, and is the * same as IDAT(2) if checking ahead for sign persistence. * IDAT(2) = -1 computation of G's is complete, return to dxrk8. * IDAT(3) gives the index of the G_i which just crossed a 0. * IDAT(2) = -2 No iteration of G at the end of this step was * necessary, simply return to DXRK8. * IDAT(2) =-3 interpolation must be set up, and an iteration * started for the G with the index in IDAT(KEOS). If * IDAT(KEOS) is -1 we are setting up to extrapolate. * IDAT(2) =-4 the stepsize needs restricting for an extrapolating * G-Stop. WORK(IDAT(LOCINT)-3) contains the multiple of the * last stepsize to use. * IDAT(1) = 3 Flags some kind of termination condition. * IDAT(2) = 0, the integration has reached TS(3) = tfinal. If * TS(3) is changed the integration can be continued. (Doing * this sort of thing only makes sense if no option that requires * interpolation has been specified.) * IDAT(2) = -1 Set by the user (when processing an output point) * to interrupt the integration at the end of the step in a state * such that the integration can be continued. At the completion * of the step, a return is made with IDAT(1) = 3, and IDAT(2) = * 1. The value of WORK need not be preserved when continuing * with IDAT(1) = 3, and IDAT(2) = 1. * IDAT(2) = 1 Continue an integration that was interupted when * IDAT(1) > 19, and IDAT(2) was set to -1. * IDAT(2) > 1 Works as does IDAT(2) = 0, except this value has * results when the user set IDAT(1) = 3 and IDAT(2) < -1. DXRK8 * then sets IDAT(2) to its absolute value and returns to the * user. The integration can not be continued after this point. * If dxrk8o is being called for output points then it will be * called with IDAT(1)=3, and IDAT(2) = 0 upon reaching the final * point. Upon return from dxrk8o, if TS(3) has not been changed, * a return to the calling program will be made with these same * values for the flags. * IDAT(1) = 4 An error has occurred. IDAT(2) identifies what. * IDAT(2) = 1 TS(1) = TS(3) at start. * = 2 Error estimate too large, and at min. stepsize. * = 3 Not used. * = 4 Looks like problem is getting stiff * = 5 Looks like too much precision has been requested. * On this return IDAT(3) is 0. If it is not changed, * the code will change the accuracy requested to what * it regards as a more reasonable value. * = 6 An extrapolating G-Stop encountered upon taking a * step that didn't change TS(1). * = 7 Precision limited by noise in computation. *# (Code has no way to set this now. Such a test would * probably be based on computing the norm of f, and * flagging condition when (E3/||f||)**(2./3.)< IDAT(4)) * = 13 A repeated spec. for the same group of equations. * = 14 An undefined option. * = 15 Entered with an invalid value for IDAT(1). * = 16 Trying to continue after a fatal error reported. * = 17 Absolute/Relative error tolerance vectors overlap. * = 18 Need to set aside more space in IDAT. * = 19 Need to set aside more space in DAT. * = 20 Need to set aside more space in WORK. * = 21 Option required inside equation group, isn't. * = 22 Bad value for IDAT(1) on entry to DXRK8A. */ /* = 23 Initial TS(2) has the wrong sign.@'')') * = 24 DXRK8G was not called when requested. * IDAT(1) = 5 Restarts the integration. Ordinarily it is better * do the restart using IDAT(2) < 0, as described for IDAT(1)>19. * IDAT(1) > 19 Values of T and Y have been set for an output point, * and IDAT(1) gives the integer index of the option that resulted * in this output point. If 21 < IDAT(1) < 24, then IDAT(2) gives * the order in which this option was declared among the options * in this range and DAT(IDAT(3)) contains the value of t * associated with the output. Thus if one has a large number of * different output points defined, one has a way to tell * specifically which is responsible for the current output, and * can change the value in DAT(IDAT(3)) (and/or DAT(IDAT(3)+1) * when IDAT(1) = 24). In the case of other values of IDAT(1)>19, * IDAT(2) is 1, except in the case of G-Stops when IDAT(2) gives * the index of the g function which has a zero. To get a return * to the program calling DXRK8 with IDAT(1) = 3, and IDAT(2) = * -1, simply set IDAT(2) = -1. If one would like to restart the * integration at the current value of T and Y, set IDAT(2) = -2 * to use the same stepsize, and IDAT(2) = -3 to get the new one * picked automatically. If one would like to interpolate * immediately to an arbitrary value of T, set T (or TS(1)) to the * interpolation point and IDAT(2) = -4. The values at the * interpolated point will be returned with IDAT(1) = 20, IDAT(2) * = 1, and IDAT(3) = 0. One can continue setting IDAT(3) = -4 as * many times as desired. * IDAT(4) [in only] Contains NEQ, the number of differential * equations. Dimension of Y and F (if used) must be at least * this large. * IDAT(5) [in/out] The space at the beginning of IDAT that may be * referenced by DXRK8. If one provides a negative number here, * its absolute value is used, and the value in IDAT(5) is * replaced by the space in IDAT actually required by DXRK8. * LOCERR is a parameter with a current value of 44, Let I8 = * number of times option 8 (XRKDIG) has been turned "on" and I9 * and I10 be defined similarly for options 9 (XRKIN) and 10 * (XRKXP). Let IE be the number of times options 2-7 (error * control) have been specified (including the implied * specification for the first group of equations if an explicit * specification is not given) and IO the number of options for * I/O (21 to 24) that have been specified. Then one should have * IDAT(5) .ge. LOCERR + 10 + 2*(I8+I9+I10+IO) + 3*IE. The * actual space needed may be smaller by 1 or 2. If you are * concerned with such details, run once with IDAT(5) < 0. If you * are concerned that minor future perturbations to your code in * the future may result in slightly more storage being required, * make these numbers a little larger than required. * IDAT(6) [in/out] As for IDAT(5), except for DAT. Make IDAT(6) .ge. * LOCTY (=30) + 2 * IDAT(4) + space for all error tolerances + * space for defining t output values and delta t for option 24 * (XRKTDT) + if using extra precision, 1 for t and for each y * that is kept in extended precision. * IDAT(7) [in/out] As for IDAT(5), except for WORK. Define LWI = * IDAT(4) if there is no interpolation being done to arbitrary * points, and LWI = max(IDAT(4), 6 + 8 * (number of components of * Y interpolated to arbitrary points (default is IDAT(4))). Make * IDAT(7) .ge. 9*IDAT(4) + LWI. * * Starting in IDAT(10) are integer variables used internally. These * locations are referenced in the code with IDAT(), where * is a parameter name. Such names are defined with IDAT under the * section "Variable Definitions" below. Beyond this area of IDAT is * space used for the processing of various options. * * DAT As for IDAT, except starting in DAT(1) are floating point * data needed during the integration. This data must be preserved if * an integration is interrupted by some other computation. As for * IDAT there are parameter names defined below for accessing scalar * data in DAT. * WORK Data that need not be saved from one step to the next if the * integration is interupted for another task. This data must be * preserved if the integration is interrupted by some other means than * using option 3. See WORK under "Variable Definitions" for more * details. * * ************************* Variable Definitions *********************** * * Note that names used in subscripts of DAT and IDAT are defined in * alphabetical order under DAT and IDAT respectively. * DAT Formal argument, see above. Below is an alphabetical listing * of parameter names that are used in referencing DAT. Thus to find * what reference to DAT() is about, one should start by looking * for in the list below, before looking for the parameter name * directly. The parameter declarations for these names gives the * values of the associated parameters. The user should only alter * values in these locations of DAT with careful reflection on all * possible consequences and after verifying that the values of the * parameters used in the code are what they expect them to be. The * end of DAT is used for data required for the processing of options. * * Contents of DAT() * HA contains the absolute values of the stepsize. * LAH The value of AH used in defining the current stepsize. * LASTE3 contains last estimate for the 3-rd order error estimate * (multiplied by a factor and not including the factor of h). * LEPS contains EPS = smallest number such that 1.D0 + EPS .ne. 1.D0. * Option XRKEPS=17 can be used to change the above value. * LERR DAT(LERR:LERR+2) contain information associated with * stepsize control. Exponential average is based on the * parameter W. */ /* DAT(LERR) = "constant" part of exponential average of errors. * DAT(LERR+1) = "linear" part of exponential average of errors. * DAT(LERR+2) = "quadratic" part of exponential average of errors. * This area is also used to store some values when computing the * initial stepsize. * LERR1 =1 if the first equation index with error control is 1, else it * is 2. * LERTRY An attempt is made to keep DAT(LERTRY) * ||estimated error / * requested error|| at 1. If this is > DAT(LERTRY+1)the step is * rejected. If there is no error control, DAT(LERTRY) = 0.D0 * LHMAX Current maximum stepsize. (= DAT(LHMAXU) unless the code * thinks there is a risk of stability problems. * LHMAXU contains the maximum stepsize (= inteval length if not set by * the user). * LHMIN contains the minimum stepsize. * LINT IDAT(IDAT(LINT)) to IDAT(IDAT(LINT)+1) contain data defining * the current T interpolation point. If IDAT(LINT) is 3, we are at * the final output point. * LINTX When interpolating, the interpolation point is allowed to be * outside the interval of the current step by DAT(LINTX) * |the * stepsize|. * LOCSAV This and following locations are used to save a few floating * point quantities when a recoverable error return is made. * Also used to save output T for user interpolation interupts. * LOCTGS contains TS(1) at G-Stop that is preceded by a T-Stop. * LOCTY The value of t at the beginning of the current step. The * vector y at the beginning of the current step follows * immediately, and the contents of f at the beginning of the * step follow immediately after y. * LOCTXG The value of T we want to use when extrapolating to check on * extrapolating G-Stops. * LPHI last value computed for PHI * LSAVEH value of stepsize saved when taking small step to final t or * for an extrapolating G-Stop. * LTOUT value of the independent variable at the next output point. * NEGPOS -1.D0 if the stepsize if < 0, and 1.D0 otherwise. * NEXTH value to be used for the stepsize on the next step. * NEXTT value to be used for the base time on the next step. * PHI The current value for phi, the coefficient expected for * the error term. */ /* ************************* End if DAT parameters ********************** * * DATERR Used for temp. storage of floating pt. data for error messages. * DIFFTS The difference TS(3) - TS(1). * EN Passed to DXRK8N. EN(1) is the square of the fifth order error * estimate, and EN(2) is the square of the third order estimate. * ERR Ratio of (error estimate / error request)**2. * ERR1 Used for storing error option values when processing options. * ERR2 Same as for ERR1, but value from next location in OPT. * ERR3 Starts as square of the 3rd order error estimate. * ERR5 Starts as squere of the 5th order error estimate. * GOTINT True if interpolation might be required, false otherwise. * I Temporary integer. * I1 Temporary integer. * I2 Temporary integer. * II Temporary integer. * IBEG Index of next location available for storing option data in * DAT. * IDAT Formal argument, see above. Below is an alphabetical listing * of parameter names that are used in referencing IDAT. Thus to find * what reference to IDAT() is about, one should start by looking * for in the list below, before looking for the parameter name * directly. The parameter declarations for these names give the * values of the associated parameters. The user should only alter * values in these locations of IDAT with careful reflection on all * possible consequences and after verifying that the values of the * parameters used in the code are what they expect them to be. The * end of IDAT is used for data required for the processing of options. * * Contents of IDAT() * INTSAV Points to the initial location for saving information about * output points * KDZERO Used in DXRK8G for flag used in call to DZERO. * KEOS Gives the state of processing at the end of a step. Set from * IDAT(LOCEOS). Possible values are defined there. This is * also used to hold the index of the current G-Stop when * iterating for a G-Stop and set to -1 when extrapolating * because of an extrapolating G-Stop. * KHMIN Number of times we have had H = DAT(LHMIN) and error to big. * KFAIL Starts at 0 and incremented by 1 on each step failure. * KSTEP Number of (accepted) steps taken. * KSTEPM Value for IDAT(KSTEP) at which a user return is called for. * On such a return, IDAT(KSTEPM) is increased by IDAT(KSTEPM+1). * KSTEPX When IDAT(KSTEP) .ge. IDAT(KSTEPX), we have enough accepted * steps so that something special needs to be done. Set to -1 * if user want to save the solution. * KSTIFF Count of times in a row we think stiffness is an issue. * KWORK Amount of space used in WORK. * LAST0 Index of last location in IDAT to set 0 initially. * LDAT Length of DAT available = | IDAT(6) | * LIDAT Length of IDAT available = | IDAT(5) | * LOCEOS Gives the initial value assigned to IDAT(KEOS) when done with * a step. * = -3 End of step then G-Stop Check * = -2 End of step then TOUT * = -1 G-Stop Check then TOUT * = 0 TOUT then finish * = 1 Finish End of step actions. */ /* LOCEPR If nonzero, points to a list: (i1, j1), (i2, j2), ..., 0. * Error estimates are printed for equations i_k to j_k for * k = 1, 2, ... * LOCERR This must be the last location specified as a constant in IDAT. * Starts a list of actions defining how to compute the square * of the norm of a vector relative to the error tolerance. This * list has the form k_1, e_1, l1, k_2, ...., e_?=NEQ, l_?. The * k_i correspond to options 2-7 for specifying error criteria, * the e_i give the indexes for the last equation for which the * specification applies, and the l_i (not used if k_i = 7) point * to where the error tolerance information is actually stored in * DAT. (When the error information is a vector, the value of * l_i is such that when the index of the equation (or 2 times * the index for the case of option 6) when added to l_i gives * the value to use for the error tolerance. * LOCINT Points to area in WORK where information required to do * interpolations is stored. Let NIY = IDAT(NUMINT), t0 and t1 * be the values of t at the beginning and end of the last step * taken, and y[t0], y[t0,t1], ... y[t0,t1,t0,t1,t0,t1,t0,t1] * denote estimated divided differences of y. The first NIY * locations contain the values of y[t0], required, the next NIY * locations contain the values of y[t0,t1] required, etc. Let * LI = IDAT(LOCINT), then WORK(LI-3:LI-1) are used to contain: * Value for t at which an interpolation is to be made, value of * t at beginning of the step with respect to which the * interpolation is being done, and the (signed) length of the * step. In addition in the case of G-Stops, WORK(LI-6) is used * to save a value of G used in the iterations, WORK(LI-5) is * used for the corresponding value of T, and WORK(LI-3) is * used to store the value of T associated with a possible zero. * LOCOUT Points to a list of actions needed to find the next output * point after informing the user that a given output point has * been reached. This list has the form: i_1, k_1, i_2, k_2, ... * 30, 0, where i_n gives the index for the option, and DAT(k_i) * contains the value of the output point. This list is always * kept ordered so that DAT(k_n) should be output before or is * the same as DAT(k_{n+1}). * LOCXP If nonzero, points to a list: k, (i1, j1), (i2, j2), ..., * NEQ, where DAT(K) contains the extra precision part of TS(1) * and succeeding locations contain the extra precision parts of * Y that are needed, and extra precision is maintained in Y for * equations i_k to j_k for k = 1, 2, ... If this option is * specified, extra precision is always maintained for TS(1). * LREVF If 0, call DXRK8F for computing f. * If -1, use reverse communication for computing f, user stores * f in WORK(IDAT(3)). * If 1, treat as for -1, except user stores f in F. * LREFO If 0, call DXRK8O for going output. * If 1, use reverse communication for output. * LSTLAS Value of IDAT(KSTEP), the last time a step was taken with * order 3 error estimate relatively small relative to the order 5 * error estimate. * LTFIN Location in the list starting with LOCOUT for the final point. * LWORK Length of WORK available = | IDAT(7) | * LWHYF Index giving the current reason for computing the derivatives. * 1-12 Identifies the stage of the Runge-Kutta process. * 13-15 Used to obtain the extra derivatives required when * computing the data needed for doing interpolation. * 16 At the initial point. * 17 Extra derivative used in automatic selection of the initial * stepsize. * 18 Just got last derivative required to complete a step. * 19 Got a G sign change, take a small step prior to iterating. * 21 Used as flag to let us know dxrk8g was called at end of step. * NERRSV = -2 on the first step, -1 on the second, 1 just after a * rejected step after getting started and 0 othewise. * NGHI Absolute value gives the highest index of a G-Stop found at a * given value for TS(1). See NGLO below. * NGLO As for NGHI, but the lowest indexed G-Stop for a given value of * TS(1). The highest indexed G-Stop is reported to the user * first. Signs on NGHI and NGLO are used as follows. * NGLO NGHI * 0 No G-stops active at this value of t. * >0 >0 Interpolation is giving values to output. * >0 <0 Iterpolation for t output followed by a G-Stop. * <0 >0 Interpolation to start iterating for interpolating G-Stop. * <0 <0 Extrapolation to start iterating for extrapolating G-Stop. * NTGS Starts equal to the number of interpolating G-Stops. Later the * index in TS of the last G-Stop, 0 if none. * NUMINT The number of y's computed when interpolating. If data needs * to be computed before doing an interpolation, IDAT(NUMINT) is * negated. * NXGS Starts equal to the number of extrapolating G-Stops. Later the * index in TS of the last extratpolating G-Stop. * ************************ End if IDAT parameters ********************** * * IERR Internal index for an error message. * IERROR Array used to map IERR to the external error messaged index. * ISTOP Array used to map IERR to a stopping index for error messages. * ITEST Index of last location in current area for saving data in * DAT. If no more space is available, this is set to 0. * J Temporary integer. * K Temporary integer. * KACT An array of integers used in processing options. The contents * of KACT(i) = 10*c + s, where * s is the space required in OPT to specify option i. * c is the category of option i, with the following possibilities: * =0 End of options -- May need final actions for some categories. * =1 Equation group -- Special actions for categories 2-5. * =2 Category for tracking error control specs. See LOCERR. */ /* =3 Category for tracking diagnostic output on error estimation * for individual equations. See LOCEPR. * =4 Category for tracking which dependent variables are to be * maintained with extra precision. See LOCXP. * =5 No longer used. Might be used later? * =6 Category for tracking actions to take to determine the next * value of the independent variable which is an output point. * See LOCOUT. * =7 Option specifying an action that takes place on the end of a * step. * =8 Option not defined. * >10 Option that sets a value in IDAT(c) (value is 1 if s=1). * >1000 Option that sets values in DAT(c-1000) & if s=3 in DAT(c-999). * Data in KACT is set as follows * Opt 10*c+s Brief Description * 0 1 End of options. * 1 12 (A1=NEQ) Index of last equation in a group of equations. * 2 22 (A2) A2 gives an absolute error test tolerance * 3 23 (A3=e**(.75D0), B3=e**(.75D0)) Relative error * 4 23 (A4, B4) Similar to option 3 * 5 22 (A5) A5 is pointer to absolute error data. * 6 22 (A6) A6 is pointer to relative error data. * 7 21 Don't accumulate anything into the measure of the error * 8 32 (A8=0) A8 = 0, turns off diagnostic output, =1 turns on * 9 52 This category is no longer used. It might be used later? * 10 42 Use extra precision in accumulating Y. * 11 81 Reserved (???) * 12 10*LREVF+1 Reverse communication for computing f(t,y). * 13 10*LREVF+1 Rev. comm. for computing f(t,y), no copy needed. * 14 10*LREVO+1 Reverse communication for output. * 15 10*(LHMAXU+100)+2 (A15=TS(3)-value in TS(1)) A15 gives hmax. * 16 10*(LHMIN+100)+2 (A16=0) A16 gives hmin. * 17 10*(LEPS+100)+2 (A17=mach_eps) Floating point precision. * 18 81 Reserved (Use arc length as independent variable?) * 19 81 Reserved (???) * 20 10*KSTEPM+2 (A20) A20 = max number of steps between output points. * 21 71 Return at end of every step. (Sets IDAT(LOCEOS) = * IDAT(LOCEOS) - 2 if IDAT(LOCEOS) > -2) * 22 62 (A22) Interpolate the solution at t = A22. * 23 62 (A23) As option 22, except return value at the end of step. * 24 63 (A24, B24) A24 gives the increment between output points * 25 72 (A25). A25 gives the number of extrapolatory G-Stops. * This and option below, subtract 1 from IDAT(LOCEOS) if * IDAT(LOCEOS) is 0 or -2 * 26 72 (A26). A26 gives the number of interpolatory G-Stops. * 27 10*(LERTRY+100)+3 (A27=100, B27=4) Params. for error control * 28 10*(LINTX+100)+2 (A28=.1D0) Fraction of step by which * interpolation can be outside an interval. * KAT Index giving the category. See KACT above for categories. * KRK??? Integer values for the parameter names starting with XRK.... * KSTATE Array indexed by KAT giving states in the processing of * options. Usually KSTATE(KAT) is 0 for an option that is off and 1 * for an option that is on. * KUSE When first checking options this accumulates the number of * times the option is set. On the second pass this points to where * data defining the option is to be stored. Indexed by KAT, the * category of the option. * L Temporary integer. * LERMSG Array used to hold the locations of error messages. * LEQ Index of last equation for a group when processing options. * LEQL Array indexed by the category that give the value of LEQ the * last time that category was used. * LLERR In initialization gives LO on the last error control option. * LO Index in OPT in the initialization. * LOCBF F(i) at he start of a step is saved in DAT(LOCBF+I). * MACT Array used to hold error messaged actions. * MACT1 -- MACT3 As for MACT. * MAXLER MAXLER(-5) contains the largest index in DAT required for * storing absolute error tolerance vectors. * MAXLOC 1 + Index of last location used in DAT for storing vectors used * for error control. If no such vectors are used, this is 0. * MEEMES Paramater used to flag the start of an error messaged. * METEXT Parameter used to output text for an error messaged. * MEFDAT Parameter giving used to set start of floating point data for * error message processor. * MEFVEC Param. used to output floating pt. vector for error messages. * MEIVEC Parameter used to output an integer vector for error messages. * MENTXT Parameter giving the value to set index of next text to output. * MERET Parameter used to flag the end of actions for error messages. * MINLER As for MAXLER, except for the smallest index in DAT that is * used. * NEQ Number of equations. * NEQNOE Keeps count of number of equations with no error control. * Used for getting DAT(LERTRY) * NEXT Track next location available in IDAT when processing options. * NIY The number of components of y that are being interpolated for. * TP Temporary floating point quantity. * TP1 Temporary floating point quantity. * TP2 Temporary floating point quantity. * XRK??? Used for parameter names of the values used to define options. * WORK Formal argument, see above. The data required in WORK depends * on where we are in the process of taking a step. Let Ki, i = 1, 2, * ... denote the derivatives computed in the process of taking a * step, where K1 is the derivative at the beginning of the step, and * is stored in DAT(LOCBF+1:LOCBF+NEQ), LOCBF = LOCTY + NEQ. (Ki) is * used to indicate that Ki is not used on the current stage, but will * be needed later. Let y denote the solution at the beginning of the * step, y* the solution at the end of the step, and dy = y* - y. Let * Si denote the state when computing the i-th y value of the step, * where the first y value is kept in DAT(LOCTY+1:LOCTY+NEQ). Let dy */ /* and err denote computations used to compute dy and the estimated * error made on the step. Let I1, I2, I3, I4 be similarly defined for * the stages involved in computing the values necessary to do * interpolation. Let "Next" denote what is needed in making the * transition to a new step. Let err1 and err2 denote the 3rd and 5th * order error estimators, and "Int." denote space used to store data * needed for interpolation. Finally let H1 and H2 denote the stages * involved in getting an initial stepsize. * * WORK(*) Ni = i*NEQ Ni+ = Ni+1 DAT * N0+: N1+: N2+: N3+: N4+: N5+: N6+: N7+: N8+: N9+: LOCTY+ LOCBF+ * N1 N2 N3 N4 N5 N6 N7 N8 N9 N10 * S2 y K1 * S3 K2 y K1 * S4 K3 y K1 * S5 K3 K4 y K1 * S6 K4 K5 y K1 * S7 K4 K5 K6 y K1 * S8 K4 K5 K6 K7 y K1 * S9 K4 K5 K6 K7 K8 y K1 * S10 K4 K5 K6 K7 K8 K9 y K1 * S11 K4 K5 K6 K7 K8 K9 K10 y K1 * S12 K11 K4 K5 K6 K7 K8 K9 K10 y K1 * dy K11 K12 K4 K5 K6 K7 K8 K9 K10 (y) K1 * errs K11 K12 K6 K7 K8 K9 K10 (y) K1 * |err| dy err2 err1 (y) * I1 K11 K12 y* K13 K6 K7 K8 K9 K10 Int. y K1 * I2 K11 K12 (y*) K13 K6 K7 K8 (K9) K14 Int. y K1 * I3 K15 (y*) K13 K6 K7 K8 K9 K14 Int. y K1 * I4 K15 (y*) K13 K16 K14 Int. * Next [K12] y* K13 * H1 K1 y (K1) * H2 K2 y K1 * * When interpolating, WORK(1) = DAT(LINTX), and 3+8*NIY locations * beyond WORK(9*NEQ) are used to save information required for * doing the interpolation, where NIY is the number of components of * y that are interpolated for. If in addition the G-Stop feature * is used an additional 3 locations are required. * * ******************** Variable Declarations *************************** * * Formal arguments */ /* Parameters used for options */ /* Integer values for these parameter used internally. */ /* Parameters used to reference DAT. * Sorted by value: * 1 LEPS * 2 LHMAXU * 3 LHMIN * 4 LOCTGS * 5 LTOUT * 6 NEGPOS * 7 NEXTH * 8 NEXTT * 9 LINTX * 10 LSAVEH * 11 HMAX * 12 LAH * 13 LASTE3 * 14 LERTRY (Two locations) * 16 LERR (Three locations) * 19 PHI * 20 HLAST * 21 Not used * 22 HA * 23 LOCTXG * 24 LPHI * 25 HCAU * 26 LOCSAV * 27-29 Not used * 30 LOCTY */ /* Parameters used to reference IDAT. (Lots flagged as availble.) * Sorted by value: * 10 LINT * 11 LOCEPR * 12 LOCXP * 13 LBF * 14 LOCOUT * 15 LSTLAS * 16 NXGS * 17 NTGS * 18 LOCEOS * 19 KEOS * 20 KSTIFF * 21 KDZERO * 22 Available * 23 LREVF * 24 LREVO * 25 Available * 26 NERRSV * 27 KSTEP * 28 KFAIL * 29 Not used * 30 KHMIN * 31 NGLO * 32 NGHI (Also LAST0) * 33 Available * 34 KSTEPM * 35 Used along with KSTEPM * 36 KSTEPX * 37 Available * 38 LOCINT * 39 NUMINT * 40 LWHYF * 41 INTSAV * 42 LTFIN * 44 LOCERR (Must be last.) */ /* Local variables */ /* Parameters for setting value in KACT. */ /* Parameters for error messages */ /* Other Stuff for error processing */ /* ********* Error message text *************** *[Last 2 letters of Param. name] [Text generating message.] *AA DXRK8$B *AB Undefined option.$E *AC Equation indexes for groups are out of order or out of range.$E *AD Option required inside equation group is outside the last group.$E *AE Option specification made twice in the same group..$E *AF Error tolerance vectors overlap in DAT.$E *AG Options must be turned on with 1.D0 and off with 0.D0, you $C * used $F.$E *AH IDAT(5) must be at least $I.$E *AI Could not find space in DAT for all options. IDAT(6) must be $C * > $I.$E *AJ IDAT(7) must be at least $I.$E *AK The paramters used in option 27, must have the first at least 4 $C * times bigger than the second.$E * $ *AL OPT up to the point where we are seeing problems.$N * OPT: $B *AM IDAT(4:7):$B */ /* **** End of automatically generated text * 1 2 3 4 5 */ /* 1 2 3 4 5 6 7 */ /* 1 2 3 4 5 6 */ /* 7 8 9 10 * 1 2 3 4 5 6 7 8 9 10 */ /* 0 1 2 3 4 5 6 7 8 9 10 11 12 */ /* ********************** Start of Executable Code ********************** * */ if (Idat[1] != 0) goto L_600; neq = Idat[4]; lidat = labs( Idat[5] ); ldat = labs( Idat[6] ); lwork = labs( Idat[7] ); Dat[LEPS] = DBL_EPSILON; /* Set IDAT(2) to signal no errors. */ Idat[2] = 0; /* Initialize IDAT and DAT. */ for (i = LOCEPR; i <= LAST0; i++) { Idat[i] = 0; } diffts = Ts[3] - Ts[1]; Dat[NEGPOS] = sign( 1.e0, diffts ); Dat[LHMAXU] = fabs( diffts ); Dat[LHMIN] = 0.e0; Idat[KFAIL] = 0; Idat[KSTEP] = 0; Idat[KSTEPM] = 100000; /* ****************** Process options, First phase ********************** * * Processing is done in two phases. First phase checks for valid input * and determines how much space in needed for internal processing of the * options. The second phase stores data for later processing of an * option. * */ gotint = FALSE; for (i = 2; i <= 6; i++) { /* Indexes, I, here, correspond to the categories mentioned above. */ kuse[i-(2)] = 0; leql[i-(2)] = 0; kstate[i-(2)] = 0; } /* Need a special flag for extra precision so we know to turn it on. */ leql[2] = -1; kuse[4] = 2; /* Need to know if no error tolerances have been set. */ kstate[0] = -1; ibeg = LOCTY + 2*neq + 1; lerr1 = 2; minler = 0; maxler = ibeg - 1; /* Flags that no extra precision is used. */ Idat[LOCXP] = -1; leq = -2; lo = 1; l = 0; L_100: lo += l; k = nint( Opt[lo] ); if ((k < 0) || (k > KRKIFS)) { /* An undefined option */ ierr = 1; goto L_700; } j = kact[k]; /* Category of the option. */ kat = j/10; /* Space required by the option. */ l = j - 10*kat; if (kat <= 5) { if (kat <= 1) { /* Reached end of a group */ i = leq; if (kat == 0) { /* End of all options */ leq = neq + 1; goto L_130; } else { /* Set LEQ = index of equation for end of a group. */ leq = nint( Opt[lo + 1] ); } if ((leq > i) && (leq <= neq)) goto L_100; /* Equation indexes not strictly increasing or too big. */ ierr = 2; goto L_700; } else { /* Option needing extra space depending on when option 1 is specified. */ if (leq > neq) { /* Options can not be specified after the last group. */ ierr = 3; goto L_700; } if (leql[kat-(2)] == leq) { /* Two specifications for item in a group. */ ierr = 4; goto L_700; } if (kat == 2) goto L_130; j = nint( Opt[lo + 1] ); if ((j != 0) && (j != 1)) { Work[1] = Opt[lo + 1]; ierr = 6; goto L_700; } if (kstate[kat-(2)] != j) { if (kat == 4) { /* Get total space in DAT that is needed into IDAT(LOCXP) */ if (j == 1) { if (leql[2] <= 0) { Idat[LOCXP] = labs( Idat[LOCXP] ) + neq; } else { Idat[LOCXP] = labs( Idat[LOCXP] ) + neq - leq + 1; } } else if (leql[2] > 0) { Idat[LOCXP] += -neq + leq; } } kstate[kat-(2)] = j; if (j == 1) { kuse[kat-(2)] += 2; } else if ((kat == 5) && (leq <= 0)) { kuse[3] -= 2; } } leql[kat-(2)] = leq; goto L_100; } } else if (kat <= 7) { /* Interp. point, no interp. needed if end of step value is returned. */ if ((k != KRKTS) && (k != KRKEOS)) gotint = TRUE; if (kat == 6) kuse[4] += 2; } else if (kat == 8) { ierr = 1; goto L_700; } else if (kat > 10) { if (kat < 100) { /* Data to store in IDAT. */ if (l == 1) { Idat[kat] = 1; if (k == KRKFR2) Idat[kat] = -1; } else { Idat[kat] = nint( Opt[lo + 1] ); } } else { /* Data to store in DAT. */ Dat[kat - 100] = Opt[lo + 1]; if (l == 3) { Dat[kat - 99] = Opt[lo + 2]; if (Dat[kat - 100] < 4.e0*Dat[kat - 99]) { ierr = 10; goto L_700; } } } } goto L_100; /* For error control * End here if processing for errors is not finished. */ L_130: if (leql[0] > 0) { j = nint( Opt[llerr] ); if ((j == KRKAEV) || (j == KRKREV)) { /* Track space used by vectors for error control */ i = nint( Opt[llerr + 1] ); if (i > maxler) { maxler = i + leq - leql[0] - 1; if (minler == 0) { minler = i; } if (maxler <= ldat) goto L_140; } /* Problem with space for error vectors. */ ierr = 5; goto L_700; } } L_140: if (kat != 0) { if (leq == 1) lerr1 = 1; leql[0] = max( leq, 1 ); kuse[0] += 3; llerr = lo; lerr1 = 1; goto L_100; } if (lerr1 != 1) kuse[0] += 3; /* Set up starting locations for the second pass * Set start of available space in DAT and WORK */ kwork = 10*neq; /* Can't use space beyond DAT(ITEST). */ itest = ldat; if (minler != 0) itest = minler - 1; /* Set pointers for internal groups. * Space for error control. */ if (kuse[0] == 0) kuse[0] = 3; next = LOCERR + kuse[0]; kuse[0] = LOCERR; if (kuse[1] != 0) { Idat[LOCEPR] = next; /* Space for diagnostic print. */ next += kuse[1] + 1; kuse[1] = Idat[LOCEPR]; /* Set first location so we can recognize an off prior to the first on. */ Idat[kuse[1]] = 0; } if (Idat[LOCXP] >= 0) { /* K is the amount of space needed in DAT */ k = Idat[LOCXP]; Idat[LOCXP] = next; /* Space for extra precision specifications. */ next += kuse[2] + 2; kuse[2] = Idat[LOCXP] + 1; if (k <= itest - ibeg) { j = ibeg; ibeg += k; } else if (k <= ldat - maxler) { j = maxler + 1; maxler += j; } else { goto L_640; } /* Set extra precision parts to 0. */ Idat[kuse[2] - 1] = j; for (i = j; i <= (j + k - 1); i++) { Dat[i] = 0.e0; } } else { Idat[LOCXP] = 0; } /* Remember the initial start for the list of output points */ Idat[INTSAV] = next; /* And the current start for the above list */ Idat[LOCOUT] = next; /* Space for remembering output points. */ next += kuse[4]; kuse[4] = Idat[LOCOUT]; if (next > lidat) { /* Error -- Dimension of IDAT is too small */ ierr = 7; Idat[1] = next - 1; goto L_700; } if (Idat[5] < 0) Idat[5] = next; /* ****** Process the options -- Second phase -- Space looks o.k. ******* * */ for (i = 2; i <= 6; i++) { /* Indexes, I, here, correspond to the categories mentioned above. */ leql[i-(2)] = 0; kstate[i-(2)] = 0; } /* Need a special flag for extra precision so we know to turn it on. */ leql[2] = -1; neqnoe = 0; leq = 0; lo = 1; l = 0; L_300: lo += l; k = nint( Opt[lo] ); j = kact[k]; /* Category of the option. */ kat = j/10; /* Space required by the option. */ l = j - 10*kat; if (kat <= 5) { if (kat <= 1) { /* Reached end of a group and if KAT is 0, end of options. */ if (kat == 0) { if (lerr1 == 1) goto L_400; /* If no error control, we need to insert the default. */ lerr1 = 1; kat = 2; l = 3; k = 0; leq = neq + 1; } else { leq = nint( Opt[lo + 1] ); goto L_300; } } if (kat == 2) { L_320: j = kuse[0]; if (kstate[0] != 0) { Idat[j - 2] = leq - 1; } else { if (leq > 1) { /* Default actions for this group */ ii = KRKRE1; err1 = pow(Dat[LEPS],.75e0); err2 = err1; n = 2; goto L_330; } kstate[0] = 1; } ii = k; if (k >= KRKAEV) { /* Want 1 less to have a base address */ Idat[j + 2] = nint( Opt[lo + 1] ) - 1; if (k == KRKNOE) neqnoe = 1; goto L_340; } n = l - 1; err1 = Opt[lo + 1]; err2 = Opt[lo + 2]; /* Getting space in DAT to store data. */ L_330: if (ibeg + n > itest) { if (itest == ldat) goto L_640; itest = ldat; ibeg = maxler + 1; goto L_330; } Idat[j + 2] = ibeg; Dat[ibeg] = err1; if (n == 2) Dat[ibeg + 1] = err2; ibeg += n; /* Get here for vector case or no error control */ L_340: ; Idat[j] = ii; Idat[j + 1] = neq; kuse[0] += 3; if (kstate[0] != 0) goto L_300; /* If we had to set the defaults, we still need to do the current. */ kstate[0] = 1; if (k == 0) goto L_400; goto L_320; } else { /* Options we can turn on or off. */ ii = nint( Opt[lo + 1] ); if (kstate[kat-(2)] != ii) { if (ii == 0) { /* We are turning this one off */ Idat[kuse[kat-(2)] - 1] = leq - 1; } else { /* Turning this one on. Last is set to NEQ, reset if turned off. */ Idat[kuse[kat-(2)]] = max( 1, leq ); Idat[kuse[kat-(2)] + 1] = neq; leql[kat-(2)] = leq; kuse[kat-(2)] += 2; } kstate[kat-(2)] = ii; } goto L_300; } } if (kat == 6) { /* Save data for tracking output points */ n = l - 1; L_360: if (ibeg + n - 1 > itest) { if (itest >= ldat) goto L_640; itest = ldat; ibeg = maxler + 1; goto L_360; } Idat[kuse[4]] = k; Idat[kuse[4] + 1] = ibeg; kuse[4] += 2; Dat[ibeg] = Opt[lo + 1]; if (n == 2) Dat[ibeg + 1] = Opt[lo + 2]; ibeg += n; goto L_300; } else if (kat == 7) { if (k == KRKEOS) { /* Return at end of every step. (0 => -2, -1 => -3) */ if (Idat[LOCEOS] > -2) Idat[LOCEOS] -= 2; } else { /* G_Stop of some type (-2 => -3, 0 => -1) */ if (labs( Idat[LOCEOS] + 1 ) == 1) Idat[LOCEOS] -= 1; /* Save the count in the proper location */ Idat[NXGS + k - KRKXG] = nint( Opt[lo + 1] ); gotint = TRUE; } } goto L_300; /* ******************* Set final termination flags, Etc. **************** * * Finish up list of output points */ L_400: Idat[kuse[4]] = 30; Idat[kuse[4] + 1] = 0; Idat[LTFIN] = kuse[4]; /* Set final 0 for diagnostic output */ if (Idat[LOCEPR] != 0) Idat[kuse[1]] = 0; /* Get space in WORK used for interpolation */ if (gotint) { kwork += 7*neq + 6; Idat[NTGS] += Idat[NXGS]; if (Idat[NTGS] != 0) { Idat[NXGS] += 3; Idat[NTGS] += 3; } } Idat[NUMINT] = -neq; if (Idat[LOCXP] != 0) { Idat[kuse[2]] = neq; } /* Set IBEG to the space actually needed (instead of next available) */ ibeg -= 1; if (kwork > lwork) { ierr = 9; Idat[1] = kwork; goto L_700; } /* If requested, tell user about space needed. (IDAT(5) is above.) */ if (Idat[6] < 1) { Idat[6] = ibeg; } if (Idat[7] < 0) { Idat[7] = kwork; } /* Initialize some stuff that depends on user input. */ Idat[LOCINT] = 9*neq + 7; Idat[KSTEPM + 1] = Idat[KSTEPM]; Idat[KSTEPX] = Idat[KSTEPM]; Idat[LBF] = neq + LOCTY; Dat[HCAU] = Dat[LHMAXU]; Dat[HMAX] = Dat[LHMAXU]; if (neqnoe != 0) { neqnoe = 0; for (i = Idat[LOCERR]; i <= 1000000; i += 3) { if (Idat[i + 1] >= neq) goto L_540; if (Idat[i] == KRKNOE) { neqnoe += Idat[i + 4] - Idat[i + 1]; } } L_540: ; } Dat[LERTRY] = 0.e0; if (neq != neqnoe) { Dat[LERTRY] = 100.e0/(double)( neq - neqnoe ); } Dat[LERTRY + 1] = 6.e0; /* Set count on derivative evaluations. * * **************************** Integrate ******************************* * */ L_600: dxrk8a( ts, y, opt, idat, dat, work ); return; /* *********************** Error Processing ***************************** * * Entry when insufficient space in DAT */ L_640: ierr = 8; Idat[1] = ibeg; L_700: Mact[3] = Ierror[ierr]; Mact[4] = Lermsg[ierr]; dmess( mact, (char*)mtxtaa,182, idat, work ); Mactv[3] = lo + l; dmess( mactv, (char*)mtxtab,71, &Idat[4], opt ); Idat[1] = 4; Idat[2] = Mact[3]; return; } /* end of function */ /* PARAMETER translations */ #define A101 4.77662536438264365890433908527e-1 #define A104 (-2.48811461997166764192642586468e0) #define A105 (-5.90290826836842996371446475743e-1) #define A106 2.12300514481811942347288949897e1 #define A107 1.52792336328824235832596922938e1 #define A108 (-3.32882109689848629194453265587e1) #define A109 (-2.03312017085086261358222928593e-2) #define A111 (-9.3714243008598732571704021658e-1) #define A1110 (-3.0467644718982195003823669022e0) #define A114 5.18637242884406370830023853209e0 #define A115 1.09143734899672957818500254654e0 #define A116 (-8.14978701074692612513997267357e0) #define A117 (-1.85200656599969598641566180701e1) #define A118 2.27394870993505042818970056734e1 #define A119 2.49360555267965238987089396762e0 #define A121 2.27331014751653820792359768449e0 #define A1210 1.23605671757943030647266201528e1 #define A1211 6.43392746015763530355970484046e-1 #define A124 (-1.05344954667372501984066689879e1) #define A125 (-2.00087205822486249909675718444e0) #define A126 (-1.79589318631187989172765950534e1) #define A127 2.79488845294199600508499808837e1 #define A128 (-2.85899827713502369474065508674e0) #define A129 (-8.87285693353062954433549289258e0) #define A141 5.61675022830479523392909219681e-2 #define A1410 1.5329179827876569731206322685e-1 #define A1411 8.20105229563468988491666602057e-3 #define A1412 7.56789766054569976138603589584e-3 #define A1413 (-8.298e-3) #define A147 2.53500210216624811088794765333e-1 #define A148 (-2.46239037470802489917441475441e-1) #define A149 (-1.24191423263816360469010140626e-1) #define A151 3.18346481635021405060768473261e-2 #define A1511 (-1.08347328697249322858509316994e-4) #define A1512 3.82571090835658412954920192323e-4 #define A1513 (-3.40465008687404560802977114492e-4) #define A1514 1.41312443674632500278074618366e-1 #define A156 2.83009096723667755288322961402e-2 #define A157 5.35419883074385676223797384372e-2 #define A158 (-5.49237485713909884646569340306e-2) #define A161 (-4.28896301583791923408573538692e-1) #define A1613 (-1.39902416515901462129418009734e-3) #define A1614 2.9475147891527723389556272149e0 #define A1615 (-9.15095847217987001081870187138e0) #define A166 (-4.69762141536116384314449447206e0) #define A167 7.68342119606259904184240953878e0 #define A168 4.06898981839711007970213554331e0 #define A169 3.56727187455281109270669543021e-1 #define A21 5.26001519587677318785587544488e-2 #define A31 1.97250569845378994544595329183e-2 #define A32 5.91751709536136983633785987549e-2 #define A41 2.95875854768068491816892993775e-2 #define A43 8.87627564304205475450678981324e-2 #define A51 2.41365134159266685502369798665e-1 #define A53 (-8.84549479328286085344864962717e-1) #define A54 9.24834003261792003115737966543e-1 #define A61 3.7037037037037037037037037037e-2 #define A64 1.70828608729473871279604482173e-1 #define A65 1.25467687566822425016691814123e-1 #define A71 3.7109375e-2 #define A74 1.70252211019544039314978060272e-1 #define A75 6.02165389804559606850219397283e-2 #define A76 (-1.7578125e-2) #define A81 3.70920001185047927108779319836e-2 #define A84 1.70383925712239993810214054705e-1 #define A85 1.07262030446373284651809199168e-1 #define A86 (-1.53194377486244017527936158236e-2) #define A87 8.27378916381402288758473766002e-3 #define A91 6.24110958716075717114429577812e-1 #define A94 (-3.36089262944694129406857109825e0) #define A95 (-8.68219346841726006818189891453e-1) #define A96 2.75920996994467083049415600797e1 #define A97 2.01540675504778934086186788979e1 #define A98 (-4.34898841810699588477366255144e1) #define B1 5.42937341165687622380535766363e-2 #define B10 (-1.52160949662516078556178806805e-1) #define B11 2.01365400804030348374776537501e-1 #define B12 4.47106157277725905176885569043e-2 #define B6 4.45031289275240888144113950566e0 #define B7 1.89151789931450038304281599044e0 #define B8 (-5.8012039600105847814672114227e0) #define B9 3.1116436695781989440891606237e-1 #define BHH1 0.244094488188976377952755905512e00 #define BHH2 0.733846688281611857341361741547e00 #define BHH3 0.220588235294117647058823529412e-01 #define C10 0.6e00 #define C11 0.857142857142857142857142857142e00 #define C14 0.1e00 #define C15 0.2e00 #define C16 0.777777777777777777777777777778e00 #define C2 0.526001519587677318785587544488e-01 #define C3 0.789002279381515978178381316732e-01 #define C4 0.118350341907227396726757197510e00 #define C5 0.281649658092772603273242802490e00 #define C6 0.333333333333333333333333333333e00 #define C7 0.25e00 #define C8 0.307692307692307692307692307692e00 #define C9 0.651282051282051282051282051282e00 #define D41 (-0.84289382761090128651353491142e01) #define D410 (-0.87139158377797299206789907490e00) #define D411 0.22404374302607882758541771650e01 #define D412 0.63157877876946881815570249290e00 #define D413 (-0.88990336451333310820698117400e-01) #define D414 0.18148505520854727256656404962e02 #define D415 (-0.91946323924783554000451984436e01) #define D416 (-0.44360363875948939664310572000e01) #define D46 0.56671495351937776962531783590e00 #define D47 (-0.30689499459498916912797304727e01) #define D48 0.23846676565120698287728149680e01 #define D49 0.21170345824450282767155149946e01 #define D51 0.10427508642579134603413151009e02 #define D510 0.77334326684722638389603898808e01 #define D511 (-0.30674084731089398182061213626e02) #define D512 (-0.93321305264302278729567221706e01) #define D513 0.15697238121770843886131091075e02 #define D514 (-0.31139403219565177677282850411e02) #define D515 (-0.93529243588444783865713862664e01) #define D516 0.35816841486394083752465898540e02 #define D56 0.24228349177525818288430175319e03 #define D57 0.16520045171727028198505394887e03 #define D58 (-0.37454675472269020279518312152e03) #define D59 (-0.22113666853125306036270938578e02) #define D61 0.19985053242002433820987653617e02 #define D610 0.68812326946963000169666922661e01 #define D611 (-0.10006050966910838403183860980e01) #define D612 0.77771377980534432092869265740e00 #define D613 (-0.27782057523535084065932004339e01) #define D614 (-0.60196695231264120758267380846e02) #define D615 0.84320405506677161018159903784e02 #define D616 0.11992291136182789328035130030e02 #define D66 (-0.38703730874935176555105901742e03) #define D67 (-0.18917813819516756882830838328e03) #define D68 0.52780815920542364900561016686e03 #define D69 (-0.11573902539959630126141871134e02) #define D71 (-0.25693933462703749003312586129e02) #define D710 (-0.37458323136451633156875139351e02) #define D711 0.10409964950896230045147246184e03 #define D712 0.29840293426660503123344363579e02 #define D713 (-0.43533456590011143754432175058e02) #define D714 0.96324553959188282948394950600e02 #define D715 (-0.39177261675615439165231486172e02) #define D716 (-0.14972683625798562581422125276e03) #define D76 (-0.15418974869023643374053993627e03) #define D77 (-0.23152937917604549567536039109e03) #define D78 0.35763911791061412378285349910e03 #define D79 0.93405324183624310003907691704e02 #define ER1 0.1312004499419488073250102996e-01 #define ER10 0.3341791187130174790297318841e00 #define ER11 0.8192320648511571246570742613e-01 #define ER12 (-0.2235530786388629525884427845e-01) #define ER6 (-0.1225156446376204440720569753e01) #define ER7 (-0.4957589496572501915214079952e00) #define ER8 0.1664377182454986536961530415e01 #define ER9 (-0.3503288487499736816886487290e00) #define HA 22 #define HLAST 20 #define KEOS 19 #define KHMIN 30 #define KSTIFF 20 #define LAH 12 #define LASTE3 13 #define LERR 16 #define LINT 10 #define LOCSAV 26 #define LOCTGS 4 #define LOCTXG 23 #define LPHI 24 #define LSAVEH 10 #define LTOUT 5 #undef LTXTAC #define LTXTAC 55 #undef LTXTAD #define LTXTAD 86 #undef LTXTAE #define LTXTAE 127 #undef LTXTAF #define LTXTAF 169 #undef LTXTAG #define LTXTAG 286 #undef LTXTAH #define LTXTAH 324 #undef LTXTAI #define LTXTAI 460 #undef LTXTAJ #define LTXTAJ 512 #undef LTXTAK #define LTXTAK 1 #define LWHYF 40 #define MEDDIG 12 #define MEFDAT 25 #define NERRSV 26 #define NEXTH 7 #define NEXTT 8 #define NGLO 31 #define P 8.e0 #define PHI 19 #define PINV (1.e0/P) #define W .1e0 #define WI1A (W/powi(1.e0 - W,2)) #define WI1B ((1.e0 - 2.e0*W)/powi(1.e0 - W,2)) #define WI2A ((2.e0*W)/powi(1.e0 - W,3)) #define WI2B ((1.e0 - 3.e0*W)/powi(1.e0 - W,3)) #define WL1A ((1.e0 - SQ(W))/W) #define WL2A (-powi(1.e0 - W,2)/W) /* end of PARAMETER translations */ void /*FUNCTION*/ dxrk8a( double ts[], double y[], double f[], long idat[], double dat[], double work[]) { long int i, i1, ierr, j, k, locbf, neq; double ah, daterr[8], h, tp, tp1, tp2; static double en[2], err; static char mtxtaa[3][203]={"DXRK8$BIDAT(1) = $I is an error on entry to DXRK8A.$EToo much precision requested.$EDXRK8G was not called as was requested.$ECan not continue after integration done.$ELAST WARNING -- At T=$F, the Error", " Est. is $F which is too big, but the stepsize has its minimum value which is $F.$EInitial TS(2)=$F has the wrong sign.$ELAST WARNING -- Looks like problem is stiff near TS(1) =$ $F, with H = $F. Consi", "der using an integrator designed for stiff equations.$EIf you are going to ignore my error flags, I quit!$EA sign change in TS($I) at TS(1)=$F is skipped due to inconsistent computation of its value.$E "}; static char mtxtab[1][99]={"DXRK8$BError tolerance increased from $(E10.3) times$ that requested to $G times that requested.$E"}; static char mtxtac[1][91]={"$NRejected step at T=$F with stepsize of $(E12.5), ERR=$(E10.3), E5/E3=$G, AH=$F, PHI=$F$E"}; static char mtxtad[1][63]={"$NAH=$F, adjusted to $F with KSTIFF=$I, Last AH=$F, E5/E3=$F$E"}; static char mtxtae[1][72]={"T=$(E15.8) H=$(E12.5) ERR=$(E10.3) PHI=$G E5/E3=$G$ HMAX=$G HCAU=$G$E"}; static char mtxtaf[1][18]={"Raw($I:)/H**2 =$B"}; static long mact[5]={MEEMES,0,0,0,MERET}; static long lermsg[10]={LTXTAB,LTXTAC,LTXTAD,LTXTAE,LTXTAF,LTXTAG, LTXTAH,LTXTAI,LTXTAJ,LTXTAK}; static long ierror[9]={22,5,24,10,2,23,4,99,3}; static long istop[9]={88,25,88,88,25,88,14,99,14}; static long mact1[2]={METEXT,MERET}; static long mact2[6]={METEXT,MEFDAT,0,MEFVEC,0,MERET}; static long mact3[5]={MEEMES,24,0,0,MERET}; static long mact4[4]={MEDDIG,4,METEXT,MERET}; /* OFFSET Vectors w/subscript range: 1 to dimension */ double *const Dat = &dat[0] - 1; double *const Daterr = &daterr[0] - 1; double *const En = &en[0] - 1; double *const F = &f[0] - 1; long *const Idat = &idat[0] - 1; long *const Ierror = &ierror[0] - 1; long *const Istop = &istop[0] - 1; long *const Lermsg = &lermsg[0] - 1; long *const Mact = &mact[0] - 1; long *const Mact1 = &mact1[0] - 1; long *const Mact2 = &mact2[0] - 1; long *const Mact3 = &mact3[0] - 1; long *const Mact4 = &mact4[0] - 1; double *const Ts = &ts[0] - 1; double *const Work = &work[0] - 1; double *const Y = &y[0] - 1; /* end of OFFSET VECTORS */ /* Main action routine for solving differential equations with the * Dormand-Prince 8th Order Runge_Kutta formulas. * * ************************* Variable Declarations ********************** * * Formal arguments */ /* Parameters used to reference DAT */ /* Parameters used to reference IDAT */ /* Parameters used in exponential averaging for the stepsize adjustment. */ /* parameter (WI3A=(3.D0*W)/(1.D0-W)**4) * parameter (WI3B=(1.D0-4.D0*W)/(1.D0-W)**4) * double precision WQ1A,WQ2A,WQ3A,WQ1B,WQ2B,WQ3B,WQ1C,WQ2C,WQ3C * parameter(WQ1A=(1.D0-W)*(W**2+W+1.D0)/W**2) * parameter(WQ2A=(1.D0-W)*(W**2+W-2.D0)/W**2) * parameter(WQ3A=(1.D0-W)**3/W**2) * parameter(WQ1B=(1.D0-W)**2*(-W-2.D0)/W**2) * parameter(WQ2B=(1.D0-W)**2*(-W**2-3.D0*W+4.D0)/W**2) * parameter(WQ3B=-2.D0 * (1.D0-W)**4/W**2) * parameter(WQ1C=((1.D0-W)**3)/W**2) * parameter(WQ2C=(-2.D0*(1.D0-W)**4)/W**2) * parameter(WQ3C=((1.D0-W)**5)/W**2) */ /* parameter (WL1B=-(1.D0-W)**2/W, WL2B=(1.D0-W)**3/W) */ /* double precision AQ, BQ, CQ, AL, BL * These don't need to be in DAT as we only need them on rejected steps. */ /* Local Variables */ /* Integration coefficients */ /* Parameters for error messages */ /* ********* Error message text *************** *[Last 2 letters of Param. name] [Text generating message.] *AA DXRK8$B *AB IDAT(1) = $I is an error on entry to DXRK8A.$E *AC Too much precision requested.$E *AD DXRK8G was not called as was requested.$E *AE Can not continue after integration done.$E *AF LAST WARNING -- At T=$F, the Error Est. is $F which is too $C * big, but the stepsize has its minimum value which is $F.$E *AG Initial TS(2)=$F has the wrong sign.$E *AH LAST WARNING -- Looks like problem is stiff near TS(1) = $F, $C * with H = $F. Consider using an integrator designed for stiff $C * equations.$E *AI If you are going to ignore my error flags, I quit!$E *AJ A sign change in TS($I) at TS(1)=$F is skipped due to $C * inconsistent computation of its value.$E * $ *AK DXRK8$B *AL Error tolerance increased from $(E10.3) times that requested to $C * $G times that requested.$E * $ *AM $NRejected step at T=$F with stepsize of $(E12.5), ERR=$(E10.3)$C * , E5/E3=$G, AH=$F, PHI=$F$E * $ *AN $NAH=$F, adjusted to $F with KSTIFF=$I, Last AH=$F, E5/E3=$F$E * $ *AO T=$(E15.8) H=$(E12.5) ERR=$(E10.3) PHI=$G E5/E3=$G HMAX=$G $C * HCAU=$G$E * $ *AP Raw($I:)/H**2 =$B */ /* **** End of automatically generated text * * 1 2 3 4 5 * Other Stuff for error processing */ /* 1 2 3 4 5 6 7 8 9 */ /* 1 2 3 4 5 6 */ /*++ Default STAT=.F. *++ Code for STAT is inactive * double precision ESUMSQ, EMAXER * save ESUMSQ, EMAXER *++ End * * * ************************* Start of Executable Code ******************* * */ neq = Idat[4]; locbf = Idat[LBF]; if (Idat[1] == 1) { if (Idat[LREVF] > 0) { /* Copy F to place wanted in WORK. */ if (Idat[2] > 0) { for (i = 1; i <= neq; i++) { Work[Idat[2] + i - 1] = F[i]; } } } h = Ts[2]; } goto L_400; /* Computing derivatives -- Starting location depends on which deriv. */ L_200: Idat[2] += neq; L_210: Idat[LWHYF] += 1; L_220: ; if (Idat[LREVF] != 0) return; dxrk8f( &Ts[1], y, &Work[Idat[2]], idat ); /* ****************** Main Loop, Branch on value of IDAT(1) ************* * */ L_400: ; h = Ts[2]; /* 1 2 3 4 */ switch (Idat[1]) { case 1: goto L_1000; case 2: goto L_4320; case 3: goto L_3000; case 4: goto L_3100; } if (Idat[1] >= 19) goto L_3700; if (Idat[1] <= 0) { if (Idat[1] == 0) goto L_500; if (Idat[1] >= -2) { Idat[KSTEPX] = Idat[1]; Idat[1] = Idat[Idat[LINT]]; if (Idat[NGLO] > 0) { if (Idat[NGHI] > 0) Idat[1] = 25; } goto L_400; } else { Idat[2] = -Idat[1]; Idat[1] = 3; return; } } ierr = 1; goto L_5000; /* ****************** Starting or Restarting an integration ************* * */ L_500: Idat[1] = 1; Idat[2] = 1; Idat[LWHYF] = 16; /* Negative IDAT(NTGS) so G-Stops know we are starting. */ Idat[NTGS] = -labs( Idat[NTGS] ); Idat[NGLO] = 0; Dat[LINTX] = .1e0; /*++ Code for STAT is inactive * ESUMSQ = 0.D0 * EMAXER = 0.D0 *++ End */ goto L_220; /* *************************** Computed derivatives ********************* * */ L_1000: ; /* 2 3 4 5 6 7 8 9 10 11 12 */ switch (Idat[LWHYF] - 1) { case 1: goto L_1200; case 2: goto L_1300; case 3: goto L_1400; case 4: goto L_1500; case 5: goto L_1600; case 6: goto L_1700; case 7: goto L_1800; case 8: goto L_1900; case 9: goto L_2000; case 10: goto L_2100; case 11: goto L_2200; case 12: goto L_2300; case 13: goto L_2400; case 14: goto L_2500; case 15: goto L_2600; case 16: goto L_2700; case 17: goto L_3580; case 18: goto L_2850; } /* 13 14 15 16 17 18, 19 * * ******************* Compute Y for the Various Stages ***************** * * Stage 2 */ L_1200: for (i = 1; i <= neq; i++) { Y[i] = Dat[LOCTY + i] + h*(A31*Dat[locbf + i] + A32*Work[i]); } Ts[1] = Dat[LOCTY] + C3*h; goto L_200; /* Stage 3 */ L_1300: for (i = 1; i <= neq; i++) { Y[i] = Dat[LOCTY + i] + h*(A41*Dat[locbf + i] + A43*Work[i + neq]); } Ts[1] = Dat[LOCTY] + C4*h; goto L_200; /* Stage 4 */ L_1400: for (i = 1; i <= neq; i++) { Y[i] = Dat[LOCTY + i] + h*(A51*Dat[locbf + i] + A53*Work[i + neq] + A54*Work[i + 2*neq]); } Ts[1] = Dat[LOCTY] + C5*h; goto L_200; /* Stage 5 */ L_1500: for (i = 1; i <= neq; i++) { Y[i] = Dat[LOCTY + i] + h*(A61*Dat[locbf + i] + A64*Work[i + 2*neq] + A65*Work[i + 3*neq]); } Ts[1] = Dat[LOCTY] + C6*h; goto L_200; /* Stage 6 */ L_1600: for (i = 1; i <= neq; i++) { Y[i] = Dat[LOCTY + i] + h*(A71*Dat[locbf + i] + A74*Work[i + 2*neq] + A75*Work[i + 3*neq] + A76*Work[i + 4*neq]); } Ts[1] = Dat[LOCTY] + C7*h; goto L_200; /* Stage 7 */ L_1700: for (i = 1; i <= neq; i++) { Y[i] = Dat[LOCTY + i] + h*(A81*Dat[locbf + i] + A84*Work[i + 2*neq] + A85*Work[i + 3*neq] + A86*Work[i + 4*neq] + A87*Work[i + 5*neq]); } Ts[1] = Dat[LOCTY] + C8*h; goto L_200; /* Stage 8 */ L_1800: for (i = 1; i <= neq; i++) { Y[i] = Dat[LOCTY + i] + h*(A91*Dat[locbf + i] + A94*Work[i + 2*neq] + A95*Work[i + 3*neq] + A96*Work[i + 4*neq] + A97*Work[i + 5*neq] + A98*Work[i + 6*neq]); } Ts[1] = Dat[LOCTY] + C9*h; goto L_200; /* Stage 9 */ L_1900: for (i = 1; i <= neq; i++) { Y[i] = Dat[LOCTY + i] + h*(A101*Dat[locbf + i] + A104*Work[i + 2*neq] + A105*Work[i + 3*neq] + A106*Work[i + 4*neq] + A107*Work[i + 5*neq] + A108*Work[i + 6*neq] + A109*Work[i + 7*neq]); } Ts[1] = Dat[LOCTY] + C10*h; goto L_200; /* Stage 10 */ L_2000: for (i = 1; i <= neq; i++) { Y[i] = Dat[LOCTY + i] + h*(A111*Dat[locbf + i] + A114*Work[i + 2*neq] + A115*Work[i + 3*neq] + A116*Work[i + 4*neq] + A117*Work[i + 5*neq] + A118*Work[i + 6*neq] + A119*Work[i + 7*neq] + A1110*Work[i + 8*neq]); } Ts[1] = Dat[LOCTY] + C11*h; Idat[2] = 1; goto L_210; /* Stage 11 */ L_2100: ; /* print '(''Start stage 11'')' */ for (i = 1; i <= neq; i++) { Y[i] = Dat[LOCTY + i] + h*(A121*Dat[locbf + i] + A124*Work[i + 2*neq] + A125*Work[i + 3*neq] + A126*Work[i + 4*neq] + A127*Work[i + 5*neq] + A128*Work[i + 6*neq] + A129*Work[i + 7*neq] + A1210*Work[i + 8*neq] + A1211*Work[i]); } /* print '(''End stage 11'')' */ Ts[1] = Dat[LOCTY] + h; goto L_200; /* ******* Step Almost Finished -- Check Errors & Adjust Stepsize ******* * */ L_2200: ; /* Set new value for Y, correction from start of step, and the * error estimate. */ for (i = 1; i <= neq; i++) { /* The correction for the final Y / h. */ Work[i + 2*neq] = B1*Dat[locbf + i] + B6*Work[i + 4*neq] + B7*Work[i + 5*neq] + B8*Work[i + 6*neq] + B9*Work[i + 7*neq] + B10*Work[i + 8*neq] + B11*Work[i] + B12*Work[i + neq]; /* The error estimate of order 3. */ Work[i + 9*neq] = Work[i + 2*neq] - BHH1*Dat[locbf + i] - BHH2*Work[7*neq + i] - BHH3*Work[i + neq]; /* The error estimate of order 5. */ Work[i + 3*neq] = ER1*Dat[locbf + i] + ER6*Work[i + 4*neq] + ER7*Work[i + 5*neq] + ER8*Work[i + 6*neq] + ER9*Work[i + 7*neq] + ER10*Work[i + 8*neq] + ER11*Work[i] + ER12*Work[i + neq]; } goto L_3260; /* *********** Processing derivatives required for interpolation ******** * * Got the first function value required only for interpolation. */ L_2300: for (i = 1; i <= neq; i++) { Y[i] = Dat[LOCTY + i] + h*(A151*Dat[locbf + i] + A156*Work[4*neq + i] + A157*Work[5*neq + i] + A158*Work[6*neq + i] + A1511*Work[i] + A1512*Work[neq + i] + A1513*Work[3*neq + i] + A1514*Work[8*neq + i]); } Ts[1] = Dat[LOCTY] + C15*h; Idat[2] = 1; goto L_210; /* Just got second derivative required for interpolation now get third */ L_2400: for (i = 1; i <= neq; i++) { Y[i] = Dat[LOCTY + i] + h*(A161*Dat[locbf + i] + A166*Work[4*neq + i] + A167*Work[5*neq + i] + A168*Work[6*neq + i] + A169*Work[7*neq + i] + A1613*Work[3*neq + i] + A1614*Work[8*neq + i] + A1615*Work[i]); } Ts[1] = Dat[LOCTY] + C16*h; Idat[2] = 4*neq + 1; goto L_210; /* Just got the last derivative required for interpolation */ L_2500: ; /* Final preparation for interpolation. */ j = Idat[LOCINT]; Work[j - 2] = Dat[LOCTY]; Work[j - 1] = h; for (i = 1; i <= neq; i++) { Work[4*neq + j] = h*(Work[4*neq + j] + D413*Work[3*neq + i] + D414*Work[8*neq + i] + D415*Work[i] + D416*Work[4*neq + i]); Work[5*neq + j] = h*(Work[5*neq + j] + D513*Work[3*neq + i] + D514*Work[8*neq + i] + D515*Work[i] + D516*Work[4*neq + i]); Work[6*neq + j] = h*(Work[6*neq + j] + D613*Work[3*neq + i] + D614*Work[8*neq + i] + D615*Work[i] + D616*Work[4*neq + i]); Work[7*neq + j] = h*(Work[7*neq + j] + D713*Work[3*neq + i] + D714*Work[8*neq + i] + D715*Work[i] + D716*Work[4*neq + i]); j += 1; } /* Go save state and then take care of interpolation. */ if (Idat[NGLO] != 0) { /* Interpolation was needed for a G-Stop */ Ts[1] = Dat[LOCTGS]; Dat[LINTX] = Dat[LSAVEH]/h; } else { Ts[1] = Dat[LTOUT]; } goto L_4200; /* ****** Got derivative at start, get initial stepsize if needed ******* * * Save initial t, y, and f in their base locations. */ L_2600: Dat[LOCTY] = Ts[1]; Dat[NEXTT] = Ts[1]; for (i = 1; i <= neq; i++) { Dat[LOCTY + i] = Y[i]; Dat[locbf + i] = Work[i]; } Idat[NERRSV] = -2; Idat[KEOS] = Idat[LOCEOS]; Dat[LASTE3] = 0.e0; /* Setting DAT(HA) so it has some initialized value. */ Dat[HA] = 1.e0; /* Get square of scaled norms of F and Y in DAT(LERR), and DAT(LERR+1) * (Most of the Y arguments do no contribute to the desired results.) */ dxrk8n( idat, dat, work, y, y, y, &Dat[LERR] ); if ((Dat[LEPS]*Dat[LERR])*Dat[LEPS] >= SQ(Dat[LERTRY + 1])) { /* Error request too small. Save data to enable reset of error allowed. */ Dat[LPHI] = (Dat[LEPS]*Dat[LERR + 1])*Dat[LEPS]; ierr = 2; Idat[3] = 0; goto L_5000; } /* If have stepsize, get next output point and then to end of step stuff. * Following "go to" so that code can convert to a case. */ goto L_3600; /* Just got an extra derivative for use in getting initial stepsize. */ L_2700: for (i = 1; i <= neq; i++) { /* F at new point - F at initial point */ Work[i + neq] = Work[i] - Dat[locbf + i]; } /* Get the error scaled norm squared of the differences in F, Y (/H**2) */ dxrk8n( idat, dat, &Work[neq + 1], &Dat[locbf + 1], &Dat[LOCTY + 1], &Dat[LOCTY + 1], en ); if (En[1] == 0.e0) { /* No change in f, try increasing stepsize by 100. */ Dat[HA] *= 100.e0; } else { if (SQ(Dat[LERR + 1]) > Dat[LERR]) { /* TP2 is sqrt(y/y'') = sqrt(y/(df/h)) = h sqrt(y / df) */ tp2 = Dat[HA]*sqrt( Dat[LERR + 1]/En[1] ); } else { /* TP2 is y'/y'' = f / (df/h) = hf/df */ tp2 = Dat[HA]*sqrt( Dat[LERR]/En[1] ); } } /* 5 picked as it works reasonably well */ tp2 *= 5.e0*pow(En[1]*Dat[HA],-.125e0); /* Reduce if concern for stiffness */ if (En[2] != 0.e0) tp2 = fmin( tp2, Dat[HA]*sqrt( En[2]/En[1] ) ); Dat[HA] = fmin( fmax( tp2, Dat[LHMIN] ), Dat[HMAX] ); Ts[2] = Dat[HA]*Dat[NEGPOS]; Dat[LERR] = 0.e0; /* Restore data as if finished a step. */ Ts[1] = Dat[LOCTY]; for (i = 1; i <= neq; i++) { Y[i] = Dat[LOCTY + i]; } ah = -1.e2; /* Use below when we have HMAX small and want a plot * open (FILE='plot.out', UNIT=31) */ goto L_3650; /* End of getting initial stepsize * End of cases for computing derivatives */ /* G sign change, Try a small step to almost reach the sign change. */ L_2850: tp = fabs( Dat[NEXTH] ); h = Ts[2]; if (fabs( tp - h ) < Dat[LEPS]*10.e0*fabs( tp )) { /* Things are messed up here. We give an error message and ignore the * sign change in g. */ Ts[Idat[NTGS] - Idat[NGHI] - 2] = Ts[-Idat[NGHI]]; h = tp; ierr = 9; Daterr[1] = Ts[1]; Idat[1] = -Idat[NGHI]; goto L_5000; } Dat[LSAVEH] = Dat[NEXTH]; Dat[LOCTGS] = Ts[1]; Dat[LOCTXG] = Dat[LOCTGS]; goto L_4140; /* ************* Just processed G-Stops * * * ******** Attempt to continue an integration that has finished. ******* * */ L_3000: if (Idat[2] == 2) goto L_4100; ierr = 4; goto L_5000; /* ***************** An error condition has been reported *************** * */ L_3100: h = Ts[2]; if (Idat[2] == 2) { /* Error estimate too large at min. stepsize */ if ((Dat[LERTRY + 1] > 0) || (Dat[HA] > Dat[LHMIN])) Idat[KHMIN] = 0; Dat[LERTRY + 1] = fabs( Dat[LERTRY + 1] ); goto L_3280; } else if (Idat[2] == 3) { /* Too many steps, just continue unless it would loop */ if (Idat[KSTEPM] != Idat[KSTEP]) goto L_4000; } else if (Idat[2] == 4) { /* Just gave stiff or too much precision warning */ Idat[1] = 1; Idat[2] = 1; goto L_3320; } else if (Idat[2] == 5) { /* Just gave diagnostic for too much precision. */ tp1 = Dat[LERR + 1]; Idat[1] = 1; Idat[2] = 1; if (Idat[3] == 0) { /* Change to loosen up the error requested. */ Daterr[1] = 1.e0/sqrt( Dat[LERTRY] ); tp = Dat[LERTRY]/SQ(Dat[LERTRY + 1]); Dat[LERTRY] *= SQ(Dat[LERTRY + 1])/Dat[LPHI]; tp1 /= Dat[LPHI]; Dat[LERR] /= Dat[LPHI]; Daterr[2] = 1.e0/sqrt( Dat[LERTRY] ); dmess( mact3, (char*)mtxtab,99, idat, daterr ); } goto L_3600; } ierr = 8; goto L_5000; /* End of cases based on input value of IDAT(1) */ /* Set up to do interpolation for a G-Stop */ L_3150: Dat[LOCTGS] = Ts[1]; /* ********************** Set up to do interpolation ******************** * */ L_3200: h = Ts[2]; j = Idat[LOCINT]; Work[j - 3] = Ts[1]; Idat[NUMINT] = neq; j = Idat[LOCINT]; for (i = 1; i <= neq; i++) { Work[j] = Dat[LOCTY + i]; Work[neq + j] = Work[2*neq + i] - Work[j]; Work[2*neq + j] = h*Dat[locbf + i] - Work[neq + j]; Work[3*neq + j] = Work[neq + j] - h*Work[3*neq + i] - Work[2*neq + j]; Work[4*neq + j] = D41*Dat[locbf + i] + D46*Work[4*neq + i] + D47*Work[5*neq + i] + D48*Work[6*neq + i] + D49*Work[7*neq + i] + D410*Work[8*neq + i] + D411*Work[i] + D412*Work[neq + i]; Work[5*neq + j] = D51*Dat[locbf + i] + D56*Work[4*neq + i] + D57*Work[5*neq + i] + D58*Work[6*neq + i] + D59*Work[7*neq + i] + D510*Work[8*neq + i] + D511*Work[i] + D512*Work[neq + i]; Work[6*neq + j] = D61*Dat[locbf + i] + D66*Work[4*neq + i] + D67*Work[5*neq + i] + D68*Work[6*neq + i] + D69*Work[7*neq + i] + D610*Work[8*neq + i] + D611*Work[i] + D612*Work[neq + i]; Work[7*neq + j] = D71*Dat[locbf + i] + D76*Work[4*neq + i] + D77*Work[5*neq + i] + D78*Work[6*neq + i] + D79*Work[7*neq + i] + D710*Work[8*neq + i] + D711*Work[i] + D712*Work[neq + i]; j += 1; } /* Set up for first function value required for interpolation. */ Ts[1] = Dat[LOCTY] + C14*h; for (i = 1; i <= neq; i++) { Y[i] = Dat[LOCTY + i] + h*(A141*Dat[locbf + i] + A147*Work[5*neq + i] + A148*Work[6*neq + i] + A149*Work[7*neq + i] + A1410*Work[8*neq + i] + A1411*Work[i] + A1412*Work[neq + i] + A1413*Work[3*neq + i]); } Idat[1] = 1; Idat[2] = 8*neq + 1; Idat[LWHYF] = 13; goto L_220; /* ************* Error Estimation and Selecting Stepsize **************** * * Skip calculation of error estimate if no error control. */ L_3260: if (Dat[LERTRY] == 0.e0) goto L_3400; dxrk8n( idat, dat, &Work[3*neq + 1], &Work[9*neq + 1], y, &Dat[LOCTY + 1], en ); /* Convert order 10 estimator, ERR, to an order 14 estimator, new ERR. */ En[2] = 1.e-2*(En[2] + 1.e-3*Dat[LASTE3]); err = Dat[HA]*En[1]/sqrt( En[2] ); Dat[PHI] = log( err ) - P*log( Dat[HA] ); if (err < Dat[LERTRY + 1]) goto L_3300; /* Error is too big -- must repeat step. */ Idat[KFAIL] += 1; /*++ Code for STAT is inactive * if (TS(1) .ne. TS(2)) then * EMAXER = max(EMAXER, ERR) * end if *++ End */ if (Idat[NERRSV] != 0) { tp1 = Dat[PHI]; Dat[HLAST] = 0.e0; } else { Idat[NERRSV] = 1; tp1 = .25e0*Dat[LPHI] + .75e0*Dat[PHI]; } Dat[HA] = exp( -PINV*tp1 ); Dat[HMAX] = fmax( Dat[HLAST], Dat[HA] ); if (Dat[HA] <= Dat[LHMIN]) { Idat[KHMIN] += 1; if (Idat[KHMIN] == 1) { ierr = 5; Daterr[1] = Ts[1]; Daterr[2] = err; Daterr[3] = Ts[2]; Dat[LERTRY + 1] = -Dat[LERTRY + 1]; goto L_5000; } } L_3280: ; if (Idat[LOCEPR] != 0) { Daterr[1] = Ts[1]; Daterr[2] = Ts[2]; Daterr[3] = err; Daterr[4] = h*sqrt( En[2] ); Daterr[5] = h*sqrt( En[1] ); dmess( mact1, (char*)mtxtac,91, idat, daterr ); } Dat[NEXTH] = Dat[NEGPOS]*Dat[HA]; Ts[2] = Dat[NEXTH]; h = Ts[2]; goto L_4140; /* *************************** Select stepsize ************************** * */ L_3300: ; Dat[HLAST] = Dat[HA]; L_3320: ; /* Usual case for adjusting stepsize. */ k = 1; /*++ Code for STAT is inactive * if (TS(1) .ne. TS(3)) then * ESUMSQ = ESUMSQ + (ERR - 1.D0)**2 * EMAXER = max(EMAXER, ERR) * end if *++ End */ Dat[LASTE3] = En[2]; if (Idat[NERRSV] != 0) { if (Idat[NERRSV] == -2) { Idat[NERRSV] = -1; /* After "first" step, just a simple rule for stepsize selection */ Dat[HA] *= pow(err,-PINV); Dat[LAH] = Dat[PHI]; Idat[KSTIFF] = 0; goto L_3400; } /* Get DAT(LERR), DAT(LERR+1), DAT(LERR+2) initialized. */ Dat[LERR] = WI1A*Dat[LPHI] + WI1B*Dat[PHI]; Dat[LERR + 1] = WI2A*Dat[LPHI] + WI2B*Dat[PHI]; /* DAT(LERR+2) = WI3A * DAT(LPHI) + WI3B * DAT(PHI) */ } else { Dat[LERR] = Dat[PHI] + W*Dat[LERR]; Dat[LERR + 1] = Dat[LERR] + W*Dat[LERR + 1]; /* DAT(LERR+2) = DAT(LERR+1) + W * DAT(LERR+2) */ } /* AQ = WQ1A*DAT(LERR) + WQ2A*DAT(LERR+1) + WQ3A*DAT(LERR+2) * BQ = WQ1B*DAT(LERR) + WQ2B*DAT(LERR+1) + WQ3B*DAT(LERR+2) * CQ = WQ1C*DAT(LERR) + WQ2C*DAT(LERR+1) + WQ3C*DAT(LERR+2) */ /* AL = WL1A * DAT(LERR) + WL2A * DAT(LERR+1) * BL = WL1B * DAT(LERR) + WL2B * DAT(LERR+1) */ /* AQ = max(AQ, AL) */ ah = WL1A*Dat[LERR] + WL2A*Dat[LERR + 1]; Idat[KSTIFF] -= 1; if (En[1] > En[2]) { tp = Dat[LAH] + .75e0*log( 1.e-2*En[1]/En[2] ); if (tp > ah) { if (Idat[LOCEPR] != 0) { Daterr[1] = ah; Daterr[2] = tp; Daterr[3] = Dat[LAH]; Daterr[4] = sqrt( En[1]/En[2] ); dmess( mact4, (char*)mtxtad,63, &Idat[KSTIFF], daterr ); } if ((tp > Dat[LAH]) && (err < 1.e-4)) tp = Dat[LAH]; ah = tp; Idat[KSTIFF] = max( 0, Idat[KSTIFF] + 2 ); if (Idat[KSTIFF] == 5) { Idat[KSTIFF] = 2001; /* Give the one time diagnostic */ ierr = 7; Daterr[1] = Ts[1]; Daterr[2] = Ts[2]; goto L_5000; } } } tp = exp( -PINV*ah ); Dat[LAH] = ah; if (Idat[NERRSV] == 1) { Dat[HMAX] = Dat[HA]; if (Dat[HCAU] == Dat[LHMAXU]) Dat[HCAU] = Dat[HMAX]; } Idat[NERRSV] = 0; if (tp > Dat[HMAX]) { if (Dat[HA] >= Dat[HMAX]) Dat[HMAX] = sqrt( tp*Dat[HMAX] ); Dat[HA] = Dat[HMAX]; } else { if (Dat[HMAX] < Dat[HCAU]) { Dat[HMAX] = Dat[HCAU]; } else { Dat[HCAU] = Dat[HMAX]; } Dat[HA] = fmax( Dat[LHMIN], tp ); } L_3400: ; Dat[NEXTH] = Dat[HA]*Dat[NEGPOS]; if (Idat[LOCEPR] != 0) { Daterr[1] = Ts[1]; Daterr[2] = Ts[2]; Daterr[3] = err; Daterr[4] = Dat[PHI]; Daterr[5] = sqrt( En[1]/En[2] ); Daterr[6] = Dat[HMAX]; Daterr[7] = Dat[HCAU]; /* Print info. about error estimation */ dmess( mact1, (char*)mtxtae,72, idat, daterr ); k = Idat[LOCEPR]; if (Idat[k + 2] != 0) { L_3420: i1 = Idat[k]; if (i1 != 0) { Mact2[3] = i1; Mact2[5] = -Idat[k + 1]; dmess( mact2, (char*)mtxtaf,18, &Idat[k], &Work[3*neq + 1] ); k += 2; goto L_3420; } } } Dat[LPHI] = Dat[PHI]; /* End of checking errors and adjusting stepsize * * ************************ Set new base t and new y ******************** * */ i1 = neq; i = 1; k = Idat[LOCXP]; if (k != 0) { /* Accumulate t in extra precision */ j = Idat[k]; Ts[1] = Dat[LOCTY] + (h + Dat[j]); Dat[j] = dxrk8x( Ts[1], Dat[LOCTY], h, Dat[j] ); Dat[NEXTT] = Ts[1]; i1 = Idat[k + 1]; if (i1 == 1) goto L_3550; } else { Dat[NEXTT] = Dat[LOCTY] + h; } L_3520: for (i = i; i <= i1; i++) { Work[9*neq + i] = Y[i]; /* The new Y */ Y[i] = Dat[LOCTY + i] + h*Work[i + 2*neq]; Work[i + 2*neq] = Y[i]; } L_3550: if (i <= neq) { k += 2; i1 = Idat[k]; for (i = i; i <= i1; i++) { Work[9*neq + i] = Y[i]; j += 1; /* First get the correction to Y. */ tp = h*Work[i + 2*neq]; Y[i] = Dat[LOCTY + i] + (tp + Dat[j]); Dat[j] = dxrk8x( Y[i], Dat[LOCTY + i], tp, Dat[j] ); Work[i + 2*neq] = Y[i]; } i1 = Idat[k + 1]; if (i <= neq) goto L_3520; } /* End of setting new base t and base y. * * Get final derivative value for the step */ Idat[2] = 3*neq + 1; Idat[LWHYF] = 18; goto L_220; /* ****************** Got last derivative on the step ******************* * */ L_3580: Idat[KSTEP] += 1; Idat[KEOS] = Idat[LOCEOS]; goto L_4000; /* **************** Get stepsize for starting the integration *********** */ L_3600: if (Ts[2] != 0.e0) { if (Dat[NEGPOS]*Ts[2] >= 0.e0) goto L_3650; Daterr[1] = Ts[2]; ierr = 6; goto L_5000; } /* Compute a stepsize for taking an Euler step. Use the information from * this step to generate a guess for the initial stepsize. Algorithm is * based (on possibly invalid) assumption that the derivatives of y * decrease (or increase) with a relatively constant ratio. To the * extent that this is not true, the initial stepsize will not be * selected very well. At worst this should just waste a few function * evaluations. * */ if (Dat[LERR] != 0.e0) { /* Set TP = est. (||y||_e / ||f||_e)^2 */ tp = fmax( sqrt( Dat[LERR + 1]/Dat[LERR] ), 1.e-2 ); Dat[HA] = .05e0*tp*pow(Dat[LERR + 1] + Dat[LERR],-.0625e0); } else { Dat[HA] = .05e0; tp = pow(Dat[LERR + 1],.0625e0); } /* Get DAT(HA) in range. */ Dat[HA] = fmin( Dat[HA], Dat[HMAX] ); /* Compute a new Y */ h = Dat[HA]*Dat[NEGPOS]; for (i = 1; i <= neq; i++) { Y[i] = Dat[LOCTY + i] + h*Work[i]; } /* Prepare for getting the derivative */ Dat[LERR + 1] = tp; Ts[1] += h; goto L_210; /* Just finished establishing the initial stepsize * */ L_3650: Dat[NEXTT] = Ts[1]; Dat[NEXTH] = Ts[2]; h = Ts[2]; /* Get next output point set. */ goto L_3750; /* *************** Just processed some kind of output point. ************ * * Got output point using some kind of interpolation. */ L_3700: if (Idat[1] < 25) { if (Idat[1] <= 21) goto L_3800; /* Processed some kind of T output point */ j = Idat[3]; k += 1; if (Dat[LOCSAV] == Dat[j]) { if (Idat[1] == 24) { if (Dat[NEGPOS]*Dat[j + 1] > 0.e0) { Dat[j] += Dat[j + 1]; goto L_3750; } } if (2*Idat[2] == Idat[LINT] - Idat[INTSAV]) { /* The first output we were examining is now inactive. */ Idat[LOCOUT] += 2; } else { /* Store a value well beyond TS(3) so this output is ignored later. */ Dat[j] = sign( 2.e0*fabs( Ts[3] ) + 1.e35, Dat[NEGPOS] ); } } } else if (Idat[1] == 30) { if (Ts[1] == Ts[3]) { /* The integration is finished. *++ Code for STAT is inactive * ESUMSQ = sqrt(ESUMSQ) * print '(''W ='', F6.3, '' KFAIL ='', I4, '' Norm Diff ='', * 1F7.1, '' Max. Diff ='', 1PE8.1)', * 2 W,IDAT(KFAIL),ESUMSQ,EMAXER/DAT(LERTRY+1) *c Returns value to the test driver. * DAT(24) = ESUMSQ * DAT(25) = EMAXER *++ End */ Idat[1] = 3; return; } Dat[NEXTH] = Dat[LSAVEH]; } else { if (Idat[2] <= Idat[3]) { /* We have reported all G-Stops at this time. */ if (Idat[2] <= Idat[NXGS]) { /* Must restart after an extrapolating G-Stop */ Ts[2] = Dat[LSAVEH]; goto L_500; } /* Interpolating G-Stop * Check at end of step for another G-Stop */ Idat[KEOS] = -1; Idat[NGLO] = 0; goto L_3800; } else { /* There should be more G-Stopa for this value of TS(1) */ goto L_4310; } } /* Get next DAT(LTOUT) and IDAT(LINT) */ L_3750: ; Dat[LTOUT] = Ts[3]; tp = Ts[1]; Idat[LINT] = Idat[LTFIN]; j = Idat[LOCOUT]; L_3780: if (Idat[j] != 30) { i1 = Idat[j + 1]; if (Dat[NEGPOS]*(Dat[i1] - Dat[LTOUT]) < 0.e0) { /* This output point precedes the current pick. */ Dat[LTOUT] = Dat[i1]; Idat[LINT] = j; } j += 2; goto L_3780; } if (Idat[NGLO] > 0) { /* An interpolating G-Stop is in the queue */ if (Idat[KSTEPX] == -1) goto L_500; Ts[1] = Dat[LOCTGS]; goto L_4200; } /* Check for user interupt after taking care of interpolation updates. */ L_3800: ; /* Check if we want a restart */ if (Idat[KSTEPX] == -1) goto L_500; /* ********************** Process end of step actions ***************** * * IDAT(KEOS) = -3 End of step then G-Stop * = -2 End of step then TOUT * = -1 G-Stop then TOUT * = 0 TOUT then finish * = 1 Finish */ L_4000: k = Idat[KEOS]; if (k < 0) { if (k < -1) { /* End of step (on every step) return */ Idat[1] = 21; Idat[2] = 1; Idat[KEOS] += 2; goto L_4250; } else { /* Process G-Stops */ Idat[KEOS] = 0; if (Idat[1] < 25) Work[Idat[LOCINT] - 4] = Dat[LOCTY]; Idat[1] = 2; Idat[3] = 0; if (Idat[NGHI] < 0) { /* We have taken a smaller step for an extrapolating G-Stop. Now * interpolate (actually extrapolate) to where we think the sign changes. */ Ts[1] = Dat[LOCTGS]; Ts[2] = h; goto L_3200; } /* Update solution to the end of the step if it is not already there. * (Needed if processing after got a G-Stop on the step already.) */ if (Ts[1] != Dat[NEXTT]) { if (Dat[NEGPOS]*(Dat[NEXTT] - Ts[1]) < 0.e0) { /* We extrapolated for a G-Stop */ Ts[1] = Dat[LOCTXG]; if (Dat[NEGPOS]*(Ts[1] - Ts[3]) > 0.e0) Ts[1] = Ts[3]; Idat[NGHI] = -Idat[NGHI]; Idat[NGLO] = Idat[NGHI]; goto L_4200; } /* Restore solution to that at the end of the step to check again. */ Ts[1] = Dat[NEXTT]; for (i = 1; i <= neq; i++) { Y[i] = Work[2*neq + i]; } } j = Idat[LOCINT]; Work[j - 2] = Dat[LOCTY]; Work[j - 1] = h; goto L_4300; } } /* Check for output point */ if (Dat[NEGPOS]*(Dat[NEXTT] - Dat[LTOUT]) >= 0.e0) { /* End of step passes an output point. * Check if reached the final point. */ Idat[1] = Idat[Idat[LINT]]; if (Idat[1] == 30) { Idat[2] = 0; if (Ts[1] != Ts[3]) { /* Restore solution to that at the end of the step to check again. */ Ts[1] = Ts[3]; for (i = 1; i <= neq; i++) { Y[i] = Work[2*neq + i]; } } goto L_4250; } if (Idat[1] == 23) { Idat[Idat[LINT] + 1] = Dat[NEXTT]; Dat[LTOUT] = Dat[NEXTT]; } /* Ready for output if same t. */ if (Ts[1] == Dat[LTOUT]) goto L_4220; /* Get ready for interpolation if not ready yet. */ if (Idat[NUMINT] < 0) goto L_3200; /* Do the interpolation */ Ts[1] = Dat[LTOUT]; goto L_4200; } Idat[KEOS] = 1; if (Dat[LOCTY] != Dat[NEXTT]) { /* Save information required to take the next step. * Set new base t, y, and f. */ Dat[LOCTY] = Dat[NEXTT]; for (i = 1; i <= neq; i++) { Dat[LOCTY + i] = Work[2*neq + i]; Dat[locbf + i] = Work[3*neq + i]; } } /* Essentially finished with the step. */ if (Idat[KSTEP] >= Idat[KSTEPX]) { if (Idat[KSTEPX] < 0) { /* User wants an interupt */ Idat[2] = -Idat[KSTEPX]; Idat[KSTEPX] = Idat[KSTEPM]; Idat[1] = 3; return; } /* Step count has reached point where something special is called for. */ if (Idat[KSTEPM] <= Idat[KSTEP]) { /* Hit maximum number of steps */ Idat[KSTEPM] += Idat[KSTEPM + 1]; Idat[KSTEPX] = Idat[KSTEPM]; Idat[1] = 20; goto L_4250; } } /* Standard value for next stepsize and continuation entry. */ L_4100: ; Ts[2] = Dat[NEXTH]; if (Dat[NEGPOS]*(Ts[3] - Dat[LOCTY] - Ts[2]) < 1.e-2*fabs( Ts[2] )) { Ts[2] = Dat[NEGPOS]*(Ts[3] - Dat[LOCTY]); Dat[HA] = fabs( Ts[2] ); if (Dat[HA] <= Dat[LEPS]*fabs( Dat[LOCTY] )) { /* Got to the final point */ Idat[1] = 3; Idat[2] = 0; Ts[1] = Ts[3]; goto L_4250; } } /* Save standard stepsize */ Dat[LSAVEH] = Dat[NEXTH]; h = Ts[2]; /* Start the next step. */ L_4140: Idat[LWHYF] = 2; Idat[NUMINT] = -neq; Idat[1] = 1; Idat[2] = 1; for (i = 1; i <= neq; i++) { Y[i] = Dat[LOCTY + i] + h*A21*Dat[locbf + i]; } Ts[1] = Dat[LOCTY] + C2*h; goto L_220; /* Do the interpolation. */ L_4200: Work[1] = Dat[LINTX]; dxrk8i( Ts[1], y, idat, work ); L_4220: if (Idat[NGLO] != 0) { if (Idat[NGLO] < 0) { /* Interpolated for a G-Stop */ Idat[1] = 2; Idat[2] = -1; Idat[3] = Idat[NGHI]; /* If interpolating for an extrapolated stop, we must start fresh. */ if (Idat[3] < 0) Idat[2] = 0; goto L_4310; } else { /* G-Stop has been found */ if (Idat[NGHI] < 0) { /* Another t output precedes the G-Stop */ Idat[NGHI] = -Idat[NGHI]; } else { Idat[2] = Idat[NGHI]; Idat[3] = Idat[NGLO]; Idat[1] = 26; if (Idat[2] <= Idat[NXGS]) { Idat[1] = 25; Dat[NEXTH] = Dat[LSAVEH]; Ts[2] = Dat[NEXTH]; } goto L_4250; } } } /* Set flags for T output points */ k = Idat[LINT]; Idat[1] = Idat[k]; Idat[2] = (k + 2 - Idat[INTSAV])/2; Dat[LOCSAV] = Ts[1]; Idat[3] = Idat[k + 1]; /* Output to user. */ L_4250: if (Idat[LREVO] != 0) return; dxrk8o( ts, y, idat, dat ); goto L_400; /* Call DXRK8F for a G-Stop */ L_4300: Idat[2] = 0; Idat[3] = 0; L_4310: Idat[LWHYF] = 1; if (Idat[LREVF] != 0) return; dxrk8f( &Ts[1], y, work, idat ); L_4320: ; if ((Idat[LWHYF] != 21) || (Idat[2] == -1)) { Daterr[1] = Ts[1]; ierr = 3; goto L_5000; } /* If IDAT(2) is 0, there was no sign change and no G-Stop */ if (Idat[2] == 0) goto L_4000; /* If IDAT(2) is -2 we need to interpolate for a G-stop. */ if (Idat[2] == -2) goto L_3150; /* Got a G-Stop */ Dat[LOCTGS] = Ts[1]; /* If G-Stop precedes next output, then output the G-Stop */ if (Dat[NEGPOS]*(Ts[1] - Dat[LTOUT]) <= 0.e0) goto L_4220; /* Output point precedes the G-Stop, save current T, solve to TOUT. */ Idat[NGHI] = -Idat[NGHI]; if (Idat[Idat[LINT]] == 30) { /* Integrate to the end point */ h = Ts[3] - Dat[NEXTT]; Ts[2] = h; goto L_4140; } else { Ts[1] = Dat[LTOUT]; } goto L_4200; /* Some kind of error. */ L_5000: Mact[2] = Istop[ierr]; Mact[3] = Ierror[ierr]; Mact[4] = Lermsg[ierr]; dmess( mact, (char*)mtxtaa,203, idat, daterr ); if (ierr == 9) goto L_4140; Idat[1] = 4; Idat[2] = Mact[3]; return; } /* end of function */ void /*FUNCTION*/ dxrk8i( double t, double y[], long idat[], double work[]) { long int i, j, niy; double daterr[2], s, s1; static char mtxtaa[1][79]={"DXRK8$BInterpolating outside of safe interval by $F of the stepsize at T=$F.$E"}; static char mtxtab[1][65]={"DXRK8$BA bug in dxrk8, dxrk8i called when not allowed at T=$F.$E"}; static long mact[5]={MEEMES,24,28,0,MERET}; /* OFFSET Vectors w/subscript range: 1 to dimension */ double *const Daterr = &daterr[0] - 1; long *const Idat = &idat[0] - 1; long *const Mact = &mact[0] - 1; double *const Work = &work[0] - 1; double *const Y = &y[0] - 1; /* end of OFFSET VECTORS */ /* Interpolate to get y(t), using data structures in IDAT and WORK. * Only certain components of y are set as defined in IDAT. * * Formal Arguments */ /* Local variables */ /* Parameters used to reference IDAT */ /* Local variables */ /* Parameters for error messages */ /* Other Stuff for error processing */ /* ********* Error message text *************** *[Last 2 letters of Param. name] [Text generating message.] *AA DXRK8$B *AB Interpolating outside of safe interval by $F of the stepsize $C * at T=$F.$E * $ *AC DXRK8$B *AD A bug in dxrk8, dxrk8i called when not allowed at T=$F.$E */ /* **** End of automatically generated text * */ /* Check if ready. */ niy = Idat[NUMINT]; if (niy <= 0) { Mact[2] = 99; Mact[3] = 29; Daterr[1] = t; dmess( mact, (char*)mtxtab,65, idat, daterr ); exit(0); } j = Idat[LOCINT]; s = (t - Work[j - 2])/Work[j - 1]; if ((s > Work[1] + 1.e0) || (s < -Work[1])) { /* Error interpolation outside allowed interval. */ Daterr[1] = s; Daterr[2] = t; dmess( mact, (char*)mtxtaa,79, idat, daterr ); } s1 = 1.e0 - s; for (i = 1; i <= niy; i++) { Y[i] = Work[j] + s*(Work[niy + j] + s1*(Work[2*niy + j] + s*(Work[3*niy + j] + s1*(Work[4*niy + j] + s*(Work[5*niy + j] + s1*(Work[6*niy + j] + s*Work[7*niy + j])))))); j += 1; } return; } /* end of function */ void /*FUNCTION*/ dxrk8n( long idat[], double dat[], double v1[], double v2[], double y[], double yold[], double en[]) { long int i, i1, i2, j, k, l; double pv1, pv2, sv1, sv2, sy, tp, tp1, tp2; /* OFFSET Vectors w/subscript range: 1 to dimension */ double *const Dat = &dat[0] - 1; double *const En = &en[0] - 1; long *const Idat = &idat[0] - 1; double *const V1 = &v1[0] - 1; double *const V2 = &v2[0] - 1; double *const Y = &y[0] - 1; double *const Yold = &yold[0] - 1; /* end of OFFSET VECTORS */ /* Computes squared norms of V1, V2, Y and F as well as Y^T F scaled by * error tolerances. * * Formal arguments */ /* Parameters used to reference DAT */ /* Parameters used to reference IDAT */ /* Local variables */ j = LOCERR; i2 = 0; sv1 = 8.e-32; sv2 = 1.e-32; goto L_20; L_10: ; /* Updates when error request is a scalar. */ sv1 += pv1*tp; sv2 += pv2*tp; L_20: if (i2 == Idat[4]) { En[1] = sv1*Dat[LERTRY]; En[2] = sv2*Dat[LERTRY]; return; } k = Idat[j]; i1 = i2 + 1; i2 = Idat[j + 1]; l = Idat[j + 2]; j += 3; if (k <= 3) { pv1 = 0.e0; pv2 = 0.e0; if (k == 2) { /* Simple absolute error test */ for (i = i1; i <= i2; i++) { pv1 += SQ(V1[i]); pv2 += SQ(V2[i]); } tp = 1.e0/SQ(Dat[l]); } else { /* Mixed absolute / relative error test */ sy = 0.e0; for (i = i1; i <= i2; i++) { pv1 += SQ(V1[i]); pv2 += SQ(V2[i]); sy += SQ(Y[i]); } tp = 1.e0/(SQ(Dat[l]) + sy*SQ(Dat[l + 1])); } goto L_10; } else if (k == 4) { /* Separate for each Y -- mixed absolute / relative error test. */ tp1 = SQ(Dat[l]); tp2 = SQ(Dat[l + 1]); for (i = i1; i <= i2; i++) { tp = 1.e0/(tp1 + tp2*(SQ(Y[i]) + SQ(Yold[i]))); sv1 += tp*SQ(V1[i]); sv2 += tp*SQ(V2[i]); } } else if (k == 5) { /* Vector absolute error test */ for (i = i1; i <= i2; i++) { tp = 1.e0/Dat[i + l]; sv1 += powi(V1[i]*tp,2); sv2 += powi(V2[i]*tp,2); } } else if (k == 6) { /* Vector mixed absolute / relative error test */ for (i = i1; i <= i2; i++) { tp = 1.e0/(SQ(Dat[l + 2*i]) + SQ(Dat[l + 2*i + 1])*(SQ(Y[i]) + SQ(Yold[i]))); sv1 += tp*SQ(V1[i]); sv1 += tp*SQ(V2[i]); } } goto L_20; /* End of subroutine DXRK8N */ } /* end of function */ double /*FUNCTION*/ dxrk8x( double vnew, double vold, double vdiff, double vlow) { double dxrk8x_v; /* Get new least significant part when carrying extra precision. * Code could be inline when floating point registers store numbers * to the same precision as memory or if one knows the compiler * stores intermediate results to memory. */ /* With perfect arithmetic this would return 0, if because of round-off * VNEW is a little too small, then DXRK8X > 0 which will give a little * boost to the Y on the next step. */ dxrk8x_v = (vold - vnew) + vdiff + vlow; return( dxrk8x_v ); } /* end of function */