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