650323268f45b56f9eb005677a95556101688f81
[scilab.git] / scilab / modules / differential_equations / src / fortran / rchek2.f
1       subroutine rchek2(job, g, neq, y, yh, nyh, g0, g1, gx, jroot, irt
2      $     ,IWORK)
3 clll. optimize
4       include 'stack.h'
5       
6       external g
7       integer job, neq, nyh, jroot, irt
8       double precision y, yh, g0, g1, gx
9       dimension neq(*), y(*), yh(nyh,*), g0(*), g1(*), gx(*), jroot(*)
10       integer iownd, iowns,
11      1   icf, ierpj, iersl, jcur, jstart, kflag, l, meth, miter,
12      2   maxord, maxcor, msbp, mxncf, n, nq, nst, nfe, nje, nqu
13       integer iownd3, iownr3, irfnd, itaskc, ngc, nge
14       integer i, iflag, jflag
15       double precision rownd, rowns,
16      1   ccmax, el0, h, hmin, hmxi, hu, rc, tn, uround
17       double precision rownr3, t0, tlast, toutc
18       double precision hming, t1, x
19       logical zroot, Mroot
20       common /ls0001/ rownd, rowns(209),
21      2   ccmax, el0, h, hmin, hmxi, hu, rc, tn, uround,
22      3   iownd(14), iowns(6),
23      4   icf, ierpj, iersl, jcur, jstart, kflag, l, meth, miter,
24      5   maxord, maxcor, msbp, mxncf, n, nq, nst, nfe, nje, nqu
25       common /lsr001/ rownr3(2), t0, tlast, toutc,
26      1   iownd3(3), iownr3(2), irfnd, itaskc, ngc, nge      
27 c     ------------------ masking ----------------
28       integer IWORK
29       dimension IWORK(*)
30 c     pointers into iwork:
31       parameter (leniw=18)
32       integer    lmask
33       data zero/0.0D0/
34 c     ------------------ masking ----------------
35
36 c!purpose
37 c this routine checks for the presence of a root in the
38 c vicinity of the current t, in a manner depending on the
39 c input flag job.  it calls subroutine roots to locate the root
40 c as precisely as possible.
41 c
42 c!calling sequence
43 c in addition to variables described previously, rchek
44 c uses the following for communication..
45 c job    = integer flag indicating type of call..
46 c          job = 1 means the problem is being initialized, and rchek
47 c                  is to look for a root at or very near the initial t.
48 c          job = 2 means a continuation call to the solver was just
49 c                  made, and rchek is to check for a root in the
50 c                  relevant part of the step last taken.
51 c          job = 3 means a successful step was just taken, and rchek
52 c                  is to look for a root in the interval of the step.
53 c g0     = array of length ng, containing the value of g at t = t0.
54 c          g0 is input for job .ge. 2 and on output in all cases.
55 c g1,gx  = arrays of length ng for work space.
56 c irt    = completion flag..
57 c          irt = 0  means no root was found.
58 c          irt = -1 means job = 1 and a root was found too near to t.
59 c          irt = 1  means a legitimate root was found (job = 2 or 3).
60 c                   on return, t0 is the root location, and y is the
61 c                   corresponding solution vector.
62 c          irt = 2 means a change from zero to a non-zero value has been
63 c          occurred, so do a cold restart to reevaluate the modes
64 c          of if-then-else, because they have mode
65 c t0     = value of t at one endpoint of interval of interest.  only
66 c          roots beyond t0 in the direction of integration are sought.
67 c          t0 is input if job .ge. 2, and output in all cases.
68 c          t0 is updated by rchek, whether a root is found or not.
69 c tlast  = last value of t returned by the solver (input only).
70 c toutc  = copy of tout (input only).
71 c irfnd  = input flag showing whether the last step taken had a root.
72 c          irfnd = 1 if it did, = 0 if not.
73 c itaskc = copy of itask (input only).
74 c ngc    = copy of ng (input only).
75 c!
76 c
77       irt = 0
78 c     -------------- masking obtaining the mask adresses-----------------
79       lmask=iwork(leniw)-ngc
80 c     -------------- masking -----------------
81       hming = (dabs(tn) + dabs(h))*uround*100.0d0
82 c
83       go to (100, 200, 300), job
84 c
85 c evaluate g at initial t, and check for zero values. ------------------
86  100  continue
87 c     -------------- masking: disabling masks in cold-major-time-step------
88       DO 103 I = 1,ngc
89          jroot(i) = 0
90  103     iwork(lmask+i)=0
91
92       t0=tn
93       call g (neq, t0, y, ngc, g0)
94       nge = nge + 1
95
96       DO 110 I = 1,ngc
97          IF (DABS(g0(I)) .EQ. ZERO) THEN
98             iwork(lmask+i)=1
99          ENDIF
100  110  CONTINUE
101       RETURN
102 c     -------------- masking -----------------
103 c
104  200  continue
105
106 c     in the previous call there was not a root, so this part can be ignored.
107 c      IF (IWORK(LIRFND) .EQ. 0) GO TO 260
108        DO 203 I = 1,ngc
109           JROOT(I) = 0
110  203      IWORK(LMASK+I)=0
111 C     If a root was found on the previous step, evaluate R0 = R(T0). -------
112        call intdy (t0, 0, yh, nyh, y, iflag)
113        call g (neq, t0, y, ngc, g0)
114        nge = nge + 1
115        DO 210 I = 1,ngc
116           IF (dABS(g0(I)) .EQ. ZERO) THEN
117              IWORK(LMASK+I)=1
118           ENDIF
119  210   CONTINUE
120 C     R0 has no zero components.  Proceed to check relevant interval. ------
121  260   if (tn .eq. tlast) return
122 c
123 c     -------------- try in manor-time-steps -----
124  300  continue
125
126 c set t1 to tn or toutc, whichever comes first, and get g at t1. -------
127       if (itaskc.eq.2 .or. itaskc.eq.3 .or. itaskc.eq.5) go to 310
128       if ((toutc - tn)*h .ge. 0.0d0) go to 310
129       t1 = toutc
130       if ((t1 - t0)*h .le. 0.0d0) go to 390
131       call intdy (t1, 0, yh, nyh, y, iflag)
132       go to 330
133  310  t1 = tn
134       do 320 i = 1,n
135  320    y(i) = yh(i,1)
136  330  call g (neq, t1, y, ngc, g1)
137       if(ierror.gt.0) return
138       nge = nge + 1
139
140 C     Call DROOTS to search for root in interval from T0 to T1. -----------
141       JFLAG = 0
142       
143       DO 340 I = 1,Ngc
144          JROOT(I)=IWORK(LMASK+I)
145  340  CONTINUE
146       
147  350  CONTINUE
148       call roots2(ngc,hming, jflag, t0, t1, g0, g1, gx, x, jroot)
149       IF (JFLAG .GT. 1) GO TO 360
150       call intdy (x, 0, yh, nyh, y, iflag)
151       call g (neq, x, y, ngc, gx)
152       if(ierror.gt.0) return
153       nge = nge + 1
154       GO TO 350
155       
156  360  CONTINUE
157       if (JFLAG.eq.2) then      ! root found         
158          ZROOT=.false.
159          MROOT=.false.
160          DO 361 I = 1,Ngc          
161             if(IWORK(LMASK+I).eq.1) then
162                if(ABS(g1(i)).ne. ZERO) THEN
163                   JROOT(I)=SIGN(2.0D0,g1(I))
164                   Mroot=.true.
165                ELSE
166                   JROOT(I)=0
167                ENDIF
168             ELSE
169                IF (ABS(g1(I)) .EQ. ZERO) THEN
170                   JROOT(I) = -SIGN(1.0D0,g0(I))
171                   zroot=.true.
172                ELSE
173                   IF (SIGN(1.0D0,g0(I)) .NE. SIGN(1.0D0,g1(I))) THEN
174                      JROOT(I) = SIGN(1.0D0,g1(I) - g0(I))
175                      zroot=.true.
176                   ELSE
177                      JROOT(I)=0
178                   ENDIF
179                ENDIF
180             ENDIF
181  361     CONTINUE
182          
183          call intdy (x, 0, yh, nyh, y, iflag)
184
185          if (Zroot) then
186             DO 380 I = 1,Ngc
187                IF(ABS(JROOT(I)).EQ.2) JROOT(I)=0
188  380        CONTINUE  
189             MROOT=.false.
190             IRT=1
191          endif
192          IF (MROOT) THEN
193             IRT=2
194          ENDIF
195       ENDIF
196       T0 = X
197       CALL DCOPY (Ngc, gx, 1, g0, 1)
198       RETURN
199 c     
200  390  continue
201       return
202 c----------------------- end of subroutine rchek -----------------------
203       end