1      SUBROUTINE SSTODE (NEQ, Y, YH, NYH, YH1, EWT, SAVF, ACOR,
2     1   WM, IWM, F, JAC, PJAC, SLVS)
3C***BEGIN PROLOGUE  SSTODE
4C***SUBSIDIARY
5C***PURPOSE  Performs one step of an ODEPACK integration.
6C***TYPE      SINGLE PRECISION (SSTODE-S, DSTODE-D)
7C***AUTHOR  Hindmarsh, Alan C., (LLNL)
8C***DESCRIPTION
9C
10C  SSTODE performs one step of the integration of an initial value
11C  problem for a system of ordinary differential equations.
12C  Note:  SSTODE is independent of the value of the iteration method
13C  indicator MITER, when this is .ne. 0, and hence is independent
14C  of the type of chord method used, or the Jacobian structure.
15C  Communication with SSTODE is done with the following variables:
16C
17C  NEQ    = integer array containing problem size in NEQ(1), and
18C           passed as the NEQ argument in all calls to F and JAC.
19C  Y      = an array of length .ge. N used as the Y argument in
20C           all calls to F and JAC.
21C  YH     = an NYH by LMAX array containing the dependent variables
22C           and their approximate scaled derivatives, where
23C           LMAX = MAXORD + 1.  YH(i,j+1) contains the approximate
24C           j-th derivative of y(i), scaled by h**j/factorial(j)
25C           (j = 0,1,...,NQ).  on entry for the first step, the first
26C           two columns of YH must be set from the initial values.
27C  NYH    = a constant integer .ge. N, the first dimension of YH.
28C  YH1    = a one-dimensional array occupying the same space as YH.
29C  EWT    = an array of length N containing multiplicative weights
30C           for local error measurements.  Local errors in Y(i) are
31C           compared to 1.0/EWT(i) in various error tests.
32C  SAVF   = an array of working storage, of length N.
33C           Also used for input of YH(*,MAXORD+2) when JSTART = -1
34C           and MAXORD .lt. the current order NQ.
35C  ACOR   = a work array of length N, used for the accumulated
36C           corrections.  On a successful return, ACOR(i) contains
37C           the estimated one-step local error in Y(i).
38C  WM,IWM = real and integer work arrays associated with matrix
39C           operations in chord iteration (MITER .ne. 0).
40C  PJAC   = name of routine to evaluate and preprocess Jacobian matrix
41C           and P = I - h*el0*JAC, if a chord method is being used.
42C  SLVS   = name of routine to solve linear system in chord iteration.
43C  CCMAX  = maximum relative change in h*el0 before PJAC is called.
44C  H      = the step size to be attempted on the next step.
45C           H is altered by the error control algorithm during the
46C           problem.  H can be either positive or negative, but its
47C           sign must remain constant throughout the problem.
48C  HMIN   = the minimum absolute value of the step size h to be used.
49C  HMXI   = inverse of the maximum absolute value of h to be used.
50C           HMXI = 0.0 is allowed and corresponds to an infinite hmax.
51C           HMIN and HMXI may be changed at any time, but will not
52C           take effect until the next change of h is considered.
53C  TN     = the independent variable. TN is updated on each step taken.
54C  JSTART = an integer used for input only, with the following
55C           values and meanings:
56C                0  perform the first step.
57C            .gt.0  take a new step continuing from the last.
58C               -1  take the next step with a new value of H, MAXORD,
59C                     N, METH, MITER, and/or matrix parameters.
60C               -2  take the next step with a new value of H,
61C                     but with other inputs unchanged.
62C           On return, JSTART is set to 1 to facilitate continuation.
63C  KFLAG  = a completion code with the following meanings:
64C                0  the step was succesful.
65C               -1  the requested error could not be achieved.
66C               -2  corrector convergence could not be achieved.
67C               -3  fatal error in PJAC or SLVS.
68C           A return with KFLAG = -1 or -2 means either
69C           abs(H) = HMIN or 10 consecutive failures occurred.
70C           On a return with KFLAG negative, the values of TN and
71C           the YH array are as of the beginning of the last
72C           step, and H is the last step size attempted.
73C  MAXORD = the maximum order of integration method to be allowed.
74C  MAXCOR = the maximum number of corrector iterations allowed.
75C  MSBP   = maximum number of steps between PJAC calls (MITER .gt. 0).
76C  MXNCF  = maximum number of convergence failures allowed.
77C  METH/MITER = the method flags.  See description in driver.
78C  N      = the number of first-order differential equations.
79C  The values of CCMAX, H, HMIN, HMXI, TN, JSTART, KFLAG, MAXORD,
80C  MAXCOR, MSBP, MXNCF, METH, MITER, and N are communicated via COMMON.
81C
82C***SEE ALSO  SLSODE
83C***ROUTINES CALLED  SCFODE, SVNORM
84C***COMMON BLOCKS    SLS001
85C***REVISION HISTORY  (YYMMDD)
86C   791129  DATE WRITTEN
87C   890501  Modified prologue to SLATEC/LDOC format.  (FNF)
88C   890503  Minor cosmetic changes.  (FNF)
89C   930809  Renamed to allow single/double precision versions. (ACH)
90C   010413  Reduced size of Common block /SLS001/. (ACH)
91C   031105  Restored 'own' variables to Common block /SLS001/, to
92C           enable interrupt/restart feature. (ACH)
93C***END PROLOGUE  SSTODE
94C**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
118C
119C***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
131C-----------------------------------------------------------------------
132C On the first call, the order is set to 1, and other variables are
133C initialized.  RMAX is the maximum ratio by which H can be increased
134C in a single step.  It is initially 1.E4 to compensate for the small
135C initial H, but then is normally equal to 10.  If a failure
136C occurs (in corrector convergence or error test), RMAX is set to 2
137C for the next increase.
138C-----------------------------------------------------------------------
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
153C-----------------------------------------------------------------------
154C The following block handles preliminaries needed when JSTART = -1.
155C IPUP is set to MITER to force a matrix update.
156C If an order increase is about to be considered (IALTH = 1),
157C IALTH is reset to 2 to postpone consideration one more step.
158C If the caller has changed METH, SCFODE is called to reset
159C the coefficients of the method.
160C If the caller has changed MAXORD to a value less than the current
161C order NQ, NQ is reduced to MAXORD, and a new H chosen accordingly.
162C If H is to be changed, YH must be rescaled.
163C If H or METH is being changed, IALTH is reset to L = NQ + 1
164C to prevent further changes in H for that many steps.
165C-----------------------------------------------------------------------
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
194C-----------------------------------------------------------------------
195C SCFODE is called to get all the integration coefficients for the
196C current METH.  Then the EL vector and related constants are reset
197C whenever the order NQ is changed, or at the start of the problem.
198C-----------------------------------------------------------------------
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
207C-----------------------------------------------------------------------
208C If H is being changed, the H ratio RH is checked against
209C RMAX, HMIN, and HMXI, and the YH array rescaled.  IALTH is set to
210C L = NQ + 1 to prevent a change of H for that many steps, unless
211C forced by a convergence or error test failure.
212C-----------------------------------------------------------------------
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
230C-----------------------------------------------------------------------
231C This section computes the predicted values by effectively
232C multiplying the YH array by the Pascal Triangle matrix.
233C RC is the ratio of new to old values of the coefficient  H*EL(1).
234C When RC differs from 1 by more than CCMAX, IPUP is set to MITER
235C to force PJAC to be called, if a Jacobian is involved.
236C In any case, PJAC is called at least every MSBP steps.
237C-----------------------------------------------------------------------
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
244Cdir$ ivdep
245        DO 210 I = I1,NQNYH
246 210      YH1(I) = YH1(I) + YH1(I+NYH)
247 215    CONTINUE
248C-----------------------------------------------------------------------
249C Up to MAXCOR corrector iterations are taken.  A convergence test is
250C made on the R.M.S. norm of each correction, weighted by the error
251C weight vector EWT.  The sum of the corrections is accumulated in the
252C vector ACOR(i).  The YH array is not altered in the corrector loop.
253C-----------------------------------------------------------------------
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
260C-----------------------------------------------------------------------
261C If indicated, the matrix P = I - h*el(1)*J is reevaluated and
262C preprocessed before starting the corrector iteration.  IPUP is set
263C to 0 as an indicator that this has been done.
264C-----------------------------------------------------------------------
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
274C-----------------------------------------------------------------------
275C In the case of functional iteration, update Y directly from
276C the result of the last function evaluation.
277C-----------------------------------------------------------------------
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
286C-----------------------------------------------------------------------
287C In the case of the chord method, compute the corrector error,
288C and solve the linear system with that as right-hand side and
289C P as coefficient matrix.
290C-----------------------------------------------------------------------
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)
300C-----------------------------------------------------------------------
301C Test for convergence.  If M.gt.0, an estimate of the convergence
302C rate constant is stored in CRATE, and this is used in the test.
303C-----------------------------------------------------------------------
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
314C-----------------------------------------------------------------------
315C The corrector iteration failed to converge.
316C If MITER .ne. 0 and the Jacobian is out of date, PJAC is called for
317C the next try.  Otherwise the YH array is retracted to its values
318C before prediction, and H is reduced, if possible.  If H cannot be
319C reduced or MXNCF failures have occurred, exit with KFLAG = -2.
320C-----------------------------------------------------------------------
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
332Cdir$ 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
343C-----------------------------------------------------------------------
344C The corrector has converged.  JCUR is set to 0
345C to signal that the Jacobian involved may need updating later.
346C The local error test is made and control passes to statement 500
347C if it fails.
348C-----------------------------------------------------------------------
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
353C-----------------------------------------------------------------------
354C After a successful step, update the YH array.
355C Consider changing H if IALTH = 1.  Otherwise decrease IALTH by 1.
356C If IALTH is then 1 and NQ .lt. MAXORD, then ACOR is saved for
357C use in a possible order increase on the next step.
358C If a change in H is considered, an increase or decrease in order
359C by one is considered also.  A change in H is made only if it is by a
360C factor of at least 1.1.  If not, IALTH is set to 3 to prevent
361C testing for that many steps.
362C-----------------------------------------------------------------------
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
378C-----------------------------------------------------------------------
379C The error test failed.  KFLAG keeps track of multiple failures.
380C Restore TN and the YH array to their previous values, and prepare
381C to try the step again.  Compute the optimum step size for this or
382C one lower order.  After 2 or more failures, H is forced to decrease
383C by a factor of 0.2 or less.
384C-----------------------------------------------------------------------
385 500  KFLAG = KFLAG - 1
386      TN = TOLD
387      I1 = NQNYH + 1
388      DO 515 JB = 1,NQ
389        I1 = I1 - NYH
390Cdir$ 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
400C-----------------------------------------------------------------------
401C Regardless of the success or failure of the step, factors
402C RHDN, RHSM, and RHUP are computed, by which H could be multiplied
403C at order NQ - 1, order NQ, or order NQ + 1, respectively.
404C In the case of failure, RHUP = 0.0 to avoid an order increase.
405C The largest of these is determined and the new order chosen
406C accordingly.  If the order is to be increased, we compute one
407C additional scaled derivative.
408C-----------------------------------------------------------------------
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)
445C-----------------------------------------------------------------------
446C If there is a change of order, reset NQ, l, and the coefficients.
447C In any case H is reset according to RH and the YH array is rescaled.
448C Then exit from 690 if the step was OK, or redo the step otherwise.
449C-----------------------------------------------------------------------
450      IF (NEWQ .EQ. NQ) GO TO 170
451 630  NQ = NEWQ
452      L = NQ + 1
453      IRET = 2
454      GO TO 150
455C-----------------------------------------------------------------------
456C Control reaches this section if 3 or more failures have occurred.
457C If 10 failures have occurred, exit with KFLAG = -1.
458C It is assumed that the derivatives that have accumulated in the
459C YH array have errors of the wrong order.  Hence the first
460C derivative is recomputed, and the order is set to 1.  Then
461C H is reduced by a factor of 10, and the step is retried,
462C until it succeeds or H reaches HMIN.
463C-----------------------------------------------------------------------
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
481C-----------------------------------------------------------------------
482C All returns are made through this section.  H is saved in HOLD
483C to allow the caller to change H on the next step.
484C-----------------------------------------------------------------------
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
498C----------------------- END OF SUBROUTINE SSTODE ----------------------
499      END
500