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