2 SUBROUTINE slsode (F, NEQ, Y, T, TOUT, ITOL, RTOL, ATOL, ITASK,
3 1 istate, iopt, rwork, lrw, iwork, liw, jac, mf)
5 INTEGER NEQ, ITOL, ITASK, ISTATE, IOPT, LRW, IWORK, LIW, MF
6 REAL Y, T, TOUT, RTOL, ATOL, RWORK
7 dimension neq(*), y(*), rtol(*), atol(*), rwork(lrw), iwork(liw)
1210 INTEGER INIT, MXSTEP, MXHNIL, NHNIL, NSLAST, NYH,
1211 1 ialth, ipup, lmax, meo, nqnyh, nslp,
1212 1 icf, ierpj, iersl, jcur, jstart, kflag, l,
1213 2 lyh, lewt, lacor, lsavf, lwm, liwm, meth, miter,
1214 3 maxord, maxcor, msbp, mxncf, n, nq, nst, nfe, nje, nqu
1215 INTEGER I, I1, I2, IFLAG, IMXER, KGO, LF0,
1216 1 leniw, lenrw, lenwm, ml, mord, mu, mxhnl0, mxstp0
1217 REAL CONIT, CRATE, EL, ELCO, HOLD, RMAX, TESCO,
1218 1 ccmax, el0, h, hmin, hmxi, hu, rc, tn, uround
1219 REAL ATOLI, AYI, BIG, EWTI, H0, HMAX, HMX, RH, RTOLI,
1220 1 tcrit, tdist, tnext, tol, tolsf, tp,
SIZE, sum, w0
1224 SAVE mord, mxstp0, mxhnl0
1235 COMMON /sls001/ conit, crate, el(13), elco(13,12),
1236 1 hold, rmax, tesco(3,12),
1237 1 ccmax, el0, h, hmin, hmxi, hu, rc, tn, uround,
1238 2 init, mxstep, mxhnil, nhnil, nslast, nyh,
1239 3 ialth, ipup, lmax, meo, nqnyh, nslp,
1240 3 icf, ierpj, iersl, jcur, jstart, kflag, l,
1241 4 lyh, lewt, lacor, lsavf, lwm, liwm, meth, miter,
1242 5 maxord, maxcor, msbp, mxncf, n, nq, nst, nfe, nje, nqu
1244 DATA mord(1),mord(2)/12,5/, mxstp0/500/, mxhnl0/10/
1255 IF (istate .LT. 1 .OR. istate .GT. 3) go
to 601
1256 IF (itask .LT. 1 .OR. itask .GT. 5) go
to 602
1257 IF (istate .EQ. 1) go
to 10
1258 IF (init .EQ. 0) go
to 603
1259 IF (istate .EQ. 2) go
to 200
1262 IF (tout .EQ. t)
RETURN
1272 20
IF (neq(1) .LE. 0) go
to 604
1273 IF (istate .EQ. 1) go
to 25
1274 IF (neq(1) .GT. n) go
to 605
1276 IF (itol .LT. 1 .OR. itol .GT. 4) go
to 606
1277 IF (iopt .LT. 0 .OR. iopt .GT. 1) go
to 607
1279 miter = mf - 10*meth
1280 IF (meth .LT. 1 .OR. meth .GT. 2) go
to 608
1281 IF (miter .LT. 0 .OR. miter .GT. 5) go
to 608
1282 IF (miter .LE. 3) go
to 30
1285 IF (ml .LT. 0 .OR. ml .GE. n) go
to 609
1286 IF (mu .LT. 0 .OR. mu .GE. n) go
to 610
1289 IF (iopt .EQ. 1) go
to 40
1293 IF (istate .EQ. 1) h0 = 0.0e0
1297 40 maxord = iwork(5)
1298 IF (maxord .LT. 0) go
to 611
1299 IF (maxord .EQ. 0) maxord = 100
1300 maxord =
min(maxord,mord(meth))
1302 IF (mxstep .LT. 0) go
to 612
1303 IF (mxstep .EQ. 0) mxstep = mxstp0
1305 IF (mxhnil .LT. 0) go
to 613
1306 IF (mxhnil .EQ. 0) mxhnil = mxhnl0
1307 IF (istate .NE. 1) go
to 50
1309 IF ((tout - t)*h0 .LT. 0.0e0) go
to 614
1311 IF (hmax .LT. 0.0e0) go
to 615
1313 IF (hmax .GT. 0.0e0) hmxi = 1.0e0/hmax
1315 IF (hmin .LT. 0.0e0) go
to 616
1323 IF (istate .EQ. 1) nyh = n
1324 lwm = lyh + (maxord + 1)*nyh
1325 IF (miter .EQ. 0) lenwm = 0
1326 IF (miter .EQ. 1 .OR. miter .EQ. 2) lenwm = n*n + 2
1327 IF (miter .EQ. 3) lenwm = n + 2
1328 IF (miter .GE. 4) lenwm = (2*ml + mu + 1)*n + 2
1332 lenrw = lacor + n - 1
1336 IF (miter .EQ. 0 .OR. miter .EQ. 3) leniw = 20
1338 IF (lenrw .GT. lrw) go
to 617
1339 IF (leniw .GT. liw) go
to 618
1344 IF (itol .GE. 3) rtoli = rtol(i)
1345 IF (itol .EQ. 2 .OR. itol .EQ. 4) atoli = atol(i)
1346 IF (rtoli .LT. 0.0e0) go
to 619
1347 IF (atoli .LT. 0.0e0) go
to 620
1349 IF (istate .EQ. 1) go
to 100
1352 IF (nq .LE. maxord) go
to 90
1355 80 rwork(i+lsavf-1) = rwork(i+lwm-1)
1357 90
IF (miter .GT. 0) rwork(lwm) =
sqrt(uround)
1358 IF (n .EQ. nyh) go
to 200
1361 i2 = lyh + (maxord + 1)*nyh - 1
1362 IF (i1 .GT. i2) go
to 200
1373 100 uround = r1mach(4)
1375 IF (itask .NE. 4 .AND. itask .NE. 5) go
to 110
1377 IF ((tcrit - tout)*(tout - t) .LT. 0.0e0) go
to 625
1378 IF (h0 .NE. 0.0e0 .AND. (t + h0 - tcrit)*h0 .GT. 0.0e0)
1381 IF (miter .GT. 0) rwork(lwm) =
sqrt(uround)
1394 CALL
f(neq, t, y, rwork(lf0))
1398 115 rwork(i+lyh-1) = y(i)
1402 CALL
sewset(n, itol, rtol, atol, rwork(lyh), rwork(lewt))
1404 IF (rwork(i+lewt-1) .LE. 0.0e0) go
to 621
1405 120 rwork(i+lewt-1) = 1.0e0/rwork(i+lewt-1)
1422 IF (h0 .NE. 0.0e0) go
to 180
1423 tdist =
abs(tout - t)
1425 IF (tdist .LT. 2.0e0*uround*w0) go
to 622
1427 IF (itol .LE. 2) go
to 140
1429 130 tol =
max(tol,rtol(i))
1430 140
IF (tol .GT. 0.0e0) go
to 160
1433 IF (itol .EQ. 2 .OR. itol .EQ. 4) atoli = atol(i)
1435 IF (ayi .NE. 0.0e0) tol =
max(tol,atoli/ayi)
1437 160 tol =
max(tol,100.0e0*uround)
1438 tol =
min(tol,0.001e0)
1439 sum = svnorm(n, rwork(lf0), rwork(lewt))
1440 sum = 1.0e0/(tol*w0*w0) + tol*sum**2
1441 h0 = 1.0e0/
sqrt(sum)
1443 h0 = sign(h0,tout-t)
1445 180 rh =
abs(h0)*hmxi
1446 IF (rh .GT. 1.0e0) h0 = h0/rh
1450 190 rwork(i+lf0-1) = h0*rwork(i+lf0-1)
1458 go
to(210, 250, 220, 230, 240), itask
1459 210
IF ((tn - tout)*h .LT. 0.0e0) go
to 250
1460 CALL
sintdy(tout, 0, rwork(lyh), nyh, y, iflag)
1461 IF (iflag .NE. 0) go
to 627
1464 220 tp = tn - hu*(1.0e0 + 100.0e0*uround)
1465 IF ((tp - tout)*h .GT. 0.0e0) go
to 623
1466 IF ((tn - tout)*h .LT. 0.0e0) go
to 250
1468 230 tcrit = rwork(1)
1469 IF ((tn - tcrit)*h .GT. 0.0e0) go
to 624
1470 IF ((tcrit - tout)*h .LT. 0.0e0) go
to 625
1471 IF ((tn - tout)*h .LT. 0.0e0) go
to 245
1472 CALL
sintdy(tout, 0, rwork(lyh), nyh, y, iflag)
1473 IF (iflag .NE. 0) go
to 627
1476 240 tcrit = rwork(1)
1477 IF ((tn - tcrit)*h .GT. 0.0e0) go
to 624
1478 245 hmx =
abs(tn) +
abs(h)
1479 ihit =
abs(tn - tcrit) .LE. 100.0e0*uround*hmx
1481 tnext = tn + h*(1.0e0 + 4.0e0*uround)
1482 IF ((tnext - tcrit)*h .LE. 0.0e0) go
to 250
1483 h = (tcrit - tn)*(1.0e0 - 4.0e0*uround)
1484 IF (istate .EQ. 2) jstart = -2
1497 IF ((nst-nslast) .GE. mxstep) go
to 500
1498 CALL
sewset(n, itol, rtol, atol, rwork(lyh), rwork(lewt))
1500 IF (rwork(i+lewt-1) .LE. 0.0e0) go
to 510
1501 260 rwork(i+lewt-1) = 1.0e0/rwork(i+lewt-1)
1502 270 tolsf = uround*svnorm(n, rwork(lyh), rwork(lewt))
1503 IF (tolsf .LE. 1.0e0) go
to 280
1505 IF (nst .EQ. 0) go
to 626
1507 280
IF ((tn + h) .NE. tn) go
to 290
1509 IF (nhnil .GT. mxhnil) go
to 290
1510 CALL
xerrwd(
'SLSODE- Warning..internal T (=R1) and H (=R2) are',
1511 1 50, 101, 0, 0, 0, 0, 0, 0.0e0, 0.0e0)
1513 1
' such that in the machine, T + H = T on the next step ',
1514 1 60, 101, 0, 0, 0, 0, 0, 0.0e0, 0.0e0)
1515 CALL
xerrwd(
' (H = step size). Solver will continue anyway',
1516 1 50, 101, 0, 0, 0, 0, 2, tn, h)
1517 IF (nhnil .LT. mxhnil) go
to 290
1518 CALL
xerrwd(
'SLSODE- Above warning has been issued I1 times. ',
1519 1 50, 102, 0, 0, 0, 0, 0, 0.0e0, 0.0e0)
1520 CALL
xerrwd(
' It will not be issued again for this problem',
1521 1 50, 102, 0, 1, mxhnil, 0, 0, 0.0e0, 0.0e0)
1526 CALL
sstode(neq, y, rwork(lyh), nyh, rwork(lyh), rwork(lewt),
1527 1 rwork(lsavf), rwork(lacor), rwork(lwm), iwork(liwm),
1530 go
to(300, 530, 540), kgo
1537 go
to(310, 400, 330, 340, 350), itask
1539 310
IF ((tn - tout)*h .LT. 0.0e0) go
to 250
1540 CALL
sintdy(tout, 0, rwork(lyh), nyh, y, iflag)
1544 330
IF ((tn - tout)*h .GE. 0.0e0) go
to 400
1547 340
IF ((tn - tout)*h .LT. 0.0e0) go
to 345
1548 CALL
sintdy(tout, 0, rwork(lyh), nyh, y, iflag)
1551 345 hmx =
abs(tn) +
abs(h)
1552 ihit =
abs(tn - tcrit) .LE. 100.0e0*uround*hmx
1554 tnext = tn + h*(1.0e0 + 4.0e0*uround)
1555 IF ((tnext - tcrit)*h .LE. 0.0e0) go
to 250
1556 h = (tcrit - tn)*(1.0e0 - 4.0e0*uround)
1560 350 hmx =
abs(tn) +
abs(h)
1561 ihit =
abs(tn - tcrit) .LE. 100.0e0*uround*hmx
1570 410 y(i) = rwork(i+lyh-1)
1572 IF (itask .NE. 4 .AND. itask .NE. 5) go
to 420
1593 500 CALL
xerrwd(
'SLSODE- At current T (=R1), MXSTEP (=I1) steps ',
1594 1 50, 201, 0, 0, 0, 0, 0, 0.0e0, 0.0e0)
1595 CALL
xerrwd(
' taken on this call before reaching TOUT ',
1596 1 50, 201, 0, 1, mxstep, 0, 1, tn, 0.0e0)
1600 510 ewti = rwork(lewt+i-1)
1601 CALL
xerrwd(.LE.
'SLSODE- At T (=R1), EWT(I1) has become R2 0.',
1602 1 50, 202, 0, 1, i, 0, 2, tn, ewti)
1606 520 CALL
xerrwd(
'SLSODE- At T (=R1), too much accuracy requested ',
1607 1 50, 203, 0, 0, 0, 0, 0, 0.0e0, 0.0e0)
1608 CALL
xerrwd(
' for precision of machine.. see TOLSF (=R2) ',
1609 1 50, 203, 0, 0, 0, 0, 2, tn, tolsf)
1614 530 CALL
xerrwd(
'SLSODE- At T(=R1) and step size H(=R2), the error',
1615 1 50, 204, 0, 0, 0, 0, 0, 0.0e0, 0.0e0)
1616 CALL
xerrwd(
' test failed repeatedly or with ABS(H) = HMIN',
1617 1 50, 204, 0, 0, 0, 0, 2, tn, h)
1621 540 CALL
xerrwd(
'SLSODE- At T (=R1) and step size H (=R2), the ',
1622 1 50, 205, 0, 0, 0, 0, 0, 0.0e0, 0.0e0)
1623 CALL
xerrwd(
' corrector convergence failed repeatedly ',
1624 1 50, 205, 0, 0, 0, 0, 0, 0.0e0, 0.0e0)
1625 CALL
xerrwd(
' or with ABS(H) = HMIN ',
1626 1 30, 205, 0, 0, 0, 0, 2, tn, h)
1632 SIZE =
abs(rwork(i+lacor-1)*rwork(i+lewt-1))
1633 IF (big .GE.
size) go
to 570
1640 590 y(i) = rwork(i+lyh-1)
1658 601 CALL
xerrwd(
'SLSODE- ISTATE (=I1) illegal ',
1659 1 30, 1, 0, 1, istate, 0, 0, 0.0e0, 0.0e0)
1660 IF (istate .LT. 0) go
to 800
1662 602 CALL
xerrwd(
'SLSODE- ITASK (=I1) illegal ',
1663 1 30, 2, 0, 1, itask, 0, 0, 0.0e0, 0.0e0)
1665 603 CALL
xerrwd(.GT.
'SLSODE- ISTATE 1 but SLSODE not initialized ',
1666 1 50, 3, 0, 0, 0, 0, 0, 0.0e0, 0.0e0)
1668 604 CALL
xerrwd(.LT.
'SLSODE- NEQ (=I1) 1 ',
1669 1 30, 4, 0, 1, neq(1), 0, 0, 0.0e0, 0.0e0)
1671 605 CALL
xerrwd(
'SLSODE- ISTATE = 3 and NEQ increased (I1 to I2) ',
1672 1 50, 5, 0, 2, n, neq(1), 0, 0.0e0, 0.0e0)
1674 606 CALL
xerrwd(
'SLSODE- ITOL (=I1) illegal ',
1675 1 30, 6, 0, 1, itol, 0, 0, 0.0e0, 0.0e0)
1677 607 CALL
xerrwd(
'SLSODE- IOPT (=I1) illegal ',
1678 1 30, 7, 0, 1, iopt, 0, 0, 0.0e0, 0.0e0)
1680 608 CALL
xerrwd(
'SLSODE- MF (=I1) illegal ',
1681 1 30, 8, 0, 1, mf, 0, 0, 0.0e0, 0.0e0)
1683 609 CALL
xerrwd(.LT..GE.
'SLSODE- ML (=I1) illegal.. 0 or NEQ (=I2)',
1684 1 50, 9, 0, 2, ml, neq(1), 0, 0.0e0, 0.0e0)
1686 610 CALL
xerrwd(.LT..GE.
'SLSODE- MU (=I1) illegal.. 0 or NEQ (=I2)',
1687 1 50, 10, 0, 2, mu, neq(1), 0, 0.0e0, 0.0e0)
1689 611 CALL
xerrwd(.LT.
'SLSODE- MAXORD (=I1) 0 ',
1690 1 30, 11, 0, 1, maxord, 0, 0, 0.0e0, 0.0e0)
1692 612 CALL
xerrwd(.LT.
'SLSODE- MXSTEP (=I1) 0 ',
1693 1 30, 12, 0, 1, mxstep, 0, 0, 0.0e0, 0.0e0)
1695 613 CALL
xerrwd(.LT.
'SLSODE- MXHNIL (=I1) 0 ',
1696 1 30, 13, 0, 1, mxhnil, 0, 0, 0.0e0, 0.0e0)
1698 614 CALL
xerrwd(
'SLSODE- TOUT (=R1) behind T (=R2) ',
1699 1 40, 14, 0, 0, 0, 0, 2, tout, t)
1700 CALL
xerrwd(
' Integration direction is given by H0 (=R1) ',
1701 1 50, 14, 0, 0, 0, 0, 1, h0, 0.0e0)
1703 615 CALL
xerrwd(.LT.
'SLSODE- HMAX (=R1) 0.0 ',
1704 1 30, 15, 0, 0, 0, 0, 1, hmax, 0.0e0)
1706 616 CALL
xerrwd(.LT.
'SLSODE- HMIN (=R1) 0.0 ',
1707 1 30, 16, 0, 0, 0, 0, 1, hmin, 0.0e0)
1710 1
'SLSODE- RWORK length needed, LENRW (=I1), exceeds LRW (=I2)',
1711 1 60, 17, 0, 2, lenrw, lrw, 0, 0.0e0, 0.0e0)
1714 1
'SLSODE- IWORK length needed, LENIW (=I1), exceeds LIW (=I2)',
1715 1 60, 18, 0, 2, leniw, liw, 0, 0.0e0, 0.0e0)
1717 619 CALL
xerrwd(.LT.
'SLSODE- RTOL(I1) is R1 0.0 ',
1718 1 40, 19, 0, 1, i, 0, 1, rtoli, 0.0e0)
1720 620 CALL
xerrwd(.LT.
'SLSODE- ATOL(I1) is R1 0.0 ',
1721 1 40, 20, 0, 1, i, 0, 1, atoli, 0.0e0)
1723 621 ewti = rwork(lewt+i-1)
1724 CALL
xerrwd(.LE.
'SLSODE- EWT(I1) is R1 0.0 ',
1725 1 40, 21, 0, 1, i, 0, 1, ewti, 0.0e0)
1728 1
'SLSODE- TOUT (=R1) too close to T(=R2) to start integration',
1729 1 60, 22, 0, 0, 0, 0, 2, tout, t)
1732 1
'SLSODE- ITASK = I1 and TOUT (=R1) behind TCUR - HU (= R2) ',
1733 1 60, 23, 0, 1, itask, 0, 2, tout, tp)
1736 1
'SLSODE- ITASK = 4 OR 5 and TCRIT (=R1) behind TCUR (=R2) ',
1737 1 60, 24, 0, 0, 0, 0, 2, tcrit, tn)
1740 1
'SLSODE- ITASK = 4 or 5 and TCRIT (=R1) behind TOUT (=R2) ',
1741 1 60, 25, 0, 0, 0, 0, 2, tcrit, tout)
1743 626 CALL
xerrwd(
'SLSODE- At start of problem, too much accuracy ',
1744 1 50, 26, 0, 0, 0, 0, 0, 0.0e0, 0.0e0)
1746 1
' requested for precision of machine.. See TOLSF (=R1) ',
1747 1 60, 26, 0, 0, 0, 0, 1, tolsf, 0.0e0)
1750 627 CALL
xerrwd(
'SLSODE- Trouble in SINTDY. ITASK = I1, TOUT = R1',
1751 1 50, 27, 0, 1, itask, 0, 1, tout, 0.0e0)
1756 800 CALL
xerrwd(
'SLSODE- Run aborted.. apparent infinite loop ',
1757 1 50, 303, 2, 0, 0, 0, 0, 0.0e0, 0.0e0)
subroutine ssolsy(WM, IWM, X, TEM)
subroutine sewset(N, ITOL, RTOL, ATOL, YCUR, EWT)
subroutine xerrwd(MSG, NMES, NERR, LEVEL, NI, I1, I2, NR, R1, R2)
F77_RET_T F77_REAL &F77_RET_T F77_DBLE &F77_RET_T F77_REAL &F77_RET_T F77_DBLE &F77_RET_T F77_REAL &F77_RET_T F77_DBLE &F77_RET_T const F77_REAL const F77_REAL F77_REAL &F77_RET_T const F77_DBLE const F77_DBLE F77_DBLE &F77_RET_T F77_REAL &F77_RET_T F77_DBLE &F77_RET_T F77_DBLE &F77_RET_T F77_REAL &F77_RET_T F77_REAL &F77_RET_T F77_DBLE &F77_RET_T const F77_DBLE F77_DBLE &F77_RET_T const F77_REAL F77_REAL &F77_RET_T F77_REAL F77_REAL &F77_RET_T F77_DBLE F77_DBLE &F77_RET_T const F77_DBLE const F77_DBLE * f
subroutine slsode(F, NEQ, Y, T, TOUT, ITOL, RTOL, ATOL, ITASK, ISTATE, IOPT, RWORK, LRW, IWORK, LIW, JAC, MF)
subroutine sintdy(T, K, YH, NYH, DKY, IFLAG)
may be zero for pure relative error test tem the relative tolerance must be greater than or equal to
OCTAVE_EXPORT octave_value_list etc The functions then dimension(columns)
charNDArray max(char d, const charNDArray &m)
OCTAVE_EXPORT octave_value_list return the value of the option it must match the dimension of the state and the relative tolerance must also be a vector of the same length tem it must match the dimension of the state and the absolute tolerance must also be a vector of the same length The local error test applied at each integration step is xample roup abs(local error in x(i))<
subroutine sprepj(NEQ, Y, YH, NYH, EWT, FTEM, SAVF, WM, IWM, F, JAC)
subroutine sstode(NEQ, Y, YH, NYH, YH1, EWT, SAVF, ACOR, WM, IWM, F, JAC, PJAC, SLVS)
octave_value sqrt(void) const
charNDArray min(char d, const charNDArray &m)