share/odepack/fortran/drchek.f

   1 *DECK DRCHEK
   2       SUBROUTINE DRCHEK (JOB, G, NEQ, Y, YH,NYH, G0, G1, GX, JROOT, IRT)
   3       EXTERNAL G
   4       INTEGER JOB, NEQ, NYH, JROOT, IRT
   5       DOUBLE PRECISION Y, YH, G0, G1, GX
   6       DIMENSION NEQ(*), Y(*), YH(NYH,*), G0(*), G1(*), GX(*), JROOT(*)
   7       INTEGER IOWND, IOWNS,
   8      1   ICF, IERPJ, IERSL, JCUR, JSTART, KFLAG, L,
   9      2   LYH, LEWT, LACOR, LSAVF, LWM, LIWM, METH, MITER,
  10      3   MAXORD, MAXCOR, MSBP, MXNCF, N, NQ, NST, NFE, NJE, NQU
  11       INTEGER IOWND3, IOWNR3, IRFND, ITASKC, NGC, NGE
  12       DOUBLE PRECISION ROWNS,
  13      1   CCMAX, EL0, H, HMIN, HMXI, HU, RC, TN, UROUND
  14       DOUBLE PRECISION ROWNR3, T0, TLAST, TOUTC
  15       COMMON /DLS001/ ROWNS(209),
  16      1   CCMAX, EL0, H, HMIN, HMXI, HU, RC, TN, UROUND,
  17      2   IOWND(6), IOWNS(6),
  18      3   ICF, IERPJ, IERSL, JCUR, JSTART, KFLAG, L,
  19      4   LYH, LEWT, LACOR, LSAVF, LWM, LIWM, METH, MITER,
  20      5   MAXORD, MAXCOR, MSBP, MXNCF, N, NQ, NST, NFE, NJE, NQU
  21       COMMON /DLSR01/ ROWNR3(2), T0, TLAST, TOUTC,
  22      1   IOWND3(3), IOWNR3(2), IRFND, ITASKC, NGC, NGE
  23       INTEGER I, IFLAG, JFLAG
  24       DOUBLE PRECISION HMING, T1, TEMP1, TEMP2, X
  25       LOGICAL ZROOT
  26 C-----------------------------------------------------------------------
  27 C This routine checks for the presence of a root in the
  28 C vicinity of the current T, in a manner depending on the
  29 C input flag JOB.  It calls Subroutine DROOTS to locate the root
  30 C as precisely as possible.
  31 C
  32 C In addition to variables described previously, DRCHEK
  33 C uses the following for communication:
  34 C JOB    = integer flag indicating type of call:
  35 C          JOB = 1 means the problem is being initialized, and DRCHEK
  36 C                  is to look for a root at or very near the initial T.
  37 C          JOB = 2 means a continuation call to the solver was just
  38 C                  made, and DRCHEK is to check for a root in the
  39 C                  relevant part of the step last taken.
  40 C          JOB = 3 means a successful step was just taken, and DRCHEK
  41 C                  is to look for a root in the interval of the step.
  42 C G0     = array of length NG, containing the value of g at T = T0.
  43 C          G0 is input for JOB .ge. 2, and output in all cases.
  44 C G1,GX  = arrays of length NG for work space.
  45 C IRT    = completion flag:
  46 C          IRT = 0  means no root was found.
  47 C          IRT = -1 means JOB = 1 and a root was found too near to T.
  48 C          IRT = 1  means a legitimate root was found (JOB = 2 or 3).
  49 C                   On return, T0 is the root location, and Y is the
  50 C                   corresponding solution vector.
  51 C T0     = value of T at one endpoint of interval of interest.  Only
  52 C          roots beyond T0 in the direction of integration are sought.
  53 C          T0 is input if JOB .ge. 2, and output in all cases.
  54 C          T0 is updated by DRCHEK, whether a root is found or not.
  55 C TLAST  = last value of T returned by the solver (input only).
  56 C TOUTC  = copy of TOUT (input only).
  57 C IRFND  = input flag showing whether the last step taken had a root.
  58 C          IRFND = 1 if it did, = 0 if not.
  59 C ITASKC = copy of ITASK (input only).
  60 C NGC    = copy of NG (input only).
  61 C-----------------------------------------------------------------------
  62       IRT = 0
  63       DO 10 I = 1,NGC
  64  10     JROOT(I) = 0
  65       HMING = (ABS(TN) + ABS(H))*UROUND*100.0D0
  66 C
  67       GO TO (100, 200, 300), JOB
  68 C
  69 C Evaluate g at initial T, and check for zero values. ------------------
  70  100  CONTINUE
  71       T0 = TN
  72       CALL G (NEQ, T0, Y, NGC, G0)
  73       NGE = 1
  74       ZROOT = .FALSE.
  75       DO 110 I = 1,NGC
  76  110    IF (ABS(G0(I)) .LE. 0.0D0) ZROOT = .TRUE.
  77       IF (.NOT. ZROOT) GO TO 190
  78 C g has a zero at T.  Look at g at T + (small increment). --------------
  79       TEMP1 = SIGN(HMING,H)
  80       T0 = T0 + TEMP1
  81       TEMP2 = TEMP1/H
  82       DO 120 I = 1,N
  83  120    Y(I) = Y(I) + TEMP2*YH(I,2)
  84       CALL G (NEQ, T0, Y, NGC, G0)
  85       NGE = NGE + 1
  86       ZROOT = .FALSE.
  87       DO 130 I = 1,NGC
  88  130    IF (ABS(G0(I)) .LE. 0.0D0) ZROOT = .TRUE.
  89       IF (.NOT. ZROOT) GO TO 190
  90 C g has a zero at T and also close to T.  Take error return. -----------
  91       IRT = -1
  92       RETURN
  93 C
  94  190  CONTINUE
  95       RETURN
  96 C
  97 C
  98  200  CONTINUE
  99       IF (IRFND .EQ. 0) GO TO 260
 100 C If a root was found on the previous step, evaluate G0 = g(T0). -------
 101       CALL DINTDY (T0, 0, YH, NYH, Y, IFLAG)
 102       CALL G (NEQ, T0, Y, NGC, G0)
 103       NGE = NGE + 1
 104       ZROOT = .FALSE.
 105       DO 210 I = 1,NGC
 106  210    IF (ABS(G0(I)) .LE. 0.0D0) ZROOT = .TRUE.
 107       IF (.NOT. ZROOT) GO TO 260
 108 C g has a zero at T0.  Look at g at T + (small increment). -------------
 109       TEMP1 = SIGN(HMING,H)
 110       T0 = T0 + TEMP1
 111       IF ((T0 - TN)*H .LT. 0.0D0) GO TO 230
 112       TEMP2 = TEMP1/H
 113       DO 220 I = 1,N
 114  220    Y(I) = Y(I) + TEMP2*YH(I,2)
 115       GO TO 240
 116  230  CALL DINTDY (T0, 0, YH, NYH, Y, IFLAG)
 117  240  CALL G (NEQ, T0, Y, NGC, G0)
 118       NGE = NGE + 1
 119       ZROOT = .FALSE.
 120       DO 250 I = 1,NGC
 121         IF (ABS(G0(I)) .GT. 0.0D0) GO TO 250
 122         JROOT(I) = 1
 123         ZROOT = .TRUE.
 124  250    CONTINUE
 125       IF (.NOT. ZROOT) GO TO 260
 126 C g has a zero at T0 and also close to T0.  Return root. ---------------
 127       IRT = 1
 128       RETURN
 129 C G0 has no zero components.  Proceed to check relevant interval. ------
 130  260  IF (TN .EQ. TLAST) GO TO 390
 131 C
 132  300  CONTINUE
 133 C Set T1 to TN or TOUTC, whichever comes first, and get g at T1. -------
 134       IF (ITASKC.EQ.2 .OR. ITASKC.EQ.3 .OR. ITASKC.EQ.5) GO TO 310
 135       IF ((TOUTC - TN)*H .GE. 0.0D0) GO TO 310
 136       T1 = TOUTC
 137       IF ((T1 - T0)*H .LE. 0.0D0) GO TO 390
 138       CALL DINTDY (T1, 0, YH, NYH, Y, IFLAG)
 139       GO TO 330
 140  310  T1 = TN
 141       DO 320 I = 1,N
 142  320    Y(I) = YH(I,1)
 143  330  CALL G (NEQ, T1, Y, NGC, G1)
 144       NGE = NGE + 1
 145 C Call DROOTS to search for root in interval from T0 to T1. ------------
 146       JFLAG = 0
 147  350  CONTINUE
 148       CALL DROOTS (NGC, HMING, JFLAG, T0, T1, G0, G1, GX, X, JROOT)
 149       IF (JFLAG .GT. 1) GO TO 360
 150       CALL DINTDY (X, 0, YH, NYH, Y, IFLAG)
 151       CALL G (NEQ, X, Y, NGC, GX)
 152       NGE = NGE + 1
 153       GO TO 350
 154  360  T0 = X
 155       CALL DCOPY (NGC, GX, 1, G0, 1)
 156       IF (JFLAG .EQ. 4) GO TO 390
 157 C Found a root.  Interpolate to X and return. --------------------------
 158       CALL DINTDY (X, 0, YH, NYH, Y, IFLAG)
 159       IRT = 1
 160       RETURN
 161 C
 162  390  CONTINUE
 163       RETURN
 164 C----------------------- End of Subroutine DRCHEK ----------------------
 165       END