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