1
subroutine roots (ng, hmin, jflag, x0, x1, g0, g1, gx, x, jroot)
3
integer ng, jflag, jroot
4
double precision hmin, x0, x1, g0, g1, gx, x
5
dimension g0(ng), g1(ng), gx(ng), jroot(ng)
6
integer iownd3, imax, last, idum3
7
double precision alpha, x2, rdum3
8
common /lsr001/ alpha, x2, rdum3(3),
9
1 iownd3(3), imax, last, idum3(4)
10
c-----------------------------------------------------------------------
11
c this subroutine finds the leftmost root of a set of arbitrary
12
c functions gi(x) (i = 1,...,ng) in an interval (x0,x1). only roots
13
c of odd multiplicity (i.e. changes of sign of the gi) are found.
14
c here the sign of x1 - x0 is arbitrary, but is constant for a given
15
c problem, and -leftmost- means nearest to x0.
16
c the values of the vector-valued function g(x) = (gi, i=1...ng)
17
c are communicated through the call sequence of roots.
18
c the method used is the illinois algorithm.
21
c kathie l. hiebert and lawrence f. shampine, implicitly defined
22
c output points for solutions of ode-s, sandia report sand80-0180,
25
c description of parameters.
27
c ng = number of functions gi, or the number of components of
28
c the vector valued function g(x). input only.
30
c hmin = resolution parameter in x. input only. when a root is
31
c found, it is located only to within an error of hmin in x.
32
c typically, hmin should be set to something on the order of
33
c 100 * uround * max(abs(x0),abs(x1)),
34
c where uround is the unit roundoff of the machine.
36
c jflag = integer flag for input and output communication.
38
c on input, set jflag = 0 on the first call for the problem,
39
c and leave it unchanged until the problem is completed.
40
c (the problem is completed when jflag .ge. 2 on return.)
42
c on output, jflag has the following values and meanings..
43
c jflag = 1 means roots needs a value of g(x). set gx = g(x)
44
c and call roots again.
45
c jflag = 2 means a root has been found. the root is
46
c at x, and gx contains g(x). (actually, x is the
47
c rightmost approximation to the root on an interval
48
c (x0,x1) of size hmin or less.)
49
c jflag = 3 means x = x1 is a root, with one or more of the gi
50
c being zero at x1 and no sign changes in (x0,x1).
51
c gx contains g(x) on output.
52
c jflag = 4 means no roots (of odd multiplicity) were
53
c found in (x0,x1) (no sign changes).
55
c x0,x1 = endpoints of the interval where roots are sought.
56
c x1 and x0 are input when jflag = 0 (first call), and
57
c must be left unchanged between calls until the problem is
58
c completed. x0 and x1 must be distinct, but x1 - x0 may be
59
c of either sign. however, the notion of -left- and -right-
60
c will be used to mean nearer to x0 or x1, respectively.
61
c when jflag .ge. 2 on return, x0 and x1 are output, and
62
c are the endpoints of the relevant interval.
64
c g0,g1 = arrays of length ng containing the vectors g(x0) and g(x1),
65
c respectively. when jflag = 0, g0 and g1 are input and
66
c none of the g0(i) should be be zero.
67
c when jflag .ge. 2 on return, g0 and g1 are output.
69
c gx = array of length ng containing g(x). gx is input
70
c when jflag = 1, and output when jflag .ge. 2.
72
c x = independent variable value. output only.
73
c when jflag = 1 on output, x is the point at which g(x)
74
c is to be evaluated and loaded into gx.
75
c when jflag = 2 or 3, x is the root.
76
c when jflag = 4, x is the right endpoint of the interval, x1.
78
c jroot = integer array of length ng. output only.
79
c when jflag = 2 or 3, jroot indicates which components
80
c of g(x) have a root at x. jroot(i) is 1 if the i-th
81
c component has a root, and jroot(i) = 0 otherwise.
83
c note.. this routine uses the common block /lsr001/ to save
84
c the values of certain variables between calls (own variables).
85
c-----------------------------------------------------------------------
86
integer i, imxold, nxlast
87
double precision t2, tmax, zero
88
logical zroot, sgnchg, xroot
91
if (jflag .eq. 1) go to 200
92
c jflag .ne. 1. check for change in sign of g or zero at x1. ----------
97
if (dabs(g1(i)) .gt. zero) go to 110
100
c at this point, g0(i) has been checked and cannot be zero. ------------
101
110 if (dsign(1.0d0,g0(i)) .eq. dsign(1.0d0,g1(i))) go to 120
102
t2 = dabs(g1(i)/(g1(i)-g0(i)))
103
if (t2 .le. tmax) go to 120
107
if (imax .gt. 0) go to 130
111
140 if (.not. sgnchg) go to 400
112
c there is a sign change. find the first root in the interval. --------
117
c repeat until the first root in the interval is found. loop point. ---
120
if (nxlast .eq. last) go to 160
123
160 if (last .eq. 0) go to 170
126
170 alpha = 2.0d0*alpha
127
180 x2 = x1 - (x1-x0)*g1(imax)/(g1(imax) - alpha*g0(imax))
128
if ((dabs(x2-x0) .lt. hmin) .and.
129
1 (dabs(x1-x0) .gt. 10.0d0*hmin)) x2 = x0 + 0.1d0*(x1-x0)
132
c return to the calling routine to get a value of gx = g(x). -----------
134
c check to see in which interval g changes sign. -----------------------
140
if (dabs(gx(i)) .gt. zero) go to 210
143
c neither g0(i) nor gx(i) can be zero at this point. -------------------
144
210 if (dsign(1.0d0,g0(i)) .eq. dsign(1.0d0,gx(i))) go to 220
145
t2 = dabs(gx(i)/(gx(i) - g0(i)))
146
if (t2 .le. tmax) go to 220
150
if (imax .gt. 0) go to 230
156
if (.not. sgnchg) go to 250
157
c sign change between x0 and x2, so replace x1 with x2. ----------------
159
call dcopy (ng, gx, 1, g1, 1)
163
250 if (.not. zroot) go to 260
164
c zero value at x2 and no sign change in (x0,x2), so x2 is a root. -----
166
call dcopy (ng, gx, 1, g1, 1)
169
c no sign change between x0 and x2. replace x0 with x2. ---------------
171
call dcopy (ng, gx, 1, g0, 1)
175
270 if (dabs(x1-x0) .le. hmin) xroot = .true.
178
c return with x1 as the root. set jroot. set x = x1 and gx = g1. -----
181
call dcopy (ng, g1, 1, gx, 1)
184
if (dabs(g1(i)) .gt. zero) go to 310
187
310 if (dsign(1.0d0,g0(i)) .ne. dsign(1.0d0,g1(i))) jroot(i) = 1
191
c no sign change in the interval. check for zero at right endpoint. ---
192
400 if (.not. zroot) go to 420
194
c zero value at x1 and no sign change in (x0,x1). return jflag = 3. ---
196
call dcopy (ng, g1, 1, gx, 1)
199
if (dabs(g1(i)) .le. zero) jroot (i) = 1
204
c no sign changes in this interval. set x = x1, return jflag = 4. -----
205
420 call dcopy (ng, g1, 1, gx, 1)
209
c----------------------- end of subroutine roots -----------------------
added
removed
Lib/integrate/odepack/roots.f
