GNU Octave  6.2.0
A high-level interpreted language, primarily intended for numerical computations, mostly compatible with Matlab
droots.f
Go to the documentation of this file.
1  SUBROUTINE droots (NG, HMIN, JFLAG, X0, X1, G0, G1, GX, X, JROOT,
2  * IMAX, LAST, ALPHA, X2)
3 C
4 C***BEGIN PROLOGUE DROOTS
5 C***REFER TO DDASRT
6 C***ROUTINES CALLED DCOPY
7 C***DATE WRITTEN 821001 (YYMMDD)
8 C***REVISION DATE 900926 (YYMMDD)
9 C***END PROLOGUE DROOTS
10 C
11  IMPLICIT DOUBLE PRECISION(a-h,o-z)
12  INTEGER NG, JFLAG, JROOT, IMAX, LAST
13  DOUBLE PRECISION HMIN, X0, X1, G0, G1, GX, X, ALPHA, X2
14  dimension g0(ng), g1(ng), gx(ng), jroot(ng)
15 C-----------------------------------------------------------------------
16 C THIS SUBROUTINE FINDS THE LEFTMOST ROOT OF A SET OF ARBITRARY
17 C FUNCTIONS GI(X) (I = 1,...,NG) IN AN INTERVAL (X0,X1). ONLY ROOTS
18 C OF ODD MULTIPLICITY (I.E. CHANGES OF SIGN OF THE GI) ARE FOUND.
19 C HERE THE SIGN OF X1 - X0 IS ARBITRARY, BUT IS CONSTANT FOR A GIVEN
20 C PROBLEM, AND -LEFTMOST- MEANS NEAREST TO X0.
21 C THE VALUES OF THE VECTOR-VALUED FUNCTION G(X) = (GI, I=1...NG)
22 C ARE COMMUNICATED THROUGH THE CALL SEQUENCE OF DROOTS.
23 C THE METHOD USED IS THE ILLINOIS ALGORITHM.
24 C
25 C REFERENCE..
26 C KATHIE L. HIEBERT AND LAWRENCE F. SHAMPINE, IMPLICITLY DEFINED
27 C OUTPUT POINTS FOR SOLUTIONS OF ODE-S, SANDIA REPORT SAND80-0180,
28 C FEBRUARY, 1980.
29 C
30 C DESCRIPTION OF PARAMETERS.
31 C
32 C NG = NUMBER OF FUNCTIONS GI, OR THE NUMBER OF COMPONENTS OF
33 C THE VECTOR VALUED FUNCTION G(X). INPUT ONLY.
34 C
35 C HMIN = RESOLUTION PARAMETER IN X. INPUT ONLY. WHEN A ROOT IS
36 C FOUND, IT IS LOCATED ONLY TO WITHIN AN ERROR OF HMIN IN X.
37 C TYPICALLY, HMIN SHOULD BE SET TO SOMETHING ON THE ORDER OF
38 C 100 * UROUND * MAX(ABS(X0),ABS(X1)),
39 C WHERE UROUND IS THE UNIT ROUNDOFF OF THE MACHINE.
40 C
41 C JFLAG = INTEGER FLAG FOR INPUT AND OUTPUT COMMUNICATION.
42 C
43 C ON INPUT, SET JFLAG = 0 ON THE FIRST CALL FOR THE PROBLEM,
44 C AND LEAVE IT UNCHANGED UNTIL THE PROBLEM IS COMPLETED.
45 C (THE PROBLEM IS COMPLETED WHEN JFLAG .GE. 2 ON RETURN.)
46 C
47 C ON OUTPUT, JFLAG HAS THE FOLLOWING VALUES AND MEANINGS..
48 C JFLAG = 1 MEANS DROOTS NEEDS A VALUE OF G(X). SET GX = G(X)
49 C AND CALL DROOTS AGAIN.
50 C JFLAG = 2 MEANS A ROOT HAS BEEN FOUND. THE ROOT IS
51 C AT X, AND GX CONTAINS G(X). (ACTUALLY, X IS THE
52 C RIGHTMOST APPROXIMATION TO THE ROOT ON AN INTERVAL
53 C (X0,X1) OF SIZE HMIN OR LESS.)
54 C JFLAG = 3 MEANS X = X1 IS A ROOT, WITH ONE OR MORE OF THE GI
55 C BEING ZERO AT X1 AND NO SIGN CHANGES IN (X0,X1).
56 C GX CONTAINS G(X) ON OUTPUT.
57 C JFLAG = 4 MEANS NO ROOTS (OF ODD MULTIPLICITY) WERE
58 C FOUND IN (X0,X1) (NO SIGN CHANGES).
59 C
60 C X0,X1 = ENDPOINTS OF THE INTERVAL WHERE ROOTS ARE SOUGHT.
61 C X1 AND X0 ARE INPUT WHEN JFLAG = 0 (FIRST CALL), AND
62 C MUST BE LEFT UNCHANGED BETWEEN CALLS UNTIL THE PROBLEM IS
63 C COMPLETED. X0 AND X1 MUST BE DISTINCT, BUT X1 - X0 MAY BE
64 C OF EITHER SIGN. HOWEVER, THE NOTION OF -LEFT- AND -RIGHT-
65 C WILL BE USED TO MEAN NEARER TO X0 OR X1, RESPECTIVELY.
66 C WHEN JFLAG .GE. 2 ON RETURN, X0 AND X1 ARE OUTPUT, AND
67 C ARE THE ENDPOINTS OF THE RELEVANT INTERVAL.
68 C
69 C G0,G1 = ARRAYS OF LENGTH NG CONTAINING THE VECTORS G(X0) AND G(X1),
70 C RESPECTIVELY. WHEN JFLAG = 0, G0 AND G1 ARE INPUT AND
71 C NONE OF THE G0(I) SHOULD BE BE ZERO.
72 C WHEN JFLAG .GE. 2 ON RETURN, G0 AND G1 ARE OUTPUT.
73 C
74 C GX = ARRAY OF LENGTH NG CONTAINING G(X). GX IS INPUT
75 C WHEN JFLAG = 1, AND OUTPUT WHEN JFLAG .GE. 2.
76 C
77 C X = INDEPENDENT VARIABLE VALUE. OUTPUT ONLY.
78 C WHEN JFLAG = 1 ON OUTPUT, X IS THE POINT AT WHICH G(X)
79 C IS TO BE EVALUATED AND LOADED INTO GX.
80 C WHEN JFLAG = 2 OR 3, X IS THE ROOT.
81 C WHEN JFLAG = 4, X IS THE RIGHT ENDPOINT OF THE INTERVAL, X1.
82 C
83 C JROOT = INTEGER ARRAY OF LENGTH NG. OUTPUT ONLY.
84 C WHEN JFLAG = 2 OR 3, JROOT INDICATES WHICH COMPONENTS
85 C OF G(X) HAVE A ROOT AT X. JROOT(I) IS 1 IF THE I-TH
86 C COMPONENT HAS A ROOT, AND JROOT(I) = 0 OTHERWISE.
87 C
88 C IMAX, LAST, ALPHA, X2 =
89 C BOOKKEEPING VARIABLES WHICH MUST BE SAVED FROM CALL
90 C TO CALL. THEY ARE SAVED INSIDE THE CALLING ROUTINE,
91 C BUT THEY ARE USED ONLY WITHIN THIS ROUTINE.
92 C-----------------------------------------------------------------------
93  INTEGER I, IMXOLD, NXLAST
94  DOUBLE PRECISION T2, TMAX, ZERO
95  LOGICAL ZROOT, SGNCHG, XROOT
96  DATA zero/0.0d0/
97 C
98  IF (jflag .EQ. 1) GO TO 200
99 C JFLAG .NE. 1. CHECK FOR CHANGE IN SIGN OF G OR ZERO AT X1. ----------
100  imax = 0
101  tmax = zero
102  zroot = .false.
103  DO 120 i = 1,ng
104  IF (dabs(g1(i)) .GT. zero) GO TO 110
105  zroot = .true.
106  GO TO 120
107 C AT THIS POINT, G0(I) HAS BEEN CHECKED AND CANNOT BE ZERO. ------------
108  110 IF (dsign(1.0d0,g0(i)) .EQ. dsign(1.0d0,g1(i))) GO TO 120
109  t2 = dabs(g1(i)/(g1(i)-g0(i)))
110  IF (t2 .LE. tmax) GO TO 120
111  tmax = t2
112  imax = i
113  120 CONTINUE
114  IF (imax .GT. 0) GO TO 130
115  sgnchg = .false.
116  GO TO 140
117  130 sgnchg = .true.
118  140 IF (.NOT. sgnchg) GO TO 400
119 C THERE IS A SIGN CHANGE. FIND THE FIRST ROOT IN THE INTERVAL. --------
120  xroot = .false.
121  nxlast = 0
122  last = 1
123 C
124 C REPEAT UNTIL THE FIRST ROOT IN THE INTERVAL IS FOUND. LOOP POINT. ---
125  150 CONTINUE
126  IF (xroot) GO TO 300
127  IF (nxlast .EQ. last) GO TO 160
128  alpha = 1.0d0
129  GO TO 180
130  160 IF (last .EQ. 0) GO TO 170
131  alpha = 0.5d0*alpha
132  GO TO 180
133  170 alpha = 2.0d0*alpha
134  180 x2 = x1 - (x1-x0)*g1(imax)/(g1(imax) - alpha*g0(imax))
135  IF ((dabs(x2-x0) .LT. hmin) .AND.
136  1 (dabs(x1-x0) .GT. 10.0d0*hmin)) x2 = x0 + 0.1d0*(x1-x0)
137  jflag = 1
138  x = x2
139 C RETURN TO THE CALLING ROUTINE TO GET A VALUE OF GX = G(X). -----------
140  RETURN
141 C CHECK TO SEE IN WHICH INTERVAL G CHANGES SIGN. -----------------------
142  200 imxold = imax
143  imax = 0
144  tmax = zero
145  zroot = .false.
146  DO 220 i = 1,ng
147  IF (dabs(gx(i)) .GT. zero) GO TO 210
148  zroot = .true.
149  GO TO 220
150 C NEITHER G0(I) NOR GX(I) CAN BE ZERO AT THIS POINT. -------------------
151  210 IF (dsign(1.0d0,g0(i)) .EQ. dsign(1.0d0,gx(i))) GO TO 220
152  t2 = dabs(gx(i)/(gx(i) - g0(i)))
153  IF (t2 .LE. tmax) GO TO 220
154  tmax = t2
155  imax = i
156  220 CONTINUE
157  IF (imax .GT. 0) GO TO 230
158  sgnchg = .false.
159  imax = imxold
160  GO TO 240
161  230 sgnchg = .true.
162  240 nxlast = last
163  IF (.NOT. sgnchg) GO TO 250
164 C SIGN CHANGE BETWEEN X0 AND X2, SO REPLACE X1 WITH X2. ----------------
165  x1 = x2
166  CALL dcopy (ng, gx, 1, g1, 1)
167  last = 1
168  xroot = .false.
169  GO TO 270
170  250 IF (.NOT. zroot) GO TO 260
171 C ZERO VALUE AT X2 AND NO SIGN CHANGE IN (X0,X2), SO X2 IS A ROOT. -----
172  x1 = x2
173  CALL dcopy (ng, gx, 1, g1, 1)
174  xroot = .true.
175  GO TO 270
176 C NO SIGN CHANGE BETWEEN X0 AND X2. REPLACE X0 WITH X2. ---------------
177  260 CONTINUE
178  CALL dcopy (ng, gx, 1, g0, 1)
179  x0 = x2
180  last = 0
181  xroot = .false.
182  270 IF (dabs(x1-x0) .LE. hmin) xroot = .true.
183  GO TO 150
184 C
185 C RETURN WITH X1 AS THE ROOT. SET JROOT. SET X = X1 AND GX = G1. -----
186  300 jflag = 2
187  x = x1
188  CALL dcopy (ng, g1, 1, gx, 1)
189  DO 320 i = 1,ng
190  jroot(i) = 0
191  IF (dabs(g1(i)) .GT. zero) GO TO 310
192  jroot(i) = 1
193  GO TO 320
194  310 IF (dsign(1.0d0,g0(i)) .NE. dsign(1.0d0,g1(i))) jroot(i) = 1
195  320 CONTINUE
196  RETURN
197 C
198 C NO SIGN CHANGE IN THE INTERVAL. CHECK FOR ZERO AT RIGHT ENDPOINT. ---
199  400 IF (.NOT. zroot) GO TO 420
200 C
201 C ZERO VALUE AT X1 AND NO SIGN CHANGE IN (X0,X1). RETURN JFLAG = 3. ---
202  x = x1
203  CALL dcopy (ng, g1, 1, gx, 1)
204  DO 410 i = 1,ng
205  jroot(i) = 0
206  IF (dabs(g1(i)) .LE. zero) jroot(i) = 1
207  410 CONTINUE
208  jflag = 3
209  RETURN
210 C
211 C NO SIGN CHANGES IN THIS INTERVAL. SET X = X1, RETURN JFLAG = 4. -----
212  420 CALL dcopy (ng, g1, 1, gx, 1)
213  x = x1
214  jflag = 4
215  RETURN
216 C---------------------- END OF SUBROUTINE DROOTS -----------------------
217  END
subroutine droots(NG, HMIN, JFLAG, X0, X1, G0, G1, GX, X, JROOT, IMAX, LAST, ALPHA, X2)
Definition: droots.f:3