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, iowns
949 INTEGER ICF, IERPJ, IERSL, JCUR, JSTART, KFLAG, L, METH, MITER,
950 1 maxord, maxcor, msbp, mxncf, n, nq, nst, nfe, nje, nqu
951 INTEGER I, I1, I2, IFLAG, IMXER, KGO, LF0,
952 1 leniw, lenrw, lenwm, ml, mord, mu, mxhnl0, mxstp0
953 DOUBLE PRECISION ROWNS,
954 1 ccmax, el0, h, hmin, hmxi, hu, rc, tn, uround
955 DOUBLE PRECISION ATOLI, AYI, BIG, EWTI, H0, HMAX, HMX, RH, RTOLI,
956 1 tcrit, tdist, tnext, tol, tolsf, tp,
SIZE, sum, w0,
974 COMMON /ls0001/ rowns(209),
975 1 ccmax, el0, h, hmin, hmxi, hu, rc, tn, uround,
976 2 illin, init, lyh, lewt, lacor, lsavf, lwm, liwm,
977 3 mxstep, mxhnil, nhnil, ntrep, nslast, nyh, iowns(6),
978 4 icf, ierpj, iersl, jcur, jstart, kflag, l, meth, miter,
979 5 maxord, maxcor, msbp, mxncf, n, nq, nst, nfe, nje, nqu
981 DATA mord(1),mord(2)/12,5/, mxstp0/500/, mxhnl0/10/
990 IF (istate .LT. 1 .OR. istate .GT. 3) go to 601
991 IF (itask .LT. 1 .OR. itask .GT. 5) go to 602
992 IF (istate .EQ. 1) go to 10
993 IF (init .EQ. 0) go to 603
994 IF (istate .EQ. 2) go to 200
997 IF (tout .EQ. t) go to 430
1008 IF (neq(1) .LE. 0) go to 604
1009 IF (istate .EQ. 1) go to 25
1010 IF (neq(1) .GT. n) go to 605
1012 IF (itol .LT. 1 .OR. itol .GT. 4) go to 606
1013 IF (iopt .LT. 0 .OR. iopt .GT. 1) go to 607
1015 miter = mf - 10*meth
1016 IF (meth .LT. 1 .OR. meth .GT. 2) go to 608
1017 IF (miter .LT. 0 .OR. miter .GT. 5) go to 608
1018 IF (miter .LE. 3) go to 30
1021 IF (ml .LT. 0 .OR. ml .GE. n) go to 609
1022 IF (mu .LT. 0 .OR. mu .GE. n) go to 610
1025 IF (iopt .EQ. 1) go to 40
1029 IF (istate .EQ. 1) h0 = 0.0d0
1033 40 maxord = iwork(5)
1034 IF (maxord .LT. 0) go to 611
1035 IF (maxord .EQ. 0) maxord = 100
1036 maxord = min0(maxord,mord(meth))
1038 IF (mxstep .LT. 0) go to 612
1039 IF (mxstep .EQ. 0) mxstep = mxstp0
1041 IF (mxhnil .LT. 0) go to 613
1042 IF (mxhnil .EQ. 0) mxhnil = mxhnl0
1043 IF (istate .NE. 1) go to 50
1045 IF ((tout - t)*h0 .LT. 0.0d0) go to 614
1047 IF (hmax .LT. 0.0d0) go to 615
1049 IF (hmax .GT. 0.0d0) hmxi = 1.0d0/hmax
1051 IF (hmin .LT. 0.0d0) go to 616
1059 IF (istate .EQ. 1) nyh = n
1060 lwm = lyh + (maxord + 1)*nyh
1061 IF (miter .EQ. 0) lenwm = 0
1062 IF (miter .EQ. 1 .OR. miter .EQ. 2) lenwm = n*n + 2
1063 IF (miter .EQ. 3) lenwm = n + 2
1064 IF (miter .GE. 4) lenwm = (2*ml + mu + 1)*n + 2
1068 lenrw = lacor + n - 1
1072 IF (miter .EQ. 0 .OR. miter .EQ. 3) leniw = 20
1074 IF (lenrw .GT. lrw) go to 617
1075 IF (leniw .GT. liw) go to 618
1080 IF (itol .GE. 3) rtoli = rtol(i)
1081 IF (itol .EQ. 2 .OR. itol .EQ. 4) atoli = atol(i)
1082 IF (rtoli .LT. 0.0d0) go to 619
1083 IF (atoli .LT. 0.0d0) go to 620
1085 IF (istate .EQ. 1) go to 100
1088 IF (nq .LE. maxord) go to 90
1091 80 rwork(i+lsavf-1) = rwork(i+lwm-1)
1093 90
IF (miter .GT. 0) rwork(lwm) = dsqrt(uround)
1094 IF (n .EQ. nyh) go to 200
1097 i2 = lyh + (maxord + 1)*nyh - 1
1098 IF (i1 .GT. i2) go to 200
1109 100 uround = d1mach(4)
1111 IF (itask .NE. 4 .AND. itask .NE. 5) go to 110
1113 IF ((tcrit - tout)*(tout - t) .LT. 0.0d0) go to 625
1114 IF (h0 .NE. 0.0d0 .AND. (t + h0 - tcrit)*h0 .GT. 0.0d0)
1117 IF (miter .GT. 0) rwork(lwm) = dsqrt(uround)
1131 CALL
f(neq, t, y, rwork(lf0), ierr)
1132 IF (ierr .LT. 0)
THEN
1139 115 rwork(i+lyh-1) = y(i)
1143 CALL
ewset(n, itol, rtol, atol, rwork(lyh), rwork(lewt))
1145 IF (rwork(i+lewt-1) .LE. 0.0d0) go to 621
1146 120 rwork(i+lewt-1) = 1.0d0/rwork(i+lewt-1)
1163 IF (h0 .NE. 0.0d0) go to 180
1164 tdist = dabs(tout - t)
1165 w0 = dmax1(dabs(t),dabs(tout))
1166 IF (tdist .LT. 2.0d0*uround*w0) go to 622
1168 IF (itol .LE. 2) go to 140
1170 130 tol = dmax1(tol,rtol(i))
1171 140
IF (tol .GT. 0.0d0) go to 160
1174 IF (itol .EQ. 2 .OR. itol .EQ. 4) atoli = atol(i)
1176 IF (ayi .NE. 0.0d0) tol = dmax1(tol,atoli/ayi)
1178 160 tol = dmax1(tol,100.0d0*uround)
1179 tol = dmin1(tol,0.001d0)
1180 sum =
vnorm(n, rwork(lf0), rwork(lewt))
1181 sum = 1.0d0/(tol*w0*w0) + tol*sum**2
1182 h0 = 1.0d0/dsqrt(sum)
1183 h0 = dmin1(h0,tdist)
1184 h0 = dsign(h0,tout-t)
1186 180 rh = dabs(h0)*hmxi
1187 IF (rh .GT. 1.0d0) h0 = h0/rh
1191 190 rwork(i+lf0-1) = h0*rwork(i+lf0-1)
1199 go to(210, 250, 220, 230, 240), itask
1200 210
IF ((tn - tout)*h .LT. 0.0d0) go to 250
1201 CALL
intdy(tout, 0, rwork(lyh), nyh, y, iflag)
1202 IF (iflag .NE. 0) go to 627
1205 220 tp = tn - hu*(1.0d0 + 100.0d0*uround)
1206 IF ((tp - tout)*h .GT. 0.0d0) go to 623
1207 IF ((tn - tout)*h .LT. 0.0d0) go to 250
1209 230 tcrit = rwork(1)
1210 IF ((tn - tcrit)*h .GT. 0.0d0) go to 624
1211 IF ((tcrit - tout)*h .LT. 0.0d0) go to 625
1212 IF ((tn - tout)*h .LT. 0.0d0) go to 245
1213 CALL
intdy(tout, 0, rwork(lyh), nyh, y, iflag)
1214 IF (iflag .NE. 0) go to 627
1217 240 tcrit = rwork(1)
1218 IF ((tn - tcrit)*h .GT. 0.0d0) go to 624
1219 245 hmx = dabs(tn) + dabs(h)
1220 ihit = dabs(tn - tcrit) .LE. 100.0d0*uround*hmx
1222 tnext = tn + h*(1.0d0 + 4.0d0*uround)
1223 IF ((tnext - tcrit)*h .LE. 0.0d0) go to 250
1224 h = (tcrit - tn)*(1.0d0 - 4.0d0*uround)
1225 IF (istate .EQ. 2) jstart = -2
1238 IF ((nst-nslast) .GE. mxstep) go to 500
1239 CALL
ewset(n, itol, rtol, atol, rwork(lyh), rwork(lewt))
1241 IF (rwork(i+lewt-1) .LE. 0.0d0) go to 510
1242 260 rwork(i+lewt-1) = 1.0d0/rwork(i+lewt-1)
1243 270 tolsf = uround*
vnorm(n, rwork(lyh), rwork(lewt))
1244 IF (tolsf .LE. 1.0d0) go to 280
1246 IF (nst .EQ. 0) go to 626
1248 280
IF ((tn + h) .NE. tn) go to 290
1250 IF (nhnil .GT. mxhnil) go to 290
1251 CALL
xerrwd(
'LSODE-- WARNING..INTERNAL T (=R1) AND H (=R2) ARE',
1252 1 50, 101, 0, 0, 0, 0, 0, 0.0d0, 0.0d0)
1254 1
' SUCH THAT IN THE MACHINE, T + H = T ON THE NEXT STEP ',
1255 1 60, 101, 0, 0, 0, 0, 0, 0.0d0, 0.0d0)
1256 CALL
xerrwd(
' (H = STEP SIZE). SOLVER WILL CONTINUE ANYWAY',
1257 1 50, 101, 0, 0, 0, 0, 2, tn, h)
1258 IF (nhnil .LT. mxhnil) go to 290
1259 CALL
xerrwd(
'LSODE-- ABOVE WARNING HAS BEEN ISSUED I1 TIMES. ',
1260 1 50, 102, 0, 0, 0, 0, 0, 0.0d0, 0.0d0)
1261 CALL
xerrwd(
' IT WILL NOT BE ISSUED AGAIN FOR THIS PROBLEM',
1262 1 50, 102, 0, 1, mxhnil, 0, 0, 0.0d0, 0.0d0)
1268 CALL
stode(neq, y, rwork(lyh), nyh, rwork(lyh), rwork(lewt),
1269 1 rwork(lsavf), rwork(lacor), rwork(lwm), iwork(liwm),
1271 IF (ierr .LT. 0)
THEN
1276 go to(300, 530, 540), kgo
1283 go to(310, 400, 330, 340, 350), itask
1285 310
IF ((tn - tout)*h .LT. 0.0d0) go to 250
1286 CALL
intdy(tout, 0, rwork(lyh), nyh, y, iflag)
1290 330
IF ((tn - tout)*h .GE. 0.0d0) go to 400
1293 340
IF ((tn - tout)*h .LT. 0.0d0) go to 345
1294 CALL
intdy(tout, 0, rwork(lyh), nyh, y, iflag)
1297 345 hmx = dabs(tn) + dabs(h)
1298 ihit = dabs(tn - tcrit) .LE. 100.0d0*uround*hmx
1300 tnext = tn + h*(1.0d0 + 4.0d0*uround)
1301 IF ((tnext - tcrit)*h .LE. 0.0d0) go to 250
1302 h = (tcrit - tn)*(1.0d0 - 4.0d0*uround)
1306 350 hmx = dabs(tn) + dabs(h)
1307 ihit = dabs(tn - tcrit) .LE. 100.0d0*uround*hmx
1318 410 y(i) = rwork(i+lyh-1)
1320 IF (itask .NE. 4 .AND. itask .NE. 5) go to 420
1334 430 ntrep = ntrep + 1
1335 IF (ntrep .LT. 5)
RETURN
1337 1
'LSODE-- REPEATED CALLS WITH ISTATE = 1 AND TOUT = T (=R1) ',
1338 1 60, 301, 0, 0, 0, 0, 1, t, 0.0d0)
1350 500 CALL
xerrwd(
'LSODE-- AT CURRENT T (=R1), MXSTEP (=I1) STEPS ',
1351 1 50, 201, 0, 0, 0, 0, 0, 0.0d0, 0.0d0)
1352 CALL
xerrwd(
' TAKEN ON THIS CALL BEFORE REACHING TOUT ',
1353 1 50, 201, 0, 1, mxstep, 0, 1, tn, 0.0d0)
1357 510 ewti = rwork(lewt+i-1)
1358 CALL
xerrwd(.LE.
'LSODE-- AT T (=R1), EWT(I1) HAS BECOME R2 0.',
1359 1 50, 202, 0, 1, i, 0, 2, tn, ewti)
1363 520 CALL
xerrwd(
'LSODE-- AT T (=R1), TOO MUCH ACCURACY REQUESTED ',
1364 1 50, 203, 0, 0, 0, 0, 0, 0.0d0, 0.0d0)
1365 CALL
xerrwd(
' FOR PRECISION OF MACHINE.. SEE TOLSF (=R2) ',
1366 1 50, 203, 0, 0, 0, 0, 2, tn, tolsf)
1371 530 CALL
xerrwd(
'LSODE-- AT T(=R1) AND STEP SIZE H(=R2), THE ERROR',
1372 1 50, 204, 0, 0, 0, 0, 0, 0.0d0, 0.0d0)
1373 CALL
xerrwd(
' TEST FAILED REPEATEDLY OR WITH ABS(H) = HMIN',
1374 1 50, 204, 0, 0, 0, 0, 2, tn, h)
1378 540 CALL
xerrwd(
'LSODE-- AT T (=R1) AND STEP SIZE H (=R2), THE ',
1379 1 50, 205, 0, 0, 0, 0, 0, 0.0d0, 0.0d0)
1380 CALL
xerrwd(
' CORRECTOR CONVERGENCE FAILED REPEATEDLY ',
1381 1 50, 205, 0, 0, 0, 0, 0, 0.0d0, 0.0d0)
1382 CALL
xerrwd(
' OR WITH ABS(H) = HMIN ',
1383 1 30, 205, 0, 0, 0, 0, 2, tn, h)
1389 SIZE = dabs(rwork(i+lacor-1)*rwork(i+lewt-1))
1390 IF (big .GE.
size) go to 570
1397 590 y(i) = rwork(i+lyh-1)
1417 601 CALL
xerrwd(
'LSODE-- ISTATE (=I1) ILLEGAL ',
1418 1 30, 1, 0, 1, istate, 0, 0, 0.0d0, 0.0d0)
1420 602 CALL
xerrwd(
'LSODE-- ITASK (=I1) ILLEGAL ',
1421 1 30, 2, 0, 1, itask, 0, 0, 0.0d0, 0.0d0)
1423 603 CALL
xerrwd(.GT.
'LSODE-- ISTATE 1 BUT LSODE NOT INITIALIZED ',
1424 1 50, 3, 0, 0, 0, 0, 0, 0.0d0, 0.0d0)
1426 604 CALL
xerrwd(.LT.
'LSODE-- NEQ (=I1) 1 ',
1427 1 30, 4, 0, 1, neq(1), 0, 0, 0.0d0, 0.0d0)
1429 605 CALL
xerrwd(
'LSODE-- ISTATE = 3 AND NEQ INCREASED (I1 TO I2) ',
1430 1 50, 5, 0, 2, n, neq(1), 0, 0.0d0, 0.0d0)
1432 606 CALL
xerrwd(
'LSODE-- ITOL (=I1) ILLEGAL ',
1433 1 30, 6, 0, 1, itol, 0, 0, 0.0d0, 0.0d0)
1435 607 CALL
xerrwd(
'LSODE-- IOPT (=I1) ILLEGAL ',
1436 1 30, 7, 0, 1, iopt, 0, 0, 0.0d0, 0.0d0)
1438 608 CALL
xerrwd(
'LSODE-- MF (=I1) ILLEGAL ',
1439 1 30, 8, 0, 1, mf, 0, 0, 0.0d0, 0.0d0)
1441 609 CALL
xerrwd(.LT..GE.
'LSODE-- ML (=I1) ILLEGAL.. 0 OR NEQ (=I2)',
1442 1 50, 9, 0, 2, ml, neq(1), 0, 0.0d0, 0.0d0)
1444 610 CALL
xerrwd(.LT..GE.
'LSODE-- MU (=I1) ILLEGAL.. 0 OR NEQ (=I2)',
1445 1 50, 10, 0, 2, mu, neq(1), 0, 0.0d0, 0.0d0)
1447 611 CALL
xerrwd(.LT.
'LSODE-- MAXORD (=I1) 0 ',
1448 1 30, 11, 0, 1, maxord, 0, 0, 0.0d0, 0.0d0)
1450 612 CALL
xerrwd(.LT.
'LSODE-- MXSTEP (=I1) 0 ',
1451 1 30, 12, 0, 1, mxstep, 0, 0, 0.0d0, 0.0d0)
1453 613 CALL
xerrwd(.LT.
'LSODE-- MXHNIL (=I1) 0 ',
1454 1 30, 13, 0, 1, mxhnil, 0, 0, 0.0d0, 0.0d0)
1456 614 CALL
xerrwd(
'LSODE-- TOUT (=R1) BEHIND T (=R2) ',
1457 1 40, 14, 0, 0, 0, 0, 2, tout, t)
1458 CALL
xerrwd(
' INTEGRATION DIRECTION IS GIVEN BY H0 (=R1) ',
1459 1 50, 14, 0, 0, 0, 0, 1, h0, 0.0d0)
1461 615 CALL
xerrwd(.LT.
'LSODE-- HMAX (=R1) 0.0 ',
1462 1 30, 15, 0, 0, 0, 0, 1, hmax, 0.0d0)
1464 616 CALL
xerrwd(.LT.
'LSODE-- HMIN (=R1) 0.0 ',
1465 1 30, 16, 0, 0, 0, 0, 1, hmin, 0.0d0)
1468 1
'LSODE-- RWORK LENGTH NEEDED, LENRW (=I1), EXCEEDS LRW (=I2)',
1469 1 60, 17, 0, 2, lenrw, lrw, 0, 0.0d0, 0.0d0)
1472 1
'LSODE-- IWORK LENGTH NEEDED, LENIW (=I1), EXCEEDS LIW (=I2)',
1473 1 60, 18, 0, 2, leniw, liw, 0, 0.0d0, 0.0d0)
1475 619 CALL
xerrwd(.LT.
'LSODE-- RTOL(I1) IS R1 0.0 ',
1476 1 40, 19, 0, 1, i, 0, 1, rtoli, 0.0d0)
1478 620 CALL
xerrwd(.LT.
'LSODE-- ATOL(I1) IS R1 0.0 ',
1479 1 40, 20, 0, 1, i, 0, 1, atoli, 0.0d0)
1481 621 ewti = rwork(lewt+i-1)
1482 CALL
xerrwd(.LE.
'LSODE-- EWT(I1) IS R1 0.0 ',
1483 1 40, 21, 0, 1, i, 0, 1, ewti, 0.0d0)
1486 1
'LSODE-- TOUT (=R1) TOO CLOSE TO T(=R2) TO START INTEGRATION',
1487 1 60, 22, 0, 0, 0, 0, 2, tout, t)
1490 1
'LSODE-- ITASK = I1 AND TOUT (=R1) BEHIND TCUR - HU (= R2) ',
1491 1 60, 23, 0, 1, itask, 0, 2, tout, tp)
1494 1
'LSODE-- ITASK = 4 OR 5 AND TCRIT (=R1) BEHIND TCUR (=R2) ',
1495 1 60, 24, 0, 0, 0, 0, 2, tcrit, tn)
1498 1
'LSODE-- ITASK = 4 OR 5 AND TCRIT (=R1) BEHIND TOUT (=R2) ',
1499 1 60, 25, 0, 0, 0, 0, 2, tcrit, tout)
1501 626 CALL
xerrwd(
'LSODE-- AT START OF PROBLEM, TOO MUCH ACCURACY ',
1502 1 50, 26, 0, 0, 0, 0, 0, 0.0d0, 0.0d0)
1504 1
' REQUESTED FOR PRECISION OF MACHINE.. SEE TOLSF (=R1) ',
1505 1 60, 26, 0, 0, 0, 0, 1, tolsf, 0.0d0)
1508 627 CALL
xerrwd(
'LSODE-- TROUBLE FROM INTDY. ITASK = I1, TOUT = R1',
1509 1 50, 27, 0, 1, itask, 0, 1, tout, 0.0d0)
1511 700
IF (illin .EQ. 5) go to 710
1515 710 CALL
xerrwd(
'LSODE-- REPEATED OCCURRENCES OF ILLEGAL INPUT ',
1516 1 50, 302, 0, 0, 0, 0, 0, 0.0d0, 0.0d0)
1518 800 CALL
xerrwd(
'LSODE-- RUN ABORTED.. APPARENT INFINITE LOOP ',
1519 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)
std::string dimension(void) const
subroutine xerrwd(MSG, NMES, NERR, LEVEL, NI, I1, I2, NR, R1, R2)
subroutine stode(NEQ, Y, YH, NYH, YH1, EWT, SAVF, ACOR, WM, IWM, F, JAC, PJAC, SLVS, IERR)
F77_RET_T const double const double * f
size_t size(T const (&)[z])
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)