1
subroutine stoda (neq, y, yh, nyh, yh1, ewt, savf, acor,
2
1 wm, iwm, f, jac, pjac, slvs)
4
external f, jac, pjac, slvs
6
integer iownd, ialth, ipup, lmax, meo, nqnyh, nslp,
7
1 icf, ierpj, iersl, jcur, jstart, kflag, l, meth, miter,
8
2 maxord, maxcor, msbp, mxncf, n, nq, nst, nfe, nje, nqu
9
integer iownd2, icount, irflag, jtyp, mused, mxordn, mxords
10
integer i, i1, iredo, iret, j, jb, m, ncf, newq
11
integer lm1, lm1p1, lm2, lm2p1, nqm1, nqm2, isav
12
double precision y, yh, yh1, ewt, savf, acor, wm, rsav
13
double precision conit, crate, el, elco, hold, rmax, tesco,
14
2 ccmax, el0, h, hmin, hmxi, hu, rc, tn, uround
15
double precision rownd2, pdest, pdlast, ratio, cm1, cm2,
17
double precision dcon, ddn, del, delp, dsm, dup, exdn, exsm, exup,
18
1 r, rh, rhdn, rhsm, rhup, told, vmnorm
19
double precision alpha, dm1, dm2, exm1, exm2, pdh, pnorm, rate,
20
1 rh1, rh1it, rh2, rm, sm1
21
dimension neq(1), y(1), yh(nyh,*), yh1(1), ewt(1), savf(1),
22
1 acor(1), wm(*), iwm(*), rsav(240), isav(50)
24
common /ls0001/ conit, crate, el(13), elco(13,12),
25
1 hold, rmax, tesco(3,12),
26
2 ccmax, el0, h, hmin, hmxi, hu, rc, tn, uround, iownd(14),
27
3 ialth, ipup, lmax, meo, nqnyh, nslp,
28
4 icf, ierpj, iersl, jcur, jstart, kflag, l, meth, miter,
29
5 maxord, maxcor, msbp, mxncf, n, nq, nst, nfe, nje, nqu
30
common /lsa001/ rownd2, pdest, pdlast, ratio, cm1(12), cm2(5),
32
2 iownd2(3), icount, irflag, jtyp, mused, mxordn, mxords
33
data sm1/0.5d0, 0.575d0, 0.55d0, 0.45d0, 0.35d0, 0.25d0,
34
1 0.20d0, 0.15d0, 0.10d0, 0.075d0, 0.050d0, 0.025d0/
35
c-----------------------------------------------------------------------
36
c stoda performs one step of the integration of an initial value
37
c problem for a system of ordinary differential equations.
38
c note.. stoda is independent of the value of the iteration method
39
c indicator miter, when this is .ne. 0, and hence is independent
40
c of the type of chord method used, or the jacobian structure.
41
c communication with stoda is done with the following variables..
43
c y = an array of length .ge. n used as the y argument in
44
c all calls to f and jac.
45
c neq = integer array containing problem size in neq(1), and
46
c passed as the neq argument in all calls to f and jac.
47
c yh = an nyh by lmax array containing the dependent variables
48
c and their approximate scaled derivatives, where
49
c lmax = maxord + 1. yh(i,j+1) contains the approximate
50
c j-th derivative of y(i), scaled by h**j/factorial(j)
51
c (j = 0,1,...,nq). on entry for the first step, the first
52
c two columns of yh must be set from the initial values.
53
c nyh = a constant integer .ge. n, the first dimension of yh.
54
c yh1 = a one-dimensional array occupying the same space as yh.
55
c ewt = an array of length n containing multiplicative weights
56
c for local error measurements. local errors in y(i) are
57
c compared to 1.0/ewt(i) in various error tests.
58
c savf = an array of working storage, of length n.
59
c acor = a work array of length n, used for the accumulated
60
c corrections. on a successful return, acor(i) contains
61
c the estimated one-step local error in y(i).
62
c wm,iwm = real and integer work arrays associated with matrix
63
c operations in chord iteration (miter .ne. 0).
64
c pjac = name of routine to evaluate and preprocess jacobian matrix
65
c and p = i - h*el0*jac, if a chord method is being used.
66
c it also returns an estimate of norm(jac) in pdnorm.
67
c slvs = name of routine to solve linear system in chord iteration.
68
c ccmax = maximum relative change in h*el0 before pjac is called.
69
c h = the step size to be attempted on the next step.
70
c h is altered by the error control algorithm during the
71
c problem. h can be either positive or negative, but its
72
c sign must remain constant throughout the problem.
73
c hmin = the minimum absolute value of the step size h to be used.
74
c hmxi = inverse of the maximum absolute value of h to be used.
75
c hmxi = 0.0 is allowed and corresponds to an infinite hmax.
76
c hmin and hmxi may be changed at any time, but will not
77
c take effect until the next change of h is considered.
78
c tn = the independent variable. tn is updated on each step taken.
79
c jstart = an integer used for input only, with the following
80
c values and meanings..
81
c 0 perform the first step.
82
c .gt.0 take a new step continuing from the last.
83
c -1 take the next step with a new value of h,
84
c n, meth, miter, and/or matrix parameters.
85
c -2 take the next step with a new value of h,
86
c but with other inputs unchanged.
87
c on return, jstart is set to 1 to facilitate continuation.
88
c kflag = a completion code with the following meanings..
89
c 0 the step was succesful.
90
c -1 the requested error could not be achieved.
91
c -2 corrector convergence could not be achieved.
92
c -3 fatal error in pjac or slvs.
93
c a return with kflag = -1 or -2 means either
94
c abs(h) = hmin or 10 consecutive failures occurred.
95
c on a return with kflag negative, the values of tn and
96
c the yh array are as of the beginning of the last
97
c step, and h is the last step size attempted.
98
c maxord = the maximum order of integration method to be allowed.
99
c maxcor = the maximum number of corrector iterations allowed.
100
c msbp = maximum number of steps between pjac calls (miter .gt. 0).
101
c mxncf = maximum number of convergence failures allowed.
102
c meth = current method.
103
c meth = 1 means adams method (nonstiff)
104
c meth = 2 means bdf method (stiff)
105
c meth may be reset by stoda.
106
c miter = corrector iteration method.
107
c miter = 0 means functional iteration.
108
c miter = jt .gt. 0 means a chord iteration corresponding
109
c to jacobian type jt. (the lsoda argument jt is
110
c communicated here as jtyp, but is not used in stoda
111
c except to load miter following a method switch.)
112
c miter may be reset by stoda.
113
c n = the number of first-order differential equations.
114
c-----------------------------------------------------------------------
123
if (jstart .gt. 0) go to 200
124
if (jstart .eq. -1) go to 100
125
if (jstart .eq. -2) go to 160
126
c-----------------------------------------------------------------------
127
c on the first call, the order is set to 1, and other variables are
128
c initialized. rmax is the maximum ratio by which h can be increased
129
c in a single step. it is initially 1.e4 to compensate for the small
130
c initial h, but then is normally equal to 10. if a failure
131
c occurs (in corrector convergence or error test), rmax is set at 2
132
c for the next increase.
133
c cfode is called to get the needed coefficients for both methods.
134
c-----------------------------------------------------------------------
147
c initialize switching parameters. meth = 1 is assumed initially. -----
153
call cfode (2, elco, tesco)
155
10 cm2(i) = tesco(2,i)*elco(i+1,i)
156
call cfode (1, elco, tesco)
158
20 cm1(i) = tesco(2,i)*elco(i+1,i)
160
c-----------------------------------------------------------------------
161
c the following block handles preliminaries needed when jstart = -1.
162
c ipup is set to miter to force a matrix update.
163
c if an order increase is about to be considered (ialth = 1),
164
c ialth is reset to 2 to postpone consideration one more step.
165
c if the caller has changed meth, cfode is called to reset
166
c the coefficients of the method.
167
c if h is to be changed, yh must be rescaled.
168
c if h or meth is being changed, ialth is reset to l = nq + 1
169
c to prevent further changes in h for that many steps.
170
c-----------------------------------------------------------------------
173
if (ialth .eq. 1) ialth = 2
174
if (meth .eq. mused) go to 160
175
call cfode (meth, elco, tesco)
178
c-----------------------------------------------------------------------
179
c the el vector and related constants are reset
180
c whenever the order nq is changed, or at the start of the problem.
181
c-----------------------------------------------------------------------
183
155 el(i) = elco(i,nq)
187
conit = 0.5d0/dfloat(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
200
170 rh = dmax1(rh,hmin/dabs(h))
201
175 rh = dmin1(rh,rmax)
202
rh = rh/dmax1(1.0d0,dabs(h)*hmxi*rh)
203
c-----------------------------------------------------------------------
204
c if meth = 1, also restrict the new step size by the stability region.
205
c if this reduces h, set irflag to 1 so that if there are roundoff
206
c problems later, we can assume that is the cause of the trouble.
207
c-----------------------------------------------------------------------
208
if (meth .eq. 2) go to 178
210
pdh = dmax1(dabs(h)*pdlast,0.000001d0)
211
if (rh*pdh*1.00001d0 .lt. sm1(nq)) go to 178
219
180 yh(i,j) = yh(i,j)*r
223
if (iredo .eq. 0) go to 690
224
c-----------------------------------------------------------------------
225
c this section computes the predicted values by effectively
226
c multiplying the yh array by the pascal triangle matrix.
227
c rc is the ratio of new to old values of the coefficient h*el(1).
228
c when rc differs from 1 by more than ccmax, ipup is set to miter
229
c to force pjac to be called, if a jacobian is involved.
230
c in any case, pjac is called at least every msbp steps.
231
c-----------------------------------------------------------------------
232
200 if (dabs(rc-1.0d0) .gt. ccmax) ipup = miter
233
if (nst .ge. nslp+msbp) ipup = miter
240
210 yh1(i) = yh1(i) + yh1(i+nyh)
242
pnorm = vmnorm (n, yh1, ewt)
243
c-----------------------------------------------------------------------
244
c up to maxcor corrector iterations are taken. a convergence test is
245
c made on the r.m.s. norm of each correction, weighted by the error
246
c weight vector ewt. the sum of the corrections is accumulated in the
247
c vector acor(i). the yh array is not altered in the corrector loop.
248
c-----------------------------------------------------------------------
254
call srcma (rsav, isav, 1)
255
call f (neq, tn, y, savf)
256
call srcma (rsav, isav, 2)
258
if (ipup .le. 0) go to 250
259
c-----------------------------------------------------------------------
260
c if indicated, the matrix p = i - h*el(1)*j is reevaluated and
261
c preprocessed before starting the corrector iteration. ipup is set
262
c to 0 as an indicator that this has been done.
263
c-----------------------------------------------------------------------
264
call pjac (neq, y, yh, nyh, ewt, acor, savf, wm, iwm, f, jac)
269
if (ierpj .ne. 0) go to 430
272
270 if (miter .ne. 0) go to 350
273
c-----------------------------------------------------------------------
274
c in the case of functional iteration, update y directly from
275
c the result of the last function evaluation.
276
c-----------------------------------------------------------------------
278
savf(i) = h*savf(i) - yh(i,2)
279
290 y(i) = savf(i) - acor(i)
280
del = vmnorm (n, y, ewt)
282
y(i) = yh(i,1) + el(1)*savf(i)
283
300 acor(i) = savf(i)
285
c-----------------------------------------------------------------------
286
c in the case of the chord method, compute the corrector error,
287
c and solve the linear system with that as right-hand side and
288
c p as coefficient matrix.
289
c-----------------------------------------------------------------------
291
360 y(i) = h*savf(i) - (yh(i,2) + acor(i))
292
call slvs (wm, iwm, y, savf)
293
if (iersl .lt. 0) go to 430
294
if (iersl .gt. 0) go to 410
295
del = vmnorm (n, y, ewt)
297
acor(i) = acor(i) + y(i)
298
380 y(i) = yh(i,1) + el(1)*acor(i)
299
c-----------------------------------------------------------------------
300
c test for convergence. if m.gt.0, an estimate of the convergence
301
c rate constant is stored in crate, and this is used in the test.
303
c we first check for a change of iterates that is the size of
304
c roundoff error. if this occurs, the iteration has converged, and a
305
c new rate estimate is not formed.
306
c in all other cases, force at least two iterations to estimate a
307
c local lipschitz constant estimate for adams methods.
308
c on convergence, form pdest = local maximum lipschitz constant
309
c estimate. pdlast is the most recent nonzero estimate.
310
c-----------------------------------------------------------------------
312
if (del .le. 100.0d0*pnorm*uround) go to 450
313
if (m .eq. 0 .and. meth .eq. 1) go to 405
314
if (m .eq. 0) go to 402
316
if (del .le. 1024.0d0*delp) rm = del/delp
317
rate = dmax1(rate,rm)
318
crate = dmax1(0.2d0*crate,rm)
319
402 dcon = del*dmin1(1.0d0,1.5d0*crate)/(tesco(2,nq)*conit)
320
if (dcon .gt. 1.0d0) go to 405
321
pdest = dmax1(pdest,rate/dabs(h*el(1)))
322
if (pdest .ne. 0.0d0) pdlast = pdest
326
if (m .eq. maxcor) go to 410
327
if (m .ge. 2 .and. del .gt. 2.0d0*delp) go to 410
329
call srcma (rsav, isav, 1)
330
call f (neq, tn, y, savf)
331
call srcma (rsav, isav, 2)
334
c-----------------------------------------------------------------------
335
c the corrector iteration failed to converge.
336
c if miter .ne. 0 and the jacobian is out of date, pjac is called for
337
c the next try. otherwise the yh array is retracted to its values
338
c before prediction, and h is reduced, if possible. if h cannot be
339
c reduced or mxncf failures have occurred, exit with kflag = -2.
340
c-----------------------------------------------------------------------
341
410 if (miter .eq. 0 .or. jcur .eq. 1) go to 430
354
440 yh1(i) = yh1(i) - yh1(i+nyh)
356
if (ierpj .lt. 0 .or. iersl .lt. 0) go to 680
357
if (dabs(h) .le. hmin*1.00001d0) go to 670
358
if (ncf .eq. mxncf) go to 670
363
c-----------------------------------------------------------------------
364
c the corrector has converged. jcur is set to 0
365
c to signal that the jacobian involved may need updating later.
366
c the local error test is made and control passes to statement 500
368
c-----------------------------------------------------------------------
370
if (m .eq. 0) dsm = del/tesco(2,nq)
371
if (m .gt. 0) dsm = vmnorm (n, acor, ewt)/tesco(2,nq)
372
if (dsm .gt. 1.0d0) go to 500
373
c-----------------------------------------------------------------------
374
c after a successful step, update the yh array.
375
c decrease icount by 1, and if it is -1, consider switching methods.
376
c if a method switch is made, reset various parameters,
377
c rescale the yh array, and exit. if there is no switch,
378
c consider changing h if ialth = 1. otherwise decrease ialth by 1.
379
c if ialth is then 1 and nq .lt. maxord, then acor is saved for
380
c use in a possible order increase on the next step.
381
c if a change in h is considered, an increase or decrease in order
382
c by one is considered also. a change in h is made only if it is by a
383
c factor of at least 1.1. if not, ialth is set to 3 to prevent
384
c testing for that many steps.
385
c-----------------------------------------------------------------------
394
460 yh(i,j) = yh(i,j) + el(j)*acor(i)
396
if (icount .ge. 0) go to 488
397
if (meth .eq. 2) go to 480
398
c-----------------------------------------------------------------------
399
c we are currently using an adams method. consider switching to bdf.
400
c if the current order is greater than 5, assume the problem is
401
c not stiff, and skip this section.
402
c if the lipschitz constant and error estimate are not polluted
403
c by roundoff, go to 470 and perform the usual test.
404
c otherwise, switch to the bdf methods if the last step was
405
c restricted to insure stability (irflag = 1), and stay with adams
406
c method if not. when switching to bdf with polluted error estimates,
407
c in the absence of other information, double the step size.
409
c when the estimates are ok, we make the usual test by computing
410
c the step size we could have (ideally) used on this step,
411
c with the current (adams) method, and also that for the bdf.
412
c if nq .gt. mxords, we consider changing to order mxords on switching.
413
c compare the two step sizes to decide whether to switch.
414
c the step size advantage must be at least ratio = 5 to switch.
415
c-----------------------------------------------------------------------
416
if (nq .gt. 5) go to 488
417
if (dsm .gt. 100.0d0*pnorm*uround .and. pdest .ne. 0.0d0)
419
if (irflag .eq. 0) go to 488
421
nqm2 = min0(nq,mxords)
424
exsm = 1.0d0/dfloat(l)
425
rh1 = 1.0d0/(1.2d0*dsm**exsm + 0.0000012d0)
428
if (pdh*rh1 .gt. 0.00001d0) rh1it = sm1(nq)/pdh
429
rh1 = dmin1(rh1,rh1it)
430
if (nq .le. mxords) go to 474
433
exm2 = 1.0d0/dfloat(lm2)
435
dm2 = vmnorm (n, yh(1,lm2p1), ewt)/cm2(mxords)
436
rh2 = 1.0d0/(1.2d0*dm2**exm2 + 0.0000012d0)
438
474 dm2 = dsm*(cm1(nq)/cm2(nq))
439
rh2 = 1.0d0/(1.2d0*dm2**exsm + 0.0000012d0)
442
if (rh2 .lt. ratio*rh1) go to 488
443
c the switch test passed. reset relevant quantities for bdf. ----------
452
c-----------------------------------------------------------------------
453
c we are currently using a bdf method. consider switching to adams.
454
c compute the step size we could have (ideally) used on this step,
455
c with the current (bdf) method, and also that for the adams.
456
c if nq .gt. mxordn, we consider changing to order mxordn on switching.
457
c compare the two step sizes to decide whether to switch.
458
c the step size advantage must be at least 5/ratio = 1 to switch.
459
c if the step size for adams would be so small as to cause
460
c roundoff pollution, we stay with bdf.
461
c-----------------------------------------------------------------------
463
exsm = 1.0d0/dfloat(l)
464
if (mxordn .ge. nq) go to 484
467
exm1 = 1.0d0/dfloat(lm1)
469
dm1 = vmnorm (n, yh(1,lm1p1), ewt)/cm1(mxordn)
470
rh1 = 1.0d0/(1.2d0*dm1**exm1 + 0.0000012d0)
472
484 dm1 = dsm*(cm2(nq)/cm1(nq))
473
rh1 = 1.0d0/(1.2d0*dm1**exsm + 0.0000012d0)
476
486 rh1it = 2.0d0*rh1
478
if (pdh*rh1 .gt. 0.00001d0) rh1it = sm1(nqm1)/pdh
479
rh1 = dmin1(rh1,rh1it)
480
rh2 = 1.0d0/(1.2d0*dsm**exsm + 0.0000012d0)
481
if (rh1*ratio .lt. 5.0d0*rh2) go to 488
482
alpha = dmax1(0.001d0,rh1)
483
dm1 = (alpha**exm1)*dm1
484
if (dm1 .le. 1000.0d0*uround*pnorm) go to 488
485
c the switch test passed. reset relevant quantities for adams. --------
495
c no method switch is being made. do the usual step/order selection. --
498
if (ialth .eq. 0) go to 520
499
if (ialth .gt. 1) go to 700
500
if (l .eq. lmax) go to 700
502
490 yh(i,lmax) = acor(i)
504
c-----------------------------------------------------------------------
505
c the error test failed. kflag keeps track of multiple failures.
506
c restore tn and the yh array to their previous values, and prepare
507
c to try the step again. compute the optimum step size for this or
508
c one lower order. after 2 or more failures, h is forced to decrease
509
c by a factor of 0.2 or less.
510
c-----------------------------------------------------------------------
511
500 kflag = kflag - 1
518
510 yh1(i) = yh1(i) - yh1(i+nyh)
521
if (dabs(h) .le. hmin*1.00001d0) go to 660
522
if (kflag .le. -3) go to 640
526
c-----------------------------------------------------------------------
527
c regardless of the success or failure of the step, factors
528
c rhdn, rhsm, and rhup are computed, by which h could be multiplied
529
c at order nq - 1, order nq, or order nq + 1, respectively.
530
c in the case of failure, rhup = 0.0 to avoid an order increase.
531
c the largest of these is determined and the new order chosen
532
c accordingly. if the order is to be increased, we compute one
533
c additional scaled derivative.
534
c-----------------------------------------------------------------------
536
if (l .eq. lmax) go to 540
538
530 savf(i) = acor(i) - yh(i,lmax)
539
dup = vmnorm (n, savf, ewt)/tesco(3,nq)
540
exup = 1.0d0/dfloat(l+1)
541
rhup = 1.0d0/(1.4d0*dup**exup + 0.0000014d0)
542
540 exsm = 1.0d0/dfloat(l)
543
rhsm = 1.0d0/(1.2d0*dsm**exsm + 0.0000012d0)
545
if (nq .eq. 1) go to 550
546
ddn = vmnorm (n, yh(1,l), ewt)/tesco(1,nq)
547
exdn = 1.0d0/dfloat(nq)
548
rhdn = 1.0d0/(1.3d0*ddn**exdn + 0.0000013d0)
549
c if meth = 1, limit rh according to the stability region also. --------
550
550 if (meth .eq. 2) go to 560
551
pdh = dmax1(dabs(h)*pdlast,0.000001d0)
552
if (l .lt. lmax) rhup = dmin1(rhup,sm1(l)/pdh)
553
rhsm = dmin1(rhsm,sm1(nq)/pdh)
554
if (nq .gt. 1) rhdn = dmin1(rhdn,sm1(nq-1)/pdh)
556
560 if (rhsm .ge. rhup) go to 570
557
if (rhup .gt. rhdn) go to 590
559
570 if (rhsm .lt. rhdn) go to 580
565
if (kflag .lt. 0 .and. rh .gt. 1.0d0) rh = 1.0d0
569
if (rh .lt. 1.1d0) go to 610
572
600 yh(i,newq+1) = acor(i)*r
576
c if meth = 1 and h is restricted by stability, bypass 10 percent test.
577
620 if (meth .eq. 2) go to 622
578
if (rh*pdh*1.00001d0 .ge. sm1(newq)) go to 625
579
622 if (kflag .eq. 0 .and. rh .lt. 1.1d0) go to 610
580
625 if (kflag .le. -2) rh = dmin1(rh,0.2d0)
581
c-----------------------------------------------------------------------
582
c if there is a change of order, reset nq, l, and the coefficients.
583
c in any case h is reset according to rh and the yh array is rescaled.
584
c then exit from 690 if the step was ok, or redo the step otherwise.
585
c-----------------------------------------------------------------------
586
if (newq .eq. nq) go to 170
591
c-----------------------------------------------------------------------
592
c control reaches this section if 3 or more failures have occured.
593
c if 10 failures have occurred, exit with kflag = -1.
594
c it is assumed that the derivatives that have accumulated in the
595
c yh array have errors of the wrong order. hence the first
596
c derivative is recomputed, and the order is set to 1. then
597
c h is reduced by a factor of 10, and the step is retried,
598
c until it succeeds or h reaches hmin.
599
c-----------------------------------------------------------------------
600
640 if (kflag .eq. -10) go to 660
602
rh = dmax1(hmin/dabs(h),rh)
606
call srcma (rsav, isav, 1)
607
call f (neq, tn, y, savf)
608
call srcma (rsav, isav, 2)
611
650 yh(i,2) = h*savf(i)
614
if (nq .eq. 1) go to 200
619
c-----------------------------------------------------------------------
620
c all returns are made through this section. h is saved in hold
621
c to allow the caller to change h on the next step.
622
c-----------------------------------------------------------------------
630
700 r = 1.0d0/tesco(2,nqu)
632
710 acor(i) = acor(i)*r
636
c----------------------- end of subroutine stoda -----------------------
added
removed
Lib/integrate/odepack/stoda.f
