GNU Octave  4.4.1
A high-level interpreted language, primarily intended for numerical computations, mostly compatible with Matlab
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 INIT, MXSTEP, MXHNIL, NHNIL, NSLAST, CNYH,
101  1 ialth, ipup, lmax, meo, nqnyh, nslp,
102  1 icf, ierpj, iersl, jcur, jstart, kflag, l,
103  2 lyh, lewt, lacor, lsavf, lwm, liwm, meth, miter,
104  3 maxord, maxcor, msbp, mxncf, n, nq, nst, nfe, nje, nqu
105  INTEGER I, I1, IREDO, IRET, J, JB, M, NCF, NEWQ
106  REAL CONIT, CRATE, EL, ELCO, HOLD, RMAX, TESCO,
107  1 ccmax, el0, h, hmin, hmxi, hu, rc, tn, uround
108  REAL DCON, DDN, DEL, DELP, DSM, DUP, EXDN, EXSM, EXUP,
109  1 r, rh, rhdn, rhsm, rhup, told, svnorm
110  COMMON /sls001/ conit, crate, el(13), elco(13,12),
111  1 hold, rmax, tesco(3,12),
112  1 ccmax, el0, h, hmin, hmxi, hu, rc, tn, uround,
113  2 init, mxstep, mxhnil, nhnil, nslast, cnyh,
114  3 ialth, ipup, lmax, meo, nqnyh, nslp,
115  3 icf, ierpj, iersl, jcur, jstart, kflag, l,
116  4 lyh, lewt, lacor, lsavf, lwm, liwm, meth, miter,
117  5 maxord, maxcor, msbp, mxncf, n, nq, nst, nfe, nje, nqu
118 C
119 C***FIRST EXECUTABLE STATEMENT SSTODE
120  kflag = 0
121  told = tn
122  ncf = 0
123  ierpj = 0
124  iersl = 0
125  jcur = 0
126  icf = 0
127  delp = 0.0e0
128  IF (jstart .GT. 0) GO TO 200
129  IF (jstart .EQ. -1) GO TO 100
130  IF (jstart .EQ. -2) GO TO 160
131 C-----------------------------------------------------------------------
132 C On the first call, the order is set to 1, and other variables are
133 C initialized. RMAX is the maximum ratio by which H can be increased
134 C in a single step. It is initially 1.E4 to compensate for the small
135 C initial H, but then is normally equal to 10. If a failure
136 C occurs (in corrector convergence or error test), RMAX is set to 2
137 C for the next increase.
138 C-----------------------------------------------------------------------
139  lmax = maxord + 1
140  nq = 1
141  l = 2
142  ialth = 2
143  rmax = 10000.0e0
144  rc = 0.0e0
145  el0 = 1.0e0
146  crate = 0.7e0
147  hold = h
148  meo = meth
149  nslp = 0
150  ipup = miter
151  iret = 3
152  GO TO 140
153 C-----------------------------------------------------------------------
154 C The following block handles preliminaries needed when JSTART = -1.
155 C IPUP is set to MITER to force a matrix update.
156 C If an order increase is about to be considered (IALTH = 1),
157 C IALTH is reset to 2 to postpone consideration one more step.
158 C If the caller has changed METH, SCFODE is called to reset
159 C the coefficients of the method.
160 C If the caller has changed MAXORD to a value less than the current
161 C order NQ, NQ is reduced to MAXORD, and a new H chosen accordingly.
162 C If H is to be changed, YH must be rescaled.
163 C If H or METH is being changed, IALTH is reset to L = NQ + 1
164 C to prevent further changes in H for that many steps.
165 C-----------------------------------------------------------------------
166  100 ipup = miter
167  lmax = maxord + 1
168  IF (ialth .EQ. 1) ialth = 2
169  IF (meth .EQ. meo) GO TO 110
170  CALL scfode (meth, elco, tesco)
171  meo = meth
172  IF (nq .GT. maxord) GO TO 120
173  ialth = l
174  iret = 1
175  GO TO 150
176  110 IF (nq .LE. maxord) GO TO 160
177  120 nq = maxord
178  l = lmax
179  DO 125 i = 1,l
180  125 el(i) = elco(i,nq)
181  nqnyh = nq*nyh
182  rc = rc*el(1)/el0
183  el0 = el(1)
184  conit = 0.5e0/(nq+2)
185  ddn = svnorm(n, savf, ewt)/tesco(1,l)
186  exdn = 1.0e0/l
187  rhdn = 1.0e0/(1.3e0*ddn**exdn + 0.0000013e0)
188  rh = min(rhdn,1.0e0)
189  iredo = 3
190  IF (h .EQ. hold) GO TO 170
191  rh = min(rh,abs(h/hold))
192  h = hold
193  GO TO 175
194 C-----------------------------------------------------------------------
195 C SCFODE is called to get all the integration coefficients for the
196 C current METH. Then the EL vector and related constants are reset
197 C whenever the order NQ is changed, or at the start of the problem.
198 C-----------------------------------------------------------------------
199  140 CALL scfode (meth, elco, tesco)
200  150 DO 155 i = 1,l
201  155 el(i) = elco(i,nq)
202  nqnyh = nq*nyh
203  rc = rc*el(1)/el0
204  el0 = el(1)
205  conit = 0.5e0/(nq+2)
206  GO TO (160, 170, 200), iret
207 C-----------------------------------------------------------------------
208 C If H is being changed, the H ratio RH is checked against
209 C RMAX, HMIN, and HMXI, and the YH array rescaled. IALTH is set to
210 C L = NQ + 1 to prevent a change of H for that many steps, unless
211 C forced by a convergence or error test failure.
212 C-----------------------------------------------------------------------
213  160 IF (h .EQ. hold) GO TO 200
214  rh = h/hold
215  h = hold
216  iredo = 3
217  GO TO 175
218  170 rh = max(rh,hmin/abs(h))
219  175 rh = min(rh,rmax)
220  rh = rh/max(1.0e0,abs(h)*hmxi*rh)
221  r = 1.0e0
222  DO 180 j = 2,l
223  r = r*rh
224  DO 180 i = 1,n
225  180 yh(i,j) = yh(i,j)*r
226  h = h*rh
227  rc = rc*rh
228  ialth = l
229  IF (iredo .EQ. 0) GO TO 690
230 C-----------------------------------------------------------------------
231 C This section computes the predicted values by effectively
232 C multiplying the YH array by the Pascal Triangle matrix.
233 C RC is the ratio of new to old values of the coefficient H*EL(1).
234 C When RC differs from 1 by more than CCMAX, IPUP is set to MITER
235 C to force PJAC to be called, if a Jacobian is involved.
236 C In any case, PJAC is called at least every MSBP steps.
237 C-----------------------------------------------------------------------
238  200 IF (abs(rc-1.0e0) .GT. ccmax) ipup = miter
239  IF (nst .GE. nslp+msbp) ipup = miter
240  tn = tn + h
241  i1 = nqnyh + 1
242  DO 215 jb = 1,nq
243  i1 = i1 - nyh
244 Cdir$ ivdep
245  DO 210 i = i1,nqnyh
246  210 yh1(i) = yh1(i) + yh1(i+nyh)
247  215 CONTINUE
248 C-----------------------------------------------------------------------
249 C Up to MAXCOR corrector iterations are taken. A convergence test is
250 C made on the R.M.S. norm of each correction, weighted by the error
251 C weight vector EWT. The sum of the corrections is accumulated in the
252 C vector ACOR(i). The YH array is not altered in the corrector loop.
253 C-----------------------------------------------------------------------
254  220 m = 0
255  DO 230 i = 1,n
256  230 y(i) = yh(i,1)
257  CALL f (neq, tn, y, savf)
258  nfe = nfe + 1
259  IF (ipup .LE. 0) GO TO 250
260 C-----------------------------------------------------------------------
261 C If indicated, the matrix P = I - h*el(1)*J is reevaluated and
262 C preprocessed before starting the corrector iteration. IPUP is set
263 C to 0 as an indicator that this has been done.
264 C-----------------------------------------------------------------------
265  CALL pjac (neq, y, yh, nyh, ewt, acor, savf, wm, iwm, f, jac)
266  ipup = 0
267  rc = 1.0e0
268  nslp = nst
269  crate = 0.7e0
270  IF (ierpj .NE. 0) GO TO 430
271  250 DO 260 i = 1,n
272  260 acor(i) = 0.0e0
273  270 IF (miter .NE. 0) GO TO 350
274 C-----------------------------------------------------------------------
275 C In the case of functional iteration, update Y directly from
276 C the result of the last function evaluation.
277 C-----------------------------------------------------------------------
278  DO 290 i = 1,n
279  savf(i) = h*savf(i) - yh(i,2)
280  290 y(i) = savf(i) - acor(i)
281  del = svnorm(n, y, ewt)
282  DO 300 i = 1,n
283  y(i) = yh(i,1) + el(1)*savf(i)
284  300 acor(i) = savf(i)
285  GO TO 400
286 C-----------------------------------------------------------------------
287 C In the case of the chord method, compute the corrector error,
288 C and solve the linear system with that as right-hand side and
289 C P as coefficient matrix.
290 C-----------------------------------------------------------------------
291  350 DO 360 i = 1,n
292  360 y(i) = h*savf(i) - (yh(i,2) + acor(i))
293  CALL slvs (wm, iwm, y, savf)
294  IF (iersl .LT. 0) GO TO 430
295  IF (iersl .GT. 0) GO TO 410
296  del = svnorm(n, y, ewt)
297  DO 380 i = 1,n
298  acor(i) = acor(i) + y(i)
299  380 y(i) = yh(i,1) + el(1)*acor(i)
300 C-----------------------------------------------------------------------
301 C Test for convergence. If M.gt.0, an estimate of the convergence
302 C rate constant is stored in CRATE, and this is used in the test.
303 C-----------------------------------------------------------------------
304  400 IF (m .NE. 0) crate = max(0.2e0*crate,del/delp)
305  dcon = del*min(1.0e0,1.5e0*crate)/(tesco(2,nq)*conit)
306  IF (dcon .LE. 1.0e0) GO TO 450
307  m = m + 1
308  IF (m .EQ. maxcor) GO TO 410
309  IF (m .GE. 2 .AND. del .GT. 2.0e0*delp) GO TO 410
310  delp = del
311  CALL f (neq, tn, y, savf)
312  nfe = nfe + 1
313  GO TO 270
314 C-----------------------------------------------------------------------
315 C The corrector iteration failed to converge.
316 C If MITER .ne. 0 and the Jacobian is out of date, PJAC is called for
317 C the next try. Otherwise the YH array is retracted to its values
318 C before prediction, and H is reduced, if possible. If H cannot be
319 C reduced or MXNCF failures have occurred, exit with KFLAG = -2.
320 C-----------------------------------------------------------------------
321  410 IF (miter .EQ. 0 .OR. jcur .EQ. 1) GO TO 430
322  icf = 1
323  ipup = miter
324  GO TO 220
325  430 icf = 2
326  ncf = ncf + 1
327  rmax = 2.0e0
328  tn = told
329  i1 = nqnyh + 1
330  DO 445 jb = 1,nq
331  i1 = i1 - nyh
332 Cdir$ ivdep
333  DO 440 i = i1,nqnyh
334  440 yh1(i) = yh1(i) - yh1(i+nyh)
335  445 CONTINUE
336  IF (ierpj .LT. 0 .OR. iersl .LT. 0) GO TO 680
337  IF (abs(h) .LE. hmin*1.00001e0) GO TO 670
338  IF (ncf .EQ. mxncf) GO TO 670
339  rh = 0.25e0
340  ipup = miter
341  iredo = 1
342  GO TO 170
343 C-----------------------------------------------------------------------
344 C The corrector has converged. JCUR is set to 0
345 C to signal that the Jacobian involved may need updating later.
346 C The local error test is made and control passes to statement 500
347 C if it fails.
348 C-----------------------------------------------------------------------
349  450 jcur = 0
350  IF (m .EQ. 0) dsm = del/tesco(2,nq)
351  IF (m .GT. 0) dsm = svnorm(n, acor, ewt)/tesco(2,nq)
352  IF (dsm .GT. 1.0e0) GO TO 500
353 C-----------------------------------------------------------------------
354 C After a successful step, update the YH array.
355 C Consider changing H if IALTH = 1. Otherwise decrease IALTH by 1.
356 C If IALTH is then 1 and NQ .lt. MAXORD, then ACOR is saved for
357 C use in a possible order increase on the next step.
358 C If a change in H is considered, an increase or decrease in order
359 C by one is considered also. A change in H is made only if it is by a
360 C factor of at least 1.1. If not, IALTH is set to 3 to prevent
361 C testing for that many steps.
362 C-----------------------------------------------------------------------
363  kflag = 0
364  iredo = 0
365  nst = nst + 1
366  hu = h
367  nqu = nq
368  DO 470 j = 1,l
369  DO 470 i = 1,n
370  470 yh(i,j) = yh(i,j) + el(j)*acor(i)
371  ialth = ialth - 1
372  IF (ialth .EQ. 0) GO TO 520
373  IF (ialth .GT. 1) GO TO 700
374  IF (l .EQ. lmax) GO TO 700
375  DO 490 i = 1,n
376  490 yh(i,lmax) = acor(i)
377  GO TO 700
378 C-----------------------------------------------------------------------
379 C The error test failed. KFLAG keeps track of multiple failures.
380 C Restore TN and the YH array to their previous values, and prepare
381 C to try the step again. Compute the optimum step size for this or
382 C one lower order. After 2 or more failures, H is forced to decrease
383 C by a factor of 0.2 or less.
384 C-----------------------------------------------------------------------
385  500 kflag = kflag - 1
386  tn = told
387  i1 = nqnyh + 1
388  DO 515 jb = 1,nq
389  i1 = i1 - nyh
390 Cdir$ ivdep
391  DO 510 i = i1,nqnyh
392  510 yh1(i) = yh1(i) - yh1(i+nyh)
393  515 CONTINUE
394  rmax = 2.0e0
395  IF (abs(h) .LE. hmin*1.00001e0) GO TO 660
396  IF (kflag .LE. -3) GO TO 640
397  iredo = 2
398  rhup = 0.0e0
399  GO TO 540
400 C-----------------------------------------------------------------------
401 C Regardless of the success or failure of the step, factors
402 C RHDN, RHSM, and RHUP are computed, by which H could be multiplied
403 C at order NQ - 1, order NQ, or order NQ + 1, respectively.
404 C In the case of failure, RHUP = 0.0 to avoid an order increase.
405 C The largest of these is determined and the new order chosen
406 C accordingly. If the order is to be increased, we compute one
407 C additional scaled derivative.
408 C-----------------------------------------------------------------------
409  520 rhup = 0.0e0
410  IF (l .EQ. lmax) GO TO 540
411  DO 530 i = 1,n
412  530 savf(i) = acor(i) - yh(i,lmax)
413  dup = svnorm(n, savf, ewt)/tesco(3,nq)
414  exup = 1.0e0/(l+1)
415  rhup = 1.0e0/(1.4e0*dup**exup + 0.0000014e0)
416  540 exsm = 1.0e0/l
417  rhsm = 1.0e0/(1.2e0*dsm**exsm + 0.0000012e0)
418  rhdn = 0.0e0
419  IF (nq .EQ. 1) GO TO 560
420  ddn = svnorm(n, yh(1,l), ewt)/tesco(1,nq)
421  exdn = 1.0e0/nq
422  rhdn = 1.0e0/(1.3e0*ddn**exdn + 0.0000013e0)
423  560 IF (rhsm .GE. rhup) GO TO 570
424  IF (rhup .GT. rhdn) GO TO 590
425  GO TO 580
426  570 IF (rhsm .LT. rhdn) GO TO 580
427  newq = nq
428  rh = rhsm
429  GO TO 620
430  580 newq = nq - 1
431  rh = rhdn
432  IF (kflag .LT. 0 .AND. rh .GT. 1.0e0) rh = 1.0e0
433  GO TO 620
434  590 newq = l
435  rh = rhup
436  IF (rh .LT. 1.1e0) GO TO 610
437  r = el(l)/l
438  DO 600 i = 1,n
439  600 yh(i,newq+1) = acor(i)*r
440  GO TO 630
441  610 ialth = 3
442  GO TO 700
443  620 IF ((kflag .EQ. 0) .AND. (rh .LT. 1.1e0)) GO TO 610
444  IF (kflag .LE. -2) rh = min(rh,0.2e0)
445 C-----------------------------------------------------------------------
446 C If there is a change of order, reset NQ, l, and the coefficients.
447 C In any case H is reset according to RH and the YH array is rescaled.
448 C Then exit from 690 if the step was OK, or redo the step otherwise.
449 C-----------------------------------------------------------------------
450  IF (newq .EQ. nq) GO TO 170
451  630 nq = newq
452  l = nq + 1
453  iret = 2
454  GO TO 150
455 C-----------------------------------------------------------------------
456 C Control reaches this section if 3 or more failures have occurred.
457 C If 10 failures have occurred, exit with KFLAG = -1.
458 C It is assumed that the derivatives that have accumulated in the
459 C YH array have errors of the wrong order. Hence the first
460 C derivative is recomputed, and the order is set to 1. Then
461 C H is reduced by a factor of 10, and the step is retried,
462 C until it succeeds or H reaches HMIN.
463 C-----------------------------------------------------------------------
464  640 IF (kflag .EQ. -10) GO TO 660
465  rh = 0.1e0
466  rh = max(hmin/abs(h),rh)
467  h = h*rh
468  DO 645 i = 1,n
469  645 y(i) = yh(i,1)
470  CALL f (neq, tn, y, savf)
471  nfe = nfe + 1
472  DO 650 i = 1,n
473  650 yh(i,2) = h*savf(i)
474  ipup = miter
475  ialth = 5
476  IF (nq .EQ. 1) GO TO 200
477  nq = 1
478  l = 2
479  iret = 3
480  GO TO 150
481 C-----------------------------------------------------------------------
482 C All returns are made through this section. H is saved in HOLD
483 C to allow the caller to change H on the next step.
484 C-----------------------------------------------------------------------
485  660 kflag = -1
486  GO TO 720
487  670 kflag = -2
488  GO TO 720
489  680 kflag = -3
490  GO TO 720
491  690 rmax = 10.0e0
492  700 r = 1.0e0/tesco(2,nqu)
493  DO 710 i = 1,n
494  710 acor(i) = acor(i)*r
495  720 hold = h
496  jstart = 1
497  RETURN
498 C----------------------- END OF SUBROUTINE SSTODE ----------------------
499  END
static T abs(T x)
Definition: pr-output.cc:1696
OCTAVE_EXPORT octave_value_list etc The functions then dimension(columns)
charNDArray max(char d, const charNDArray &m)
Definition: chNDArray.cc:227
charNDArray min(char d, const charNDArray &m)
Definition: chNDArray.cc:204