2 SUBROUTINE dlsode (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 DOUBLE PRECISION y, t, tout, rtol, atol, rwork
7 dimension neq(*), y(*), rtol(*), atol(*), rwork(lrw), iwork(liw)
1206 EXTERNAL dprepj, dsolsy
1207 DOUBLE PRECISION dumach, dvnorm
1210 INTEGER init, mxstep, mxhnil, nhnil, nslast, nyh, iowns,
1211 1 icf, ierpj, iersl, jcur, jstart, kflag, l,
1212 2 lyh, lewt, lacor, lsavf, lwm, liwm, meth, miter,
1213 3 maxord, maxcor, msbp, mxncf, n, nq, nst, nfe, nje, nqu
1214 INTEGER i, i1, i2, iflag, imxer, kgo, lf0,
1215 1 leniw, lenrw, lenwm, ml, mord, mu, mxhnl0, mxstp0
1216 DOUBLE PRECISION rowns,
1217 1 ccmax, el0, h, hmin, hmxi, hu, rc, tn, uround
1218 DOUBLE PRECISION atoli, ayi, big, ewti, h0, hmax, hmx, rh, rtoli,
1219 1 tcrit, tdist, tnext, tol, tolsf, tp,
SIZE, sum, w0
1223 SAVE mord, mxstp0, mxhnl0
1234 COMMON /dls001/ rowns(209),
1235 1 ccmax, el0, h, hmin, hmxi, hu, rc, tn, uround,
1236 2 init, mxstep, mxhnil, nhnil, nslast, nyh, iowns(6),
1237 3 icf, ierpj, iersl, jcur, jstart, kflag, l,
1238 4 lyh, lewt, lacor, lsavf, lwm, liwm, meth, miter,
1239 5 maxord, maxcor, msbp, mxncf, n, nq, nst, nfe, nje, nqu
1241 DATA mord(1),mord(2)/12,5/, mxstp0/500/, mxhnl0/10/
1252 IF (istate .LT. 1 .OR. istate .GT. 3)
GO TO 601
1253 IF (itask .LT. 1 .OR. itask .GT. 5)
GO TO 602
1254 IF (istate .EQ. 1)
GO TO 10
1255 IF (init .EQ. 0)
GO TO 603
1256 IF (istate .EQ. 2)
GO TO 200
1259 IF (tout .EQ. t)
RETURN 1269 20
IF (neq(1) .LE. 0)
GO TO 604
1270 IF (istate .EQ. 1)
GO TO 25
1271 IF (neq(1) .GT. n)
GO TO 605
1273 IF (itol .LT. 1 .OR. itol .GT. 4)
GO TO 606
1274 IF (iopt .LT. 0 .OR. iopt .GT. 1)
GO TO 607
1276 miter = mf - 10*meth
1277 IF (meth .LT. 1 .OR. meth .GT. 2)
GO TO 608
1278 IF (miter .LT. 0 .OR. miter .GT. 5)
GO TO 608
1279 IF (miter .LE. 3)
GO TO 30
1282 IF (ml .LT. 0 .OR. ml .GE. n)
GO TO 609
1283 IF (mu .LT. 0 .OR. mu .GE. n)
GO TO 610
1286 IF (iopt .EQ. 1)
GO TO 40
1290 IF (istate .EQ. 1) h0 = 0.0d0
1294 40 maxord = iwork(5)
1295 IF (maxord .LT. 0)
GO TO 611
1296 IF (maxord .EQ. 0) maxord = 100
1297 maxord = min(maxord,mord(meth))
1299 IF (mxstep .LT. 0)
GO TO 612
1300 IF (mxstep .EQ. 0) mxstep = mxstp0
1302 IF (mxhnil .LT. 0)
GO TO 613
1303 IF (mxhnil .EQ. 0) mxhnil = mxhnl0
1304 IF (istate .NE. 1)
GO TO 50
1306 IF ((tout - t)*h0 .LT. 0.0d0)
GO TO 614
1308 IF (hmax .LT. 0.0d0)
GO TO 615
1310 IF (hmax .GT. 0.0d0) hmxi = 1.0d0/hmax
1312 IF (hmin .LT. 0.0d0)
GO TO 616
1320 IF (istate .EQ. 1) nyh = n
1321 lwm = lyh + (maxord + 1)*nyh
1322 IF (miter .EQ. 0) lenwm = 0
1323 IF (miter .EQ. 1 .OR. miter .EQ. 2) lenwm = n*n + 2
1324 IF (miter .EQ. 3) lenwm = n + 2
1325 IF (miter .GE. 4) lenwm = (2*ml + mu + 1)*n + 2
1329 lenrw = lacor + n - 1
1333 IF (miter .EQ. 0 .OR. miter .EQ. 3) leniw = 20
1335 IF (lenrw .GT. lrw)
GO TO 617
1336 IF (leniw .GT. liw)
GO TO 618
1341 IF (itol .GE. 3) rtoli = rtol(i)
1342 IF (itol .EQ. 2 .OR. itol .EQ. 4) atoli = atol(i)
1343 IF (rtoli .LT. 0.0d0)
GO TO 619
1344 IF (atoli .LT. 0.0d0)
GO TO 620
1346 IF (istate .EQ. 1)
GO TO 100
1349 IF (nq .LE. maxord)
GO TO 90
1352 80 rwork(i+lsavf-1) = rwork(i+lwm-1)
1354 90
IF (miter .GT. 0) rwork(lwm) = sqrt(uround)
1355 IF (n .EQ. nyh)
GO TO 200
1358 i2 = lyh + (maxord + 1)*nyh - 1
1359 IF (i1 .GT. i2)
GO TO 200
1370 100 uround = dumach()
1372 IF (itask .NE. 4 .AND. itask .NE. 5)
GO TO 110
1374 IF ((tcrit - tout)*(tout - t) .LT. 0.0d0)
GO TO 625
1375 IF (h0 .NE. 0.0d0 .AND. (t + h0 - tcrit)*h0 .GT. 0.0d0)
1378 IF (miter .GT. 0) rwork(lwm) = sqrt(uround)
1391 CALL f (neq, t, y, rwork(lf0))
1395 115 rwork(i+lyh-1) = y(i)
1399 CALL dewset (n, itol, rtol, atol, rwork(lyh), rwork(lewt))
1401 IF (rwork(i+lewt-1) .LE. 0.0d0)
GO TO 621
1402 120 rwork(i+lewt-1) = 1.0d0/rwork(i+lewt-1)
1419 IF (h0 .NE. 0.0d0)
GO TO 180
1420 tdist = abs(tout - t)
1421 w0 = max(abs(t),abs(tout))
1422 IF (tdist .LT. 2.0d0*uround*w0)
GO TO 622
1424 IF (itol .LE. 2)
GO TO 140
1426 130 tol = max(tol,rtol(i))
1427 140
IF (tol .GT. 0.0d0)
GO TO 160
1430 IF (itol .EQ. 2 .OR. itol .EQ. 4) atoli = atol(i)
1432 IF (ayi .NE. 0.0d0) tol = max(tol,atoli/ayi)
1434 160 tol = max(tol,100.0d0*uround)
1435 tol = min(tol,0.001d0)
1436 sum = dvnorm(n, rwork(lf0), rwork(lewt))
1437 sum = 1.0d0/(tol*w0*w0) + tol*sum**2
1438 h0 = 1.0d0/sqrt(sum)
1440 h0 = sign(h0,tout-t)
1442 180 rh = abs(h0)*hmxi
1443 IF (rh .GT. 1.0d0) h0 = h0/rh
1447 190 rwork(i+lf0-1) = h0*rwork(i+lf0-1)
1455 GO TO (210, 250, 220, 230, 240), itask
1456 210
IF ((tn - tout)*h .LT. 0.0d0)
GO TO 250
1457 CALL dintdy (tout, 0, rwork(lyh), nyh, y, iflag)
1458 IF (iflag .NE. 0)
GO TO 627
1461 220 tp = tn - hu*(1.0d0 + 100.0d0*uround)
1462 IF ((tp - tout)*h .GT. 0.0d0)
GO TO 623
1463 IF ((tn - tout)*h .LT. 0.0d0)
GO TO 250
1465 230 tcrit = rwork(1)
1466 IF ((tn - tcrit)*h .GT. 0.0d0)
GO TO 624
1467 IF ((tcrit - tout)*h .LT. 0.0d0)
GO TO 625
1468 IF ((tn - tout)*h .LT. 0.0d0)
GO TO 245
1469 CALL dintdy (tout, 0, rwork(lyh), nyh, y, iflag)
1470 IF (iflag .NE. 0)
GO TO 627
1473 240 tcrit = rwork(1)
1474 IF ((tn - tcrit)*h .GT. 0.0d0)
GO TO 624
1475 245 hmx = abs(tn) + abs(h)
1476 ihit = abs(tn - tcrit) .LE. 100.0d0*uround*hmx
1478 tnext = tn + h*(1.0d0 + 4.0d0*uround)
1479 IF ((tnext - tcrit)*h .LE. 0.0d0)
GO TO 250
1480 h = (tcrit - tn)*(1.0d0 - 4.0d0*uround)
1481 IF (istate .EQ. 2) jstart = -2
1494 IF ((nst-nslast) .GE. mxstep)
GO TO 500
1495 CALL dewset (n, itol, rtol, atol, rwork(lyh), rwork(lewt))
1497 IF (rwork(i+lewt-1) .LE. 0.0d0)
GO TO 510
1498 260 rwork(i+lewt-1) = 1.0d0/rwork(i+lewt-1)
1499 270 tolsf = uround*dvnorm(n, rwork(lyh), rwork(lewt))
1500 IF (tolsf .LE. 1.0d0)
GO TO 280
1502 IF (nst .EQ. 0)
GO TO 626
1504 280
IF ((tn + h) .NE. tn)
GO TO 290
1506 IF (nhnil .GT. mxhnil)
GO TO 290
1507 msg =
'DLSODE- Warning..internal T (=R1) and H (=R2) are' 1508 CALL xerrwd (msg, 50, 101, 0, 0, 0, 0, 0, 0.0d0, 0.0d0)
1509 msg=
' such that in the machine, T + H = T on the next step ' 1510 CALL xerrwd (msg, 60, 101, 0, 0, 0, 0, 0, 0.0d0, 0.0d0)
1511 msg =
' (H = step size). Solver will continue anyway' 1512 CALL xerrwd (msg, 50, 101, 0, 0, 0, 0, 2, tn, h)
1513 IF (nhnil .LT. mxhnil)
GO TO 290
1514 msg =
'DLSODE- Above warning has been issued I1 times. ' 1515 CALL xerrwd (msg, 50, 102, 0, 0, 0, 0, 0, 0.0d0, 0.0d0)
1516 msg =
' It will not be issued again for this problem' 1517 CALL xerrwd (msg, 50, 102, 0, 1, mxhnil, 0, 0, 0.0d0, 0.0d0)
1522 CALL dstode (neq, y, rwork(lyh), nyh, rwork(lyh), rwork(lewt),
1523 1 rwork(lsavf), rwork(lacor), rwork(lwm), iwork(liwm),
1524 2 f, jac, dprepj, dsolsy)
1526 GO TO (300, 530, 540), kgo
1533 GO TO (310, 400, 330, 340, 350), itask
1535 310
IF ((tn - tout)*h .LT. 0.0d0)
GO TO 250
1536 CALL dintdy (tout, 0, rwork(lyh), nyh, y, iflag)
1540 330
IF ((tn - tout)*h .GE. 0.0d0)
GO TO 400
1543 340
IF ((tn - tout)*h .LT. 0.0d0)
GO TO 345
1544 CALL dintdy (tout, 0, rwork(lyh), nyh, y, iflag)
1547 345 hmx = abs(tn) + abs(h)
1548 ihit = abs(tn - tcrit) .LE. 100.0d0*uround*hmx
1550 tnext = tn + h*(1.0d0 + 4.0d0*uround)
1551 IF ((tnext - tcrit)*h .LE. 0.0d0)
GO TO 250
1552 h = (tcrit - tn)*(1.0d0 - 4.0d0*uround)
1556 350 hmx = abs(tn) + abs(h)
1557 ihit = abs(tn - tcrit) .LE. 100.0d0*uround*hmx
1566 410 y(i) = rwork(i+lyh-1)
1568 IF (itask .NE. 4 .AND. itask .NE. 5)
GO TO 420
1589 500 msg =
'DLSODE- At current T (=R1), MXSTEP (=I1) steps ' 1590 CALL xerrwd (msg, 50, 201, 0, 0, 0, 0, 0, 0.0d0, 0.0d0)
1591 msg =
' taken on this call before reaching TOUT ' 1592 CALL xerrwd (msg, 50, 201, 0, 1, mxstep, 0, 1, tn, 0.0d0)
1596 510 ewti = rwork(lewt+i-1)
1597 msg = .LE.
'DLSODE- At T (=R1), EWT(I1) has become R2 0.' 1598 CALL xerrwd (msg, 50, 202, 0, 1, i, 0, 2, tn, ewti)
1602 520 msg =
'DLSODE- At T (=R1), too much accuracy requested ' 1603 CALL xerrwd (msg, 50, 203, 0, 0, 0, 0, 0, 0.0d0, 0.0d0)
1604 msg =
' for precision of machine.. see TOLSF (=R2) ' 1605 CALL xerrwd (msg, 50, 203, 0, 0, 0, 0, 2, tn, tolsf)
1610 530 msg =
'DLSODE- At T(=R1) and step size H(=R2), the error' 1611 CALL xerrwd (msg, 50, 204, 0, 0, 0, 0, 0, 0.0d0, 0.0d0)
1612 msg =
' test failed repeatedly or with ABS(H) = HMIN' 1613 CALL xerrwd (msg, 50, 204, 0, 0, 0, 0, 2, tn, h)
1617 540 msg =
'DLSODE- At T (=R1) and step size H (=R2), the ' 1618 CALL xerrwd (msg, 50, 205, 0, 0, 0, 0, 0, 0.0d0, 0.0d0)
1619 msg =
' corrector convergence failed repeatedly ' 1620 CALL xerrwd (msg, 50, 205, 0, 0, 0, 0, 0, 0.0d0, 0.0d0)
1621 msg =
' or with ABS(H) = HMIN ' 1622 CALL xerrwd (msg, 30, 205, 0, 0, 0, 0, 2, tn, h)
1628 SIZE = abs(rwork(i+lacor-1)*rwork(i+lewt-1))
1629 IF (big .GE. size)
GO TO 570
1636 590 y(i) = rwork(i+lyh-1)
1654 601 msg =
'DLSODE- ISTATE (=I1) illegal ' 1655 CALL xerrwd (msg, 30, 1, 0, 1, istate, 0, 0, 0.0d0, 0.0d0)
1656 IF (istate .LT. 0)
GO TO 800
1658 602 msg =
'DLSODE- ITASK (=I1) illegal ' 1659 CALL xerrwd (msg, 30, 2, 0, 1, itask, 0, 0, 0.0d0, 0.0d0)
1661 603 msg = .GT.
'DLSODE- ISTATE 1 but DLSODE not initialized ' 1662 CALL xerrwd (msg, 50, 3, 0, 0, 0, 0, 0, 0.0d0, 0.0d0)
1664 604 msg = .LT.
'DLSODE- NEQ (=I1) 1 ' 1665 CALL xerrwd (msg, 30, 4, 0, 1, neq(1), 0, 0, 0.0d0, 0.0d0)
1667 605 msg =
'DLSODE- ISTATE = 3 and NEQ increased (I1 to I2) ' 1668 CALL xerrwd (msg, 50, 5, 0, 2, n, neq(1), 0, 0.0d0, 0.0d0)
1670 606 msg =
'DLSODE- ITOL (=I1) illegal ' 1671 CALL xerrwd (msg, 30, 6, 0, 1, itol, 0, 0, 0.0d0, 0.0d0)
1673 607 msg =
'DLSODE- IOPT (=I1) illegal ' 1674 CALL xerrwd (msg, 30, 7, 0, 1, iopt, 0, 0, 0.0d0, 0.0d0)
1676 608 msg =
'DLSODE- MF (=I1) illegal ' 1677 CALL xerrwd (msg, 30, 8, 0, 1, mf, 0, 0, 0.0d0, 0.0d0)
1679 609 msg = .LT..GE.
'DLSODE- ML (=I1) illegal.. 0 or NEQ (=I2)' 1680 CALL xerrwd (msg, 50, 9, 0, 2, ml, neq(1), 0, 0.0d0, 0.0d0)
1682 610 msg = .LT..GE.
'DLSODE- MU (=I1) illegal.. 0 or NEQ (=I2)' 1683 CALL xerrwd (msg, 50, 10, 0, 2, mu, neq(1), 0, 0.0d0, 0.0d0)
1685 611 msg = .LT.
'DLSODE- MAXORD (=I1) 0 ' 1686 CALL xerrwd (msg, 30, 11, 0, 1, maxord, 0, 0, 0.0d0, 0.0d0)
1688 612 msg = .LT.
'DLSODE- MXSTEP (=I1) 0 ' 1689 CALL xerrwd (msg, 30, 12, 0, 1, mxstep, 0, 0, 0.0d0, 0.0d0)
1691 613 msg = .LT.
'DLSODE- MXHNIL (=I1) 0 ' 1692 CALL xerrwd (msg, 30, 13, 0, 1, mxhnil, 0, 0, 0.0d0, 0.0d0)
1694 614 msg =
'DLSODE- TOUT (=R1) behind T (=R2) ' 1695 CALL xerrwd (msg, 40, 14, 0, 0, 0, 0, 2, tout, t)
1696 msg =
' Integration direction is given by H0 (=R1) ' 1697 CALL xerrwd (msg, 50, 14, 0, 0, 0, 0, 1, h0, 0.0d0)
1699 615 msg = .LT.
'DLSODE- HMAX (=R1) 0.0 ' 1700 CALL xerrwd (msg, 30, 15, 0, 0, 0, 0, 1, hmax, 0.0d0)
1702 616 msg = .LT.
'DLSODE- HMIN (=R1) 0.0 ' 1703 CALL xerrwd (msg, 30, 16, 0, 0, 0, 0, 1, hmin, 0.0d0)
1706 msg=
'DLSODE- RWORK length needed, LENRW (=I1), exceeds LRW (=I2)' 1707 CALL xerrwd (msg, 60, 17, 0, 2, lenrw, lrw, 0, 0.0d0, 0.0d0)
1710 msg=
'DLSODE- IWORK length needed, LENIW (=I1), exceeds LIW (=I2)' 1711 CALL xerrwd (msg, 60, 18, 0, 2, leniw, liw, 0, 0.0d0, 0.0d0)
1713 619 msg = .LT.
'DLSODE- RTOL(I1) is R1 0.0 ' 1714 CALL xerrwd (msg, 40, 19, 0, 1, i, 0, 1, rtoli, 0.0d0)
1716 620 msg = .LT.
'DLSODE- ATOL(I1) is R1 0.0 ' 1717 CALL xerrwd (msg, 40, 20, 0, 1, i, 0, 1, atoli, 0.0d0)
1719 621 ewti = rwork(lewt+i-1)
1720 msg = .LE.
'DLSODE- EWT(I1) is R1 0.0 ' 1721 CALL xerrwd (msg, 40, 21, 0, 1, i, 0, 1, ewti, 0.0d0)
1724 msg=
'DLSODE- TOUT (=R1) too close to T(=R2) to start integration' 1725 CALL xerrwd (msg, 60, 22, 0, 0, 0, 0, 2, tout, t)
1728 msg=
'DLSODE- ITASK = I1 and TOUT (=R1) behind TCUR - HU (= R2) ' 1729 CALL xerrwd (msg, 60, 23, 0, 1, itask, 0, 2, tout, tp)
1732 msg=
'DLSODE- ITASK = 4 OR 5 and TCRIT (=R1) behind TCUR (=R2) ' 1733 CALL xerrwd (msg, 60, 24, 0, 0, 0, 0, 2, tcrit, tn)
1736 msg=
'DLSODE- ITASK = 4 or 5 and TCRIT (=R1) behind TOUT (=R2) ' 1737 CALL xerrwd (msg, 60, 25, 0, 0, 0, 0, 2, tcrit, tout)
1739 626 msg =
'DLSODE- At start of problem, too much accuracy ' 1740 CALL xerrwd (msg, 50, 26, 0, 0, 0, 0, 0, 0.0d0, 0.0d0)
1741 msg=
' requested for precision of machine.. See TOLSF (=R1) ' 1742 CALL xerrwd (msg, 60, 26, 0, 0, 0, 0, 1, tolsf, 0.0d0)
1745 627 msg =
'DLSODE- Trouble in DINTDY. ITASK = I1, TOUT = R1' 1746 CALL xerrwd (msg, 50, 27, 0, 1, itask, 0, 1, tout, 0.0d0)
1751 800 msg =
'DLSODE- Run aborted.. apparent infinite loop ' 1752 CALL xerrwd (msg, 50, 303, 2, 0, 0, 0, 0, 0.0d0, 0.0d0)
1757 SUBROUTINE dlsodes (F, NEQ, Y, T, TOUT, ITOL, RTOL, ATOL, ITASK,
1758 1 ISTATE, IOPT, RWORK, LRW, IWORK, LIW, JAC, MF)
1760 INTEGER neq, itol, itask, istate, iopt, lrw, iwork, liw, mf
1761 DOUBLE PRECISION y, t, tout, rtol, atol, rwork
1762 dimension neq(*), y(*), rtol(*), atol(*), rwork(lrw), iwork(liw)
3006 EXTERNAL dprjs, dsolss
3007 DOUBLE PRECISION dumach, dvnorm
3008 INTEGER init, mxstep, mxhnil, nhnil, nslast, nyh, iowns,
3009 1 icf, ierpj, iersl, jcur, jstart, kflag, l,
3010 2 lyh, lewt, lacor, lsavf, lwm, liwm, meth, miter,
3011 3 maxord, maxcor, msbp, mxncf, n, nq, nst, nfe, nje, nqu
3012 INTEGER iplost, iesp, istatc, iys, iba, ibian, ibjan, ibjgp,
3013 1 ipian, ipjan, ipjgp, ipigp, ipr, ipc, ipic, ipisp, iprsp, ipa,
3014 2 lenyh, lenyhm, lenwk, lreq, lrat, lrest, lwmin, moss, msbj,
3015 3 nslj, ngp, nlu, nnz, nsp, nzl, nzu
3016 INTEGER i, i1, i2, iflag, imax, imul, imxer, ipflag, ipgo, irem,
3017 1 j, kgo, lenrat, lenyht, leniw, lenrw, lf0, lia, lja,
3018 2 lrtem, lwtem, lyhd, lyhn, mf1, mord, mxhnl0, mxstp0, ncolm
3019 DOUBLE PRECISION rowns,
3020 1 ccmax, el0, h, hmin, hmxi, hu, rc, tn, uround
3021 DOUBLE PRECISION con0, conmin, ccmxj, psmall, rbig, seth
3022 DOUBLE PRECISION atoli, ayi, big, ewti, h0, hmax, hmx, rh, rtoli,
3023 1 tcrit, tdist, tnext, tol, tolsf, tp,
SIZE, sum, w0
3027 SAVE lenrat, mord, mxstp0, mxhnl0
3040 COMMON /dls001/ rowns(209),
3041 1 ccmax, el0, h, hmin, hmxi, hu, rc, tn, uround,
3042 2 init, mxstep, mxhnil, nhnil, nslast, nyh, iowns(6),
3043 3 icf, ierpj, iersl, jcur, jstart, kflag, l,
3044 4 lyh, lewt, lacor, lsavf, lwm, liwm, meth, miter,
3045 5 maxord, maxcor, msbp, mxncf, n, nq, nst, nfe, nje, nqu
3047 COMMON /dlss01/ con0, conmin, ccmxj, psmall, rbig, seth,
3048 1 iplost, iesp, istatc, iys, iba, ibian, ibjan, ibjgp,
3049 2 ipian, ipjan, ipjgp, ipigp, ipr, ipc, ipic, ipisp, iprsp, ipa,
3050 3 lenyh, lenyhm, lenwk, lreq, lrat, lrest, lwmin, moss, msbj,
3051 4 nslj, ngp, nlu, nnz, nsp, nzl, nzu
3053 DATA mord(1),mord(2)/12,5/, mxstp0/500/, mxhnl0/10/
3069 IF (istate .LT. 1 .OR. istate .GT. 3)
GO TO 601
3070 IF (itask .LT. 1 .OR. itask .GT. 5)
GO TO 602
3071 IF (istate .EQ. 1)
GO TO 10
3072 IF (init .EQ. 0)
GO TO 603
3073 IF (istate .EQ. 2)
GO TO 200
3076 IF (tout .EQ. t)
RETURN 3088 20
IF (neq(1) .LE. 0)
GO TO 604
3089 IF (istate .EQ. 1)
GO TO 25
3090 IF (neq(1) .GT. n)
GO TO 605
3092 IF (itol .LT. 1 .OR. itol .GT. 4)
GO TO 606
3093 IF (iopt .LT. 0 .OR. iopt .GT. 1)
GO TO 607
3097 miter = mf1 - 10*meth
3098 IF (moss .LT. 0 .OR. moss .GT. 2)
GO TO 608
3099 IF (meth .LT. 1 .OR. meth .GT. 2)
GO TO 608
3100 IF (miter .LT. 0 .OR. miter .GT. 3)
GO TO 608
3101 IF (miter .EQ. 0 .OR. miter .EQ. 3) moss = 0
3103 IF (iopt .EQ. 1)
GO TO 40
3107 IF (istate .EQ. 1) h0 = 0.0d0
3112 40 maxord = iwork(5)
3113 IF (maxord .LT. 0)
GO TO 611
3114 IF (maxord .EQ. 0) maxord = 100
3115 maxord = min(maxord,mord(meth))
3117 IF (mxstep .LT. 0)
GO TO 612
3118 IF (mxstep .EQ. 0) mxstep = mxstp0
3120 IF (mxhnil .LT. 0)
GO TO 613
3121 IF (mxhnil .EQ. 0) mxhnil = mxhnl0
3122 IF (istate .NE. 1)
GO TO 50
3124 IF ((tout - t)*h0 .LT. 0.0d0)
GO TO 614
3126 IF (hmax .LT. 0.0d0)
GO TO 615
3128 IF (hmax .GT. 0.0d0) hmxi = 1.0d0/hmax
3130 IF (hmin .LT. 0.0d0)
GO TO 616
3132 IF (seth .LT. 0.0d0)
GO TO 609
3137 IF (itol .GE. 3) rtoli = rtol(i)
3138 IF (itol .EQ. 2 .OR. itol .EQ. 4) atoli = atol(i)
3139 IF (rtoli .LT. 0.0d0)
GO TO 619
3140 IF (atoli .LT. 0.0d0)
GO TO 620
3159 IF (istate .EQ. 1) nyh = n
3161 IF (miter .EQ. 1) lwmin = 4*n + 10*n/lrat
3162 IF (miter .EQ. 2) lwmin = 4*n + 11*n/lrat
3163 IF (miter .EQ. 3) lwmin = n + 2
3164 lenyh = (maxord+1)*nyh
3166 lenrw = 20 + lwmin + lrest
3169 IF (moss .EQ. 0 .AND. miter .NE. 0 .AND. miter .NE. 3)
3170 1 leniw = leniw + n + 1
3172 IF (lenrw .GT. lrw)
GO TO 617
3173 IF (leniw .GT. liw)
GO TO 618
3175 IF (moss .EQ. 0 .AND. miter .NE. 0 .AND. miter .NE. 3)
3176 1 leniw = leniw + iwork(lia+n) - 1
3178 IF (leniw .GT. liw)
GO TO 618
3183 IF (istate .EQ. 1) nq = 1
3184 ncolm = min(nq+1,maxord+2)
3187 IF (miter .EQ. 1 .OR. miter .EQ. 2) lenyht = lenyhm
3189 IF (istate .EQ. 3) imul = moss
3190 IF (moss .EQ. 2) imul = 3
3191 lrtem = lenyht + imul*n
3193 IF (miter .EQ. 1 .OR. miter .EQ. 2) lwtem = lrw - 20 - lrtem
3196 lsavf = lyhn + lenyht
3200 IF (istate .EQ. 1)
GO TO 100
3211 imax = lyhn - 1 + lenyhm
3213 IF (lyhd .LT. 0)
THEN 3216 72 rwork(j) = rwork(j+lyhd)
3218 IF (lyhd .GT. 0)
THEN 3220 76 rwork(i) = rwork(i+lyhd)
3224 IF (miter .EQ. 0 .OR. miter .EQ. 3)
GO TO 92
3225 IF (moss .NE. 2)
GO TO 85
3227 CALL dewset (n, itol, rtol, atol, rwork(lyh), rwork(lewt))
3229 IF (rwork(i+lewt-1) .LE. 0.0d0)
GO TO 621
3230 82 rwork(i+lewt-1) = 1.0d0/rwork(i+lewt-1)
3233 lsavf = min(lsavf,lrw)
3234 lewt = min(lewt,lrw)
3235 lacor = min(lacor,lrw)
3236 CALL diprep (neq, y, rwork, iwork(lia),iwork(lja), ipflag, f, jac)
3237 lenrw = lwm - 1 + lenwk + lrest
3239 IF (ipflag .NE. -1) iwork(23) = ipian
3240 IF (ipflag .NE. -1) iwork(24) = ipjan
3242 GO TO (90, 628, 629, 630, 631, 632, 633), ipgo
3244 IF (lenrw .GT. lrw)
GO TO 617
3247 IF (n .EQ. nyh)
GO TO 200
3250 i2 = lyh + (maxord + 1)*nyh - 1
3251 IF (i1 .GT. i2)
GO TO 200
3275 105 rwork(i+lyh-1) = y(i)
3278 CALL f (neq, t, y, rwork(lf0))
3281 CALL dewset (n, itol, rtol, atol, rwork(lyh), rwork(lewt))
3283 IF (rwork(i+lewt-1) .LE. 0.0d0)
GO TO 621
3284 110 rwork(i+lewt-1) = 1.0d0/rwork(i+lewt-1)
3285 IF (miter .EQ. 0 .OR. miter .EQ. 3)
GO TO 120
3287 lacor = min(lacor,lrw)
3288 CALL diprep (neq, y, rwork, iwork(lia),iwork(lja), ipflag, f, jac)
3289 lenrw = lwm - 1 + lenwk + lrest
3291 IF (ipflag .NE. -1) iwork(23) = ipian
3292 IF (ipflag .NE. -1) iwork(24) = ipjan
3294 GO TO (115, 628, 629, 630, 631, 632, 633), ipgo
3296 IF (lenrw .GT. lrw)
GO TO 617
3299 IF (itask .NE. 4 .AND. itask .NE. 5)
GO TO 125
3301 IF ((tcrit - tout)*(tout - t) .LT. 0.0d0)
GO TO 625
3302 IF (h0 .NE. 0.0d0 .AND. (t + h0 - tcrit)*h0 .GT. 0.0d0)
3305 125 uround = dumach()
3307 IF (miter .NE. 0) rwork(lwm) = sqrt(uround)
3311 psmall = 1000.0d0*uround
3312 rbig = 0.01d0/psmall
3341 IF (h0 .NE. 0.0d0)
GO TO 180
3342 tdist = abs(tout - t)
3343 w0 = max(abs(t),abs(tout))
3344 IF (tdist .LT. 2.0d0*uround*w0)
GO TO 622
3346 IF (itol .LE. 2)
GO TO 140
3348 130 tol = max(tol,rtol(i))
3349 140
IF (tol .GT. 0.0d0)
GO TO 160
3352 IF (itol .EQ. 2 .OR. itol .EQ. 4) atoli = atol(i)
3354 IF (ayi .NE. 0.0d0) tol = max(tol,atoli/ayi)
3356 160 tol = max(tol,100.0d0*uround)
3357 tol = min(tol,0.001d0)
3358 sum = dvnorm(n, rwork(lf0), rwork(lewt))
3359 sum = 1.0d0/(tol*w0*w0) + tol*sum**2
3360 h0 = 1.0d0/sqrt(sum)
3362 h0 = sign(h0,tout-t)
3364 180 rh = abs(h0)*hmxi
3365 IF (rh .GT. 1.0d0) h0 = h0/rh
3369 190 rwork(i+lf0-1) = h0*rwork(i+lf0-1)
3377 GO TO (210, 250, 220, 230, 240), itask
3378 210
IF ((tn - tout)*h .LT. 0.0d0)
GO TO 250
3379 CALL dintdy (tout, 0, rwork(lyh), nyh, y, iflag)
3380 IF (iflag .NE. 0)
GO TO 627
3383 220 tp = tn - hu*(1.0d0 + 100.0d0*uround)
3384 IF ((tp - tout)*h .GT. 0.0d0)
GO TO 623
3385 IF ((tn - tout)*h .LT. 0.0d0)
GO TO 250
3387 230 tcrit = rwork(1)
3388 IF ((tn - tcrit)*h .GT. 0.0d0)
GO TO 624
3389 IF ((tcrit - tout)*h .LT. 0.0d0)
GO TO 625
3390 IF ((tn - tout)*h .LT. 0.0d0)
GO TO 245
3391 CALL dintdy (tout, 0, rwork(lyh), nyh, y, iflag)
3392 IF (iflag .NE. 0)
GO TO 627
3395 240 tcrit = rwork(1)
3396 IF ((tn - tcrit)*h .GT. 0.0d0)
GO TO 624
3397 245 hmx = abs(tn) + abs(h)
3398 ihit = abs(tn - tcrit) .LE. 100.0d0*uround*hmx
3400 tnext = tn + h*(1.0d0 + 4.0d0*uround)
3401 IF ((tnext - tcrit)*h .LE. 0.0d0)
GO TO 250
3402 h = (tcrit - tn)*(1.0d0 - 4.0d0*uround)
3403 IF (istate .EQ. 2) jstart = -2
3416 IF ((nst-nslast) .GE. mxstep)
GO TO 500
3417 CALL dewset (n, itol, rtol, atol, rwork(lyh), rwork(lewt))
3419 IF (rwork(i+lewt-1) .LE. 0.0d0)
GO TO 510
3420 260 rwork(i+lewt-1) = 1.0d0/rwork(i+lewt-1)
3421 270 tolsf = uround*dvnorm(n, rwork(lyh), rwork(lewt))
3422 IF (tolsf .LE. 1.0d0)
GO TO 280
3424 IF (nst .EQ. 0)
GO TO 626
3426 280
IF ((tn + h) .NE. tn)
GO TO 290
3428 IF (nhnil .GT. mxhnil)
GO TO 290
3429 msg =
'DLSODES- Warning..Internal T (=R1) and H (=R2) are' 3430 CALL xerrwd (msg, 50, 101, 0, 0, 0, 0, 0, 0.0d0, 0.0d0)
3431 msg=
' such that in the machine, T + H = T on the next step ' 3432 CALL xerrwd (msg, 60, 101, 0, 0, 0, 0, 0, 0.0d0, 0.0d0)
3433 msg =
' (H = step size). Solver will continue anyway.' 3434 CALL xerrwd (msg, 50, 101, 0, 0, 0, 0, 2, tn, h)
3435 IF (nhnil .LT. mxhnil)
GO TO 290
3436 msg =
'DLSODES- Above warning has been issued I1 times. ' 3437 CALL xerrwd (msg, 50, 102, 0, 0, 0, 0, 0, 0.0d0, 0.0d0)
3438 msg =
' It will not be issued again for this problem.' 3439 CALL xerrwd (msg, 50, 102, 0, 1, mxhnil, 0, 0, 0.0d0, 0.0d0)
3444 CALL dstode (neq, y, rwork(lyh), nyh, rwork(lyh), rwork(lewt),
3445 1 rwork(lsavf), rwork(lacor), rwork(lwm), rwork(lwm),
3446 2 f, jac, dprjs, dsolss)
3448 GO TO (300, 530, 540, 550), kgo
3455 GO TO (310, 400, 330, 340, 350), itask
3457 310
IF ((tn - tout)*h .LT. 0.0d0)
GO TO 250
3458 CALL dintdy (tout, 0, rwork(lyh), nyh, y, iflag)
3462 330
IF ((tn - tout)*h .GE. 0.0d0)
GO TO 400
3465 340
IF ((tn - tout)*h .LT. 0.0d0)
GO TO 345
3466 CALL dintdy (tout, 0, rwork(lyh), nyh, y, iflag)
3469 345 hmx = abs(tn) + abs(h)
3470 ihit = abs(tn - tcrit) .LE. 100.0d0*uround*hmx
3472 tnext = tn + h*(1.0d0 + 4.0d0*uround)
3473 IF ((tnext - tcrit)*h .LE. 0.0d0)
GO TO 250
3474 h = (tcrit - tn)*(1.0d0 - 4.0d0*uround)
3478 350 hmx = abs(tn) + abs(h)
3479 ihit = abs(tn - tcrit) .LE. 100.0d0*uround*hmx
3488 410 y(i) = rwork(i+lyh-1)
3490 IF (itask .NE. 4 .AND. itask .NE. 5)
GO TO 420
3516 500 msg =
'DLSODES- At current T (=R1), MXSTEP (=I1) steps ' 3517 CALL xerrwd (msg, 50, 201, 0, 0, 0, 0, 0, 0.0d0, 0.0d0)
3518 msg =
' taken on this call before reaching TOUT ' 3519 CALL xerrwd (msg, 50, 201, 0, 1, mxstep, 0, 1, tn, 0.0d0)
3523 510 ewti = rwork(lewt+i-1)
3524 msg = .le.
'DLSODES- At T (=R1), EWT(I1) has become R2 0.' 3525 CALL xerrwd (msg, 50, 202, 0, 1, i, 0, 2, tn, ewti)
3529 520 msg =
'DLSODES- At T (=R1), too much accuracy requested ' 3530 CALL xerrwd (msg, 50, 203, 0, 0, 0, 0, 0, 0.0d0, 0.0d0)
3531 msg =
' for precision of machine.. See TOLSF (=R2) ' 3532 CALL xerrwd (msg, 50, 203, 0, 0, 0, 0, 2, tn, tolsf)
3537 530 msg =
'DLSODES- At T(=R1) and step size H(=R2), the error' 3538 CALL xerrwd (msg, 50, 204, 0, 0, 0, 0, 0, 0.0d0, 0.0d0)
3539 msg =
' test failed repeatedly or with ABS(H) = HMIN' 3540 CALL xerrwd (msg, 50, 204, 0, 0, 0, 0, 2, tn, h)
3544 540 msg =
'DLSODES- At T (=R1) and step size H (=R2), the ' 3545 CALL xerrwd (msg, 50, 205, 0, 0, 0, 0, 0, 0.0d0, 0.0d0)
3546 msg =
' corrector convergence failed repeatedly ' 3547 CALL xerrwd (msg, 50, 205, 0, 0, 0, 0, 0, 0.0d0, 0.0d0)
3548 msg =
' or with ABS(H) = HMIN ' 3549 CALL xerrwd (msg, 30, 205, 0, 0, 0, 0, 2, tn, h)
3553 550 msg =
'DLSODES- At T (=R1) and step size H (=R2), a fatal' 3554 CALL xerrwd (msg, 50, 207, 0, 0, 0, 0, 0, 0.0d0, 0.0d0)
3555 msg =
' error flag was returned by CDRV (by way of ' 3556 CALL xerrwd (msg, 50, 207, 0, 0, 0, 0, 0, 0.0d0, 0.0d0)
3557 msg =
' Subroutine DPRJS or DSOLSS) ' 3558 CALL xerrwd (msg, 40, 207, 0, 0, 0, 0, 2, tn, h)
3565 SIZE = abs(rwork(i+lacor-1)*rwork(i+lewt-1))
3566 IF (big .GE. size)
GO TO 570
3573 590 y(i) = rwork(i+lyh-1)
3596 601 msg =
'DLSODES- ISTATE (=I1) illegal.' 3597 CALL xerrwd (msg, 30, 1, 0, 1, istate, 0, 0, 0.0d0, 0.0d0)
3598 IF (istate .LT. 0)
GO TO 800
3600 602 msg =
'DLSODES- ITASK (=I1) illegal. ' 3601 CALL xerrwd (msg, 30, 2, 0, 1, itask, 0, 0, 0.0d0, 0.0d0)
3603 603 msg = .gt.
'DLSODES- ISTATE1 but DLSODES not initialized. ' 3604 CALL xerrwd (msg, 50, 3, 0, 0, 0, 0, 0, 0.0d0, 0.0d0)
3606 604 msg = .lt.
'DLSODES- NEQ (=I1) 1 ' 3607 CALL xerrwd (msg, 30, 4, 0, 1, neq(1), 0, 0, 0.0d0, 0.0d0)
3609 605 msg =
'DLSODES- ISTATE = 3 and NEQ increased (I1 to I2). ' 3610 CALL xerrwd (msg, 50, 5, 0, 2, n, neq(1), 0, 0.0d0, 0.0d0)
3612 606 msg =
'DLSODES- ITOL (=I1) illegal. ' 3613 CALL xerrwd (msg, 30, 6, 0, 1, itol, 0, 0, 0.0d0, 0.0d0)
3615 607 msg =
'DLSODES- IOPT (=I1) illegal. ' 3616 CALL xerrwd (msg, 30, 7, 0, 1, iopt, 0, 0, 0.0d0, 0.0d0)
3618 608 msg =
'DLSODES- MF (=I1) illegal. ' 3619 CALL xerrwd (msg, 30, 8, 0, 1, mf, 0, 0, 0.0d0, 0.0d0)
3621 609 msg = .lt.
'DLSODES- SETH (=R1) 0.0 ' 3622 CALL xerrwd (msg, 30, 9, 0, 0, 0, 0, 1, seth, 0.0d0)
3624 611 msg = .lt.
'DLSODES- MAXORD (=I1) 0 ' 3625 CALL xerrwd (msg, 30, 11, 0, 1, maxord, 0, 0, 0.0d0, 0.0d0)
3627 612 msg = .lt.
'DLSODES- MXSTEP (=I1) 0 ' 3628 CALL xerrwd (msg, 30, 12, 0, 1, mxstep, 0, 0, 0.0d0, 0.0d0)
3630 613 msg = .lt.
'DLSODES- MXHNIL (=I1) 0 ' 3631 CALL xerrwd (msg, 30, 13, 0, 1, mxhnil, 0, 0, 0.0d0, 0.0d0)
3633 614 msg =
'DLSODES- TOUT (=R1) behind T (=R2) ' 3634 CALL xerrwd (msg, 40, 14, 0, 0, 0, 0, 2, tout, t)
3635 msg =
' Integration direction is given by H0 (=R1) ' 3636 CALL xerrwd (msg, 50, 14, 0, 0, 0, 0, 1, h0, 0.0d0)
3638 615 msg = .lt.
'DLSODES- HMAX (=R1) 0.0 ' 3639 CALL xerrwd (msg, 30, 15, 0, 0, 0, 0, 1, hmax, 0.0d0)
3641 616 msg = .lt.
'DLSODES- HMIN (=R1) 0.0 ' 3642 CALL xerrwd (msg, 30, 16, 0, 0, 0, 0, 1, hmin, 0.0d0)
3644 617 msg =
'DLSODES- RWORK length is insufficient to proceed. ' 3645 CALL xerrwd (msg, 50, 17, 0, 0, 0, 0, 0, 0.0d0, 0.0d0)
3646 msg=.ge.
' Length needed is LENRW (=I1), exceeds LRW (=I2)' 3647 CALL xerrwd (msg, 60, 17, 0, 2, lenrw, lrw, 0, 0.0d0, 0.0d0)
3649 618 msg =
'DLSODES- IWORK length is insufficient to proceed. ' 3650 CALL xerrwd (msg, 50, 18, 0, 0, 0, 0, 0, 0.0d0, 0.0d0)
3651 msg=.ge.
' Length needed is LENIW (=I1), exceeds LIW (=I2)' 3652 CALL xerrwd (msg, 60, 18, 0, 2, leniw, liw, 0, 0.0d0, 0.0d0)
3654 619 msg = .lt.
'DLSODES- RTOL(I1) is R1 0.0 ' 3655 CALL xerrwd (msg, 40, 19, 0, 1, i, 0, 1, rtoli, 0.0d0)
3657 620 msg = .lt.
'DLSODES- ATOL(I1) is R1 0.0 ' 3658 CALL xerrwd (msg, 40, 20, 0, 1, i, 0, 1, atoli, 0.0d0)
3660 621 ewti = rwork(lewt+i-1)
3661 msg = .le.
'DLSODES- EWT(I1) is R1 0.0 ' 3662 CALL xerrwd (msg, 40, 21, 0, 1, i, 0, 1, ewti, 0.0d0)
3664 622 msg=
'DLSODES- TOUT(=R1) too close to T(=R2) to start integration.' 3665 CALL xerrwd (msg, 60, 22, 0, 0, 0, 0, 2, tout, t)
3667 623 msg=
'DLSODES- ITASK = I1 and TOUT (=R1) behind TCUR - HU (= R2) ' 3668 CALL xerrwd (msg, 60, 23, 0, 1, itask, 0, 2, tout, tp)
3670 624 msg=
'DLSODES- ITASK = 4 or 5 and TCRIT (=R1) behind TCUR (=R2) ' 3671 CALL xerrwd (msg, 60, 24, 0, 0, 0, 0, 2, tcrit, tn)
3673 625 msg=
'DLSODES- ITASK = 4 or 5 and TCRIT (=R1) behind TOUT (=R2) ' 3674 CALL xerrwd (msg, 60, 25, 0, 0, 0, 0, 2, tcrit, tout)
3676 626 msg =
'DLSODES- At start of problem, too much accuracy ' 3677 CALL xerrwd (msg, 50, 26, 0, 0, 0, 0, 0, 0.0d0, 0.0d0)
3678 msg=
' requested for precision of machine.. See TOLSF (=R1) ' 3679 CALL xerrwd (msg, 60, 26, 0, 0, 0, 0, 1, tolsf, 0.0d0)
3682 627 msg =
'DLSODES- Trouble in DINTDY. ITASK = I1, TOUT = R1' 3683 CALL xerrwd (msg, 50, 27, 0, 1, itask, 0, 1, tout, 0.0d0)
3685 628 msg=
'DLSODES- RWORK length insufficient (for Subroutine DPREP). ' 3686 CALL xerrwd (msg, 60, 28, 0, 0, 0, 0, 0, 0.0d0, 0.0d0)
3687 msg=.ge.
' Length needed is LENRW (=I1), exceeds LRW (=I2)' 3688 CALL xerrwd (msg, 60, 28, 0, 2, lenrw, lrw, 0, 0.0d0, 0.0d0)
3690 629 msg=
'DLSODES- RWORK length insufficient (for Subroutine JGROUP). ' 3691 CALL xerrwd (msg, 60, 29, 0, 0, 0, 0, 0, 0.0d0, 0.0d0)
3692 msg=.ge.
' Length needed is LENRW (=I1), exceeds LRW (=I2)' 3693 CALL xerrwd (msg, 60, 29, 0, 2, lenrw, lrw, 0, 0.0d0, 0.0d0)
3695 630 msg=
'DLSODES- RWORK length insufficient (for Subroutine ODRV). ' 3696 CALL xerrwd (msg, 60, 30, 0, 0, 0, 0, 0, 0.0d0, 0.0d0)
3697 msg=.ge.
' Length needed is LENRW (=I1), exceeds LRW (=I2)' 3698 CALL xerrwd (msg, 60, 30, 0, 2, lenrw, lrw, 0, 0.0d0, 0.0d0)
3700 631 msg=
'DLSODES- Error from ODRV in Yale Sparse Matrix Package. ' 3701 CALL xerrwd (msg, 60, 31, 0, 0, 0, 0, 0, 0.0d0, 0.0d0)
3704 msg=
' At T (=R1), ODRV returned error flag = I1*NEQ + I2. ' 3705 CALL xerrwd (msg, 60, 31, 0, 2, imul, irem, 1, tn, 0.0d0)
3707 632 msg=
'DLSODES- RWORK length insufficient (for Subroutine CDRV). ' 3708 CALL xerrwd (msg, 60, 32, 0, 0, 0, 0, 0, 0.0d0, 0.0d0)
3709 msg=.ge.
' Length needed is LENRW (=I1), exceeds LRW (=I2)' 3710 CALL xerrwd (msg, 60, 32, 0, 2, lenrw, lrw, 0, 0.0d0, 0.0d0)
3712 633 msg=
'DLSODES- Error from CDRV in Yale Sparse Matrix Package. ' 3713 CALL xerrwd (msg, 60, 33, 0, 0, 0, 0, 0, 0.0d0, 0.0d0)
3716 msg=
' At T (=R1), CDRV returned error flag = I1*NEQ + I2. ' 3717 CALL xerrwd (msg, 60, 33, 0, 2, imul, irem, 1, tn, 0.0d0)
3718 IF (imul .EQ. 2)
THEN 3719 msg=
' Duplicate entry in sparsity structure descriptors. ' 3720 CALL xerrwd (msg, 60, 33, 0, 0, 0, 0, 0, 0.0d0, 0.0d0)
3722 IF (imul .EQ. 3 .OR. imul .EQ. 6)
THEN 3723 msg=
' Insufficient storage for NSFC (called by CDRV). ' 3724 CALL xerrwd (msg, 60, 33, 0, 0, 0, 0, 0, 0.0d0, 0.0d0)
3730 800 msg =
'DLSODES- Run aborted.. apparent infinite loop. ' 3731 CALL xerrwd (msg, 50, 303, 2, 0, 0, 0, 0, 0.0d0, 0.0d0)
3736 SUBROUTINE dlsoda (F, NEQ, Y, T, TOUT, ITOL, RTOL, ATOL, ITASK,
3737 1 ISTATE, IOPT, RWORK, LRW, IWORK, LIW, JAC, JT)
3739 INTEGER neq, itol, itask, istate, iopt, lrw, iwork, liw, jt
3740 DOUBLE PRECISION y, t, tout, rtol, atol, rwork
3741 dimension neq(*), y(*), rtol(*), atol(*), rwork(lrw), iwork(liw)
4721 EXTERNAL dprja, dsolsy
4722 DOUBLE PRECISION dumach, dmnorm
4723 INTEGER init, mxstep, mxhnil, nhnil, nslast, nyh, iowns,
4724 1 icf, ierpj, iersl, jcur, jstart, kflag, l,
4725 2 lyh, lewt, lacor, lsavf, lwm, liwm, meth, miter,
4726 3 maxord, maxcor, msbp, mxncf, n, nq, nst, nfe, nje, nqu
4727 INTEGER insufr, insufi, ixpr, iowns2, jtyp, mused, mxordn, mxords
4728 INTEGER i, i1, i2, iflag, imxer, kgo, lf0,
4729 1 leniw, lenrw, lenwm, ml, mord, mu, mxhnl0, mxstp0
4730 INTEGER len1, len1c, len1n, len1s, len2, leniwc, lenrwc
4731 DOUBLE PRECISION rowns,
4732 1 ccmax, el0, h, hmin, hmxi, hu, rc, tn, uround
4733 DOUBLE PRECISION tsw, rowns2, pdnorm
4734 DOUBLE PRECISION atoli, ayi, big, ewti, h0, hmax, hmx, rh, rtoli,
4735 1 tcrit, tdist, tnext, tol, tolsf, tp,
SIZE, sum, w0
4739 SAVE mord, mxstp0, mxhnl0
4751 COMMON /dls001/ rowns(209),
4752 1 ccmax, el0, h, hmin, hmxi, hu, rc, tn, uround,
4753 2 init, mxstep, mxhnil, nhnil, nslast, nyh, iowns(6),
4754 3 icf, ierpj, iersl, jcur, jstart, kflag, l,
4755 4 lyh, lewt, lacor, lsavf, lwm, liwm, meth, miter,
4756 5 maxord, maxcor, msbp, mxncf, n, nq, nst, nfe, nje, nqu
4758 COMMON /dlsa01/ tsw, rowns2(20), pdnorm,
4759 1 insufr, insufi, ixpr, iowns2(2), jtyp, mused, mxordn, mxords
4761 DATA mord(1),mord(2)/12,5/, mxstp0/500/, mxhnl0/10/
4770 IF (istate .LT. 1 .OR. istate .GT. 3)
GO TO 601
4771 IF (itask .LT. 1 .OR. itask .GT. 5)
GO TO 602
4772 IF (istate .EQ. 1)
GO TO 10
4773 IF (init .EQ. 0)
GO TO 603
4774 IF (istate .EQ. 2)
GO TO 200
4777 IF (tout .EQ. t)
RETURN 4787 20
IF (neq(1) .LE. 0)
GO TO 604
4788 IF (istate .EQ. 1)
GO TO 25
4789 IF (neq(1) .GT. n)
GO TO 605
4791 IF (itol .LT. 1 .OR. itol .GT. 4)
GO TO 606
4792 IF (iopt .LT. 0 .OR. iopt .GT. 1)
GO TO 607
4793 IF (jt .EQ. 3 .OR. jt .LT. 1 .OR. jt .GT. 5)
GO TO 608
4795 IF (jt .LE. 2)
GO TO 30
4798 IF (ml .LT. 0 .OR. ml .GE. n)
GO TO 609
4799 IF (mu .LT. 0 .OR. mu .GE. n)
GO TO 610
4802 IF (iopt .EQ. 1)
GO TO 40
4808 IF (istate .NE. 1)
GO TO 60
4814 IF (ixpr .LT. 0 .OR. ixpr .GT. 1)
GO TO 611
4816 IF (mxstep .LT. 0)
GO TO 612
4817 IF (mxstep .EQ. 0) mxstep = mxstp0
4819 IF (mxhnil .LT. 0)
GO TO 613
4820 IF (mxhnil .EQ. 0) mxhnil = mxhnl0
4821 IF (istate .NE. 1)
GO TO 50
4824 IF (mxordn .LT. 0)
GO TO 628
4825 IF (mxordn .EQ. 0) mxordn = 100
4826 mxordn = min(mxordn,mord(1))
4828 IF (mxords .LT. 0)
GO TO 629
4829 IF (mxords .EQ. 0) mxords = 100
4830 mxords = min(mxords,mord(2))
4831 IF ((tout - t)*h0 .LT. 0.0d0)
GO TO 614
4833 IF (hmax .LT. 0.0d0)
GO TO 615
4835 IF (hmax .GT. 0.0d0) hmxi = 1.0d0/hmax
4837 IF (hmin .LT. 0.0d0)
GO TO 616
4851 60
IF (istate .EQ. 1) meth = 1
4852 IF (istate .EQ. 1) nyh = n
4854 len1n = 20 + (mxordn + 1)*nyh
4855 len1s = 20 + (mxords + 1)*nyh
4857 IF (jt .LE. 2) lenwm = n*n + 2
4858 IF (jt .GE. 4) lenwm = (2*ml + mu + 1)*n + 2
4859 len1s = len1s + lenwm
4861 IF (meth .EQ. 2) len1c = len1s
4862 len1 = max(len1n,len1s)
4865 lenrwc = len1c + len2
4870 IF (meth .EQ. 2) leniwc = leniw
4872 IF (istate .EQ. 1 .AND. lrw .LT. lenrwc)
GO TO 617
4873 IF (istate .EQ. 1 .AND. liw .LT. leniwc)
GO TO 618
4874 IF (istate .EQ. 3 .AND. lrw .LT. lenrwc)
GO TO 550
4875 IF (istate .EQ. 3 .AND. liw .LT. leniwc)
GO TO 555
4878 IF (lrw .GE. lenrw)
GO TO 65
4881 msg=
'DLSODA- Warning.. RWORK length is sufficient for now, but ' 4882 CALL xerrwd (msg, 60, 103, 0, 0, 0, 0, 0, 0.0d0, 0.0d0)
4883 msg=
' may not be later. Integration will proceed anyway. ' 4884 CALL xerrwd (msg, 60, 103, 0, 0, 0, 0, 0, 0.0d0, 0.0d0)
4885 msg =
' Length needed is LENRW = I1, while LRW = I2.' 4886 CALL xerrwd (msg, 50, 103, 0, 2, lenrw, lrw, 0, 0.0d0, 0.0d0)
4890 IF (liw .GE. leniw)
GO TO 70
4892 msg=
'DLSODA- Warning.. IWORK length is sufficient for now, but ' 4893 CALL xerrwd (msg, 60, 104, 0, 0, 0, 0, 0, 0.0d0, 0.0d0)
4894 msg=
' may not be later. Integration will proceed anyway. ' 4895 CALL xerrwd (msg, 60, 104, 0, 0, 0, 0, 0, 0.0d0, 0.0d0)
4896 msg =
' Length needed is LENIW = I1, while LIW = I2.' 4897 CALL xerrwd (msg, 50, 104, 0, 2, leniw, liw, 0, 0.0d0, 0.0d0)
4903 IF (itol .GE. 3) rtoli = rtol(i)
4904 IF (itol .EQ. 2 .OR. itol .EQ. 4) atoli = atol(i)
4905 IF (rtoli .LT. 0.0d0)
GO TO 619
4906 IF (atoli .LT. 0.0d0)
GO TO 620
4908 IF (istate .EQ. 1)
GO TO 100
4911 IF (n .EQ. nyh)
GO TO 200
4914 i2 = lyh + (maxord + 1)*nyh - 1
4915 IF (i1 .GT. i2)
GO TO 200
4926 100 uround = dumach()
4930 IF (itask .NE. 4 .AND. itask .NE. 5)
GO TO 110
4932 IF ((tcrit - tout)*(tout - t) .LT. 0.0d0)
GO TO 625
4933 IF (h0 .NE. 0.0d0 .AND. (t + h0 - tcrit)*h0 .GT. 0.0d0)
4950 CALL f (neq, t, y, rwork(lf0))
4954 115 rwork(i+lyh-1) = y(i)
4958 CALL dewset (n, itol, rtol, atol, rwork(lyh), rwork(lewt))
4960 IF (rwork(i+lewt-1) .LE. 0.0d0)
GO TO 621
4961 120 rwork(i+lewt-1) = 1.0d0/rwork(i+lewt-1)
4981 IF (h0 .NE. 0.0d0)
GO TO 180
4982 tdist = abs(tout - t)
4983 w0 = max(abs(t),abs(tout))
4984 IF (tdist .LT. 2.0d0*uround*w0)
GO TO 622
4986 IF (itol .LE. 2)
GO TO 140
4988 130 tol = max(tol,rtol(i))
4989 140
IF (tol .GT. 0.0d0)
GO TO 160
4992 IF (itol .EQ. 2 .OR. itol .EQ. 4) atoli = atol(i)
4994 IF (ayi .NE. 0.0d0) tol = max(tol,atoli/ayi)
4996 160 tol = max(tol,100.0d0*uround)
4997 tol = min(tol,0.001d0)
4998 sum = dmnorm(n, rwork(lf0), rwork(lewt))
4999 sum = 1.0d0/(tol*w0*w0) + tol*sum**2
5000 h0 = 1.0d0/sqrt(sum)
5002 h0 = sign(h0,tout-t)
5004 180 rh = abs(h0)*hmxi
5005 IF (rh .GT. 1.0d0) h0 = h0/rh
5009 190 rwork(i+lf0-1) = h0*rwork(i+lf0-1)
5017 GO TO (210, 250, 220, 230, 240), itask
5018 210
IF ((tn - tout)*h .LT. 0.0d0)
GO TO 250
5019 CALL dintdy (tout, 0, rwork(lyh), nyh, y, iflag)
5020 IF (iflag .NE. 0)
GO TO 627
5023 220 tp = tn - hu*(1.0d0 + 100.0d0*uround)
5024 IF ((tp - tout)*h .GT. 0.0d0)
GO TO 623
5025 IF ((tn - tout)*h .LT. 0.0d0)
GO TO 250
5028 230 tcrit = rwork(1)
5029 IF ((tn - tcrit)*h .GT. 0.0d0)
GO TO 624
5030 IF ((tcrit - tout)*h .LT. 0.0d0)
GO TO 625
5031 IF ((tn - tout)*h .LT. 0.0d0)
GO TO 245
5032 CALL dintdy (tout, 0, rwork(lyh), nyh, y, iflag)
5033 IF (iflag .NE. 0)
GO TO 627
5036 240 tcrit = rwork(1)
5037 IF ((tn - tcrit)*h .GT. 0.0d0)
GO TO 624
5038 245 hmx = abs(tn) + abs(h)
5039 ihit = abs(tn - tcrit) .LE. 100.0d0*uround*hmx
5042 tnext = tn + h*(1.0d0 + 4.0d0*uround)
5043 IF ((tnext - tcrit)*h .LE. 0.0d0)
GO TO 250
5044 h = (tcrit - tn)*(1.0d0 - 4.0d0*uround)
5045 IF (istate .EQ. 2 .AND. jstart .GE. 0) jstart = -2
5058 IF (meth .EQ. mused)
GO TO 255
5059 IF (insufr .EQ. 1)
GO TO 550
5060 IF (insufi .EQ. 1)
GO TO 555
5061 255
IF ((nst-nslast) .GE. mxstep)
GO TO 500
5062 CALL dewset (n, itol, rtol, atol, rwork(lyh), rwork(lewt))
5064 IF (rwork(i+lewt-1) .LE. 0.0d0)
GO TO 510
5065 260 rwork(i+lewt-1) = 1.0d0/rwork(i+lewt-1)
5066 270 tolsf = uround*dmnorm(n, rwork(lyh), rwork(lewt))
5067 IF (tolsf .LE. 1.0d0)
GO TO 280
5069 IF (nst .EQ. 0)
GO TO 626
5071 280
IF ((tn + h) .NE. tn)
GO TO 290
5073 IF (nhnil .GT. mxhnil)
GO TO 290
5074 msg =
'DLSODA- Warning..Internal T (=R1) and H (=R2) are' 5075 CALL xerrwd (msg, 50, 101, 0, 0, 0, 0, 0, 0.0d0, 0.0d0)
5076 msg=
' such that in the machine, T + H = T on the next step ' 5077 CALL xerrwd (msg, 60, 101, 0, 0, 0, 0, 0, 0.0d0, 0.0d0)
5078 msg =
' (H = step size). Solver will continue anyway.' 5079 CALL xerrwd (msg, 50, 101, 0, 0, 0, 0, 2, tn, h)
5080 IF (nhnil .LT. mxhnil)
GO TO 290
5081 msg =
'DLSODA- Above warning has been issued I1 times. ' 5082 CALL xerrwd (msg, 50, 102, 0, 0, 0, 0, 0, 0.0d0, 0.0d0)
5083 msg =
' It will not be issued again for this problem.' 5084 CALL xerrwd (msg, 50, 102, 0, 1, mxhnil, 0, 0, 0.0d0, 0.0d0)
5089 CALL dstoda (neq, y, rwork(lyh), nyh, rwork(lyh), rwork(lewt),
5090 1 rwork(lsavf), rwork(lacor), rwork(lwm), iwork(liwm),
5091 2 f, jac, dprja, dsolsy)
5093 GO TO (300, 530, 540), kgo
5104 IF (meth .EQ. mused)
GO TO 310
5107 IF (meth .EQ. 2) maxord = mxords
5108 IF (meth .EQ. 2) rwork(lwm) = sqrt(uround)
5109 insufr = min(insufr,1)
5110 insufi = min(insufi,1)
5112 IF (ixpr .EQ. 0)
GO TO 310
5113 IF (meth .EQ. 2)
THEN 5114 msg=
'DLSODA- A switch to the BDF (stiff) method has occurred ' 5115 CALL xerrwd (msg, 60, 105, 0, 0, 0, 0, 0, 0.0d0, 0.0d0)
5117 IF (meth .EQ. 1)
THEN 5118 msg=
'DLSODA- A switch to the Adams (nonstiff) method has occurred' 5119 CALL xerrwd (msg, 60, 106, 0, 0, 0, 0, 0, 0.0d0, 0.0d0)
5121 msg=
' at T = R1, tentative step size H = R2, step NST = I1 ' 5122 CALL xerrwd (msg, 60, 107, 0, 1, nst, 0, 2, tn, h)
5123 310
GO TO (320, 400, 330, 340, 350), itask
5125 320
IF ((tn - tout)*h .LT. 0.0d0)
GO TO 250
5126 CALL dintdy (tout, 0, rwork(lyh), nyh, y, iflag)
5130 330
IF ((tn - tout)*h .GE. 0.0d0)
GO TO 400
5133 340
IF ((tn - tout)*h .LT. 0.0d0)
GO TO 345
5134 CALL dintdy (tout, 0, rwork(lyh), nyh, y, iflag)
5137 345 hmx = abs(tn) + abs(h)
5138 ihit = abs(tn - tcrit) .LE. 100.0d0*uround*hmx
5140 tnext = tn + h*(1.0d0 + 4.0d0*uround)
5141 IF ((tnext - tcrit)*h .LE. 0.0d0)
GO TO 250
5142 h = (tcrit - tn)*(1.0d0 - 4.0d0*uround)
5143 IF (jstart .GE. 0) jstart = -2
5146 350 hmx = abs(tn) + abs(h)
5147 ihit = abs(tn - tcrit) .LE. 100.0d0*uround*hmx
5156 410 y(i) = rwork(i+lyh-1)
5158 IF (itask .NE. 4 .AND. itask .NE. 5)
GO TO 420
5182 500 msg =
'DLSODA- At current T (=R1), MXSTEP (=I1) steps ' 5183 CALL xerrwd (msg, 50, 201, 0, 0, 0, 0, 0, 0.0d0, 0.0d0)
5184 msg =
' taken on this call before reaching TOUT ' 5185 CALL xerrwd (msg, 50, 201, 0, 1, mxstep, 0, 1, tn, 0.0d0)
5189 510 ewti = rwork(lewt+i-1)
5190 msg = .le.
'DLSODA- At T (=R1), EWT(I1) has become R2 0.' 5191 CALL xerrwd (msg, 50, 202, 0, 1, i, 0, 2, tn, ewti)
5195 520 msg =
'DLSODA- At T (=R1), too much accuracy requested ' 5196 CALL xerrwd (msg, 50, 203, 0, 0, 0, 0, 0, 0.0d0, 0.0d0)
5197 msg =
' for precision of machine.. See TOLSF (=R2) ' 5198 CALL xerrwd (msg, 50, 203, 0, 0, 0, 0, 2, tn, tolsf)
5203 530 msg =
'DLSODA- At T(=R1) and step size H(=R2), the error' 5204 CALL xerrwd (msg, 50, 204, 0, 0, 0, 0, 0, 0.0d0, 0.0d0)
5205 msg =
' test failed repeatedly or with ABS(H) = HMIN' 5206 CALL xerrwd (msg, 50, 204, 0, 0, 0, 0, 2, tn, h)
5210 540 msg =
'DLSODA- At T (=R1) and step size H (=R2), the ' 5211 CALL xerrwd (msg, 50, 205, 0, 0, 0, 0, 0, 0.0d0, 0.0d0)
5212 msg =
' corrector convergence failed repeatedly ' 5213 CALL xerrwd (msg, 50, 205, 0, 0, 0, 0, 0, 0.0d0, 0.0d0)
5214 msg =
' or with ABS(H) = HMIN ' 5215 CALL xerrwd (msg, 30, 205, 0, 0, 0, 0, 2, tn, h)
5219 550 msg =
'DLSODA- At current T(=R1), RWORK length too small' 5220 CALL xerrwd (msg, 50, 206, 0, 0, 0, 0, 0, 0.0d0, 0.0d0)
5221 msg=
' to proceed. The integration was otherwise successful.' 5222 CALL xerrwd (msg, 60, 206, 0, 0, 0, 0, 1, tn, 0.0d0)
5226 555 msg =
'DLSODA- At current T(=R1), IWORK length too small' 5227 CALL xerrwd (msg, 50, 207, 0, 0, 0, 0, 0, 0.0d0, 0.0d0)
5228 msg=
' to proceed. The integration was otherwise successful.' 5229 CALL xerrwd (msg, 60, 207, 0, 0, 0, 0, 1, tn, 0.0d0)
5236 SIZE = abs(rwork(i+lacor-1)*rwork(i+lewt-1))
5237 IF (big .GE. size)
GO TO 570
5244 590 y(i) = rwork(i+lyh-1)
5265 601 msg =
'DLSODA- ISTATE (=I1) illegal.' 5266 CALL xerrwd (msg, 30, 1, 0, 1, istate, 0, 0, 0.0d0, 0.0d0)
5267 IF (istate .LT. 0)
GO TO 800
5269 602 msg =
'DLSODA- ITASK (=I1) illegal. ' 5270 CALL xerrwd (msg, 30, 2, 0, 1, itask, 0, 0, 0.0d0, 0.0d0)
5272 603 msg = .gt.
'DLSODA- ISTATE 1 but DLSODA not initialized.' 5273 CALL xerrwd (msg, 50, 3, 0, 0, 0, 0, 0, 0.0d0, 0.0d0)
5275 604 msg = .lt.
'DLSODA- NEQ (=I1) 1 ' 5276 CALL xerrwd (msg, 30, 4, 0, 1, neq(1), 0, 0, 0.0d0, 0.0d0)
5278 605 msg =
'DLSODA- ISTATE = 3 and NEQ increased (I1 to I2). ' 5279 CALL xerrwd (msg, 50, 5, 0, 2, n, neq(1), 0, 0.0d0, 0.0d0)
5281 606 msg =
'DLSODA- ITOL (=I1) illegal. ' 5282 CALL xerrwd (msg, 30, 6, 0, 1, itol, 0, 0, 0.0d0, 0.0d0)
5284 607 msg =
'DLSODA- IOPT (=I1) illegal. ' 5285 CALL xerrwd (msg, 30, 7, 0, 1, iopt, 0, 0, 0.0d0, 0.0d0)
5287 608 msg =
'DLSODA- JT (=I1) illegal. ' 5288 CALL xerrwd (msg, 30, 8, 0, 1, jt, 0, 0, 0.0d0, 0.0d0)
5290 609 msg = .lt..ge.
'DLSODA- ML (=I1) illegal: 0 or NEQ (=I2) ' 5291 CALL xerrwd (msg, 50, 9, 0, 2, ml, neq(1), 0, 0.0d0, 0.0d0)
5293 610 msg = .lt..ge.
'DLSODA- MU (=I1) illegal: 0 or NEQ (=I2) ' 5294 CALL xerrwd (msg, 50, 10, 0, 2, mu, neq(1), 0, 0.0d0, 0.0d0)
5296 611 msg =
'DLSODA- IXPR (=I1) illegal. ' 5297 CALL xerrwd (msg, 30, 11, 0, 1, ixpr, 0, 0, 0.0d0, 0.0d0)
5299 612 msg = .lt.
'DLSODA- MXSTEP (=I1) 0 ' 5300 CALL xerrwd (msg, 30, 12, 0, 1, mxstep, 0, 0, 0.0d0, 0.0d0)
5302 613 msg = .lt.
'DLSODA- MXHNIL (=I1) 0 ' 5303 CALL xerrwd (msg, 30, 13, 0, 1, mxhnil, 0, 0, 0.0d0, 0.0d0)
5305 614 msg =
'DLSODA- TOUT (=R1) behind T (=R2) ' 5306 CALL xerrwd (msg, 40, 14, 0, 0, 0, 0, 2, tout, t)
5307 msg =
' Integration direction is given by H0 (=R1) ' 5308 CALL xerrwd (msg, 50, 14, 0, 0, 0, 0, 1, h0, 0.0d0)
5310 615 msg = .lt.
'DLSODA- HMAX (=R1) 0.0 ' 5311 CALL xerrwd (msg, 30, 15, 0, 0, 0, 0, 1, hmax, 0.0d0)
5313 616 msg = .lt.
'DLSODA- HMIN (=R1) 0.0 ' 5314 CALL xerrwd (msg, 30, 16, 0, 0, 0, 0, 1, hmin, 0.0d0)
5316 617 msg=
'DLSODA- RWORK length needed, LENRW (=I1), exceeds LRW (=I2)' 5317 CALL xerrwd (msg, 60, 17, 0, 2, lenrw, lrw, 0, 0.0d0, 0.0d0)
5319 618 msg=
'DLSODA- IWORK length needed, LENIW (=I1), exceeds LIW (=I2)' 5320 CALL xerrwd (msg, 60, 18, 0, 2, leniw, liw, 0, 0.0d0, 0.0d0)
5322 619 msg = .lt.
'DLSODA- RTOL(I1) is R1 0.0 ' 5323 CALL xerrwd (msg, 40, 19, 0, 1, i, 0, 1, rtoli, 0.0d0)
5325 620 msg = .lt.
'DLSODA- ATOL(I1) is R1 0.0 ' 5326 CALL xerrwd (msg, 40, 20, 0, 1, i, 0, 1, atoli, 0.0d0)
5328 621 ewti = rwork(lewt+i-1)
5329 msg = .le.
'DLSODA- EWT(I1) is R1 0.0 ' 5330 CALL xerrwd (msg, 40, 21, 0, 1, i, 0, 1, ewti, 0.0d0)
5332 622 msg=
'DLSODA- TOUT(=R1) too close to T(=R2) to start integration.' 5333 CALL xerrwd (msg, 60, 22, 0, 0, 0, 0, 2, tout, t)
5335 623 msg=
'DLSODA- ITASK = I1 and TOUT (=R1) behind TCUR - HU (= R2) ' 5336 CALL xerrwd (msg, 60, 23, 0, 1, itask, 0, 2, tout, tp)
5338 624 msg=
'DLSODA- ITASK = 4 or 5 and TCRIT (=R1) behind TCUR (=R2) ' 5339 CALL xerrwd (msg, 60, 24, 0, 0, 0, 0, 2, tcrit, tn)
5341 625 msg=
'DLSODA- ITASK = 4 or 5 and TCRIT (=R1) behind TOUT (=R2) ' 5342 CALL xerrwd (msg, 60, 25, 0, 0, 0, 0, 2, tcrit, tout)
5344 626 msg =
'DLSODA- At start of problem, too much accuracy ' 5345 CALL xerrwd (msg, 50, 26, 0, 0, 0, 0, 0, 0.0d0, 0.0d0)
5346 msg=
' requested for precision of machine.. See TOLSF (=R1) ' 5347 CALL xerrwd (msg, 60, 26, 0, 0, 0, 0, 1, tolsf, 0.0d0)
5350 627 msg =
'DLSODA- Trouble in DINTDY. ITASK = I1, TOUT = R1' 5351 CALL xerrwd (msg, 50, 27, 0, 1, itask, 0, 1, tout, 0.0d0)
5353 628 msg = .lt.
'DLSODA- MXORDN (=I1) 0 ' 5354 CALL xerrwd (msg, 30, 28, 0, 1, mxordn, 0, 0, 0.0d0, 0.0d0)
5356 629 msg = .lt.
'DLSODA- MXORDS (=I1) 0 ' 5357 CALL xerrwd (msg, 30, 29, 0, 1, mxords, 0, 0, 0.0d0, 0.0d0)
5362 800 msg =
'DLSODA- Run aborted.. apparent infinite loop. ' 5363 CALL xerrwd (msg, 50, 303, 2, 0, 0, 0, 0, 0.0d0, 0.0d0)
5368 SUBROUTINE dlsodar (F, NEQ, Y, T, TOUT, ITOL, RTOL, ATOL, ITASK,
5369 1 ISTATE, IOPT, RWORK, LRW, IWORK, LIW, JAC, JT,
5372 INTEGER neq, itol, itask, istate, iopt, lrw, iwork, liw, jt,
5374 DOUBLE PRECISION y, t, tout, rtol, atol, rwork
5375 dimension neq(*), y(*), rtol(*), atol(*), rwork(lrw), iwork(liw),
6475 EXTERNAL dprja, dsolsy
6476 DOUBLE PRECISION dumach, dmnorm
6477 INTEGER init, mxstep, mxhnil, nhnil, nslast, nyh, iowns,
6478 1 icf, ierpj, iersl, jcur, jstart, kflag, l,
6479 2 lyh, lewt, lacor, lsavf, lwm, liwm, meth, miter,
6480 3 maxord, maxcor, msbp, mxncf, n, nq, nst, nfe, nje, nqu
6481 INTEGER insufr, insufi, ixpr, iowns2, jtyp, mused, mxordn, mxords
6482 INTEGER lg0, lg1, lgx, iownr3, irfnd, itaskc, ngc, nge
6483 INTEGER i, i1, i2, iflag, imxer, kgo, leniw,
6484 1 lenrw, lenwm, lf0, ml, mord, mu, mxhnl0, mxstp0
6485 INTEGER len1, len1c, len1n, len1s, len2, leniwc, lenrwc
6486 INTEGER irfp, irt, lenyh, lyhnew
6487 DOUBLE PRECISION rowns,
6488 1 ccmax, el0, h, hmin, hmxi, hu, rc, tn, uround
6489 DOUBLE PRECISION tsw, rowns2, pdnorm
6490 DOUBLE PRECISION rownr3, t0, tlast, toutc
6491 DOUBLE PRECISION atoli, ayi, big, ewti, h0, hmax, hmx, rh, rtoli,
6492 1 tcrit, tdist, tnext, tol, tolsf, tp,
SIZE, sum, w0
6496 SAVE mord, mxstp0, mxhnl0
6509 COMMON /dls001/ rowns(209),
6510 1 ccmax, el0, h, hmin, hmxi, hu, rc, tn, uround,
6511 2 init, mxstep, mxhnil, nhnil, nslast, nyh, iowns(6),
6512 3 icf, ierpj, iersl, jcur, jstart, kflag, l,
6513 4 lyh, lewt, lacor, lsavf, lwm, liwm, meth, miter,
6514 5 maxord, maxcor, msbp, mxncf, n, nq, nst, nfe, nje, nqu
6516 COMMON /dlsa01/ tsw, rowns2(20), pdnorm,
6517 1 insufr, insufi, ixpr, iowns2(2), jtyp, mused, mxordn, mxords
6519 COMMON /dlsr01/ rownr3(2), t0, tlast, toutc,
6520 1 lg0, lg1, lgx, iownr3(2), irfnd, itaskc, ngc, nge
6522 DATA mord(1),mord(2)/12,5/, mxstp0/500/, mxhnl0/10/
6531 IF (istate .LT. 1 .OR. istate .GT. 3)
GO TO 601
6532 IF (itask .LT. 1 .OR. itask .GT. 5)
GO TO 602
6534 IF (istate .EQ. 1)
GO TO 10
6535 IF (init .EQ. 0)
GO TO 603
6536 IF (istate .EQ. 2)
GO TO 200
6539 IF (tout .EQ. t)
RETURN 6549 20
IF (neq(1) .LE. 0)
GO TO 604
6550 IF (istate .EQ. 1)
GO TO 25
6551 IF (neq(1) .GT. n)
GO TO 605
6553 IF (itol .LT. 1 .OR. itol .GT. 4)
GO TO 606
6554 IF (iopt .LT. 0 .OR. iopt .GT. 1)
GO TO 607
6555 IF (jt .EQ. 3 .OR. jt .LT. 1 .OR. jt .GT. 5)
GO TO 608
6557 IF (jt .LE. 2)
GO TO 30
6560 IF (ml .LT. 0 .OR. ml .GE. n)
GO TO 609
6561 IF (mu .LT. 0 .OR. mu .GE. n)
GO TO 610
6563 IF (ng .LT. 0)
GO TO 630
6564 IF (istate .EQ. 1)
GO TO 35
6565 IF (irfnd .EQ. 0 .AND. ng .NE. ngc)
GO TO 631
6568 IF (iopt .EQ. 1)
GO TO 40
6574 IF (istate .NE. 1)
GO TO 60
6580 IF (ixpr .LT. 0 .OR. ixpr .GT. 1)
GO TO 611
6582 IF (mxstep .LT. 0)
GO TO 612
6583 IF (mxstep .EQ. 0) mxstep = mxstp0
6585 IF (mxhnil .LT. 0)
GO TO 613
6586 IF (mxhnil .EQ. 0) mxhnil = mxhnl0
6587 IF (istate .NE. 1)
GO TO 50
6590 IF (mxordn .LT. 0)
GO TO 628
6591 IF (mxordn .EQ. 0) mxordn = 100
6592 mxordn = min(mxordn,mord(1))
6594 IF (mxords .LT. 0)
GO TO 629
6595 IF (mxords .EQ. 0) mxords = 100
6596 mxords = min(mxords,mord(2))
6597 IF ((tout - t)*h0 .LT. 0.0d0)
GO TO 614
6599 IF (hmax .LT. 0.0d0)
GO TO 615
6601 IF (hmax .GT. 0.0d0) hmxi = 1.0d0/hmax
6603 IF (hmin .LT. 0.0d0)
GO TO 616
6618 60
IF (istate .EQ. 1) meth = 1
6619 IF (istate .EQ. 1) nyh = n
6624 IF (istate .EQ. 1) lyh = lyhnew
6625 IF (lyhnew .EQ. lyh)
GO TO 62
6628 IF (lrw .LT. lyhnew-1+lenyh)
GO TO 62
6630 IF (lyhnew .GT. lyh) i1 = -1
6631 CALL dcopy (lenyh, rwork(lyh), i1, rwork(lyhnew), i1)
6634 len1n = lyhnew - 1 + (mxordn + 1)*nyh
6635 len1s = lyhnew - 1 + (mxords + 1)*nyh
6637 IF (jt .LE. 2) lenwm = n*n + 2
6638 IF (jt .GE. 4) lenwm = (2*ml + mu + 1)*n + 2
6639 len1s = len1s + lenwm
6641 IF (meth .EQ. 2) len1c = len1s
6642 len1 = max(len1n,len1s)
6645 lenrwc = len1c + len2
6650 IF (meth .EQ. 2) leniwc = leniw
6652 IF (istate .EQ. 1 .AND. lrw .LT. lenrwc)
GO TO 617
6653 IF (istate .EQ. 1 .AND. liw .LT. leniwc)
GO TO 618
6654 IF (istate .EQ. 3 .AND. lrw .LT. lenrwc)
GO TO 550
6655 IF (istate .EQ. 3 .AND. liw .LT. leniwc)
GO TO 555
6658 IF (lrw .GE. lenrw)
GO TO 65
6661 msg=
'DLSODAR- Warning.. RWORK length is sufficient for now, but ' 6662 CALL xerrwd (msg, 60, 103, 0, 0, 0, 0, 0, 0.0d0, 0.0d0)
6663 msg=
' may not be later. Integration will proceed anyway. ' 6664 CALL xerrwd (msg, 60, 103, 0, 0, 0, 0, 0, 0.0d0, 0.0d0)
6665 msg =
' Length needed is LENRW = I1, while LRW = I2.' 6666 CALL xerrwd (msg, 50, 103, 0, 2, lenrw, lrw, 0, 0.0d0, 0.0d0)
6670 IF (liw .GE. leniw)
GO TO 70
6672 msg=
'DLSODAR- Warning.. IWORK length is sufficient for now, but ' 6673 CALL xerrwd (msg, 60, 104, 0, 0, 0, 0, 0, 0.0d0, 0.0d0)
6674 msg=
' may not be later. Integration will proceed anyway. ' 6675 CALL xerrwd (msg, 60, 104, 0, 0, 0, 0, 0, 0.0d0, 0.0d0)
6676 msg =
' Length needed is LENIW = I1, while LIW = I2.' 6677 CALL xerrwd (msg, 50, 104, 0, 2, leniw, liw, 0, 0.0d0, 0.0d0)
6683 IF (itol .GE. 3) rtoli = rtol(i)
6684 IF (itol .EQ. 2 .OR. itol .EQ. 4) atoli = atol(i)
6685 IF (rtoli .LT. 0.0d0)
GO TO 619
6686 IF (atoli .LT. 0.0d0)
GO TO 620
6688 IF (istate .EQ. 1)
GO TO 100
6691 IF (n .EQ. nyh)
GO TO 200
6694 i2 = lyh + (maxord + 1)*nyh - 1
6695 IF (i1 .GT. i2)
GO TO 200
6706 100 uround = dumach()
6710 IF (itask .NE. 4 .AND. itask .NE. 5)
GO TO 110
6712 IF ((tcrit - tout)*(tout - t) .LT. 0.0d0)
GO TO 625
6713 IF (h0 .NE. 0.0d0 .AND. (t + h0 - tcrit)*h0 .GT. 0.0d0)
6730 CALL f (neq, t, y, rwork(lf0))
6734 115 rwork(i+lyh-1) = y(i)
6738 CALL dewset (n, itol, rtol, atol, rwork(lyh), rwork(lewt))
6740 IF (rwork(i+lewt-1) .LE. 0.0d0)
GO TO 621
6741 120 rwork(i+lewt-1) = 1.0d0/rwork(i+lewt-1)
6761 IF (h0 .NE. 0.0d0)
GO TO 180
6762 tdist = abs(tout - t)
6763 w0 = max(abs(t),abs(tout))
6764 IF (tdist .LT. 2.0d0*uround*w0)
GO TO 622
6766 IF (itol .LE. 2)
GO TO 140
6768 130 tol = max(tol,rtol(i))
6769 140
IF (tol .GT. 0.0d0)
GO TO 160
6772 IF (itol .EQ. 2 .OR. itol .EQ. 4) atoli = atol(i)
6774 IF (ayi .NE. 0.0d0) tol = max(tol,atoli/ayi)
6776 160 tol = max(tol,100.0d0*uround)
6777 tol = min(tol,0.001d0)
6778 sum = dmnorm(n, rwork(lf0), rwork(lewt))
6779 sum = 1.0d0/(tol*w0*w0) + tol*sum**2
6780 h0 = 1.0d0/sqrt(sum)
6782 h0 = sign(h0,tout-t)
6784 180 rh = abs(h0)*hmxi
6785 IF (rh .GT. 1.0d0) h0 = h0/rh
6789 190 rwork(i+lf0-1) = h0*rwork(i+lf0-1)
6794 IF (ngc .EQ. 0)
GO TO 270
6795 CALL drchek (1, g, neq, y, rwork(lyh), nyh,
6796 1 rwork(lg0), rwork(lg1), rwork(lgx), jroot, irt)
6797 IF (irt .EQ. 0)
GO TO 270
6811 IF (ngc .EQ. 0)
GO TO 205
6812 IF (itask .EQ. 1 .OR. itask .EQ. 4) toutc = tout
6813 CALL drchek (2, g, neq, y, rwork(lyh), nyh,
6814 1 rwork(lg0), rwork(lg1), rwork(lgx), jroot, irt)
6815 IF (irt .NE. 1)
GO TO 205
6822 IF (irfp .EQ. 1 .AND. tlast .NE. tn .AND. itask .EQ. 2)
GO TO 400
6824 GO TO (210, 250, 220, 230, 240), itask
6825 210
IF ((tn - tout)*h .LT. 0.0d0)
GO TO 250
6826 CALL dintdy (tout, 0, rwork(lyh), nyh, y, iflag)
6827 IF (iflag .NE. 0)
GO TO 627
6830 220 tp = tn - hu*(1.0d0 + 100.0d0*uround)
6831 IF ((tp - tout)*h .GT. 0.0d0)
GO TO 623
6832 IF ((tn - tout)*h .LT. 0.0d0)
GO TO 250
6835 230 tcrit = rwork(1)
6836 IF ((tn - tcrit)*h .GT. 0.0d0)
GO TO 624
6837 IF ((tcrit - tout)*h .LT. 0.0d0)
GO TO 625
6838 IF ((tn - tout)*h .LT. 0.0d0)
GO TO 245
6839 CALL dintdy (tout, 0, rwork(lyh), nyh, y, iflag)
6840 IF (iflag .NE. 0)
GO TO 627
6843 240 tcrit = rwork(1)
6844 IF ((tn - tcrit)*h .GT. 0.0d0)
GO TO 624
6845 245 hmx = abs(tn) + abs(h)
6846 ihit = abs(tn - tcrit) .LE. 100.0d0*uround*hmx
6848 IF (irfp .EQ. 1 .AND. tlast .NE. tn .AND. itask .EQ. 5)
GO TO 400
6850 tnext = tn + h*(1.0d0 + 4.0d0*uround)
6851 IF ((tnext - tcrit)*h .LE. 0.0d0)
GO TO 250
6852 h = (tcrit - tn)*(1.0d0 - 4.0d0*uround)
6853 IF (istate .EQ. 2 .AND. jstart .GE. 0) jstart = -2
6866 IF (meth .EQ. mused)
GO TO 255
6867 IF (insufr .EQ. 1)
GO TO 550
6868 IF (insufi .EQ. 1)
GO TO 555
6869 255
IF ((nst-nslast) .GE. mxstep)
GO TO 500
6870 CALL dewset (n, itol, rtol, atol, rwork(lyh), rwork(lewt))
6872 IF (rwork(i+lewt-1) .LE. 0.0d0)
GO TO 510
6873 260 rwork(i+lewt-1) = 1.0d0/rwork(i+lewt-1)
6874 270 tolsf = uround*dmnorm(n, rwork(lyh), rwork(lewt))
6875 IF (tolsf .LE. 1.0d0)
GO TO 280
6877 IF (nst .EQ. 0)
GO TO 626
6879 280
IF ((tn + h) .NE. tn)
GO TO 290
6881 IF (nhnil .GT. mxhnil)
GO TO 290
6882 msg =
'DLSODAR- Warning..Internal T(=R1) and H(=R2) are ' 6883 CALL xerrwd (msg, 50, 101, 0, 0, 0, 0, 0, 0.0d0, 0.0d0)
6884 msg=
' such that in the machine, T + H = T on the next step ' 6885 CALL xerrwd (msg, 60, 101, 0, 0, 0, 0, 0, 0.0d0, 0.0d0)
6886 msg =
' (H = step size). Solver will continue anyway.' 6887 CALL xerrwd (msg, 50, 101, 0, 0, 0, 0, 2, tn, h)
6888 IF (nhnil .LT. mxhnil)
GO TO 290
6889 msg =
'DLSODAR- Above warning has been issued I1 times. ' 6890 CALL xerrwd (msg, 50, 102, 0, 0, 0, 0, 0, 0.0d0, 0.0d0)
6891 msg =
' It will not be issued again for this problem.' 6892 CALL xerrwd (msg, 50, 102, 0, 1, mxhnil, 0, 0, 0.0d0, 0.0d0)
6897 CALL dstoda (neq, y, rwork(lyh), nyh, rwork(lyh), rwork(lewt),
6898 1 rwork(lsavf), rwork(lacor), rwork(lwm), iwork(liwm),
6899 2 f, jac, dprja, dsolsy)
6901 GO TO (300, 530, 540), kgo
6913 IF (meth .EQ. mused)
GO TO 310
6916 IF (meth .EQ. 2) maxord = mxords
6917 IF (meth .EQ. 2) rwork(lwm) = sqrt(uround)
6918 insufr = min(insufr,1)
6919 insufi = min(insufi,1)
6921 IF (ixpr .EQ. 0)
GO TO 310
6922 IF (meth .EQ. 2)
THEN 6923 msg=
'DLSODAR- A switch to the BDF (stiff) method has occurred ' 6924 CALL xerrwd (msg, 60, 105, 0, 0, 0, 0, 0, 0.0d0, 0.0d0)
6926 IF (meth .EQ. 1)
THEN 6927 msg=
'DLSODAR- A switch to the Adams (nonstiff) method occurred ' 6928 CALL xerrwd (msg, 60, 106, 0, 0, 0, 0, 0, 0.0d0, 0.0d0)
6930 msg=
' at T = R1, tentative step size H = R2, step NST = I1 ' 6931 CALL xerrwd (msg, 60, 107, 0, 1, nst, 0, 2, tn, h)
6934 IF (ngc .EQ. 0)
GO TO 315
6935 CALL drchek (3, g, neq, y, rwork(lyh), nyh,
6936 1 rwork(lg0), rwork(lg1), rwork(lgx), jroot, irt)
6937 IF (irt .NE. 1)
GO TO 315
6944 GO TO (320, 400, 330, 340, 350), itask
6946 320
IF ((tn - tout)*h .LT. 0.0d0)
GO TO 250
6947 CALL dintdy (tout, 0, rwork(lyh), nyh, y, iflag)
6951 330
IF ((tn - tout)*h .GE. 0.0d0)
GO TO 400
6954 340
IF ((tn - tout)*h .LT. 0.0d0)
GO TO 345
6955 CALL dintdy (tout, 0, rwork(lyh), nyh, y, iflag)
6958 345 hmx = abs(tn) + abs(h)
6959 ihit = abs(tn - tcrit) .LE. 100.0d0*uround*hmx
6961 tnext = tn + h*(1.0d0 + 4.0d0*uround)
6962 IF ((tnext - tcrit)*h .LE. 0.0d0)
GO TO 250
6963 h = (tcrit - tn)*(1.0d0 - 4.0d0*uround)
6964 IF (jstart .GE. 0) jstart = -2
6967 350 hmx = abs(tn) + abs(h)
6968 ihit = abs(tn - tcrit) .LE. 100.0d0*uround*hmx
6977 410 y(i) = rwork(i+lyh-1)
6979 IF (itask .NE. 4 .AND. itask .NE. 5)
GO TO 420
7006 500 msg =
'DLSODAR- At current T (=R1), MXSTEP (=I1) steps ' 7007 CALL xerrwd (msg, 50, 201, 0, 0, 0, 0, 0, 0.0d0, 0.0d0)
7008 msg =
' taken on this call before reaching TOUT ' 7009 CALL xerrwd (msg, 50, 201, 0, 1, mxstep, 0, 1, tn, 0.0d0)
7013 510 ewti = rwork(lewt+i-1)
7014 msg = .le.
'DLSODAR- At T(=R1), EWT(I1) has become R2 0.' 7015 CALL xerrwd (msg, 50, 202, 0, 1, i, 0, 2, tn, ewti)
7019 520 msg =
'DLSODAR- At T (=R1), too much accuracy requested ' 7020 CALL xerrwd (msg, 50, 203, 0, 0, 0, 0, 0, 0.0d0, 0.0d0)
7021 msg =
' for precision of machine.. See TOLSF (=R2) ' 7022 CALL xerrwd (msg, 50, 203, 0, 0, 0, 0, 2, tn, tolsf)
7027 530 msg =
'DLSODAR- At T(=R1), step size H(=R2), the error ' 7028 CALL xerrwd (msg, 50, 204, 0, 0, 0, 0, 0, 0.0d0, 0.0d0)
7029 msg =
' test failed repeatedly or with ABS(H) = HMIN' 7030 CALL xerrwd (msg, 50, 204, 0, 0, 0, 0, 2, tn, h)
7034 540 msg =
'DLSODAR- At T (=R1) and step size H (=R2), the ' 7035 CALL xerrwd (msg, 50, 205, 0, 0, 0, 0, 0, 0.0d0, 0.0d0)
7036 msg =
' corrector convergence failed repeatedly ' 7037 CALL xerrwd (msg, 50, 205, 0, 0, 0, 0, 0, 0.0d0, 0.0d0)
7038 msg =
' or with ABS(H) = HMIN ' 7039 CALL xerrwd (msg, 30, 205, 0, 0, 0, 0, 2, tn, h)
7043 550 msg =
'DLSODAR- At current T(=R1), RWORK length too small' 7044 CALL xerrwd (msg, 50, 206, 0, 0, 0, 0, 0, 0.0d0, 0.0d0)
7045 msg=
' to proceed. The integration was otherwise successful.' 7046 CALL xerrwd (msg, 60, 206, 0, 0, 0, 0, 1, tn, 0.0d0)
7050 555 msg =
'DLSODAR- At current T(=R1), IWORK length too small' 7051 CALL xerrwd (msg, 50, 207, 0, 0, 0, 0, 0, 0.0d0, 0.0d0)
7052 msg=
' to proceed. The integration was otherwise successful.' 7053 CALL xerrwd (msg, 60, 207, 0, 0, 0, 0, 1, tn, 0.0d0)
7060 SIZE = abs(rwork(i+lacor-1)*rwork(i+lewt-1))
7061 IF (big .GE. size)
GO TO 570
7068 590 y(i) = rwork(i+lyh-1)
7091 601 msg =
'DLSODAR- ISTATE(=I1) illegal.' 7092 CALL xerrwd (msg, 30, 1, 0, 1, istate, 0, 0, 0.0d0, 0.0d0)
7093 IF (istate .LT. 0)
GO TO 800
7095 602 msg =
'DLSODAR- ITASK (=I1) illegal.' 7096 CALL xerrwd (msg, 30, 2, 0, 1, itask, 0, 0, 0.0d0, 0.0d0)
7098 603 msg = .gt.
'DLSODAR- ISTATE1 but DLSODAR not initialized.' 7099 CALL xerrwd (msg, 50, 3, 0, 0, 0, 0, 0, 0.0d0, 0.0d0)
7101 604 msg = .lt.
'DLSODAR- NEQ (=I1) 1 ' 7102 CALL xerrwd (msg, 30, 4, 0, 1, neq(1), 0, 0, 0.0d0, 0.0d0)
7104 605 msg =
'DLSODAR- ISTATE = 3 and NEQ increased (I1 to I2).' 7105 CALL xerrwd (msg, 50, 5, 0, 2, n, neq(1), 0, 0.0d0, 0.0d0)
7107 606 msg =
'DLSODAR- ITOL (=I1) illegal. ' 7108 CALL xerrwd (msg, 30, 6, 0, 1, itol, 0, 0, 0.0d0, 0.0d0)
7110 607 msg =
'DLSODAR- IOPT (=I1) illegal. ' 7111 CALL xerrwd (msg, 30, 7, 0, 1, iopt, 0, 0, 0.0d0, 0.0d0)
7113 608 msg =
'DLSODAR- JT (=I1) illegal. ' 7114 CALL xerrwd (msg, 30, 8, 0, 1, jt, 0, 0, 0.0d0, 0.0d0)
7116 609 msg = .lt..ge.
'DLSODAR- ML (=I1) illegal: 0 or NEQ (=I2)' 7117 CALL xerrwd (msg, 50, 9, 0, 2, ml, neq(1), 0, 0.0d0, 0.0d0)
7119 610 msg = .lt..ge.
'DLSODAR- MU (=I1) illegal: 0 or NEQ (=I2)' 7120 CALL xerrwd (msg, 50, 10, 0, 2, mu, neq(1), 0, 0.0d0, 0.0d0)
7122 611 msg =
'DLSODAR- IXPR (=I1) illegal. ' 7123 CALL xerrwd (msg, 30, 11, 0, 1, ixpr, 0, 0, 0.0d0, 0.0d0)
7125 612 msg = .lt.
'DLSODAR- MXSTEP (=I1) 0 ' 7126 CALL xerrwd (msg, 30, 12, 0, 1, mxstep, 0, 0, 0.0d0, 0.0d0)
7128 613 msg = .lt.
'DLSODAR- MXHNIL (=I1) 0 ' 7129 CALL xerrwd (msg, 30, 13, 0, 1, mxhnil, 0, 0, 0.0d0, 0.0d0)
7131 614 msg =
'DLSODAR- TOUT (=R1) behind T (=R2) ' 7132 CALL xerrwd (msg, 40, 14, 0, 0, 0, 0, 2, tout, t)
7133 msg =
' Integration direction is given by H0 (=R1) ' 7134 CALL xerrwd (msg, 50, 14, 0, 0, 0, 0, 1, h0, 0.0d0)
7136 615 msg = .lt.
'DLSODAR- HMAX (=R1) 0.0 ' 7137 CALL xerrwd (msg, 30, 15, 0, 0, 0, 0, 1, hmax, 0.0d0)
7139 616 msg = .lt.
'DLSODAR- HMIN (=R1) 0.0 ' 7140 CALL xerrwd (msg, 30, 16, 0, 0, 0, 0, 1, hmin, 0.0d0)
7142 617 msg=
'DLSODAR- RWORK length needed, LENRW(=I1), exceeds LRW(=I2) ' 7143 CALL xerrwd (msg, 60, 17, 0, 2, lenrw, lrw, 0, 0.0d0, 0.0d0)
7145 618 msg=
'DLSODAR- IWORK length needed, LENIW(=I1), exceeds LIW(=I2) ' 7146 CALL xerrwd (msg, 60, 18, 0, 2, leniw, liw, 0, 0.0d0, 0.0d0)
7148 619 msg = .lt.
'DLSODAR- RTOL(I1) is R1 0.0 ' 7149 CALL xerrwd (msg, 40, 19, 0, 1, i, 0, 1, rtoli, 0.0d0)
7151 620 msg = .lt.
'DLSODAR- ATOL(I1) is R1 0.0 ' 7152 CALL xerrwd (msg, 40, 20, 0, 1, i, 0, 1, atoli, 0.0d0)
7154 621 ewti = rwork(lewt+i-1)
7155 msg = .le.
'DLSODAR- EWT(I1) is R1 0.0 ' 7156 CALL xerrwd (msg, 40, 21, 0, 1, i, 0, 1, ewti, 0.0d0)
7158 622 msg=
'DLSODAR- TOUT(=R1) too close to T(=R2) to start integration.' 7159 CALL xerrwd (msg, 60, 22, 0, 0, 0, 0, 2, tout, t)
7161 623 msg=
'DLSODAR- ITASK = I1 and TOUT (=R1) behind TCUR - HU (= R2) ' 7162 CALL xerrwd (msg, 60, 23, 0, 1, itask, 0, 2, tout, tp)
7164 624 msg=
'DLSODAR- ITASK = 4 or 5 and TCRIT (=R1) behind TCUR (=R2) ' 7165 CALL xerrwd (msg, 60, 24, 0, 0, 0, 0, 2, tcrit, tn)
7167 625 msg=
'DLSODAR- ITASK = 4 or 5 and TCRIT (=R1) behind TOUT (=R2) ' 7168 CALL xerrwd (msg, 60, 25, 0, 0, 0, 0, 2, tcrit, tout)
7170 626 msg =
'DLSODAR- At start of problem, too much accuracy ' 7171 CALL xerrwd (msg, 50, 26, 0, 0, 0, 0, 0, 0.0d0, 0.0d0)
7172 msg=
' requested for precision of machine.. See TOLSF (=R1) ' 7173 CALL xerrwd (msg, 60, 26, 0, 0, 0, 0, 1, tolsf, 0.0d0)
7176 627 msg =
'DLSODAR- Trouble in DINTDY. ITASK = I1, TOUT = R1' 7177 CALL xerrwd (msg, 50, 27, 0, 1, itask, 0, 1, tout, 0.0d0)
7179 628 msg = .lt.
'DLSODAR- MXORDN (=I1) 0 ' 7180 CALL xerrwd (msg, 30, 28, 0, 1, mxordn, 0, 0, 0.0d0, 0.0d0)
7182 629 msg = .lt.
'DLSODAR- MXORDS (=I1) 0 ' 7183 CALL xerrwd (msg, 30, 29, 0, 1, mxords, 0, 0, 0.0d0, 0.0d0)
7185 630 msg = .lt.
'DLSODAR- NG (=I1) 0 ' 7186 CALL xerrwd (msg, 30, 30, 0, 1, ng, 0, 0, 0.0d0, 0.0d0)
7188 631 msg =
'DLSODAR- NG changed (from I1 to I2) illegally, ' 7189 CALL xerrwd (msg, 50, 31, 0, 0, 0, 0, 0, 0.0d0, 0.0d0)
7190 msg =
' i.e. not immediately after a root was found.' 7191 CALL xerrwd (msg, 50, 31, 0, 2, ngc, ng, 0, 0.0d0, 0.0d0)
7193 632 msg =
'DLSODAR- One or more components of g has a root ' 7194 CALL xerrwd (msg, 50, 32, 0, 0, 0, 0, 0, 0.0d0, 0.0d0)
7195 msg =
' too near to the initial point. ' 7196 CALL xerrwd (msg, 40, 32, 0, 0, 0, 0, 0, 0.0d0, 0.0d0)
7201 800 msg =
'DLSODAR- Run aborted.. apparent infinite loop. ' 7202 CALL xerrwd (msg, 50, 303, 2, 0, 0, 0, 0, 0.0d0, 0.0d0)
7207 SUBROUTINE dlsodpk (F, NEQ, Y, T, TOUT, ITOL, RTOL, ATOL, ITASK,
7208 1 ISTATE, IOPT, RWORK, LRW, IWORK, LIW, JAC, PSOL, MF)
7209 EXTERNAL f, jac, psol
7210 INTEGER neq, itol, itask, istate, iopt, lrw, iwork, liw, mf
7211 DOUBLE PRECISION y, t, tout, rtol, atol, rwork
7212 dimension neq(*), y(*), rtol(*), atol(*), rwork(lrw), iwork(liw)
8312 DOUBLE PRECISION dumach, dvnorm
8313 INTEGER init, mxstep, mxhnil, nhnil, nslast, nyh, iowns,
8314 1 icf, ierpj, iersl, jcur, jstart, kflag, l,
8315 2 lyh, lewt, lacor, lsavf, lwm, liwm, meth, miter,
8316 3 maxord, maxcor, msbp, mxncf, n, nq, nst, nfe, nje, nqu
8317 INTEGER jpre, jacflg, locwp, lociwp, lsavx, kmp, maxl, mnewt,
8318 1 nni, nli, nps, ncfn, ncfl
8319 INTEGER i, i1, i2, iflag, imxer, kgo, lf0, leniw,
8320 1 leniwk, lenrw, lenwm, lenwk, liwp, lwp, mord, mxhnl0, mxstp0,
8321 2 ncfn0, ncfl0, nli0, nni0, nnid, nstd, nwarn
8322 DOUBLE PRECISION rowns,
8323 1 ccmax, el0, h, hmin, hmxi, hu, rc, tn, uround
8324 DOUBLE PRECISION delt, epcon, sqrtn, rsqrtn
8325 DOUBLE PRECISION atoli, avdim, ayi, big, ewti, h0, hmax, hmx,
8326 1 rcfl, rcfn, rh, rtoli, tcrit,
8327 2 tdist, tnext, tol, tolsf, tp,
SIZE, sum, w0
8329 LOGICAL ihit, lavd, lcfn, lcfl, lwarn
8331 SAVE mord, mxstp0, mxhnl0
8344 COMMON /dls001/ rowns(209),
8345 1 ccmax, el0, h, hmin, hmxi, hu, rc, tn, uround,
8346 2 init, mxstep, mxhnil, nhnil, nslast, nyh, iowns(6),
8347 3 icf, ierpj, iersl, jcur, jstart, kflag, l,
8348 4 lyh, lewt, lacor, lsavf, lwm, liwm, meth, miter,
8349 5 maxord, maxcor, msbp, mxncf, n, nq, nst, nfe, nje, nqu
8351 COMMON /dlpk01/ delt, epcon, sqrtn, rsqrtn,
8352 1 jpre, jacflg, locwp, lociwp, lsavx, kmp, maxl, mnewt,
8353 2 nni, nli, nps, ncfn, ncfl
8355 DATA mord(1),mord(2)/12,5/, mxstp0/500/, mxhnl0/10/
8364 IF (istate .LT. 1 .OR. istate .GT. 3)
GO TO 601
8365 IF (itask .LT. 1 .OR. itask .GT. 5)
GO TO 602
8366 IF (istate .EQ. 1)
GO TO 10
8367 IF (init .EQ. 0)
GO TO 603
8368 IF (istate .EQ. 2)
GO TO 200
8371 IF (tout .EQ. t)
RETURN 8380 20
IF (neq(1) .LE. 0)
GO TO 604
8381 IF (istate .EQ. 1)
GO TO 25
8382 IF (neq(1) .GT. n)
GO TO 605
8384 IF (itol .LT. 1 .OR. itol .GT. 4)
GO TO 606
8385 IF (iopt .LT. 0 .OR. iopt .GT. 1)
GO TO 607
8387 miter = mf - 10*meth
8388 IF (meth .LT. 1 .OR. meth .GT. 2)
GO TO 608
8389 IF (miter .LT. 0)
GO TO 608
8390 IF (miter .GT. 4 .AND. miter .LT. 9)
GO TO 608
8391 IF (miter .GE. 1) jpre = iwork(3)
8393 IF (miter .GE. 1) jacflg = iwork(4)
8395 IF (iopt .EQ. 1)
GO TO 40
8399 IF (istate .EQ. 1) h0 = 0.0d0
8406 40 maxord = iwork(5)
8407 IF (maxord .LT. 0)
GO TO 611
8408 IF (maxord .EQ. 0) maxord = 100
8409 maxord = min(maxord,mord(meth))
8411 IF (mxstep .LT. 0)
GO TO 612
8412 IF (mxstep .EQ. 0) mxstep = mxstp0
8414 IF (mxhnil .LT. 0)
GO TO 613
8415 IF (mxhnil .EQ. 0) mxhnil = mxhnl0
8416 IF (istate .NE. 1)
GO TO 50
8418 IF ((tout - t)*h0 .LT. 0.0d0)
GO TO 614
8420 IF (hmax .LT. 0.0d0)
GO TO 615
8422 IF (hmax .GT. 0.0d0) hmxi = 1.0d0/hmax
8424 IF (hmin .LT. 0.0d0)
GO TO 616
8426 IF (maxl .EQ. 0) maxl = 5
8429 IF (kmp .EQ. 0 .OR. kmp .GT. maxl) kmp = maxl
8431 IF (delt .EQ. 0.0d0) delt = 0.05d0
8439 IF (istate .EQ. 1) nyh = n
8440 lwm = lyh + (maxord + 1)*nyh
8441 IF (miter .EQ. 0) lenwk = 0
8442 IF (miter .EQ. 1) lenwk = n*(maxl+2) + maxl*maxl
8444 1 lenwk = n*(maxl+2+min(1,maxl-kmp)) + (maxl+3)*maxl + 1
8445 IF (miter .EQ. 3 .OR. miter .EQ. 4) lenwk = 5*n
8446 IF (miter .EQ. 9) lenwk = 2*n
8448 IF (miter .GE. 1) lwp = iwork(1)
8455 IF (miter .EQ. 0) lacor = lsavf + n
8456 lenrw = lacor + n - 1
8460 IF (miter .EQ. 1) leniwk = maxl
8462 IF (miter .GE. 1) liwp = iwork(2)
8463 leniw = 30 + leniwk + liwp
8466 IF (lenrw .GT. lrw)
GO TO 617
8467 IF (leniw .GT. liw)
GO TO 618
8472 IF (itol .GE. 3) rtoli = rtol(i)
8473 IF (itol .EQ. 2 .OR. itol .EQ. 4) atoli = atol(i)
8474 IF (rtoli .LT. 0.0d0)
GO TO 619
8475 IF (atoli .LT. 0.0d0)
GO TO 620
8478 sqrtn = sqrt(
REAL(n))
8479 rsqrtn = 1.0d0/sqrtn
8480 IF (istate .EQ. 1)
GO TO 100
8483 IF (nq .LE. maxord)
GO TO 90
8486 80 rwork(i+lsavf-1) = rwork(i+lwm-1)
8488 IF (n .EQ. nyh)
GO TO 200
8491 i2 = lyh + (maxord + 1)*nyh - 1
8492 IF (i1 .GT. i2)
GO TO 200
8503 100 uround = dumach()
8505 IF (itask .NE. 4 .AND. itask .NE. 5)
GO TO 110
8507 IF ((tcrit - tout)*(tout - t) .LT. 0.0d0)
GO TO 625
8508 IF (h0 .NE. 0.0d0 .AND. (t + h0 - tcrit)*h0 .GT. 0.0d0)
8533 CALL f (neq, t, y, rwork(lf0))
8537 115 rwork(i+lyh-1) = y(i)
8541 CALL dewset (n, itol, rtol, atol, rwork(lyh), rwork(lewt))
8543 IF (rwork(i+lewt-1) .LE. 0.0d0)
GO TO 621
8544 120 rwork(i+lewt-1) = 1.0d0/rwork(i+lewt-1)
8561 IF (h0 .NE. 0.0d0)
GO TO 180
8562 tdist = abs(tout - t)
8563 w0 = max(abs(t),abs(tout))
8564 IF (tdist .LT. 2.0d0*uround*w0)
GO TO 622
8566 IF (itol .LE. 2)
GO TO 140
8568 130 tol = max(tol,rtol(i))
8569 140
IF (tol .GT. 0.0d0)
GO TO 160
8572 IF (itol .EQ. 2 .OR. itol .EQ. 4) atoli = atol(i)
8574 IF (ayi .NE. 0.0d0) tol = max(tol,atoli/ayi)
8576 160 tol = max(tol,100.0d0*uround)
8577 tol = min(tol,0.001d0)
8578 sum = dvnorm(n, rwork(lf0), rwork(lewt))
8579 sum = 1.0d0/(tol*w0*w0) + tol*sum**2
8580 h0 = 1.0d0/sqrt(sum)
8582 h0 = sign(h0,tout-t)
8584 180 rh = abs(h0)*hmxi
8585 IF (rh .GT. 1.0d0) h0 = h0/rh
8589 190 rwork(i+lf0-1) = h0*rwork(i+lf0-1)
8602 GO TO (210, 250, 220, 230, 240), itask
8603 210
IF ((tn - tout)*h .LT. 0.0d0)
GO TO 250
8604 CALL dintdy (tout, 0, rwork(lyh), nyh, y, iflag)
8605 IF (iflag .NE. 0)
GO TO 627
8608 220 tp = tn - hu*(1.0d0 + 100.0d0*uround)
8609 IF ((tp - tout)*h .GT. 0.0d0)
GO TO 623
8610 IF ((tn - tout)*h .LT. 0.0d0)
GO TO 250
8612 230 tcrit = rwork(1)
8613 IF ((tn - tcrit)*h .GT. 0.0d0)
GO TO 624
8614 IF ((tcrit - tout)*h .LT. 0.0d0)
GO TO 625
8615 IF ((tn - tout)*h .LT. 0.0d0)
GO TO 245
8616 CALL dintdy (tout, 0, rwork(lyh), nyh, y, iflag)
8617 IF (iflag .NE. 0)
GO TO 627
8620 240 tcrit = rwork(1)
8621 IF ((tn - tcrit)*h .GT. 0.0d0)
GO TO 624
8622 245 hmx = abs(tn) + abs(h)
8623 ihit = abs(tn - tcrit) .LE. 100.0d0*uround*hmx
8625 tnext = tn + h*(1.0d0 + 4.0d0*uround)
8626 IF ((tnext - tcrit)*h .LE. 0.0d0)
GO TO 250
8627 h = (tcrit - tn)*(1.0d0 - 4.0d0*uround)
8628 IF (istate .EQ. 2) jstart = -2
8642 IF ((nst-nslast) .GE. mxstep)
GO TO 500
8645 IF (nstd .LT. 10 .OR. nnid .EQ. 0)
GO TO 255
8646 avdim =
REAL(nli - nli0)/
REAL(nnid)
8647 rcfn =
REAL(ncfn - ncfn0)/
REAL(nstd)
8648 rcfl =
REAL(ncfl - ncfl0)/
REAL(nnid)
8649 lavd = avdim .GT. (maxl - 0.05d0)
8650 lcfn = rcfn .GT. 0.9d0
8651 lcfl = rcfl .GT. 0.9d0
8652 lwarn = lavd .OR. lcfn .OR. lcfl
8653 IF (.NOT.lwarn)
GO TO 255
8655 IF (nwarn .GT. 10)
GO TO 255
8657 msg=
'DLSODPK- Warning. Poor iterative algorithm performance seen ' 8658 CALL xerrwd (msg, 60, 111, 0, 0, 0, 0, 0, 0.0d0, 0.0d0)
8661 msg=
' at T = R1 by average no. of linear iterations = R2 ' 8662 CALL xerrwd (msg, 60, 111, 0, 0, 0, 0, 2, tn, avdim)
8665 msg=
'DLSODPK- Warning. Poor iterative algorithm performance seen ' 8666 CALL xerrwd (msg, 60, 112, 0, 0, 0, 0, 0, 0.0d0, 0.0d0)
8669 msg=
' at T = R1 by nonlinear convergence failure rate = R2 ' 8670 CALL xerrwd (msg, 60, 112, 0, 0, 0, 0, 2, tn, rcfn)
8673 msg=
'DLSODPK- Warning. Poor iterative algorithm performance seen ' 8674 CALL xerrwd (msg, 60, 113, 0, 0, 0, 0, 0, 0.0d0, 0.0d0)
8677 msg=
' at T = R1 by linear convergence failure rate = R2 ' 8678 CALL xerrwd (msg, 60, 113, 0, 0, 0, 0, 2, tn, rcfl)
8681 CALL dewset (n, itol, rtol, atol, rwork(lyh), rwork(lewt))
8683 IF (rwork(i+lewt-1) .LE. 0.0d0)
GO TO 510
8684 260 rwork(i+lewt-1) = 1.0d0/rwork(i+lewt-1)
8685 270 tolsf = uround*dvnorm(n, rwork(lyh), rwork(lewt))
8686 IF (tolsf .LE. 1.0d0)
GO TO 280
8688 IF (nst .EQ. 0)
GO TO 626
8690 280
IF ((tn + h) .NE. tn)
GO TO 290
8692 IF (nhnil .GT. mxhnil)
GO TO 290
8693 msg =
'DLSODPK- Warning..Internal T(=R1) and H(=R2) are ' 8694 CALL xerrwd (msg, 50, 101, 0, 0, 0, 0, 0, 0.0d0, 0.0d0)
8695 msg=
' such that in the machine, T + H = T on the next step ' 8696 CALL xerrwd (msg, 60, 101, 0, 0, 0, 0, 0, 0.0d0, 0.0d0)
8697 msg =
' (H = step size). Solver will continue anyway.' 8698 CALL xerrwd (msg, 50, 101, 0, 0, 0, 0, 2, tn, h)
8699 IF (nhnil .LT. mxhnil)
GO TO 290
8700 msg =
'DLSODPK- Above warning has been issued I1 times. ' 8701 CALL xerrwd (msg, 50, 102, 0, 0, 0, 0, 0, 0.0d0, 0.0d0)
8702 msg =
' It will not be issued again for this problem.' 8703 CALL xerrwd (msg, 50, 102, 0, 1, mxhnil, 0, 0, 0.0d0, 0.0d0)
8708 CALL dstodpk (neq, y, rwork(lyh), nyh, rwork(lyh), rwork(lewt),
8709 1 rwork(lsavf), rwork(lsavx), rwork(lacor), rwork(lwm),
8710 2 iwork(liwm), f, jac, psol)
8712 GO TO (300, 530, 540, 550), kgo
8719 GO TO (310, 400, 330, 340, 350), itask
8721 310
IF ((tn - tout)*h .LT. 0.0d0)
GO TO 250
8722 CALL dintdy (tout, 0, rwork(lyh), nyh, y, iflag)
8726 330
IF ((tn - tout)*h .GE. 0.0d0)
GO TO 400
8729 340
IF ((tn - tout)*h .LT. 0.0d0)
GO TO 345
8730 CALL dintdy (tout, 0, rwork(lyh), nyh, y, iflag)
8733 345 hmx = abs(tn) + abs(h)
8734 ihit = abs(tn - tcrit) .LE. 100.0d0*uround*hmx
8736 tnext = tn + h*(1.0d0 + 4.0d0*uround)
8737 IF ((tnext - tcrit)*h .LE. 0.0d0)
GO TO 250
8738 h = (tcrit - tn)*(1.0d0 - 4.0d0*uround)
8742 350 hmx = abs(tn) + abs(h)
8743 ihit = abs(tn - tcrit) .LE. 100.0d0*uround*hmx
8752 410 y(i) = rwork(i+lyh-1)
8754 IF (itask .NE. 4 .AND. itask .NE. 5)
GO TO 420
8780 500 msg =
'DLSODPK- At current T (=R1), MXSTEP (=I1) steps ' 8781 CALL xerrwd (msg, 50, 201, 0, 0, 0, 0, 0, 0.0d0, 0.0d0)
8782 msg =
' taken on this call before reaching TOUT ' 8783 CALL xerrwd (msg, 50, 201, 0, 1, mxstep, 0, 1, tn, 0.0d0)
8787 510 ewti = rwork(lewt+i-1)
8788 msg = .le.
'DLSODPK- At T (=R1), EWT(I1) has become R20. ' 8789 CALL xerrwd (msg, 50, 202, 0, 1, i, 0, 2, tn, ewti)
8793 520 msg =
'DLSODPK- At T (=R1), too much accuracy requested ' 8794 CALL xerrwd (msg, 50, 203, 0, 0, 0, 0, 0, 0.0d0, 0.0d0)
8795 msg =
' for precision of machine.. See TOLSF (=R2) ' 8796 CALL xerrwd (msg, 50, 203, 0, 0, 0, 0, 2, tn, tolsf)
8801 530 msg =
'DLSODPK- At T(=R1), step size H(=R2), the error ' 8802 CALL xerrwd (msg, 50, 204, 0, 0, 0, 0, 0, 0.0d0, 0.0d0)
8803 msg =
' test failed repeatedly or with ABS(H) = HMIN' 8804 CALL xerrwd (msg, 50, 204, 0, 0, 0, 0, 2, tn, h)
8808 540 msg =
'DLSODPK- At T (=R1) and step size H (=R2), the ' 8809 CALL xerrwd (msg, 50, 205, 0, 0, 0, 0, 0, 0.0d0, 0.0d0)
8810 msg =
' corrector convergence failed repeatedly ' 8811 CALL xerrwd (msg, 50, 205, 0, 0, 0, 0, 0, 0.0d0, 0.0d0)
8812 msg =
' or with ABS(H) = HMIN ' 8813 CALL xerrwd (msg, 30, 205, 0, 0, 0, 0, 2, tn, h)
8817 550 msg =
'DLSODPK- At T (=R1) an unrecoverable error return' 8818 CALL xerrwd (msg, 50, 205, 0, 0, 0, 0, 0, 0.0d0, 0.0d0)
8819 msg =
' was made from Subroutine PSOL ' 8820 CALL xerrwd (msg, 40, 205, 0, 0, 0, 0, 1, tn, 0.0d0)
8827 SIZE = abs(rwork(i+lacor-1)*rwork(i+lewt-1))
8828 IF (big .GE. size)
GO TO 570
8835 590 y(i) = rwork(i+lyh-1)
8858 601 msg =
'DLSODPK- ISTATE(=I1) illegal.' 8859 CALL xerrwd (msg, 30, 1, 0, 1, istate, 0, 0, 0.0d0, 0.0d0)
8860 IF (istate .LT. 0)
GO TO 800
8862 602 msg =
'DLSODPK- ITASK (=I1) illegal.' 8863 CALL xerrwd (msg, 30, 2, 0, 1, itask, 0, 0, 0.0d0, 0.0d0)
8865 603 msg = .gt.
'DLSODPK- ISTATE1 but DLSODPK not initialized.' 8866 CALL xerrwd (msg, 50, 3, 0, 0, 0, 0, 0, 0.0d0, 0.0d0)
8868 604 msg = .lt.
'DLSODPK- NEQ (=I1) 1 ' 8869 CALL xerrwd (msg, 30, 4, 0, 1, neq(1), 0, 0, 0.0d0, 0.0d0)
8871 605 msg =
'DLSODPK- ISTATE = 3 and NEQ increased (I1 to I2).' 8872 CALL xerrwd (msg, 50, 5, 0, 2, n, neq(1), 0, 0.0d0, 0.0d0)
8874 606 msg =
'DLSODPK- ITOL (=I1) illegal. ' 8875 CALL xerrwd (msg, 30, 6, 0, 1, itol, 0, 0, 0.0d0, 0.0d0)
8877 607 msg =
'DLSODPK- IOPT (=I1) illegal. ' 8878 CALL xerrwd (msg, 30, 7, 0, 1, iopt, 0, 0, 0.0d0, 0.0d0)
8880 608 msg =
'DLSODPK- MF (=I1) illegal. ' 8881 CALL xerrwd (msg, 30, 8, 0, 1, mf, 0, 0, 0.0d0, 0.0d0)
8883 611 msg = .lt.
'DLSODPK- MAXORD (=I1) 0 ' 8884 CALL xerrwd (msg, 30, 11, 0, 1, maxord, 0, 0, 0.0d0, 0.0d0)
8886 612 msg = .lt.
'DLSODPK- MXSTEP (=I1) 0 ' 8887 CALL xerrwd (msg, 30, 12, 0, 1, mxstep, 0, 0, 0.0d0, 0.0d0)
8889 613 msg = .lt.
'DLSODPK- MXHNIL (=I1) 0 ' 8890 CALL xerrwd (msg, 30, 13, 0, 1, mxhnil, 0, 0, 0.0d0, 0.0d0)
8892 614 msg =
'DLSODPK- TOUT (=R1) behind T (=R2) ' 8893 CALL xerrwd (msg, 40, 14, 0, 0, 0, 0, 2, tout, t)
8894 msg =
' Integration direction is given by H0 (=R1) ' 8895 CALL xerrwd (msg, 50, 14, 0, 0, 0, 0, 1, h0, 0.0d0)
8897 615 msg = .lt.
'DLSODPK- HMAX (=R1) 0.0 ' 8898 CALL xerrwd (msg, 30, 15, 0, 0, 0, 0, 1, hmax, 0.0d0)
8900 616 msg = .lt.
'DLSODPK- HMIN (=R1) 0.0 ' 8901 CALL xerrwd (msg, 30, 16, 0, 0, 0, 0, 1, hmin, 0.0d0)
8903 617 msg=
'DLSODPK- RWORK length needed, LENRW(=I1), exceeds LRW(=I2) ' 8904 CALL xerrwd (msg, 60, 17, 0, 2, lenrw, lrw, 0, 0.0d0, 0.0d0)
8906 618 msg=
'DLSODPK- IWORK length needed, LENIW(=I1), exceeds LIW(=I2) ' 8907 CALL xerrwd (msg, 60, 18, 0, 2, leniw, liw, 0, 0.0d0, 0.0d0)
8909 619 msg = .lt.
'DLSODPK- RTOL(I1) is R1 0.0 ' 8910 CALL xerrwd (msg, 40, 19, 0, 1, i, 0, 1, rtoli, 0.0d0)
8912 620 msg = .lt.
'DLSODPK- ATOL(I1) is R1 0.0 ' 8913 CALL xerrwd (msg, 40, 20, 0, 1, i, 0, 1, atoli, 0.0d0)
8915 621 ewti = rwork(lewt+i-1)
8916 msg = .le.
'DLSODPK- EWT(I1) is R1 0.0 ' 8917 CALL xerrwd (msg, 40, 21, 0, 1, i, 0, 1, ewti, 0.0d0)
8919 622 msg=
'DLSODPK- TOUT(=R1) too close to T(=R2) to start integration.' 8920 CALL xerrwd (msg, 60, 22, 0, 0, 0, 0, 2, tout, t)
8922 623 msg=
'DLSODPK- ITASK = I1 and TOUT (=R1) behind TCUR - HU (= R2) ' 8923 CALL xerrwd (msg, 60, 23, 0, 1, itask, 0, 2, tout, tp)
8925 624 msg=
'DLSODPK- ITASK = 4 or 5 and TCRIT (=R1) behind TCUR (=R2) ' 8926 CALL xerrwd (msg, 60, 24, 0, 0, 0, 0, 2, tcrit, tn)
8928 625 msg=
'DLSODPK- ITASK = 4 or 5 and TCRIT (=R1) behind TOUT (=R2) ' 8929 CALL xerrwd (msg, 60, 25, 0, 0, 0, 0, 2, tcrit, tout)
8931 626 msg =
'DLSODPK- At start of problem, too much accuracy ' 8932 CALL xerrwd (msg, 50, 26, 0, 0, 0, 0, 0, 0.0d0, 0.0d0)
8933 msg=
' requested for precision of machine.. See TOLSF (=R1) ' 8934 CALL xerrwd (msg, 60, 26, 0, 0, 0, 0, 1, tolsf, 0.0d0)
8937 627 msg =
'DLSODPK- Trouble in DINTDY. ITASK = I1, TOUT = R1' 8938 CALL xerrwd (msg, 50, 27, 0, 1, itask, 0, 1, tout, 0.0d0)
8943 800 msg =
'DLSODPK- Run aborted.. apparent infinite loop. ' 8944 CALL xerrwd (msg, 50, 303, 2, 0, 0, 0, 0, 0.0d0, 0.0d0)
8949 SUBROUTINE dlsodkr (F, NEQ, Y, T, TOUT, ITOL, RTOL, ATOL, ITASK,
8950 1 ISTATE, IOPT, RWORK, LRW, IWORK, LIW, JAC, PSOL,
8952 EXTERNAL f, jac, psol, g
8953 INTEGER neq, itol, itask, istate, iopt, lrw, iwork, liw, mf,
8955 DOUBLE PRECISION y, t, tout, rtol, atol, rwork
8956 dimension neq(*), y(*), rtol(*), atol(*), rwork(lrw), iwork(liw),
10153 DOUBLE PRECISION dumach, dvnorm
10154 INTEGER init, mxstep, mxhnil, nhnil, nslast, nyh, iowns,
10155 1 icf, ierpj, iersl, jcur, jstart, kflag, l,
10156 2 lyh, lewt, lacor, lsavf, lwm, liwm, meth, miter,
10157 3 maxord, maxcor, msbp, mxncf, n, nq, nst, nfe, nje, nqu
10158 INTEGER newt, nsfi, nslj, njev
10159 INTEGER lg0, lg1, lgx, iownr3, irfnd, itaskc, ngc, nge
10160 INTEGER jpre, jacflg, locwp, lociwp, lsavx, kmp, maxl, mnewt,
10161 1 nni, nli, nps, ncfn, ncfl
10162 INTEGER i, i1, i2, ier, iflag, imxer, kgo, lf0,
10163 1 leniw, leniwk, lenrw, lenwm, lenwk, liwp, lwp, mord, mxhnl0,
10164 2 mxstp0, ncfn0, ncfl0, niter, nli0, nni0, nnid, nstd, nwarn
10165 INTEGER irfp, irt, lenyh, lyhnew
10166 DOUBLE PRECISION rowns,
10167 1 ccmax, el0, h, hmin, hmxi, hu, rc, tn, uround
10168 DOUBLE PRECISION stifr
10169 DOUBLE PRECISION rownr3, t0, tlast, toutc
10170 DOUBLE PRECISION delt, epcon, sqrtn, rsqrtn
10171 DOUBLE PRECISION atoli, avdim, big, ewti, h0, hmax, hmx, rcfl,
10172 1 rcfn, rh, rtoli, tcrit, tnext, tolsf, tp,
SIZE 10174 LOGICAL ihit, lavd, lcfn, lcfl, lwarn
10176 SAVE mord, mxstp0, mxhnl0
10191 COMMON /dls001/ rowns(209),
10192 1 ccmax, el0, h, hmin, hmxi, hu, rc, tn, uround,
10193 2 init, mxstep, mxhnil, nhnil, nslast, nyh, iowns(6),
10194 3 icf, ierpj, iersl, jcur, jstart, kflag, l,
10195 4 lyh, lewt, lacor, lsavf, lwm, liwm, meth, miter,
10196 5 maxord, maxcor, msbp, mxncf, n, nq, nst, nfe, nje, nqu
10198 COMMON /dls002/ stifr, newt, nsfi, nslj, njev
10200 COMMON /dlsr01/ rownr3(2), t0, tlast, toutc,
10201 1 lg0, lg1, lgx, iownr3(2), irfnd, itaskc, ngc, nge
10203 COMMON /dlpk01/ delt, epcon, sqrtn, rsqrtn,
10204 1 jpre, jacflg, locwp, lociwp, lsavx, kmp, maxl, mnewt,
10205 2 nni, nli, nps, ncfn, ncfl
10207 DATA mord(1),mord(2)/12,5/, mxstp0/500/, mxhnl0/10/
10216 IF (istate .LT. 1 .OR. istate .GT. 3)
GO TO 601
10217 IF (itask .LT. 1 .OR. itask .GT. 5)
GO TO 602
10219 IF (istate .EQ. 1)
GO TO 10
10220 IF (init .EQ. 0)
GO TO 603
10221 IF (istate .EQ. 2)
GO TO 200
10224 IF (tout .EQ. t)
RETURN 10234 20
IF (neq(1) .LE. 0)
GO TO 604
10235 IF (istate .EQ. 1)
GO TO 25
10236 IF (neq(1) .GT. n)
GO TO 605
10238 IF (itol .LT. 1 .OR. itol .GT. 4)
GO TO 606
10239 IF (iopt .LT. 0 .OR. iopt .GT. 1)
GO TO 607
10241 miter = mf - 10*meth
10242 IF (meth .LT. 1 .OR. meth .GT. 2)
GO TO 608
10243 IF (miter .LT. 0)
GO TO 608
10244 IF (miter .GT. 4 .AND. miter .LT. 9)
GO TO 608
10245 IF (miter .GE. 1) jpre = iwork(3)
10247 IF (miter .GE. 1) jacflg = iwork(4)
10248 IF (ng .LT. 0)
GO TO 630
10249 IF (istate .EQ. 1)
GO TO 35
10250 IF (irfnd .EQ. 0 .AND. ng .NE. ngc)
GO TO 631
10253 IF (iopt .EQ. 1)
GO TO 40
10254 maxord = mord(meth)
10257 IF (istate .EQ. 1) h0 = 0.0d0
10264 40 maxord = iwork(5)
10265 IF (maxord .LT. 0)
GO TO 611
10266 IF (maxord .EQ. 0) maxord = 100
10267 maxord = min(maxord,mord(meth))
10269 IF (mxstep .LT. 0)
GO TO 612
10270 IF (mxstep .EQ. 0) mxstep = mxstp0
10272 IF (mxhnil .LT. 0)
GO TO 613
10273 IF (mxhnil .EQ. 0) mxhnil = mxhnl0
10274 IF (istate .NE. 1)
GO TO 50
10276 IF ((tout - t)*h0 .LT. 0.0d0)
GO TO 614
10278 IF (hmax .LT. 0.0d0)
GO TO 615
10280 IF (hmax .GT. 0.0d0) hmxi = 1.0d0/hmax
10282 IF (hmin .LT. 0.0d0)
GO TO 616
10284 IF (maxl .EQ. 0) maxl = 5
10287 IF (kmp .EQ. 0 .OR. kmp .GT. maxl) kmp = maxl
10289 IF (delt .EQ. 0.0d0) delt = 0.05d0
10297 60
IF (istate .EQ. 1) nyh = n
10302 IF (istate .EQ. 1) lyh = lyhnew
10303 IF (lyhnew .EQ. lyh)
GO TO 62
10306 IF (lrw .LT. lyhnew-1+lenyh)
GO TO 62
10308 IF (lyhnew .GT. lyh) i1 = -1
10309 CALL dcopy (lenyh, rwork(lyh), i1, rwork(lyhnew), i1)
10312 lwm = lyh + (maxord + 1)*nyh
10313 IF (miter .EQ. 0) lenwk = 0
10314 IF (miter .EQ. 1) lenwk = n*(maxl+2) + maxl*maxl
10316 1 lenwk = n*(maxl+2+min(1,maxl-kmp)) + (maxl+3)*maxl + 1
10317 IF (miter .EQ. 3 .OR. miter .EQ. 4) lenwk = 5*n
10318 IF (miter .EQ. 9) lenwk = 2*n
10320 IF (miter .GE. 1) lwp = iwork(1)
10321 lenwm = lenwk + lwp
10327 IF (miter .EQ. 0) lacor = lsavf + n
10328 lenrw = lacor + n - 1
10332 IF (miter .EQ. 1) leniwk = maxl
10334 IF (miter .GE. 1) liwp = iwork(2)
10335 leniw = 30 + leniwk + liwp
10336 lociwp = leniwk + 1
10338 IF (lenrw .GT. lrw)
GO TO 617
10339 IF (leniw .GT. liw)
GO TO 618
10344 IF (itol .GE. 3) rtoli = rtol(i)
10345 IF (itol .EQ. 2 .OR. itol .EQ. 4) atoli = atol(i)
10346 IF (rtoli .LT. 0.0d0)
GO TO 619
10347 IF (atoli .LT. 0.0d0)
GO TO 620
10350 sqrtn = sqrt(
REAL(n))
10351 rsqrtn = 1.0d0/sqrtn
10352 IF (istate .EQ. 1)
GO TO 100
10355 IF (nq .LE. maxord)
GO TO 90
10358 80 rwork(i+lsavf-1) = rwork(i+lwm-1)
10360 IF (n .EQ. nyh)
GO TO 200
10363 i2 = lyh + (maxord + 1)*nyh - 1
10364 IF (i1 .GT. i2)
GO TO 200
10366 95 rwork(i) = 0.0d0
10375 100 uround = dumach()
10377 IF (itask .NE. 4 .AND. itask .NE. 5)
GO TO 110
10379 IF ((tcrit - tout)*(tout - t) .LT. 0.0d0)
GO TO 625
10380 IF (h0 .NE. 0.0d0 .AND. (t + h0 - tcrit)*h0 .GT. 0.0d0)
10407 CALL f (neq, t, y, rwork(lf0))
10411 115 rwork(i+lyh-1) = y(i)
10415 CALL dewset (n, itol, rtol, atol, rwork(lyh), rwork(lewt))
10417 IF (rwork(i+lewt-1) .LE. 0.0d0)
GO TO 621
10418 120 rwork(i+lewt-1) = 1.0d0/rwork(i+lewt-1)
10419 IF (h0 .NE. 0.0d0)
GO TO 180
10421 CALL dlhin (neq, n, t, rwork(lyh), rwork(lf0), f, tout, uround,
10422 1 rwork(lewt), itol, atol, y, rwork(lacor), h0, niter, ier)
10424 IF (ier .NE. 0)
GO TO 622
10426 180 rh = abs(h0)*hmxi
10427 IF (rh .GT. 1.0d0) h0 = h0/rh
10431 190 rwork(i+lf0-1) = h0*rwork(i+lf0-1)
10435 IF (ngc .EQ. 0)
GO TO 270
10436 CALL drchek (1, g, neq, y, rwork(lyh), nyh,
10437 1 rwork(lg0), rwork(lg1), rwork(lgx), jroot, irt)
10438 IF (irt .EQ. 0)
GO TO 270
10452 IF (ngc .EQ. 0)
GO TO 205
10453 IF (itask .EQ. 1 .OR. itask .EQ. 4) toutc = tout
10454 CALL drchek (2, g, neq, y, rwork(lyh), nyh,
10455 1 rwork(lg0), rwork(lg1), rwork(lgx), jroot, irt)
10456 IF (irt .NE. 1)
GO TO 205
10463 IF (irfp .EQ. 1 .AND. tlast .NE. tn .AND. itask .EQ. 2)
GO TO 400
10470 GO TO (210, 250, 220, 230, 240), itask
10471 210
IF ((tn - tout)*h .LT. 0.0d0)
GO TO 250
10472 CALL dintdy (tout, 0, rwork(lyh), nyh, y, iflag)
10473 IF (iflag .NE. 0)
GO TO 627
10476 220 tp = tn - hu*(1.0d0 + 100.0d0*uround)
10477 IF ((tp - tout)*h .GT. 0.0d0)
GO TO 623
10478 IF ((tn - tout)*h .LT. 0.0d0)
GO TO 250
10480 230 tcrit = rwork(1)
10481 IF ((tn - tcrit)*h .GT. 0.0d0)
GO TO 624
10482 IF ((tcrit - tout)*h .LT. 0.0d0)
GO TO 625
10483 IF ((tn - tout)*h .LT. 0.0d0)
GO TO 245
10484 CALL dintdy (tout, 0, rwork(lyh), nyh, y, iflag)
10485 IF (iflag .NE. 0)
GO TO 627
10488 240 tcrit = rwork(1)
10489 IF ((tn - tcrit)*h .GT. 0.0d0)
GO TO 624
10490 245 hmx = abs(tn) + abs(h)
10491 ihit = abs(tn - tcrit) .LE. 100.0d0*uround*hmx
10492 IF (ihit) t = tcrit
10493 IF (irfp .EQ. 1 .AND. tlast .NE. tn .AND. itask .EQ. 5)
GO TO 400
10494 IF (ihit)
GO TO 400
10495 tnext = tn + h*(1.0d0 + 4.0d0*uround)
10496 IF ((tnext - tcrit)*h .LE. 0.0d0)
GO TO 250
10497 h = (tcrit - tn)*(1.0d0 - 4.0d0*uround)
10498 IF (istate .EQ. 2) jstart = -2
10512 IF ((nst-nslast) .GE. mxstep)
GO TO 500
10513 nstd = nst - nslast
10515 IF (nstd .LT. 10 .OR. nnid .EQ. 0)
GO TO 255
10516 avdim =
REAL(nli - nli0)/
REAL(nnid)
10517 rcfn =
REAL(ncfn - ncfn0)/
REAL(nstd)
10518 rcfl =
REAL(ncfl - ncfl0)/
REAL(nnid)
10519 lavd = avdim .GT. (maxl - 0.05d0)
10520 lcfn = rcfn .GT. 0.9d0
10521 lcfl = rcfl .GT. 0.9d0
10522 lwarn = lavd .OR. lcfn .OR. lcfl
10523 IF (.NOT.lwarn)
GO TO 255
10525 IF (nwarn .GT. 10)
GO TO 255
10527 msg=
'DLSODKR- Warning. Poor iterative algorithm performance seen ' 10528 CALL xerrwd (msg, 60, 111, 0, 0, 0, 0, 0, 0.0d0, 0.0d0)
10531 msg=
' at T = R1 by average no. of linear iterations = R2 ' 10532 CALL xerrwd (msg, 60, 111, 0, 0, 0, 0, 2, tn, avdim)
10535 msg=
'DLSODKR- Warning. Poor iterative algorithm performance seen ' 10536 CALL xerrwd (msg, 60, 112, 0, 0, 0, 0, 0, 0.0d0, 0.0d0)
10539 msg=
' at T = R1 by nonlinear convergence failure rate = R2 ' 10540 CALL xerrwd (msg, 60, 112, 0, 0, 0, 0, 2, tn, rcfn)
10543 msg=
'DLSODKR- Warning. Poor iterative algorithm performance seen ' 10544 CALL xerrwd (msg, 60, 113, 0, 0, 0, 0, 0, 0.0d0, 0.0d0)
10547 msg=
' at T = R1 by linear convergence failure rate = R2 ' 10548 CALL xerrwd (msg, 60, 113, 0, 0, 0, 0, 2, tn, rcfl)
10551 CALL dewset (n, itol, rtol, atol, rwork(lyh), rwork(lewt))
10553 IF (rwork(i+lewt-1) .LE. 0.0d0)
GO TO 510
10554 260 rwork(i+lewt-1) = 1.0d0/rwork(i+lewt-1)
10555 270 tolsf = uround*dvnorm(n, rwork(lyh), rwork(lewt))
10556 IF (tolsf .LE. 1.0d0)
GO TO 280
10557 tolsf = tolsf*2.0d0
10558 IF (nst .EQ. 0)
GO TO 626
10560 280
IF ((tn + h) .NE. tn)
GO TO 290
10562 IF (nhnil .GT. mxhnil)
GO TO 290
10563 msg =
'DLSODKR- Warning.. Internal T(=R1) and H(=R2) are' 10564 CALL xerrwd (msg, 50, 101, 0, 0, 0, 0, 0, 0.0d0, 0.0d0)
10565 msg=
' such that in the machine, T + H = T on the next step ' 10566 CALL xerrwd (msg, 60, 101, 0, 0, 0, 0, 0, 0.0d0, 0.0d0)
10567 msg =
' (H = step size). Solver will continue anyway.' 10568 CALL xerrwd (msg, 50, 101, 0, 0, 0, 0, 2, tn, h)
10569 IF (nhnil .LT. mxhnil)
GO TO 290
10570 msg =
'DLSODKR- Above warning has been issued I1 times. ' 10571 CALL xerrwd (msg, 50, 102, 0, 0, 0, 0, 0, 0.0d0, 0.0d0)
10572 msg =
' It will not be issued again for this problem.' 10573 CALL xerrwd (msg, 50, 102, 0, 1, mxhnil, 0, 0, 0.0d0, 0.0d0)
10578 CALL dstoka (neq, y, rwork(lyh), nyh, rwork(lyh), rwork(lewt),
10579 1 rwork(lsavf), rwork(lsavx), rwork(lacor), rwork(lwm),
10580 2 iwork(liwm), f, jac, psol)
10582 GO TO (300, 530, 540, 550), kgo
10592 IF (ngc .EQ. 0)
GO TO 315
10593 CALL drchek (3, g, neq, y, rwork(lyh), nyh,
10594 1 rwork(lg0), rwork(lg1), rwork(lgx), jroot, irt)
10595 IF (irt .NE. 1)
GO TO 315
10602 GO TO (310, 400, 330, 340, 350), itask
10604 310
IF ((tn - tout)*h .LT. 0.0d0)
GO TO 250
10605 CALL dintdy (tout, 0, rwork(lyh), nyh, y, iflag)
10609 330
IF ((tn - tout)*h .GE. 0.0d0)
GO TO 400
10612 340
IF ((tn - tout)*h .LT. 0.0d0)
GO TO 345
10613 CALL dintdy (tout, 0, rwork(lyh), nyh, y, iflag)
10616 345 hmx = abs(tn) + abs(h)
10617 ihit = abs(tn - tcrit) .LE. 100.0d0*uround*hmx
10618 IF (ihit)
GO TO 400
10619 tnext = tn + h*(1.0d0 + 4.0d0*uround)
10620 IF ((tnext - tcrit)*h .LE. 0.0d0)
GO TO 250
10621 h = (tcrit - tn)*(1.0d0 - 4.0d0*uround)
10625 350 hmx = abs(tn) + abs(h)
10626 ihit = abs(tn - tcrit) .LE. 100.0d0*uround*hmx
10635 410 y(i) = rwork(i+lyh-1)
10637 IF (itask .NE. 4 .AND. itask .NE. 5)
GO TO 420
10638 IF (ihit) t = tcrit
10668 500 msg =
'DLSODKR- At current T (=R1), MXSTEP (=I1) steps ' 10669 CALL xerrwd (msg, 50, 201, 0, 0, 0, 0, 0, 0.0d0, 0.0d0)
10670 msg =
' taken on this call before reaching TOUT ' 10671 CALL xerrwd (msg, 50, 201, 0, 1, mxstep, 0, 1, tn, 0.0d0)
10675 510 ewti = rwork(lewt+i-1)
10676 msg = .le.
'DLSODKR- At T(=R1), EWT(I1) has become R2 0.' 10677 CALL xerrwd (msg, 50, 202, 0, 1, i, 0, 2, tn, ewti)
10681 520 msg =
'DLSODKR- At T (=R1), too much accuracy requested ' 10682 CALL xerrwd (msg, 50, 203, 0, 0, 0, 0, 0, 0.0d0, 0.0d0)
10683 msg =
' for precision of machine.. See TOLSF (=R2) ' 10684 CALL xerrwd (msg, 50, 203, 0, 0, 0, 0, 2, tn, tolsf)
10689 530 msg =
'DLSODKR- At T(=R1) and step size H(=R2), the error' 10690 CALL xerrwd (msg, 50, 204, 0, 0, 0, 0, 0, 0.0d0, 0.0d0)
10691 msg =
' test failed repeatedly or with ABS(H) = HMIN' 10692 CALL xerrwd (msg, 50, 204, 0, 0, 0, 0, 2, tn, h)
10696 540 msg =
'DLSODKR- At T (=R1) and step size H (=R2), the ' 10697 CALL xerrwd (msg, 50, 205, 0, 0, 0, 0, 0, 0.0d0, 0.0d0)
10698 msg =
' corrector convergence failed repeatedly ' 10699 CALL xerrwd (msg, 50, 205, 0, 0, 0, 0, 0, 0.0d0, 0.0d0)
10700 msg =
' or with ABS(H) = HMIN ' 10701 CALL xerrwd (msg, 30, 205, 0, 0, 0, 0, 2, tn, h)
10705 550 msg =
'DLSODKR- At T (=R1) an unrecoverable error return' 10706 CALL xerrwd (msg, 50, 206, 0, 0, 0, 0, 0, 0.0d0, 0.0d0)
10707 msg =
' was made from Subroutine PSOL ' 10708 CALL xerrwd (msg, 40, 206, 0, 0, 0, 0, 1, tn, 0.0d0)
10715 SIZE = abs(rwork(i+lacor-1)*rwork(i+lewt-1))
10716 IF (big .GE. size)
GO TO 570
10723 590 y(i) = rwork(i+lyh-1)
10750 601 msg =
'DLSODKR- ISTATE(=I1) illegal.' 10751 CALL xerrwd (msg, 30, 1, 0, 1, istate, 0, 0, 0.0d0, 0.0d0)
10752 IF (istate .LT. 0)
GO TO 800
10754 602 msg =
'DLSODKR- ITASK (=I1) illegal.' 10755 CALL xerrwd (msg, 30, 2, 0, 1, itask, 0, 0, 0.0d0, 0.0d0)
10757 603 msg = .gt.
'DLSODKR- ISTATE1 but DLSODKR not initialized. ' 10758 CALL xerrwd (msg, 50, 3, 0, 0, 0, 0, 0, 0.0d0, 0.0d0)
10760 604 msg = .lt.
'DLSODKR- NEQ (=I1) 1 ' 10761 CALL xerrwd (msg, 30, 4, 0, 1, neq(1), 0, 0, 0.0d0, 0.0d0)
10763 605 msg =
'DLSODKR- ISTATE = 3 and NEQ increased (I1 to I2).' 10764 CALL xerrwd (msg, 50, 5, 0, 2, n, neq(1), 0, 0.0d0, 0.0d0)
10766 606 msg =
'DLSODKR- ITOL (=I1) illegal. ' 10767 CALL xerrwd (msg, 30, 6, 0, 1, itol, 0, 0, 0.0d0, 0.0d0)
10769 607 msg =
'DLSODKR- IOPT (=I1) illegal. ' 10770 CALL xerrwd (msg, 30, 7, 0, 1, iopt, 0, 0, 0.0d0, 0.0d0)
10772 608 msg =
'DLSODKR- MF (=I1) illegal. ' 10773 CALL xerrwd (msg, 30, 8, 0, 1, mf, 0, 0, 0.0d0, 0.0d0)
10775 611 msg = .lt.
'DLSODKR- MAXORD (=I1) 0 ' 10776 CALL xerrwd (msg, 30, 11, 0, 1, maxord, 0, 0, 0.0d0, 0.0d0)
10778 612 msg = .lt.
'DLSODKR- MXSTEP (=I1) 0 ' 10779 CALL xerrwd (msg, 30, 12, 0, 1, mxstep, 0, 0, 0.0d0, 0.0d0)
10781 613 msg = .lt.
'DLSODKR- MXHNIL (=I1) 0 ' 10782 CALL xerrwd (msg, 30, 13, 0, 1, mxhnil, 0, 0, 0.0d0, 0.0d0)
10784 614 msg =
'DLSODKR- TOUT (=R1) behind T (=R2) ' 10785 CALL xerrwd (msg, 40, 14, 0, 0, 0, 0, 2, tout, t)
10786 msg =
' Integration direction is given by H0 (=R1) ' 10787 CALL xerrwd (msg, 50, 14, 0, 0, 0, 0, 1, h0, 0.0d0)
10789 615 msg = .lt.
'DLSODKR- HMAX (=R1) 0.0 ' 10790 CALL xerrwd (msg, 30, 15, 0, 0, 0, 0, 1, hmax, 0.0d0)
10792 616 msg = .lt.
'DLSODKR- HMIN (=R1) 0.0 ' 10793 CALL xerrwd (msg, 30, 16, 0, 0, 0, 0, 1, hmin, 0.0d0)
10795 617 msg=
'DLSODKR- RWORK length needed, LENRW(=I1), exceeds LRW(=I2) ' 10796 CALL xerrwd (msg, 60, 17, 0, 2, lenrw, lrw, 0, 0.0d0, 0.0d0)
10798 618 msg=
'DLSODKR- IWORK length needed, LENIW(=I1), exceeds LIW(=I2) ' 10799 CALL xerrwd (msg, 60, 18, 0, 2, leniw, liw, 0, 0.0d0, 0.0d0)
10801 619 msg = .lt.
'DLSODKR- RTOL(I1) is R1 0.0 ' 10802 CALL xerrwd (msg, 40, 19, 0, 1, i, 0, 1, rtoli, 0.0d0)
10804 620 msg = .lt.
'DLSODKR- ATOL(I1) is R1 0.0 ' 10805 CALL xerrwd (msg, 40, 20, 0, 1, i, 0, 1, atoli, 0.0d0)
10807 621 ewti = rwork(lewt+i-1)
10808 msg = .le.
'DLSODKR- EWT(I1) is R1 0.0 ' 10809 CALL xerrwd (msg, 40, 21, 0, 1, i, 0, 1, ewti, 0.0d0)
10811 622 msg=
'DLSODKR- TOUT(=R1) too close to T(=R2) to start integration.' 10812 CALL xerrwd (msg, 60, 22, 0, 0, 0, 0, 2, tout, t)
10814 623 msg=
'DLSODKR- ITASK = I1 and TOUT (=R1) behind TCUR - HU (= R2) ' 10815 CALL xerrwd (msg, 60, 23, 0, 1, itask, 0, 2, tout, tp)
10817 624 msg=
'DLSODKR- ITASK = 4 or 5 and TCRIT (=R1) behind TCUR (=R2) ' 10818 CALL xerrwd (msg, 60, 24, 0, 0, 0, 0, 2, tcrit, tn)
10820 625 msg=
'DLSODKR- ITASK = 4 or 5 and TCRIT (=R1) behind TOUT (=R2) ' 10821 CALL xerrwd (msg, 60, 25, 0, 0, 0, 0, 2, tcrit, tout)
10823 626 msg =
'DLSODKR- At start of problem, too much accuracy ' 10824 CALL xerrwd (msg, 50, 26, 0, 0, 0, 0, 0, 0.0d0, 0.0d0)
10825 msg=
' requested for precision of machine.. See TOLSF (=R1) ' 10826 CALL xerrwd (msg, 60, 26, 0, 0, 0, 0, 1, tolsf, 0.0d0)
10829 627 msg =
'DLSODKR- Trouble in DINTDY. ITASK = I1, TOUT = R1' 10830 CALL xerrwd (msg, 50, 27, 0, 1, itask, 0, 1, tout, 0.0d0)
10832 630 msg = .lt.
'DLSODKR- NG (=I1) 0 ' 10833 CALL xerrwd (msg, 30, 30, 0, 1, ng, 0, 0, 0.0d0, 0.0d0)
10835 631 msg =
'DLSODKR- NG changed (from I1 to I2) illegally, ' 10836 CALL xerrwd (msg, 50, 31, 0, 0, 0, 0, 0, 0.0d0, 0.0d0)
10837 msg =
' i.e. not immediately after a root was found.' 10838 CALL xerrwd (msg, 50, 31, 0, 2, ngc, ng, 0, 0.0d0, 0.0d0)
10840 632 msg =
'DLSODKR- One or more components of g has a root ' 10841 CALL xerrwd (msg, 50, 32, 0, 0, 0, 0, 0, 0.0d0, 0.0d0)
10842 msg =
' too near to the initial point. ' 10843 CALL xerrwd (msg, 40, 32, 0, 0, 0, 0, 0, 0.0d0, 0.0d0)
10848 800 msg =
'DLSODKR- Run aborted.. apparent infinite loop. ' 10849 CALL xerrwd (msg, 50, 303, 2, 0, 0, 0, 0, 0.0d0, 0.0d0)
10854 SUBROUTINE dlsodi (RES, ADDA, JAC, NEQ, Y, YDOTI, T, TOUT, ITOL,
10855 1 RTOL, ATOL, ITASK, ISTATE, IOPT, RWORK, LRW, IWORK, LIW, MF )
10856 EXTERNAL res, adda, jac
10857 INTEGER neq, itol, itask, istate, iopt, lrw, iwork, liw, mf
10858 DOUBLE PRECISION y, ydoti, t, tout, rtol, atol, rwork
10859 dimension neq(*), y(*), ydoti(*), rtol(*), atol(*), rwork(lrw),
12025 EXTERNAL dprepji, dsolsy
12026 DOUBLE PRECISION dumach, dvnorm
12027 INTEGER init, mxstep, mxhnil, nhnil, nslast, nyh, iowns,
12028 1 icf, ierpj, iersl, jcur, jstart, kflag, l,
12029 2 lyh, lewt, lacor, lsavf, lwm, liwm, meth, miter,
12030 3 maxord, maxcor, msbp, mxncf, n, nq, nst, nfe, nje, nqu
12031 INTEGER i, i1, i2, ier, iflag, imxer, ires, kgo,
12032 1 leniw, lenrw, lenwm, lp, lyd0, ml, mord, mu, mxhnl0, mxstp0
12033 DOUBLE PRECISION rowns,
12034 1 ccmax, el0, h, hmin, hmxi, hu, rc, tn, uround
12035 DOUBLE PRECISION atoli, ayi, big, ewti, h0, hmax, hmx, rh, rtoli,
12036 1 tcrit, tdist, tnext, tol, tolsf, tp,
SIZE, sum, w0
12040 SAVE mord, mxstp0, mxhnl0
12051 COMMON /dls001/ rowns(209),
12052 1 ccmax, el0, h, hmin, hmxi, hu, rc, tn, uround,
12053 2 init, mxstep, mxhnil, nhnil, nslast, nyh, iowns(6),
12054 3 icf, ierpj, iersl, jcur, jstart, kflag, l,
12055 4 lyh, lewt, lacor, lsavf, lwm, liwm, meth, miter,
12056 5 maxord, maxcor, msbp, mxncf, n, nq, nst, nfe, nje, nqu
12058 DATA mord(1),mord(2)/12,5/, mxstp0/500/, mxhnl0/10/
12067 IF (istate .LT. 0 .OR. istate .GT. 3)
GO TO 601
12068 IF (itask .LT. 1 .OR. itask .GT. 5)
GO TO 602
12069 IF (istate .LE. 1)
GO TO 10
12070 IF (init .EQ. 0)
GO TO 603
12071 IF (istate .EQ. 2)
GO TO 200
12074 IF (tout .EQ. t)
RETURN 12084 20
IF (neq(1) .LE. 0)
GO TO 604
12085 IF (istate .LE. 1)
GO TO 25
12086 IF (neq(1) .GT. n)
GO TO 605
12088 IF (itol .LT. 1 .OR. itol .GT. 4)
GO TO 606
12089 IF (iopt .LT. 0 .OR. iopt .GT. 1)
GO TO 607
12091 miter = mf - 10*meth
12092 IF (meth .LT. 1 .OR. meth .GT. 2)
GO TO 608
12093 IF (miter .LE. 0 .OR. miter .GT. 5)
GO TO 608
12094 IF (miter .EQ. 3)
GO TO 608
12095 IF (miter .LT. 3)
GO TO 30
12098 IF (ml .LT. 0 .OR. ml .GE. n)
GO TO 609
12099 IF (mu .LT. 0 .OR. mu .GE. n)
GO TO 610
12102 IF (iopt .EQ. 1)
GO TO 40
12103 maxord = mord(meth)
12106 IF (istate .LE. 1) h0 = 0.0d0
12110 40 maxord = iwork(5)
12111 IF (maxord .LT. 0)
GO TO 611
12112 IF (maxord .EQ. 0) maxord = 100
12113 maxord = min(maxord,mord(meth))
12115 IF (mxstep .LT. 0)
GO TO 612
12116 IF (mxstep .EQ. 0) mxstep = mxstp0
12118 IF (mxhnil .LT. 0)
GO TO 613
12119 IF (mxhnil .EQ. 0) mxhnil = mxhnl0
12120 IF (istate .GT. 1)
GO TO 50
12122 IF ((tout - t)*h0 .LT. 0.0d0)
GO TO 614
12124 IF (hmax .LT. 0.0d0)
GO TO 615
12126 IF (hmax .GT. 0.0d0) hmxi = 1.0d0/hmax
12128 IF (hmin .LT. 0.0d0)
GO TO 616
12136 IF (istate .LE. 1) nyh = n
12137 lwm = lyh + (maxord + 1)*nyh
12138 IF (miter .LE. 2) lenwm = n*n + 2
12139 IF (miter .GE. 4) lenwm = (2*ml + mu + 1)*n + 2
12143 lenrw = lacor + n - 1
12148 IF (lenrw .GT. lrw)
GO TO 617
12149 IF (leniw .GT. liw)
GO TO 618
12154 IF (itol .GE. 3) rtoli = rtol(i)
12155 IF (itol .EQ. 2 .OR. itol .EQ. 4) atoli = atol(i)
12156 IF (rtoli .LT. 0.0d0)
GO TO 619
12157 IF (atoli .LT. 0.0d0)
GO TO 620
12159 IF (istate .LE. 1)
GO TO 100
12162 IF (nq .LE. maxord)
GO TO 90
12165 80 ydoti(i) = rwork(i+lwm-1)
12167 90 rwork(lwm) = sqrt(uround)
12168 IF (n .EQ. nyh)
GO TO 200
12171 i2 = lyh + (maxord + 1)*nyh - 1
12172 IF (i1 .GT. i2)
GO TO 200
12174 95 rwork(i) = 0.0d0
12183 100 uround = dumach()
12185 IF (itask .NE. 4 .AND. itask .NE. 5)
GO TO 105
12187 IF ((tcrit - tout)*(tout - t) .LT. 0.0d0)
GO TO 625
12188 IF (h0 .NE. 0.0d0 .AND. (t + h0 - tcrit)*h0 .GT. 0.0d0)
12191 rwork(lwm) = sqrt(uround)
12206 IF (istate .EQ. 1)
GO TO 120
12208 CALL dainvg( res, adda, neq, t, y, rwork(lyd0), miter,
12209 1 ml, mu, rwork(lp), iwork(21), ier )
12211 IF (ier .LT. 0)
GO TO 560
12212 IF (ier .GT. 0)
GO TO 565
12214 115 rwork(i+lyh-1) = y(i)
12218 rwork(i+lyh-1) = y(i)
12219 125 rwork(i+lyd0-1) = ydoti(i)
12224 CALL dewset (n, itol, rtol, atol, rwork(lyh), rwork(lewt))
12226 IF (rwork(i+lewt-1) .LE. 0.0d0)
GO TO 621
12227 135 rwork(i+lewt-1) = 1.0d0/rwork(i+lewt-1)
12244 IF (h0 .NE. 0.0d0)
GO TO 180
12245 tdist = abs(tout - t)
12246 w0 = max(abs(t),abs(tout))
12247 IF (tdist .LT. 2.0d0*uround*w0)
GO TO 622
12249 IF (itol .LE. 2)
GO TO 145
12251 140 tol = max(tol,rtol(i))
12252 145
IF (tol .GT. 0.0d0)
GO TO 160
12255 IF (itol .EQ. 2 .OR. itol .EQ. 4) atoli = atol(i)
12257 IF (ayi .NE. 0.0d0) tol = max(tol,atoli/ayi)
12259 160 tol = max(tol,100.0d0*uround)
12260 tol = min(tol,0.001d0)
12261 sum = dvnorm(n, rwork(lyd0), rwork(lewt))
12262 sum = 1.0d0/(tol*w0*w0) + tol*sum**2
12263 h0 = 1.0d0/sqrt(sum)
12265 h0 = sign(h0,tout-t)
12267 180 rh = abs(h0)*hmxi
12268 IF (rh .GT. 1.0d0) h0 = h0/rh
12272 190 rwork(i+lyd0-1) = h0*rwork(i+lyd0-1)
12280 GO TO (210, 250, 220, 230, 240), itask
12281 210
IF ((tn - tout)*h .LT. 0.0d0)
GO TO 250
12282 CALL dintdy (tout, 0, rwork(lyh), nyh, y, iflag)
12283 IF (iflag .NE. 0)
GO TO 627
12286 220 tp = tn - hu*(1.0d0 + 100.0d0*uround)
12287 IF ((tp - tout)*h .GT. 0.0d0)
GO TO 623
12288 IF ((tn - tout)*h .LT. 0.0d0)
GO TO 250
12290 230 tcrit = rwork(1)
12291 IF ((tn - tcrit)*h .GT. 0.0d0)
GO TO 624
12292 IF ((tcrit - tout)*h .LT. 0.0d0)
GO TO 625
12293 IF ((tn - tout)*h .LT. 0.0d0)
GO TO 245
12294 CALL dintdy (tout, 0, rwork(lyh), nyh, y, iflag)
12295 IF (iflag .NE. 0)
GO TO 627
12298 240 tcrit = rwork(1)
12299 IF ((tn - tcrit)*h .GT. 0.0d0)
GO TO 624
12300 245 hmx = abs(tn) + abs(h)
12301 ihit = abs(tn - tcrit) .LE. 100.0d0*uround*hmx
12302 IF (ihit)
GO TO 400
12303 tnext = tn + h*(1.0d0 + 4.0d0*uround)
12304 IF ((tnext - tcrit)*h .LE. 0.0d0)
GO TO 250
12305 h = (tcrit - tn)*(1.0d0 - 4.0d0*uround)
12306 IF (istate .EQ. 2) jstart = -2
12319 IF ((nst-nslast) .GE. mxstep)
GO TO 500
12320 CALL dewset (n, itol, rtol, atol, rwork(lyh), rwork(lewt))
12322 IF (rwork(i+lewt-1) .LE. 0.0d0)
GO TO 510
12323 260 rwork(i+lewt-1) = 1.0d0/rwork(i+lewt-1)
12324 270 tolsf = uround*dvnorm(n, rwork(lyh), rwork(lewt))
12325 IF (tolsf .LE. 1.0d0)
GO TO 280
12326 tolsf = tolsf*2.0d0
12327 IF (nst .EQ. 0)
GO TO 626
12329 280
IF ((tn + h) .NE. tn)
GO TO 290
12331 IF (nhnil .GT. mxhnil)
GO TO 290
12332 msg =
'DLSODI- Warning..Internal T (=R1) and H (=R2) are' 12333 CALL xerrwd (msg, 50, 101, 0, 0, 0, 0, 0, 0.0d0, 0.0d0)
12334 msg=
' such that in the machine, T + H = T on the next step ' 12335 CALL xerrwd (msg, 60, 101, 0, 0, 0, 0, 0, 0.0d0, 0.0d0)
12336 msg =
' (H = step size). Solver will continue anyway.' 12337 CALL xerrwd (msg, 50, 101, 0, 0, 0, 0, 2, tn, h)
12338 IF (nhnil .LT. mxhnil)
GO TO 290
12339 msg =
'DLSODI- Above warning has been issued I1 times. ' 12340 CALL xerrwd (msg, 50, 102, 0, 0, 0, 0, 0, 0.0d0, 0.0d0)
12341 msg =
' It will not be issued again for this problem.' 12342 CALL xerrwd (msg, 50, 102, 0, 1, mxhnil, 0, 0, 0.0d0, 0.0d0)
12349 CALL dstodi (neq, y, rwork(lyh), nyh, rwork(lyh), rwork(lewt),
12350 1 ydoti, rwork(lsavf), rwork(lacor), rwork(lwm),
12351 2 iwork(liwm), res, adda, jac, dprepji, dsolsy )
12353 GO TO (300, 530, 540, 400, 550), kgo
12363 GO TO (310, 400, 330, 340, 350), itask
12365 310
IF ((tn - tout)*h .LT. 0.0d0)
GO TO 250
12366 CALL dintdy (tout, 0, rwork(lyh), nyh, y, iflag)
12370 330
IF ((tn - tout)*h .GE. 0.0d0)
GO TO 400
12373 340
IF ((tn - tout)*h .LT. 0.0d0)
GO TO 345
12374 CALL dintdy (tout, 0, rwork(lyh), nyh, y, iflag)
12377 345 hmx = abs(tn) + abs(h)
12378 ihit = abs(tn - tcrit) .LE. 100.0d0*uround*hmx
12379 IF (ihit)
GO TO 400
12380 tnext = tn + h*(1.0d0 + 4.0d0*uround)
12381 IF ((tnext - tcrit)*h .LE. 0.0d0)
GO TO 250
12382 h = (tcrit - tn)*(1.0d0 - 4.0d0*uround)
12386 350 hmx = abs(tn) + abs(h)
12387 ihit = abs(tn - tcrit) .LE. 100.0d0*uround*hmx
12396 410 y(i) = rwork(i+lyh-1)
12398 IF (itask .NE. 4 .AND. itask .NE. 5)
GO TO 420
12399 IF (ihit) t = tcrit
12401 IF (kflag .EQ. -3) istate = 3
12420 500 msg =
'DLSODI- At current T (=R1), MXSTEP (=I1) steps ' 12421 CALL xerrwd (msg, 50, 201, 0, 0, 0, 0, 0, 0.0d0, 0.0d0)
12422 msg =
' taken on this call before reaching TOUT ' 12423 CALL xerrwd (msg, 50, 201, 0, 1, mxstep, 0, 1, tn, 0.0d0)
12427 510 ewti = rwork(lewt+i-1)
12428 msg = .le.
'DLSODI- At T (=R1), EWT(I1) has become R2 0.' 12429 CALL xerrwd (msg, 50, 202, 0, 1, i, 0, 2, tn, ewti)
12433 520 msg =
'DLSODI- At T (=R1), too much accuracy requested ' 12434 CALL xerrwd (msg, 50, 203, 0, 0, 0, 0, 0, 0.0d0, 0.0d0)
12435 msg =
' for precision of machine.. See TOLSF (=R2) ' 12436 CALL xerrwd (msg, 50, 203, 0, 0, 0, 0, 2, tn, tolsf)
12441 530 msg =
'DLSODI- At T(=R1) and step size H(=R2), the error' 12442 CALL xerrwd (msg, 50, 204, 0, 0, 0, 0, 0, 0.0d0, 0.0d0)
12443 msg =
' test failed repeatedly or with ABS(H) = HMIN' 12444 CALL xerrwd (msg, 50, 204, 0, 0, 0, 0, 2, tn, h)
12448 540 msg =
'DLSODI- At T (=R1) and step size H (=R2), the ' 12449 CALL xerrwd (msg, 50, 205, 0, 0, 0, 0, 0, 0.0d0, 0.0d0)
12450 msg =
' corrector convergence failed repeatedly ' 12451 CALL xerrwd (msg, 50, 205, 0, 0, 0, 0, 0, 0.0d0, 0.0d0)
12452 msg =
' or with ABS(H) = HMIN ' 12453 CALL xerrwd (msg, 30, 205, 0, 0, 0, 0, 2, tn, h)
12457 550 msg =
'DLSODI- At T (=R1) residual routine returned ' 12458 CALL xerrwd (msg, 50, 206, 0, 0, 0, 0, 0, 0.0d0, 0.0d0)
12459 msg =
' error IRES = 3 repeatedly. ' 12460 CALL xerrwd (msg, 40, 206, 0, 0, 0, 0, 1, tn, 0.0d0)
12465 msg=
'DLSODI- Attempt to initialize dy/dt failed: Matrix A is ' 12466 CALL xerrwd (msg, 60, 207, 0, 0, 0, 0, 0, 0.0d0, 0.0d0)
12467 msg =
' singular. DGEFA or DGBFA returned INFO = I1' 12468 CALL xerrwd (msg, 50, 207, 0, 1, ier, 0, 0, 0.0d0, 0.0d0)
12472 565 msg =
'DLSODI- Attempt to initialize dy/dt failed ' 12473 CALL xerrwd (msg, 50, 208, 0, 0, 0, 0, 0, 0.0d0, 0.0d0)
12474 msg =
' because residual routine set its error flag ' 12475 CALL xerrwd (msg, 50, 208, 0, 0, 0, 0, 0, 0.0d0, 0.0d0)
12476 msg =
' to IRES = (I1)' 12477 CALL xerrwd (msg, 20, 208, 0, 1, ier, 0, 0, 0.0d0, 0.0d0)
12484 SIZE = abs(rwork(i+lacor-1)*rwork(i+lewt-1))
12485 IF (big .GE. size)
GO TO 575
12491 580 lyd0 = lyh + nyh
12493 rwork(i+lsavf-1) = rwork(i+lyd0-1)/h
12494 585 y(i) = rwork(i+lyh-1)
12496 CALL res (neq, tn, y, rwork(lsavf), ydoti, ires )
12498 IF (ires .LE. 1)
GO TO 595
12499 msg =
'DLSODI- Residual routine set its flag IRES ' 12500 CALL xerrwd (msg, 50, 210, 0, 0, 0, 0, 0, 0.0d0, 0.0d0)
12501 msg =
' to (I1) when called for final output. ' 12502 CALL xerrwd (msg, 50, 210, 0, 1, ires, 0, 0, 0.0d0, 0.0d0)
12506 592 y(i) = rwork(i+lyh-1)
12524 601 msg =
'DLSODI- ISTATE (=I1) illegal.' 12525 CALL xerrwd (msg, 30, 1, 0, 1, istate, 0, 0, 0.0d0, 0.0d0)
12526 IF (istate .LT. 0)
GO TO 800
12528 602 msg =
'DLSODI- ITASK (=I1) illegal. ' 12529 CALL xerrwd (msg, 30, 2, 0, 1, itask, 0, 0, 0.0d0, 0.0d0)
12531 603 msg = .gt.
'DLSODI- ISTATE 1 but DLSODI not initialized.' 12532 CALL xerrwd (msg, 50, 3, 0, 0, 0, 0, 0, 0.0d0, 0.0d0)
12534 604 msg = .lt.
'DLSODI- NEQ (=I1) 1 ' 12535 CALL xerrwd (msg, 30, 4, 0, 1, neq(1), 0, 0, 0.0d0, 0.0d0)
12537 605 msg =
'DLSODI- ISTATE = 3 and NEQ increased (I1 to I2). ' 12538 CALL xerrwd (msg, 50, 5, 0, 2, n, neq(1), 0, 0.0d0, 0.0d0)
12540 606 msg =
'DLSODI- ITOL (=I1) illegal. ' 12541 CALL xerrwd (msg, 30, 6, 0, 1, itol, 0, 0, 0.0d0, 0.0d0)
12543 607 msg =
'DLSODI- IOPT (=I1) illegal. ' 12544 CALL xerrwd (msg, 30, 7, 0, 1, iopt, 0, 0, 0.0d0, 0.0d0)
12546 608 msg =
'DLSODI- MF (=I1) illegal. ' 12547 CALL xerrwd (msg, 30, 8, 0, 1, mf, 0, 0, 0.0d0, 0.0d0)
12549 609 msg = .lt..ge.
'DLSODI- ML(=I1) illegal: 0 or NEQ(=I2) ' 12550 CALL xerrwd (msg, 50, 9, 0, 2, ml, neq(1), 0, 0.0d0, 0.0d0)
12552 610 msg = .lt..ge.
'DLSODI- MU(=I1) illegal: 0 or NEQ(=I2) ' 12553 CALL xerrwd (msg, 50, 10, 0, 2, mu, neq(1), 0, 0.0d0, 0.0d0)
12555 611 msg = .lt.
'DLSODI- MAXORD (=I1) 0 ' 12556 CALL xerrwd (msg, 30, 11, 0, 1, maxord, 0, 0, 0.0d0, 0.0d0)
12558 612 msg = .lt.
'DLSODI- MXSTEP (=I1) 0 ' 12559 CALL xerrwd (msg, 30, 12, 0, 1, mxstep, 0, 0, 0.0d0, 0.0d0)
12561 613 msg = .lt.
'DLSODI- MXHNIL (=I1) 0 ' 12562 CALL xerrwd (msg, 30, 13, 0, 1, mxhnil, 0, 0, 0.0d0, 0.0d0)
12564 614 msg =
'DLSODI- TOUT (=R1) behind T (=R2) ' 12565 CALL xerrwd (msg, 40, 14, 0, 0, 0, 0, 2, tout, t)
12566 msg =
' Integration direction is given by H0 (=R1) ' 12567 CALL xerrwd (msg, 50, 14, 0, 0, 0, 0, 1, h0, 0.0d0)
12569 615 msg = .lt.
'DLSODI- HMAX (=R1) 0.0 ' 12570 CALL xerrwd (msg, 30, 15, 0, 0, 0, 0, 1, hmax, 0.0d0)
12572 616 msg = .lt.
'DLSODI- HMIN (=R1) 0.0 ' 12573 CALL xerrwd (msg, 30, 16, 0, 0, 0, 0, 1, hmin, 0.0d0)
12575 617 msg=
'DLSODI- RWORK length needed, LENRW (=I1), exceeds LRW (=I2)' 12576 CALL xerrwd (msg, 60, 17, 0, 2, lenrw, lrw, 0, 0.0d0, 0.0d0)
12578 618 msg=
'DLSODI- IWORK length needed, LENIW (=I1), exceeds LIW (=I2)' 12579 CALL xerrwd (msg, 60, 18, 0, 2, leniw, liw, 0, 0.0d0, 0.0d0)
12581 619 msg = .lt.
'DLSODI- RTOL(=I1) is R1 0.0 ' 12582 CALL xerrwd (msg, 40, 19, 0, 1, i, 0, 1, rtoli, 0.0d0)
12584 620 msg = .lt.
'DLSODI- ATOL(=I1) is R1 0.0 ' 12585 CALL xerrwd (msg, 40, 20, 0, 1, i, 0, 1, atoli, 0.0d0)
12587 621 ewti = rwork(lewt+i-1)
12588 msg = .le.
'DLSODI- EWT(I1) is R1 0.0 ' 12589 CALL xerrwd (msg, 40, 21, 0, 1, i, 0, 1, ewti, 0.0d0)
12591 622 msg=
'DLSODI- TOUT(=R1) too close to T(=R2) to start integration.' 12592 CALL xerrwd (msg, 60, 22, 0, 0, 0, 0, 2, tout, t)
12594 623 msg=
'DLSODI- ITASK = I1 and TOUT (=R1) behind TCUR - HU (= R2) ' 12595 CALL xerrwd (msg, 60, 23, 0, 1, itask, 0, 2, tout, tp)
12597 624 msg=
'DLSODI- ITASK = 4 or 5 and TCRIT (=R1) behind TCUR (=R2) ' 12598 CALL xerrwd (msg, 60, 24, 0, 0, 0, 0, 2, tcrit, tn)
12600 625 msg=
'DLSODI- ITASK = 4 or 5 and TCRIT (=R1) behind TOUT (=R2) ' 12601 CALL xerrwd (msg, 60, 25, 0, 0, 0, 0, 2, tcrit, tout)
12603 626 msg =
'DLSODI- At start of problem, too much accuracy ' 12604 CALL xerrwd (msg, 50, 26, 0, 0, 0, 0, 0, 0.0d0, 0.0d0)
12605 msg=
' requested for precision of machine.. See TOLSF (=R1) ' 12606 CALL xerrwd (msg, 60, 26, 0, 0, 0, 0, 1, tolsf, 0.0d0)
12609 627 msg =
'DLSODI- Trouble in DINTDY. ITASK = I1, TOUT = R1' 12610 CALL xerrwd (msg, 50, 27, 0, 1, itask, 0, 1, tout, 0.0d0)
12615 800 msg =
'DLSODI- Run aborted.. apparent infinite loop. ' 12616 CALL xerrwd (msg, 50, 303, 2, 0, 0, 0, 0, 0.0d0, 0.0d0)
12621 SUBROUTINE dlsoibt (RES, ADDA, JAC, NEQ, Y, YDOTI, T, TOUT, ITOL,
12622 1 RTOL, ATOL, ITASK, ISTATE, IOPT, RWORK, LRW, IWORK, LIW, MF )
12623 EXTERNAL res, adda, jac
12624 INTEGER neq, itol, itask, istate, iopt, lrw, iwork, liw, mf
12625 DOUBLE PRECISION y, ydoti, t, tout, rtol, atol, rwork
12626 dimension neq(*), y(*), ydoti(*), rtol(*), atol(*), rwork(lrw),
13827 EXTERNAL dpjibt, dslsbt
13828 DOUBLE PRECISION dumach, dvnorm
13829 INTEGER init, mxstep, mxhnil, nhnil, nslast, nyh, iowns,
13830 1 icf, ierpj, iersl, jcur, jstart, kflag, l,
13831 2 lyh, lewt, lacor, lsavf, lwm, liwm, meth, miter,
13832 3 maxord, maxcor, msbp, mxncf, n, nq, nst, nfe, nje, nqu
13833 INTEGER i, i1, i2, ier, iflag, imxer, ires, kgo,
13834 1 leniw, lenrw, lenwm, lp, lyd0, mb, mord, mxhnl0, mxstp0, nb
13835 DOUBLE PRECISION rowns,
13836 1 ccmax, el0, h, hmin, hmxi, hu, rc, tn, uround
13837 DOUBLE PRECISION atoli, ayi, big, ewti, h0, hmax, hmx, rh, rtoli,
13838 1 tcrit, tdist, tnext, tol, tolsf, tp,
SIZE, sum, w0
13842 SAVE mord, mxstp0, mxhnl0
13853 COMMON /dls001/ rowns(209),
13854 1 ccmax, el0, h, hmin, hmxi, hu, rc, tn, uround,
13855 2 init, mxstep, mxhnil, nhnil, nslast, nyh, iowns(6),
13856 3 icf, ierpj, iersl, jcur, jstart, kflag, l,
13857 4 lyh, lewt, lacor, lsavf, lwm, liwm, meth, miter,
13858 5 maxord, maxcor, msbp, mxncf, n, nq, nst, nfe, nje, nqu
13860 DATA mord(1),mord(2)/12,5/, mxstp0/500/, mxhnl0/10/
13869 IF (istate .LT. 0 .OR. istate .GT. 3)
GO TO 601
13870 IF (itask .LT. 1 .OR. itask .GT. 5)
GO TO 602
13871 IF (istate .LE. 1)
GO TO 10
13872 IF (init .EQ. 0)
GO TO 603
13873 IF (istate .EQ. 2)
GO TO 200
13876 IF (tout .EQ. t)
RETURN 13886 20
IF (neq(1) .LE. 0)
GO TO 604
13887 IF (istate .LE. 1)
GO TO 25
13888 IF (neq(1) .GT. n)
GO TO 605
13890 IF (itol .LT. 1 .OR. itol .GT. 4)
GO TO 606
13891 IF (iopt .LT. 0 .OR. iopt .GT. 1)
GO TO 607
13893 miter = mf - 10*meth
13894 IF (meth .LT. 1 .OR. meth .GT. 2)
GO TO 608
13895 IF (miter .LT. 1 .OR. miter .GT. 2)
GO TO 608
13898 IF (mb .LT. 1 .OR. mb .GT. n)
GO TO 609
13899 IF (nb .LT. 4)
GO TO 610
13900 IF (mb*nb .NE. n)
GO TO 609
13902 IF (iopt .EQ. 1)
GO TO 40
13903 maxord = mord(meth)
13906 IF (istate .LE. 1) h0 = 0.0d0
13910 40 maxord = iwork(5)
13911 IF (maxord .LT. 0)
GO TO 611
13912 IF (maxord .EQ. 0) maxord = 100
13913 maxord = min(maxord,mord(meth))
13915 IF (mxstep .LT. 0)
GO TO 612
13916 IF (mxstep .EQ. 0) mxstep = mxstp0
13918 IF (mxhnil .LT. 0)
GO TO 613
13919 IF (mxhnil .EQ. 0) mxhnil = mxhnl0
13920 IF (istate .GT. 1)
GO TO 50
13922 IF ((tout - t)*h0 .LT. 0.0d0)
GO TO 614
13924 IF (hmax .LT. 0.0d0)
GO TO 615
13926 IF (hmax .GT. 0.0d0) hmxi = 1.0d0/hmax
13928 IF (hmin .LT. 0.0d0)
GO TO 616
13936 IF (istate .LE. 1) nyh = n
13937 lwm = lyh + (maxord + 1)*nyh
13938 lenwm = 3*mb*mb*nb + 2
13942 lenrw = lacor + n - 1
13947 IF (lenrw .GT. lrw)
GO TO 617
13948 IF (leniw .GT. liw)
GO TO 618
13953 IF (itol .GE. 3) rtoli = rtol(i)
13954 IF (itol .EQ. 2 .OR. itol .EQ. 4) atoli = atol(i)
13955 IF (rtoli .LT. 0.0d0)
GO TO 619
13956 IF (atoli .LT. 0.0d0)
GO TO 620
13958 IF (istate .LE. 1)
GO TO 100
13961 IF (nq .LE. maxord)
GO TO 90
13964 80 ydoti(i) = rwork(i+lwm-1)
13966 90 rwork(lwm) = sqrt(uround)
13967 IF (n .EQ. nyh)
GO TO 200
13970 i2 = lyh + (maxord + 1)*nyh - 1
13971 IF (i1 .GT. i2)
GO TO 200
13973 95 rwork(i) = 0.0d0
13982 100 uround = dumach()
13984 IF (itask .NE. 4 .AND. itask .NE. 5)
GO TO 105
13986 IF ((tcrit - tout)*(tout - t) .LT. 0.0d0)
GO TO 625
13987 IF (h0 .NE. 0.0d0 .AND. (t + h0 - tcrit)*h0 .GT. 0.0d0)
13990 rwork(lwm) = sqrt(uround)
14005 IF ( istate .EQ. 1 )
GO TO 120
14007 CALL daigbt( res, adda, neq, t, y, rwork(lyd0),
14008 1 mb, nb, rwork(lp), iwork(21), ier )
14010 IF (ier .LT. 0)
GO TO 560
14011 IF (ier .GT. 0)
GO TO 565
14013 115 rwork(i+lyh-1) = y(i)
14017 rwork(i+lyh-1) = y(i)
14018 125 rwork(i+lyd0-1) = ydoti(i)
14023 CALL dewset (n, itol, rtol, atol, rwork(lyh), rwork(lewt))
14025 IF (rwork(i+lewt-1) .LE. 0.0d0)
GO TO 621
14026 135 rwork(i+lewt-1) = 1.0d0/rwork(i+lewt-1)
14043 IF (h0 .NE. 0.0d0)
GO TO 180
14044 tdist = abs(tout - t)
14045 w0 = max(abs(t),abs(tout))
14046 IF (tdist .LT. 2.0d0*uround*w0)
GO TO 622
14048 IF (itol .LE. 2)
GO TO 145
14050 140 tol = max(tol,rtol(i))
14051 145
IF (tol .GT. 0.0d0)
GO TO 160
14054 IF (itol .EQ. 2 .OR. itol .EQ. 4) atoli = atol(i)
14056 IF (ayi .NE. 0.0d0) tol = max(tol,atoli/ayi)
14058 160 tol = max(tol,100.0d0*uround)
14059 tol = min(tol,0.001d0)
14060 sum = dvnorm(n, rwork(lyd0), rwork(lewt))
14061 sum = 1.0d0/(tol*w0*w0) + tol*sum**2
14062 h0 = 1.0d0/sqrt(sum)
14064 h0 = sign(h0,tout-t)
14066 180 rh = abs(h0)*hmxi
14067 IF (rh .GT. 1.0d0) h0 = h0/rh
14071 190 rwork(i+lyd0-1) = h0*rwork(i+lyd0-1)
14079 GO TO (210, 250, 220, 230, 240), itask
14080 210
IF ((tn - tout)*h .LT. 0.0d0)
GO TO 250
14081 CALL dintdy (tout, 0, rwork(lyh), nyh, y, iflag)
14082 IF (iflag .NE. 0)
GO TO 627
14085 220 tp = tn - hu*(1.0d0 + 100.0d0*uround)
14086 IF ((tp - tout)*h .GT. 0.0d0)
GO TO 623
14087 IF ((tn - tout)*h .LT. 0.0d0)
GO TO 250
14089 230 tcrit = rwork(1)
14090 IF ((tn - tcrit)*h .GT. 0.0d0)
GO TO 624
14091 IF ((tcrit - tout)*h .LT. 0.0d0)
GO TO 625
14092 IF ((tn - tout)*h .LT. 0.0d0)
GO TO 245
14093 CALL dintdy (tout, 0, rwork(lyh), nyh, y, iflag)
14094 IF (iflag .NE. 0)
GO TO 627
14097 240 tcrit = rwork(1)
14098 IF ((tn - tcrit)*h .GT. 0.0d0)
GO TO 624
14099 245 hmx = abs(tn) + abs(h)
14100 ihit = abs(tn - tcrit) .LE. 100.0d0*uround*hmx
14101 IF (ihit)
GO TO 400
14102 tnext = tn + h*(1.0d0 + 4.0d0*uround)
14103 IF ((tnext - tcrit)*h .LE. 0.0d0)
GO TO 250
14104 h = (tcrit - tn)*(1.0d0 - 4.0d0*uround)
14105 IF (istate .EQ. 2) jstart = -2
14118 IF ((nst-nslast) .GE. mxstep)
GO TO 500
14119 CALL dewset (n, itol, rtol, atol, rwork(lyh), rwork(lewt))
14121 IF (rwork(i+lewt-1) .LE. 0.0d0)
GO TO 510
14122 260 rwork(i+lewt-1) = 1.0d0/rwork(i+lewt-1)
14123 270 tolsf = uround*dvnorm(n, rwork(lyh), rwork(lewt))
14124 IF (tolsf .LE. 1.0d0)
GO TO 280
14125 tolsf = tolsf*2.0d0
14126 IF (nst .EQ. 0)
GO TO 626
14128 280
IF ((tn + h) .NE. tn)
GO TO 290
14130 IF (nhnil .GT. mxhnil)
GO TO 290
14131 msg =
'DLSOIBT- Warning..Internal T (=R1) and H (=R2) are' 14132 CALL xerrwd (msg, 50, 101, 0, 0, 0, 0, 0, 0.0d0, 0.0d0)
14133 msg=
' such that in the machine, T + H = T on the next step ' 14134 CALL xerrwd (msg, 60, 101, 0, 0, 0, 0, 0, 0.0d0, 0.0d0)
14135 msg =
' (H = step size). Solver will continue anyway.' 14136 CALL xerrwd (msg, 50, 101, 0, 0, 0, 0, 2, tn, h)
14137 IF (nhnil .LT. mxhnil)
GO TO 290
14138 msg =
'DLSOIBT- Above warning has been issued I1 times. ' 14139 CALL xerrwd (msg, 50, 102, 0, 0, 0, 0, 0, 0.0d0, 0.0d0)
14140 msg =
' It will not be issued again for this problem.' 14141 CALL xerrwd (msg, 50, 102, 0, 1, mxhnil, 0, 0, 0.0d0, 0.0d0)
14148 CALL dstodi (neq, y, rwork(lyh), nyh, rwork(lyh), rwork(lewt),
14149 1 ydoti, rwork(lsavf), rwork(lacor), rwork(lwm),
14150 2 iwork(liwm), res, adda, jac, dpjibt, dslsbt )
14152 GO TO (300, 530, 540, 400, 550), kgo
14162 GO TO (310, 400, 330, 340, 350), itask
14164 310
IF ((tn - tout)*h .LT. 0.0d0)
GO TO 250
14165 CALL dintdy (tout, 0, rwork(lyh), nyh, y, iflag)
14169 330
IF ((tn - tout)*h .GE. 0.0d0)
GO TO 400
14172 340
IF ((tn - tout)*h .LT. 0.0d0)
GO TO 345
14173 CALL dintdy (tout, 0, rwork(lyh), nyh, y, iflag)
14176 345 hmx = abs(tn) + abs(h)
14177 ihit = abs(tn - tcrit) .LE. 100.0d0*uround*hmx
14178 IF (ihit)
GO TO 400
14179 tnext = tn + h*(1.0d0 + 4.0d0*uround)
14180 IF ((tnext - tcrit)*h .LE. 0.0d0)
GO TO 250
14181 h = (tcrit - tn)*(1.0d0 - 4.0d0*uround)
14185 350 hmx = abs(tn) + abs(h)
14186 ihit = abs(tn - tcrit) .LE. 100.0d0*uround*hmx
14195 410 y(i) = rwork(i+lyh-1)
14197 IF (itask .NE. 4 .AND. itask .NE. 5)
GO TO 420
14198 IF (ihit) t = tcrit
14200 IF ( kflag .EQ. -3 ) istate = 3
14219 500 msg =
'DLSOIBT- At current T (=R1), MXSTEP (=I1) steps ' 14220 CALL xerrwd (msg, 50, 201, 0, 0, 0, 0, 0, 0.0d0, 0.0d0)
14221 msg =
' taken on this call before reaching TOUT ' 14222 CALL xerrwd (msg, 50, 201, 0, 1, mxstep, 0, 1, tn, 0.0d0)
14226 510 ewti = rwork(lewt+i-1)
14227 msg = .le.
'DLSOIBT- At T (=R1), EWT(I1) has become R2 0.' 14228 CALL xerrwd (msg, 50, 202, 0, 1, i, 0, 2, tn, ewti)
14232 520 msg =
'DLSOIBT- At T (=R1), too much accuracy requested ' 14233 CALL xerrwd (msg, 50, 203, 0, 0, 0, 0, 0, 0.0d0, 0.0d0)
14234 msg =
' for precision of machine.. See TOLSF (=R2) ' 14235 CALL xerrwd (msg, 50, 203, 0, 0, 0, 0, 2, tn, tolsf)
14240 530 msg =
'DLSOIBT- At T (=R1) and step size H (=R2), the ' 14241 CALL xerrwd (msg, 50, 204, 0, 0, 0, 0, 0, 0.0d0, 0.0d0)
14242 msg =
'error test failed repeatedly or with ABS(H) = HMIN' 14243 CALL xerrwd (msg, 50, 204, 0, 0, 0, 0, 2, tn, h)
14247 540 msg =
'DLSOIBT- At T (=R1) and step size H (=R2), the ' 14248 CALL xerrwd (msg, 50, 205, 0, 0, 0, 0, 0, 0.0d0, 0.0d0)
14249 msg =
' corrector convergence failed repeatedly ' 14250 CALL xerrwd (msg, 50, 205, 0, 0, 0, 0, 0, 0.0d0, 0.0d0)
14251 msg =
' or with ABS(H) = HMIN ' 14252 CALL xerrwd (msg, 30, 205, 0, 0, 0, 0, 2, tn, h)
14256 550 msg =
'DLSOIBT- At T (=R1) residual routine returned ' 14257 CALL xerrwd (msg, 50, 206, 0, 0, 0, 0, 0, 0.0d0, 0.0d0)
14258 msg =
' error IRES = 3 repeatedly. ' 14259 CALL xerrwd (msg, 40, 206, 0, 0, 0, 0, 1, tn, 0.0d0)
14264 msg=
'DLSOIBT- Attempt to initialize dy/dt failed: Matrix A has a' 14265 CALL xerrwd (msg, 60, 207, 0, 0, 0, 0, 0, 0.0d0, 0.0d0)
14266 msg =
' singular diagonal block, block no. = (I1) ' 14267 CALL xerrwd (msg, 50, 207, 0, 1, ier, 0, 0, 0.0d0, 0.0d0)
14271 565 msg =
'DLSOIBT- Attempt to initialize dy/dt failed ' 14272 CALL xerrwd (msg, 50, 208, 0, 0, 0, 0, 0, 0.0d0, 0.0d0)
14273 msg =
' because residual routine set its error flag ' 14274 CALL xerrwd (msg, 50, 208, 0, 0, 0, 0, 0, 0.0d0, 0.0d0)
14275 msg =
' to IRES = (I1)' 14276 CALL xerrwd (msg, 20, 208, 0, 1, ier, 0, 0, 0.0d0, 0.0d0)
14283 SIZE = abs(rwork(i+lacor-1)*rwork(i+lewt-1))
14284 IF (big .GE. size)
GO TO 575
14290 580 lyd0 = lyh + nyh
14292 rwork(i+lsavf-1) = rwork(i+lyd0-1)/h
14293 585 y(i) = rwork(i+lyh-1)
14295 CALL res (neq, tn, y, rwork(lsavf), ydoti, ires)
14297 IF (ires .LE. 1)
GO TO 595
14298 msg =
'DLSOIBT- Residual routine set its flag IRES ' 14299 CALL xerrwd (msg, 50, 210, 0, 0, 0, 0, 0, 0.0d0, 0.0d0)
14300 msg =
' to (I1) when called for final output. ' 14301 CALL xerrwd (msg, 50, 210, 0, 1, ires, 0, 0, 0.0d0, 0.0d0)
14305 592 y(i) = rwork(i+lyh-1)
14323 601 msg =
'DLSOIBT- ISTATE (=I1) illegal.' 14324 CALL xerrwd (msg, 30, 1, 0, 1, istate, 0, 0, 0.0d0, 0.0d0)
14325 IF (istate .LT. 0)
GO TO 800
14327 602 msg =
'DLSOIBT- ITASK (=I1) illegal. ' 14328 CALL xerrwd (msg, 30, 2, 0, 1, itask, 0, 0, 0.0d0, 0.0d0)
14330 603 msg = .gt.
'DLSOIBT- ISTATE1 but DLSOIBT not initialized. ' 14331 CALL xerrwd (msg, 50, 3, 0, 0, 0, 0, 0, 0.0d0, 0.0d0)
14333 604 msg = .lt.
'DLSOIBT- NEQ (=I1) 1 ' 14334 CALL xerrwd (msg, 30, 4, 0, 1, neq(1), 0, 0, 0.0d0, 0.0d0)
14336 605 msg =
'DLSOIBT- ISTATE = 3 and NEQ increased (I1 to I2). ' 14337 CALL xerrwd (msg, 50, 5, 0, 2, n, neq(1), 0, 0.0d0, 0.0d0)
14339 606 msg =
'DLSOIBT- ITOL (=I1) illegal. ' 14340 CALL xerrwd (msg, 30, 6, 0, 1, itol, 0, 0, 0.0d0, 0.0d0)
14342 607 msg =
'DLSOIBT- IOPT (=I1) illegal. ' 14343 CALL xerrwd (msg, 30, 7, 0, 1, iopt, 0, 0, 0.0d0, 0.0d0)
14345 608 msg =
'DLSOIBT- MF (=I1) illegal. ' 14346 CALL xerrwd (msg, 30, 8, 0, 1, mf, 0, 0, 0.0d0, 0.0d0)
14348 609 msg =
'DLSOIBT- MB (=I1) or NB (=I2) illegal. ' 14349 CALL xerrwd (msg, 40, 9, 0, 2, mb, nb, 0, 0.0d0, 0.0d0)
14351 610 msg = .lt.
'DLSOIBT- NB (=I1) 4 illegal. ' 14352 CALL xerrwd (msg, 40, 10, 0, 1, nb, 0, 0, 0.0d0, 0.0d0)
14354 611 msg = .lt.
'DLSOIBT- MAXORD (=I1) 0 ' 14355 CALL xerrwd (msg, 30, 11, 0, 1, maxord, 0, 0, 0.0d0, 0.0d0)
14357 612 msg = .lt.
'DLSOIBT- MXSTEP (=I1) 0 ' 14358 CALL xerrwd (msg, 30, 12, 0, 1, mxstep, 0, 0, 0.0d0, 0.0d0)
14360 613 msg = .lt.
'DLSOIBT- MXHNIL (=I1) 0 ' 14361 CALL xerrwd (msg, 30, 13, 0, 1, mxhnil, 0, 0, 0.0d0, 0.0d0)
14363 614 msg =
'DLSOIBT- TOUT (=R1) behind T (=R2) ' 14364 CALL xerrwd (msg, 40, 14, 0, 0, 0, 0, 2, tout, t)
14365 msg =
' Integration direction is given by H0 (=R1) ' 14366 CALL xerrwd (msg, 50, 14, 0, 0, 0, 0, 1, h0, 0.0d0)
14368 615 msg = .lt.
'DLSOIBT- HMAX (=R1) 0.0 ' 14369 CALL xerrwd (msg, 30, 15, 0, 0, 0, 0, 1, hmax, 0.0d0)
14371 616 msg = .lt.
'DLSOIBT- HMIN (=R1) 0.0 ' 14372 CALL xerrwd (msg, 30, 16, 0, 0, 0, 0, 1, hmin, 0.0d0)
14374 617 msg=
'DLSOIBT- RWORK length needed, LENRW (=I1), exceeds LRW (=I2)' 14375 CALL xerrwd (msg, 60, 17, 0, 2, lenrw, lrw, 0, 0.0d0, 0.0d0)
14377 618 msg=
'DLSOIBT- IWORK length needed, LENIW (=I1), exceeds LIW (=I2)' 14378 CALL xerrwd (msg, 60, 18, 0, 2, leniw, liw, 0, 0.0d0, 0.0d0)
14380 619 msg = .lt.
'DLSOIBT- RTOL(=I1) is R1 0.0 ' 14381 CALL xerrwd (msg, 40, 19, 0, 1, i, 0, 1, rtoli, 0.0d0)
14383 620 msg = .lt.
'DLSOIBT- ATOL(=I1) is R1 0.0 ' 14384 CALL xerrwd (msg, 40, 20, 0, 1, i, 0, 1, atoli, 0.0d0)
14386 621 ewti = rwork(lewt+i-1)
14387 msg = .le.
'DLSOIBT- EWT(I1) is R1 0.0 ' 14388 CALL xerrwd (msg, 40, 21, 0, 1, i, 0, 1, ewti, 0.0d0)
14390 622 msg=
'DLSOIBT- TOUT(=R1) too close to T(=R2) to start integration.' 14391 CALL xerrwd (msg, 60, 22, 0, 0, 0, 0, 2, tout, t)
14393 623 msg=
'DLSOIBT- ITASK = I1 and TOUT (=R1) behind TCUR - HU (= R2) ' 14394 CALL xerrwd (msg, 60, 23, 0, 1, itask, 0, 2, tout, tp)
14396 624 msg=
'DLSOIBT- ITASK = 4 or 5 and TCRIT (=R1) behind TCUR (=R2) ' 14397 CALL xerrwd (msg, 60, 24, 0, 0, 0, 0, 2, tcrit, tn)
14399 625 msg=
'DLSOIBT- ITASK = 4 or 5 and TCRIT (=R1) behind TOUT (=R2) ' 14400 CALL xerrwd (msg, 60, 25, 0, 0, 0, 0, 2, tcrit, tout)
14402 626 msg =
'DLSOIBT- At start of problem, too much accuracy ' 14403 CALL xerrwd (msg, 50, 26, 0, 0, 0, 0, 0, 0.0d0, 0.0d0)
14404 msg=
' requested for precision of machine.. See TOLSF (=R1) ' 14405 CALL xerrwd (msg, 60, 26, 0, 0, 0, 0, 1, tolsf, 0.0d0)
14408 627 msg =
'DLSOIBT- Trouble in DINTDY. ITASK = I1, TOUT = R1' 14409 CALL xerrwd (msg, 50, 27, 0, 1, itask, 0, 1, tout, 0.0d0)
14414 800 msg =
'DLSOIBT- Run aborted.. apparent infinite loop. ' 14415 CALL xerrwd (msg, 50, 303, 2, 0, 0, 0, 0, 0.0d0, 0.0d0)
14420 SUBROUTINE dlsodis (RES, ADDA, JAC, NEQ, Y, YDOTI, T, TOUT, ITOL,
14421 1 RTOL, ATOL, ITASK, ISTATE, IOPT, RWORK, LRW, IWORK, LIW, MF )
14422 EXTERNAL res, adda, jac
14423 INTEGER neq, itol, itask, istate, iopt, lrw, iwork, liw, mf
14424 DOUBLE PRECISION y, ydoti, t, tout, rtol, atol, rwork
14425 dimension neq(*), y(*), ydoti(*), rtol(*), atol(*), rwork(lrw),
15792 EXTERNAL dprjis, dsolss
15793 DOUBLE PRECISION dumach, dvnorm
15794 INTEGER init, mxstep, mxhnil, nhnil, nslast, nyh, iowns,
15795 1 icf, ierpj, iersl, jcur, jstart, kflag, l,
15796 2 lyh, lewt, lacor, lsavf, lwm, liwm, meth, miter,
15797 3 maxord, maxcor, msbp, mxncf, n, nq, nst, nfe, nje, nqu
15798 INTEGER iplost, iesp, istatc, iys, iba, ibian, ibjan, ibjgp,
15799 1 ipian, ipjan, ipjgp, ipigp, ipr, ipc, ipic, ipisp, iprsp, ipa,
15800 2 lenyh, lenyhm, lenwk, lreq, lrat, lrest, lwmin, moss, msbj,
15801 3 nslj, ngp, nlu, nnz, nsp, nzl, nzu
15802 INTEGER i, i1, i2, ier, igo, iflag, imax, imul, imxer, ipflag,
15803 1 ipgo, irem, ires, j, kgo, lenrat, lenyht, leniw, lenrw,
15804 2 lia, lic, lja, ljc, lrtem, lwtem, lyd0, lyhd, lyhn, mf1,
15805 3 mord, mxhnl0, mxstp0, ncolm
15806 DOUBLE PRECISION rowns,
15807 1 ccmax, el0, h, hmin, hmxi, hu, rc, tn, uround
15808 DOUBLE PRECISION con0, conmin, ccmxj, psmall, rbig, seth
15809 DOUBLE PRECISION atoli, ayi, big, ewti, h0, hmax, hmx, rh, rtoli,
15810 1 tcrit, tdist, tnext, tol, tolsf, tp,
SIZE, sum, w0
15814 SAVE lenrat, mord, mxstp0, mxhnl0
15827 COMMON /dls001/ rowns(209),
15828 1 ccmax, el0, h, hmin, hmxi, hu, rc, tn, uround,
15829 2 init, mxstep, mxhnil, nhnil, nslast, nyh, iowns(6),
15830 3 icf, ierpj, iersl, jcur, jstart, kflag, l,
15831 4 lyh, lewt, lacor, lsavf, lwm, liwm, meth, miter,
15832 5 maxord, maxcor, msbp, mxncf, n, nq, nst, nfe, nje, nqu
15834 COMMON /dlss01/ con0, conmin, ccmxj, psmall, rbig, seth,
15835 1 iplost, iesp, istatc, iys, iba, ibian, ibjan, ibjgp,
15836 2 ipian, ipjan, ipjgp, ipigp, ipr, ipc, ipic, ipisp, iprsp, ipa,
15837 3 lenyh, lenyhm, lenwk, lreq, lrat, lrest, lwmin, moss, msbj,
15838 4 nslj, ngp, nlu, nnz, nsp, nzl, nzu
15840 DATA mord(1),mord(2)/12,5/, mxstp0/500/, mxhnl0/10/
15856 IF (istate .LT. 0 .OR. istate .GT. 3)
GO TO 601
15857 IF (itask .LT. 1 .OR. itask .GT. 5)
GO TO 602
15858 IF (istate .LE. 1)
GO TO 10
15859 IF (init .EQ. 0)
GO TO 603
15860 IF (istate .EQ. 2)
GO TO 200
15863 IF (tout .EQ. t)
RETURN 15875 20
IF (neq(1) .LE. 0)
GO TO 604
15876 IF (istate .LE. 1)
GO TO 25
15877 IF (neq(1) .GT. n)
GO TO 605
15879 IF (itol .LT. 1 .OR. itol .GT. 4)
GO TO 606
15880 IF (iopt .LT. 0 .OR. iopt .GT. 1)
GO TO 607
15882 mf1 = mf - 100*moss
15884 miter = mf1 - 10*meth
15885 IF (moss .LT. 0 .OR. moss .GT. 4)
GO TO 608
15886 IF (miter .EQ. 2 .AND. moss .EQ. 1) moss = moss + 1
15887 IF (miter .EQ. 2 .AND. moss .EQ. 3) moss = moss + 1
15888 IF (miter .EQ. 1 .AND. moss .EQ. 2) moss = moss - 1
15889 IF (miter .EQ. 1 .AND. moss .EQ. 4) moss = moss - 1
15890 IF (meth .LT. 1 .OR. meth .GT. 2)
GO TO 608
15891 IF (miter .LT. 1 .OR. miter .GT. 2)
GO TO 608
15893 IF (iopt .EQ. 1)
GO TO 40
15894 maxord = mord(meth)
15897 IF (istate .LE. 1) h0 = 0.0d0
15901 40 maxord = iwork(5)
15902 IF (maxord .LT. 0)
GO TO 611
15903 IF (maxord .EQ. 0) maxord = 100
15904 maxord = min(maxord,mord(meth))
15906 IF (mxstep .LT. 0)
GO TO 612
15907 IF (mxstep .EQ. 0) mxstep = mxstp0
15909 IF (mxhnil .LT. 0)
GO TO 613
15910 IF (mxhnil .EQ. 0) mxhnil = mxhnl0
15911 IF (istate .GT. 1)
GO TO 50
15913 IF ((tout - t)*h0 .LT. 0.0d0)
GO TO 614
15915 IF (hmax .LT. 0.0d0)
GO TO 615
15917 IF (hmax .GT. 0.0d0) hmxi = 1.0d0/hmax
15919 IF (hmin .LT. 0.0d0)
GO TO 616
15924 IF (itol .GE. 3) rtoli = rtol(i)
15925 IF (itol .EQ. 2 .OR. itol .EQ. 4) atoli = atol(i)
15926 IF (rtoli .LT. 0.0d0)
GO TO 619
15927 IF (atoli .LT. 0.0d0)
GO TO 620
15946 IF (istate .LE. 1) nyh = n
15947 IF (miter .EQ. 1) lwmin = 4*n + 10*n/lrat
15948 IF (miter .EQ. 2) lwmin = 4*n + 11*n/lrat
15949 lenyh = (maxord+1)*nyh
15950 lrest = lenyh + 3*n
15951 lenrw = 20 + lwmin + lrest
15954 IF (moss .NE. 1 .AND. moss .NE. 2) leniw = leniw + n + 1
15956 IF (lenrw .GT. lrw)
GO TO 617
15957 IF (leniw .GT. liw)
GO TO 618
15959 IF (moss .NE. 1 .AND. moss .NE. 2)
15960 1 leniw = leniw + iwork(lia+n) - 1
15962 IF (leniw .GT. liw)
GO TO 618
15967 IF (moss .EQ. 0) leniw = leniw + n + 1
15969 IF (leniw .GT. liw)
GO TO 618
15970 IF (moss .EQ. 0) leniw = leniw + iwork(lic+n) - 1
15972 IF (leniw .GT. liw)
GO TO 618
15977 IF (istate .LE. 1) nq = istate
15978 ncolm = min(nq+1,maxord+2)
15982 IF (istate .EQ. 3) imul = moss
15983 IF (istate .EQ. 3 .AND. moss .EQ. 3) imul = 1
15984 IF (moss .EQ. 2 .OR. moss .EQ. 4) imul = 3
15985 lrtem = lenyht + imul*n
15986 lwtem = lrw - 20 - lrtem
15989 lsavf = lyhn + lenyht
15993 IF (istate .LE. 1)
GO TO 100
16004 imax = lyhn - 1 + lenyhm
16006 IF (lyhd .LT. 0)
THEN 16007 DO 72 i = lyhn,imax
16008 j = imax + lyhn - i
16009 72 rwork(j) = rwork(j+lyhd)
16011 IF (lyhd .GT. 0)
THEN 16012 DO 76 i = lyhn,imax
16013 76 rwork(i) = rwork(i+lyhd)
16017 IF (moss .NE. 2 .AND. moss .NE. 4)
GO TO 85
16019 CALL dewset (n,itol,rtol,atol,rwork(lyh),rwork(lewt))
16021 IF (rwork(i+lewt-1) .LE. 0.0d0)
GO TO 621
16022 82 rwork(i+lewt-1) = 1.0d0/rwork(i+lewt-1)
16025 lsavf = min(lsavf,lrw)
16026 lewt = min(lewt,lrw)
16027 lacor = min(lacor,lrw)
16028 CALL diprepi (neq, y, ydoti, rwork, iwork(lia), iwork(lja),
16029 1 iwork(lic), iwork(ljc), ipflag, res, jac, adda)
16030 lenrw = lwm - 1 + lenwk + lrest
16032 IF (ipflag .NE. -1) iwork(23) = ipian
16033 IF (ipflag .NE. -1) iwork(24) = ipjan
16035 GO TO (90, 628, 629, 630, 631, 632, 633, 634, 634), ipgo
16038 IF (lenrw .GT. lrw)
GO TO 617
16041 IF (nq .LE. maxord)
GO TO 94
16044 92 ydoti(i) = rwork(i+lsavf-1)
16045 94
IF (n .EQ. nyh)
GO TO 200
16048 i2 = lyh + (maxord + 1)*nyh - 1
16049 IF (i1 .GT. i2)
GO TO 200
16051 95 rwork(i) = 0.0d0
16074 105 rwork(i+lyh-1) = y(i)
16075 IF (istate .NE. 1)
GO TO 108
16079 106 rwork(i+lyd0-1) = ydoti(i)
16082 CALL dewset (n,itol,rtol,atol,rwork(lyh),rwork(lewt))
16084 IF (rwork(i+lewt-1) .LE. 0.0d0)
GO TO 621
16085 110 rwork(i+lewt-1) = 1.0d0/rwork(i+lewt-1)
16087 lacor = min(lacor,lrw)
16088 CALL diprepi (neq, y, ydoti, rwork, iwork(lia), iwork(lja),
16089 1 iwork(lic), iwork(ljc), ipflag, res, jac, adda)
16090 lenrw = lwm - 1 + lenwk + lrest
16092 IF (ipflag .NE. -1) iwork(23) = ipian
16093 IF (ipflag .NE. -1) iwork(24) = ipjan
16095 GO TO (115, 628, 629, 630, 631, 632, 633, 634, 634), ipgo
16096 115 iwork(22) = lyh
16097 IF (lenrw .GT. lrw)
GO TO 617
16100 IF (istate .NE. 0)
GO TO 120
16101 CALL dainvgs (neq, t, y, rwork(lwm), rwork(lwm), rwork(lacor),
16102 1 rwork(lyd0), ier, res, adda)
16105 GO TO (120, 565, 560, 560), igo
16108 IF (itask .NE. 4 .AND. itask .NE. 5)
GO TO 125
16110 IF ((tcrit - tout)*(tout - t) .LT. 0.0d0)
GO TO 625
16111 IF (h0 .NE. 0.0d0 .AND. (t + h0 - tcrit)*h0 .GT. 0.0d0)
16114 125 uround = dumach()
16116 rwork(lwm) = sqrt(uround)
16143 IF (h0 .NE. 0.0d0)
GO TO 180
16144 tdist = abs(tout - t)
16145 w0 = max(abs(t),abs(tout))
16146 IF (tdist .LT. 2.0d0*uround*w0)
GO TO 622
16148 IF (itol .LE. 2)
GO TO 145
16150 140 tol = max(tol,rtol(i))
16151 145
IF (tol .GT. 0.0d0)
GO TO 160
16154 IF (itol .EQ. 2 .OR. itol .EQ. 4) atoli = atol(i)
16156 IF (ayi .NE. 0.0d0) tol = max(tol,atoli/ayi)
16158 160 tol = max(tol,100.0d0*uround)
16159 tol = min(tol,0.001d0)
16160 sum = dvnorm(n, rwork(lyd0), rwork(lewt))
16161 sum = 1.0d0/(tol*w0*w0) + tol*sum**2
16162 h0 = 1.0d0/sqrt(sum)
16164 h0 = sign(h0,tout-t)
16166 180 rh = abs(h0)*hmxi
16167 IF (rh .GT. 1.0d0) h0 = h0/rh
16171 190 rwork(i+lyd0-1) = h0*rwork(i+lyd0-1)
16179 GO TO (210, 250, 220, 230, 240), itask
16180 210
IF ((tn - tout)*h .LT. 0.0d0)
GO TO 250
16181 CALL dintdy (tout, 0, rwork(lyh), nyh, y, iflag)
16182 IF (iflag .NE. 0)
GO TO 627
16185 220 tp = tn - hu*(1.0d0 + 100.0d0*uround)
16186 IF ((tp - tout)*h .GT. 0.0d0)
GO TO 623
16187 IF ((tn - tout)*h .LT. 0.0d0)
GO TO 250
16189 230 tcrit = rwork(1)
16190 IF ((tn - tcrit)*h .GT. 0.0d0)
GO TO 624
16191 IF ((tcrit - tout)*h .LT. 0.0d0)
GO TO 625
16192 IF ((tn - tout)*h .LT. 0.0d0)
GO TO 245
16193 CALL dintdy (tout, 0, rwork(lyh), nyh, y, iflag)
16194 IF (iflag .NE. 0)
GO TO 627
16197 240 tcrit = rwork(1)
16198 IF ((tn - tcrit)*h .GT. 0.0d0)
GO TO 624
16199 245 hmx = abs(tn) + abs(h)
16200 ihit = abs(tn - tcrit) .LE. 100.0d0*uround*hmx
16201 IF (ihit)
GO TO 400
16202 tnext = tn + h*(1.0d0 + 4.0d0*uround)
16203 IF ((tnext - tcrit)*h .LE. 0.0d0)
GO TO 250
16204 h = (tcrit - tn)*(1.0d0 - 4.0d0*uround)
16205 IF (istate .EQ. 2) jstart = -2
16218 IF ((nst-nslast) .GE. mxstep)
GO TO 500
16219 CALL dewset (n, itol, rtol, atol, rwork(lyh), rwork(lewt))
16221 IF (rwork(i+lewt-1) .LE. 0.0d0)
GO TO 510
16222 260 rwork(i+lewt-1) = 1.0d0/rwork(i+lewt-1)
16223 270 tolsf = uround*dvnorm(n, rwork(lyh), rwork(lewt))
16224 IF (tolsf .LE. 1.0d0)
GO TO 280
16225 tolsf = tolsf*2.0d0
16226 IF (nst .EQ. 0)
GO TO 626
16228 280
IF ((tn + h) .NE. tn)
GO TO 290
16230 IF (nhnil .GT. mxhnil)
GO TO 290
16231 msg =
'DLSODIS- Warning..Internal T (=R1) and H (=R2) are' 16232 CALL xerrwd (msg, 50, 101, 0, 0, 0, 0, 0, 0.0d0, 0.0d0)
16233 msg=
' such that in the machine, T + H = T on the next step ' 16234 CALL xerrwd (msg, 60, 101, 0, 0, 0, 0, 0, 0.0d0, 0.0d0)
16235 msg =
' (H = step size). Solver will continue anyway.' 16236 CALL xerrwd (msg, 50, 101, 0, 0, 0, 0, 2, tn, h)
16237 IF (nhnil .LT. mxhnil)
GO TO 290
16238 msg =
'DLSODIS- Above warning has been issued I1 times. ' 16239 CALL xerrwd (msg, 50, 102, 0, 0, 0, 0, 0, 0.0d0, 0.0d0)
16240 msg =
' It will not be issued again for this problem.' 16241 CALL xerrwd (msg, 50, 102, 0, 1, mxhnil, 0, 0, 0.0d0, 0.0d0)
16248 CALL dstodi (neq, y, rwork(lyh), nyh, rwork(lyh), rwork(lewt),
16249 1 ydoti, rwork(lsavf), rwork(lacor), rwork(lwm),
16250 2 rwork(lwm), res, adda, jac, dprjis, dsolss )
16252 GO TO (300, 530, 540, 400, 550, 555), kgo
16263 GO TO (310, 400, 330, 340, 350), itask
16265 310
iF ((tn - tout)*h .LT. 0.0d0)
GO TO 250
16266 CALL dintdy (tout, 0, rwork(lyh), nyh, y, iflag)
16270 330
IF ((tn - tout)*h .GE. 0.0d0)
GO TO 400
16273 340
IF ((tn - tout)*h .LT. 0.0d0)
GO TO 345
16274 CALL dintdy (tout, 0, rwork(lyh), nyh, y, iflag)
16277 345 hmx = abs(tn) + abs(h)
16278 ihit = abs(tn - tcrit) .LE. 100.0d0*uround*hmx
16279 IF (ihit)
GO TO 400
16280 tnext = tn + h*(1.0d0 + 4.0d0*uround)
16281 IF ((tnext - tcrit)*h .LE. 0.0d0)
GO TO 250
16282 h = (tcrit - tn)*(1.0d0 - 4.0d0*uround)
16286 350 hmx = abs(tn) + abs(h)
16287 ihit = abs(tn - tcrit) .LE. 100.0d0*uround*hmx
16296 410 y(i) = rwork(i+lyh-1)
16298 IF (itask .NE. 4 .AND. itask .NE. 5)
GO TO 420
16299 IF (ihit) t = tcrit
16301 IF ( kflag .EQ. -3 ) istate = 3
16325 500 msg =
'DLSODIS- At current T (=R1), MXSTEP (=I1) steps ' 16326 CALL xerrwd (msg, 50, 201, 0, 0, 0, 0, 0, 0.0d0, 0.0d0)
16327 msg =
' taken on this call before reaching TOUT ' 16328 CALL xerrwd (msg, 50, 201, 0, 1, mxstep, 0, 1, tn, 0.0d0)
16332 510 ewti = rwork(lewt+i-1)
16333 msg = .le.
'DLSODIS- At T (=R1), EWT(I1) has become R2 0.' 16334 CALL xerrwd (msg, 50, 202, 0, 1, i, 0, 2, tn, ewti)
16338 520 msg =
'DLSODIS- At T (=R1), too much accuracy requested ' 16339 CALL xerrwd (msg, 50, 203, 0, 0, 0, 0, 0, 0.0d0, 0.0d0)
16340 msg =
' for precision of machine.. See TOLSF (=R2) ' 16341 CALL xerrwd (msg, 50, 203, 0, 0, 0, 0, 2, tn, tolsf)
16346 530 msg =
'DLSODIS- At T (=R1) and step size H (=R2), the ' 16347 CALL xerrwd (msg, 50, 204, 0, 0, 0, 0, 0, 0.0d0, 0.0d0)
16348 msg=
' error test failed repeatedly or with ABS(H) = HMIN ' 16349 CALL xerrwd (msg, 60, 204, 0, 0, 0, 0, 2, tn, h)
16353 540 msg =
'DLSODIS- At T (=R1) and step size H (=R2), the ' 16354 CALL xerrwd (msg, 50, 205, 0, 0, 0, 0, 0, 0.0d0, 0.0d0)
16355 msg =
' corrector convergence failed repeatedly ' 16356 CALL xerrwd (msg, 50, 205, 0, 0, 0, 0, 0, 0.0d0, 0.0d0)
16357 msg =
' or with ABS(H) = HMIN ' 16358 CALL xerrwd (msg, 30, 205, 0, 0, 0, 0, 2, tn, h)
16362 550 msg =
'DLSODIS- At T (=R1) residual routine returned ' 16363 CALL xerrwd (msg, 50, 206, 0, 0, 0, 0, 0, 0.0d0, 0.0d0)
16364 msg =
' error IRES = 3 repeatedly.' 16365 CALL xerrwd (msg, 30, 206, 1, 0, 0, 0, 0, tn, 0.0d0)
16369 555 msg =
'DLSODIS- At T (=R1) and step size H (=R2), a fatal' 16370 CALL xerrwd (msg, 50, 207, 0, 0, 0, 0, 0, 0.0d0, 0.0d0)
16371 msg =
' error flag was returned by CDRV (by way of ' 16372 CALL xerrwd (msg, 50, 207, 0, 0, 0, 0, 0, 0.0d0, 0.0d0)
16373 msg =
' Subroutine DPRJIS or DSOLSS) ' 16374 CALL xerrwd (msg, 40, 207, 0, 0, 0, 0, 2, tn, h)
16378 560 msg=
'DLSODIS- Attempt to initialize dy/dt failed because matrix A' 16379 CALL xerrwd (msg, 60, 208, 0, 0, 0, 0, 0, 0.0d0, 0.0d0)
16380 msg=
' was singular. CDRV returned zero pivot error flag. ' 16381 CALL xerrwd (msg, 60, 208, 0, 0, 0, 0, 0, 0.0d0, 0.0d0)
16382 msg =
'DAINVGS set its error flag to IER = (I1)' 16383 CALL xerrwd (msg, 40, 208, 0, 1, ier, 0, 0, 0.0d0, 0.0d0)
16387 565 msg =
'DLSODIS- Attempt to initialize dy/dt failed ' 16388 CALL xerrwd (msg, 50, 209, 0, 0, 0, 0, 0, 0.0d0, 0.0d0)
16389 msg =
' because residual routine set its error flag ' 16390 CALL xerrwd (msg, 50, 209, 0, 0, 0, 0, 0, 0.0d0, 0.0d0)
16391 msg =
' to IRES = (I1)' 16392 CALL xerrwd (msg, 20, 209, 0, 1, ier, 0, 0, 0.0d0, 0.0d0)
16399 SIZE = abs(rwork(i+lacor-1)*rwork(i+lewt-1))
16400 IF (big .GE. size)
GO TO 575
16406 580 lyd0 = lyh + nyh
16408 rwork(i+lsavf-1) = rwork(i+lyd0-1) / h
16409 585 y(i) = rwork(i+lyh-1)
16411 CALL res (neq, tn, y, rwork(lsavf), ydoti, ires)
16413 IF ( ires .LE. 1 )
GO TO 595
16414 msg =
'DLSODIS- Residual routine set its flag IRES ' 16415 CALL xerrwd (msg, 50, 210, 0, 0, 0, 0, 0, 0.0d0, 0.0d0)
16416 msg =
' to (I1) when called for final output. ' 16417 CALL xerrwd (msg, 50, 210, 0, 1, ires, 0, 0, 0.0d0, 0.0d0)
16421 592 y(i) = rwork(i+lyh-1)
16444 601 msg =
'DLSODIS- ISTATE (=I1) illegal.' 16445 CALL xerrwd (msg, 30, 1, 0, 1, istate, 0, 0, 0.0d0, 0.0d0)
16446 IF (istate .LT. 0)
GO TO 800
16448 602 msg =
'DLSODIS- ITASK (=I1) illegal. ' 16449 CALL xerrwd (msg, 30, 2, 0, 1, itask, 0, 0, 0.0d0, 0.0d0)
16451 603 msg = .gt.
'DLSODIS-ISTATE 1 but DLSODIS not initialized.' 16452 CALL xerrwd (msg, 50, 3, 0, 0, 0, 0, 0, 0.0d0, 0.0d0)
16454 604 msg = .lt.
'DLSODIS- NEQ (=I1) 1 ' 16455 CALL xerrwd (msg, 30, 4, 0, 1, neq(1), 0, 0, 0.0d0, 0.0d0)
16457 605 msg =
'DLSODIS- ISTATE = 3 and NEQ increased (I1 to I2). ' 16458 CALL xerrwd (msg, 50, 5, 0, 2, n, neq(1), 0, 0.0d0, 0.0d0)
16460 606 msg =
'DLSODIS- ITOL (=I1) illegal. ' 16461 CALL xerrwd (msg, 30, 6, 0, 1, itol, 0, 0, 0.0d0, 0.0d0)
16463 607 msg =
'DLSODIS- IOPT (=I1) illegal. ' 16464 CALL xerrwd (msg, 30, 7, 0, 1, iopt, 0, 0, 0.0d0, 0.0d0)
16466 608 msg =
'DLSODIS- MF (=I1) illegal. ' 16467 CALL xerrwd (msg, 30, 8, 0, 1, mf, 0, 0, 0.0d0, 0.0d0)
16469 611 msg = .lt.
'DLSODIS- MAXORD (=I1) 0 ' 16470 CALL xerrwd (msg, 30, 11, 0, 1, maxord, 0, 0, 0.0d0, 0.0d0)
16472 612 msg = .lt.
'DLSODIS- MXSTEP (=I1) 0 ' 16473 CALL xerrwd (msg, 30, 12, 0, 1, mxstep, 0, 0, 0.0d0, 0.0d0)
16475 613 msg = .lt.
'DLSODIS- MXHNIL (=I1) 0 ' 16476 CALL xerrwd (msg, 30, 13, 0, 1, mxhnil, 0, 0, 0.0d0, 0.0d0)
16478 614 msg =
'DLSODIS- TOUT (=R1) behind T (=R2) ' 16479 CALL xerrwd (msg, 40, 14, 0, 0, 0, 0, 2, tout, t)
16480 msg =
' Integration direction is given by H0 (=R1) ' 16481 CALL xerrwd (msg, 50, 14, 0, 0, 0, 0, 1, h0, 0.0d0)
16483 615 msg = .lt.
'DLSODIS- HMAX (=R1) 0.0 ' 16484 CALL xerrwd (msg, 30, 15, 0, 0, 0, 0, 1, hmax, 0.0d0)
16486 616 msg = .lt.
'DLSODIS- HMIN (=R1) 0.0 ' 16487 CALL xerrwd (msg, 30, 16, 0, 0, 0, 0, 1, hmin, 0.0d0)
16489 617 msg =
'DLSODIS- RWORK length is insufficient to proceed. ' 16490 CALL xerrwd (msg, 50, 17, 0, 0, 0, 0, 0, 0.0d0, 0.0d0)
16491 msg=.ge.
' Length needed is LENRW (=I1), exceeds LRW (=I2)' 16492 CALL xerrwd (msg, 60, 17, 0, 2, lenrw, lrw, 0, 0.0d0, 0.0d0)
16494 618 msg =
'DLSODIS- IWORK length is insufficient to proceed. ' 16495 CALL xerrwd (msg, 50, 18, 0, 0, 0, 0, 0, 0.0d0, 0.0d0)
16496 msg=.ge.
' Length needed is LENIW (=I1), exceeds LIW (=I2)' 16497 CALL xerrwd (msg, 60, 18, 0, 2, leniw, liw, 0, 0.0d0, 0.0d0)
16499 619 msg = .lt.
'DLSODIS- RTOL(=I1) is R1 0.0 ' 16500 CALL xerrwd (msg, 40, 19, 0, 1, i, 0, 1, rtoli, 0.0d0)
16502 620 msg = .lt.
'DLSODIS- ATOL(=I1) is R1 0.0 ' 16503 CALL xerrwd (msg, 40, 20, 0, 1, i, 0, 1, atoli, 0.0d0)
16505 621 ewti = rwork(lewt+i-1)
16506 msg = .le.
'DLSODIS- EWT(I1) is R1 0.0 ' 16507 CALL xerrwd (msg, 40, 21, 0, 1, i, 0, 1, ewti, 0.0d0)
16509 622 msg=
'DLSODIS- TOUT(=R1) too close to T(=R2) to start integration.' 16510 CALL xerrwd (msg, 60, 22, 0, 0, 0, 0, 2, tout, t)
16512 623 msg=
'DLSODIS- ITASK = I1 and TOUT (=R1) behind TCUR - HU (= R2) ' 16513 CALL xerrwd (msg, 60, 23, 0, 1, itask, 0, 2, tout, tp)
16515 624 msg=
'DLSODIS- ITASK = 4 or 5 and TCRIT (=R1) behind TCUR (=R2) ' 16516 CALL xerrwd (msg, 60, 24, 0, 0, 0, 0, 2, tcrit, tn)
16518 625 msg=
'DLSODIS- ITASK = 4 or 5 and TCRIT (=R1) behind TOUT (=R2) ' 16519 CALL xerrwd (msg, 60, 25, 0, 0, 0, 0, 2, tcrit, tout)
16521 626 msg =
'DLSODIS- At start of problem, too much accuracy ' 16522 CALL xerrwd (msg, 50, 26, 0, 0, 0, 0, 0, 0.0d0, 0.0d0)
16523 msg=
' requested for precision of machine.. See TOLSF (=R1) ' 16524 CALL xerrwd (msg, 60, 26, 0, 0, 0, 0, 1, tolsf, 0.0d0)
16527 627 msg =
'DLSODIS- Trouble in DINTDY. ITASK = I1, TOUT = R1' 16528 CALL xerrwd (msg, 50, 27, 0, 1, itask, 0, 1, tout, 0.0d0)
16530 628 msg=
'DLSODIS- RWORK length insufficient (for Subroutine DPREPI). ' 16531 CALL xerrwd (msg, 60, 28, 0, 0, 0, 0, 0, 0.0d0, 0.0d0)
16532 msg=.ge.
' Length needed is LENRW (=I1), exceeds LRW (=I2)' 16533 CALL xerrwd (msg, 60, 28, 0, 2, lenrw, lrw, 0, 0.0d0, 0.0d0)
16535 629 msg=
'DLSODIS- RWORK length insufficient (for Subroutine JGROUP). ' 16536 CALL xerrwd (msg, 60, 29, 0, 0, 0, 0, 0, 0.0d0, 0.0d0)
16537 msg=.ge.
' Length needed is LENRW (=I1), exceeds LRW (=I2)' 16538 CALL xerrwd (msg, 60, 29, 0, 2, lenrw, lrw, 0, 0.0d0, 0.0d0)
16540 630 msg=
'DLSODIS- RWORK length insufficient (for Subroutine ODRV). ' 16541 CALL xerrwd (msg, 60, 30, 0, 0, 0, 0, 0, 0.0d0, 0.0d0)
16542 msg=.ge.
' Length needed is LENRW (=I1), exceeds LRW (=I2)' 16543 CALL xerrwd (msg, 60, 30, 0, 2, lenrw, lrw, 0, 0.0d0, 0.0d0)
16545 631 msg=
'DLSODIS- Error from ODRV in Yale Sparse Matrix Package. ' 16546 CALL xerrwd (msg, 60, 31, 0, 0, 0, 0, 0, 0.0d0, 0.0d0)
16548 irem = iys - imul*n
16549 msg=
' At T (=R1), ODRV returned error flag = I1*NEQ + I2. ' 16550 CALL xerrwd (msg, 60, 31, 0, 2, imul, irem, 1, tn, 0.0d0)
16552 632 msg=
'DLSODIS- RWORK length insufficient (for Subroutine CDRV). ' 16553 CALL xerrwd (msg, 60, 32, 0, 0, 0, 0, 0, 0.0d0, 0.0d0)
16554 msg=.ge.
' Length needed is LENRW (=I1), exceeds LRW (=I2)' 16555 CALL xerrwd (msg, 60, 32, 0, 2, lenrw, lrw, 0, 0.0d0, 0.0d0)
16557 633 msg=
'DLSODIS- Error from CDRV in Yale Sparse Matrix Package. ' 16558 CALL xerrwd (msg, 60, 33, 0, 0, 0, 0, 0, 0.0d0, 0.0d0)
16560 irem = iys - imul*n
16561 msg=
' At T (=R1), CDRV returned error flag = I1*NEQ + I2. ' 16562 CALL xerrwd (msg, 60, 33, 0, 2, imul, irem, 1, tn, 0.0d0)
16563 IF (imul .EQ. 2)
THEN 16564 msg=
' Duplicate entry in sparsity structure descriptors. ' 16565 CALL xerrwd (msg, 60, 33, 0, 0, 0, 0, 0, 0.0d0, 0.0d0)
16567 IF (imul .EQ. 3 .OR. imul .EQ. 6)
THEN 16568 msg=
' Insufficient storage for NSFC (called by CDRV). ' 16569 CALL xerrwd (msg, 60, 33, 0, 0, 0, 0, 0, 0.0d0, 0.0d0)
16572 634 msg=
'DLSODIS- At T (=R1) residual routine (called by DPREPI) ' 16573 CALL xerrwd (msg, 60, 34, 0, 0, 0, 0, 0, 0.0d0, 0.0d0)
16575 msg =
' returned error IRES (=I1)' 16576 CALL xerrwd (msg, 30, 34, 0, 1, ier, 0, 1, tn, 0.0d0)
16581 800 msg =
'DLSODIS- Run aborted.. apparent infinite loop. ' 16582 CALL xerrwd (msg, 50, 303, 2, 0, 0, 0, 0, 0.0d0, 0.0d0)