1 subroutine dqag(f,a,b,epsabs,epsrel,key,result,abserr,neval,ier,
2 * limit,lenw,last,iwork,work)
150 double precision a,abserr,b,epsabs,epsrel,f,result,work
151 integer ier,iwork,key,last,lenw,limit,lvl,l1,l2,l3,neval
153 dimension iwork(limit),work(lenw)
165 if(limit.lt.1.or.lenw.lt.limit*4)
go to 10
173 call dqage(f,a,b,epsabs,epsrel,key,limit,result,abserr,neval,
174 * ier,work(1),work(l1),work(l2),work(l3),iwork,last)
179 10
if(ier.eq.6) lvl = 1
184 subroutine dqage(f,a,b,epsabs,epsrel,key,limit,result,abserr,
185 * neval,ier,alist,blist,rlist,elist,iord,last)
331 double precision a,abserr,alist,area,area1,area12,area2,a1,a2,b,
332 * blist,b1,b2,dabs,defabs,defab1,defab2,dmax1,d1mach,elist,epmach,
333 * epsabs,epsrel,errbnd,errmax,error1,error2,erro12,errsum,f,
334 * resabs,result,rlist,uflow
335 integer ier,iord,iroff1,iroff2,k,key,keyf,last,limit,maxerr,neval,
338 dimension alist(limit),blist(limit),elist(limit),iord(limit),
388 if(epsabs.le.0.0d+00.and.
389 * epsrel.lt.dmax1(0.5d+02*epmach,0.5d-28)) ier = 6
390 if(ier.eq.6)
go to 999
396 if(key.le.0) keyf = 1
397 if(key.ge.7) keyf = 6
399 if(keyf.eq.1)
call dqk15(f,a,b,result,abserr,defabs,resabs)
400 if(keyf.eq.2)
call dqk21(f,a,b,result,abserr,defabs,resabs)
401 if(keyf.eq.3)
call dqk31(f,a,b,result,abserr,defabs,resabs)
402 if(keyf.eq.4)
call dqk41(f,a,b,result,abserr,defabs,resabs)
403 if(keyf.eq.5)
call dqk51(f,a,b,result,abserr,defabs,resabs)
404 if(keyf.eq.6)
call dqk61(f,a,b,result,abserr,defabs,resabs)
412 errbnd = dmax1(epsabs,epsrel*dabs(result))
413 if(abserr.le.0.5d+02*epmach*defabs.and.abserr.gt.errbnd) ier = 2
414 if(limit.eq.1) ier = 1
415 if(ier.ne.0.or.(abserr.le.errbnd.and.abserr.ne.resabs)
416 * .or.abserr.eq.0.0d+00)
go to 60
438 b1 = 0.5d+00*(alist(maxerr)+blist(maxerr))
441 if(keyf.eq.1)
call dqk15(f,a1,b1,area1,error1,resabs,defab1)
442 if(keyf.eq.2)
call dqk21(f,a1,b1,area1,error1,resabs,defab1)
443 if(keyf.eq.3)
call dqk31(f,a1,b1,area1,error1,resabs,defab1)
444 if(keyf.eq.4)
call dqk41(f,a1,b1,area1,error1,resabs,defab1)
445 if(keyf.eq.5)
call dqk51(f,a1,b1,area1,error1,resabs,defab1)
446 if(keyf.eq.6)
call dqk61(f,a1,b1,area1,error1,resabs,defab1)
447 if(keyf.eq.1)
call dqk15(f,a2,b2,area2,error2,resabs,defab2)
448 if(keyf.eq.2)
call dqk21(f,a2,b2,area2,error2,resabs,defab2)
449 if(keyf.eq.3)
call dqk31(f,a2,b2,area2,error2,resabs,defab2)
450 if(keyf.eq.4)
call dqk41(f,a2,b2,area2,error2,resabs,defab2)
451 if(keyf.eq.5)
call dqk51(f,a2,b2,area2,error2,resabs,defab2)
452 if(keyf.eq.6)
call dqk61(f,a2,b2,area2,error2,resabs,defab2)
459 erro12 = error1+error2
460 errsum = errsum+erro12-errmax
461 area = area+area12-rlist(maxerr)
462 if(defab1.eq.error1.or.defab2.eq.error2)
go to 5
463 if(dabs(rlist(maxerr)-area12).le.0.1d-04*dabs(area12)
464 * .and.erro12.ge.0.99d+00*errmax) iroff1 = iroff1+1
465 if(last.gt.10.and.erro12.gt.errmax) iroff2 = iroff2+1
466 5 rlist(maxerr) = area1
468 errbnd = dmax1(epsabs,epsrel*dabs(area))
469 if(errsum.le.errbnd)
go to 8
473 if(iroff1.ge.6.or.iroff2.ge.20) ier = 2
478 if(last.eq.limit) ier = 1
483 if(dmax1(dabs(a1),dabs(b2)).le.(0.1d+01+0.1d+03*
484 * epmach)*(dabs(a2)+0.1d+04*uflow)) ier = 3
488 8
if(error2.gt.error1)
go to 10
492 elist(maxerr) = error1
495 10 alist(maxerr) = a2
498 rlist(maxerr) = area2
500 elist(maxerr) = error2
507 20
call dqpsrt(limit,last,maxerr,errmax,elist,iord,nrmax)
509 if(ier.ne.0.or.errsum.le.errbnd)
go to 40
517 result = result+rlist(k)
520 60
if(keyf.ne.1) neval = (10*keyf+1)*(2*neval+1)
521 if(keyf.eq.1) neval = 30*neval+15
524 subroutine dqk15(f,a,b,result,abserr,resabs,resasc)
576 double precision a,absc,abserr,b,centr,dabs,dhlgth,dmax1,dmin1,
577 * d1mach,epmach,f,fc,fsum,fval1,fval2,fv1,fv2,hlgth,resabs,resasc,
578 * resg,resk,reskh,result,uflow,wg,wgk,xgk
582 dimension fv1(7),fv2(7),wg(4),wgk(8),xgk(8)
603 data wg( 1) / 0.1294849661 6886969327 0611432679 082 d0 /
604 data wg( 2) / 0.2797053914 8927666790 1467771423 780 d0 /
605 data wg( 3) / 0.3818300505 0511894495 0369775488 975 d0 /
606 data wg( 4) / 0.4179591836 7346938775 5102040816 327 d0 /
608 data xgk( 1) / 0.9914553711 2081263920 6854697526 329 d0 /
609 data xgk( 2) / 0.9491079123 4275852452 6189684047 851 d0 /
610 data xgk( 3) / 0.8648644233 5976907278 9712788640 926 d0 /
611 data xgk( 4) / 0.7415311855 9939443986 3864773280 788 d0 /
612 data xgk( 5) / 0.5860872354 6769113029 4144838258 730 d0 /
613 data xgk( 6) / 0.4058451513 7739716690 6606412076 961 d0 /
614 data xgk( 7) / 0.2077849550 0789846760 0689403773 245 d0 /
615 data xgk( 8) / 0.0000000000 0000000000 0000000000 000 d0 /
617 data wgk( 1) / 0.0229353220 1052922496 3732008058 970 d0 /
618 data wgk( 2) / 0.0630920926 2997855329 0700663189 204 d0 /
619 data wgk( 3) / 0.1047900103 2225018383 9876322541 518 d0 /
620 data wgk( 4) / 0.1406532597 1552591874 5189590510 238 d0 /
621 data wgk( 5) / 0.1690047266 3926790282 6583426598 550 d0 /
622 data wgk( 6) / 0.1903505780 6478540991 3256402421 014 d0 /
623 data wgk( 7) / 0.2044329400 7529889241 4161999234 649 d0 /
624 data wgk( 8) / 0.2094821410 8472782801 2999174891 714 d0 /
649 centr = 0.5d+00*(a+b)
650 hlgth = 0.5d+00*(b-a)
662 absc = hlgth*xgk(jtw)
663 fval1 = f(centr-absc)
664 fval2 = f(centr+absc)
668 resg = resg+wg(j)*fsum
669 resk = resk+wgk(jtw)*fsum
670 resabs = resabs+wgk(jtw)*(dabs(fval1)+dabs(fval2))
674 absc = hlgth*xgk(jtwm1)
675 fval1 = f(centr-absc)
676 fval2 = f(centr+absc)
680 resk = resk+wgk(jtwm1)*fsum
681 resabs = resabs+wgk(jtwm1)*(dabs(fval1)+dabs(fval2))
684 resasc = wgk(8)*dabs(fc-reskh)
686 resasc = resasc+wgk(j)*(dabs(fv1(j)-reskh)+dabs(fv2(j)-reskh))
689 resabs = resabs*dhlgth
690 resasc = resasc*dhlgth
691 abserr = dabs((resk-resg)*hlgth)
692 if(resasc.ne.0.0d+00.and.abserr.ne.0.0d+00)
693 * abserr = resasc*dmin1(0.1d+01,(0.2d+03*abserr/resasc)**1.5d+00)
694 if(resabs.gt.uflow/(0.5d+02*epmach)) abserr = dmax1
695 * ((epmach*0.5d+02)*resabs,abserr)
698 subroutine dqk21(f,a,b,result,abserr,resabs,resasc)
750 double precision a,absc,abserr,b,centr,dabs,dhlgth,dmax1,dmin1,
751 * d1mach,epmach,f,fc,fsum,fval1,fval2,fv1,fv2,hlgth,resabs,resasc,
752 * resg,resk,reskh,result,uflow,wg,wgk,xgk
756 dimension fv1(10),fv2(10),wg(5),wgk(11),xgk(11)
777 data wg( 1) / 0.0666713443 0868813759 3568809893 332 d0 /
778 data wg( 2) / 0.1494513491 5058059314 5776339657 697 d0 /
779 data wg( 3) / 0.2190863625 1598204399 5534934228 163 d0 /
780 data wg( 4) / 0.2692667193 0999635509 1226921569 469 d0 /
781 data wg( 5) / 0.2955242247 1475287017 3892994651 338 d0 /
783 data xgk( 1) / 0.9956571630 2580808073 5527280689 003 d0 /
784 data xgk( 2) / 0.9739065285 1717172007 7964012084 452 d0 /
785 data xgk( 3) / 0.9301574913 5570822600 1207180059 508 d0 /
786 data xgk( 4) / 0.8650633666 8898451073 2096688423 493 d0 /
787 data xgk( 5) / 0.7808177265 8641689706 3717578345 042 d0 /
788 data xgk( 6) / 0.6794095682 9902440623 4327365114 874 d0 /
789 data xgk( 7) / 0.5627571346 6860468333 9000099272 694 d0 /
790 data xgk( 8) / 0.4333953941 2924719079 9265943165 784 d0 /
791 data xgk( 9) / 0.2943928627 0146019813 1126603103 866 d0 /
792 data xgk( 10) / 0.1488743389 8163121088 4826001129 720 d0 /
793 data xgk( 11) / 0.0000000000 0000000000 0000000000 000 d0 /
795 data wgk( 1) / 0.0116946388 6737187427 8064396062 192 d0 /
796 data wgk( 2) / 0.0325581623 0796472747 8818972459 390 d0 /
797 data wgk( 3) / 0.0547558965 7435199603 1381300244 580 d0 /
798 data wgk( 4) / 0.0750396748 1091995276 7043140916 190 d0 /
799 data wgk( 5) / 0.0931254545 8369760553 5065465083 366 d0 /
800 data wgk( 6) / 0.1093871588 0229764189 9210590325 805 d0 /
801 data wgk( 7) / 0.1234919762 6206585107 7958109831 074 d0 /
802 data wgk( 8) / 0.1347092173 1147332592 8054001771 707 d0 /
803 data wgk( 9) / 0.1427759385 7706008079 7094273138 717 d0 /
804 data wgk( 10) / 0.1477391049 0133849137 4841515972 068 d0 /
805 data wgk( 11) / 0.1494455540 0291690566 4936468389 821 d0 /
831 centr = 0.5d+00*(a+b)
832 hlgth = 0.5d+00*(b-a)
844 absc = hlgth*xgk(jtw)
845 fval1 = f(centr-absc)
846 fval2 = f(centr+absc)
850 resg = resg+wg(j)*fsum
851 resk = resk+wgk(jtw)*fsum
852 resabs = resabs+wgk(jtw)*(dabs(fval1)+dabs(fval2))
856 absc = hlgth*xgk(jtwm1)
857 fval1 = f(centr-absc)
858 fval2 = f(centr+absc)
862 resk = resk+wgk(jtwm1)*fsum
863 resabs = resabs+wgk(jtwm1)*(dabs(fval1)+dabs(fval2))
866 resasc = wgk(11)*dabs(fc-reskh)
868 resasc = resasc+wgk(j)*(dabs(fv1(j)-reskh)+dabs(fv2(j)-reskh))
871 resabs = resabs*dhlgth
872 resasc = resasc*dhlgth
873 abserr = dabs((resk-resg)*hlgth)
874 if(resasc.ne.0.0d+00.and.abserr.ne.0.0d+00)
875 * abserr = resasc*dmin1(0.1d+01,(0.2d+03*abserr/resasc)**1.5d+00)
876 if(resabs.gt.uflow/(0.5d+02*epmach)) abserr = dmax1
877 * ((epmach*0.5d+02)*resabs,abserr)
880 subroutine dqk31(f,a,b,result,abserr,resabs,resasc)
932 double precision a,absc,abserr,b,centr,dabs,dhlgth,dmax1,dmin1,
933 * d1mach,epmach,f,fc,fsum,fval1,fval2,fv1,fv2,hlgth,resabs,resasc,
934 * resg,resk,reskh,result,uflow,wg,wgk,xgk
938 dimension fv1(15),fv2(15),xgk(16),wgk(16),wg(8)
959 data wg( 1) / 0.0307532419 9611726835 4628393577 204 d0 /
960 data wg( 2) / 0.0703660474 8810812470 9267416450 667 d0 /
961 data wg( 3) / 0.1071592204 6717193501 1869546685 869 d0 /
962 data wg( 4) / 0.1395706779 2615431444 7804794511 028 d0 /
963 data wg( 5) / 0.1662692058 1699393355 3200860481 209 d0 /
964 data wg( 6) / 0.1861610000 1556221102 6800561866 423 d0 /
965 data wg( 7) / 0.1984314853 2711157645 6118326443 839 d0 /
966 data wg( 8) / 0.2025782419 2556127288 0620199967 519 d0 /
968 data xgk( 1) / 0.9980022986 9339706028 5172840152 271 d0 /
969 data xgk( 2) / 0.9879925180 2048542848 9565718586 613 d0 /
970 data xgk( 3) / 0.9677390756 7913913425 7347978784 337 d0 /
971 data xgk( 4) / 0.9372733924 0070590430 7758947710 209 d0 /
972 data xgk( 5) / 0.8972645323 4408190088 2509656454 496 d0 /
973 data xgk( 6) / 0.8482065834 1042721620 0648320774 217 d0 /
974 data xgk( 7) / 0.7904185014 4246593296 7649294817 947 d0 /
975 data xgk( 8) / 0.7244177313 6017004741 6186054613 938 d0 /
976 data xgk( 9) / 0.6509967412 9741697053 3735895313 275 d0 /
977 data xgk( 10) / 0.5709721726 0853884753 7226737253 911 d0 /
978 data xgk( 11) / 0.4850818636 4023968069 3655740232 351 d0 /
979 data xgk( 12) / 0.3941513470 7756336989 7207370981 045 d0 /
980 data xgk( 13) / 0.2991800071 5316881216 6780024266 389 d0 /
981 data xgk( 14) / 0.2011940939 9743452230 0628303394 596 d0 /
982 data xgk( 15) / 0.1011420669 1871749902 7074231447 392 d0 /
983 data xgk( 16) / 0.0000000000 0000000000 0000000000 000 d0 /
985 data wgk( 1) / 0.0053774798 7292334898 7792051430 128 d0 /
986 data wgk( 2) / 0.0150079473 2931612253 8374763075 807 d0 /
987 data wgk( 3) / 0.0254608473 2671532018 6874001019 653 d0 /
988 data wgk( 4) / 0.0353463607 9137584622 2037948478 360 d0 /
989 data wgk( 5) / 0.0445897513 2476487660 8227299373 280 d0 /
990 data wgk( 6) / 0.0534815246 9092808726 5343147239 430 d0 /
991 data wgk( 7) / 0.0620095678 0067064028 5139230960 803 d0 /
992 data wgk( 8) / 0.0698541213 1872825870 9520077099 147 d0 /
993 data wgk( 9) / 0.0768496807 5772037889 4432777482 659 d0 /
994 data wgk( 10) / 0.0830805028 2313302103 8289247286 104 d0 /
995 data wgk( 11) / 0.0885644430 5621177064 7275443693 774 d0 /
996 data wgk( 12) / 0.0931265981 7082532122 5486872747 346 d0 /
997 data wgk( 13) / 0.0966427269 8362367850 5179907627 589 d0 /
998 data wgk( 14) / 0.0991735987 2179195933 2393173484 603 d0 /
999 data wgk( 15) / 0.1007698455 2387559504 4946662617 570 d0 /
1000 data wgk( 16) / 0.1013300070 1479154901 7374792767 493 d0 /
1022 centr = 0.5d+00*(a+b)
1023 hlgth = 0.5d+00*(b-a)
1024 dhlgth = dabs(hlgth)
1035 absc = hlgth*xgk(jtw)
1036 fval1 = f(centr-absc)
1037 fval2 = f(centr+absc)
1041 resg = resg+wg(j)*fsum
1042 resk = resk+wgk(jtw)*fsum
1043 resabs = resabs+wgk(jtw)*(dabs(fval1)+dabs(fval2))
1047 absc = hlgth*xgk(jtwm1)
1048 fval1 = f(centr-absc)
1049 fval2 = f(centr+absc)
1053 resk = resk+wgk(jtwm1)*fsum
1054 resabs = resabs+wgk(jtwm1)*(dabs(fval1)+dabs(fval2))
1056 reskh = resk*0.5d+00
1057 resasc = wgk(16)*dabs(fc-reskh)
1059 resasc = resasc+wgk(j)*(dabs(fv1(j)-reskh)+dabs(fv2(j)-reskh))
1062 resabs = resabs*dhlgth
1063 resasc = resasc*dhlgth
1064 abserr = dabs((resk-resg)*hlgth)
1065 if(resasc.ne.0.0d+00.and.abserr.ne.0.0d+00)
1066 * abserr = resasc*dmin1(0.1d+01,(0.2d+03*abserr/resasc)**1.5d+00)
1067 if(resabs.gt.uflow/(0.5d+02*epmach)) abserr = dmax1
1068 * ((epmach*0.5d+02)*resabs,abserr)
1071 subroutine dqk41(f,a,b,result,abserr,resabs,resasc)
1124 double precision a,absc,abserr,b,centr,dabs,dhlgth,dmax1,dmin1,
1125 * d1mach,epmach,f,fc,fsum,fval1,fval2,fv1,fv2,hlgth,resabs,resasc,
1126 * resg,resk,reskh,result,uflow,wg,wgk,xgk
1130 dimension fv1(20),fv2(20),xgk(21),wgk(21),wg(10)
1151 data wg( 1) / 0.0176140071 3915211831 1861962351 853 d0 /
1152 data wg( 2) / 0.0406014298 0038694133 1039952274 932 d0 /
1153 data wg( 3) / 0.0626720483 3410906356 9506535187 042 d0 /
1154 data wg( 4) / 0.0832767415 7670474872 4758143222 046 d0 /
1155 data wg( 5) / 0.1019301198 1724043503 6750135480 350 d0 /
1156 data wg( 6) / 0.1181945319 6151841731 2377377711 382 d0 /
1157 data wg( 7) / 0.1316886384 4917662689 8494499748 163 d0 /
1158 data wg( 8) / 0.1420961093 1838205132 9298325067 165 d0 /
1159 data wg( 9) / 0.1491729864 7260374678 7828737001 969 d0 /
1160 data wg( 10) / 0.1527533871 3072585069 8084331955 098 d0 /
1162 data xgk( 1) / 0.9988590315 8827766383 8315576545 863 d0 /
1163 data xgk( 2) / 0.9931285991 8509492478 6122388471 320 d0 /
1164 data xgk( 3) / 0.9815078774 5025025919 3342994720 217 d0 /
1165 data xgk( 4) / 0.9639719272 7791379126 7666131197 277 d0 /
1166 data xgk( 5) / 0.9408226338 3175475351 9982722212 443 d0 /
1167 data xgk( 6) / 0.9122344282 5132590586 7752441203 298 d0 /
1168 data xgk( 7) / 0.8782768112 5228197607 7442995113 078 d0 /
1169 data xgk( 8) / 0.8391169718 2221882339 4529061701 521 d0 /
1170 data xgk( 9) / 0.7950414288 3755119835 0638833272 788 d0 /
1171 data xgk( 10) / 0.7463319064 6015079261 4305070355 642 d0 /
1172 data xgk( 11) / 0.6932376563 3475138480 5490711845 932 d0 /
1173 data xgk( 12) / 0.6360536807 2651502545 2836696226 286 d0 /
1174 data xgk( 13) / 0.5751404468 1971031534 2946036586 425 d0 /
1175 data xgk( 14) / 0.5108670019 5082709800 4364050955 251 d0 /
1176 data xgk( 15) / 0.4435931752 3872510319 9992213492 640 d0 /
1177 data xgk( 16) / 0.3737060887 1541956067 2548177024 927 d0 /
1178 data xgk( 17) / 0.3016278681 1491300432 0555356858 592 d0 /
1179 data xgk( 18) / 0.2277858511 4164507808 0496195368 575 d0 /
1180 data xgk( 19) / 0.1526054652 4092267550 5220241022 678 d0 /
1181 data xgk( 20) / 0.0765265211 3349733375 4640409398 838 d0 /
1182 data xgk( 21) / 0.0000000000 0000000000 0000000000 000 d0 /
1184 data wgk( 1) / 0.0030735837 1852053150 1218293246 031 d0 /
1185 data wgk( 2) / 0.0086002698 5564294219 8661787950 102 d0 /
1186 data wgk( 3) / 0.0146261692 5697125298 3787960308 868 d0 /
1187 data wgk( 4) / 0.0203883734 6126652359 8010231432 755 d0 /
1188 data wgk( 5) / 0.0258821336 0495115883 4505067096 153 d0 /
1189 data wgk( 6) / 0.0312873067 7703279895 8543119323 801 d0 /
1190 data wgk( 7) / 0.0366001697 5820079803 0557240707 211 d0 /
1191 data wgk( 8) / 0.0416688733 2797368626 3788305936 895 d0 /
1192 data wgk( 9) / 0.0464348218 6749767472 0231880926 108 d0 /
1193 data wgk( 10) / 0.0509445739 2372869193 2707670050 345 d0 /
1194 data wgk( 11) / 0.0551951053 4828599474 4832372419 777 d0 /
1195 data wgk( 12) / 0.0591114008 8063957237 4967220648 594 d0 /
1196 data wgk( 13) / 0.0626532375 5478116802 5870122174 255 d0 /
1197 data wgk( 14) / 0.0658345971 3361842211 1563556969 398 d0 /
1198 data wgk( 15) / 0.0686486729 2852161934 5623411885 368 d0 /
1199 data wgk( 16) / 0.0710544235 5344406830 5790361723 210 d0 /
1200 data wgk( 17) / 0.0730306903 3278666749 5189417658 913 d0 /
1201 data wgk( 18) / 0.0745828754 0049918898 6581418362 488 d0 /
1202 data wgk( 19) / 0.0757044976 8455667465 9542775376 617 d0 /
1203 data wgk( 20) / 0.0763778676 7208073670 5502835038 061 d0 /
1204 data wgk( 21) / 0.0766007119 1799965644 5049901530 102 d0 /
1229 centr = 0.5d+00*(a+b)
1230 hlgth = 0.5d+00*(b-a)
1231 dhlgth = dabs(hlgth)
1242 absc = hlgth*xgk(jtw)
1243 fval1 = f(centr-absc)
1244 fval2 = f(centr+absc)
1248 resg = resg+wg(j)*fsum
1249 resk = resk+wgk(jtw)*fsum
1250 resabs = resabs+wgk(jtw)*(dabs(fval1)+dabs(fval2))
1254 absc = hlgth*xgk(jtwm1)
1255 fval1 = f(centr-absc)
1256 fval2 = f(centr+absc)
1260 resk = resk+wgk(jtwm1)*fsum
1261 resabs = resabs+wgk(jtwm1)*(dabs(fval1)+dabs(fval2))
1263 reskh = resk*0.5d+00
1264 resasc = wgk(21)*dabs(fc-reskh)
1266 resasc = resasc+wgk(j)*(dabs(fv1(j)-reskh)+dabs(fv2(j)-reskh))
1269 resabs = resabs*dhlgth
1270 resasc = resasc*dhlgth
1271 abserr = dabs((resk-resg)*hlgth)
1272 if(resasc.ne.0.0d+00.and.abserr.ne.0.d+00)
1273 * abserr = resasc*dmin1(0.1d+01,(0.2d+03*abserr/resasc)**1.5d+00)
1274 if(resabs.gt.uflow/(0.5d+02*epmach)) abserr = dmax1
1275 * ((epmach*0.5d+02)*resabs,abserr)
1278 subroutine dqk51(f,a,b,result,abserr,resabs,resasc)
1330 double precision a,absc,abserr,b,centr,dabs,dhlgth,dmax1,dmin1,
1331 * d1mach,epmach,f,fc,fsum,fval1,fval2,fv1,fv2,hlgth,resabs,resasc,
1332 * resg,resk,reskh,result,uflow,wg,wgk,xgk
1336 dimension fv1(25),fv2(25),xgk(26),wgk(26),wg(13)
1357 data wg( 1) / 0.0113937985 0102628794 7902964113 235 d0 /
1358 data wg( 2) / 0.0263549866 1503213726 1901815295 299 d0 /
1359 data wg( 3) / 0.0409391567 0130631265 5623487711 646 d0 /
1360 data wg( 4) / 0.0549046959 7583519192 5936891540 473 d0 /
1361 data wg( 5) / 0.0680383338 1235691720 7187185656 708 d0 /
1362 data wg( 6) / 0.0801407003 3500101801 3234959669 111 d0 /
1363 data wg( 7) / 0.0910282619 8296364981 1497220702 892 d0 /
1364 data wg( 8) / 0.1005359490 6705064420 2206890392 686 d0 /
1365 data wg( 9) / 0.1085196244 7426365311 6093957050 117 d0 /
1366 data wg( 10) / 0.1148582591 4571164833 9325545869 556 d0 /
1367 data wg( 11) / 0.1194557635 3578477222 8178126512 901 d0 /
1368 data wg( 12) / 0.1222424429 9031004168 8959518945 852 d0 /
1369 data wg( 13) / 0.1231760537 2671545120 3902873079 050 d0 /
1371 data xgk( 1) / 0.9992621049 9260983419 3457486540 341 d0 /
1372 data xgk( 2) / 0.9955569697 9049809790 8784946893 902 d0 /
1373 data xgk( 3) / 0.9880357945 3407724763 7331014577 406 d0 /
1374 data xgk( 4) / 0.9766639214 5951751149 8315386479 594 d0 /
1375 data xgk( 5) / 0.9616149864 2584251241 8130033660 167 d0 /
1376 data xgk( 6) / 0.9429745712 2897433941 4011169658 471 d0 /
1377 data xgk( 7) / 0.9207471152 8170156174 6346084546 331 d0 /
1378 data xgk( 8) / 0.8949919978 7827536885 1042006782 805 d0 /
1379 data xgk( 9) / 0.8658470652 9327559544 8996969588 340 d0 /
1380 data xgk( 10) / 0.8334426287 6083400142 1021108693 570 d0 /
1381 data xgk( 11) / 0.7978737979 9850005941 0410904994 307 d0 /
1382 data xgk( 12) / 0.7592592630 3735763057 7282865204 361 d0 /
1383 data xgk( 13) / 0.7177664068 1308438818 6654079773 298 d0 /
1384 data xgk( 14) / 0.6735663684 7346836448 5120633247 622 d0 /
1385 data xgk( 15) / 0.6268100990 1031741278 8122681624 518 d0 /
1386 data xgk( 16) / 0.5776629302 4122296772 3689841612 654 d0 /
1387 data xgk( 17) / 0.5263252843 3471918259 9623778158 010 d0 /
1388 data xgk( 18) / 0.4730027314 4571496052 2182115009 192 d0 /
1389 data xgk( 19) / 0.4178853821 9303774885 1814394594 572 d0 /
1390 data xgk( 20) / 0.3611723058 0938783773 5821730127 641 d0 /
1391 data xgk( 21) / 0.3030895389 3110783016 7478909980 339 d0 /
1392 data xgk( 22) / 0.2438668837 2098843204 5190362797 452 d0 /
1393 data xgk( 23) / 0.1837189394 2104889201 5969888759 528 d0 /
1394 data xgk( 24) / 0.1228646926 1071039638 7359818808 037 d0 /
1395 data xgk( 25) / 0.0615444830 0568507888 6546392366 797 d0 /
1396 data xgk( 26) / 0.0000000000 0000000000 0000000000 000 d0 /
1398 data wgk( 1) / 0.0019873838 9233031592 6507851882 843 d0 /
1399 data wgk( 2) / 0.0055619321 3535671375 8040236901 066 d0 /
1400 data wgk( 3) / 0.0094739733 8617415160 7207710523 655 d0 /
1401 data wgk( 4) / 0.0132362291 9557167481 3656405846 976 d0 /
1402 data wgk( 5) / 0.0168478177 0912829823 1516667536 336 d0 /
1403 data wgk( 6) / 0.0204353711 4588283545 6568292235 939 d0 /
1404 data wgk( 7) / 0.0240099456 0695321622 0092489164 881 d0 /
1405 data wgk( 8) / 0.0274753175 8785173780 2948455517 811 d0 /
1406 data wgk( 9) / 0.0307923001 6738748889 1109020215 229 d0 /
1407 data wgk( 10) / 0.0340021302 7432933783 6748795229 551 d0 /
1408 data wgk( 11) / 0.0371162714 8341554356 0330625367 620 d0 /
1409 data wgk( 12) / 0.0400838255 0403238207 4839284467 076 d0 /
1410 data wgk( 13) / 0.0428728450 2017004947 6895792439 495 d0 /
1411 data wgk( 14) / 0.0455029130 4992178890 9870584752 660 d0 /
1412 data wgk( 15) / 0.0479825371 3883671390 6392255756 915 d0 /
1413 data wgk( 16) / 0.0502776790 8071567196 3325259433 440 d0 /
1414 data wgk( 17) / 0.0523628858 0640747586 4366712137 873 d0 /
1415 data wgk( 18) / 0.0542511298 8854549014 4543370459 876 d0 /
1416 data wgk( 19) / 0.0559508112 2041231730 8240686382 747 d0 /
1417 data wgk( 20) / 0.0574371163 6156783285 3582693939 506 d0 /
1418 data wgk( 21) / 0.0586896800 2239420796 1974175856 788 d0 /
1419 data wgk( 22) / 0.0597203403 2417405997 9099291932 562 d0 /
1420 data wgk( 23) / 0.0605394553 7604586294 5360267517 565 d0 /
1421 data wgk( 24) / 0.0611285097 1705304830 5859030416 293 d0 /
1422 data wgk( 25) / 0.0614711898 7142531666 1544131965 264 d0 /
1424 data wgk( 26) / 0.0615808180 6783293507 8759824240 066 d0 /
1449 centr = 0.5d+00*(a+b)
1450 hlgth = 0.5d+00*(b-a)
1451 dhlgth = dabs(hlgth)
1462 absc = hlgth*xgk(jtw)
1463 fval1 = f(centr-absc)
1464 fval2 = f(centr+absc)
1468 resg = resg+wg(j)*fsum
1469 resk = resk+wgk(jtw)*fsum
1470 resabs = resabs+wgk(jtw)*(dabs(fval1)+dabs(fval2))
1474 absc = hlgth*xgk(jtwm1)
1475 fval1 = f(centr-absc)
1476 fval2 = f(centr+absc)
1480 resk = resk+wgk(jtwm1)*fsum
1481 resabs = resabs+wgk(jtwm1)*(dabs(fval1)+dabs(fval2))
1483 reskh = resk*0.5d+00
1484 resasc = wgk(26)*dabs(fc-reskh)
1486 resasc = resasc+wgk(j)*(dabs(fv1(j)-reskh)+dabs(fv2(j)-reskh))
1489 resabs = resabs*dhlgth
1490 resasc = resasc*dhlgth
1491 abserr = dabs((resk-resg)*hlgth)
1492 if(resasc.ne.0.0d+00.and.abserr.ne.0.0d+00)
1493 * abserr = resasc*dmin1(0.1d+01,(0.2d+03*abserr/resasc)**1.5d+00)
1494 if(resabs.gt.uflow/(0.5d+02*epmach)) abserr = dmax1
1495 * ((epmach*0.5d+02)*resabs,abserr)
1498 subroutine dqk61(f,a,b,result,abserr,resabs,resasc)
1551 double precision a,dabsc,abserr,b,centr,dabs,dhlgth,dmax1,dmin1,
1552 * d1mach,epmach,f,fc,fsum,fval1,fval2,fv1,fv2,hlgth,resabs,resasc,
1553 * resg,resk,reskh,result,uflow,wg,wgk,xgk
1557 dimension fv1(30),fv2(30),xgk(31),wgk(31),wg(15)
1578 data wg( 1) / 0.0079681924 9616660561 5465883474 674 d0 /
1579 data wg( 2) / 0.0184664683 1109095914 2302131912 047 d0 /
1580 data wg( 3) / 0.0287847078 8332336934 9719179611 292 d0 /
1581 data wg( 4) / 0.0387991925 6962704959 6801936446 348 d0 /
1582 data wg( 5) / 0.0484026728 3059405290 2938140422 808 d0 /
1583 data wg( 6) / 0.0574931562 1761906648 1721689402 056 d0 /
1584 data wg( 7) / 0.0659742298 8218049512 8128515115 962 d0 /
1585 data wg( 8) / 0.0737559747 3770520626 8243850022 191 d0 /
1586 data wg( 9) / 0.0807558952 2942021535 4694938460 530 d0 /
1587 data wg( 10) / 0.0868997872 0108297980 2387530715 126 d0 /
1588 data wg( 11) / 0.0921225222 3778612871 7632707087 619 d0 /
1589 data wg( 12) / 0.0963687371 7464425963 9468626351 810 d0 /
1590 data wg( 13) / 0.0995934205 8679526706 2780282103 569 d0 /
1591 data wg( 14) / 0.1017623897 4840550459 6428952168 554 d0 /
1592 data wg( 15) / 0.1028526528 9355884034 1285636705 415 d0 /
1594 data xgk( 1) / 0.9994844100 5049063757 1325895705 811 d0 /
1595 data xgk( 2) / 0.9968934840 7464954027 1630050918 695 d0 /
1596 data xgk( 3) / 0.9916309968 7040459485 8628366109 486 d0 /
1597 data xgk( 4) / 0.9836681232 7974720997 0032581605 663 d0 /
1598 data xgk( 5) / 0.9731163225 0112626837 4693868423 707 d0 /
1599 data xgk( 6) / 0.9600218649 6830751221 6871025581 798 d0 /
1600 data xgk( 7) / 0.9443744447 4855997941 5831324037 439 d0 /
1601 data xgk( 8) / 0.9262000474 2927432587 9324277080 474 d0 /
1602 data xgk( 9) / 0.9055733076 9990779854 6522558925 958 d0 /
1603 data xgk( 10) / 0.8825605357 9205268154 3116462530 226 d0 /
1604 data xgk( 11) / 0.8572052335 4606109895 8658510658 944 d0 /
1605 data xgk( 12) / 0.8295657623 8276839744 2898119732 502 d0 /
1606 data xgk( 13) / 0.7997278358 2183908301 3668942322 683 d0 /
1607 data xgk( 14) / 0.7677774321 0482619491 7977340974 503 d0 /
1608 data xgk( 15) / 0.7337900624 5322680472 6171131369 528 d0 /
1609 data xgk( 16) / 0.6978504947 9331579693 2292388026 640 d0 /
1610 data xgk( 17) / 0.6600610641 2662696137 0053668149 271 d0 /
1611 data xgk( 18) / 0.6205261829 8924286114 0477556431 189 d0 /
1612 data xgk( 19) / 0.5793452358 2636169175 6024932172 540 d0 /
1613 data xgk( 20) / 0.5366241481 4201989926 4169793311 073 d0 /
1614 data xgk( 21) / 0.4924804678 6177857499 3693061207 709 d0 /
1615 data xgk( 22) / 0.4470337695 3808917678 0609900322 854 d0 /
1616 data xgk( 23) / 0.4004012548 3039439253 5476211542 661 d0 /
1617 data xgk( 24) / 0.3527047255 3087811347 1037207089 374 d0 /
1618 data xgk( 25) / 0.3040732022 7362507737 2677107199 257 d0 /
1619 data xgk( 26) / 0.2546369261 6788984643 9805129817 805 d0 /
1620 data xgk( 27) / 0.2045251166 8230989143 8957671002 025 d0 /
1621 data xgk( 28) / 0.1538699136 0858354696 3794672743 256 d0 /
1622 data xgk( 29) / 0.1028069379 6673703014 7096751318 001 d0 /
1623 data xgk( 30) / 0.0514718425 5531769583 3025213166 723 d0 /
1624 data xgk( 31) / 0.0000000000 0000000000 0000000000 000 d0 /
1626 data wgk( 1) / 0.0013890136 9867700762 4551591226 760 d0 /
1627 data wgk( 2) / 0.0038904611 2709988405 1267201844 516 d0 /
1628 data wgk( 3) / 0.0066307039 1593129217 3319826369 750 d0 /
1629 data wgk( 4) / 0.0092732796 5951776342 8441146892 024 d0 /
1630 data wgk( 5) / 0.0118230152 5349634174 2232898853 251 d0 /
1631 data wgk( 6) / 0.0143697295 0704580481 2451432443 580 d0 /
1632 data wgk( 7) / 0.0169208891 8905327262 7572289420 322 d0 /
1633 data wgk( 8) / 0.0194141411 9394238117 3408951050 128 d0 /
1634 data wgk( 9) / 0.0218280358 2160919229 7167485738 339 d0 /
1635 data wgk( 10) / 0.0241911620 7808060136 5686370725 232 d0 /
1636 data wgk( 11) / 0.0265099548 8233310161 0601709335 075 d0 /
1637 data wgk( 12) / 0.0287540487 6504129284 3978785354 334 d0 /
1638 data wgk( 13) / 0.0309072575 6238776247 2884252943 092 d0 /
1639 data wgk( 14) / 0.0329814470 5748372603 1814191016 854 d0 /
1640 data wgk( 15) / 0.0349793380 2806002413 7499670731 468 d0 /
1641 data wgk( 16) / 0.0368823646 5182122922 3911065617 136 d0 /
1642 data wgk( 17) / 0.0386789456 2472759295 0348651532 281 d0 /
1643 data wgk( 18) / 0.0403745389 5153595911 1995279752 468 d0 /
1644 data wgk( 19) / 0.0419698102 1516424614 7147541285 970 d0 /
1645 data wgk( 20) / 0.0434525397 0135606931 6831728117 073 d0 /
1646 data wgk( 21) / 0.0448148001 3316266319 2355551616 723 d0 /
1647 data wgk( 22) / 0.0460592382 7100698811 6271735559 374 d0 /
1648 data wgk( 23) / 0.0471855465 6929915394 5261478181 099 d0 /
1649 data wgk( 24) / 0.0481858617 5708712914 0779492298 305 d0 /
1650 data wgk( 25) / 0.0490554345 5502977888 7528165367 238 d0 /
1651 data wgk( 26) / 0.0497956834 2707420635 7811569379 942 d0 /
1652 data wgk( 27) / 0.0504059214 0278234684 0893085653 585 d0 /
1653 data wgk( 28) / 0.0508817958 9874960649 2297473049 805 d0 /
1654 data wgk( 29) / 0.0512215478 4925877217 0656282604 944 d0 /
1655 data wgk( 30) / 0.0514261285 3745902593 3862879215 781 d0 /
1656 data wgk( 31) / 0.0514947294 2945156755 8340433647 099 d0 /
1679 centr = 0.5d+00*(b+a)
1680 hlgth = 0.5d+00*(b-a)
1681 dhlgth = dabs(hlgth)
1693 dabsc = hlgth*xgk(jtw)
1694 fval1 = f(centr-dabsc)
1695 fval2 = f(centr+dabsc)
1699 resg = resg+wg(j)*fsum
1700 resk = resk+wgk(jtw)*fsum
1701 resabs = resabs+wgk(jtw)*(dabs(fval1)+dabs(fval2))
1705 dabsc = hlgth*xgk(jtwm1)
1706 fval1 = f(centr-dabsc)
1707 fval2 = f(centr+dabsc)
1711 resk = resk+wgk(jtwm1)*fsum
1712 resabs = resabs+wgk(jtwm1)*(dabs(fval1)+dabs(fval2))
1714 reskh = resk*0.5d+00
1715 resasc = wgk(31)*dabs(fc-reskh)
1717 resasc = resasc+wgk(j)*(dabs(fv1(j)-reskh)+dabs(fv2(j)-reskh))
1720 resabs = resabs*dhlgth
1721 resasc = resasc*dhlgth
1722 abserr = dabs((resk-resg)*hlgth)
1723 if(resasc.ne.0.0d+00.and.abserr.ne.0.0d+00)
1724 * abserr = resasc*dmin1(0.1d+01,(0.2d+03*abserr/resasc)**1.5d+00)
1725 if(resabs.gt.uflow/(0.5d+02*epmach)) abserr = dmax1
1726 * ((epmach*0.5d+02)*resabs,abserr)
1729 subroutine dqpsrt(limit,last,maxerr,ermax,elist,iord,nrmax)
1783 double precision elist,ermax,errmax,errmin
1784 integer i,ibeg,ido,iord,isucc,j,jbnd,jupbn,k,last,limit,maxerr,
1786 dimension elist(last),iord(last)
1792 if(last.gt.2)
go to 10
1802 10 errmax = elist(maxerr)
1803 if(nrmax.eq.1)
go to 30
1806 isucc = iord(nrmax-1)
1808 if(errmax.le.elist(isucc))
go to 30
1818 if(last.gt.(limit/2+2)) jupbn = limit+3-last
1819 errmin = elist(last)
1826 if(ibeg.gt.jbnd)
go to 50
1830 if(errmax.ge.elist(isucc))
go to 60
1833 50 iord(jbnd) = maxerr
1839 60 iord(i-1) = maxerr
1844 if(errmin.lt.elist(isucc))
go to 80
1854 90 maxerr = iord(nrmax)
1855 ermax = elist(maxerr)
1859 DOUBLE PRECISION FUNCTION d1mach(I)
1874 INTEGER SC, CRAY1(38), J
1875 COMMON /d9mach/ cray1
1876 SAVE small, large, right, diver, log10, sc
1877 DOUBLE PRECISION DMACH(5)
1878 equivalence(dmach(1),small(1))
1879 equivalence(dmach(2),large(1))
1880 equivalence(dmach(3),right(1))
1881 equivalence(dmach(4),diver(1))
1882 equivalence(dmach(5),log10(1))
1919 IF (sc .NE. 987)
THEN
1921 IF ( small(1) .EQ. 1117925532
1922 * .AND. small(2) .EQ. -448790528)
THEN
1926 large(1) = 2146435071
1928 right(1) = 1017118720
1930 diver(1) = 1018167296
1932 log10(1) = 1070810131
1933 log10(2) = 1352628735
1934 ELSE IF ( small(2) .EQ. 1117925532
1935 * .AND. small(1) .EQ. -448790528)
THEN
1939 large(2) = 2146435071
1941 right(2) = 1017118720
1943 diver(2) = 1018167296
1945 log10(2) = 1070810131
1946 log10(1) = 1352628735
1947 ELSE IF ( small(1) .EQ. -2065213935
1948 * .AND. small(2) .EQ. 10752)
THEN
1958 log10(1) = 546979738
1959 log10(2) = -805796613
1960 ELSE IF ( small(1) .EQ. 1267827943
1961 * .AND. small(2) .EQ. 704643072)
THEN
1965 large(1) = 2147483647
1967 right(1) = 856686592
1969 diver(1) = 873463808
1971 log10(1) = 1091781651
1972 log10(2) = 1352628735
1973 ELSE IF ( small(1) .EQ. 1120022684
1974 * .AND. small(2) .EQ. -448790528)
THEN
1978 large(1) = 2147483647
1980 right(1) = 1019215872
1982 diver(1) = 1020264448
1984 log10(1) = 1072907283
1985 log10(2) = 1352628735
1986 ELSE IF ( small(1) .EQ. 815547074
1987 * .AND. small(2) .EQ. 58688)
THEN
1997 log10(1) = 1142112243
1998 log10(2) = 2046775455
2000 dmach(2) = 1.d27 + 1
2002 large(2) = large(2) - right(2)
2003 IF (large(2) .EQ. 64 .AND. small(2) .EQ. 0)
THEN
2006 cray1(j+1) = cray1(j) + cray1(j)
2008 cray1(22) = cray1(21) + 321322
2010 cray1(j+1) = cray1(j) + cray1(j)
2012 IF (cray1(38) .EQ. small(1))
THEN
2014 CALL i1mcry(small(1), j, 8285, 8388608, 0)
2016 CALL i1mcry(large(1), j, 24574, 16777215, 16777215)
2017 CALL i1mcry(large(2), j, 0, 16777215, 16777214)
2018 CALL i1mcry(right(1), j, 16291, 8388608, 0)
2020 CALL i1mcry(diver(1), j, 16292, 8388608, 0)
2022 CALL i1mcry(log10(1), j, 16383, 10100890, 8715215)
2023 CALL i1mcry(log10(2), j, 0, 16226447, 9001388)
2036 IF (dmach(4) .GE. 1.0d0) stop 778
2037 IF (i .LT. 1 .OR. i .GT. 5)
THEN
2038 WRITE(*,*)
'D1MACH(I): I =',i,
' is out of bounds.'
2043 9000
FORMAT(/
' Adjust D1MACH by uncommenting data statements'/
2044 *
' appropriate for your machine.')
2047 SUBROUTINE i1mcry(A, A1, B, C, D)
2049 INTEGER A, A1, B, C, D