subroutine rchek (job, g, neq, y, yh, nyh, g0, g1, gx, jroot, irt)
clll. optimize
external g
integer job, neq, nyh, jroot, irt
double precision y, yh, g0, g1, gx
dimension neq(1), y(1), yh(nyh,1), g0(1), g1(1), gx(1), jroot(1)
integer iownd, iowns,
1 icf, ierpj, iersl, jcur, jstart, kflag, l, meth, miter,
2 maxord, maxcor, msbp, mxncf, n, nq, nst, nfe, nje, nqu
integer iownd3, iownr3, irfnd, itaskc, ngc, nge
integer i, iflag, jflag
double precision rownd, rowns,
1 ccmax, el0, h, hmin, hmxi, hu, rc, tn, uround
double precision rownr3, t0, tlast, toutc
double precision hming, t1, temp1, temp2, x
logical zroot
common /ls0001/ rownd, rowns(209),
2 ccmax, el0, h, hmin, hmxi, hu, rc, tn, uround,
3 iownd(14), iowns(6),
4 icf, ierpj, iersl, jcur, jstart, kflag, l, meth, miter,
5 maxord, maxcor, msbp, mxncf, n, nq, nst, nfe, nje, nqu
common /lsr001/ rownr3(2), t0, tlast, toutc,
1 iownd3(3), iownr3(2), irfnd, itaskc, ngc, nge
c-----------------------------------------------------------------------
c this routine checks for the presence of a root in the
c vicinity of the current t, in a manner depending on the
c input flag job. it calls subroutine roots to locate the root
c as precisely as possible.
c
c in addition to variables described previously, rchek
c uses the following for communication..
c job = integer flag indicating type of call..
c job = 1 means the problem is being initialized, and rchek
c is to look for a root at or very near the initial t.
c job = 2 means a continuation call to the solver was just
c made, and rchek is to check for a root in the
c relevant part of the step last taken.
c job = 3 means a successful step was just taken, and rchek
c is to look for a root in the interval of the step.
c g0 = array of length ng, containing the value of g at t = t0.
c g0 is input for job .ge. 2 and on output in all cases.
c g1,gx = arrays of length ng for work space.
c irt = completion flag..
c irt = 0 means no root was found.
c irt = -1 means job = 1 and a root was found too near to t.
c irt = 1 means a legitimate root was found (job = 2 or 3).
c on return, t0 is the root location, and y is the
c corresponding solution vector.
c t0 = value of t at one endpoint of interval of interest. only
c roots beyond t0 in the direction of integration are sought.
c t0 is input if job .ge. 2, and output in all cases.
c t0 is updated by rchek, whether a root is found or not.
c tlast = last value of t returned by the solver (input only).
c toutc = copy of tout (input only).
c irfnd = input flag showing whether the last step taken had a root.
c irfnd = 1 if it did, = 0 if not.
c itaskc = copy of itask (input only).
c ngc = copy of ng (input only).
c-----------------------------------------------------------------------
c
irt = 0
do 10 i = 1,ngc
10 jroot(i) = 0
hming = (dabs(tn) + dabs(h))*uround*100.0d0
c
go to (100, 200, 300), job
c
c evaluate g at initial t, and check for zero values. ------------------
100 continue
t0 = tn
call g (neq, t0, y, ngc, g0)
nge = 1
zroot = .false.
do 110 i = 1,ngc
110 if (dabs(g0(i)) .le. 0.0d0) zroot = .true.
if (.not. zroot) go to 190
c g has a zero at t. look at g at t + (small increment). --------------
temp1 = dsign(hming,h)
t0 = t0 + temp1
temp2 = temp1/h
do 120 i = 1,n
120 y(i) = y(i) + temp2*yh(i,2)
call g (neq, t0, y, ngc, g0)
nge = nge + 1
zroot = .false.
do 130 i = 1,ngc
130 if (dabs(g0(i)) .le. 0.0d0) zroot = .true.
if (.not. zroot) go to 190
c g has a zero at t and also close to t. take error return. -----------
irt = -1
return
c
190 continue
return
c
c
200 continue
if (irfnd .eq. 0) go to 260
c if a root was found on the previous step, evaluate g0 = g(t0). -------
call intdy (t0, 0, yh, nyh, y, iflag)
call g (neq, t0, y, ngc, g0)
nge = nge + 1
zroot = .false.
do 210 i = 1,ngc
210 if (dabs(g0(i)) .le. 0.0d0) zroot = .true.
if (.not. zroot) go to 260
c g has a zero at t0. look at g at t + (small increment). -------------
temp1 = dsign(hming,h)
t0 = t0 + temp1
if ((t0 - tn)*h .lt. 0.0d0) go to 230
temp2 = temp1/h
do 220 i = 1,n
220 y(i) = y(i) + temp2*yh(i,2)
go to 240
230 call intdy (t0, 0, yh, nyh, y, iflag)
240 call g (neq, t0, y, ngc, g0)
nge = nge + 1
zroot = .false.
do 250 i = 1,ngc
if (dabs(g0(i)) .gt. 0.0d0) go to 250
jroot(i) = 1
zroot = .true.
250 continue
if (.not. zroot) go to 260
c g has a zero at t0 and also close to t0. return root. ---------------
irt = 1
return
c here, g0 does not have a root
c g0 has no zero components. proceed to check relevant interval. ------
260 if (tn .eq. tlast) go to 390
c
300 continue
c set t1 to tn or toutc, whichever comes first, and get g at t1. -------
if (itaskc.eq.2 .or. itaskc.eq.3 .or. itaskc.eq.5) go to 310
if ((toutc - tn)*h .ge. 0.0d0) go to 310
t1 = toutc
if ((t1 - t0)*h .le. 0.0d0) go to 390
call intdy (t1, 0, yh, nyh, y, iflag)
go to 330
310 t1 = tn
do 320 i = 1,n
320 y(i) = yh(i,1)
330 call g (neq, t1, y, ngc, g1)
nge = nge + 1
c call roots to search for root in interval from t0 to t1. -------------
jflag = 0
350 continue
call roots (ngc, hming, jflag, t0, t1, g0, g1, gx, x, jroot)
if (jflag .gt. 1) go to 360
call intdy (x, 0, yh, nyh, y, iflag)
call g (neq, x, y, ngc, gx)
nge = nge + 1
go to 350
360 t0 = x
call dcopy (ngc, gx, 1, g0, 1)
if (jflag .eq. 4) go to 390
c found a root. interpolate to x and return. --------------------------
call intdy (x, 0, yh, nyh, y, iflag)
irt = 1
return
c
390 continue
return
c----------------------- end of subroutine rchek -----------------------
end