/*Translated by FOR_C, v3.4.2 (-), on 07/09/115 at 08:31:44 */ /*FOR_C Options SET: ftn=u io=c no=p op=aimnv s=dbov str=l x=f - prototypes */ #include #include "fcrt.h" #include "dxrk8g.h" #include #include /* PARAMETER translations */ #define KDZERO 21 #define LOCINT 38 #define LTXTAB 9 #define LTXTAC 103 #define LWHYF 40 #define MEEMES 52 #define MERET 51 #define NGHI 32 #define NGLO 31 #define NTGS 17 #define NUMINT 39 #define NXGS 16 /* end of PARAMETER translations */ void /*FUNCTION*/ dxrk8g( double ts[], double y[], double work[], long idat[]) { long int i, lg, li, ng; double daterr[1], tl; static char mtxtaa[2][144]={"DXRK8G$BCalling DXRK8G when IDAT(1) is not 1 or 2, and no G-Stop was found for IDAT(1) = $I, T = $F.$ETS($I) has changed sign in an interval wh", "ere previous checks showed no sign change. The location of the stop returned, $F, is at best questionable, and is certainly of low accuracy.$E"}; static long mact[5]={MEEMES,36,0,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 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. *>> 2008-02-25 DXRK8G Krogh Lots of changes in usage, bug fixes. *>> 1997-02-18 DXRK8G Krogh Initial code. *--D replaces "?": ?XRK8G, ?XRK8I, ?XRK8, ?MESS, ?ZERO */ /* Checks and finds G-Stops for the integrator DXRK8. * From dxrk8f: call dxrk8g(TS, Y, F, IDAT) * Reverse Comm: call dxrk8g(TS, Y, WORK, IDAT) * * User code from dxrk8f looks something like: * if (IDAT(1) .eq. 1) then * if (There are extrapolating G-Stops) then * TS(3+i) = value of the i-th extrapolating G_i, i = 1, ... * call DXRKG(TS, Y, F, IDAT) * if (IDAT(3) .ne. 0) return * end if * Compute F and return * end if * do * Compute TS(3+i) for all G_i, i = 1, ... * if (IDAT(3) .ne. -1) return * end do * * * ************************ Calling Sequence Arguments ****************** * * T The current value of the independent variable. * Y The current value of the dependent variable. * WORK Contains information about the stepsize, previous T, last T for * which G's were evaluated, and when interpolating, contains other * data required for the interpolation. * IDAT Used to hold various integer values used in the integration. * * ************************* Variable Definitions *********************** * * Most names are defined in DXRK8. Just those local to DXRK8G here. * * LI IDAT(LOCINT) = Location in WORK where information connected with * interpolation is saved. See IDAT(LOCINT) in DXRK8. * LG Difference between indexes for g and the old g. * NG Index for current G. * TL Value of T where G's were last computed on call to DZERO. * * ******************** Variable Declarations *************************** * */ /* 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 DXRK8G$B *AB Calling DXRK8G when IDAT(1) is not 1 or 2, and no G-Stop was $C * found for IDAT(1) = $I, T = $F.$E *AC TS($I) has changed sign in an interval where previous checks $C * showed no sign change. The location of the stop returned, $F, $C * is at best questionable, and is certainly of low accuracy.$E */ /* **** End of automatically generated text * */ /* ********************** Start of Executable Code ********************** * */ lg = labs( Idat[NTGS] ) - 2; if (Idat[1] != 2) { if (Idat[1] == 1) { Idat[3] = 0; if (Idat[NTGS] < 0) { /* Just started an integraion, just Copy G to GOLD */ for (ng = 4; ng <= Idat[NXGS]; ng++) { Ts[ng + lg] = Ts[ng]; } /* Remember t at the saved values for G/ */ Ts[lg + 3] = Ts[1]; return; } /* Checking for sign changes. */ for (ng = 4; ng <= Idat[NXGS]; ng++) { /* Check extrapolating G's when computing F. */ if (Ts[ng]*Ts[ng + lg] <= 0) { if (Ts[ng + lg] != 0.e0) { if (Idat[LWHYF] >= 13) { /* Sign change occurred while getting derivatives for interpolation. */ Mact[3] = 25; Mact[4] = LTXTAC; Idat[1] = 2; Idat[2] = ng; Idat[3] = ng; i = 2; goto L_510; } /* Take a step that doesn't quite reach this sign change. */ Ts[2] = .9e0*(Ts[1] - Ts[lg + 3])*fabs( Ts[lg + ng]/ (Ts[lg + ng] - Ts[ng]) ); Ts[1] = Ts[lg + 3] + 1.8*Ts[2]; Idat[LWHYF] = 19; Idat[3] = ng; Idat[NGLO] = -ng; Idat[NGHI] = -ng; return; } } } /* Found no stops */ return; } /* Set IDAT(LWHYF) so we know G's were computed. */ Idat[LWHYF] = 21; if (((Idat[3] >= 4) && (Idat[3] < Idat[2])) && (Idat[2] < lg + 3)) { /* Apparently there were multiple stops athe same value of t, get next. */ for (ng = Idat[2] - 1; ng >= Idat[3]; ng--) { if ((Ts[ng]*Ts[ng + lg] <= 0.e0) && (Ts[ng + lg] != 0.e0)) { Idat[2] = ng; Idat[NGHI] = ng; Ts[ng + lg] = 0.e0; return; } } } Mact[3] = 7; Mact[4] = LTXTAB; /* No good way to get error flag back to users main program * without setting IDAT(1) which could snowball the problem. */ goto L_500; } /* End of step actions -- Set IDAT(LWHYF) so we know G's were computed. */ Idat[LWHYF] = 21; if (Idat[NTGS] < 0) { /* Initial setup for TOLD and GOLD */ Idat[NTGS] = -Idat[NTGS]; for (ng = 4; ng <= (lg + 2); ng++) { Ts[ng + lg] = Ts[ng]; } Ts[ng] = Ts[1]; return; } li = Idat[LOCINT]; /* If IDAT(2) is -1, then we are in the middle of an iteration */ if (Idat[2] == -1) goto L_320; /* User just computed a full set of G's. Set flag assuming no stop. */ if (Idat[3] > 0) { ng = Idat[3] + 1; } else { ng = 4; } if (fabs( Ts[1] - Ts[lg + 3] ) <= .25e0*fabs( Ts[2] )) { /* We want to leave GOLD_i as 0 if it is 0 when close to last 0. * (This is prevent noise from resulting in multiple G-Stops.) */ for (ng = ng; ng <= (lg + 2); ng++) { if (Ts[ng + lg] != 0) { if (Ts[ng]*Ts[ng + lg] <= 0.e0) { if (Ts[ng] != 0.e0) goto L_300; } Ts[ng + lg] = Ts[ng]; } } goto L_100; } /* Continue here after finding a stop. */ L_70: ; /* t is not close to location of previously computed or just found a * sign change */ for (ng = ng; ng <= (lg + 2); ng++) { if (Ts[ng]*Ts[ng + lg] <= 0.e0) { if (Ts[ng + lg] != 0.e0) goto L_300; } Ts[ng + lg] = Ts[ng]; } L_100: if (Idat[3] == 0.e0) { /* No sign change was seen. */ Ts[ng] = Ts[1]; Idat[2] = 0; return; } ng = Idat[3]; Idat[NGHI] = ng; Idat[NGLO] = ng; Idat[2] = ng; /* Set to 0 so rounding errors don't give us an extra stop. */ Ts[ng + lg] = 0.e0; /* Check if there is a G-Stop preceding the one found. * (This can only happen if there are multiple stops at a given time.) */ for (ng = Idat[2] - 1; ng >= 4; ng--) { if (Ts[ng]*Ts[ng + lg] <= 0.e0) { if (Ts[ng + lg] != 0.e0) { /* Found a stop remember the first. */ Idat[NGLO] = ng; } else { Ts[ng + lg] = Ts[ng]; } } } /* Set last TS(1) where G-Stops were computed */ Ts[lg + 3] = Ts[1]; return; /* Start search for a G -Stop */ L_300: Idat[2] = -1; Idat[KDZERO] = 0; Idat[3] = -ng; Work[li - 6] = Ts[ng + lg]; Work[li - 5] = Ts[lg + 3]; goto L_350; /* In the middle of an iteration */ L_320: ng = Idat[3]; if (ng > 0) goto L_420; ng = -ng; /*## Note DZERO has save variables as does MESS and DMESS. * At some point we need to set aside space in WORK for all saves that * might be required. In DZERO this is a good size set: DIV, FL, FLMFB, * FO, KNKP, KS, KTYP, LCHG, LMODE, LNLP, MACT, NP, RND, SMALL, XL, * XLMXB, XO, XX, XXMXOL. (Even more in MESS and DMESS!) A simpler code * would need less. Or we could allow save variable here since the * "save" is only over the duration of the iteration for the 0??? * */ L_350: tl = Ts[1]; dzero( &Ts[1], &Ts[ng], &Work[li - 5], &Work[li - 6], &Idat[KDZERO], 0.e0 ); if (Idat[NUMINT] <= 0) { /* Not ready to interpolate */ Idat[2] = -2; Idat[3] = -ng; Idat[NGHI] = ng; Idat[NGLO] = -ng; return; } if (Idat[KDZERO] != 1) goto L_400; L_380: dxrk8i( Ts[1], y, idat, work ); return; /* Some kind of convergence */ L_400: if (Ts[2]*(Work[li - 5] - Ts[1]) > 0.e0) Ts[1] = Work[li - 5]; Idat[3] = ng; if (Ts[2]*(Ts[1] - tl) > 0) goto L_380; Ts[1] = tl; /* Get here after a final g evaluation or on multiple stops. */ L_420: Idat[2] = ng; ng += 1; goto L_70; /* Process an error. */ L_500: i = 1; L_510: Daterr[1] = Ts[1]; dmess( mact, (char*)mtxtaa,144, &Idat[i], daterr ); return; } /* end of function */