drchek.f

Go to the documentation of this file.

      SUBROUTINE drchek (JOB, G, NG, NEQ, TN, TOUT, Y, YP, PHI, PSI,
     *  kold, g0, g1, gx, jroot, irt, uround, info3, rwork, iwork,
     *  rpar, ipar)
C
C***BEGIN PROLOGUE  DRCHEK
C***REFER TO DDASRT
C***ROUTINES CALLED  DDATRP, DROOTS, DCOPY
C***DATE WRITTEN   821001   (YYMMDD)
C***REVISION DATE  900926   (YYMMDD)
C***END PROLOGUE  DRCHEK
C
      IMPLICIT DOUBLE PRECISION(a-h,o-z)
      parameter(lnge=16, lirfnd=18, llast=19, limax=20,
     *           lt0=41, ltlast=42, lalphr=43, lx2=44)
      EXTERNAL g
      INTEGER JOB, NG, NEQ, KOLD, JROOT, IRT, INFO3, IWORK, IPAR
      DOUBLE PRECISION TN, TOUT, Y, YP, PHI, PSI, G0, G1, GX, UROUND,
     *  rwork, rpar
      dimension  y(*), yp(*), phi(neq,*), psi(*),
     1  g0(*), g1(*), gx(*), jroot(*), rwork(*), iwork(*)
      INTEGER I, JFLAG
      DOUBLE PRECISION H
      DOUBLE PRECISION HMING, T1, TEMP1, TEMP2, X
      LOGICAL ZROOT
C-----------------------------------------------------------------------
C THIS ROUTINE CHECKS FOR THE PRESENCE OF A ROOT IN THE
C VICINITY OF THE CURRENT T, IN A MANNER DEPENDING ON THE
C INPUT FLAG JOB.  IT CALLS SUBROUTINE DROOTS TO LOCATE THE ROOT
C AS PRECISELY AS POSSIBLE.
C
C IN ADDITION TO VARIABLES DESCRIBED PREVIOUSLY, DRCHEK
C USES THE FOLLOWING FOR COMMUNICATION..
C JOB    = INTEGER FLAG INDICATING TYPE OF CALL..
C          JOB = 1 MEANS THE PROBLEM IS BEING INITIALIZED, AND DRCHEK
C                  IS TO LOOK FOR A ROOT AT OR VERY NEAR THE INITIAL T.
C          JOB = 2 MEANS A CONTINUATION CALL TO THE SOLVER WAS JUST
C                  MADE, AND DRCHEK IS TO CHECK FOR A ROOT IN THE
C                  RELEVANT PART OF THE STEP LAST TAKEN.
C          JOB = 3 MEANS A SUCCESSFUL STEP WAS JUST TAKEN, AND DRCHEK
C                  IS TO LOOK FOR A ROOT IN THE INTERVAL OF THE STEP.
C G0     = ARRAY OF LENGTH NG, CONTAINING THE VALUE OF G AT T = T0.
C          G0 IS INPUT FOR JOB .GE. 2 AND ON OUTPUT IN ALL CASES.
C G1,GX  = ARRAYS OF LENGTH NG FOR WORK SPACE.
C IRT    = COMPLETION FLAG..
C          IRT = 0  MEANS NO ROOT WAS FOUND.
C          IRT = -1 MEANS JOB = 1 AND A ROOT WAS FOUND TOO NEAR TO T.
C          IRT = 1  MEANS A LEGITIMATE ROOT WAS FOUND (JOB = 2 OR 3).
C                   ON RETURN, T0 IS THE ROOT LOCATION, AND Y IS THE
C                   CORRESPONDING SOLUTION VECTOR.
C T0     = VALUE OF T AT ONE ENDPOINT OF INTERVAL OF INTEREST.  ONLY
C          ROOTS BEYOND T0 IN THE DIRECTION OF INTEGRATION ARE SOUGHT.
C          T0 IS INPUT IF JOB .GE. 2, AND OUTPUT IN ALL CASES.
C          T0 IS UPDATED BY DRCHEK, WHETHER A ROOT IS FOUND OR NOT.
C          STORED IN THE GLOBAL ARRAY RWORK.
C TLAST  = LAST VALUE OF T RETURNED BY THE SOLVER (INPUT ONLY).
C          STORED IN THE GLOBAL ARRAY RWORK.
C TOUT   = FINAL OUTPUT TIME FOR THE SOLVER.
C IRFND  = INPUT FLAG SHOWING WHETHER THE LAST STEP TAKEN HAD A ROOT.
C          IRFND = 1 IF IT DID, = 0 IF NOT.
C          STORED IN THE GLOBAL ARRAY IWORK.
C INFO3  = COPY OF INFO(3) (INPUT ONLY).
C-----------------------------------------------------------------------
C
      h = psi(1)
      irt = 0
      DO 10 i = 1,ng
 10     jroot(i) = 0
      hming = (dabs(tn) + dabs(h))*uround*100.0d0
C
      go to(100, 200, 300), job
C
C EVALUATE G AT INITIAL T (STORED IN RWORK(LT0)), AND CHECK FOR
C ZERO VALUES.----------------------------------------------------------
 100  CONTINUE
      CALL ddatrp(tn,rwork(lt0),y,yp,neq,kold,phi,psi)
      CALL g(neq, rwork(lt0), y, ng, g0, rpar, ipar)
      iwork(lnge) = 1
      zroot = .false.
      DO 110 i = 1,ng
 110    IF (dabs(g0(i)) .LE. 0.0d0) zroot = .true.
      IF (.NOT. zroot) go to 190
C G HAS A ZERO AT T.  LOOK AT G AT T + (SMALL INCREMENT). --------------
      temp1 = dsign(hming,h)
      rwork(lt0) = rwork(lt0) + temp1
      temp2 = temp1/h
      DO 120 i = 1,neq
 120    y(i) = y(i) + temp2*phi(i,2)
      CALL g(neq, rwork(lt0), y, ng, g0, rpar, ipar)
      iwork(lnge) = iwork(lnge) + 1
      zroot = .false.
      DO 130 i = 1,ng
 130    IF (dabs(g0(i)) .LE. 0.0d0) zroot = .true.
      IF (.NOT. zroot) go to 190
C G HAS A ZERO AT T AND ALSO CLOSE TO T.  TAKE ERROR RETURN. -----------
      irt = -1
      RETURN
C
 190  CONTINUE
      RETURN
C
C
 200  CONTINUE
      IF (iwork(lirfnd) .EQ. 0) go to 260
C IF A ROOT WAS FOUND ON THE PREVIOUS STEP, EVALUATE G0 = G(T0). -------
      CALL ddatrp(tn, rwork(lt0), y, yp, neq, kold, phi, psi)
      CALL g(neq, rwork(lt0), y, ng, g0, rpar, ipar)
      iwork(lnge) = iwork(lnge) + 1
      zroot = .false.
      DO 210 i = 1,ng
 210    IF (dabs(g0(i)) .LE. 0.0d0) zroot = .true.
      IF (.NOT. zroot) go to 260
C G HAS A ZERO AT T0.  LOOK AT G AT T + (SMALL INCREMENT). -------------
      temp1 = dsign(hming,h)
      rwork(lt0) = rwork(lt0) + temp1
      IF ((rwork(lt0) - tn)*h .LT. 0.0d0) go to 230
      temp2 = temp1/h
      DO 220 i = 1,neq
 220    y(i) = y(i) + temp2*phi(i,2)
      go to 240
 230  CALL ddatrp(tn, rwork(lt0), y, yp, neq, kold, phi, psi)
 240  CALL g(neq, rwork(lt0), y, ng, g0, rpar, ipar)
      iwork(lnge) = iwork(lnge) + 1
      zroot = .false.
      DO 250 i = 1,ng
        IF (dabs(g0(i)) .GT. 0.0d0) go to 250
        jroot(i) = 1
        zroot = .true.
 250    CONTINUE
      IF (.NOT. zroot) go to 260
C G HAS A ZERO AT T0 AND ALSO CLOSE TO T0.  RETURN ROOT. ---------------
      irt = 1
      RETURN
C     HERE, G0 DOES NOT HAVE A ROOT
C G0 HAS NO ZERO COMPONENTS.  PROCEED TO CHECK RELEVANT INTERVAL. ------
 260  IF (tn .EQ. rwork(ltlast)) go to 390
C
 300  CONTINUE
C SET T1 TO TN OR TOUT, WHICHEVER COMES FIRST, AND GET G AT T1. --------
      IF (info3 .EQ. 1) go to 310
      IF ((tout - tn)*h .GE. 0.0d0) go to 310
      t1 = tout
      IF ((t1 - rwork(lt0))*h .LE. 0.0d0) go to 390
      CALL ddatrp(tn, t1, y, yp, neq, kold, phi, psi)
      go to 330
 310  t1 = tn
      DO 320 i = 1,neq
 320    y(i) = phi(i,1)
 330  CALL g(neq, t1, y, ng, g1, rpar, ipar)
      iwork(lnge) = iwork(lnge) + 1
C CALL DROOTS TO SEARCH FOR ROOT IN INTERVAL FROM T0 TO T1. ------------
      jflag = 0
 350  CONTINUE
      CALL droots(ng, hming, jflag, rwork(lt0), t1, g0, g1, gx, x,
     *             jroot, iwork(limax), iwork(llast), rwork(lalphr),
     *             rwork(lx2))
      IF (jflag .GT. 1) go to 360
      CALL ddatrp(tn, x, y, yp, neq, kold, phi, psi)
      CALL g(neq, x, y, ng, gx, rpar, ipar)
      iwork(lnge) = iwork(lnge) + 1
      go to 350
 360  rwork(lt0) = x
      CALL dcopy(ng, gx, 1, g0, 1)
      IF (jflag .EQ. 4) go to 390
C FOUND A ROOT.  INTERPOLATE TO X AND RETURN. --------------------------
      CALL ddatrp(tn, x, y, yp, neq, kold, phi, psi)
      irt = 1
      RETURN
C
 390  CONTINUE
      RETURN
C---------------------- END OF SUBROUTINE DRCHEK -----------------------
      END