fix import/export of ierode common on Windows
[scilab.git] / scilab / modules / differential_equations / src / fortran / rchek.f
1       subroutine rchek (job, g, neq, y, yh, nyh, g0, g1, gx, jroot, irt)
2 clll. optimize
3       include 'stack.h'
4       
5       external g
6       integer job, neq, nyh, jroot, irt
7       double precision y, yh, g0, g1, gx
8       dimension neq(*), y(*), yh(nyh,*), g0(*), g1(*), gx(*), jroot(*)
9       integer iownd, iowns,
10      1   icf, ierpj, iersl, jcur, jstart, kflag, l, meth, miter,
11      2   maxord, maxcor, msbp, mxncf, n, nq, nst, nfe, nje, nqu
12       integer iownd3, iownr3, irfnd, itaskc, ngc, nge
13       integer i, iflag, jflag
14       double precision rownd, rowns,
15      1   ccmax, el0, h, hmin, hmxi, hu, rc, tn, uround
16       double precision rownr3, t0, tlast, toutc
17       double precision hming, t1, temp1, temp2, x
18       logical zroot
19       common /ls0001/ rownd, rowns(209),
20      2   ccmax, el0, h, hmin, hmxi, hu, rc, tn, uround,
21      3   iownd(14), iowns(6),
22      4   icf, ierpj, iersl, jcur, jstart, kflag, l, meth, miter,
23      5   maxord, maxcor, msbp, mxncf, n, nq, nst, nfe, nje, nqu
24       common /lsr001/ rownr3(2), t0, tlast, toutc,
25      1   iownd3(3), iownr3(2), irfnd, itaskc, ngc, nge      
26
27
28 c!purpose
29 c this routine checks for the presence of a root in the
30 c vicinity of the current t, in a manner depending on the
31 c input flag job.  it calls subroutine roots to locate the root
32 c as precisely as possible.
33 c
34 c!calling sequence
35 c in addition to variables described previously, rchek
36 c uses the following for communication..
37 c job    = integer flag indicating type of call..
38 c          job = 1 means the problem is being initialized, and rchek
39 c                  is to look for a root at or very near the initial t.
40 c          job = 2 means a continuation call to the solver was just
41 c                  made, and rchek is to check for a root in the
42 c                  relevant part of the step last taken.
43 c          job = 3 means a successful step was just taken, and rchek
44 c                  is to look for a root in the interval of the step.
45 c g0     = array of length ng, containing the value of g at t = t0.
46 c          g0 is input for job .ge. 2 and on output in all cases.
47 c g1,gx  = arrays of length ng for work space.
48 c irt    = completion flag..
49 c          irt = 0  means no root was found.
50 c          irt = -1 means job = 1 and a root was found too near to t.
51 c          irt = 1  means a legitimate root was found (job = 2 or 3).
52 c                   on return, t0 is the root location, and y is the
53 c                   corresponding solution vector.
54 c t0     = value of t at one endpoint of interval of interest.  only
55 c          roots beyond t0 in the direction of integration are sought.
56 c          t0 is input if job .ge. 2, and output in all cases.
57 c          t0 is updated by rchek, whether a root is found or not.
58 c tlast  = last value of t returned by the solver (input only).
59 c toutc  = copy of tout (input only).
60 c irfnd  = input flag showing whether the last step taken had a root.
61 c          irfnd = 1 if it did, = 0 if not.
62 c itaskc = copy of itask (input only).
63 c ngc    = copy of ng (input only).
64 c!
65 c
66       irt = 0
67       do 10 i = 1,ngc
68  10     jroot(i) = 0
69       hming = (dabs(tn) + dabs(h))*uround*100.0d0
70 c
71       go to (100, 200, 300), job
72 c
73 c 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
83 c 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
96 c g has a zero at t and also close to t.  take error return. -----------
97       irt = -1
98       return
99 c
100  190  continue
101       return
102 c
103 c
104  200  continue
105       if (irfnd .eq. 0) go to 260
106 c 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
119 c 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      
145 c g0 has no zero components.  proceed to check relevant interval. ------
146  260  if (tn .eq. tlast) go to 390
147 c
148  300  continue
149 c 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
162 c 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
175 c found a root.  interpolate to x and return. --------------------------
176       call intdy (x, 0, yh, nyh, y, iflag)
177       irt = 1
178       return
179 c
180  390  continue
181       return
182 c----------------------- end of subroutine rchek -----------------------
183       end