GNU Octave  3.8.0
A high-level interpreted language, primarily intended for numerical computations, mostly compatible with Matlab
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
sstode.f
Go to the documentation of this file.
1  SUBROUTINE sstode (NEQ, Y, YH, NYH, YH1, EWT, SAVF, ACOR,
2  1 wm, iwm, f, jac, pjac, slvs)
3 C***BEGIN PROLOGUE SSTODE
4 C***SUBSIDIARY
5 C***PURPOSE Performs one step of an ODEPACK integration.
6 C***TYPE SINGLE PRECISION (SSTODE-S, DSTODE-D)
7 C***AUTHOR Hindmarsh, Alan C., (LLNL)
8 C***DESCRIPTION
9 C
10 C SSTODE performs one step of the integration of an initial value
11 C problem for a system of ordinary differential equations.
12 C Note: SSTODE is independent of the value of the iteration method
13 C indicator MITER, when this is .ne. 0, and hence is independent
14 C of the type of chord method used, or the Jacobian structure.
15 C Communication with SSTODE is done with the following variables:
16 C
17 C NEQ = integer array containing problem size in NEQ(1), and
18 C passed as the NEQ argument in all calls to F and JAC.
19 C Y = an array of length .ge. N used as the Y argument in
20 C all calls to F and JAC.
21 C YH = an NYH by LMAX array containing the dependent variables
22 C and their approximate scaled derivatives, where
23 C LMAX = MAXORD + 1. YH(i,j+1) contains the approximate
24 C j-th derivative of y(i), scaled by h**j/factorial(j)
25 C (j = 0,1,...,NQ). on entry for the first step, the first
26 C two columns of YH must be set from the initial values.
27 C NYH = a constant integer .ge. N, the first dimension of YH.
28 C YH1 = a one-dimensional array occupying the same space as YH.
29 C EWT = an array of length N containing multiplicative weights
30 C for local error measurements. Local errors in Y(i) are
31 C compared to 1.0/EWT(i) in various error tests.
32 C SAVF = an array of working storage, of length N.
33 C Also used for input of YH(*,MAXORD+2) when JSTART = -1
34 C and MAXORD .lt. the current order NQ.
35 C ACOR = a work array of length N, used for the accumulated
36 C corrections. On a successful return, ACOR(i) contains
37 C the estimated one-step local error in Y(i).
38 C WM,IWM = real and integer work arrays associated with matrix
39 C operations in chord iteration (MITER .ne. 0).
40 C PJAC = name of routine to evaluate and preprocess Jacobian matrix
41 C and P = I - h*el0*JAC, if a chord method is being used.
42 C SLVS = name of routine to solve linear system in chord iteration.
43 C CCMAX = maximum relative change in h*el0 before PJAC is called.
44 C H = the step size to be attempted on the next step.
45 C H is altered by the error control algorithm during the
46 C problem. H can be either positive or negative, but its
47 C sign must remain constant throughout the problem.
48 C HMIN = the minimum absolute value of the step size h to be used.
49 C HMXI = inverse of the maximum absolute value of h to be used.
50 C HMXI = 0.0 is allowed and corresponds to an infinite hmax.
51 C HMIN and HMXI may be changed at any time, but will not
52 C take effect until the next change of h is considered.
53 C TN = the independent variable. TN is updated on each step taken.
54 C JSTART = an integer used for input only, with the following
55 C values and meanings:
56 C 0 perform the first step.
57 C .gt.0 take a new step continuing from the last.
58 C -1 take the next step with a new value of H, MAXORD,
59 C N, METH, MITER, and/or matrix parameters.
60 C -2 take the next step with a new value of H,
61 C but with other inputs unchanged.
62 C On return, JSTART is set to 1 to facilitate continuation.
63 C KFLAG = a completion code with the following meanings:
64 C 0 the step was succesful.
65 C -1 the requested error could not be achieved.
66 C -2 corrector convergence could not be achieved.
67 C -3 fatal error in PJAC or SLVS.
68 C A return with KFLAG = -1 or -2 means either
69 C abs(H) = HMIN or 10 consecutive failures occurred.
70 C On a return with KFLAG negative, the values of TN and
71 C the YH array are as of the beginning of the last
72 C step, and H is the last step size attempted.
73 C MAXORD = the maximum order of integration method to be allowed.
74 C MAXCOR = the maximum number of corrector iterations allowed.
75 C MSBP = maximum number of steps between PJAC calls (MITER .gt. 0).
76 C MXNCF = maximum number of convergence failures allowed.
77 C METH/MITER = the method flags. See description in driver.
78 C N = the number of first-order differential equations.
79 C The values of CCMAX, H, HMIN, HMXI, TN, JSTART, KFLAG, MAXORD,
80 C MAXCOR, MSBP, MXNCF, METH, MITER, and N are communicated via COMMON.
81 C
82 C***SEE ALSO SLSODE
83 C***ROUTINES CALLED SCFODE, SVNORM
84 C***COMMON BLOCKS SLS001
85 C***REVISION HISTORY (YYMMDD)
86 C 791129 DATE WRITTEN
87 C 890501 Modified prologue to SLATEC/LDOC format. (FNF)
88 C 890503 Minor cosmetic changes. (FNF)
89 C 930809 Renamed to allow single/double precision versions. (ACH)
90 C 010413 Reduced size of Common block /SLS001/. (ACH)
91 C 031105 Restored 'own' variables to Common block /SLS001/, to
92 C enable interrupt/restart feature. (ACH)
93 C***END PROLOGUE SSTODE
94 C**End
95  EXTERNAL f, jac, pjac, slvs
96  INTEGER neq, nyh, iwm
97  REAL y, yh, yh1, ewt, savf, acor, wm
98  dimension neq(*), y(*), yh(nyh,*), yh1(*), ewt(*), savf(*),
99  1 acor(*), wm(*), iwm(*)
100  INTEGER iownd, ialth, ipup, lmax, meo, nqnyh, nslp,
101  1 icf, ierpj, iersl, jcur, jstart, kflag, l,
102  2 lyh, lewt, lacor, lsavf, lwm, liwm, meth, miter,
103  3 maxord, maxcor, msbp, mxncf, n, nq, nst, nfe, nje, nqu
104  INTEGER i, i1, iredo, iret, j, jb, m, ncf, newq
105  REAL conit, crate, el, elco, hold, rmax, tesco,
106  2 ccmax, el0, h, hmin, hmxi, hu, rc, tn, uround
107  REAL dcon, ddn, del, delp, dsm, dup, exdn, exsm, exup,
108  1 r, rh, rhdn, rhsm, rhup, told, svnorm
109  COMMON /sls001/ conit, crate, el(13), elco(13,12),
110  1 hold, rmax, tesco(3,12),
111  2 ccmax, el0, h, hmin, hmxi, hu, rc, tn, uround,
112  3 iownd(6), ialth, ipup, lmax, meo, nqnyh, nslp,
113  3 icf, ierpj, iersl, jcur, jstart, kflag, l,
114  4 lyh, lewt, lacor, lsavf, lwm, liwm, meth, miter,
115  5 maxord, maxcor, msbp, mxncf, n, nq, nst, nfe, nje, nqu
116 C
117 C***FIRST EXECUTABLE STATEMENT SSTODE
118  kflag = 0
119  told = tn
120  ncf = 0
121  ierpj = 0
122  iersl = 0
123  jcur = 0
124  icf = 0
125  delp = 0.0e0
126  IF (jstart .GT. 0) go to 200
127  IF (jstart .EQ. -1) go to 100
128  IF (jstart .EQ. -2) go to 160
129 C-----------------------------------------------------------------------
130 C On the first call, the order is set to 1, and other variables are
131 C initialized. RMAX is the maximum ratio by which H can be increased
132 C in a single step. It is initially 1.E4 to compensate for the small
133 C initial H, but then is normally equal to 10. If a failure
134 C occurs (in corrector convergence or error test), RMAX is set to 2
135 C for the next increase.
136 C-----------------------------------------------------------------------
137  lmax = maxord + 1
138  nq = 1
139  l = 2
140  ialth = 2
141  rmax = 10000.0e0
142  rc = 0.0e0
143  el0 = 1.0e0
144  crate = 0.7e0
145  hold = h
146  meo = meth
147  nslp = 0
148  ipup = miter
149  iret = 3
150  go to 140
151 C-----------------------------------------------------------------------
152 C The following block handles preliminaries needed when JSTART = -1.
153 C IPUP is set to MITER to force a matrix update.
154 C If an order increase is about to be considered (IALTH = 1),
155 C IALTH is reset to 2 to postpone consideration one more step.
156 C If the caller has changed METH, SCFODE is called to reset
157 C the coefficients of the method.
158 C If the caller has changed MAXORD to a value less than the current
159 C order NQ, NQ is reduced to MAXORD, and a new H chosen accordingly.
160 C If H is to be changed, YH must be rescaled.
161 C If H or METH is being changed, IALTH is reset to L = NQ + 1
162 C to prevent further changes in H for that many steps.
163 C-----------------------------------------------------------------------
164  100 ipup = miter
165  lmax = maxord + 1
166  IF (ialth .EQ. 1) ialth = 2
167  IF (meth .EQ. meo) go to 110
168  CALL scfode(meth, elco, tesco)
169  meo = meth
170  IF (nq .GT. maxord) go to 120
171  ialth = l
172  iret = 1
173  go to 150
174  110 IF (nq .LE. maxord) go to 160
175  120 nq = maxord
176  l = lmax
177  DO 125 i = 1,l
178  125 el(i) = elco(i,nq)
179  nqnyh = nq*nyh
180  rc = rc*el(1)/el0
181  el0 = el(1)
182  conit = 0.5e0/(nq+2)
183  ddn = svnorm(n, savf, ewt)/tesco(1,l)
184  exdn = 1.0e0/l
185  rhdn = 1.0e0/(1.3e0*ddn**exdn + 0.0000013e0)
186  rh = min(rhdn,1.0e0)
187  iredo = 3
188  IF (h .EQ. hold) go to 170
189  rh = min(rh,abs(h/hold))
190  h = hold
191  go to 175
192 C-----------------------------------------------------------------------
193 C SCFODE is called to get all the integration coefficients for the
194 C current METH. Then the EL vector and related constants are reset
195 C whenever the order NQ is changed, or at the start of the problem.
196 C-----------------------------------------------------------------------
197  140 CALL scfode(meth, elco, tesco)
198  150 DO 155 i = 1,l
199  155 el(i) = elco(i,nq)
200  nqnyh = nq*nyh
201  rc = rc*el(1)/el0
202  el0 = el(1)
203  conit = 0.5e0/(nq+2)
204  go to(160, 170, 200), iret
205 C-----------------------------------------------------------------------
206 C If H is being changed, the H ratio RH is checked against
207 C RMAX, HMIN, and HMXI, and the YH array rescaled. IALTH is set to
208 C L = NQ + 1 to prevent a change of H for that many steps, unless
209 C forced by a convergence or error test failure.
210 C-----------------------------------------------------------------------
211  160 IF (h .EQ. hold) go to 200
212  rh = h/hold
213  h = hold
214  iredo = 3
215  go to 175
216  170 rh = max(rh,hmin/abs(h))
217  175 rh = min(rh,rmax)
218  rh = rh/max(1.0e0,abs(h)*hmxi*rh)
219  r = 1.0e0
220  DO 180 j = 2,l
221  r = r*rh
222  DO 180 i = 1,n
223  180 yh(i,j) = yh(i,j)*r
224  h = h*rh
225  rc = rc*rh
226  ialth = l
227  IF (iredo .EQ. 0) go to 690
228 C-----------------------------------------------------------------------
229 C This section computes the predicted values by effectively
230 C multiplying the YH array by the Pascal Triangle matrix.
231 C RC is the ratio of new to old values of the coefficient H*EL(1).
232 C When RC differs from 1 by more than CCMAX, IPUP is set to MITER
233 C to force PJAC to be called, if a Jacobian is involved.
234 C In any case, PJAC is called at least every MSBP steps.
235 C-----------------------------------------------------------------------
236  200 IF (abs(rc-1.0e0) .GT. ccmax) ipup = miter
237  IF (nst .GE. nslp+msbp) ipup = miter
238  tn = tn + h
239  i1 = nqnyh + 1
240  DO 215 jb = 1,nq
241  i1 = i1 - nyh
242 Cdir$ ivdep
243  DO 210 i = i1,nqnyh
244  210 yh1(i) = yh1(i) + yh1(i+nyh)
245  215 CONTINUE
246 C-----------------------------------------------------------------------
247 C Up to MAXCOR corrector iterations are taken. A convergence test is
248 C made on the R.M.S. norm of each correction, weighted by the error
249 C weight vector EWT. The sum of the corrections is accumulated in the
250 C vector ACOR(i). The YH array is not altered in the corrector loop.
251 C-----------------------------------------------------------------------
252  220 m = 0
253  DO 230 i = 1,n
254  230 y(i) = yh(i,1)
255  CALL f(neq, tn, y, savf)
256  nfe = nfe + 1
257  IF (ipup .LE. 0) go to 250
258 C-----------------------------------------------------------------------
259 C If indicated, the matrix P = I - h*el(1)*J is reevaluated and
260 C preprocessed before starting the corrector iteration. IPUP is set
261 C to 0 as an indicator that this has been done.
262 C-----------------------------------------------------------------------
263  CALL pjac(neq, y, yh, nyh, ewt, acor, savf, wm, iwm, f, jac)
264  ipup = 0
265  rc = 1.0e0
266  nslp = nst
267  crate = 0.7e0
268  IF (ierpj .NE. 0) go to 430
269  250 DO 260 i = 1,n
270  260 acor(i) = 0.0e0
271  270 IF (miter .NE. 0) go to 350
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  DO 290 i = 1,n
277  savf(i) = h*savf(i) - yh(i,2)
278  290 y(i) = savf(i) - acor(i)
279  del = svnorm(n, y, ewt)
280  DO 300 i = 1,n
281  y(i) = yh(i,1) + el(1)*savf(i)
282  300 acor(i) = savf(i)
283  go to 400
284 C-----------------------------------------------------------------------
285 C In the case of the chord method, compute the corrector error,
286 C and solve the linear system with that as right-hand side and
287 C P as coefficient matrix.
288 C-----------------------------------------------------------------------
289  350 DO 360 i = 1,n
290  360 y(i) = h*savf(i) - (yh(i,2) + acor(i))
291  CALL slvs(wm, iwm, y, savf)
292  IF (iersl .LT. 0) go to 430
293  IF (iersl .GT. 0) go to 410
294  del = svnorm(n, y, ewt)
295  DO 380 i = 1,n
296  acor(i) = acor(i) + y(i)
297  380 y(i) = yh(i,1) + el(1)*acor(i)
298 C-----------------------------------------------------------------------
299 C Test for convergence. If M.gt.0, an estimate of the convergence
300 C rate constant is stored in CRATE, and this is used in the test.
301 C-----------------------------------------------------------------------
302  400 IF (m .NE. 0) crate = max(0.2e0*crate,del/delp)
303  dcon = del*min(1.0e0,1.5e0*crate)/(tesco(2,nq)*conit)
304  IF (dcon .LE. 1.0e0) go to 450
305  m = m + 1
306  IF (m .EQ. maxcor) go to 410
307  IF (m .GE. 2 .AND. del .GT. 2.0e0*delp) go to 410
308  delp = del
309  CALL f(neq, tn, y, savf)
310  nfe = nfe + 1
311  go to 270
312 C-----------------------------------------------------------------------
313 C The corrector iteration failed to converge.
314 C If MITER .ne. 0 and the Jacobian is out of date, PJAC is called for
315 C the next try. Otherwise the YH array is retracted to its values
316 C before prediction, and H is reduced, if possible. If H cannot be
317 C reduced or MXNCF failures have occurred, exit with KFLAG = -2.
318 C-----------------------------------------------------------------------
319  410 IF (miter .EQ. 0 .OR. jcur .EQ. 1) go to 430
320  icf = 1
321  ipup = miter
322  go to 220
323  430 icf = 2
324  ncf = ncf + 1
325  rmax = 2.0e0
326  tn = told
327  i1 = nqnyh + 1
328  DO 445 jb = 1,nq
329  i1 = i1 - nyh
330 Cdir$ ivdep
331  DO 440 i = i1,nqnyh
332  440 yh1(i) = yh1(i) - yh1(i+nyh)
333  445 CONTINUE
334  IF (ierpj .LT. 0 .OR. iersl .LT. 0) go to 680
335  IF (abs(h) .LE. hmin*1.00001e0) go to 670
336  IF (ncf .EQ. mxncf) go to 670
337  rh = 0.25e0
338  ipup = miter
339  iredo = 1
340  go to 170
341 C-----------------------------------------------------------------------
342 C The corrector has converged. JCUR is set to 0
343 C to signal that the Jacobian involved may need updating later.
344 C The local error test is made and control passes to statement 500
345 C if it fails.
346 C-----------------------------------------------------------------------
347  450 jcur = 0
348  IF (m .EQ. 0) dsm = del/tesco(2,nq)
349  IF (m .GT. 0) dsm = svnorm(n, acor, ewt)/tesco(2,nq)
350  IF (dsm .GT. 1.0e0) go to 500
351 C-----------------------------------------------------------------------
352 C After a successful step, update the YH array.
353 C Consider changing H if IALTH = 1. Otherwise decrease IALTH by 1.
354 C If IALTH is then 1 and NQ .lt. MAXORD, then ACOR is saved for
355 C use in a possible order increase on the next step.
356 C If a change in H is considered, an increase or decrease in order
357 C by one is considered also. A change in H is made only if it is by a
358 C factor of at least 1.1. If not, IALTH is set to 3 to prevent
359 C testing for that many steps.
360 C-----------------------------------------------------------------------
361  kflag = 0
362  iredo = 0
363  nst = nst + 1
364  hu = h
365  nqu = nq
366  DO 470 j = 1,l
367  DO 470 i = 1,n
368  470 yh(i,j) = yh(i,j) + el(j)*acor(i)
369  ialth = ialth - 1
370  IF (ialth .EQ. 0) go to 520
371  IF (ialth .GT. 1) go to 700
372  IF (l .EQ. lmax) go to 700
373  DO 490 i = 1,n
374  490 yh(i,lmax) = acor(i)
375  go to 700
376 C-----------------------------------------------------------------------
377 C The error test failed. KFLAG keeps track of multiple failures.
378 C Restore TN and the YH array to their previous values, and prepare
379 C to try the step again. Compute the optimum step size for this or
380 C one lower order. After 2 or more failures, H is forced to decrease
381 C by a factor of 0.2 or less.
382 C-----------------------------------------------------------------------
383  500 kflag = kflag - 1
384  tn = told
385  i1 = nqnyh + 1
386  DO 515 jb = 1,nq
387  i1 = i1 - nyh
388 Cdir$ ivdep
389  DO 510 i = i1,nqnyh
390  510 yh1(i) = yh1(i) - yh1(i+nyh)
391  515 CONTINUE
392  rmax = 2.0e0
393  IF (abs(h) .LE. hmin*1.00001e0) go to 660
394  IF (kflag .LE. -3) go to 640
395  iredo = 2
396  rhup = 0.0e0
397  go to 540
398 C-----------------------------------------------------------------------
399 C Regardless of the success or failure of the step, factors
400 C RHDN, RHSM, and RHUP are computed, by which H could be multiplied
401 C at order NQ - 1, order NQ, or order NQ + 1, respectively.
402 C In the case of failure, RHUP = 0.0 to avoid an order increase.
403 C The largest of these is determined and the new order chosen
404 C accordingly. If the order is to be increased, we compute one
405 C additional scaled derivative.
406 C-----------------------------------------------------------------------
407  520 rhup = 0.0e0
408  IF (l .EQ. lmax) go to 540
409  DO 530 i = 1,n
410  530 savf(i) = acor(i) - yh(i,lmax)
411  dup = svnorm(n, savf, ewt)/tesco(3,nq)
412  exup = 1.0e0/(l+1)
413  rhup = 1.0e0/(1.4e0*dup**exup + 0.0000014e0)
414  540 exsm = 1.0e0/l
415  rhsm = 1.0e0/(1.2e0*dsm**exsm + 0.0000012e0)
416  rhdn = 0.0e0
417  IF (nq .EQ. 1) go to 560
418  ddn = svnorm(n, yh(1,l), ewt)/tesco(1,nq)
419  exdn = 1.0e0/nq
420  rhdn = 1.0e0/(1.3e0*ddn**exdn + 0.0000013e0)
421  560 IF (rhsm .GE. rhup) go to 570
422  IF (rhup .GT. rhdn) go to 590
423  go to 580
424  570 IF (rhsm .LT. rhdn) go to 580
425  newq = nq
426  rh = rhsm
427  go to 620
428  580 newq = nq - 1
429  rh = rhdn
430  IF (kflag .LT. 0 .AND. rh .GT. 1.0e0) rh = 1.0e0
431  go to 620
432  590 newq = l
433  rh = rhup
434  IF (rh .LT. 1.1e0) go to 610
435  r = el(l)/l
436  DO 600 i = 1,n
437  600 yh(i,newq+1) = acor(i)*r
438  go to 630
439  610 ialth = 3
440  go to 700
441  620 IF ((kflag .EQ. 0) .AND. (rh .LT. 1.1e0)) go to 610
442  IF (kflag .LE. -2) rh = min(rh,0.2e0)
443 C-----------------------------------------------------------------------
444 C If there is a change of order, reset NQ, l, and the coefficients.
445 C In any case H is reset according to RH and the YH array is rescaled.
446 C Then exit from 690 if the step was OK, or redo the step otherwise.
447 C-----------------------------------------------------------------------
448  IF (newq .EQ. nq) go to 170
449  630 nq = newq
450  l = nq + 1
451  iret = 2
452  go to 150
453 C-----------------------------------------------------------------------
454 C Control reaches this section if 3 or more failures have occured.
455 C If 10 failures have occurred, exit with KFLAG = -1.
456 C It is assumed that the derivatives that have accumulated in the
457 C YH array have errors of the wrong order. Hence the first
458 C derivative is recomputed, and the order is set to 1. Then
459 C H is reduced by a factor of 10, and the step is retried,
460 C until it succeeds or H reaches HMIN.
461 C-----------------------------------------------------------------------
462  640 IF (kflag .EQ. -10) go to 660
463  rh = 0.1e0
464  rh = max(hmin/abs(h),rh)
465  h = h*rh
466  DO 645 i = 1,n
467  645 y(i) = yh(i,1)
468  CALL f(neq, tn, y, savf)
469  nfe = nfe + 1
470  DO 650 i = 1,n
471  650 yh(i,2) = h*savf(i)
472  ipup = miter
473  ialth = 5
474  IF (nq .EQ. 1) go to 200
475  nq = 1
476  l = 2
477  iret = 3
478  go to 150
479 C-----------------------------------------------------------------------
480 C All returns are made through this section. H is saved in HOLD
481 C to allow the caller to change H on the next step.
482 C-----------------------------------------------------------------------
483  660 kflag = -1
484  go to 720
485  670 kflag = -2
486  go to 720
487  680 kflag = -3
488  go to 720
489  690 rmax = 10.0e0
490  700 r = 1.0e0/tesco(2,nqu)
491  DO 710 i = 1,n
492  710 acor(i) = acor(i)*r
493  720 hold = h
494  jstart = 1
495  RETURN
496 C----------------------- END OF SUBROUTINE SSTODE ----------------------
497  END