1 SUBROUTINE dcuhre(NDIM,NUMFUN,A,B,MINPTS,MAXPTS,FUNSUB,EPSABS,
2 + EPSREL,KEY,NW,RESTAR,RESULT,ABSERR,NEVAL,IFAIL,
328 INTEGER NDIM,NUMFUN,MINPTS,MAXPTS,KEY,NW,RESTAR
330 DOUBLE PRECISION A(NDIM),B(NDIM),EPSABS,EPSREL
331 DOUBLE PRECISION RESULT(NUMFUN),ABSERR(NUMFUN),WORK(NW)
363 INTEGER MDIV,MAXWT,WTLENG,MAXDIM,LENW2,MAXSUB,MINSUB
364 INTEGER NUM,NSUB,LENW,KEYF
368 parameter(lenw2=2*mdiv*maxdim* (maxwt+1)+12*maxwt+2*maxdim)
369 INTEGER WRKSUB,I1,I2,I3,I4,I5,I6,I7,I8,K1,K2,K3,K4,K5,K6,K7,K8
370 DOUBLE PRECISION WORK2(LENW2)
377 CALL dchhre(maxdim,ndim,numfun,mdiv,a,b,minpts,maxpts,epsabs,
378 + epsrel,key,nw,restar,num,maxsub,minsub,keyf,
380 wrksub = (nw - 1 - 17*mdiv*numfun)/(2*ndim + 2*numfun + 2)
388 i2 = i1 + wrksub*numfun
389 i3 = i2 + wrksub*numfun
390 i4 = i3 + wrksub*ndim
391 i5 = i4 + wrksub*ndim
394 i8 = i7 + numfun*mdiv
396 k2 = k1 + 2*mdiv*wtleng*ndim
401 k7 = k6 + 2*mdiv*ndim
407 IF (restar.EQ.1)
THEN
413 lenw = 16*mdiv*numfun
414 CALL dadhre(ndim,numfun,mdiv,a,b,minsub,maxsub,funsub,epsabs,
415 + epsrel,keyf,restar,num,lenw,wtleng,
416 + result,abserr,neval,nsub,ifail,work(i1),work(i2),
417 + work(i3),work(i4),work(i5),work(i6),work(i7),work(i8),
418 + work2(k1),work2(k2),work2(k3),work2(k4),work2(k5),
419 + work2(k6),work2(k7),work2(k8))
426 SUBROUTINE dchhre(MAXDIM,NDIM,NUMFUN,MDIV,A,B,MINPTS,MAXPTS,
427 + EPSABS,EPSREL,KEY,NW,RESTAR,NUM,MAXSUB,MINSUB,
545 INTEGER NDIM,NUMFUN,MDIV,MINPTS,MAXPTS,KEY,NW,MINSUB,MAXSUB
546 INTEGER RESTAR,NUM,KEYF,IFAIL,MAXDIM,WTLENG
547 DOUBLE PRECISION A(NDIM),B(NDIM),EPSABS,EPSREL
559 IF (key.LT.0 .OR. key.GT.4)
THEN
566 IF (ndim.LT.2 .OR. ndim.GT.maxdim)
THEN
573 IF (key.EQ.1 .AND. ndim.NE.2)
THEN
580 IF (key.EQ.2 .AND. ndim.NE.3)
THEN
590 ELSE IF (ndim.EQ.3)
THEN
604 ELSE IF (keyf.EQ.2)
THEN
607 ELSE IF (keyf.EQ.3)
THEN
608 num = 1 + 4*2*ndim + 2*ndim* (ndim-1) + 4*ndim* (ndim-1) +
609 + 4*ndim* (ndim-1)* (ndim-2)/3 + 2**ndim
611 IF (ndim.EQ.2) wtleng = 8
612 ELSE IF (keyf.EQ.4)
THEN
613 num = 1 + 3*2*ndim + 2*ndim* (ndim-1) + 2**ndim
619 maxsub = (maxpts-num)/ (2*num) + 1
623 minsub = (minpts-num)/ (2*num) + 1
624 IF (mod(minpts-num,2*num).NE.0)
THEN
627 minsub = max(2,minsub)
631 IF (numfun.LT.1)
THEN
639 IF (a(j)-b(j).EQ.0)
THEN
647 IF (maxpts.LT.3*num)
THEN
654 IF (maxpts.LT.minpts)
THEN
661 IF (epsabs.LT.0 .AND. epsrel.LT.0)
THEN
668 limit = maxsub* (2*ndim+2*numfun+2) + 17*mdiv*numfun + 1
669 IF (nw.LT.limit)
THEN
676 IF (restar.NE.0 .AND. restar.NE.1)
THEN
685 SUBROUTINE dadhre(NDIM,NUMFUN,MDIV,A,B,MINSUB,MAXSUB,FUNSUB,
686 + EPSABS,EPSREL,KEY,RESTAR,NUM,LENW,WTLENG,
687 + RESULT,ABSERR,NEVAL,NSUB,IFAIL,VALUES,
688 + ERRORS,CENTRS,HWIDTS,GREATE,DIR,OLDRES,WORK,G,W,
689 + RULPTS,CENTER,HWIDTH,X,SCALES,NORMS)
862 INTEGER NDIM,NUMFUN,MDIV,MINSUB,MAXSUB,KEY,LENW,RESTAR
863 INTEGER NUM,NEVAL,NSUB,IFAIL,WTLENG
864 DOUBLE PRECISION A(NDIM),B(NDIM),EPSABS,EPSREL
865 DOUBLE PRECISION RESULT(NUMFUN),ABSERR(NUMFUN)
866 DOUBLE PRECISION VALUES(NUMFUN,MAXSUB),ERRORS(NUMFUN,MAXSUB)
867 DOUBLE PRECISION CENTRS(NDIM,MAXSUB)
868 DOUBLE PRECISION HWIDTS(NDIM,MAXSUB)
869 DOUBLE PRECISION GREATE(MAXSUB),DIR(MAXSUB)
870 DOUBLE PRECISION OLDRES(NUMFUN,MDIV)
871 DOUBLE PRECISION WORK(LENW),RULPTS(WTLENG)
872 DOUBLE PRECISION G(NDIM,WTLENG,2*MDIV),W(5,WTLENG)
873 DOUBLE PRECISION CENTER(NDIM),HWIDTH(NDIM),X(NDIM,2*MDIV)
874 DOUBLE PRECISION SCALES(3,WTLENG),NORMS(3,WTLENG)
888 INTEGER INTSGN,SBRGNS
890 INTEGER NDIV,POINTR,DIRECT,INDEX
891 DOUBLE PRECISION OLDCEN,EST1,EST2,ERRCOF(6)
899 IF (b(j).LT.a(j))
THEN
907 CALL dinhre(ndim,key,wtleng,w,g,errcof,rulpts,scales,norms)
911 IF (restar.EQ.1)
THEN
920 centrs(j,1) = (a(j)+b(j))/2
921 hwidts(j,1) = abs(b(j)-a(j))/2
934 CALL drlhre(ndim,centrs(1,1),hwidts(1,1),wtleng,g,w,errcof,numfun,
935 + funsub,scales,norms,x,work,values(1,1),errors(1,1),
942 result(j) = result(j) + values(j,1)
945 abserr(j) = abserr(j) + errors(j,1)
951 CALL dtrhre(2,ndim,numfun,index,values,errors,centrs,hwidts,
952 + greate,work(1),work(numfun+1),center,hwidth,dir)
959 110
IF (sbrgns+1.LE.maxsub)
THEN
967 ndiv = maxsub - sbrgns
968 ndiv = min(ndiv,mdiv,sbrgns)
977 pointr = sbrgns + ndiv + 1 - i
982 result(j) = result(j) - values(j,1)
983 abserr(j) = abserr(j) - errors(j,1)
989 centrs(j,pointr) = centrs(j,1)
990 hwidts(j,pointr) = hwidts(j,1)
994 hwidts(direct,pointr) = hwidts(direct,1)/2
995 oldcen = centrs(direct,1)
996 centrs(direct,pointr) = oldcen - hwidts(direct,pointr)
1001 oldres(j,ndiv-i+1) = values(j,1)
1006 CALL dtrhre(1,ndim,numfun,sbrgns,values,errors,centrs,
1007 + hwidts,greate,work(1),work(numfun+1),center,
1013 centrs(j,pointr-1) = centrs(j,pointr)
1014 hwidts(j,pointr-1) = hwidts(j,pointr)
1016 centrs(direct,pointr-1) = oldcen + hwidts(direct,pointr)
1017 hwidts(direct,pointr-1) = hwidts(direct,pointr)
1018 dir(pointr-1) = direct
1034 l1 = 1 + (i-1)*8*numfun
1035 CALL drlhre(ndim,centrs(1,index),hwidts(1,index),wtleng,
1036 + g(1,1,i),w,errcof,numfun,funsub,scales,norms,
1037 + x(1,i),work(l1),values(1,index),
1038 + errors(1,index),dir(index))
1040 neval = neval + 2*ndiv*num
1046 result(j) = result(j) + values(j,sbrgns+i)
1054 greate(sbrgns+2*i-1) = 0
1055 greate(sbrgns+2*i) = 0
1057 est1 = abs(oldres(j,i)- (values(j,
1058 + sbrgns+2*i-1)+values(j,sbrgns+2*i)))
1059 est2 = errors(j,sbrgns+2*i-1) + errors(j,sbrgns+2*i)
1061 errors(j,sbrgns+2*i-1) = errors(j,sbrgns+2*i-1)*
1062 + (1+errcof(5)*est1/est2)
1063 errors(j,sbrgns+2*i) = errors(j,sbrgns+2*i)*
1064 + (1+errcof(5)*est1/est2)
1066 errors(j,sbrgns+2*i-1) = errors(j,sbrgns+2*i-1) +
1068 errors(j,sbrgns+2*i) = errors(j,sbrgns+2*i) +
1070 IF (errors(j,sbrgns+2*i-1).GT.
1071 + greate(sbrgns+2*i-1))
THEN
1072 greate(sbrgns+2*i-1) = errors(j,sbrgns+2*i-1)
1074 IF (errors(j,sbrgns+2*i).GT.greate(sbrgns+2*i))
THEN
1075 greate(sbrgns+2*i) = errors(j,sbrgns+2*i)
1077 abserr(j) = abserr(j) + errors(j,sbrgns+2*i-1) +
1078 + errors(j,sbrgns+2*i)
1086 CALL dtrhre(2,ndim,numfun,index,values,errors,centrs,
1087 + hwidts,greate,work(1),work(numfun+1),center,
1090 sbrgns = sbrgns + 2*ndiv
1094 IF (sbrgns.LT.minsub)
THEN
1098 IF (abserr(j).GT.epsrel*abs(result(j)) .AND.
1099 + abserr(j).GT.epsabs)
THEN
1122 result(j) = result(j) + values(j,i)
1123 abserr(j) = abserr(j) + errors(j,i)
1130 result(j) = result(j)*intsgn
1138 SUBROUTINE dinhre(NDIM,KEY,WTLENG,W,G,ERRCOF,RULPTS,SCALES,NORMS)
1201 INTEGER NDIM,KEY,WTLENG
1202 DOUBLE PRECISION G(NDIM,WTLENG),W(5,WTLENG),ERRCOF(6)
1203 DOUBLE PRECISION RULPTS(WTLENG),SCALES(3,WTLENG)
1204 DOUBLE PRECISION NORMS(3,WTLENG)
1209 DOUBLE PRECISION WE(14)
1216 CALL d132re(wtleng,w,g,errcof,rulpts)
1217 ELSE IF (key.EQ.2)
THEN
1218 CALL d113re(wtleng,w,g,errcof,rulpts)
1219 ELSE IF (key.EQ.3)
THEN
1220 CALL d09hre(ndim,wtleng,w,g,errcof,rulpts)
1221 ELSE IF (key.EQ.4)
THEN
1222 CALL d07hre(ndim,wtleng,w,g,errcof,rulpts)
1229 IF (w(k+1,i).NE.0)
THEN
1230 scales(k,i) = - w(k+2,i)/w(k+1,i)
1235 we(j) = w(k+2,j) + scales(k,i)*w(k+1,j)
1239 norms(k,i) = norms(k,i) + rulpts(j)*abs(we(j))
1241 norms(k,i) = 2**ndim/norms(k,i)
1249 SUBROUTINE d132re(WTLENG,W,G,ERRCOF,RULPTS)
1296 DOUBLE PRECISION W(5,WTLENG),G(2,WTLENG),ERRCOF(6)
1297 DOUBLE PRECISION RULPTS(WTLENG)
1302 DOUBLE PRECISION DIM2G(16)
1303 DOUBLE PRECISION DIM2W(14,5)
1305 DATA (dim2g(i),i=1,16)/0.2517129343453109d+00,
1306 + 0.7013933644534266d+00,0.9590960631619962d+00,
1307 + 0.9956010478552127d+00,0.5000000000000000d+00,
1308 + 0.1594544658297559d+00,0.3808991135940188d+00,
1309 + 0.6582769255267192d+00,0.8761473165029315d+00,
1310 + 0.9982431840531980d+00,0.9790222658168462d+00,
1311 + 0.6492284325645389d+00,0.8727421201131239d+00,
1312 + 0.3582614645881228d+00,0.5666666666666666d+00,
1313 + 0.2077777777777778d+00/
1315 DATA (dim2w(i,1),i=1,14)/0.3379692360134460d-01,
1316 + 0.9508589607597761d-01,0.1176006468056962d+00,
1317 + 0.2657774586326950d-01,0.1701441770200640d-01,
1318 + 0.0000000000000000d+00,0.1626593098637410d-01,
1319 + 0.1344892658526199d+00,0.1328032165460149d+00,
1320 + 0.5637474769991870d-01,0.3908279081310500d-02,
1321 + 0.3012798777432150d-01,0.1030873234689166d+00,
1322 + 0.6250000000000000d-01/
1324 DATA (dim2w(i,2),i=1,14)/0.3213775489050763d+00,
1325 + - .1767341636743844d+00,0.7347600537466072d-01,
1326 + - .3638022004364754d-01,0.2125297922098712d-01,
1327 + 0.1460984204026913d+00,0.1747613286152099d-01,
1328 + 0.1444954045641582d+00,0.1307687976001325d-03,
1329 + 0.5380992313941161d-03,0.1042259576889814d-03,
1330 + - .1401152865045733d-02,0.8041788181514763d-02,
1331 + - .1420416552759383d+00/
1333 DATA (dim2w(i,3),i=1,14)/0.3372900883288987d+00,
1334 + - .1644903060344491d+00,0.7707849911634622d-01,
1335 + - .3804478358506310d-01,0.2223559940380806d-01,
1336 + 0.1480693879765931d+00,0.4467143702185814d-05,
1337 + 0.1508944767074130d+00,0.3647200107516215d-04,
1338 + 0.5777198999013880d-03,0.1041757313688177d-03,
1339 + - .1452822267047819d-02,0.8338339968783705d-02,
1340 + - .1472796329231960d+00/
1342 DATA (dim2w(i,4),i=1,14)/ - .8264123822525677d+00,
1343 + 0.3065838614094360d+00,0.2389292538329435d-02,
1344 + - .1343024157997222d+00,0.8833366840533900d-01,
1345 + 0.0000000000000000d+00,0.9786283074168292d-03,
1346 + - .1319227889147519d+00,0.7990012200150630d-02,
1347 + 0.3391747079760626d-02,0.2294915718283264d-02,
1348 + - .1358584986119197d-01,0.4025866859057809d-01,
1349 + 0.3760268580063992d-02/
1351 DATA (dim2w(i,5),i=1,14)/0.6539094339575232d+00,
1352 + - .2041614154424632d+00, - .1746981515794990d+00,
1353 + 0.3937939671417803d-01,0.6974520545933992d-02,
1354 + 0.0000000000000000d+00,0.6667702171778258d-02,
1355 + 0.5512960621544304d-01,0.5443846381278607d-01,
1356 + 0.2310903863953934d-01,0.1506937747477189d-01,
1357 + - .6057021648901890d-01,0.4225737654686337d-01,
1358 + 0.2561989142123099d-01/
1420 SUBROUTINE d113re(WTLENG,W,G,ERRCOF,RULPTS)
1472 DOUBLE PRECISION W(5,WTLENG),G(3,WTLENG),ERRCOF(6)
1473 DOUBLE PRECISION RULPTS(WTLENG)
1478 DOUBLE PRECISION DIM3G(14)
1479 DOUBLE PRECISION DIM3W(13,5)
1481 DATA (dim3g(i),i=1,14)/0.1900000000000000d+00,
1482 + 0.5000000000000000d+00,0.7500000000000000d+00,
1483 + 0.8000000000000000d+00,0.9949999999999999d+00,
1484 + 0.9987344998351400d+00,0.7793703685672423d+00,
1485 + 0.9999698993088767d+00,0.7902637224771788d+00,
1486 + 0.4403396687650737d+00,0.4378478459006862d+00,
1487 + 0.9549373822794593d+00,0.9661093133630748d+00,
1488 + 0.4577105877763134d+00/
1490 DATA (dim3w(i,1),i=1,13)/0.7923078151105734d-02,
1491 + 0.6797177392788080d-01,0.1086986538805825d-02,
1492 + 0.1838633662212829d+00,0.3362119777829031d-01,
1493 + 0.1013751123334062d-01,0.1687648683985235d-02,
1494 + 0.1346468564512807d+00,0.1750145884600386d-02,
1495 + 0.7752336383837454d-01,0.2461864902770251d+00,
1496 + 0.6797944868483039d-01,0.1419962823300713d-01/
1498 DATA (dim3w(i,2),i=1,13)/0.1715006248224684d+01,
1499 + - .3755893815889209d+00,0.1488632145140549d+00,
1500 + - .2497046640620823d+00,0.1792501419135204d+00,
1501 + 0.3446126758973890d-02, - .5140483185555825d-02,
1502 + 0.6536017839876425d-02, - .6513454939229700d-03,
1503 + - .6304672433547204d-02,0.1266959399788263d-01,
1504 + - .5454241018647931d-02,0.4826995274768427d-02/
1506 DATA (dim3w(i,3),i=1,13)/0.1936014978949526d+01,
1507 + - .3673449403754268d+00,0.2929778657898176d-01,
1508 + - .1151883520260315d+00,0.5086658220872218d-01,
1509 + 0.4453911087786469d-01, - .2287828257125900d-01,
1510 + 0.2908926216345833d-01, - .2898884350669207d-02,
1511 + - .2805963413307495d-01,0.5638741361145884d-01,
1512 + - .2427469611942451d-01,0.2148307034182882d-01/
1514 DATA (dim3w(i,4),i=1,13)/0.5170828195605760d+00,
1515 + 0.1445269144914044d-01, - .3601489663995932d+00,
1516 + 0.3628307003418485d+00,0.7148802650872729d-02,
1517 + - .9222852896022966d-01,0.1719339732471725d-01,
1518 + - .1021416537460350d+00, - .7504397861080493d-02,
1519 + 0.1648362537726711d-01,0.5234610158469334d-01,
1520 + 0.1445432331613066d-01,0.3019236275367777d-02/
1522 DATA (dim3w(i,5),i=1,13)/0.2054404503818520d+01,
1523 + 0.1377759988490120d-01, - .5768062917904410d+00,
1524 + 0.3726835047700328d-01,0.6814878939777219d-02,
1525 + 0.5723169733851849d-01, - .4493018743811285d-01,
1526 + 0.2729236573866348d-01,0.3547473950556990d-03,
1527 + 0.1571366799739551d-01,0.4990099219278567d-01,
1528 + 0.1377915552666770d-01,0.2878206423099872d-02/
1599 SUBROUTINE d09hre(NDIM,WTLENG,W,G,ERRCOF,RULPTS)
1642 DOUBLE PRECISION W(5,WTLENG),G(NDIM,WTLENG),ERRCOF(6)
1643 DOUBLE PRECISION RULPTS(WTLENG)
1647 DOUBLE PRECISION RATIO,LAM0,LAM1,LAM2,LAM3,LAMP,TWONDM
1665 rulpts(wtleng) = twondm
1666 IF (ndim.GT.2) rulpts(8) = (4*ndim* (ndim-1)* (ndim-2))/3
1667 rulpts(7) = 4*ndim* (ndim-1)
1668 rulpts(6) = 2*ndim* (ndim-1)
1674 lam1 = 4/ (15-5/lam0)
1675 ratio = (1-lam1/lam0)/27
1676 lam2 = (5-7*lam1-35*ratio)/ (7-35*lam1/3-35*ratio/lam0)
1677 ratio = ratio* (1-lam2/lam0)/3
1678 lam3 = (7-9* (lam2+lam1)+63*lam2*lam1/5-63*ratio)/
1679 + (9-63* (lam2+lam1)/5+21*lam2*lam1-63*ratio/lam0)
1684 w(1,wtleng) = 1/ (3*lam0)**4/twondm
1685 IF (ndim.GT.2) w(1,8) = (1-1/ (3*lam0))/ (6*lam1)**3
1686 w(1,7) = (1-7* (lam0+lam1)/5+7*lam0*lam1/3)/
1687 + (84*lam1*lam2* (lam2-lam0)* (lam2-lam1))
1688 w(1,6) = (1-7* (lam0+lam2)/5+7*lam0*lam2/3)/
1689 + (84*lam1*lam1* (lam1-lam0)* (lam1-lam2)) -
1690 + w(1,7)*lam2/lam1 - 2* (ndim-2)*w(1,8)
1691 w(1,4) = (1-9* ((lam0+lam1+lam2)/7- (lam0*lam1+lam0*lam2+
1692 + lam1*lam2)/5)-3*lam0*lam1*lam2)/
1693 + (18*lam3* (lam3-lam0)* (lam3-lam1)* (lam3-lam2))
1694 w(1,3) = (1-9* ((lam0+lam1+lam3)/7- (lam0*lam1+lam0*lam3+
1695 + lam1*lam3)/5)-3*lam0*lam1*lam3)/
1696 + (18*lam2* (lam2-lam0)* (lam2-lam1)* (lam2-lam3)) -
1697 + 2* (ndim-1)*w(1,7)
1698 w(1,2) = (1-9* ((lam0+lam2+lam3)/7- (lam0*lam2+lam0*lam3+
1699 + lam2*lam3)/5)-3*lam0*lam2*lam3)/
1700 + (18*lam1* (lam1-lam0)* (lam1-lam2)* (lam1-lam3)) -
1701 + 2* (ndim-1)* (w(1,7)+w(1,6)+ (ndim-2)*w(1,8))
1705 w(2,wtleng) = 1/ (108*lam0**4)/twondm
1706 IF (ndim.GT.2) w(2,8) = (1-27*twondm*w(2,9)*lam0**3)/ (6*lam1)**3
1707 w(2,7) = (1-5*lam1/3-15*twondm*w(2,wtleng)*lam0**2* (lam0-lam1))/
1708 + (60*lam1*lam2* (lam2-lam1))
1709 w(2,6) = (1-9* (8*lam1*lam2*w(2,7)+twondm*w(2,wtleng)*lam0**2))/
1710 + (36*lam1*lam1) - 2*w(2,8)* (ndim-2)
1711 w(2,4) = (1-7* ((lam1+lam2)/5-lam1*lam2/3+twondm*w(2,
1712 + wtleng)*lam0* (lam0-lam1)* (lam0-lam2)))/
1713 + (14*lam3* (lam3-lam1)* (lam3-lam2))
1714 w(2,3) = (1-7* ((lam1+lam3)/5-lam1*lam3/3+twondm*w(2,
1715 + wtleng)*lam0* (lam0-lam1)* (lam0-lam3)))/
1716 + (14*lam2* (lam2-lam1)* (lam2-lam3)) - 2* (ndim-1)*w(2,7)
1717 w(2,2) = (1-7* ((lam2+lam3)/5-lam2*lam3/3+twondm*w(2,
1718 + wtleng)*lam0* (lam0-lam2)* (lam0-lam3)))/
1719 + (14*lam1* (lam1-lam2)* (lam1-lam3)) -
1720 + 2* (ndim-1)* (w(2,7)+w(2,6)+ (ndim-2)*w(2,8))
1721 w(3,wtleng) = 5/ (324*lam0**4)/twondm
1722 IF (ndim.GT.2) w(3,8) = (1-27*twondm*w(3,9)*lam0**3)/ (6*lam1)**3
1723 w(3,7) = (1-5*lam1/3-15*twondm*w(3,wtleng)*lam0**2* (lam0-lam1))/
1724 + (60*lam1*lam2* (lam2-lam1))
1725 w(3,6) = (1-9* (8*lam1*lam2*w(3,7)+twondm*w(3,wtleng)*lam0**2))/
1726 + (36*lam1*lam1) - 2*w(3,8)* (ndim-2)
1727 w(3,5) = (1-7* ((lam1+lam2)/5-lam1*lam2/3+twondm*w(3,
1728 + wtleng)*lam0* (lam0-lam1)* (lam0-lam2)))/
1729 + (14*lamp* (lamp-lam1)* (lamp-lam2))
1730 w(3,3) = (1-7* ((lam1+lamp)/5-lam1*lamp/3+twondm*w(3,
1731 + wtleng)*lam0* (lam0-lam1)* (lam0-lamp)))/
1732 + (14*lam2* (lam2-lam1)* (lam2-lamp)) - 2* (ndim-1)*w(3,7)
1733 w(3,2) = (1-7* ((lam2+lamp)/5-lam2*lamp/3+twondm*w(3,
1734 + wtleng)*lam0* (lam0-lam2)* (lam0-lamp)))/
1735 + (14*lam1* (lam1-lam2)* (lam1-lamp)) -
1736 + 2* (ndim-1)* (w(3,7)+w(3,6)+ (ndim-2)*w(3,8))
1737 w(4,wtleng) = 2/ (81*lam0**4)/twondm
1738 IF (ndim.GT.2) w(4,8) = (2-27*twondm*w(4,9)*lam0**3)/ (6*lam1)**3
1739 w(4,7) = (2-15*lam1/9-15*twondm*w(4,wtleng)*lam0* (lam0-lam1))/
1740 + (60*lam1*lam2* (lam2-lam1))
1741 w(4,6) = (1-9* (8*lam1*lam2*w(4,7)+twondm*w(4,wtleng)*lam0**2))/
1742 + (36*lam1*lam1) - 2*w(4,8)* (ndim-2)
1743 w(4,4) = (2-7* ((lam1+lam2)/5-lam1*lam2/3+twondm*w(4,
1744 + wtleng)*lam0* (lam0-lam1)* (lam0-lam2)))/
1745 + (14*lam3* (lam3-lam1)* (lam3-lam2))
1746 w(4,3) = (2-7* ((lam1+lam3)/5-lam1*lam3/3+twondm*w(4,
1747 + wtleng)*lam0* (lam0-lam1)* (lam0-lam3)))/
1748 + (14*lam2* (lam2-lam1)* (lam2-lam3)) - 2* (ndim-1)*w(4,7)
1749 w(4,2) = (2-7* ((lam2+lam3)/5-lam2*lam3/3+twondm*w(4,
1750 + wtleng)*lam0* (lam0-lam2)* (lam0-lam3)))/
1751 + (14*lam1* (lam1-lam2)* (lam1-lam3)) -
1752 + 2* (ndim-1)* (w(4,7)+w(4,6)+ (ndim-2)*w(4,8))
1753 w(5,2) = 1/ (6*lam1)
1786 w(j,i) = w(j,i) - w(1,i)
1787 w(j,1) = w(j,1) - rulpts(i)*w(j,i)
1791 w(1,i) = twondm*w(1,i)
1792 w(1,1) = w(1,1) - rulpts(i)*w(1,i)
1807 SUBROUTINE d07hre(NDIM,WTLENG,W,G,ERRCOF,RULPTS)
1851 DOUBLE PRECISION W(5,WTLENG),G(NDIM,WTLENG),ERRCOF(6)
1852 DOUBLE PRECISION RULPTS(WTLENG)
1856 DOUBLE PRECISION RATIO,LAM0,LAM1,LAM2,LAMP,TWONDM
1874 rulpts(wtleng) = twondm
1875 rulpts(wtleng-1) = 2*ndim* (ndim-1)
1882 lam1 = 4/ (15-5/lam0)
1883 ratio = (1-lam1/lam0)/27
1884 lam2 = (5-7*lam1-35*ratio)/ (7-35*lam1/3-35*ratio/lam0)
1888 w(1,6) = 1/ (3*lam0)**3/twondm
1889 w(1,5) = (1-5*lam0/3)/ (60* (lam1-lam0)*lam1**2)
1890 w(1,3) = (1-5*lam2/3-5*twondm*w(1,6)*lam0* (lam0-lam2))/
1891 + (10*lam1* (lam1-lam2)) - 2* (ndim-1)*w(1,5)
1892 w(1,2) = (1-5*lam1/3-5*twondm*w(1,6)*lam0* (lam0-lam1))/
1893 + (10*lam2* (lam2-lam1))
1897 w(2,6) = 1/ (36*lam0**3)/twondm
1898 w(2,5) = (1-9*twondm*w(2,6)*lam0**2)/ (36*lam1**2)
1899 w(2,3) = (1-5*lam2/3-5*twondm*w(2,6)*lam0* (lam0-lam2))/
1900 + (10*lam1* (lam1-lam2)) - 2* (ndim-1)*w(2,5)
1901 w(2,2) = (1-5*lam1/3-5*twondm*w(2,6)*lam0* (lam0-lam1))/
1902 + (10*lam2* (lam2-lam1))
1903 w(3,6) = 5/ (108*lam0**3)/twondm
1904 w(3,5) = (1-9*twondm*w(3,6)*lam0**2)/ (36*lam1**2)
1905 w(3,3) = (1-5*lamp/3-5*twondm*w(3,6)*lam0* (lam0-lamp))/
1906 + (10*lam1* (lam1-lamp)) - 2* (ndim-1)*w(3,5)
1907 w(3,4) = (1-5*lam1/3-5*twondm*w(3,6)*lam0* (lam0-lam1))/
1908 + (10*lamp* (lamp-lam1))
1909 w(4,6) = 1/ (54*lam0**3)/twondm
1910 w(4,5) = (1-18*twondm*w(4,6)*lam0**2)/ (72*lam1**2)
1911 w(4,3) = (1-10*lam2/3-10*twondm*w(4,6)*lam0* (lam0-lam2))/
1912 + (20*lam1* (lam1-lam2)) - 2* (ndim-1)*w(4,5)
1913 w(4,2) = (1-10*lam1/3-10*twondm*w(4,6)*lam0* (lam0-lam1))/
1914 + (20*lam2* (lam2-lam1))
1925 g(1,wtleng-1) = lam1
1926 g(2,wtleng-1) = lam1
1927 g(1,wtleng-4) = lam2
1928 g(1,wtleng-3) = lam1
1929 g(1,wtleng-2) = lamp
1938 w(j,i) = w(j,i) - w(1,i)
1939 w(j,1) = w(j,1) - rulpts(i)*w(j,i)
1943 w(1,i) = twondm*w(1,i)
1944 w(1,1) = w(1,1) - rulpts(i)*w(1,i)
1959 SUBROUTINE drlhre(NDIM,CENTER,HWIDTH,WTLENG,G,W,ERRCOF,NUMFUN,
1960 + FUNSUB,SCALES,NORMS,X,NULL,BASVAL,RGNERR,DIRECT)
2055 INTEGER WTLENG,NUMFUN,NDIM
2056 DOUBLE PRECISION CENTER(NDIM),X(NDIM),HWIDTH(NDIM),BASVAL(NUMFUN),
2057 + RGNERR(NUMFUN),NULL(NUMFUN,8),W(5,WTLENG),
2058 + G(NDIM,WTLENG),ERRCOF(6),DIRECT,SCALES(3,WTLENG),
2063 DOUBLE PRECISION RGNVOL,DIFSUM,DIFMAX,FRTHDF
2064 INTEGER I,J,K,DIVAXN
2065 DOUBLE PRECISION SEARCH,RATIO
2078 rgnvol = rgnvol*hwidth(i)
2080 IF (hwidth(i).GT.hwidth(divaxn)) divaxn = i
2082 CALL funsub(ndim,x,numfun,rgnerr)
2084 basval(j) = w(1,1)*rgnerr(j)
2086 null(j,k) = w(k+1,1)*rgnerr(j)
2090 ratio = (g(1,3)/g(1,2))**2
2092 x(i) = center(i) - hwidth(i)*g(1,2)
2093 CALL funsub(ndim,x,numfun,null(1,5))
2094 x(i) = center(i) + hwidth(i)*g(1,2)
2095 CALL funsub(ndim,x,numfun,null(1,6))
2096 x(i) = center(i) - hwidth(i)*g(1,3)
2097 CALL funsub(ndim,x,numfun,null(1,7))
2098 x(i) = center(i) + hwidth(i)*g(1,3)
2099 CALL funsub(ndim,x,numfun,null(1,8))
2103 frthdf = 2* (1-ratio)*rgnerr(j) - (null(j,7)+null(j,8)) +
2104 + ratio* (null(j,5)+null(j,6))
2108 IF (rgnerr(j)+frthdf/4.NE.rgnerr(j)) difsum = difsum +
2111 null(j,k) = null(j,k) + w(k+1,2)*
2112 + (null(j,5)+null(j,6)) +
2113 + w(k+1,3)* (null(j,7)+null(j,8))
2115 basval(j) = basval(j) + w(1,2)* (null(j,5)+null(j,6)) +
2116 + w(1,3)* (null(j,7)+null(j,8))
2118 IF (difsum.GT.difmax)
THEN
2128 CALL dfshre(ndim,center,hwidth,x,g(1,i),numfun,funsub,rgnerr,
2131 basval(j) = basval(j) + w(1,i)*rgnerr(j)
2133 null(j,k) = null(j,k) + w(k+1,i)*rgnerr(j)
2150 search = max(search,abs(null(j,i+1)+scales(i,
2151 + k)*null(j,i))*norms(i,k))
2155 IF (errcof(1)*null(j,1).LE.null(j,2) .AND.
2156 + errcof(2)*null(j,2).LE.null(j,3))
THEN
2157 rgnerr(j) = errcof(3)*null(j,1)
2159 rgnerr(j) = errcof(4)*max(null(j,1),null(j,2),null(j,3))
2161 rgnerr(j) = rgnvol*rgnerr(j)
2162 basval(j) = rgnvol*basval(j)
2168 SUBROUTINE dfshre(NDIM,CENTER,HWIDTH,X,G,NUMFUN,FUNSUB,FULSMS,
2226 DOUBLE PRECISION CENTER(NDIM),HWIDTH(NDIM),X(NDIM),G(NDIM),
2227 + FULSMS(NUMFUN),FUNVLS(NUMFUN)
2231 INTEGER IXCHNG,LXCHNG,I,J,L
2232 DOUBLE PRECISION GL,GI
2243 x(i) = center(i) + g(i)*hwidth(i)
2245 40
CALL funsub(ndim,x,numfun,funvls)
2247 fulsms(j) = fulsms(j) + funvls(j)
2251 x(i) = center(i) + g(i)*hwidth(i)
2252 IF (g(i).LT.0)
GO TO 40
2259 IF (g(i-1).GT.g(i))
THEN
2262 DO 70 l = 1, (i-1)/2
2266 IF (gl.LE.gi) ixchng = ixchng - 1
2267 IF (g(l).GT.gi) lxchng = l
2269 IF (g(ixchng).LE.gi) ixchng = lxchng
2287 SUBROUTINE dtrhre(DVFLAG,NDIM,NUMFUN,SBRGNS,VALUES,ERRORS,CENTRS,
2288 + HWIDTS,GREATE,ERROR,VALUE,CENTER,HWIDTH,DIR)
2339 INTEGER DVFLAG,NDIM,NUMFUN,SBRGNS
2340 DOUBLE PRECISION VALUES(NUMFUN,*),ERRORS(NUMFUN,*)
2341 DOUBLE PRECISION CENTRS(NDIM,*)
2342 DOUBLE PRECISION HWIDTS(NDIM,*)
2343 DOUBLE PRECISION GREATE(*)
2344 DOUBLE PRECISION ERROR(NUMFUN),
VALUE(NUMFUN)
2345 DOUBLE PRECISION CENTER(NDIM),HWIDTH(NDIM)
2346 DOUBLE PRECISION DIR(*)
2357 INTEGER J,SUBRGN,SUBTMP
2358 DOUBLE PRECISION GREAT,DIRECT
2364 GREAT = greate(sbrgns)
2365 direct = dir(sbrgns)
2367 error(j) = errors(j,sbrgns)
2368 value(j) = values(j,sbrgns)
2371 center(j) = centrs(j,sbrgns)
2372 hwidth(j) = hwidts(j,sbrgns)
2378 IF (dvflag.EQ.1)
THEN
2381 20 subtmp = 2*subrgn
2382 IF (subtmp.LE.sbrgns)
THEN
2383 IF (subtmp.NE.sbrgns)
THEN
2387 IF (greate(subtmp).LT.greate(subtmp+1))
THEN
2395 IF (great.LT.greate(subtmp))
THEN
2399 greate(subrgn) = greate(subtmp)
2401 errors(j,subrgn) = errors(j,subtmp)
2402 values(j,subrgn) = values(j,subtmp)
2404 dir(subrgn) = dir(subtmp)
2406 centrs(j,subrgn) = centrs(j,subtmp)
2407 hwidts(j,subrgn) = hwidts(j,subtmp)
2413 ELSE IF (dvflag.EQ.2)
THEN
2418 40 subtmp = subrgn/2
2419 IF (subtmp.GE.1)
THEN
2424 IF (great.GT.greate(subtmp))
THEN
2428 greate(subrgn) = greate(subtmp)
2430 errors(j,subrgn) = errors(j,subtmp)
2431 values(j,subrgn) = values(j,subtmp)
2433 dir(subrgn) = dir(subtmp)
2435 centrs(j,subrgn) = centrs(j,subtmp)
2436 hwidts(j,subrgn) = hwidts(j,subtmp)
2446 IF (sbrgns.GT.0)
THEN
2447 greate(subrgn) = great
2449 errors(j,subrgn) = error(j)
2450 values(j,subrgn) = value(j)
2452 dir(subrgn) = direct
2454 centrs(j,subrgn) = center(j)
2455 hwidts(j,subrgn) = hwidth(j)