1 SUBROUTINE ddasrt (RES,NEQ,T,Y,YPRIME,TOUT,
2 * info,rtol,atol,idid,rwork,lrw,iwork,liw,rpar,ipar,jac,
867 IMPLICIT DOUBLE PRECISION(a-h,o-z)
878 parameter(lml=1, lmu=2, lmxord=3, lmtype=4, lnst=11,
879 * lnre=12, lnje=13, letf=14, lctf=15, lnge=16, lnpd=17,
880 * lirfnd=18, lmxstp=21, lipvt=22, ljcalc=5, lphase=6, lk=7,
881 * lkold=8, lns=9, lnstl=10, liwm=1)
887 parameter(ltstop=1, lhmax=2, lh=3, ltn=4,
888 * lcj=5, lcjold=6, lhold=7, ls=8, lround=9,
889 * lalpha=11, lbeta=17,
lgamma=23,
890 * lpsi=29, lsigma=35, lt0=41, ltlast=42, lalphr=43, lx2=44,
894 IF(info(1).NE.0)go to 100
904 IF(info(i).NE.0.AND.info(i).NE.1)go to 701
907 IF(neq.LE.0)go to 702
911 IF(info(9).EQ.0)go to 20
913 IF(mxord.LT.1.OR.mxord.GT.5)go to 703
914 20 iwork(lmxord)=mxord
917 IF(info(6).NE.0)go to 40
919 lenrw=50+(iwork(lmxord)+4)*neq+lenpd+3*ng
920 IF(info(5).NE.0)go to 30
925 40
IF(iwork(lml).LT.0.OR.iwork(lml).GE.neq)go to 717
926 IF(iwork(lmu).LT.0.OR.iwork(lmu).GE.neq)go to 718
927 lenpd=(2*iwork(lml)+iwork(lmu)+1)*neq
928 IF(info(5).NE.0)go to 50
930 mband=iwork(lml)+iwork(lmu)+1
932 lenrw=50+(iwork(lmxord)+4)*neq+lenpd+2*msave+3*ng
935 lenrw=50+(iwork(lmxord)+4)*neq+lenpd+3*ng
940 IF(lrw.LT.lenrw)go to 704
941 IF(liw.LT.leniw)go to 705
945 IF(tout .EQ. t)go to 719
946 IF(ng .LT. 0) go to 730
949 IF(info(7).EQ.0)go to 70
951 IF(hmax.LE.0.0d0)go to 710
956 IF(info(12).EQ.0)go to 80
958 IF(mxstp.LT.0)go to 716
959 80 iwork(lmxstp)=mxstp
979 IF(info(1).EQ.1)go to 110
980 IF(info(1).NE.-1)go to 701
985 msg =
'DASRT-- THE LAST STEP TERMINATED WITH A NEGATIVE'
986 CALL
xerrwd(msg,49,201,0,0,0,0,0,0.0d0,0.0d0)
987 msg =
'DASRT-- VALUE (=I1) OF IDID AND NO APPROPRIATE'
988 CALL
xerrwd(msg,47,202,0,1,idid,0,0,0.0d0,0.0d0)
989 msg =
'DASRT-- ACTION WAS TAKEN. RUN TERMINATED'
990 CALL
xerrwd(msg,41,203,1,0,0,0,0,0.0d0,0.0d0)
993 iwork(lnstl)=iwork(lnst)
1008 IF(info(2).EQ.1)rtoli=rtol(i)
1009 IF(info(2).EQ.1)atoli=atol(i)
1010 IF(rtoli.GT.0.0d0.OR.atoli.GT.0.0d0)nzflg=1
1011 IF(rtoli.LT.0.0d0)go to 706
1012 IF(atoli.LT.0.0d0)go to 707
1014 IF(nzflg.EQ.0)go to 708
1024 lpd=lphi+(iwork(lmxord)+1)*neq
1026 ntemp=npd+iwork(lnpd)
1027 IF(info(1).EQ.1)go to 400
1041 CALL
ddawts(neq,info(2),rtol,atol,y,rwork(lwt),rpar,ipar)
1043 IF(rwork(lwt+i-1).LE.0.0d0) go to 713
1048 rwork(lround) = uround
1049 hmin = 4.0d0*uround*dmax1(dabs(t),dabs(tout))
1052 tdist = dabs(tout - t)
1053 IF(tdist .LT. hmin) go to 714
1056 IF (info(8) .EQ. 0) go to 310
1058 IF ((tout - t)*ho .LT. 0.0d0) go to 711
1059 IF (ho .EQ. 0.0d0) go to 712
1066 ypnorm =
ddanrm(neq,yprime,rwork(lwt),rpar,ipar)
1067 IF (ypnorm .GT. 0.5d0/ho) ho = 0.5d0/ypnorm
1068 ho = dsign(ho,tout-t)
1070 320
IF (info(7) .EQ. 0) go to 330
1071 rh = dabs(ho)/rwork(lhmax)
1072 IF (rh .GT. 1.0d0) ho = ho/rh
1074 330
IF (info(4) .EQ. 0) go to 340
1075 tstop = rwork(ltstop)
1076 IF ((tstop - t)*ho .LT. 0.0d0) go to 715
1077 IF ((t + ho - tstop)*ho .GT. 0.0d0) ho = tstop - t
1078 IF ((tstop - tout)*ho .LT. 0.0d0) go to 709
1081 340
IF (info(11) .EQ. 0) go to 350
1082 CALL
ddaini(tn,y,yprime,neq,
1083 * res,jac,ho,rwork(lwt),idid,rpar,ipar,
1084 * rwork(lphi),rwork(ldelta),rwork(le),
1085 * rwork(lwm),iwork(liwm),hmin,rwork(lround),
1087 IF (idid .LT. 0) go to 390
1094 360 itemp = lphi + neq
1096 rwork(lphi + i - 1) = y(i)
1097 370 rwork(itemp + i - 1) = h*yprime(i)
1105 rwork(lpsi+1)=2.0d0*h
1107 IF(ng .EQ. 0) go to 390
1108 CALL
drchek(1,g,ng,neq,t,tout,y,rwork(le),rwork(lphi),
1109 * rwork(lpsi),iwork(lkold),rwork(lg0),rwork(lg1),
1110 * rwork(lgx),jroot,irt,rwork(lround),info(3),
1111 * rwork,iwork,rpar,ipar)
1112 IF(irt .NE. 0) go to 732
1117 IF(ng .EQ. 0 .OR. info(11) .EQ. 0) go to 390
1118 CALL
drchek(3,g,ng,neq,tn,tout,y,rwork(le),rwork(lphi),
1119 * rwork(lpsi),iwork(lkold),rwork(lg0),rwork(lg1),
1120 * rwork(lgx),jroot,irt,rwork(lround),info(3),
1121 * rwork,iwork,rpar,ipar)
1122 IF(irt .NE. 1) go to 390
1138 uround=rwork(lround)
1142 IF(ng .EQ. 0) go to 405
1146 CALL
drchek(2,g,ng,neq,tn,tout,y,rwork(le),rwork(lphi),
1147 * rwork(lpsi),iwork(lkold),rwork(lg0),rwork(lg1),
1148 * rwork(lgx),jroot,irt,rwork(lround),info(3),
1149 * rwork,iwork,rpar,ipar)
1150 IF(irt .NE. 1) go to 405
1158 IF(info(7) .EQ. 0) go to 410
1159 rh = dabs(h)/rwork(lhmax)
1160 IF(rh .GT. 1.0d0) h = h/rh
1162 IF(t .EQ. tout) go to 719
1163 IF((t - tout)*h .GT. 0.0d0) go to 711
1164 IF(info(4) .EQ. 1) go to 430
1165 IF(info(3) .EQ. 1) go to 420
1166 IF((tn-tout)*h.LT.0.0d0)go to 490
1167 CALL
ddatrp(tn,tout,y,yprime,neq,iwork(lkold),
1168 * rwork(lphi),rwork(lpsi))
1173 420
IF((tn-t)*h .LE. 0.0d0) go to 490
1174 IF((tn - tout)*h .GT. 0.0d0) go to 425
1175 CALL
ddatrp(tn,tn,y,yprime,neq,iwork(lkold),
1176 * rwork(lphi),rwork(lpsi))
1182 CALL
ddatrp(tn,tout,y,yprime,neq,iwork(lkold),
1183 * rwork(lphi),rwork(lpsi))
1188 430
IF(info(3) .EQ. 1) go to 440
1190 IF((tn-tstop)*h.GT.0.0d0) go to 715
1191 IF((tstop-tout)*h.LT.0.0d0)go to 709
1192 IF((tn-tout)*h.LT.0.0d0)go to 450
1193 CALL
ddatrp(tn,tout,y,yprime,neq,iwork(lkold),
1194 * rwork(lphi),rwork(lpsi))
1199 440 tstop = rwork(ltstop)
1200 IF((tn-tstop)*h .GT. 0.0d0) go to 715
1201 IF((tstop-tout)*h .LT. 0.0d0) go to 709
1202 IF((tn-t)*h .LE. 0.0d0) go to 450
1203 IF((tn - tout)*h .GT. 0.0d0) go to 445
1204 CALL
ddatrp(tn,tn,y,yprime,neq,iwork(lkold),
1205 * rwork(lphi),rwork(lpsi))
1211 CALL
ddatrp(tn,tout,y,yprime,neq,iwork(lkold),
1212 * rwork(lphi),rwork(lpsi))
1219 IF(dabs(tn-tstop).GT.100.0d0*uround*
1220 * (dabs(tn)+dabs(h)))go to 460
1221 CALL
ddatrp(tn,tstop,y,yprime,neq,iwork(lkold),
1222 * rwork(lphi),rwork(lpsi))
1228 IF((tnext-tstop)*h.LE.0.0d0)go to 490
1232 490
IF (done) go to 590
1246 IF (idid .EQ. -12) go to 527
1249 IF((iwork(lnst)-iwork(lnstl)).LT.iwork(lmxstp))
1255 510 CALL
ddawts(neq,info(2),rtol,atol,rwork(lphi),
1256 * rwork(lwt),rpar,ipar)
1258 IF(rwork(i+lwt-1).GT.0.0d0)go to 520
1264 r=
ddanrm(neq,rwork(lphi),rwork(lwt),rpar,ipar)*
1266 IF(r.LE.1.0d0)go to 525
1268 IF(info(2).EQ.1)go to 523
1275 524 atol(i)=r*atol(i)
1281 hmin=4.0d0*uround*dmax1(dabs(tn),dabs(tout))
1284 IF (info(7) .EQ. 0) go to 526
1285 rh =
abs(h)/rwork(lhmax)
1286 IF (rh .GT. 1.0d0) h = h/rh
1289 CALL
ddastp(tn,y,yprime,neq,
1290 * res,jac,h,rwork(lwt),info(1),idid,rpar,ipar,
1291 * rwork(lphi),rwork(ldelta),rwork(le),
1292 * rwork(lwm),iwork(liwm),
1293 * rwork(lalpha),rwork(lbeta),rwork(
lgamma),
1294 * rwork(lpsi),rwork(lsigma),
1295 * rwork(lcj),rwork(lcjold),rwork(lhold),
1296 * rwork(ls),hmin,rwork(lround),
1297 * iwork(lphase),iwork(ljcalc),iwork(lk),
1298 * iwork(lkold),iwork(lns),info(10),ntemp)
1299 527
IF(idid.LT.0)go to 600
1306 IF(ng .EQ. 0) go to 529
1310 CALL
drchek(3,g,ng,neq,tn,tout,y,rwork(le),rwork(lphi),
1311 * rwork(lpsi),iwork(lkold),rwork(lg0),rwork(lg1),
1312 * rwork(lgx),jroot,irt,rwork(lround),info(3),
1313 * rwork,iwork,rpar,ipar)
1314 IF(irt .NE. 1) go to 529
1321 IF(info(4).NE.0)go to 540
1322 IF(info(3).NE.0)go to 530
1323 IF((tn-tout)*h.LT.0.0d0)go to 500
1324 CALL
ddatrp(tn,tout,y,yprime,neq,
1325 * iwork(lkold),rwork(lphi),rwork(lpsi))
1329 530
IF((tn-tout)*h.GE.0.0d0)go to 535
1333 535 CALL
ddatrp(tn,tout,y,yprime,neq,
1334 * iwork(lkold),rwork(lphi),rwork(lpsi))
1338 540
IF(info(3).NE.0)go to 550
1339 IF((tn-tout)*h.LT.0.0d0)go to 542
1340 CALL
ddatrp(tn,tout,y,yprime,neq,
1341 * iwork(lkold),rwork(lphi),rwork(lpsi))
1345 542
IF(dabs(tn-tstop).LE.100.0d0*uround*
1346 * (dabs(tn)+dabs(h)))go to 545
1348 IF((tnext-tstop)*h.LE.0.0d0)go to 500
1351 545 CALL
ddatrp(tn,tstop,y,yprime,neq,
1352 * iwork(lkold),rwork(lphi),rwork(lpsi))
1356 550
IF((tn-tout)*h.GE.0.0d0)go to 555
1357 IF(dabs(tn-tstop).LE.100.0d0*uround*(dabs(tn)+dabs(h)))go to 552
1361 552 CALL
ddatrp(tn,tstop,y,yprime,neq,
1362 * iwork(lkold),rwork(lphi),rwork(lpsi))
1366 555 CALL
ddatrp(tn,tout,y,yprime,neq,
1367 * iwork(lkold),rwork(lphi),rwork(lpsi))
1390 go to(610,620,630,690,690,640,650,660,670,675,
1395 610 msg =
'DASRT-- AT CURRENT T (=R1) 500 STEPS'
1396 CALL
xerrwd(msg,38,610,0,0,0,0,1,tn,0.0d0)
1397 msg =
'DASRT-- TAKEN ON THIS CALL BEFORE REACHING TOUT'
1398 CALL
xerrwd(msg,48,611,0,0,0,0,0,0.0d0,0.0d0)
1402 620 msg =
'DASRT-- AT T (=R1) TOO MUCH ACCURACY REQUESTED'
1403 CALL
xerrwd(msg,47,620,0,0,0,0,1,tn,0.0d0)
1404 msg =
'DASRT-- FOR PRECISION OF MACHINE. RTOL AND ATOL'
1405 CALL
xerrwd(msg,48,621,0,0,0,0,0,0.0d0,0.0d0)
1406 msg =
'DASRT-- WERE INCREASED TO APPROPRIATE VALUES'
1407 CALL
xerrwd(msg,45,622,0,0,0,0,0,0.0d0,0.0d0)
1411 630 msg =
'DASRT-- AT T (=R1) SOME ELEMENT OF WT'
1412 CALL
xerrwd(msg,38,630,0,0,0,0,1,tn,0.0d0)
1413 msg = .LE.
'DASRT-- HAS BECOME 0.0'
1414 CALL
xerrwd(msg,28,631,0,0,0,0,0,0.0d0,0.0d0)
1418 640 msg =
'DASRT-- AT T (=R1) AND STEPSIZE H (=R2) THE'
1419 CALL
xerrwd(msg,44,640,0,0,0,0,2,tn,h)
1420 msg=
'DASRT-- ERROR TEST FAILED REPEATEDLY OR WITH ABS(H)=HMIN'
1421 CALL
xerrwd(msg,57,641,0,0,0,0,0,0.0d0,0.0d0)
1425 650 msg =
'DASRT-- AT T (=R1) AND STEPSIZE H (=R2) THE'
1426 CALL
xerrwd(msg,44,650,0,0,0,0,2,tn,h)
1427 msg =
'DASRT-- CORRECTOR FAILED TO CONVERGE REPEATEDLY'
1428 CALL
xerrwd(msg,48,651,0,0,0,0,0,0.0d0,0.0d0)
1429 msg =
'DASRT-- OR WITH ABS(H)=HMIN'
1430 CALL
xerrwd(msg,28,652,0,0,0,0,0,0.0d0,0.0d0)
1434 660 msg =
'DASRT-- AT T (=R1) AND STEPSIZE H (=R2) THE'
1435 CALL
xerrwd(msg,44,660,0,0,0,0,2,tn,h)
1436 msg =
'DASRT-- ITERATION MATRIX IS SINGULAR'
1437 CALL
xerrwd(msg,37,661,0,0,0,0,0,0.0d0,0.0d0)
1441 670 msg =
'DASRT-- AT T (=R1) AND STEPSIZE H (=R2) THE'
1442 CALL
xerrwd(msg,44,670,0,0,0,0,2,tn,h)
1443 msg =
'DASRT-- CORRECTOR COULD NOT CONVERGE. ALSO, THE'
1444 CALL
xerrwd(msg,49,671,0,0,0,0,0,0.0d0,0.0d0)
1445 msg =
'DASRT-- ERROR TEST FAILED REPEATEDLY.'
1446 CALL
xerrwd(msg,38,672,0,0,0,0,0,0.0d0,0.0d0)
1450 675 msg =
'DASRT-- AT T (=R1) AND STEPSIZE H (=R2) THE'
1451 CALL
xerrwd(msg,44,675,0,0,0,0,2,tn,h)
1452 msg =
'DASRT-- CORRECTOR COULD NOT CONVERGE BECAUSE'
1453 CALL
xerrwd(msg,45,676,0,0,0,0,0,0.0d0,0.0d0)
1454 msg =
'DASRT-- IRES WAS EQUAL TO MINUS ONE'
1455 CALL
xerrwd(msg,36,677,0,0,0,0,0,0.0d0,0.0d0)
1459 680 msg =
'DASRT-- AT T (=R1) AND STEPSIZE H (=R2)'
1460 CALL
xerrwd(msg,40,680,0,0,0,0,2,tn,h)
1461 msg =
'DASRT-- IRES WAS EQUAL TO MINUS TWO'
1462 CALL
xerrwd(msg,36,681,0,0,0,0,0,0.0d0,0.0d0)
1466 685 msg =
'DASRT-- AT T (=R1) AND STEPSIZE H (=R2) THE'
1467 CALL
xerrwd(msg,44,685,0,0,0,0,2,tn,ho)
1468 msg =
'DASRT-- INITIAL YPRIME COULD NOT BE COMPUTED'
1469 CALL
xerrwd(msg,45,686,0,0,0,0,0,0.0d0,0.0d0)
1485 701 msg =
'DASRT-- SOME ELEMENT OF INFO VECTOR IS NOT ZERO OR ONE'
1486 CALL
xerrwd(msg,55,1,0,0,0,0,0,0.0d0,0.0d0)
1488 702 msg = .LE.
'DASRT-- NEQ (=I1) 0'
1489 CALL
xerrwd(msg,25,2,0,1,neq,0,0,0.0d0,0.0d0)
1491 703 msg =
'DASRT-- MAXORD (=I1) NOT IN RANGE'
1492 CALL
xerrwd(msg,34,3,0,1,mxord,0,0,0.0d0,0.0d0)
1494 704 msg=
'DASRT-- RWORK LENGTH NEEDED, LENRW (=I1), EXCEEDS LRW (=I2)'
1495 CALL
xerrwd(msg,60,4,0,2,lenrw,lrw,0,0.0d0,0.0d0)
1497 705 msg=
'DASRT-- IWORK LENGTH NEEDED, LENIW (=I1), EXCEEDS LIW (=I2)'
1498 CALL
xerrwd(msg,60,5,0,2,leniw,liw,0,0.0d0,0.0d0)
1500 706 msg = .LT.
'DASRT-- SOME ELEMENT OF RTOL IS 0'
1501 CALL
xerrwd(msg,39,6,0,0,0,0,0,0.0d0,0.0d0)
1503 707 msg = .LT.
'DASRT-- SOME ELEMENT OF ATOL IS 0'
1504 CALL
xerrwd(msg,39,7,0,0,0,0,0,0.0d0,0.0d0)
1506 708 msg =
'DASRT-- ALL ELEMENTS OF RTOL AND ATOL ARE ZERO'
1507 CALL
xerrwd(msg,47,8,0,0,0,0,0,0.0d0,0.0d0)
1509 709 msg=
'DASRT-- INFO(4) = 1 AND TSTOP (=R1) BEHIND TOUT (=R2)'
1510 CALL
xerrwd(msg,54,9,0,0,0,0,2,tstop,tout)
1512 710 msg = .LT.
'DASRT-- HMAX (=R1) 0.0'
1513 CALL
xerrwd(msg,28,10,0,0,0,0,1,hmax,0.0d0)
1515 711 msg =
'DASRT-- TOUT (=R1) BEHIND T (=R2)'
1516 CALL
xerrwd(msg,34,11,0,0,0,0,2,tout,t)
1518 712 msg =
'DASRT-- INFO(8)=1 AND H0=0.0'
1519 CALL
xerrwd(msg,29,12,0,0,0,0,0,0.0d0,0.0d0)
1521 713 msg = .LE.
'DASRT-- SOME ELEMENT OF WT IS 0.0'
1522 CALL
xerrwd(msg,39,13,0,0,0,0,0,0.0d0,0.0d0)
1524 714 msg=
'DASRT-- TOUT (=R1) TOO CLOSE TO T (=R2) TO START INTEGRATION'
1525 CALL
xerrwd(msg,60,14,0,0,0,0,2,tout,t)
1527 715 msg =
'DASRT-- INFO(4)=1 AND TSTOP (=R1) BEHIND T (=R2)'
1528 CALL
xerrwd(msg,49,15,0,0,0,0,2,tstop,t)
1530 716 msg = .LT.
'DASRT-- INFO(12)=1 AND MXSTP (=I1) 0'
1531 CALL
xerrwd(msg,42,16,0,1,iwork(lmxstp),0,0,0.0d0,0.0d0)
1533 717 msg = .LT..GT.
'DASRT-- ML (=I1) ILLEGAL. EITHER 0 OR NEQ'
1534 CALL
xerrwd(msg,52,17,0,1,iwork(lml),0,0,0.0d0,0.0d0)
1536 718 msg = .LT..GT.
'DASRT-- MU (=I1) ILLEGAL. EITHER 0 OR NEQ'
1537 CALL
xerrwd(msg,52,18,0,1,iwork(lmu),0,0,0.0d0,0.0d0)
1539 719 msg =
'DASRT-- TOUT (=R1) IS EQUAL TO T (=R2)'
1540 CALL
xerrwd(msg,39,19,0,0,0,0,2,tout,t)
1542 730 msg = .LT.
'DASRT-- NG (=I1) 0'
1543 CALL
xerrwd(msg,24,30,1,1,ng,0,0,0.0d0,0.0d0)
1545 732 msg =
'DASRT-- ONE OR MORE COMPONENTS OF G HAS A ROOT'
1546 CALL
xerrwd(msg,47,32,1,0,0,0,0,0.0d0,0.0d0)
1547 msg =
' TOO NEAR TO THE INITIAL POINT'
1548 CALL
xerrwd(msg,38,32,1,0,0,0,0,0.0d0,0.0d0)
1549 750
IF(info(1).EQ.-1) go to 760
1553 760 msg =
'DASRT-- REPEATED OCCURRENCES OF ILLEGAL INPUT'
1554 CALL
xerrwd(msg,46,801,0,0,0,0,0,0.0d0,0.0d0)
1555 770 msg =
'DASRT-- RUN TERMINATED. APPARENT INFINITE LOOP'
1556 CALL
xerrwd(msg,47,802,1,0,0,0,0,0.0d0,0.0d0)
subroutine ddasrt(RES, NEQ, T, Y, YPRIME, TOUT, INFO, RTOL, ATOL, IDID, RWORK, LRW, IWORK, LIW, RPAR, IPAR, JAC, G, NG, JROOT)
std::string dimension(void) const
subroutine xerrwd(MSG, NMES, NERR, LEVEL, NI, I1, I2, NR, R1, R2)
subroutine ddatrp(X, XOUT, YOUT, YPOUT, NEQ, KOLD, PHI, PSI)
double precision function d1mach(i)
subroutine ddawts(NEQ, IWT, RTOL, ATOL, Y, WT, RPAR, IPAR)
double precision function ddanrm(NEQ, V, WT, RPAR, IPAR)
subroutine ddastp(X, Y, YPRIME, NEQ, RES, JAC, H, WT, JSTART, IDID, RPAR, IPAR, PHI, DELTA, E, WM, IWM, ALPHA, BETA, GAMMA, PSI, SIGMA, CJ, CJOLD, HOLD, S, HMIN, UROUND, IPHASE, JCALC, K, KOLD, NS, NONNEG, NTEMP)
subroutine ddaini(X, Y, YPRIME, NEQ, RES, JAC, H, WT, IDID, RPAR, IPAR, PHI, DELTA, E, WM, IWM, HMIN, UROUND, NONNEG, NTEMP)
octave_value lgamma(void) const
subroutine drchek(JOB, G, NG, NEQ, TN, TOUT, Y, YP, PHI, PSI, KOLD, G0, G1, GX, JROOT, IRT, UROUND, INFO3, RWORK, IWORK, RPAR, IPAR)