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