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