GNU Octave  3.8.0
A high-level interpreted language, primarily intended for numerical computations, mostly compatible with Matlab
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
drchek.f
Go to the documentation of this file.
1  SUBROUTINE drchek (JOB, G, NG, NEQ, TN, TOUT, Y, YP, PHI, PSI,
2  * kold, g0, g1, gx, jroot, irt, uround, info3, rwork, iwork,
3  * rpar, ipar)
4 C
5 C***BEGIN PROLOGUE DRCHEK
6 C***REFER TO DDASRT
7 C***ROUTINES CALLED DDATRP, DROOTS, DCOPY
8 C***DATE WRITTEN 821001 (YYMMDD)
9 C***REVISION DATE 900926 (YYMMDD)
10 C***END PROLOGUE DRCHEK
11 C
12  IMPLICIT DOUBLE PRECISION(a-h,o-z)
13  parameter(lnge=16, lirfnd=18, llast=19, limax=20,
14  * lt0=41, ltlast=42, lalphr=43, lx2=44)
15  EXTERNAL g
16  INTEGER job, ng, neq, kold, jroot, irt, info3, iwork, ipar
17  DOUBLE PRECISION tn, tout, y, yp, phi, psi, g0, g1, gx, uround,
18  * rwork, rpar
19  dimension y(*), yp(*), phi(neq,*), psi(*),
20  1 g0(*), g1(*), gx(*), jroot(*), rwork(*), iwork(*)
21  INTEGER i, jflag
22  DOUBLE PRECISION h
23  DOUBLE PRECISION hming, t1, temp1, temp2, x
24  LOGICAL zroot
25 C-----------------------------------------------------------------------
26 C THIS ROUTINE CHECKS FOR THE PRESENCE OF A ROOT IN THE
27 C VICINITY OF THE CURRENT T, IN A MANNER DEPENDING ON THE
28 C INPUT FLAG JOB. IT CALLS SUBROUTINE DROOTS TO LOCATE THE ROOT
29 C AS PRECISELY AS POSSIBLE.
30 C
31 C IN ADDITION TO VARIABLES DESCRIBED PREVIOUSLY, DRCHEK
32 C USES THE FOLLOWING FOR COMMUNICATION..
33 C JOB = INTEGER FLAG INDICATING TYPE OF CALL..
34 C JOB = 1 MEANS THE PROBLEM IS BEING INITIALIZED, AND DRCHEK
35 C IS TO LOOK FOR A ROOT AT OR VERY NEAR THE INITIAL T.
36 C JOB = 2 MEANS A CONTINUATION CALL TO THE SOLVER WAS JUST
37 C MADE, AND DRCHEK IS TO CHECK FOR A ROOT IN THE
38 C RELEVANT PART OF THE STEP LAST TAKEN.
39 C JOB = 3 MEANS A SUCCESSFUL STEP WAS JUST TAKEN, AND DRCHEK
40 C IS TO LOOK FOR A ROOT IN THE INTERVAL OF THE STEP.
41 C G0 = ARRAY OF LENGTH NG, CONTAINING THE VALUE OF G AT T = T0.
42 C G0 IS INPUT FOR JOB .GE. 2 AND ON OUTPUT IN ALL CASES.
43 C G1,GX = ARRAYS OF LENGTH NG FOR WORK SPACE.
44 C IRT = COMPLETION FLAG..
45 C IRT = 0 MEANS NO ROOT WAS FOUND.
46 C IRT = -1 MEANS JOB = 1 AND A ROOT WAS FOUND TOO NEAR TO T.
47 C IRT = 1 MEANS A LEGITIMATE ROOT WAS FOUND (JOB = 2 OR 3).
48 C ON RETURN, T0 IS THE ROOT LOCATION, AND Y IS THE
49 C CORRESPONDING SOLUTION VECTOR.
50 C T0 = VALUE OF T AT ONE ENDPOINT OF INTERVAL OF INTEREST. ONLY
51 C ROOTS BEYOND T0 IN THE DIRECTION OF INTEGRATION ARE SOUGHT.
52 C T0 IS INPUT IF JOB .GE. 2, AND OUTPUT IN ALL CASES.
53 C T0 IS UPDATED BY DRCHEK, WHETHER A ROOT IS FOUND OR NOT.
54 C STORED IN THE GLOBAL ARRAY RWORK.
55 C TLAST = LAST VALUE OF T RETURNED BY THE SOLVER (INPUT ONLY).
56 C STORED IN THE GLOBAL ARRAY RWORK.
57 C TOUT = FINAL OUTPUT TIME FOR THE SOLVER.
58 C IRFND = INPUT FLAG SHOWING WHETHER THE LAST STEP TAKEN HAD A ROOT.
59 C IRFND = 1 IF IT DID, = 0 IF NOT.
60 C STORED IN THE GLOBAL ARRAY IWORK.
61 C INFO3 = COPY OF INFO(3) (INPUT ONLY).
62 C-----------------------------------------------------------------------
63 C
64  h = psi(1)
65  irt = 0
66  DO 10 i = 1,ng
67  10 jroot(i) = 0
68  hming = (dabs(tn) + dabs(h))*uround*100.0d0
69 C
70  go to(100, 200, 300), job
71 C
72 C EVALUATE G AT INITIAL T (STORED IN RWORK(LT0)), AND CHECK FOR
73 C ZERO VALUES.----------------------------------------------------------
74  100 CONTINUE
75  CALL ddatrp(tn,rwork(lt0),y,yp,neq,kold,phi,psi)
76  CALL g(neq, rwork(lt0), y, ng, g0, rpar, ipar)
77  iwork(lnge) = 1
78  zroot = .false.
79  DO 110 i = 1,ng
80  110 IF (dabs(g0(i)) .LE. 0.0d0) zroot = .true.
81  IF (.NOT. zroot) go to 190
82 C G HAS A ZERO AT T. LOOK AT G AT T + (SMALL INCREMENT). --------------
83  temp1 = dsign(hming,h)
84  rwork(lt0) = rwork(lt0) + temp1
85  temp2 = temp1/h
86  DO 120 i = 1,neq
87  120 y(i) = y(i) + temp2*phi(i,2)
88  CALL g(neq, rwork(lt0), y, ng, g0, rpar, ipar)
89  iwork(lnge) = iwork(lnge) + 1
90  zroot = .false.
91  DO 130 i = 1,ng
92  130 IF (dabs(g0(i)) .LE. 0.0d0) zroot = .true.
93  IF (.NOT. zroot) go to 190
94 C G HAS A ZERO AT T AND ALSO CLOSE TO T. TAKE ERROR RETURN. -----------
95  irt = -1
96  RETURN
97 C
98  190 CONTINUE
99  RETURN
100 C
101 C
102  200 CONTINUE
103  IF (iwork(lirfnd) .EQ. 0) go to 260
104 C IF A ROOT WAS FOUND ON THE PREVIOUS STEP, EVALUATE G0 = G(T0). -------
105  CALL ddatrp(tn, rwork(lt0), y, yp, neq, kold, phi, psi)
106  CALL g(neq, rwork(lt0), y, ng, g0, rpar, ipar)
107  iwork(lnge) = iwork(lnge) + 1
108  zroot = .false.
109  DO 210 i = 1,ng
110  210 IF (dabs(g0(i)) .LE. 0.0d0) zroot = .true.
111  IF (.NOT. zroot) go to 260
112 C G HAS A ZERO AT T0. LOOK AT G AT T + (SMALL INCREMENT). -------------
113  temp1 = dsign(hming,h)
114  rwork(lt0) = rwork(lt0) + temp1
115  IF ((rwork(lt0) - tn)*h .LT. 0.0d0) go to 230
116  temp2 = temp1/h
117  DO 220 i = 1,neq
118  220 y(i) = y(i) + temp2*phi(i,2)
119  go to 240
120  230 CALL ddatrp(tn, rwork(lt0), y, yp, neq, kold, phi, psi)
121  240 CALL g(neq, rwork(lt0), y, ng, g0, rpar, ipar)
122  iwork(lnge) = iwork(lnge) + 1
123  zroot = .false.
124  DO 250 i = 1,ng
125  IF (dabs(g0(i)) .GT. 0.0d0) go to 250
126  jroot(i) = 1
127  zroot = .true.
128  250 CONTINUE
129  IF (.NOT. zroot) go to 260
130 C G HAS A ZERO AT T0 AND ALSO CLOSE TO T0. RETURN ROOT. ---------------
131  irt = 1
132  RETURN
133 C HERE, G0 DOES NOT HAVE A ROOT
134 C G0 HAS NO ZERO COMPONENTS. PROCEED TO CHECK RELEVANT INTERVAL. ------
135  260 IF (tn .EQ. rwork(ltlast)) go to 390
136 C
137  300 CONTINUE
138 C SET T1 TO TN OR TOUT, WHICHEVER COMES FIRST, AND GET G AT T1. --------
139  IF (info3 .EQ. 1) go to 310
140  IF ((tout - tn)*h .GE. 0.0d0) go to 310
141  t1 = tout
142  IF ((t1 - rwork(lt0))*h .LE. 0.0d0) go to 390
143  CALL ddatrp(tn, t1, y, yp, neq, kold, phi, psi)
144  go to 330
145  310 t1 = tn
146  DO 320 i = 1,neq
147  320 y(i) = phi(i,1)
148  330 CALL g(neq, t1, y, ng, g1, rpar, ipar)
149  iwork(lnge) = iwork(lnge) + 1
150 C CALL DROOTS TO SEARCH FOR ROOT IN INTERVAL FROM T0 TO T1. ------------
151  jflag = 0
152  350 CONTINUE
153  CALL droots(ng, hming, jflag, rwork(lt0), t1, g0, g1, gx, x,
154  * jroot, iwork(limax), iwork(llast), rwork(lalphr),
155  * rwork(lx2))
156  IF (jflag .GT. 1) go to 360
157  CALL ddatrp(tn, x, y, yp, neq, kold, phi, psi)
158  CALL g(neq, x, y, ng, gx, rpar, ipar)
159  iwork(lnge) = iwork(lnge) + 1
160  go to 350
161  360 rwork(lt0) = x
162  CALL dcopy(ng, gx, 1, g0, 1)
163  IF (jflag .EQ. 4) go to 390
164 C FOUND A ROOT. INTERPOLATE TO X AND RETURN. --------------------------
165  CALL ddatrp(tn, x, y, yp, neq, kold, phi, psi)
166  irt = 1
167  RETURN
168 C
169  390 CONTINUE
170  RETURN
171 C---------------------- END OF SUBROUTINE DRCHEK -----------------------
172  END