1 subroutine rchek (job, g, neq, y, yh, nyh, g0, g1, gx, jroot, irt) 2clll. optimize 3 external g 4 integer job, neq, nyh, jroot, irt 5 double precision y, yh, g0, g1, gx 6 dimension neq(*), y(*), yh(nyh,*), g0(*), g1(*), gx(*), jroot(*) 7 integer iownd, iowns, 8 1 icf, ierpj, iersl, jcur, jstart, kflag, l, meth, miter, 9 2 maxord, maxcor, msbp, mxncf, n, nq, nst, nfe, nje, nqu 10 integer iownd3, iownr3, irfnd, itaskc, ngc, nge 11 integer i, iflag, jflag 12 double precision rownd, rowns, 13 1 ccmax, el0, h, hmin, hmxi, hu, rc, tn, uround 14 double precision rownr3, t0, tlast, toutc 15 double precision hming, t1, temp1, temp2, x 16 logical zroot 17cDEC$ ATTRIBUTES DLLIMPORT:: /ls0001/ 18 common /ls0001/ rownd, rowns(209), 19 2 ccmax, el0, h, hmin, hmxi, hu, rc, tn, uround, 20 3 iownd(14), iowns(6), 21 4 icf, ierpj, iersl, jcur, jstart, kflag, l, meth, miter, 22 5 maxord, maxcor, msbp, mxncf, n, nq, nst, nfe, nje, nqu 23cDEC$ ATTRIBUTES DLLIMPORT:: /lsr001/ 24 common /lsr001/ rownr3(2), t0, tlast, toutc, 25 1 iownd3(3), iownr3(2), irfnd, itaskc, ngc, nge 26 27 28c!purpose 29c this routine checks for the presence of a root in the 30c vicinity of the current t, in a manner depending on the 31c input flag job. it calls subroutine roots to locate the root 32c as precisely as possible. 33c 34c!calling sequence 35c in addition to variables described previously, rchek 36c uses the following for communication.. 37c job = integer flag indicating type of call.. 38c job = 1 means the problem is being initialized, and rchek 39c is to look for a root at or very near the initial t. 40c job = 2 means a continuation call to the solver was just 41c made, and rchek is to check for a root in the 42c relevant part of the step last taken. 43c job = 3 means a successful step was just taken, and rchek 44c is to look for a root in the interval of the step. 45c g0 = array of length ng, containing the value of g at t = t0. 46c g0 is input for job .ge. 2 and on output in all cases. 47c g1,gx = arrays of length ng for work space. 48c irt = completion flag.. 49c irt = 0 means no root was found. 50c irt = -1 means job = 1 and a root was found too near to t. 51c irt = 1 means a legitimate root was found (job = 2 or 3). 52c on return, t0 is the root location, and y is the 53c corresponding solution vector. 54c t0 = value of t at one endpoint of interval of interest. only 55c roots beyond t0 in the direction of integration are sought. 56c t0 is input if job .ge. 2, and output in all cases. 57c t0 is updated by rchek, whether a root is found or not. 58c tlast = last value of t returned by the solver (input only). 59c toutc = copy of tout (input only). 60c irfnd = input flag showing whether the last step taken had a root. 61c irfnd = 1 if it did, = 0 if not. 62c itaskc = copy of itask (input only). 63c ngc = copy of ng (input only). 64c! 65c 66 irt = 0 67 do 10 i = 1,ngc 68 10 jroot(i) = 0 69 hming = (dabs(tn) + dabs(h))*uround*100.0d0 70c 71 go to (100, 200, 300), job 72c 73c evaluate g at initial t, and check for zero values. ------------------ 74 100 continue 75 t0 = tn 76 call g (neq, t0, y, ngc, g0) 77 if(ierror.gt.0) return 78 nge = 1 79 zroot = .false. 80 do 110 i = 1,ngc 81 110 if (dabs(g0(i)) .le. 0.0d0) zroot = .true. 82 if (.not. zroot) go to 190 83c g has a zero at t. look at g at t + (small increment). -------------- 84 temp1 = dsign(hming,h) 85 t0 = t0 + temp1 86 temp2 = temp1/h 87 do 120 i = 1,n 88 120 y(i) = y(i) + temp2*yh(i,2) 89 call g (neq, t0, y, ngc, g0) 90 if(ierror.gt.0) return 91 nge = nge + 1 92 zroot = .false. 93 do 130 i = 1,ngc 94 130 if (dabs(g0(i)) .le. 0.0d0) zroot = .true. 95 if (.not. zroot) go to 190 96c g has a zero at t and also close to t. take error return. ----------- 97 irt = -1 98 return 99c 100 190 continue 101 return 102c 103c 104 200 continue 105 if (irfnd .eq. 0) go to 260 106c if a root was found on the previous step, evaluate g0 = g(t0). ------- 107 call intdy (t0, 0, yh, nyh, y, iflag) 108 call g (neq, t0, y, ngc, g0) 109 if(ierror.gt.0) return 110 nge = nge + 1 111 zroot = .false. 112 do 210 i = 1,ngc 113 if (dabs(g0(i)) .le. 0.0d0) then 114 jroot(i)=1 115 zroot = .true. 116 endif 117 210 continue 118 if (.not. zroot) go to 260 119c g has a zero at t0. look at g at t + (small increment). ------------- 120 temp1 = dsign(hming,h) 121 t0 = t0 + temp1 122 if ((t0 - tn)*h .lt. 0.0d0) go to 230 123 temp2 = temp1/h 124 do 220 i = 1,n 125 220 y(i) = y(i) + temp2*yh(i,2) 126 go to 240 127 230 call intdy (t0, 0, yh, nyh, y, iflag) 128 240 call g (neq, t0, y, ngc, g0) 129 if(ierror.gt.0) return 130 nge = nge + 1 131 zroot = .false. 132 do 250 i = 1,ngc 133 if (dabs(g0(i)) .gt. 0.0d0) go to 250 134 135 if (jroot(i) .eq. 1) then 136 irt = -1 137 return 138 else 139 jroot(i) = -sign(1.0d0,g0(i)) 140 irt = 1 141 endif 142 143 250 continue 144 if (irt .eq. 1) return 145c g0 has no zero components. proceed to check relevant interval. ------ 146 260 if (tn .eq. tlast) go to 390 147c 148 300 continue 149c set t1 to tn or toutc, whichever comes first, and get g at t1. ------- 150 if (itaskc.eq.2 .or. itaskc.eq.3 .or. itaskc.eq.5) go to 310 151 if ((toutc - tn)*h .ge. 0.0d0) go to 310 152 t1 = toutc 153 if ((t1 - t0)*h .le. 0.0d0) go to 390 154 call intdy (t1, 0, yh, nyh, y, iflag) 155 go to 330 156 310 t1 = tn 157 do 320 i = 1,n 158 320 y(i) = yh(i,1) 159 330 call g (neq, t1, y, ngc, g1) 160 if(ierror.gt.0) return 161 nge = nge + 1 162c call roots to search for root in interval from t0 to t1. ------------- 163 jflag = 0 164 350 continue 165 call roots (ngc, hming, jflag, t0, t1, g0, g1, gx, x, jroot) 166 if (jflag .gt. 1) go to 360 167 call intdy (x, 0, yh, nyh, y, iflag) 168 call g (neq, x, y, ngc, gx) 169 if(ierror.gt.0) return 170 nge = nge + 1 171 go to 350 172 360 t0 = x 173 call dcopy (ngc, gx, 1, g0, 1) 174 if (jflag .eq. 4) go to 390 175c found a root. interpolate to x and return. -------------------------- 176 call intdy (x, 0, yh, nyh, y, iflag) 177 irt = 1 178 return 179c 180 390 continue 181 return 182c----------------------- end of subroutine rchek ----------------------- 183 end 184