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