PROGRAM TESTEX C C MULTST Test Package with example test program and early ADAPT C EXTERNAL ADAPT CHARACTER *6 SBNAME REAL DIF(6), EXP(6), EPS INTEGER TST(6), NDIMS(5), I, MAXCLS DATA SBNAME/' ADAPT'/ DATA (DIF(I),I=1,6)/110.0,600.0,600.0,100.0,150.0,100.0/ DATA (EXP(I),I=1,6)/1.5,2.0,2.0,1.0,2.0,2.0/ DATA (TST(I),I=1,6)/1,2,3,4,5,6/ DATA (NDIMS(I),I=1,5)/2,3,4,6,8/ EPS = 1E-6 MAXCLS = 20000 CALL MULTST(20,6,TST,6,DIF,EXP,5,NDIMS,SBNAME,ADAPT,EPS,10000) END C SUBROUTINE MULTST(NSAMP, TSTLIM, TSTFNS, TSTMAX, DIFCLT, EXPNTS, * NDIML, NDIMS, SBNAME, SUBRTN, EPS, MAXCLS) C C SUBROUTINE FOR TESTING MULTIDIMENSIONAL INTEGRATION ROUTINES C C AS DESCRIBED IN C Genz, A. (1987), A Package for Testing Multiple Integration C Subroutines, in {\it Numerical Integration}, P. Keast and C G. Fairweather (Eds.), D. Riedel, pp. 337--340. C BY ALAN GENZ, DEPARTMENT of MATHEMATICS, WASHINGTON STATE C UNIVERSITY, PULLMAN, WASHINGTON 99164-3113, U.S.A. C alangenz@wsu.edu C C************** PARAMETERS ********************************************* C NSAMP INTEGER NUMBER OF SAMPLES, MUST NOT EXCEED 50. C TSTLIM INTEGER NUMBER OF TEST INTEGRANDS TO BE USED. C TSTFNS INTEGER ARRAY OF DIMENSION(TSTLIM) OF TEST INTEGRAND INDICES. C TSTFNS(I) SPECIFIES THE INDEX OF THE ITH INTEGRAND TO BE C TESTED. TSTFNS(I) MUST NOT EXCEED 6. C TSTMAX INTEGER C DIFCLT REAL ARRAY OF DIMENSION(TSTMAX) OF DIFFICULTY LEVELS. C DIFCLT(I) SPECIFIES THE DIFFICULTY LEVEL FOR THE INTEGRAND C WITH INDEX I. C EXPNTS REAL ARRAY OF DIMENSION(TSTMAX) OF DIFFICULTY EXPONENTS. C EXPNTS(I) SPECIFIES THE DIFFICULTY EXPONENT FOR THE INTEGRAND C WITH INDEX I. C NDIML INTEGER C NDIMS INTEGER ARRAY OF DIMENSION(NDIML). NDIMS(J) SPECIFIES THE C NUMBER OF VARIABLES FOR THE INTEGRALS IN THE JTH SERIES OF C TESTS. NDIMS(J) MUST NOT EXCEED 20. C SBNAME CHARACTER STRING OF LENGTH SIX USED TO STORE THE C NAME OF THE INTEGRATION SUBROUTINE TO BE TESTED. C SUBRTN EXTERNALLY DECLARED INTEGRATION SUBROUTINE TO BE TESTED. C THE CALLING SEQUENCE SHOULD BE THE SAME AS NAG ROUTINE C D01FCF(MARK10). C EPS REAL REQUIRED RELATIVE ACCURACY PARAMETER FOR ALL TESTS. C MAXCLS INTEGER MAXIMUM NUMBER OF INTEGRAND CALLS FOR ALL TESTS. C*********************************************************************** C EXTERNAL FUNCTN, SUBRTN INTEGER I, ITST, IFAIL, ITEST, J, K, LENWRK, MAXCLS, N, * NCONF, NDIML, NSAMP, NDIM, NDIMV, MXTSFN, MAXDIM, * RCALSA, RCALSB, TSTLIM, TSTFNS(TSTLIM), TSTMAX, MAXSMP, IT, * DIGITS, IFAILS, MINCLS PARAMETER (MXTSFN = 6, MAXDIM = 20, MAXSMP = 50, LENWRK = 50000) INTEGER IDFCLT(MXTSFN), NDIMS(NDIML) REAL MEDACT(MAXSMP), MEDEST(MAXSMP), MEDDSC(MAXSMP), * MEDCLA(MXTSFN), MEDCLB(MXTSFN), MEDCLS(MAXSMP), CONCOF, * CALLSA(MXTSFN,MXTSFN), CALLSB(MXTSFN,MXTSFN), DFCLT, EXN, EPS, * ERREST, EXPONS(MXTSFN), ERRLOG, ESTLOG, MEDREL, MEDRLL(MAXSMP), * ERSREL(MXTSFN,MXTSFN), ERSDSC(MXTSFN,MXTSFN), * ERSACT(MXTSFN,MXTSFN), ERSEST(MXTSFN,MXTSFN), ONE, * FINEST, RANDOM, RELERR, SMALL, TRELIB(MXTSFN), * TACTRS(MXTSFN), TESTRS(MXTSFN), VALUE, TERDSC(MXTSFN), * WRKSTR(LENWRK), ZERO, TQUALT(MXTSFN), QUALTY(MXTSFN,MXTSFN), * TCALSA(MXTSFN), TCALSB(MXTSFN), QALITY, * QALLTY(MAXSMP), TACTRB(MXTSFN), TESTRB(MXTSFN), TERDSB(MXTSFN), * ERSESB(MXTSFN,MXTSFN), ERSACB(MXTSFN,MXTSFN), * ERSDSB(MXTSFN,MXTSFN), MEDACB(MXTSFN), MEDESB(MXTSFN), * MEDDSB(MXTSFN), DIFCLT(TSTMAX), EXPNTS(TSTMAX), VALINT REAL A(MAXDIM), B(MAXDIM), ALPHA(MAXDIM), BETA(MAXDIM), FUNCTN CHARACTER*6 SBNAME CHARACTER*14 TYPE(MXTSFN) COMMON /PRMBLK/ NDIM, A, B COMMON /TSTBLK/ ITST, ALPHA, BETA, DFCLT, EXN, VALUE SAVE TYPE DATA (TYPE(I), I = 1,MXTSFN)/ 'OSCILLATORY', 'PRODUCT PEAK', * 'CORNER PEAK', 'GAUSSIAN', * 'C0 FUNCTION', 'DISCONTINUOUS'/ C C INITIALISE AND COMPUTE CONFIDENCE COEFFICIENT C ZERO = 0 ONE = 1 CONCOF = 0 NCONF = MAX( 1, (2*NSAMP)/5-2 ) DO 10 I = 1,NCONF CONCOF = 1 + (NSAMP-NCONF+I)*CONCOF/(NCONF-I+1) 10 CONTINUE CONCOF = 1 - CONCOF/2**FLOAT(NSAMP-1) SMALL = 1 20 SMALL = SMALL/2 IF ( 1 + SMALL .NE. 1 ) GO TO 20 SMALL = 2*SMALL DO 30 IT = 1,TSTLIM ITEST = TSTFNS(IT) IDFCLT(IT) = DIFCLT(ITEST) EXPONS(IT) = EXPNTS(ITEST) 30 CONTINUE C C BEGIN MAIN LOOP FOR DIFFERENT NUMBERS OF VARIABLES C DO 140 NDIMV = 1,NDIML NDIM = NDIMS(NDIMV) IF ( MOD(NDIMV-1,6) .EQ. 0 ) THEN WRITE (*,99999) NSAMP, (IDFCLT(J),J=1,TSTLIM) 99999 FORMAT (9X, 'TEST RESULTS WITH', I4,' SAMPLES PER TEST'/ * ' DIFFICULTY LEVELS', 10I6) WRITE (*,99998) (EXPONS(J),J=1,TSTLIM) 99998 FORMAT (' EXPONENTS ', 10F6.1) DIGITS = -LOG10(EPS) WRITE (*,99997) DIGITS, MAXCLS 99997 FORMAT (' REQUESTED DIGITS =',I3, ', MAXIMUM VALUES =',I8) WRITE (*,99996) SBNAME, CONCOF 99996 FORMAT (/9X, A6, ' TESTS, VARIABLE RESULTS WITH CONFIDENCE', * F5.2/' VARI- INTEGRAND CORRECT DIGITS RELIA-', * ' WRONG INTEGRAND QUALITY TOTAL'/' ABLES TYPE', * 6X, 'ESTIMATED ACTUAL BILITY DIGITS VALUES ',12X, * 'FAILS') ENDIF C C BEGIN LOOP FOR DIFFERENT TEST INTEGRANDS C DO 130 IT = 1,TSTLIM ITEST = TSTFNS(IT) ITST = ITEST EXN = EXPNTS(ITEST) DFCLT = DIFCLT(ITEST) DO 40 J = 1,NDIM A(J) = 0 B(J) = 1 40 CONTINUE IFAILS = 0 MEDREL = 0 C C BEGIN LOOP FOR DIFFERENT SAMPLES C DO 120 K = 1,NSAMP IFAIL = 1 DO 50 N = 1,NDIM ALPHA(N) = RANDOM() BETA(N) = RANDOM() 50 CONTINUE VALUE = VALINT(ITEST) C C CALL INTEGRATION SUBROUTINE C MINCLS = 4*2**NDIM CALL SUBRTN(NDIM, A, B, MINCLS, MAXCLS, FUNCTN, EPS, * ERREST, LENWRK, WRKSTR, FINEST, IFAIL) RELERR = ABS( (FINEST-VALUE)/VALUE ) IFAILS = IFAILS + MIN( IFAIL, 1 ) RELERR = MAX( MIN(ONE, RELERR), SMALL ) ERRLOG = MAX( ZERO, -LOG10(RELERR) ) ERREST = MAX( MIN(ONE, ERREST), SMALL) ESTLOG = MAX( ZERO, -LOG10(ERREST) ) MEDDSC(K) = MAX( ZERO, ESTLOG-ERRLOG ) MEDEST(K) = ESTLOG MEDACT(K) = ERRLOG MEDCLS(K) = MINCLS IF ( ERREST .GE. RELERR ) MEDREL = MEDREL + 1 120 CONTINUE C C END LOOP FOR DIFFERENT SAMPLES AND COMPUTE MEDIANS C CALL MEDIAN(NSAMP, MEDEST) CALL MEDIAN(NSAMP, MEDACT) CALL MEDIAN(NSAMP, MEDCLS) CALL MEDIAN(NSAMP, MEDDSC) MEDREL = MEDREL/NSAMP TRELIB(IT) = MEDREL TACTRS(IT) = MEDACT(2) TESTRS(IT) = MEDEST(2) TERDSC(IT) = MEDDSC(2) TCALSA(IT) = MEDCLS(2) TCALSB(IT) = MEDCLS(3) TACTRB(IT) = MEDACT(3) TESTRB(IT) = MEDEST(3) TERDSB(IT) = MEDDSC(3) ERSREL(ITEST,NDIMV) = MEDREL ERSEST(ITEST,NDIMV) = MEDEST(2) ERSACT(ITEST,NDIMV) = MEDACT(2) ERSDSC(ITEST,NDIMV) = MEDDSC(2) ERSESB(ITEST,NDIMV) = MEDEST(3) ERSACB(ITEST,NDIMV) = MEDACT(3) ERSDSB(ITEST,NDIMV) = MEDDSC(3) CALLSA(ITEST,NDIMV) = MEDCLS(2) CALLSB(ITEST,NDIMV) = MEDCLS(3) QALITY = 0 IF ( MEDCLS(1) .NE. ZERO ) QALITY = (MEDACT(1)+1)* * ( MEDEST(1) + 1 - MEDDSC(1) )/LOG(MEDCLS(1)) TQUALT(IT) = QALITY QUALTY(ITEST,NDIMV) = QALITY RCALSA = MEDCLS(2) RCALSB = MEDCLS(3) WRITE (*,99995) NDIM, TYPE(ITEST), MEDEST(2), * MEDEST(3), MEDACT(2), MEDACT(3), MEDREL, MEDDSC(2), * MEDDSC(3), RCALSA, RCALSB, QALITY, IFAILS 99995 FORMAT (2X, I2, 2X, A14, 2(F4.1, ',', F4.1, 1X), F4.2, * F4.1, ',', F3.1, I7, ',', I7, F6.2, I5) 130 CONTINUE C C END LOOP FOR DIFFERENT TEST INTEGRANDS C CALL MEDIAN(TSTLIM, TACTRS) CALL MEDIAN(TSTLIM, TRELIB) CALL MEDIAN(TSTLIM, TESTRS) CALL MEDIAN(TSTLIM, TERDSC) CALL MEDIAN(TSTLIM, TACTRB) CALL MEDIAN(TSTLIM, TESTRB) CALL MEDIAN(TSTLIM, TERDSB) CALL MEDIAN(TSTLIM, TQUALT) CALL MEDIAN(TSTLIM, TCALSA) CALL MEDIAN(TSTLIM, TCALSB) RCALSA = TCALSA(1) RCALSB = TCALSB(1) WRITE (*,99994) NDIM, TESTRS(1), TESTRB(1), TACTRS(1), * TACTRB( 1), TRELIB(1), TERDSC(1), TERDSB(1), * RCALSA, RCALSB, TQUALT(1) 99994 FORMAT (2X, I2, ' MEDIANS', 6X, 2(F4.1, ',', F4.1, 1X), * F4.2, F4.1, ',', F3.1, I7, ',', I7, F6.2/1X) 140 CONTINUE C C END LOOP FOR DIFFERENT NUMBERS OF VARIABLES C IF (NDIML.GT.1) THEN WRITE (*,99993) SBNAME, (NDIMS(NDIMV),NDIMV=1,NDIML) 99993 FORMAT (/6X, A6, ' TEST INTEGRAND MEDIANS FOR VARIABLES', 12I3) WRITE (*,99992) 99992 FORMAT ( * 8X, 'INTEGRAND CORRECT DIGITS RELIA- WRONG', * ' INTEGRAND QUALITY'/8X, ' TYPE ESTIMATED', * ' ACTUAL BILITY DIGITS VALUES') DO 160 IT = 1,TSTLIM ITEST = TSTFNS(IT) DO 150 NDIMV = 1,NDIML MEDACT(NDIMV) = ERSACT(ITEST,NDIMV) MEDEST(NDIMV) = ERSEST(ITEST,NDIMV) MEDDSC(NDIMV) = ERSDSC(ITEST,NDIMV) MEDACB(NDIMV) = ERSACB(ITEST,NDIMV) MEDESB(NDIMV) = ERSESB(ITEST,NDIMV) MEDDSB(NDIMV) = ERSDSB(ITEST,NDIMV) MEDRLL(NDIMV) = ERSREL(ITEST,NDIMV) QALLTY(NDIMV) = QUALTY(ITEST,NDIMV) MEDCLA(NDIMV) = CALLSA(ITEST,NDIMV) MEDCLB(NDIMV) = CALLSB(ITEST,NDIMV) 150 CONTINUE CALL MEDIAN(NDIML, MEDRLL) CALL MEDIAN(NDIML, MEDACT) CALL MEDIAN(NDIML, MEDEST) CALL MEDIAN(NDIML, MEDDSC) CALL MEDIAN(NDIML, MEDACB) CALL MEDIAN(NDIML, MEDESB) CALL MEDIAN(NDIML, MEDDSB) CALL MEDIAN(NDIML, QALLTY) CALL MEDIAN(NDIML, MEDCLA) CALL MEDIAN(NDIML, MEDCLB) RCALSA = MEDCLA(1) RCALSB = MEDCLB(1) WRITE (*,99991) TYPE(ITEST), MEDEST(1), MEDESB(1), * MEDACT(1), MEDACB(1),MEDRLL(1), MEDDSC(1), MEDDSB(1), * RCALSA, RCALSB, QALLTY(1) 99991 FORMAT (6X, A14, 2(F4.1, ',', F4.1, 1X), F4.2, F4.1, ',', * F3.1, I7, ',', I7, F6.2) TACTRS(IT) = MEDACT(1) TESTRS(IT) = MEDEST(1) TERDSC(IT) = MEDDSC(1) TACTRB(IT) = MEDACB(1) TESTRB(IT) = MEDESB(1) TERDSB(IT) = MEDDSB(1) TCALSA(IT) = MEDCLA(1) TCALSB(IT) = MEDCLB(1) TRELIB(IT) = MEDRLL(1) TQUALT(IT) = QALLTY(1) 160 CONTINUE CALL MEDIAN(TSTLIM, TACTRS) CALL MEDIAN(TSTLIM, TESTRS) CALL MEDIAN(TSTLIM, TERDSC) CALL MEDIAN(TSTLIM, TACTRB) CALL MEDIAN(TSTLIM, TESTRB) CALL MEDIAN(TSTLIM, TERDSB) CALL MEDIAN(TSTLIM, TRELIB) CALL MEDIAN(TSTLIM, TQUALT) CALL MEDIAN(TSTLIM, TCALSA) CALL MEDIAN(TSTLIM, TCALSB) RCALSA = TCALSA(1) RCALSB = TCALSB(1) WRITE (*,99990) TESTRS(1), TESTRB(1), TACTRS(1), TACTRB(1), * TRELIB(1), TERDSC(1), TERDSB(1), RCALSA, RCALSB, TQUALT(1) 99990 FORMAT (5X, 'GLOBAL MEDIANS ', 2(F4.1, ',', F4.1, 1X), F4.2, * F4.1, ',', F3.1, I7, ',', I7, F6.2/1X) ENDIF END SUBROUTINE MEDIAN(N, R) C C SUBROUTINE TO COMPUTE MEDIANS FOR ARRAY R OF REALS. C ON EXIT R(1) CONTAINS THE MEDIANS, (R(2),R(3)) SPECIFIES C THE CONFIDENCE INTERVAL C INTEGER J, K, KMAX, N, NCONF, ND REAL R(*), RMAX, RMED DO 20 J = 1,N KMAX = J DO 10 K = J , N IF ( R(K) .GT. R(KMAX) ) KMAX = K 10 CONTINUE RMAX = R(KMAX) R(KMAX) = R(J) R(J) = RMAX 20 CONTINUE ND = N/2 IF ( MOD(N, 2) .EQ. 0 ) THEN RMED = (R(ND)+R(ND+1))/2 ELSE RMED = R(ND+1) ENDIF NCONF = MAX( 1, (2*N)/5-2 ) RMAX = R(N-NCONF+1) R(3) = R(NCONF) R(2) = RMAX R(1) = RMED END REAL FUNCTION VALINT(ITEST) C C FUNCTION TO COMPUTE CORRECT VALUES FOR INTEGRALS C INTEGER ITEST, ISUM, MAXDIM, NDIM, ITST, J PARAMETER ( MAXDIM = 20 ) REAL A(MAXDIM), B(MAXDIM), ALPHA(MAXDIM), BETA(MAXDIM) REAL TWOPI, FNDIM, SUM, DFCLT, EXN, VALUE, * SIGN, DFACT, SGNDM, AB, PI, PHI INTEGER IC(MAXDIM) PARAMETER ( PI = 3.14159 26535 89793 23844, TWOPI = 2*PI ) COMMON /PRMBLK/ NDIM, A, B COMMON /TSTBLK/ ITST, ALPHA, BETA, DFCLT, EXN, VALUE FNDIM = NDIM SUM = 0 DO 10 J=1,NDIM SUM = SUM + ALPHA(J) 10 CONTINUE DFACT = SUM*FNDIM**EXN/DFCLT DO 20 J=1,NDIM ALPHA(J) = ALPHA(J)/DFACT 20 CONTINUE IF ( ITEST .EQ. 1 ) THEN DO 40 J=1,NDIM B(J) = ALPHA(J) IC(J) = 0 40 CONTINUE VALUE = 0 50 ISUM = 0 SUM = TWOPI*BETA(1) DO 70 J = 1,NDIM IF ( IC(J) .NE. 1 ) SUM = SUM + ALPHA(J) ISUM = ISUM + IC(J) 70 CONTINUE SIGN = 1 + 2*((ISUM/2)*2-ISUM) IF ( (NDIM/2)*2 .EQ. NDIM ) VALUE = VALUE + SIGN*COS(SUM) IF ( (NDIM/2)*2 .NE. NDIM ) VALUE = VALUE + SIGN*SIN(SUM) DO 80 J = 1,NDIM IC(J) = IC(J) + 1 IF ( IC(J) .LT. 2 ) GO TO 50 IC(J) = 0 80 CONTINUE VALUE = VALUE*FLOAT((-1)**(NDIM/2)) ELSE IF ( ITEST .EQ. 2 ) THEN VALUE = 1 DO 100 J=1,NDIM VALUE = VALUE*(ATAN((1-BETA(J))*ALPHA(J))-ATAN(-BETA(J)* * ALPHA(J)))*ALPHA(J) 100 CONTINUE ELSE IF ( ITEST .EQ. 3 ) THEN VALUE = 0 SGNDM = 1 DO 120 J=1,NDIM SGNDM = -SGNDM/FLOAT(J) B(J) = ALPHA(J) IC(J) = 0 120 CONTINUE 130 SUM = 1 ISUM = 0 DO 150 J=1,NDIM IF ( IC(J) .NE. 1 ) SUM = SUM + ALPHA(J) ISUM = ISUM + IC(J) 150 CONTINUE SIGN = 1 + 2*((ISUM/2)*2-ISUM) VALUE = VALUE + SIGN/SUM DO 160 J=1,NDIM IC(J) = IC(J) + 1 IF ( IC(J) .LT. 2 ) GO TO 130 IC(J) = 0 160 CONTINUE VALUE = VALUE*SGNDM ELSE IF ( ITEST .EQ. 4 ) THEN VALUE = 1 AB = SQRT(2D0) DO 180 J = 1, NDIM VALUE = VALUE*(SQRT(PI)/ALPHA(J)) * * ( PHI( (1-BETA(J))*AB*ALPHA(J) ) * - PHI( -BETA(J) *AB*ALPHA(J) ) ) 180 CONTINUE ELSE IF ( ITEST .EQ. 5 ) THEN VALUE = 1 DO 200 J=1,NDIM AB = ALPHA(J)*BETA(J) VALUE = VALUE*(2-EXP(-AB)-EXP(AB-ALPHA(J)))/ALPHA(J) 200 CONTINUE ELSE IF ( ITEST .EQ. 6 ) THEN VALUE = 1 DO 220 J=1,NDIM IF ( J.GT.2 ) BETA(J) = 1 VALUE = VALUE*(EXP(ALPHA(J)*BETA(J))-1)/ALPHA(J) 220 CONTINUE ENDIF VALINT = VALUE END REAL FUNCTION FUNCTN(NDIM, Z) C C TEST INTEGRAND FAMILY FUNCTION C INTEGER NDIM, MAXDIM, ITST, J PARAMETER ( MAXDIM = 20 ) REAL SUM, DFCLT, EXN, Z(NDIM), TWOPI, EXPMAX, * ALPHA(MAXDIM), BETA(MAXDIM), VALUE, F COMMON /TSTBLK/ ITST, ALPHA, BETA, DFCLT, EXN, VALUE PARAMETER ( TWOPI = 6.28318 53072, EXPMAX = 100 ) F = 0 IF ( ITST .EQ. 1 ) THEN SUM = TWOPI*BETA(1) DO 20 J = 1,NDIM SUM = SUM + Z(J) 20 CONTINUE F = COS(SUM) ELSE IF ( ITST .EQ. 2 ) THEN SUM = 1 DO 40 J = 1,NDIM SUM = SUM/( 1/ALPHA(J)**2 + (Z(J)-BETA(J))**2 ) 40 CONTINUE F = SUM ELSE IF ( ITST .EQ. 3 ) THEN C C Note: for this case, the BETA's are used to randomly select C (from 2**NDIM possibilities) a corner for the peak C SUM = 1 DO 60 J = 1,NDIM IF ( BETA(J) .LT. 0.5 ) THEN SUM = SUM + Z(J) ELSE SUM = SUM + ALPHA(J) - Z(J) ENDIF 60 CONTINUE F = SUM**(-NDIM-1) ELSE IF ( ITST .EQ. 4 ) THEN SUM = 0 DO 80 J = 1,NDIM SUM = SUM + (ALPHA(J)*(Z(J)-BETA(J)))**2 80 CONTINUE IF ( SUM .LT. EXPMAX ) F = EXP(-SUM) ELSE IF ( ITST .EQ. 5 ) THEN SUM = 0 DO 100 J = 1,NDIM SUM = SUM + ALPHA(J)*ABS(Z(J)-BETA(J)) 100 CONTINUE IF ( SUM .LT. EXPMAX ) F = EXP(-SUM) ELSE IF ( ITST .EQ. 6 ) THEN DO 120 J = 1,NDIM IF ( Z(J) .GT. BETA(J) ) GO TO 140 120 CONTINUE SUM = 0 DO 130 J = 1,NDIM SUM = SUM + ALPHA(J)*Z(J) 130 CONTINUE F = EXP(SUM) ENDIF 140 FUNCTN = F END REAL FUNCTION RANDOM() C C PORTABLE RANDOM NUMBER GENERATOR C INTEGER SEED, A, P, B15, B16, XHI, XALO, LEFTLO, FHI, K REAL PRIME PARAMETER (A = 16807, B15 = 32768, B16 = 65536, P = 2147483647) PARAMETER (PRIME = P) SAVE SEED DATA SEED/123456/ XHI = SEED/B16 XALO = (SEED-XHI*B16)*A LEFTLO = XALO/B16 FHI = XHI*A + LEFTLO K = FHI/B15 SEED = (((XALO-LEFTLO*B16)-P)+(FHI-K*B15)*B16) + K IF (SEED.LT.0) SEED = SEED + P RANDOM = SEED/PRIME END REAL FUNCTION PHI(Z) C C Normal distribution probabilities accurate to 1.e-7. C Z = no. of standard deviations from the mean. C PHI = probability to the left of Z. C EXPNTL = the probability density. C C Based upon algorithm 5666 for the error function, from: C Hart, J.F. et al, 'Computer Approximations', Wiley 1968 C C Programmer: Alan Miller C C Latest revision - 30 March 1986 C REAL P0, P1, P2, P3, P4, P5, P6, * Q0, Q1, Q2, Q3, Q4, Q5, Q6, Q7, * Z, P, EXPNTL, ROOTPI, ZABS PARAMETER(P0 = 220.20 68679 12376 1, * P1 = 221.21 35961 69931 1, * P2 = 112.07 92914 97870 9, * P3 = 33.912 86607 83830 0, * P4 = 6.3739 62203 53165 0, * P5 = .70038 30644 43688 1, * P6 = .035262 49659 98910 9) PARAMETER(Q0 = 440.41 37358 24752 2, * Q1 = 793.82 65125 19948 4, * Q2 = 637.33 36333 78831 1, * Q3 = 296.56 42487 79673 7, * Q4 = 86.780 73220 29460 8, * Q5 = 16.064 17757 92069 5, * Q6 = 1.7556 67163 18264 2, * Q7 = .088388 34764 83184 4) PARAMETER(ROOTPI = 2.5066 28274 63100 1) C ZABS = ABS(Z) C C |Z| > 12 C IF ( ZABS .GT. 12 ) THEN P = 0 ELSE C C |Z| <= 12 C EXPNTL = EXP(-ZABS**2/2) C C |Z| < 7 C IF ( ZABS .LT. 7 ) THEN P = EXPNTL*((((((P6*ZABS + P5)*ZABS + P4)*ZABS + P3)*ZABS * + P2)*ZABS + P1)*ZABS + P0)/(((((((Q7*ZABS + Q6)*ZABS * + Q5)*ZABS + Q4)*ZABS + Q3)*ZABS + Q2)*ZABS + Q1)*ZABS * + Q0) C C |Z| >= CUTOFF. C ELSE P = EXPNTL/(ZABS + 1/(ZABS + 2/(ZABS + 3/(ZABS + 4/ * (ZABS + 0.65)))))/ROOTPI END IF END IF IF (Z .GT. 0) P = 1 - P PHI = P END C SUBROUTINE ADAPT(NDIM,A,B,MINPTS,MAXPTS,FUNCTN,EPS,RELERR,LENWRK, * WRKSTR,FINEST,IFAIL) C ADAPTIVE MULTIDIMENSIONAL INTEGRATION SUBROUTINE C C AUTHOR: ALAN GENZ C MATHEMATICS DEPARTMENT C WASHINGTON STATE UNIVERSITY C PULMAN, WA 99164-3130 C C************** PARAMETERS FOR ADAPT ******************************** C***** INPUT PARAMETERS C NDIM NUMBER OF VARIABLES, MUST EXCEED 1, BUT NOT EXCEED 100 C A REAL ARRAY OF LOWER LIMITS, WITH DIMENSION NDIM C B REAL ARRAY OF UPPER LIMITS, WITH DIMENSION NDIM C MINPTS MINIMUM NUMBER OF FUNCTION EVALUATIONS TO BE ALLOWED, C MINPTS MUST NOT EXCEED MAXPTS. IF MINPTS<0 THEN THE C ROUTINE ASSUMES A PREVIOUS CALL HAS BEEN MADE WITH THE C SAME INTEGRAND AND CONTINUES THAT CALCULATION. C MAXPTS MAXIMUM NUMBER OF FUNCTION EVALUATIONS TO BE ALLOWED, C WHICH MUST BE AT LEAST RULCLS, WHERE C RULCLS = 2**NDIM+2*NDIM**2+2*NDIM+1, WHEN NDIM <16 AND C RULCLS = (NDIM*(14-NDIM*(6-4*NDIM))/3+1, WHEN NDIM >15. C FOR NDIM = 2 3 4 5 6 7 8 9 10 11 12 C RULCLS = 17 33 57 93 149 241 401 693 1245 2313 4409 C A SUGGESTED STARTING VALUE FOR MAXPTS IS 100*RULCLS. IF C THIS IS NOT LARGE ENOUGH FOR THE REQUIRED ACCURACY, THEN C MAXPTS (AND LENWRK) SHOULD BE INCREASED ACCORDINGLY. C FUNCTN EXTERNALLY DECLARED USER DEFINED FUNCTION TO BE INTEGRATED. C IT MUST HAVE PARAMETERS (NDIM,Z), WHERE Z IS A REAL ARRAY C OF DIMENSION NDIM. C EPS REQUIRED RELATIVE ACCURACY C LENWRK LENGTH OF ARRAY WRKSTR OF WORKING STORAGE, THE ROUTINE C NEEDS (2*NDIM+3)*(1+MAXPTS/RULCLS)/2 FOR LENWRK IF C MAXPTS FUNCTION CALLS ARE USED. C***** OUTPUT PARAMETERS C MINPTS ACTUAL NUMBER OF FUNCTION EVALUATIONS USED BY ADAPT C WRKSTR REAL ARRAY OF WORKING STORAGE OF DIMENSION (LENWRK). C RELERR ESTIMATED RELATIVE ACCURACY OF FINEST C FINEST ESTIMATED VALUE OF INTEGRAL C IFAIL IFAIL = 0 FOR NORMAL EXIT, WHEN ESTIMATED RELATIVE ACCURACY C RELERR IS LESS THAN EPS WITH MAXPTS OR LESS FUNCTION C CALLS MADE. C IFAIL = 1 IF MAXPTS WAS TOO SMALL FOR ADAPT TO OBTAIN THE C REQUIRED RELATIVE ACCURACY EPS. IN THIS CASE ADAPT C RETURNS A VALUE OF FINEST WITH ESTIMATED RELATIVE C ACCURACY RELERR. C IFAIL = 2 IF LENWRK TOO SMALL FOR MAXPTS FUNCTION CALLS. IN C THIS CASE ADAPT RETURNS A VALUE OF FINEST WITH C ESTIMATED ACCURACY RELERR USING THE WORKING STORAGE C AVAILABLE, BUT RELERR WILL BE GREATER THAN EPS. C IFAIL = 3 IF NDIM < 2, NDIM > 100, MINPTS > MAXPTS, C OR MAXPTS < RULCLS. C*********************************************************************** C***** FOR DOUBLE PRECISION CHANGE REAL TO DOUBLE PRECISION IN THE C NEXT STATEMENT. INTEGER DIVAXO,DIVAXN,DIVFLG,FUNCLS,IFAIL,INDEX1,INDEX2,I,J,K, * L,M,N,LENWRK,MAXPTS,MINPTS,NDIM,RGNSTR,RULCLS,SBRGNS,SBTMPP, * SUBRGN,SUBTMP REAL A(NDIM),B(NDIM),CENTER(100),DF1,DF2,DIF, * DIFMAX,EPS,FINEST,FUNCTN,F1,F2,F3,F4,LAMDA2,LAMDA4,LAMDA5, * NINE, ONE,RATIO,RELERR,RGNCMP,RGNERR,RGNVAL,RGNVOL, * SUM1,SUM2,SUM3,SUM4,SUM5,TWO,TWONDM,WEITP1,WEITP2, * WEITP3,WEITP4,WEIT1,WEIT2,WEIT3,WEIT4,WEIT5,WIDTH(100), * WIDTHL(100),WRKSTR(LENWRK),Z(100) EXTERNAL FUNCTN IFAIL = 3 RELERR = 1 FUNCLS = 0 IF(NDIM.LT.2.OR.NDIM.GT.100) GOTO 300 IF(MINPTS.GT.MAXPTS) GOTO 300 C C***** INITIALISATION OF SUBROUTINE C ONE = 1 TWO = 2 NINE = 9 TWONDM = TWO**NDIM RGNSTR = 2*NDIM+3 DIVAXO = 0 C C***** END SUBROUTINE INITIALISATION C***** BASIC RULE INITIALISATION C LAMDA5 = NINE/19 IF(NDIM.GT.15) GOTO 10 RULCLS = 2**NDIM+2*NDIM*NDIM+2*NDIM+1 LAMDA4 = NINE/10 LAMDA2 = NINE/70 WEIT5 = 1/(3*LAMDA5)**3/TWONDM GOTO 20 10 RULCLS = 1+(NDIM*(12+(NDIM-1)*(6+(NDIM-2)*4)))/3 RATIO = (NDIM-2)/NINE LAMDA4 = (ONE/5-RATIO)/(ONE/3-RATIO/LAMDA5) RATIO = (1-LAMDA4/LAMDA5)*(NDIM-1)*RATIO/6 LAMDA2 = (ONE/7-LAMDA4/5-RATIO)/(ONE/5-LAMDA4/3-RATIO/LAMDA5) WEIT5 = 1/(6*LAMDA5)**3 20 WEIT4 = (ONE/15-LAMDA5/9)/(4*(LAMDA4-LAMDA5)*LAMDA4**2) WEIT3 = (ONE/7-(LAMDA5+LAMDA2)/5+LAMDA5*LAMDA2/3)/ * (2*LAMDA4*(LAMDA4-LAMDA5)*(LAMDA4-LAMDA2))-2*(NDIM-1)*WEIT4 WEIT2 = (ONE/7-(LAMDA5+LAMDA4)/5+LAMDA5*LAMDA4/3)/ * (2*LAMDA2*(LAMDA2-LAMDA5)*(LAMDA2-LAMDA4)) IF(NDIM.GT.15) WEIT1 = 1-2*NDIM*(WEIT2+WEIT3+(NDIM-1)* * (WEIT4+2*(NDIM-2)*WEIT5/3)) IF(NDIM.LT.16) WEIT1 = 1-2*NDIM*(WEIT2+WEIT3+(NDIM-1)* * WEIT4)-TWONDM*WEIT5 WEITP4 = 1/(6*LAMDA4)**2 WEITP3 = (ONE/5-LAMDA2/3)/(2*LAMDA4*(LAMDA4-LAMDA2))- * 2*(NDIM-1)*WEITP4 WEITP2 = (ONE/5-LAMDA4/3)/(TWO*LAMDA2*(LAMDA2-LAMDA4)) WEITP1 = 1-2*NDIM*(WEITP2+WEITP3+(NDIM-1)*WEITP4) RATIO = LAMDA2/LAMDA4 LAMDA5 = SQRT(LAMDA5) LAMDA4 = SQRT(LAMDA4) LAMDA2 = SQRT(LAMDA2) IF(MAXPTS.LT.RULCLS) RETURN C C***** END BASIC RULE INITIALISATION IF(MINPTS.LT.0) SBRGNS = WRKSTR(LENWRK-1) IF(MINPTS.LT.0) GOTO 280 DO 30 J = 1,NDIM WIDTH(J) = (B(J)-A(J))/2 30 CENTER(J) = A(J)+WIDTH(J) FINEST = 0 WRKSTR(LENWRK) = 0 DIVFLG = 1 SUBRGN = RGNSTR SBRGNS = RGNSTR C C***** BEGIN BASIC RULE 40 RGNVOL = TWONDM DO 50 J = 1,NDIM RGNVOL = RGNVOL*WIDTH(J) 50 Z(J) = CENTER(J) SUM1 = FUNCTN(NDIM,Z) C***** COMPUTE SYMMETRIC SUMS OF FUNCTN(LAMDA2,0,0,...,0) AND C FUNCTN(LAMDA4,0,0,...,0), AND MAXIMUM FOURTH DIFFERENCE DIFMAX = -1 SUM2 = 0 SUM3 = 0 DO 60 J = 1,NDIM Z(J) = CENTER(J)-LAMDA2*WIDTH(J) F1 = FUNCTN(NDIM,Z) Z(J) = CENTER(J)+LAMDA2*WIDTH(J) F2 = FUNCTN(NDIM,Z) WIDTHL(J) = LAMDA4*WIDTH(J) Z(J) = CENTER(J)-WIDTHL(J) F3 = FUNCTN(NDIM,Z) Z(J) = CENTER(J)+WIDTHL(J) F4 = FUNCTN(NDIM,Z) SUM2 = SUM2+F1+F2 SUM3 = SUM3+F3+F4 DF1 = F1+F2-2*SUM1 DF2 = F3+F4-2*SUM1 DIF = ABS(DF1-RATIO*DF2) IF(DIF.LE.DIFMAX) GO TO 60 DIFMAX = DIF DIVAXN = J 60 Z(J) = CENTER(J) IF(SUM1.EQ.SUM1+DIFMAX/8) DIVAXN = MOD(DIVAXO,NDIM)+1 C***** COMPUTE SYMMETRIC SUM OF FUNCTN(LAMDA4,LAMDA4,0,0,...,0) SUM4 = 0 DO 90 J = 2,NDIM DO 80 K = J,NDIM DO 70 L = 1,2 WIDTHL(J-1) = -WIDTHL(J-1) Z(J-1) = CENTER(J-1)+WIDTHL(J-1) DO 70 M = 1,2 WIDTHL(K) = -WIDTHL(K) Z(K) = CENTER(K)+WIDTHL(K) 70 SUM4 = SUM4+FUNCTN(NDIM,Z) 80 Z(K) = CENTER(K) 90 Z(J-1) = CENTER(J-1) C***** IF NDIM < 16 COMPUTE SYMMETRIC SUM OF C***** FUNCTN(LAMDA5,LAMDA5,...,LAMDA5) SUM5 = 0 IF(NDIM.GT.15) GOTO 130 DO 100 J = 1,NDIM WIDTHL(J) = -LAMDA5*WIDTH(J) 100 Z(J) = CENTER(J)+WIDTHL(J) 110 SUM5 = SUM5+FUNCTN(NDIM,Z) DO 120 J = 1,NDIM WIDTHL(J) = -WIDTHL(J) Z(J) = CENTER(J)+WIDTHL(J) IF(WIDTHL(J).GT.0) GOTO 110 120 CONTINUE GOTO 190 C***** IF NDIM > 15 COMPUTE SYMMETRIC SUM OF C***** FUNCTN(LAMDA5,LAMDA5,LAMDA5,0,0,...,0) 130 DO 140 J = 1,NDIM 140 WIDTHL(J) = LAMDA5*WIDTH(J) DO 180 I = 3,NDIM DO 170 J = I,NDIM DO 160 K = J,NDIM DO 150 L = 1,2 WIDTHL(I-2) = -WIDTHL(I-2) Z(I-2) = CENTER(I-2)+WIDTHL(I-2) DO 150 M = 1,2 WIDTHL(J-1) = -WIDTHL(J-1) Z(J-1) = CENTER(J-1)+WIDTHL(J-1) DO 150 N = 1,2 WIDTHL(K) = -WIDTHL(K) Z(K) = CENTER(K)+WIDTHL(K) 150 SUM5 = SUM5+FUNCTN(NDIM,Z) 160 Z(K) = CENTER(K) 170 Z(J-1) = CENTER(J-1) 180 Z(I-2) = CENTER(I-2) C***** COMPUTE FIFTH AND SEVENTH DEGREE RULES AND ERROR 190 RGNCMP = RGNVOL*(WEITP1*SUM1+WEITP2*SUM2+WEITP3*SUM3+WEITP4*SUM4) RGNVAL = RGNVOL*(WEIT1*SUM1+WEIT2*SUM2+WEIT3*SUM3+WEIT4*SUM4+ * WEIT5*SUM5) RGNERR = ABS(RGNVAL-RGNCMP) C C***** END BASIC RULE FINEST = FINEST+RGNVAL WRKSTR(LENWRK) = WRKSTR(LENWRK)+RGNERR FUNCLS = FUNCLS+RULCLS C C***** PLACE RESULTS OF BASIC RULE INTO PARTIALLY ORDERED LIST C***** ACCORDING TO SUBREGION ERROR IF(DIVFLG.EQ.1) GO TO 230 C C***** WHEN DIVFLG = 0 START AT TOP OF LIST AND MOVE DOWN LIST TREE TO C FIND CORRECT POSITION FOR RESULTS FROM FIRST HALF OF RECENTLY C DIVIDED SUBREGION 200 SUBTMP = 2*SUBRGN IF(SUBTMP.GT.SBRGNS) GO TO 250 IF(SUBTMP.EQ.SBRGNS) GO TO 210 SBTMPP = SUBTMP+RGNSTR IF(WRKSTR(SUBTMP).LT.WRKSTR(SBTMPP)) SUBTMP = SBTMPP 210 IF(RGNERR.GE.WRKSTR(SUBTMP)) GO TO 250 DO 220 K = 1,RGNSTR INDEX1 = SUBRGN-K+1 INDEX2 = SUBTMP-K+1 220 WRKSTR(INDEX1) = WRKSTR(INDEX2) SUBRGN = SUBTMP GOTO 200 C C***** WHEN DIVFLG = 1 START AT BOTTOM RIGHT BRANCH AND MOVE UP LIST C TREE TO FIND CORRECT POSITION FOR RESULTS FROM SECOND HALF OF C RECENTLY DIVIDED SUBREGION 230 SUBTMP = (SUBRGN/(RGNSTR*2))*RGNSTR IF(SUBTMP.LT.RGNSTR) GO TO 250 IF(RGNERR.LE.WRKSTR(SUBTMP)) GO TO 250 DO 240 K = 1,RGNSTR INDEX1 = SUBRGN-K+1 INDEX2 = SUBTMP-K+1 240 WRKSTR(INDEX1) = WRKSTR(INDEX2) SUBRGN = SUBTMP GOTO 230 C***** STORE RESULTS OF BASIC RULE IN CORRECT POSITION IN LIST 250 WRKSTR(SUBRGN) = RGNERR WRKSTR(SUBRGN-1) = RGNVAL WRKSTR(SUBRGN-2) = DIVAXN DO 260 J = 1,NDIM SUBTMP = SUBRGN-2*(J+1) WRKSTR(SUBTMP+1) = CENTER(J) 260 WRKSTR(SUBTMP) = WIDTH(J) IF(DIVFLG.EQ.1) GO TO 270 C***** WHEN DIVFLG = 0 PREPARE FOR SECOND APPLICATION OF BASIC RULE CENTER(DIVAXO) = CENTER(DIVAXO)+2*WIDTH(DIVAXO) SBRGNS = SBRGNS+RGNSTR SUBRGN = SBRGNS DIVFLG = 1 C***** LOOP BACK TO APPLY BASIC RULE TO OTHER HALF OF SUBREGION GO TO 40 C C***** END ORDERING AND STORAGE OF BASIC RULE RESULTS C***** MAKE CHECKS FOR POSSIBLE TERMINATION OF ROUTINE C 270 RELERR = ONE IF(WRKSTR(LENWRK).LE.0) WRKSTR(LENWRK) = 0 IF(ABS(FINEST).NE.0) RELERR = WRKSTR(LENWRK)/ABS(FINEST) IF(RELERR.GT.1) RELERR = 1 IF(SBRGNS+RGNSTR.GT.LENWRK-2) IFAIL = 2 IF(FUNCLS+2*RULCLS.GT.MAXPTS) IFAIL = 1 IF(RELERR.LT.EPS.AND.FUNCLS.GE.MINPTS) IFAIL = 0 IF(IFAIL.LT.3) GOTO 300 C C***** PREPARE TO USE BASIC RULE ON EACH HALF OF SUBREGION WITH LARGEST C ERROR 280 DIVFLG = 0 SUBRGN = RGNSTR WRKSTR(LENWRK) = WRKSTR(LENWRK)-WRKSTR(SUBRGN) FINEST = FINEST-WRKSTR(SUBRGN-1) DIVAXO = WRKSTR(SUBRGN-2) DO 290 J = 1,NDIM SUBTMP = SUBRGN-2*(J+1) CENTER(J) = WRKSTR(SUBTMP+1) 290 WIDTH(J) = WRKSTR(SUBTMP) WIDTH(DIVAXO) = WIDTH(DIVAXO)/2 CENTER(DIVAXO) = CENTER(DIVAXO)-WIDTH(DIVAXO) C C***** LOOP BACK TO APPLY BASIC RULE C GOTO 40 C C***** TERMINATION POINT C 300 MINPTS = FUNCLS WRKSTR(LENWRK-1) = SBRGNS END