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
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