1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
|
subroutine roots (ng, hmin, jflag, x0, x1, g0, g1, gx, x, jroot)
clll. optimize
integer ng, jflag, jroot
double precision hmin, x0, x1, g0, g1, gx, x
dimension g0(ng), g1(ng), gx(ng), jroot(ng)
integer iownd3, imax, last, idum3
double precision alpha, x2, rdum3
common /lsr001/ alpha, x2, rdum3(3),
1 iownd3(3), imax, last, idum3(4)
c-----------------------------------------------------------------------
c this subroutine finds the leftmost root of a set of arbitrary
c functions gi(x) (i = 1,...,ng) in an interval (x0,x1). only roots
c of odd multiplicity (i.e. changes of sign of the gi) are found.
c here the sign of x1 - x0 is arbitrary, but is constant for a given
c problem, and -leftmost- means nearest to x0.
c the values of the vector-valued function g(x) = (gi, i=1...ng)
c are communicated through the call sequence of roots.
c the method used is the illinois algorithm.
c
c reference..
c kathie l. hiebert and lawrence f. shampine, implicitly defined
c output points for solutions of ode-s, sandia report sand80-0180,
c february, 1980.
c
c description of parameters.
c
c ng = number of functions gi, or the number of components of
c the vector valued function g(x). input only.
c
c hmin = resolution parameter in x. input only. when a root is
c found, it is located only to within an error of hmin in x.
c typically, hmin should be set to something on the order of
c 100 * uround * max(abs(x0),abs(x1)),
c where uround is the unit roundoff of the machine.
c
c jflag = integer flag for input and output communication.
c
c on input, set jflag = 0 on the first call for the problem,
c and leave it unchanged until the problem is completed.
c (the problem is completed when jflag .ge. 2 on return.)
c
c on output, jflag has the following values and meanings..
c jflag = 1 means roots needs a value of g(x). set gx = g(x)
c and call roots again.
c jflag = 2 means a root has been found. the root is
c at x, and gx contains g(x). (actually, x is the
c rightmost approximation to the root on an interval
c (x0,x1) of size hmin or less.)
c jflag = 3 means x = x1 is a root, with one or more of the gi
c being zero at x1 and no sign changes in (x0,x1).
c gx contains g(x) on output.
c jflag = 4 means no roots (of odd multiplicity) were
c found in (x0,x1) (no sign changes).
c
c x0,x1 = endpoints of the interval where roots are sought.
c x1 and x0 are input when jflag = 0 (first call), and
c must be left unchanged between calls until the problem is
c completed. x0 and x1 must be distinct, but x1 - x0 may be
c of either sign. however, the notion of -left- and -right-
c will be used to mean nearer to x0 or x1, respectively.
c when jflag .ge. 2 on return, x0 and x1 are output, and
c are the endpoints of the relevant interval.
c
c g0,g1 = arrays of length ng containing the vectors g(x0) and g(x1),
c respectively. when jflag = 0, g0 and g1 are input and
c none of the g0(i) should be be zero.
c when jflag .ge. 2 on return, g0 and g1 are output.
c
c gx = array of length ng containing g(x). gx is input
c when jflag = 1, and output when jflag .ge. 2.
c
c x = independent variable value. output only.
c when jflag = 1 on output, x is the point at which g(x)
c is to be evaluated and loaded into gx.
c when jflag = 2 or 3, x is the root.
c when jflag = 4, x is the right endpoint of the interval, x1.
c
c jroot = integer array of length ng. output only.
c when jflag = 2 or 3, jroot indicates which components
c of g(x) have a root at x. jroot(i) is 1 if the i-th
c component has a root, and jroot(i) = 0 otherwise.
c
c note.. this routine uses the common block /lsr001/ to save
c the values of certain variables between calls (own variables).
c-----------------------------------------------------------------------
integer i, imxold, nxlast
double precision t2, tmax, zero
logical zroot, sgnchg, xroot
data zero/0.0d0/
c
if (jflag .eq. 1) go to 200
c jflag .ne. 1. check for change in sign of g or zero at x1. ----------
imax = 0
tmax = zero
zroot = .false.
do 120 i = 1,ng
if (dabs(g1(i)) .gt. zero) go to 110
zroot = .true.
go to 120
c at this point, g0(i) has been checked and cannot be zero. ------------
110 if (dsign(1.0d0,g0(i)) .eq. dsign(1.0d0,g1(i))) go to 120
t2 = dabs(g1(i)/(g1(i)-g0(i)))
if (t2 .le. tmax) go to 120
tmax = t2
imax = i
120 continue
if (imax .gt. 0) go to 130
sgnchg = .false.
go to 140
130 sgnchg = .true.
140 if (.not. sgnchg) go to 400
c there is a sign change. find the first root in the interval. --------
xroot = .false.
nxlast = 0
last = 1
c
c repeat until the first root in the interval is found. loop point. ---
150 continue
if (xroot) go to 300
if (nxlast .eq. last) go to 160
alpha = 1.0d0
go to 180
160 if (last .eq. 0) go to 170
alpha = 0.5d0*alpha
go to 180
170 alpha = 2.0d0*alpha
180 x2 = x1 - (x1-x0)*g1(imax)/(g1(imax) - alpha*g0(imax))
if ((dabs(x2-x0) .lt. hmin) .and.
1 (dabs(x1-x0) .gt. 10.0d0*hmin)) x2 = x0 + 0.1d0*(x1-x0)
jflag = 1
x = x2
c return to the calling routine to get a value of gx = g(x). -----------
return
c check to see in which interval g changes sign. -----------------------
200 imxold = imax
imax = 0
tmax = zero
zroot = .false.
do 220 i = 1,ng
if (dabs(gx(i)) .gt. zero) go to 210
zroot = .true.
go to 220
c neither g0(i) nor gx(i) can be zero at this point. -------------------
210 if (dsign(1.0d0,g0(i)) .eq. dsign(1.0d0,gx(i))) go to 220
t2 = dabs(gx(i)/(gx(i) - g0(i)))
if (t2 .le. tmax) go to 220
tmax = t2
imax = i
220 continue
if (imax .gt. 0) go to 230
sgnchg = .false.
imax = imxold
go to 240
230 sgnchg = .true.
240 nxlast = last
if (.not. sgnchg) go to 250
c sign change between x0 and x2, so replace x1 with x2. ----------------
x1 = x2
call dcopy (ng, gx, 1, g1, 1)
last = 1
xroot = .false.
go to 270
250 if (.not. zroot) go to 260
c zero value at x2 and no sign change in (x0,x2), so x2 is a root. -----
x1 = x2
call dcopy (ng, gx, 1, g1, 1)
xroot = .true.
go to 270
c no sign change between x0 and x2. replace x0 with x2. ---------------
260 continue
call dcopy (ng, gx, 1, g0, 1)
x0 = x2
last = 0
xroot = .false.
270 if (dabs(x1-x0) .le. hmin) xroot = .true.
go to 150
c
c return with x1 as the root. set jroot. set x = x1 and gx = g1. -----
300 jflag = 2
x = x1
call dcopy (ng, g1, 1, gx, 1)
do 320 i = 1,ng
jroot(i) = 0
if (dabs(g1(i)) .gt. zero) go to 310
jroot(i) = 1
go to 320
310 if (dsign(1.0d0,g0(i)) .ne. dsign(1.0d0,g1(i))) jroot(i) = 1
320 continue
return
c
c no sign change in the interval. check for zero at right endpoint. ---
400 if (.not. zroot) go to 420
c
c zero value at x1 and no sign change in (x0,x1). return jflag = 3. ---
x = x1
call dcopy (ng, g1, 1, gx, 1)
do 410 i = 1,ng
jroot(i) = 0
if (dabs(g1(i)) .le. zero) jroot (i) = 1
410 continue
jflag = 3
return
c
c no sign changes in this interval. set x = x1, return jflag = 4. -----
420 call dcopy (ng, g1, 1, gx, 1)
x = x1
jflag = 4
return
c----------------------- end of subroutine roots -----------------------
end
|
