C ***** BENCHMARK SP1111 ***** IMPLICIT DOUBLE PRECISION(A-H,O-Z) PARAMETER NN=5000 COMMON/MISC/NGANGB COMMON/GEOM/RAB,RABSQ COMMON/PGEOM/GP(36),EP(36),DP00P(36),DP01P(36),DP10P(36),APP(36), $ BPP(36) COMMON/GINF/GCD COMMON/EABECD/EAB,ECD COMMON/QGEOM/AQZ,QPERP,QPERP2 COMMON/CONST/CONP(36) COMMON/MAXC/ISMLP(36),ISMLQ,ISML COMMON/H/H(160) COMMON/PICON/PITO52,PIDIV4 COMMON/IO/IN,IOUT,IPUNCH C DATA TWO/2.0D0/,TWOPT5/2.5D0/,FOUR/4.0D0/,TEN/10.0D0/,F326/3.26077 $D0/,F409/4.091326D0/,ONE/1.0D0/,PT2/0.2D0/,PT1/0.1D0/ DATA ZERO/0.0D0/ DATA F160/160.0D0/ DATA CHEK/4.590475501D0/,THR/1.0D-08/ C 2001 FORMAT(' FOR I = ',I6, ' SUM = ',D20.13) 2002 FORMAT(' THERE WERE A TOTAL OF ',I6,' FAILURES.') C IN=5 IOUT=6 IPUNCH=7 C C CALL HEADER('BENCHMARK SP1111',16) CALL TIMRB C C PI=FOUR*DATAN(ONE) PITO52=TWO*(PI**TWOPT5) PIDIV4=PI/FOUR C NGANGB=36 C RABSQ=PI/TEN RAB=DSQRT(RABSQ) C EXXA=F326 EXXB=F409 GCD=EXXA+EXXB C ECD=ONE/GCD C DO 10 I=1,36 GP(I)=TWO*F326 EP(I)=ONE/GP(I) DP00P(I)=F326/TEN DP01P(I)=F326/TEN DP10P(I)=F326/TEN APP(I)=GP(I) BPP(I)=EP(I) CONP(I)=PITO52*GP(I) 10 ISMLP(I)=0 C AQZ=PT1*PI QPERP2=PT2*PI QPERP=DSQRT(QPERP2) C ISMLQ=0 C CALL GAMGEN(0) C IFAIL=0 DO 40 I=1,NN DO 20 J=1,160 20 H(J)=ZERO CALL SP1111(I) SUM=ZERO DO 30 J=1,160 30 SUM=SUM+DABS(H(J)) SUM=SUM/F160 IF(DABS(SUM-CHEK).LE.THR)GO TO 40 WRITE(IOUT,2001)I,SUM IFAIL=IFAIL+1 40 CONTINUE C CALL TIMRE C C WRITE(IOUT,2002)IFAIL C STOP END SUBROUTINE SP1111(IQXX) 0010 IMPLICIT DOUBLE PRECISION(A-H,O-Z) 0020 C 0030 C -------------------------- 0040 C GAUSSIAN 74 (QCPE VERSION) 0050 C DECEMBER 1974 0060 C UNIVAC 1108/CDC 7600 0070 C -------------------------- 0080 C 0090 C 0100 COMMON/MISC/NGANGB COMMON/GEOM/RAB,RABSQ COMMON/PGEOM/GP(36),EP(36),DP00P(36),DP01P(36),DP10P(36) *,APP(36),BPP(36) 0190 COMMON/PQGEOM/AP,BP COMMON/GINF/GCD COMMON/EABECD/EAB,ECD 0230 COMMON/DPQ/DP00,DP01,DP10 COMMON/QGEOM/AQZ,QPERP,QPERP2 COMMON/CONST/CONP(36) COMMON/MAXC/ISMLP(36),ISMLQ,ISML COMMON/ASTORE/QQ,THETA,N 0290 COMMON/TABLE/ 0300 *AA(400),BA(400),CA(400), 0310 *AB(400),BB(400),CB(400), 0320 *AC(400),BC(400),CC(400), 0330 *AD(400),BD(400),CD(400), 0340 *AE(400),BE(400),CE(400) 0350 COMMON/AUXVAR/AUXVAR,VAR1,VAR2 0360 COMMON/H/ 0370 *H0000,H0001,H0002,H0003,H0011,H0012,H0013,H0022,H0023,H0033, 0380 *H0100,H0101,H0102,H0103,H0111,H0112,H0113,H0122,H0123,H0133, 0390 *H0200,H0201,H0202,H0203,H0211,H0212,H0213,H0222,H0223,H0233, 0400 *H0300,H0301,H0302,H0303,H0311,H0312,H0313,H0322,H0323,H0333, 0410 *H1000,H1001,H1002,H1003,H1011,H1012,H1013,H1022,H1023,H1033, 0420 *H1100,H1101,H1102,H1103,H1111,H1112,H1113,H1122,H1123,H1133, 0430 *H1200,H1201,H1202,H1203,H1211,H1212,H1213,H1222,H1223,H1233, 0440 *H1300,H1301,H1302,H1303,H1311,H1312,H1313,H1322,H1323,H1333, 0450 *H2000,H2001,H2002,H2003,H2011,H2012,H2013,H2022,H2023,H2033, 0460 *H2100,H2101,H2102,H2103,H2111,H2112,H2113,H2122,H2123,H2133, 0470 *H2200,H2201,H2202,H2203,H2211,H2212,H2213,H2222,H2223,H2233, 0480 *H2300,H2301,H2302,H2303,H2311,H2312,H2313,H2322,H2323,H2333, 0490 *H3000,H3001,H3002,H3003,H3011,H3012,H3013,H3022,H3023,H3033, 0500 *H3100,H3101,H3102,H3103,H3111,H3112,H3113,H3122,H3123,H3133, 0510 *H3200,H3201,H3202,H3203,H3211,H3212,H3213,H3222,H3223,H3233, 0520 *H3300,H3301,H3302,H3303,H3311,H3312,H3313,H3322,H3323,H3333 0530 COMMON/PICON/PITO52,PIDIV4 C 0540 DATA ZERO/0.0D0/,PT25/0.25D0/,PT5/0.5D0/,ONE/1.0D0/,ONEPT5/1.5D0/, 0550 *TWO/2.0D0/,TWOPT5/2.5D0/,THREE/3.0D0/,THRPT5/3.5D0/,TWENTY/20.0D0/ 0560 C 0570 X1=ZERO 0580 X2=ZERO 0590 X3=ZERO 0600 X4=ZERO 0610 X5=ZERO 0620 X6=ZERO 0630 Y1=ZERO 0640 Y2=ZERO 0650 Y3=ZERO 0660 Y4=ZERO 0670 Y5=ZERO 0680 Y6=ZERO 0690 Z1=ZERO 0700 Z2=ZERO 0710 Z3=ZERO 0720 Z4=ZERO 0730 Z5=ZERO 0740 Z6=ZERO 0750 Z7=ZERO 0760 Z8=ZERO 0770 Z9=ZERO 0780 V1=ZERO 0790 V2=ZERO 0800 V3=ZERO 0810 V4=ZERO 0820 V5=ZERO 0830 V6=ZERO 0840 W1=ZERO 0850 W2=ZERO 0860 W3=ZERO 0870 W4=ZERO 0880 W5=ZERO 0890 W6=ZERO 0900 W7=ZERO 0910 W8=ZERO 0920 W9=ZERO 0930 S1=ZERO 0940 S2=ZERO 0950 S3=ZERO 0960 S4=ZERO 0970 S5=ZERO 0980 S6=ZERO 0990 S7=ZERO 1000 S8=ZERO 1010 S9=ZERO 1020 S10=ZERO 1030 S11=ZERO 1040 S12=ZERO 1050 S13=ZERO 1060 S14=ZERO 1070 T1=ZERO 1080 T2=ZERO 1090 T3=ZERO 1100 T4=ZERO 1110 T5=ZERO 1120 T6=ZERO 1130 T7=ZERO 1140 T8=ZERO 1150 T9=ZERO 1160 T10=ZERO 1170 T11=ZERO 1180 T12=ZERO 1190 T13=ZERO 1200 T14=ZERO 1210 C1=ZERO 1220 C2=ZERO 1230 C3=ZERO 1240 C4=ZERO 1250 C5=ZERO 1260 C6=ZERO 1270 DO 90 IND=1,NGANGB 1280 ISML=ISMLQ+ISMLP(IND) IF(ISML-2)10,90,90 1300 10 IF(ISML-1)20,30,30 20 AUXVAR=VAR1 1320 GO TO 40 1330 30 AUXVAR=VAR2 1340 40 EAB=EP(IND) 1350 DP00=DP00P(IND) 1360 DP01=DP01P(IND) 1370 DP10=DP10P(IND) 1380 AP=APP(IND) 1390 BP=BPP(IND) 1400 PQAB=AQZ-AP 1410 PQAB2=PQAB*PQAB 1420 G=ONE/(EAB+ECD) 1430 X=G*(QPERP2+PQAB2) 1440 IF(X-AUXVAR)60,60,50 1450 50 F0=CONP(IND)*DSQRT(PIDIV4/(X*(GP(IND)+GCD))) GTX=G/X 1470 F1=PT5*F0*GTX 1480 F2=ONEPT5*F1*GTX 1490 F3=TWOPT5*F2*GTX 1500 F4=THRPT5*F3*GTX 1510 GO TO 70 1520 60 Y=CONP(IND)/DSQRT(GP(IND)+GCD) 1530 GY=G*Y 1540 GGY=G*GY 1550 GGGY=G*GGY 1560 QQ=X*TWENTY 1570 THETA=QQ-DINT(QQ) 1580 N=QQ-THETA 1590 THETA2=THETA*(THETA-ONE) 1600 THETA3=THETA2*(THETA-TWO) 1610 THETA4=THETA2*(THETA+ONE) 1620 F0=(AA(N+1)+THETA*BA(N+1)-THETA3*CA(N+1)+THETA4*CA(N+2))*Y 1630 F1=(AB(N+1)+THETA*BB(N+1)-THETA3*CB(N+1)+THETA4*CB(N+2))*GY 1640 F2=(AC(N+1)+THETA*BC(N+1)-THETA3*CC(N+1)+THETA4*CC(N+2))*GGY 1650 F3=(AD(N+1)+THETA*BD(N+1)-THETA3*CD(N+1)+THETA4*CD(N+2))*GGGY 1660 F4=(AE(N+1)+THETA*BE(N+1)-THETA3*CE(N+1)+THETA4*CE(N+2))*GGGY*G 1670 70 APBP=AP*BP 1680 EAB2=EAB*EAB 1690 BPDP01=BP*DP01 1700 APDP10=AP*DP10 1710 EDP01=EAB*DP01 1720 EDP10=EAB*DP10 1730 F1PQAB=F1*PQAB 1740 F2PQAB=F2*PQAB 1750 F3PQAB=F3*PQAB 1760 F4PQAB=F4*PQAB 1770 F1PQA2=F1*PQAB2 1780 F2PQA2=F2*PQAB2 1790 F3PQA2=F3*PQAB2 1800 F4PQA2=F4*PQAB2 1810 F2PQA3=F2PQA2*PQAB 1820 F3PQA3=F3PQA2*PQAB 1830 F4PQA3=F4PQA2*PQAB 1840 X1=X1+F0 *DP00 1850 X2=X2+F1 *DP00 1860 X3=X3+F2 *DP00 1870 X4=X4+F1PQAB*DP00 1880 X5=X5+F2PQAB*DP00 1890 X6=X6+F2PQA2*DP00 1900 Z1=Z1+F1 *EDP01 1910 Z2=Z2+F2 *EDP01 1920 Z3=Z3+F3 *EDP01 1930 Z4=Z4+F1PQAB *EDP01 Z5=Z5+F2PQAB *EDP01 1950 Z6=Z6+F3PQAB *EDP01 1960 Z7=Z7+F2PQA2 *EDP01 1970 Z8=Z8+F3PQA2 *EDP01 1980 Z9=Z9+F3PQA3 *EDP01 1990 W1=W1+F1 *EDP10 2000 W2=W2+F2 *EDP10 2010 W3=W3+F3 *EDP10 2020 W4=W4+F1PQAB *EDP10 2030 W5=W5+F2PQAB *EDP10 2040 W6=W6+F3PQAB *EDP10 2050 W7=W7+F2PQA2 *EDP10 2060 W8=W8+F3PQA2 *EDP10 2070 W9=W9+F3PQA3 *EDP10 2080 S1=S1+F0 *EAB 2090 S2=S2+F1 *EAB 2100 S3=S3+F2 *EAB 2110 S4=S4+F3 *EAB 2120 S6=S6+F1PQAB*EAB 2130 S7=S7+F2PQAB*EAB 2140 S8=S8+F3PQAB*EAB 2150 S9=S9+F1PQA2*EAB 2160 S10=S10+F2PQA2*EAB 2170 S11=S11+F3PQA2*EAB 2180 S12=S12+F2PQA3*EAB 2190 S13=S13+F3PQA3*EAB 2200 S14=S14+F3PQA3*PQAB*EAB 2210 T1=T1+F0 *EAB2 2220 T2=T2+F1 *EAB2 2230 T3=T3+F2 *EAB2 2240 T4=T4+F3 *EAB2 2250 T5=T5+F4 *EAB2 2260 T6=T6+F2PQAB*EAB2 2270 T7=T7+F3PQAB*EAB2 2280 T8=T8+F4PQAB*EAB2 2290 T9=T9+F2PQA2*EAB2 2300 T10=T10+F3PQA2*EAB2 2310 T11=T11+F4PQA2*EAB2 2320 T12=T12+F3PQA3*EAB2 2330 T13=T13+F4PQA3*EAB2 2340 T14=T14+F4PQA3*PQAB*EAB2 2350 IF (RABSQ) 80,90,80 2360 80 CONTINUE 2370 Y1=Y1+F0 *BPDP01 2380 Y2=Y2+F1 *BPDP01 2390 Y3=Y3+F2 *BPDP01 2400 Y4=Y4+F1PQAB*BPDP01 2410 Y5=Y5+F2PQAB*BPDP01 2420 Y6=Y6+F2PQA2*BPDP01 2430 V1=V1+F0 *APDP10 2440 V2=V2+F1 *APDP10 2450 V3=V3+F2 *APDP10 2460 V4=V4+F1PQAB*APDP10 2470 V5=V5+F2PQAB*APDP10 2480 V6=V6+F2PQA2*APDP10 2490 C1=C1+F0 *APBP 2500 C2=C2+F1 *APBP 2510 C3=C3+F2 *APBP 2520 C4=C4+F1PQAB*APBP 2530 C5=C5+F2PQAB*APBP 2540 C6=C6+F2PQA2*APBP 2550 90 CONTINUE 2560 A1=AQZ*S2-S6 2570 A2=AQZ*S3-S7 2580 A3=AQZ*S4-S8 2590 A4=AQZ*S6-S9 2600 A5=AQZ*S7-S10 2610 A6=AQZ*S8-S11 2620 A8=AQZ*S10-S12 2630 A9=AQZ*S11-S13 2640 A10=AQZ*S13-S14 2650 BQZ=AQZ-RAB 2660 B1=BQZ*S2-S6 2670 B2=BQZ*S3-S7 2680 B3=BQZ*S4-S8 2690 B4=BQZ*S6-S9 2700 B5=BQZ*S7-S10 2710 B6=BQZ*S8-S11 2720 B8=BQZ*S10-S12 2730 B9=BQZ*S11-S13 2740 B10=BQZ*S13-S14 2750 HECD=PT5*ECD 2760 ECD2=ECD*ECD 2770 HECD2=PT5*ECD2 2780 QECD=QPERP*ECD 2790 HQECD=PT5*QECD 2800 QECD2=QPERP*ECD2 2810 HQECD2=PT5*QECD2 2820 Q2ECD=QPERP2*ECD 2830 Q3ECD=QPERP*Q2ECD 2840 Q2ECD2=QPERP2*ECD2 2850 Q3ECD2=Q2ECD2*QPERP 2860 H0000=X1 2870 H0001=QECD*X2 2880 H0003=-ECD*X4 2890 H0022=HECD*(X1-ECD*X2) 2900 H0011=H0022+Q2ECD2*X3 2910 H0013=-QECD2*X5 2920 H0033=H0022+ECD2*X6 2930 H0100=-QPERP*Z1 2940 H0300=Z4+Y1 2950 H0202=HECD*Z1 2960 H0101=H0202-Q2ECD*Z2 2970 H0103=QECD*Z5 2980 H0301=H0103+QECD*Y2 2990 H0303=H0202-ECD*Z7-ECD*Y4 3000 H0212=HQECD2*Z2 3010 H0223=-HECD2*Z5 3020 H0122=H0212-QPERP*H0202 3030 H0322=H0223+HECD*(H0300-ECD*Y2) 3040 H0113=H0223+Q2ECD2*Z6 3050 H0313=H0212-QECD2*(Z8+Y5) 3060 H0111=H0122+H0212+H0212-Q3ECD2*Z3 3070 H0133=H0122-QECD2*Z8 3080 H0311=H0322+Q2ECD2*(Z6+Y3) 3090 H0333=H0322+H0223+H0223+ECD2*(Z9+Y6) 3100 H1000=-QPERP*W1 3110 H3000=W4+V1 3120 H2002=HECD*W1 3130 H1001=H2002-Q2ECD*W2 3140 H1003=QECD*W5 3150 H3001=H1003+QECD*V2 3160 H3003=H2002-ECD*W7-ECD*V4 3170 H2012=HQECD2*W2 3180 H2023=-HECD2*W5 3190 H1022=H2012-QPERP*H2002 3200 H3022=H2023+HECD*(H3000-ECD*V2) 3210 H1013=H2023+Q2ECD2*W6 3220 H3013=H2012-QECD2*(W8+V5) 3230 H1011=H1022+H2012+H2012-Q3ECD2*W3 3240 H1033=H1022-QECD2*W8 3250 H3011=H3022+Q2ECD2*(W6+V3) 3260 H3033=H3022+H2023+H2023+ECD2*(W9+V6) 3270 H2200=PT5*(S1-T2) 3280 H1100=H2200+QPERP2*T3 3290 H1300=-QPERP*(T6+B1) 3300 H3100=-QPERP*(T6+A1) 3310 H3300=H2200+T9+A4+B4+C1 3320 H2201=HQECD*(S2-T3) 3330 H1101=H2201-QECD*T3+Q3ECD*T4 3340 TEMP =HECD*T6-Q2ECD*T7 3350 H1301=TEMP+HECD*B1-Q2ECD*B2 3360 H3101=TEMP+HECD*A1-Q2ECD*A2 3370 H3301=H2201+QECD*(T10+A5+B5+C2) 3380 H1202=-HQECD*T3 3390 H2102=H1202 3400 H2302=HECD*(T6+B1) 3410 H3202=HECD*(T6+A1) 3420 H2203=HECD*(T6-S6) 3430 H1103=H2203-Q2ECD*T7 3440 TEMP =-HQECD*T3+QECD*T10 3450 H1303=TEMP+QECD*B5 3460 H3103=TEMP+QECD*A5 3470 H3303=H2203+ECD*(T6-T12-A8-B8-C4)+HECD*(A1+B1) 3480 H1212=PT25*ECD2*T3-PT5*Q2ECD2*T4 3490 H2112=H1212 3500 H2312=HQECD2*(T7+B2) 3510 H3212=HQECD2*(T7+A2) 3520 H1223=HQECD2*T7 3530 H2123=H1223 3540 H2323=HECD2*(PT5*T3-T10-B5) 3550 H3223=HECD2*(PT5*T3-T10-A5) 3560 HXXYY=PT25*(ECD*(S1-T2)-ECD2*(S2-T3)) 3570 H2222=HXXYY+HECD2*T3 3580 H1122=HXXYY+PT5*(Q2ECD*T3-Q2ECD2*T4) 3590 TEMP =HQECD*(ECD*T7-T6) 3600 H1322=TEMP+HQECD*(ECD*B2-B1) 3610 H3122=TEMP+HQECD*(ECD*A2-A1) 3620 H3322=HXXYY+HECD*(T9+A4+B4+C1)-HECD2*(T10+A5+B5+C2) 3630 H2211=HXXYY+PT5*Q2ECD2*(S3-T4) 3640 H1111=HXXYY+(HECD2+PT5*Q2ECD)*T3+Q2ECD2*(-THREE*T4+PT5*S3+QPERP2*T 3650 15) 3660 H1311=ONEPT5*QECD2*(T7+B2)-HQECD*(T6+B1)-Q3ECD2*(B3+T8) 3670 H3111=ONEPT5*QECD2*(T7+A2)-HQECD*(T6+A1)-Q3ECD2*(A3+T8) 3680 H3311=HXXYY-HECD2*(QPERP2*T4+T10+A5+B5)+HECD*(T9+A4+B4+C1-ECD*C2)+ 3690 1Q2ECD2*(T11+PT5*S3+A6+B6+C3) 3700 H2213=HQECD2*(T7-S7) 3710 H1113=ONEPT5*QECD2*T7-HQECD2*S7-Q3ECD2*T8 3720 TEMP =HECD2*(PT5*T3-T10)+Q2ECD2*(T11-PT5*T4) 3730 H1313=TEMP-HECD2*B5+Q2ECD2*B6 3740 H3113=TEMP-HECD2*A5+Q2ECD2*A6 3750 H3313=QECD2*(ONEPT5*T7-T13-A9-B9-C5)-HQECD2*(S7-A2-B2) 3760 H2233=HXXYY+HECD2*(S10-T10) 3770 H1133=HXXYY-HECD2*(QPERP2*T4+T10-S10)+PT5*Q2ECD*T3+Q2ECD2*T11 3780 H1333=QECD2*(ONEPT5*T7-T13-B9)-HQECD*(T6+B1)+HQECD2*B2 3790 H3133=QECD2*(ONEPT5*T7-T13-A9)-HQECD*(T6+A1)+HQECD2*A2 3800 H3333=HXXYY+HECD2*(-THREE*(A5+B5)+T3+S10-C2)+ECD2*(-THREE*T10+T14+ 3810 1A10+ B10+C6)+HECD*(T9+A4+B4+C1) 3820 RETURN 3830 END 3840 SUBROUTINE GAMGEN(KOP) IMPLICIT DOUBLE PRECISION(A-H,O-Z) C C -------------------------- C GAUSSIAN 74 (QCPE VERSION) C DECEMBER 1974 C UNIVAC 1108/CDC 7600 C -------------------------- C COMMON/TABLE/C(1200,5) DIMENSION Y(410),F(9) DATA PT15/0.15D0/,PT05/0.05D0/ DATA PT184/0.184D0/ DATA SIX/6.0D0/,FOUR/4.0D0/ DATA IOGAM/3/ C C INITIALIZE THE F(M,T) ROUTINE. CALL FMTSET(0,0,0) C C GENERATE THE DESIRED F(M,T) FOR THE COMPLETE RANGE. C WE WILL COMPUTE F(0,T) THROUGH F(4,T) IN THIS SECTION. T=-PT15 DO 10 I=1,404 T=T+PT05 CALL FMTGEN(F,T,5,ICK) C COPY THE RETURNED VALUES INTO PLACES WHERE THEY CAN BE C REACHED LATER. C(I,2)=F(1) C(I,3)=F(2) C(I,4)=F(3) C(I,5)=F(4) 10 Y(I )=F(5) C C COMPUTE THE INTERPOLATION TABLE WITH THE VALUES AVAILABLE. C THIS IS SOMEWHAT COMPLICATED BY THE FACT THAT SOME OF THE C INTEGRALS ARE IN Y RATHER THAN IN C. DO 40 K=1,5 DO 40 I=1,400 J=I+2 C K INDEXES THE VALUE OF M IN F(M,T). C I INDEXES THE INTERPOLATION TABLE (C). C J INDEXES THE STORED VALUES OF F(M,T). C IF K .EQ. 5, WE MUST USE ALTERNATE CODE BECAUSE THE INTEGRALS C ARE STORED IN Y. IF(K-5)20,30,20 C INTEGRALS IN C, PROCEED AS NORMAL. 20 TEMP1=C(J+1,K+1)+C(J-1,K+1)-(C(J,K+1)+C(J,K+1)) TEMP2=SIX*C(J,K+1)-FOUR*(C(J+1,K+1)+C(J-1,K+1))+C(J-2,K+1)+C(J+2,K $+1) C(I ,K)=C(J,K+1) C(I+400,K)=C(J+1,K+1)-C(J,K+1) C(I+800,K)=(TEMP1-PT184*TEMP2)/SIX GO TO 40 C ALTERNATE CODE USING Y. 30 TEMP1=Y(J+1)+Y(J-1)-(Y(J)+Y(J)) TEMP2=SIX*Y(J)-FOUR*(Y(J+1)+Y(J-1))+Y(J-2)+Y(J+2) C(I ,K)=Y(J) C(I+400,K)=Y(J+1)-Y(J) C(I+800,K)=(TEMP1-PT184*TEMP2)/SIX 40 CONTINUE RETURN END SUBROUTINE FMTSET(KOP1,KOP2,KOP3) IMPLICIT DOUBLE PRECISION(A-H,O-Z) C C ROUTINE TO PRE-SET CUTOFFS FOR FMTGEN. C COMMON/FM/GA(15),RPITWO,FMZERO(15),TOL,CUT0S,CUTSM,CUTML COMMON/IO/IN,IOUT,IPUNCH COMMON/FMCONS/FOUR,ONE,HALF,TWO,ZERO,TEN,TENM9,F20,F42,F500 DIMENSION F(15) DATA FOUR/4.0D0/,ONE/1.0D0/,HALF/0.5D0/,TWO/2.0D0/,ZERO/0.0D0/ DATA TEN/10.0D0/,TENM9/1.0D-9/,F20/20.0D0/,F42/42.0D0/,F500/500.0D $0/ C 2001 FORMAT(21H FROM FMTSET, KOPS = ,3I2,12H AND CUTS = ,3E10.3) C C COMPUTE PI RELATED CONSTANTS. PI=FOUR*DATAN(ONE) GA(1)=DSQRT(PI) RPITWO=HALF*GA(1) C C FILL GAMMA FUNCTION ARRAY. TOL=HALF DO 10 I=2,15 GA(I)=GA(I-1)*TOL 10 TOL=TOL+ONE C C FILL FMZERO (FOR ARGUMENT OF ZERO). TOL=ONE FMZERO(1)=ONE DO 20 I=2,15 TOL=TOL+TWO 20 FMZERO(I)=ONE/TOL C C OBTAIN CUTOFFS. C C CUT0S = CUTOFF, ZERO TO SMALL. C CUTSM = CUTOFF, SMALL TO MEDIUM. C CUTML = CUTOFF, MEDIUM TO LARGE. C C CUTOS IS OBTAINED FROM KOP1. CUT0S=ZERO IF(KOP1)30,40,30 30 CUT0S=TEN**(-2*KOP1) C C CUTSM IS OBTAINED FROM KOP2. C KOP2 SPECIFIES THE DESIRED ACCURACY OF THE ASSYMPTOTIC C EXPANSION. ROUTINE FMTGEN IS CALLED TO CHECK THE VIABILITY C OF THE ASSYMPTOTIC EXPANSION FOR VARIOUS VALUES OF CUTSM. 40 TOL=TENM9 CUTSM=TEN IF(KOP2)50,90,50 50 TOL=TEN**(-6-KOP2) T=F20 60 CALL FMTGEN(F,T,1,ICK) IF(ICK)80,70,80 70 T=T-ONE IF(T-TEN)90,60,60 80 CUT=T+ONE C C CUTML IS TAKEN FROM KOP3. 90 CUTML=F42 IF(KOP3)100,130,100 100 IF(KOP3-7)120,110,120 110 CUTML=F500 GO TO 130 120 CUTML=DFLOAT(KOP3)*FOUR 130 IF((KOP1+KOP2+KOP3).NE.0)WRITE(IOUT,2001)KOP1,KOP2,KOP3,CUT0S, $ CUTSM,CUTML RETURN END SUBROUTINE FMTGEN(F,T,M,ICK) IMPLICIT DOUBLE PRECISION(A-H,O-Z) C C -------------------------- C GAUSSIAN 74 (QCPE VERSION) C DECEMBER 1974 C UNIVAC 1108/CDC 7600 C -------------------------- C DIMENSION F(M) COMMON/IO/IN,IOUT,IPUNCH COMMON/FMCONS/FOUR,ONE,HALF,TWO,ZERO,TEN,TENM9,F20,F42,F500 COMMON/FM/GA(15),RPITWO,FMZERO(15),TOL,CUT0S,CUTSM,CUTML EQUIVALENCE(APPROX,OLDSUM) DATA ZERO/0.0D0/,TEN/10.0D0/,HALF/0.5D0/,ONE/1.0D0/,F42/42.0D0/ DATA TWO/2.0D0/ 2001 FORMAT(41H0FAILURE IN FMGEN FOR SMALL T: IX ' 50, / $ 6H IX = ,I3,7H, T = ,E20.14) 2002 FORMAT(37H0FAILURE IN FMGEN FOR INTERMEDIATE T:/ $ 6H T = ,E20.14) C C ICK IS AN ERROR INDICATOR. C ON RETURN, ICK=0 IMPLIES THAT ALL IS WELL. C IF ON RETURN, ICK IS NON-ZERO, THE ASSYMPTOTIC EXPANSION C HAS FAILED. ICK=0 C TEST FOR TYPE OF ALGORITHM. IF(DABS(T)-CUT0S)1,1,3 C*********************************************************************** C FILL F(M) FOR ARGUMENT OF ZERO. C************************************************************************ 1 DO 2 I=1,M 2 F(I)=FMZERO(I) RETURN C TEST FOR EVALUATION OF THE EXP. 3 TEXP=ZERO IF(DABS(T)-CUTML)4,150,150 4 TEXP=DEXP(-T) IF(DABS(T)-CUTSM)10,80,80 C*********************************************************************** C 0 .LT. T .LT. 10 C*********************************************************************** 10 A=DFLOAT(M-1)+HALF TERM=ONE/A SUM=TERM DO 20 IX=2,200 A=A+ONE TERM=TERM*T/A SUM=SUM+TERM IF(DABS(TERM/SUM)-TOL)30,20,20 20 CONTINUE WRITE(IOUT,2001)IX,T STOP 30 F(M)=HALF*SUM*TEXP GO TO 160 C C*********************************************************************** C 10 .LE. T .LT. 42 C*********************************************************************** 80 A=DFLOAT(M-1) B=A+HALF A=A-HALF TX=ONE/T MM1=M-1 APPROX=RPITWO*DSQRT(TX)*(TX**MM1) IF(MM1)90,110,90 90 DO 100 IX=1,MM1 B=B-ONE 100 APPROX=APPROX*B 110 FIMULT=HALF*TEXP*TX SUM=ZERO IF(FIMULT)120,140,120 120 FIPROP=FIMULT/APPROX TERM=ONE SUM =ONE NOTRMS=IDINT(T)+MM1 DO 130 IX=2,NOTRMS TERM=TERM*A*TX SUM=SUM+TERM IF(DABS(TERM*FIPROP/SUM)-TOL)140,140,130 130 A=A-ONE WRITE(IOUT,2002)T ICK=1 RETURN 140 F(M)=APPROX-FIMULT*SUM GO TO 160 C*********************************************************************** C T .GE. 42 C*********************************************************************** 150 TX=DFLOAT(M)-HALF F(M)=HALF*GA(M)/(T**TX) C*********************************************************************** C RECUR DOWNWARDS TO F(1) C*********************************************************************** 160 TX=T+T SUM=DFLOAT(M+M-3) MM1=M-1 IF(MM1)170,190,170 170 DO 180 IX=1,MM1 F(M-IX)=(TX*F(M-IX+1)+TEXP)/SUM 180 SUM=SUM-TWO 190 RETURN END IF(MM1)170,190,170 170 DO 180 IX=1,MM1 F(M-IX)=(TX*F(M-IX+1)+TEXP)/SUM 180 SUM=SUM-TWO 190 RETURN END