fix import/export of ierode common on Windows
[scilab.git] / scilab / modules / differential_equations / src / fortran / stode.f
1 C/MEMBR ADD NAME=STODE,SSI=0
2       subroutine stode (neq, y, yh, nyh, yh1, ewt, savf, acor,
3      1   wm, iwm, f, jac, pjac, slvs)
4 clll. optimize
5       
6       include 'stack.h'
7       
8       external f, jac, pjac, slvs
9       integer neq, nyh, iwm
10       integer iownd, ialth, ipup, lmax, meo, nqnyh, nslp,
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 i, i1, iredo, iret, j, jb, m, ncf, newq
14       double precision y, yh, yh1, ewt, savf, acor, wm
15       double precision rownd,
16      1   conit, crate, el, elco, hold, rmax, tesco,
17      2   ccmax, el0, h, hmin, hmxi, hu, rc, tn, uround
18       double precision dcon, ddn, del, delp, dsm, dup, exdn, exsm, exup,
19      1   r, rh, rhdn, rhsm, rhup, told, vnorm
20       dimension neq(*), y(*), yh(nyh,*), yh1(*), ewt(*), savf(*),
21      1   acor(*), wm(*), iwm(*)
22       common /ls0001/ rownd, conit, crate, el(13), elco(13,12),
23      1   hold, rmax, tesco(3,12),
24      2   ccmax, el0, h, hmin, hmxi, hu, rc, tn, uround, iownd(14),
25      3   ialth, ipup, lmax, meo, nqnyh, nslp,
26      4   icf, ierpj, iersl, jcur, jstart, kflag, l, meth, miter,
27      5   maxord, maxcor, msbp, mxncf, n, nq, nst, nfe, nje, nqu
28 c-----------------------------------------------------------------------
29 c%purpose
30 c stode performs one step of the integration of an initial value
31 c problem for a system of ordinary differential equations.
32 c note.. stode is independent of the value of the iteration method
33 c indicator miter, when this is .ne. 0, and hence is independent
34 c of the type of chord method used, or the jacobian structure.
35 c%calling sequence
36 c communication with stode is done with the following variables..
37 c
38 c neq    = integer array containing problem size in neq(1), and
39 c          passed as the neq argument in all calls to f and jac.
40 c y      = an array of length .ge. n used as the y argument in
41 c          all calls to f and jac.
42 c yh     = an nyh by lmax array containing the dependent variables
43 c          and their approximate scaled derivatives, where
44 c          lmax = maxord + 1.  yh(i,j+1) contains the approximate
45 c          j-th derivative of y(i), scaled by h**j/factorial(j)
46 c          (j = 0,1,...,nq).  on entry for the first step, the first
47 c          two columns of yh must be set from the initial values.
48 c nyh    = a constant integer .ge. n, the first dimension of yh.
49 c yh1    = a one-dimensional array occupying the same space as yh.
50 c ewt    = an array of length n containing multiplicative weights
51 c          for local error measurements.  local errors in y(i) are
52 c          compared to 1.0/ewt(i) in various error tests.
53 c savf   = an array of working storage, of length n.
54 c          also used for input of yh(*,maxord+2) when jstart = -1
55 c          and maxord .lt. the current order nq.
56 c acor   = a work array of length n, used for the accumulated
57 c          corrections.  on a successful return, acor(i) contains
58 c          the estimated one-step local error in y(i).
59 c wm,iwm = real and integer work arrays associated with matrix
60 c          operations in chord iteration (miter .ne. 0).
61 c pjac   = name of routine to evaluate and preprocess jacobian matrix
62 c          and p = i - h*el0*jac, if a chord method is being used.
63 c slvs   = name of routine to solve linear system in chord iteration.
64 c ccmax  = maximum relative change in h*el0 before pjac is called.
65 c h      = the step size to be attempted on the next step.
66 c          h is altered by the error control algorithm during the
67 c          problem.  h can be either positive or negative, but its
68 c          sign must remain constant throughout the problem.
69 c hmin   = the minimum absolute value of the step size h to be used.
70 c hmxi   = inverse of the maximum absolute value of h to be used.
71 c          hmxi = 0.0 is allowed and corresponds to an infinite hmax.
72 c          hmin and hmxi may be changed at any time, but will not
73 c          take effect until the next change of h is considered.
74 c tn     = the independent variable. tn is updated on each step taken.
75 c jstart = an integer used for input only, with the following
76 c          values and meanings..
77 c               0  perform the first step.
78 c           .gt.0  take a new step continuing from the last.
79 c              -1  take the next step with a new value of h, maxord,
80 c                    n, meth, miter, and/or matrix parameters.
81 c              -2  take the next step with a new value of h,
82 c                    but with other inputs unchanged.
83 c          on return, jstart is set to 1 to facilitate continuation.
84 c kflag  = a completion code with the following meanings..
85 c               0  the step was succesful.
86 c              -1  the requested error could not be achieved.
87 c              -2  corrector convergence could not be achieved.
88 c              -3  fatal error in pjac or slvs.
89 c          a return with kflag = -1 or -2 means either
90 c          abs(h) = hmin or 10 consecutive failures occurred.
91 c          on a return with kflag negative, the values of tn and
92 c          the yh array are as of the beginning of the last
93 c          step, and h is the last step size attempted.
94 c maxord = the maximum order of integration method to be allowed.
95 c maxcor = the maximum number of corrector iterations allowed.
96 c msbp   = maximum number of steps between pjac calls (miter .gt. 0).
97 c mxncf  = maximum number of convergence failures allowed.
98 c meth/miter = the method flags.  see description in driver.
99 c n      = the number of first-order differential equations.
100 c!
101 c-----------------------------------------------------------------------
102       save/ls0001/
103       kflag = 0
104       told = tn
105       ncf = 0
106       ierpj = 0
107       iersl = 0
108       jcur = 0
109       icf = 0
110       if (jstart .gt. 0) go to 200
111       if (jstart .eq. -1) go to 100
112       if (jstart .eq. -2) go to 160
113 c-----------------------------------------------------------------------
114 c on the first call, the order is set to 1, and other variables are
115 c initialized.  rmax is the maximum ratio by which h can be increased
116 c in a single step.  it is initially 1.e4 to compensate for the small
117 c initial h, but then is normally equal to 10.  if a failure
118 c occurs (in corrector convergence or error test), rmax is set at 2
119 c for the next increase.
120 c-----------------------------------------------------------------------
121       lmax = maxord + 1
122       nq = 1
123       l = 2
124       ialth = 2
125       rmax = 10000.0d+0
126       rc = 0.0d+0
127       el0 = 1.0d+0
128       crate = 0.70d+0
129       delp = 0.0d+0
130       hold = h
131       meo = meth
132       nslp = 0
133       ipup = miter
134       iret = 3
135       go to 140
136 c-----------------------------------------------------------------------
137 c the following block handles preliminaries needed when jstart = -1.
138 c ipup is set to miter to force a matrix update.
139 c if an order increase is about to be considered (ialth = 1),
140 c ialth is reset to 2 to postpone consideration one more step.
141 c if the caller has changed meth, cfode is called to reset
142 c the coefficients of the method.
143 c if the caller has changed maxord to a value less than the current
144 c order nq, nq is reduced to maxord, and a new h chosen accordingly.
145 c if h is to be changed, yh must be rescaled.
146 c if h or meth is being changed, ialth is reset to l = nq + 1
147 c to prevent further changes in h for that many steps.
148 c-----------------------------------------------------------------------
149  100  ipup = miter
150       lmax = maxord + 1
151       if (ialth .eq. 1) ialth = 2
152       if (meth .eq. meo) go to 110
153       call cfode (meth, elco(1,1), tesco(1,1))
154       meo = meth
155       if (nq .gt. maxord) go to 120
156       ialth = l
157       iret = 1
158       go to 150
159  110  if (nq .le. maxord) go to 160
160  120  nq = maxord
161       l = lmax
162       do 125 i = 1,l
163  125    el(i) = elco(i,nq)
164       nqnyh = nq*nyh
165       rc = rc*el(1)/el0
166       el0 = el(1)
167       conit = 0.50d+0/dble(nq+2)
168       ddn = vnorm (n, savf, ewt)/tesco(1,l)
169       exdn = 1.0d+0/dble(l)
170       rhdn = 1.0d+0/(1.30d+0*ddn**exdn + 0.00000130d+0)
171       rh = min(rhdn,1.0d+0)
172       iredo = 3
173       if (h .eq. hold) go to 170
174       rh = min(rh,abs(h/hold))
175       h = hold
176       go to 175
177 c-----------------------------------------------------------------------
178 c cfode is called to get all the integration coefficients for the
179 c current meth.  then the el vector and related constants are reset
180 c whenever the order nq is changed, or at the start of the problem.
181 c-----------------------------------------------------------------------
182  140  call cfode (meth, elco(1,1), tesco(1,1))
183  150  do 155 i = 1,l
184  155    el(i) = elco(i,nq)
185       nqnyh = nq*nyh
186       rc = rc*el(1)/el0
187       el0 = el(1)
188       conit = 0.50d+0/dble(nq+2)
189       go to (160, 170, 200), iret
190 c-----------------------------------------------------------------------
191 c if h is being changed, the h ratio rh is checked against
192 c rmax, hmin, and hmxi, and the yh array rescaled.  ialth is set to
193 c l = nq + 1 to prevent a change of h for that many steps, unless
194 c forced by a convergence or error test failure.
195 c-----------------------------------------------------------------------
196  160  if (h .eq. hold) go to 200
197       rh = h/hold
198       h = hold
199       iredo = 3
200       go to 175
201  170  rh = max(rh,hmin/abs(h))
202  175  rh = min(rh,rmax)
203       rh = rh/max(1.0d+0,abs(h)*hmxi*rh)
204       r = 1.0d+0
205       do 180 j = 2,l
206         r = r*rh
207         do 180 i = 1,n
208  180      yh(i,j) = yh(i,j)*r
209       h = h*rh
210       rc = rc*rh
211       ialth = l
212       if (iredo .eq. 0) go to 690
213 c-----------------------------------------------------------------------
214 c this section computes the predicted values by effectively
215 c multiplying the yh array by the pascal triangle matrix.
216 c rc is the ratio of new to old values of the coefficient  h*el(1).
217 c when rc differs from 1 by more than ccmax, ipup is set to miter
218 c to force pjac to be called, if a jacobian is involved.
219 c in any case, pjac is called at least every msbp steps.
220 c-----------------------------------------------------------------------
221  200  if (abs(rc-1.0d+0) .gt. ccmax) ipup = miter
222       if (nst .ge. nslp+msbp) ipup = miter
223       tn = tn + h
224       i1 = nqnyh + 1
225       do 215 jb = 1,nq
226         i1 = i1 - nyh
227         do 210 i = i1,nqnyh
228  210      yh1(i) = yh1(i) + yh1(i+nyh)
229  215    continue
230 c-----------------------------------------------------------------------
231 c up to maxcor corrector iterations are taken.  a convergence test is
232 c made on the r.m.s. norm of each correction, weighted by the error
233 c weight vector ewt.  the sum of the corrections is accumulated in the
234 c vector acor(i).  the yh array is not altered in the corrector loop.
235 c-----------------------------------------------------------------------
236  220  m = 0
237       do 230 i = 1,n
238  230    y(i) = yh(i,1)
239       call f (neq, tn, y, savf)
240       if(ierror.gt.0) return
241       nfe = nfe + 1
242       if (ipup .le. 0) go to 250
243 c-----------------------------------------------------------------------
244 c if indicated, the matrix p = i - h*el(1)*j is reevaluated and
245 c preprocessed before starting the corrector iteration.  ipup is set
246 c to 0 as an indicator that this has been done.
247 c-----------------------------------------------------------------------
248       ipup = 0
249       rc = 1.0d+0
250       nslp = nst
251       crate = 0.70d+0
252       call pjac (neq, y, yh, nyh, ewt, acor, savf, wm, iwm, f, jac)
253       if(ierror.gt.0) return
254       if (ierpj .ne. 0) go to 430
255  250  do 260 i = 1,n
256  260    acor(i) = 0.0d+0
257  270  if (miter .ne. 0) go to 350
258 c-----------------------------------------------------------------------
259 c in the case of functional iteration, update y directly from
260 c the result of the last function evaluation.
261 c-----------------------------------------------------------------------
262       do 290 i = 1,n
263         savf(i) = h*savf(i) - yh(i,2)
264  290    y(i) = savf(i) - acor(i)
265       del = vnorm (n, y, ewt)
266       do 300 i = 1,n
267         y(i) = yh(i,1) + el(1)*savf(i)
268  300    acor(i) = savf(i)
269       go to 400
270 c-----------------------------------------------------------------------
271 c in the case of the chord method, compute the corrector error,
272 c and solve the linear system with that as right-hand side and
273 c p as coefficient matrix.
274 c-----------------------------------------------------------------------
275  350  do 360 i = 1,n
276  360    y(i) = h*savf(i) - (yh(i,2) + acor(i))
277       call slvs (wm, iwm, y, savf)
278       if (iersl .lt. 0) go to 430
279       if (iersl .gt. 0) go to 410
280       del = vnorm (n, y, ewt)
281       do 380 i = 1,n
282         acor(i) = acor(i) + y(i)
283  380    y(i) = yh(i,1) + el(1)*acor(i)
284 c-----------------------------------------------------------------------
285 c test for convergence.  if m.gt.0, an estimate of the convergence
286 c rate constant is stored in crate, and this is used in the test.
287 c-----------------------------------------------------------------------
288  400  if (m .ne. 0) crate = max(0.20d+0*crate,del/delp)
289       dcon = del*min(1.0d+0,1.50d+0*crate)/(tesco(2,nq)*conit)
290       if (dcon .le. 1.0d+0) go to 450
291       m = m + 1
292       if (m .eq. maxcor) go to 410
293       if (m .ge. 2 .and. del .gt. 2.0d+0*delp) go to 410
294       delp = del
295       call f (neq, tn, y, savf)
296       if(ierror.gt.0) return
297       nfe = nfe + 1
298       go to 270
299 c-----------------------------------------------------------------------
300 c the corrector iteration failed to converge in maxcor tries.
301 c if miter .ne. 0 and the jacobian is out of date, pjac is called for
302 c the next try.  otherwise the yh array is retracted to its values
303 c before prediction, and h is reduced, if possible.  if h cannot be
304 c reduced or mxncf failures have occurred, exit with kflag = -2.
305 c-----------------------------------------------------------------------
306  410  if (miter .eq. 0 .or. jcur .eq. 1) go to 430
307       icf = 1
308       ipup = miter
309       go to 220
310  430  icf = 2
311       ncf = ncf + 1
312       rmax = 2.0d+0
313       tn = told
314       i1 = nqnyh + 1
315       do 445 jb = 1,nq
316         i1 = i1 - nyh
317         do 440 i = i1,nqnyh
318  440      yh1(i) = yh1(i) - yh1(i+nyh)
319  445    continue
320       if (ierpj .lt. 0 .or. iersl .lt. 0) go to 680
321       if (abs(h) .le. hmin*1.000010d+0) go to 670
322       if (ncf .eq. mxncf) go to 670
323       rh = 0.250d+0
324       ipup = miter
325       iredo = 1
326       go to 170
327 c-----------------------------------------------------------------------
328 c the corrector has converged.  jcur is set to 0
329 c to signal that the jacobian involved may need updating later.
330 c the local error test is made and control passes to statement 500
331 c if it fails.
332 c-----------------------------------------------------------------------
333  450  jcur = 0
334       if (m .eq. 0) dsm = del/tesco(2,nq)
335       if (m .gt. 0) dsm = vnorm (n, acor, ewt)/tesco(2,nq)
336       if (dsm .gt. 1.0d+0) go to 500
337 c-----------------------------------------------------------------------
338 c after a successful step, update the yh array.
339 c consider changing h if ialth = 1.  otherwise decrease ialth by 1.
340 c if ialth is then 1 and nq .lt. maxord, then acor is saved for
341 c use in a possible order increase on the next step.
342 c if a change in h is considered, an increase or decrease in order
343 c by one is considered also.  a change in h is made only if it is by a
344 c factor of at least 1.1.  if not, ialth is set to 3 to prevent
345 c testing for that many steps.
346 c-----------------------------------------------------------------------
347       kflag = 0
348       iredo = 0
349       nst = nst + 1
350       hu = h
351       nqu = nq
352       do 470 j = 1,l
353         do 470 i = 1,n
354  470      yh(i,j) = yh(i,j) + el(j)*acor(i)
355       ialth = ialth - 1
356       if (ialth .eq. 0) go to 520
357       if (ialth .gt. 1) go to 700
358       if (l .eq. lmax) go to 700
359       do 490 i = 1,n
360  490    yh(i,lmax) = acor(i)
361       go to 700
362 c-----------------------------------------------------------------------
363 c the error test failed.  kflag keeps track of multiple failures.
364 c restore tn and the yh array to their previous values, and prepare
365 c to try the step again.  compute the optimum step size for this or
366 c one lower order.  after 2 or more failures, h is forced to decrease
367 c by a factor of 0.2 or less.
368 c-----------------------------------------------------------------------
369  500  kflag = kflag - 1
370       tn = told
371       i1 = nqnyh + 1
372       do 515 jb = 1,nq
373         i1 = i1 - nyh
374         do 510 i = i1,nqnyh
375  510      yh1(i) = yh1(i) - yh1(i+nyh)
376  515    continue
377       rmax = 2.0d+0
378       if (abs(h) .le. hmin*1.000010d+0) go to 660
379       if (kflag .le. -3) go to 640
380       iredo = 2
381       rhup = 0.0d+0
382       go to 540
383 c-----------------------------------------------------------------------
384 c regardless of the success or failure of the step, factors
385 c rhdn, rhsm, and rhup are computed, by which h could be multiplied
386 c at order nq - 1, order nq, or order nq + 1, respectively.
387 c in the case of failure, rhup = 0.0 to avoid an order increase.
388 c the largest of these is determined and the new order chosen
389 c accordingly.  if the order is to be increased, we compute one
390 c additional scaled derivative.
391 c-----------------------------------------------------------------------
392  520  rhup = 0.0d+0
393       if (l .eq. lmax) go to 540
394       do 530 i = 1,n
395  530    savf(i) = acor(i) - yh(i,lmax)
396       dup = vnorm (n, savf, ewt)/tesco(3,nq)
397       exup = 1.0d+0/dble(l+1)
398       rhup = 1.0d+0/(1.40d+0*dup**exup + 0.00000140d+0)
399  540  exsm = 1.0d+0/dble(l)
400       rhsm = 1.0d+0/(1.20d+0*dsm**exsm + 0.00000120d+0)
401       rhdn = 0.0d+0
402       if (nq .eq. 1) go to 560
403       ddn = vnorm (n, yh(1,l), ewt)/tesco(1,nq)
404       exdn = 1.0d+0/dble(nq)
405       rhdn = 1.0d+0/(1.30d+0*ddn**exdn + 0.00000130d+0)
406  560  if (rhsm .ge. rhup) go to 570
407       if (rhup .gt. rhdn) go to 590
408       go to 580
409  570  if (rhsm .lt. rhdn) go to 580
410       newq = nq
411       rh = rhsm
412       go to 620
413  580  newq = nq - 1
414       rh = rhdn
415       if (kflag .lt. 0 .and. rh .gt. 1.0d+0) rh = 1.0d+0
416       go to 620
417  590  newq = l
418       rh = rhup
419       if (rh .lt. 1.10d+0) go to 610
420       r = el(l)/dble(l)
421       do 600 i = 1,n
422  600    yh(i,newq+1) = acor(i)*r
423       go to 630
424  610  ialth = 3
425       go to 700
426  620  if ((kflag .eq. 0) .and. (rh .lt. 1.10d+0)) go to 610
427       if (kflag .le. -2) rh = min(rh,0.20d+0)
428 c-----------------------------------------------------------------------
429 c if there is a change of order, reset nq, l, and the coefficients.
430 c in any case h is reset according to rh and the yh array is rescaled.
431 c then exit from 690 if the step was ok, or redo the step otherwise.
432 c-----------------------------------------------------------------------
433       if (newq .eq. nq) go to 170
434  630  nq = newq
435       l = nq + 1
436       iret = 2
437       go to 150
438 c-----------------------------------------------------------------------
439 c control reaches this section if 3 or more failures have occurred.
440 c if 10 failures have occurred, exit with kflag = -1.
441 c it is assumed that the derivatives that have accumulated in the
442 c yh array have errors of the wrong order.  hence the first
443 c derivative is recomputed, and the order is set to 1.  then
444 c h is reduced by a factor of 10, and the step is retried,
445 c until it succeeds or h reaches hmin.
446 c-----------------------------------------------------------------------
447  640  if (kflag .eq. -10) go to 660
448       rh = 0.10d+0
449       rh = max(hmin/abs(h),rh)
450       h = h*rh
451       do 645 i = 1,n
452  645    y(i) = yh(i,1)
453       call f (neq, tn, y, savf)
454       if(ierror.gt.0) return
455       nfe = nfe + 1
456       do 650 i = 1,n
457  650    yh(i,2) = h*savf(i)
458       ipup = miter
459       ialth = 5
460       if (nq .eq. 1) go to 200
461       nq = 1
462       l = 2
463       iret = 3
464       go to 150
465 c-----------------------------------------------------------------------
466 c all returns are made through this section.  h is saved in hold
467 c to allow the caller to change h on the next step.
468 c-----------------------------------------------------------------------
469  660  kflag = -1
470       go to 720
471  670  kflag = -2
472       go to 720
473  680  kflag = -3
474       go to 720
475  690  rmax = 10.0d+0
476  700  r = 1.0d+0/tesco(2,nqu)
477       do 710 i = 1,n
478  710    acor(i) = acor(i)*r
479  720  hold = h
480       jstart = 1
481       return
482 c----------------------- end of subroutine stode -----------------------
483       end