1 SUBROUTINE dlsode (F, NEQ, Y, T, TOUT, ITOL, RTOL, ATOL, ITASK,
2 1 istate, iopt, rwork, lrw, iwork, liw, jac, mf)
4 INTEGER NEQ, ITOL, ITASK, ISTATE, IOPT, LRW, IWORK, LIW, MF
5 DOUBLE PRECISION Y, T, TOUT, RTOL, ATOL, RWORK
6 dimension neq(*), y(*), rtol(*), atol(*), rwork(lrw), iwork(liw)
947 INTEGER ILLIN, INIT, LYH, LEWT, LACOR, LSAVF, LWM, LIWM,
948 1 mxstep, mxhnil, nhnil, ntrep, nslast, nyh,
949 2 ialth, ipup, lmax, meo, nqnyh, nslp
950 INTEGER ICF, IERPJ, IERSL, JCUR, JSTART, KFLAG, L, METH, MITER,
951 1 maxord, maxcor, msbp, mxncf, n, nq, nst, nfe, nje, nqu
952 INTEGER I, I1, I2, IFLAG, IMXER, KGO, LF0,
953 1 leniw, lenrw, lenwm, ml, mord, mu, mxhnl0, mxstp0
954 DOUBLE PRECISION CONIT, CRATE, EL, ELCO, HOLD, RMAX, TESCO,
955 1 ccmax, el0, h, hmin, hmxi, hu, rc, tn, uround
956 DOUBLE PRECISION ATOLI, AYI, BIG, EWTI, H0, HMAX, HMX, RH, RTOLI,
957 1 tcrit, tdist, tnext, tol, tolsf, tp,
SIZE, sum, w0,
975 COMMON /ls0001/ conit, crate, el(13), elco(13,12),
976 1 hold, rmax, tesco(3,12),
977 1 ccmax, el0, h, hmin, hmxi, hu, rc, tn, uround,
978 2 illin, init, lyh, lewt, lacor, lsavf, lwm, liwm,
979 3 mxstep, mxhnil, nhnil, ntrep, nslast, nyh,
980 3 ialth, ipup, lmax, meo, nqnyh, nslp,
981 4 icf, ierpj, iersl, jcur, jstart, kflag, l, meth, miter,
982 5 maxord, maxcor, msbp, mxncf, n, nq, nst, nfe, nje, nqu
984 DATA mord(1),mord(2)/12,5/, mxstp0/500/, mxhnl0/10/
993 IF (istate .LT. 1 .OR. istate .GT. 3) go
to 601
994 IF (itask .LT. 1 .OR. itask .GT. 5) go
to 602
995 IF (istate .EQ. 1) go
to 10
996 IF (init .EQ. 0) go
to 603
997 IF (istate .EQ. 2) go
to 200
1000 IF (tout .EQ. t) go
to 430
1011 IF (neq(1) .LE. 0) go
to 604
1012 IF (istate .EQ. 1) go
to 25
1013 IF (neq(1) .GT. n) go
to 605
1015 IF (itol .LT. 1 .OR. itol .GT. 4) go
to 606
1016 IF (iopt .LT. 0 .OR. iopt .GT. 1) go
to 607
1018 miter = mf - 10*meth
1019 IF (meth .LT. 1 .OR. meth .GT. 2) go
to 608
1020 IF (miter .LT. 0 .OR. miter .GT. 5) go
to 608
1021 IF (miter .LE. 3) go
to 30
1024 IF (ml .LT. 0 .OR. ml .GE. n) go
to 609
1025 IF (mu .LT. 0 .OR. mu .GE. n) go
to 610
1028 IF (iopt .EQ. 1) go
to 40
1032 IF (istate .EQ. 1) h0 = 0.0d0
1036 40 maxord = iwork(5)
1037 IF (maxord .LT. 0) go
to 611
1038 IF (maxord .EQ. 0) maxord = 100
1039 maxord = min0(maxord,mord(meth))
1041 IF (mxstep .LT. 0) go
to 612
1042 IF (mxstep .EQ. 0) mxstep = mxstp0
1044 IF (mxhnil .LT. 0) go
to 613
1045 IF (mxhnil .EQ. 0) mxhnil = mxhnl0
1046 IF (istate .NE. 1) go
to 50
1048 IF ((tout - t)*h0 .LT. 0.0d0) go
to 614
1050 IF (hmax .LT. 0.0d0) go
to 615
1052 IF (hmax .GT. 0.0d0) hmxi = 1.0d0/hmax
1054 IF (hmin .LT. 0.0d0) go
to 616
1062 IF (istate .EQ. 1) nyh = n
1063 lwm = lyh + (maxord + 1)*nyh
1064 IF (miter .EQ. 0) lenwm = 0
1065 IF (miter .EQ. 1 .OR. miter .EQ. 2) lenwm = n*n + 2
1066 IF (miter .EQ. 3) lenwm = n + 2
1067 IF (miter .GE. 4) lenwm = (2*ml + mu + 1)*n + 2
1071 lenrw = lacor + n - 1
1075 IF (miter .EQ. 0 .OR. miter .EQ. 3) leniw = 20
1077 IF (lenrw .GT. lrw) go
to 617
1078 IF (leniw .GT. liw) go
to 618
1083 IF (itol .GE. 3) rtoli = rtol(i)
1084 IF (itol .EQ. 2 .OR. itol .EQ. 4) atoli = atol(i)
1085 IF (rtoli .LT. 0.0d0) go
to 619
1086 IF (atoli .LT. 0.0d0) go
to 620
1088 IF (istate .EQ. 1) go
to 100
1091 IF (nq .LE. maxord) go
to 90
1094 80 rwork(i+lsavf-1) = rwork(i+lwm-1)
1096 90
IF (miter .GT. 0) rwork(lwm) = dsqrt(uround)
1097 IF (n .EQ. nyh) go
to 200
1100 i2 = lyh + (maxord + 1)*nyh - 1
1101 IF (i1 .GT. i2) go
to 200
1112 100 uround = d1mach(4)
1114 IF (itask .NE. 4 .AND. itask .NE. 5) go
to 110
1116 IF ((tcrit - tout)*(tout - t) .LT. 0.0d0) go
to 625
1117 IF (h0 .NE. 0.0d0 .AND. (t + h0 - tcrit)*h0 .GT. 0.0d0)
1120 IF (miter .GT. 0) rwork(lwm) = dsqrt(uround)
1134 CALL
f(neq, t, y, rwork(lf0), ierr)
1135 IF (ierr .LT. 0)
THEN
1142 115 rwork(i+lyh-1) = y(i)
1146 CALL
ewset(n, itol, rtol, atol, rwork(lyh), rwork(lewt))
1148 IF (rwork(i+lewt-1) .LE. 0.0d0) go
to 621
1149 120 rwork(i+lewt-1) = 1.0d0/rwork(i+lewt-1)
1166 IF (h0 .NE. 0.0d0) go
to 180
1167 tdist = dabs(tout - t)
1168 w0 = dmax1(dabs(t),dabs(tout))
1169 IF (tdist .LT. 2.0d0*uround*w0) go
to 622
1171 IF (itol .LE. 2) go
to 140
1173 130 tol = dmax1(tol,rtol(i))
1174 140
IF (tol .GT. 0.0d0) go
to 160
1177 IF (itol .EQ. 2 .OR. itol .EQ. 4) atoli = atol(i)
1179 IF (ayi .NE. 0.0d0) tol = dmax1(tol,atoli/ayi)
1181 160 tol = dmax1(tol,100.0d0*uround)
1182 tol = dmin1(tol,0.001d0)
1183 sum =
vnorm(n, rwork(lf0), rwork(lewt))
1184 sum = 1.0d0/(tol*w0*w0) + tol*sum**2
1185 h0 = 1.0d0/dsqrt(sum)
1186 h0 = dmin1(h0,tdist)
1187 h0 = dsign(h0,tout-t)
1189 180 rh = dabs(h0)*hmxi
1190 IF (rh .GT. 1.0d0) h0 = h0/rh
1194 190 rwork(i+lf0-1) = h0*rwork(i+lf0-1)
1202 go
to(210, 250, 220, 230, 240), itask
1203 210
IF ((tn - tout)*h .LT. 0.0d0) go
to 250
1204 CALL
intdy(tout, 0, rwork(lyh), nyh, y, iflag)
1205 IF (iflag .NE. 0) go
to 627
1208 220 tp = tn - hu*(1.0d0 + 100.0d0*uround)
1209 IF ((tp - tout)*h .GT. 0.0d0) go
to 623
1210 IF ((tn - tout)*h .LT. 0.0d0) go
to 250
1212 230 tcrit = rwork(1)
1213 IF ((tn - tcrit)*h .GT. 0.0d0) go
to 624
1214 IF ((tcrit - tout)*h .LT. 0.0d0) go
to 625
1215 IF ((tn - tout)*h .LT. 0.0d0) go
to 245
1216 CALL
intdy(tout, 0, rwork(lyh), nyh, y, iflag)
1217 IF (iflag .NE. 0) go
to 627
1220 240 tcrit = rwork(1)
1221 IF ((tn - tcrit)*h .GT. 0.0d0) go
to 624
1222 245 hmx = dabs(tn) + dabs(h)
1223 ihit = dabs(tn - tcrit) .LE. 100.0d0*uround*hmx
1225 tnext = tn + h*(1.0d0 + 4.0d0*uround)
1226 IF ((tnext - tcrit)*h .LE. 0.0d0) go
to 250
1227 h = (tcrit - tn)*(1.0d0 - 4.0d0*uround)
1228 IF (istate .EQ. 2) jstart = -2
1241 IF ((nst-nslast) .GE. mxstep) go
to 500
1242 CALL
ewset(n, itol, rtol, atol, rwork(lyh), rwork(lewt))
1244 IF (rwork(i+lewt-1) .LE. 0.0d0) go
to 510
1245 260 rwork(i+lewt-1) = 1.0d0/rwork(i+lewt-1)
1246 270 tolsf = uround*
vnorm(n, rwork(lyh), rwork(lewt))
1247 IF (tolsf .LE. 1.0d0) go
to 280
1249 IF (nst .EQ. 0) go
to 626
1251 280
IF ((tn + h) .NE. tn) go
to 290
1253 IF (nhnil .GT. mxhnil) go
to 290
1254 CALL
xerrwd(
'LSODE-- WARNING..INTERNAL T (=R1) AND H (=R2) ARE',
1255 1 50, 101, 0, 0, 0, 0, 0, 0.0d0, 0.0d0)
1257 1
' SUCH THAT IN THE MACHINE, T + H = T ON THE NEXT STEP ',
1258 1 60, 101, 0, 0, 0, 0, 0, 0.0d0, 0.0d0)
1259 CALL
xerrwd(
' (H = STEP SIZE). SOLVER WILL CONTINUE ANYWAY',
1260 1 50, 101, 0, 0, 0, 0, 2, tn, h)
1261 IF (nhnil .LT. mxhnil) go
to 290
1262 CALL
xerrwd(
'LSODE-- ABOVE WARNING HAS BEEN ISSUED I1 TIMES. ',
1263 1 50, 102, 0, 0, 0, 0, 0, 0.0d0, 0.0d0)
1264 CALL
xerrwd(
' IT WILL NOT BE ISSUED AGAIN FOR THIS PROBLEM',
1265 1 50, 102, 0, 1, mxhnil, 0, 0, 0.0d0, 0.0d0)
1271 CALL
stode(neq, y, rwork(lyh), nyh, rwork(lyh), rwork(lewt),
1272 1 rwork(lsavf), rwork(lacor), rwork(lwm), iwork(liwm),
1274 IF (ierr .LT. 0)
THEN
1279 go
to(300, 530, 540), kgo
1286 go
to(310, 400, 330, 340, 350), itask
1288 310
IF ((tn - tout)*h .LT. 0.0d0) go
to 250
1289 CALL
intdy(tout, 0, rwork(lyh), nyh, y, iflag)
1293 330
IF ((tn - tout)*h .GE. 0.0d0) go
to 400
1296 340
IF ((tn - tout)*h .LT. 0.0d0) go
to 345
1297 CALL
intdy(tout, 0, rwork(lyh), nyh, y, iflag)
1300 345 hmx = dabs(tn) + dabs(h)
1301 ihit = dabs(tn - tcrit) .LE. 100.0d0*uround*hmx
1303 tnext = tn + h*(1.0d0 + 4.0d0*uround)
1304 IF ((tnext - tcrit)*h .LE. 0.0d0) go
to 250
1305 h = (tcrit - tn)*(1.0d0 - 4.0d0*uround)
1309 350 hmx = dabs(tn) + dabs(h)
1310 ihit = dabs(tn - tcrit) .LE. 100.0d0*uround*hmx
1321 410 y(i) = rwork(i+lyh-1)
1323 IF (itask .NE. 4 .AND. itask .NE. 5) go
to 420
1337 430 ntrep = ntrep + 1
1338 IF (ntrep .LT. 5)
RETURN
1340 1
'LSODE-- REPEATED CALLS WITH ISTATE = 1 AND TOUT = T (=R1) ',
1341 1 60, 301, 0, 0, 0, 0, 1, t, 0.0d0)
1353 500 CALL
xerrwd(
'LSODE-- AT CURRENT T (=R1), MXSTEP (=I1) STEPS ',
1354 1 50, 201, 0, 0, 0, 0, 0, 0.0d0, 0.0d0)
1355 CALL
xerrwd(
' TAKEN ON THIS CALL BEFORE REACHING TOUT ',
1356 1 50, 201, 0, 1, mxstep, 0, 1, tn, 0.0d0)
1360 510 ewti = rwork(lewt+i-1)
1361 CALL
xerrwd(.LE.
'LSODE-- AT T (=R1), EWT(I1) HAS BECOME R2 0.',
1362 1 50, 202, 0, 1, i, 0, 2, tn, ewti)
1366 520 CALL
xerrwd(
'LSODE-- AT T (=R1), TOO MUCH ACCURACY REQUESTED ',
1367 1 50, 203, 0, 0, 0, 0, 0, 0.0d0, 0.0d0)
1368 CALL
xerrwd(
' FOR PRECISION OF MACHINE.. SEE TOLSF (=R2) ',
1369 1 50, 203, 0, 0, 0, 0, 2, tn, tolsf)
1374 530 CALL
xerrwd(
'LSODE-- AT T(=R1) AND STEP SIZE H(=R2), THE ERROR',
1375 1 50, 204, 0, 0, 0, 0, 0, 0.0d0, 0.0d0)
1376 CALL
xerrwd(
' TEST FAILED REPEATEDLY OR WITH ABS(H) = HMIN',
1377 1 50, 204, 0, 0, 0, 0, 2, tn, h)
1381 540 CALL
xerrwd(
'LSODE-- AT T (=R1) AND STEP SIZE H (=R2), THE ',
1382 1 50, 205, 0, 0, 0, 0, 0, 0.0d0, 0.0d0)
1383 CALL
xerrwd(
' CORRECTOR CONVERGENCE FAILED REPEATEDLY ',
1384 1 50, 205, 0, 0, 0, 0, 0, 0.0d0, 0.0d0)
1385 CALL
xerrwd(
' OR WITH ABS(H) = HMIN ',
1386 1 30, 205, 0, 0, 0, 0, 2, tn, h)
1392 SIZE = dabs(rwork(i+lacor-1)*rwork(i+lewt-1))
1393 IF (big .GE.
size) go
to 570
1400 590 y(i) = rwork(i+lyh-1)
1420 601 CALL
xerrwd(
'LSODE-- ISTATE (=I1) ILLEGAL ',
1421 1 30, 1, 0, 1, istate, 0, 0, 0.0d0, 0.0d0)
1423 602 CALL
xerrwd(
'LSODE-- ITASK (=I1) ILLEGAL ',
1424 1 30, 2, 0, 1, itask, 0, 0, 0.0d0, 0.0d0)
1426 603 CALL
xerrwd(.GT.
'LSODE-- ISTATE 1 BUT LSODE NOT INITIALIZED ',
1427 1 50, 3, 0, 0, 0, 0, 0, 0.0d0, 0.0d0)
1429 604 CALL
xerrwd(.LT.
'LSODE-- NEQ (=I1) 1 ',
1430 1 30, 4, 0, 1, neq(1), 0, 0, 0.0d0, 0.0d0)
1432 605 CALL
xerrwd(
'LSODE-- ISTATE = 3 AND NEQ INCREASED (I1 TO I2) ',
1433 1 50, 5, 0, 2, n, neq(1), 0, 0.0d0, 0.0d0)
1435 606 CALL
xerrwd(
'LSODE-- ITOL (=I1) ILLEGAL ',
1436 1 30, 6, 0, 1, itol, 0, 0, 0.0d0, 0.0d0)
1438 607 CALL
xerrwd(
'LSODE-- IOPT (=I1) ILLEGAL ',
1439 1 30, 7, 0, 1, iopt, 0, 0, 0.0d0, 0.0d0)
1441 608 CALL
xerrwd(
'LSODE-- MF (=I1) ILLEGAL ',
1442 1 30, 8, 0, 1, mf, 0, 0, 0.0d0, 0.0d0)
1444 609 CALL
xerrwd(.LT..GE.
'LSODE-- ML (=I1) ILLEGAL.. 0 OR NEQ (=I2)',
1445 1 50, 9, 0, 2, ml, neq(1), 0, 0.0d0, 0.0d0)
1447 610 CALL
xerrwd(.LT..GE.
'LSODE-- MU (=I1) ILLEGAL.. 0 OR NEQ (=I2)',
1448 1 50, 10, 0, 2, mu, neq(1), 0, 0.0d0, 0.0d0)
1450 611 CALL
xerrwd(.LT.
'LSODE-- MAXORD (=I1) 0 ',
1451 1 30, 11, 0, 1, maxord, 0, 0, 0.0d0, 0.0d0)
1453 612 CALL
xerrwd(.LT.
'LSODE-- MXSTEP (=I1) 0 ',
1454 1 30, 12, 0, 1, mxstep, 0, 0, 0.0d0, 0.0d0)
1456 613 CALL
xerrwd(.LT.
'LSODE-- MXHNIL (=I1) 0 ',
1457 1 30, 13, 0, 1, mxhnil, 0, 0, 0.0d0, 0.0d0)
1459 614 CALL
xerrwd(
'LSODE-- TOUT (=R1) BEHIND T (=R2) ',
1460 1 40, 14, 0, 0, 0, 0, 2, tout, t)
1461 CALL
xerrwd(
' INTEGRATION DIRECTION IS GIVEN BY H0 (=R1) ',
1462 1 50, 14, 0, 0, 0, 0, 1, h0, 0.0d0)
1464 615 CALL
xerrwd(.LT.
'LSODE-- HMAX (=R1) 0.0 ',
1465 1 30, 15, 0, 0, 0, 0, 1, hmax, 0.0d0)
1467 616 CALL
xerrwd(.LT.
'LSODE-- HMIN (=R1) 0.0 ',
1468 1 30, 16, 0, 0, 0, 0, 1, hmin, 0.0d0)
1471 1
'LSODE-- RWORK LENGTH NEEDED, LENRW (=I1), EXCEEDS LRW (=I2)',
1472 1 60, 17, 0, 2, lenrw, lrw, 0, 0.0d0, 0.0d0)
1475 1
'LSODE-- IWORK LENGTH NEEDED, LENIW (=I1), EXCEEDS LIW (=I2)',
1476 1 60, 18, 0, 2, leniw, liw, 0, 0.0d0, 0.0d0)
1478 619 CALL
xerrwd(.LT.
'LSODE-- RTOL(I1) IS R1 0.0 ',
1479 1 40, 19, 0, 1, i, 0, 1, rtoli, 0.0d0)
1481 620 CALL
xerrwd(.LT.
'LSODE-- ATOL(I1) IS R1 0.0 ',
1482 1 40, 20, 0, 1, i, 0, 1, atoli, 0.0d0)
1484 621 ewti = rwork(lewt+i-1)
1485 CALL
xerrwd(.LE.
'LSODE-- EWT(I1) IS R1 0.0 ',
1486 1 40, 21, 0, 1, i, 0, 1, ewti, 0.0d0)
1489 1
'LSODE-- TOUT (=R1) TOO CLOSE TO T(=R2) TO START INTEGRATION',
1490 1 60, 22, 0, 0, 0, 0, 2, tout, t)
1493 1
'LSODE-- ITASK = I1 AND TOUT (=R1) BEHIND TCUR - HU (= R2) ',
1494 1 60, 23, 0, 1, itask, 0, 2, tout, tp)
1497 1
'LSODE-- ITASK = 4 OR 5 AND TCRIT (=R1) BEHIND TCUR (=R2) ',
1498 1 60, 24, 0, 0, 0, 0, 2, tcrit, tn)
1501 1
'LSODE-- ITASK = 4 OR 5 AND TCRIT (=R1) BEHIND TOUT (=R2) ',
1502 1 60, 25, 0, 0, 0, 0, 2, tcrit, tout)
1504 626 CALL
xerrwd(
'LSODE-- AT START OF PROBLEM, TOO MUCH ACCURACY ',
1505 1 50, 26, 0, 0, 0, 0, 0, 0.0d0, 0.0d0)
1507 1
' REQUESTED FOR PRECISION OF MACHINE.. SEE TOLSF (=R1) ',
1508 1 60, 26, 0, 0, 0, 0, 1, tolsf, 0.0d0)
1511 627 CALL
xerrwd(
'LSODE-- TROUBLE FROM INTDY. ITASK = I1, TOUT = R1',
1512 1 50, 27, 0, 1, itask, 0, 1, tout, 0.0d0)
1514 700
IF (illin .EQ. 5) go
to 710
1518 710 CALL
xerrwd(
'LSODE-- REPEATED OCCURRENCES OF ILLEGAL INPUT ',
1519 1 50, 302, 0, 0, 0, 0, 0, 0.0d0, 0.0d0)
1521 800 CALL
xerrwd(
'LSODE-- RUN ABORTED.. APPARENT INFINITE LOOP ',
1522 1 50, 303, 2, 0, 0, 0, 0, 0.0d0, 0.0d0)
subroutine prepj(NEQ, Y, YH, NYH, EWT, FTEM, SAVF, WM, IWM, F, JAC, IERR)
subroutine dlsode(F, NEQ, Y, T, TOUT, ITOL, RTOL, ATOL, ITASK, ISTATE, IOPT, RWORK, LRW, IWORK, LIW, JAC, MF)
subroutine xerrwd(MSG, NMES, NERR, LEVEL, NI, I1, I2, NR, R1, R2)
F77_RET_T F77_REAL &F77_RET_T F77_DBLE &F77_RET_T F77_REAL &F77_RET_T F77_DBLE &F77_RET_T F77_REAL &F77_RET_T F77_DBLE &F77_RET_T const F77_REAL const F77_REAL F77_REAL &F77_RET_T const F77_DBLE const F77_DBLE F77_DBLE &F77_RET_T F77_REAL &F77_RET_T F77_DBLE &F77_RET_T F77_DBLE &F77_RET_T F77_REAL &F77_RET_T F77_REAL &F77_RET_T F77_DBLE &F77_RET_T const F77_DBLE F77_DBLE &F77_RET_T const F77_REAL F77_REAL &F77_RET_T F77_REAL F77_REAL &F77_RET_T F77_DBLE F77_DBLE &F77_RET_T const F77_DBLE const F77_DBLE * f
subroutine stode(NEQ, Y, YH, NYH, YH1, EWT, SAVF, ACOR, WM, IWM, F, JAC, PJAC, SLVS, IERR)
may be zero for pure relative error test tem the relative tolerance must be greater than or equal to
OCTAVE_EXPORT octave_value_list etc The functions then dimension(columns)
subroutine ewset(N, ITOL, RTOL, ATOL, YCUR, EWT)
subroutine intdy(T, K, YH, NYH, DKY, IFLAG)
subroutine solsy(WM, IWM, X, TEM)
double precision function vnorm(N, V, W)