1 SUBROUTINE sstode (NEQ, Y, YH, NYH, YH1, EWT, SAVF, ACOR,
2 1 wm, iwm,
f, jac, pjac, slvs)
95 EXTERNAL f, jac, pjac, slvs
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
126 IF (jstart .GT. 0) go to 200
127 IF (jstart .EQ. -1) go to 100
128 IF (jstart .EQ. -2) go to 160
166 IF (ialth .EQ. 1) ialth = 2
167 IF (meth .EQ. meo) go to 110
168 CALL
scfode(meth, elco, tesco)
170 IF (nq .GT. maxord) go to 120
174 110
IF (nq .LE. maxord) go to 160
178 125 el(i) = elco(i,nq)
183 ddn = svnorm(n, savf, ewt)/tesco(1,l)
185 rhdn = 1.0e0/(1.3e0*ddn**exdn + 0.0000013e0)
188 IF (h .EQ. hold) go to 170
197 140 CALL
scfode(meth, elco, tesco)
199 155 el(i) = elco(i,nq)
204 go to(160, 170, 200), iret
211 160
IF (h .EQ. hold) go to 200
216 170 rh =
max(rh,hmin/
abs(h))
217 175 rh =
min(rh,rmax)
218 rh = rh/
max(1.0e0,
abs(h)*hmxi*rh)
223 180 yh(i,j) = yh(i,j)*r
227 IF (iredo .EQ. 0) go to 690
236 200
IF (
abs(rc-1.0e0) .GT. ccmax) ipup = miter
237 IF (nst .GE. nslp+msbp) ipup = miter
244 210 yh1(i) = yh1(i) + yh1(i+nyh)
255 CALL
f(neq, tn, y, savf)
257 IF (ipup .LE. 0) go to 250
263 CALL pjac(neq, y, yh, nyh, ewt, acor, savf, wm, iwm,
f, jac)
268 IF (ierpj .NE. 0) go to 430
271 270
IF (miter .NE. 0) go to 350
277 savf(i) = h*savf(i) - yh(i,2)
278 290 y(i) = savf(i) - acor(i)
279 del = svnorm(n, y, ewt)
281 y(i) = yh(i,1) + el(1)*savf(i)
282 300 acor(i) = savf(i)
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)
296 acor(i) = acor(i) + y(i)
297 380 y(i) = yh(i,1) + el(1)*acor(i)
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
306 IF (m .EQ. maxcor) go to 410
307 IF (m .GE. 2 .AND. del .GT. 2.0e0*delp) go to 410
309 CALL
f(neq, tn, y, savf)
319 410
IF (miter .EQ. 0 .OR. jcur .EQ. 1) go to 430
332 440 yh1(i) = yh1(i) - yh1(i+nyh)
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
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
368 470 yh(i,j) = yh(i,j) + el(j)*acor(i)
370 IF (ialth .EQ. 0) go to 520
371 IF (ialth .GT. 1) go to 700
372 IF (l .EQ. lmax) go to 700
374 490 yh(i,lmax) = acor(i)
383 500 kflag = kflag - 1
390 510 yh1(i) = yh1(i) - yh1(i+nyh)
393 IF (
abs(h) .LE. hmin*1.00001e0) go to 660
394 IF (kflag .LE. -3) go to 640
408 IF (l .EQ. lmax) go to 540
410 530 savf(i) = acor(i) - yh(i,lmax)
411 dup = svnorm(n, savf, ewt)/tesco(3,nq)
413 rhup = 1.0e0/(1.4e0*dup**exup + 0.0000014e0)
415 rhsm = 1.0e0/(1.2e0*dsm**exsm + 0.0000012e0)
417 IF (nq .EQ. 1) go to 560
418 ddn = svnorm(n, yh(1,l), ewt)/tesco(1,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
424 570
IF (rhsm .LT. rhdn) go to 580
430 IF (kflag .LT. 0 .AND. rh .GT. 1.0e0) rh = 1.0e0
434 IF (rh .LT. 1.1e0) go to 610
437 600 yh(i,newq+1) = acor(i)*r
441 620
IF ((kflag .EQ. 0) .AND. (rh .LT. 1.1e0)) go to 610
442 IF (kflag .LE. -2) rh =
min(rh,0.2e0)
448 IF (newq .EQ. nq) go to 170
462 640
IF (kflag .EQ. -10) go to 660
468 CALL
f(neq, tn, y, savf)
471 650 yh(i,2) = h*savf(i)
474 IF (nq .EQ. 1) go to 200
490 700 r = 1.0e0/tesco(2,nqu)
492 710 acor(i) = acor(i)*r
std::string dimension(void) const
subroutine scfode(METH, ELCO, TESCO)
F77_RET_T const double const double * f
charNDArray max(char d, const charNDArray &m)
subroutine sstode(NEQ, Y, YH, NYH, YH1, EWT, SAVF, ACOR, WM, IWM, F, JAC, PJAC, SLVS)
charNDArray min(char d, const charNDArray &m)