GNU Octave: libcruft/dasrt/drchek.f Source File

Go to the documentation of this file.
00001       SUBROUTINE DRCHEK (JOB, G, NG, NEQ, TN, TOUT, Y, YP, PHI, PSI,
00002      *  KOLD, G0, G1, GX, JROOT, IRT, UROUND, INFO3, RWORK, IWORK,
00003      *  RPAR, IPAR)
00004 C
00005 C***BEGIN PROLOGUE  DRCHEK
00006 C***REFER TO DDASRT
00007 C***ROUTINES CALLED  DDATRP, DROOTS, DCOPY
00008 C***DATE WRITTEN   821001   (YYMMDD)
00009 C***REVISION DATE  900926   (YYMMDD)
00010 C***END PROLOGUE  DRCHEK
00011 C
00012       IMPLICIT DOUBLE PRECISION(A-H,O-Z)
00013       PARAMETER (LNGE=16, LIRFND=18, LLAST=19, LIMAX=20,
00014      *           LT0=41, LTLAST=42, LALPHR=43, LX2=44)
00015       EXTERNAL G
00016       INTEGER JOB, NG, NEQ, KOLD, JROOT, IRT, INFO3, IWORK, IPAR
00017       DOUBLE PRECISION TN, TOUT, Y, YP, PHI, PSI, G0, G1, GX, UROUND,
00018      *  RWORK, RPAR
00019       DIMENSION  Y(*), YP(*), PHI(NEQ,*), PSI(*),
00020      1  G0(*), G1(*), GX(*), JROOT(*), RWORK(*), IWORK(*)
00021       INTEGER I, JFLAG
00022       DOUBLE PRECISION H
00023       DOUBLE PRECISION HMING, T1, TEMP1, TEMP2, X
00024       LOGICAL ZROOT
00025 C-----------------------------------------------------------------------
00026 C THIS ROUTINE CHECKS FOR THE PRESENCE OF A ROOT IN THE
00027 C VICINITY OF THE CURRENT T, IN A MANNER DEPENDING ON THE
00028 C INPUT FLAG JOB.  IT CALLS SUBROUTINE DROOTS TO LOCATE THE ROOT
00029 C AS PRECISELY AS POSSIBLE.
00030 C
00031 C IN ADDITION TO VARIABLES DESCRIBED PREVIOUSLY, DRCHEK
00032 C USES THE FOLLOWING FOR COMMUNICATION..
00033 C JOB    = INTEGER FLAG INDICATING TYPE OF CALL..
00034 C          JOB = 1 MEANS THE PROBLEM IS BEING INITIALIZED, AND DRCHEK
00035 C                  IS TO LOOK FOR A ROOT AT OR VERY NEAR THE INITIAL T.
00036 C          JOB = 2 MEANS A CONTINUATION CALL TO THE SOLVER WAS JUST
00037 C                  MADE, AND DRCHEK IS TO CHECK FOR A ROOT IN THE
00038 C                  RELEVANT PART OF THE STEP LAST TAKEN.
00039 C          JOB = 3 MEANS A SUCCESSFUL STEP WAS JUST TAKEN, AND DRCHEK
00040 C                  IS TO LOOK FOR A ROOT IN THE INTERVAL OF THE STEP.
00041 C G0     = ARRAY OF LENGTH NG, CONTAINING THE VALUE OF G AT T = T0.
00042 C          G0 IS INPUT FOR JOB .GE. 2 AND ON OUTPUT IN ALL CASES.
00043 C G1,GX  = ARRAYS OF LENGTH NG FOR WORK SPACE.
00044 C IRT    = COMPLETION FLAG..
00045 C          IRT = 0  MEANS NO ROOT WAS FOUND.
00046 C          IRT = -1 MEANS JOB = 1 AND A ROOT WAS FOUND TOO NEAR TO T.
00047 C          IRT = 1  MEANS A LEGITIMATE ROOT WAS FOUND (JOB = 2 OR 3).
00048 C                   ON RETURN, T0 IS THE ROOT LOCATION, AND Y IS THE
00049 C                   CORRESPONDING SOLUTION VECTOR.
00050 C T0     = VALUE OF T AT ONE ENDPOINT OF INTERVAL OF INTEREST.  ONLY
00051 C          ROOTS BEYOND T0 IN THE DIRECTION OF INTEGRATION ARE SOUGHT.
00052 C          T0 IS INPUT IF JOB .GE. 2, AND OUTPUT IN ALL CASES.
00053 C          T0 IS UPDATED BY DRCHEK, WHETHER A ROOT IS FOUND OR NOT.
00054 C          STORED IN THE GLOBAL ARRAY RWORK.
00055 C TLAST  = LAST VALUE OF T RETURNED BY THE SOLVER (INPUT ONLY).
00056 C          STORED IN THE GLOBAL ARRAY RWORK.
00057 C TOUT   = FINAL OUTPUT TIME FOR THE SOLVER.
00058 C IRFND  = INPUT FLAG SHOWING WHETHER THE LAST STEP TAKEN HAD A ROOT.
00059 C          IRFND = 1 IF IT DID, = 0 IF NOT.
00060 C          STORED IN THE GLOBAL ARRAY IWORK.
00061 C INFO3  = COPY OF INFO(3) (INPUT ONLY).
00062 C-----------------------------------------------------------------------
00063 C     
00064       H = PSI(1)
00065       IRT = 0
00066       DO 10 I = 1,NG
00067  10     JROOT(I) = 0
00068       HMING = (DABS(TN) + DABS(H))*UROUND*100.0D0
00069 C
00070       GO TO (100, 200, 300), JOB
00071 C
00072 C EVALUATE G AT INITIAL T (STORED IN RWORK(LT0)), AND CHECK FOR
00073 C ZERO VALUES.----------------------------------------------------------
00074  100  CONTINUE
00075       CALL DDATRP(TN,RWORK(LT0),Y,YP,NEQ,KOLD,PHI,PSI)
00076       CALL G (NEQ, RWORK(LT0), Y, NG, G0, RPAR, IPAR)
00077       IWORK(LNGE) = 1
00078       ZROOT = .FALSE.
00079       DO 110 I = 1,NG
00080  110    IF (DABS(G0(I)) .LE. 0.0D0) ZROOT = .TRUE.
00081       IF (.NOT. ZROOT) GO TO 190
00082 C G HAS A ZERO AT T.  LOOK AT G AT T + (SMALL INCREMENT). --------------
00083       TEMP1 = DSIGN(HMING,H)
00084       RWORK(LT0) = RWORK(LT0) + TEMP1
00085       TEMP2 = TEMP1/H
00086       DO 120 I = 1,NEQ
00087  120    Y(I) = Y(I) + TEMP2*PHI(I,2)
00088       CALL G (NEQ, RWORK(LT0), Y, NG, G0, RPAR, IPAR)
00089       IWORK(LNGE) = IWORK(LNGE) + 1
00090       ZROOT = .FALSE.
00091       DO 130 I = 1,NG
00092  130    IF (DABS(G0(I)) .LE. 0.0D0) ZROOT = .TRUE.
00093       IF (.NOT. ZROOT) GO TO 190
00094 C G HAS A ZERO AT T AND ALSO CLOSE TO T.  TAKE ERROR RETURN. -----------
00095       IRT = -1
00096       RETURN
00097 C
00098  190  CONTINUE
00099       RETURN
00100 C
00101 C
00102  200  CONTINUE
00103       IF (IWORK(LIRFND) .EQ. 0) GO TO 260
00104 C IF A ROOT WAS FOUND ON THE PREVIOUS STEP, EVALUATE G0 = G(T0). -------
00105       CALL DDATRP (TN, RWORK(LT0), Y, YP, NEQ, KOLD, PHI, PSI)
00106       CALL G (NEQ, RWORK(LT0), Y, NG, G0, RPAR, IPAR)
00107       IWORK(LNGE) = IWORK(LNGE) + 1
00108       ZROOT = .FALSE.
00109       DO 210 I = 1,NG
00110  210    IF (DABS(G0(I)) .LE. 0.0D0) ZROOT = .TRUE.
00111       IF (.NOT. ZROOT) GO TO 260
00112 C G HAS A ZERO AT T0.  LOOK AT G AT T + (SMALL INCREMENT). -------------
00113       TEMP1 = DSIGN(HMING,H)
00114       RWORK(LT0) = RWORK(LT0) + TEMP1
00115       IF ((RWORK(LT0) - TN)*H .LT. 0.0D0) GO TO 230
00116       TEMP2 = TEMP1/H
00117       DO 220 I = 1,NEQ
00118  220    Y(I) = Y(I) + TEMP2*PHI(I,2)
00119       GO TO 240
00120  230  CALL DDATRP (TN, RWORK(LT0), Y, YP, NEQ, KOLD, PHI, PSI)
00121  240  CALL G (NEQ, RWORK(LT0), Y, NG, G0, RPAR, IPAR)
00122       IWORK(LNGE) = IWORK(LNGE) + 1
00123       ZROOT = .FALSE.
00124       DO 250 I = 1,NG
00125         IF (DABS(G0(I)) .GT. 0.0D0) GO TO 250
00126         JROOT(I) = 1
00127         ZROOT = .TRUE.
00128  250    CONTINUE
00129       IF (.NOT. ZROOT) GO TO 260
00130 C G HAS A ZERO AT T0 AND ALSO CLOSE TO T0.  RETURN ROOT. ---------------
00131       IRT = 1
00132       RETURN
00133 C     HERE, G0 DOES NOT HAVE A ROOT
00134 C G0 HAS NO ZERO COMPONENTS.  PROCEED TO CHECK RELEVANT INTERVAL. ------
00135  260  IF (TN .EQ. RWORK(LTLAST)) GO TO 390
00136 C
00137  300  CONTINUE
00138 C SET T1 TO TN OR TOUT, WHICHEVER COMES FIRST, AND GET G AT T1. --------
00139       IF (INFO3 .EQ. 1) GO TO 310
00140       IF ((TOUT - TN)*H .GE. 0.0D0) GO TO 310
00141       T1 = TOUT
00142       IF ((T1 - RWORK(LT0))*H .LE. 0.0D0) GO TO 390
00143       CALL DDATRP (TN, T1, Y, YP, NEQ, KOLD, PHI, PSI)
00144       GO TO 330
00145  310  T1 = TN
00146       DO 320 I = 1,NEQ
00147  320    Y(I) = PHI(I,1)
00148  330  CALL G (NEQ, T1, Y, NG, G1, RPAR, IPAR)
00149       IWORK(LNGE) = IWORK(LNGE) + 1
00150 C CALL DROOTS TO SEARCH FOR ROOT IN INTERVAL FROM T0 TO T1. ------------
00151       JFLAG = 0
00152  350  CONTINUE
00153       CALL DROOTS (NG, HMING, JFLAG, RWORK(LT0), T1, G0, G1, GX, X,
00154      *             JROOT, IWORK(LIMAX), IWORK(LLAST), RWORK(LALPHR),
00155      *             RWORK(LX2))
00156       IF (JFLAG .GT. 1) GO TO 360
00157       CALL DDATRP (TN, X, Y, YP, NEQ, KOLD, PHI, PSI)
00158       CALL G (NEQ, X, Y, NG, GX, RPAR, IPAR)
00159       IWORK(LNGE) = IWORK(LNGE) + 1
00160       GO TO 350
00161  360  RWORK(LT0) = X
00162       CALL DCOPY (NG, GX, 1, G0, 1)
00163       IF (JFLAG .EQ. 4) GO TO 390
00164 C FOUND A ROOT.  INTERPOLATE TO X AND RETURN. --------------------------
00165       CALL DDATRP (TN, X, Y, YP, NEQ, KOLD, PHI, PSI)
00166       IRT = 1
00167       RETURN
00168 C
00169  390  CONTINUE
00170       RETURN
00171 C---------------------- END OF SUBROUTINE DRCHEK -----------------------
00172       END