C ALGORITHM 612 COLLECTED ALGORITHMS FROM ACM. C ALGORITHM APPEARED IN ACM-TRANS. MATH. SOFTWARE, VOL.10, NO. 1, C MAR., 1984, P. 1-16, P. 17-22. DOUBLE PRECISION FUNCTION D1MACH(I) C***BEGIN PROLOGUE D1MACH C***DATE WRITTEN 750101 (YYMMDD) C***REVISION DATE 820801 (YYMMDD) C***CATEGORY NO. R1 C***KEYWORDS MACHINE CONSTANTS C***AUTHOR FOX, P. A., (BELL LABS) C HALL, A. D., (BELL LABS) C SCHRYER, N. L., (BELL LABS) C***PURPOSE RETURNS DOUBLE PRECISION MACHINE DEPENDENT CONSTANTS C***DESCRIPTION C D1MACH CAN BE USED TO OBTAIN MACHINE-DEPENDENT PARAMETERS C FOR THE LOCAL MACHINE ENVIRONMENT. IT IS A FUNCTION C SUBPROGRAM WITH ONE (INPUT) ARGUMENT, AND CAN BE CALLED C AS FOLLOWS, FOR EXAMPLE C C D = D1MACH(I) C C WHERE I=1,...,5. THE (OUTPUT) VALUE OF D ABOVE IS C DETERMINED BY THE (INPUT) VALUE OF I. THE RESULTS FOR C VARIOUS VALUES OF I ARE DISCUSSED BELOW. C C DOUBLE-PRECISION MACHINE CONSTANTS C D1MACH( 1) = B**(EMIN-1), THE SMALLEST POSITIVE MAGNITUDE. C D1MACH( 2) = B**EMAX*(1 - B**(-T)), THE LARGEST MAGNITUDE. C D1MACH( 3) = B**(-T), THE SMALLEST RELATIVE SPACING. C D1MACH( 4) = B**(1-T), THE LARGEST RELATIVE SPACING. C D1MACH( 5) = LOG10(B) C***REFERENCES FOX P.A., HALL A.D., SCHRYER N.L.,*FRAMEWORK FOR A C PORTABLE LIBRARY*, ACM TRANSACTIONS ON MATHEMATICAL C SOFTWARE, VOL. 4, NO. 2, JUNE 1978, PP. 177-188. C***ROUTINES CALLED XERROR C***END PROLOGUE D1MACH C INTEGER SMALL(4) INTEGER LARGE(4) INTEGER RIGHT(4) INTEGER DIVER(4) INTEGER LOG10(4) C DOUBLE PRECISION DMACH(5) C EQUIVALENCE (DMACH(1),SMALL(1)) EQUIVALENCE (DMACH(2),LARGE(1)) EQUIVALENCE (DMACH(3),RIGHT(1)) EQUIVALENCE (DMACH(4),DIVER(1)) EQUIVALENCE (DMACH(5),LOG10(1)) C C MACHINE CONSTANTS FOR THE BURROUGHS 1700 SYSTEM. C C DATA SMALL(1) / ZC00800000 / C DATA SMALL(2) / Z000000000 / C C DATA LARGE(1) / ZDFFFFFFFF / C DATA LARGE(2) / ZFFFFFFFFF / C C DATA RIGHT(1) / ZCC5800000 / C DATA RIGHT(2) / Z000000000 / C C DATA DIVER(1) / ZCC6800000 / C DATA DIVER(2) / Z000000000 / C C DATA LOG10(1) / ZD00E730E7 / C DATA LOG10(2) / ZC77800DC0 / C C MACHINE CONSTANTS FOR THE BURROUGHS 5700 SYSTEM. C C DATA SMALL(1) / O1771000000000000 / C DATA SMALL(2) / O0000000000000000 / C C DATA LARGE(1) / O0777777777777777 / C DATA LARGE(2) / O0007777777777777 / C C DATA RIGHT(1) / O1461000000000000 / C DATA RIGHT(2) / O0000000000000000 / C C DATA DIVER(1) / O1451000000000000 / C DATA DIVER(2) / O0000000000000000 / C C DATA LOG10(1) / O1157163034761674 / C DATA LOG10(2) / O0006677466732724 / C C MACHINE CONSTANTS FOR THE BURROUGHS 6700/7700 SYSTEMS. C C DATA SMALL(1) / O1771000000000000 / C DATA SMALL(2) / O7770000000000000 / C C DATA LARGE(1) / O0777777777777777 / C DATA LARGE(2) / O7777777777777777 / C C DATA RIGHT(1) / O1461000000000000 / C DATA RIGHT(2) / O0000000000000000 / C C DATA DIVER(1) / O1451000000000000 / C DATA DIVER(2) / O0000000000000000 / C C DATA LOG10(1) / O1157163034761674 / C DATA LOG10(2) / O0006677466732724 / C C MACHINE CONSTANTS FOR THE CDC 6000/7000 SERIES. C C DATA SMALL(1) / 00604000000000000000B / C DATA SMALL(2) / 00000000000000000000B / C C DATA LARGE(1) / 37767777777777777777B / C DATA LARGE(2) / 37167777777777777777B / C C DATA RIGHT(1) / 15604000000000000000B / C DATA RIGHT(2) / 15000000000000000000B / C C DATA DIVER(1) / 15614000000000000000B / C DATA DIVER(2) / 15010000000000000000B / C C DATA LOG10(1) / 17164642023241175717B / C DATA LOG10(2) / 16367571421742254654B / C C MACHINE CONSTANTS FOR THE CRAY 1 C C DATA SMALL(1) / 200004000000000000000B / C DATA SMALL(2) / 00000000000000000000B / C C DATA LARGE(1) / 577777777777777777777B / C DATA LARGE(2) / 000007777777777777777B / C C DATA RIGHT(1) / 377214000000000000000B / C DATA RIGHT(2) / 000000000000000000000B / C C DATA DIVER(1) / 377224000000000000000B / C DATA DIVER(2) / 000000000000000000000B / C C DATA LOG10(1) / 377774642023241175717B / C DATA LOG10(2) / 000007571421742254654B / C C MACHINE CONSTANTS FOR THE DATA GENERAL ECLIPSE S/200 C C NOTE - IT MAY BE APPROPRIATE TO INCLUDE THE FOLLOWING CARD - C STATIC DMACH(5) C C DATA SMALL/20K,3*0/,LARGE/77777K,3*177777K/ C DATA RIGHT/31420K,3*0/,DIVER/32020K,3*0/ C DATA LOG10/40423K,42023K,50237K,74776K/ C C MACHINE CONSTANTS FOR THE HARRIS 220 C C DATA SMALL(1),SMALL(2) / [20000000, [00000201 / C DATA LARGE(1),LARGE(2) / [37777777, [37777577 / C DATA RIGHT(1),RIGHT(2) / [20000000, [00000333 / C DATA DIVER(1),DIVER(2) / [20000000, [00000334 / C DATA LOG10(1),LOG10(2) / [23210115, [10237777 / C C MACHINE CONSTANTS FOR THE HONEYWELL 600/6000 SERIES. C C DATA SMALL(1),SMALL(2) / O402400000000, O000000000000 / C DATA LARGE(1),LARGE(2) / O376777777777, O777777777777 / C DATA RIGHT(1),RIGHT(2) / O604400000000, O000000000000 / C DATA DIVER(1),DIVER(2) / O606400000000, O000000000000 / C DATA LOG10(1),LOG10(2) / O776464202324, O117571775714 / C C MACHINE CONSTANTS FOR THE HP 2100 C THREE WORD DOUBLE PRECISION OPTION WITH FTN4 C C DATA SMALL(1), SMALL(2), SMALL(3) / 40000B, 0, 1 / C DATA LARGE(1), LARGE(2), LARGE(3) / 77777B, 177777B, 177776B / C DATA RIGHT(1), RIGHT(2), RIGHT(3) / 40000B, 0, 265B / C DATA DIVER(1), DIVER(2), DIVER(3) / 40000B, 0, 276B / C DATA LOG10(1), LOG10(2), LOG10(3) / 46420B, 46502B, 77777B / C C C MACHINE CONSTANTS FOR THE HP 2100 C FOUR WORD DOUBLE PRECISION OPTION WITH FTN4 C C DATA SMALL(1), SMALL(2) / 40000B, 0 / C DATA SMALL(3), SMALL(4) / 0, 1 / C DATA LARGE(1), LARGE(2) / 77777B, 177777B / C DATA LARGE(3), LARGE(4) / 177777B, 177776B / C DATA RIGHT(1), RIGHT(2) / 40000B, 0 / C DATA RIGHT(3), RIGHT(4) / 0, 225B / C DATA DIVER(1), DIVER(2) / 40000B, 0 / C DATA DIVER(3), DIVER(4) / 0, 227B / C DATA LOG10(1), LOG10(2) / 46420B, 46502B / C DATA LOG10(3), LOG10(4) / 76747B, 176377B / C C C MACHINE CONSTANTS FOR THE IBM 360/370 SERIES, C THE XEROX SIGMA 5/7/9, THE SEL SYSTEMS 85/86, AND C THE PERKIN ELMER (INTERDATA) 7/32. C C DATA SMALL(1),SMALL(2) / Z00100000, Z00000000 / C DATA LARGE(1),LARGE(2) / Z7FFFFFFF, ZFFFFFFFF / C DATA RIGHT(1),RIGHT(2) / Z33100000, Z00000000 / C DATA DIVER(1),DIVER(2) / Z34100000, Z00000000 / C DATA LOG10(1),LOG10(2) / Z41134413, Z509F79FF / C C MACHINE CONSTANTS FOR THE PDP-10 (KA PROCESSOR). C C DATA SMALL(1),SMALL(2) / "033400000000, "000000000000 / C DATA LARGE(1),LARGE(2) / "377777777777, "344777777777 / C DATA RIGHT(1),RIGHT(2) / "113400000000, "000000000000 / C DATA DIVER(1),DIVER(2) / "114400000000, "000000000000 / C DATA LOG10(1),LOG10(2) / "177464202324, "144117571776 / C C MACHINE CONSTANTS FOR THE PDP-10 (KI PROCESSOR). C C DATA SMALL(1),SMALL(2) / "000400000000, "000000000000 / C DATA LARGE(1),LARGE(2) / "377777777777, "377777777777 / C DATA RIGHT(1),RIGHT(2) / "103400000000, "000000000000 / C DATA DIVER(1),DIVER(2) / "104400000000, "000000000000 / C DATA LOG10(1),LOG10(2) / "177464202324, "476747767461 / C C MACHINE CONSTANTS FOR PDP-11 FORTRAN SUPPORTING C 32-BIT INTEGERS (EXPRESSED IN INTEGER AND OCTAL). C C DATA SMALL(1),SMALL(2) / 8388608, 0 / C DATA LARGE(1),LARGE(2) / 2147483647, -1 / C DATA RIGHT(1),RIGHT(2) / 612368384, 0 / C DATA DIVER(1),DIVER(2) / 620756992, 0 / C DATA LOG10(1),LOG10(2) / 1067065498, -2063872008 / C C DATA SMALL(1),SMALL(2) / O00040000000, O00000000000 / C DATA LARGE(1),LARGE(2) / O17777777777, O37777777777 / C DATA RIGHT(1),RIGHT(2) / O04440000000, O00000000000 / C DATA DIVER(1),DIVER(2) / O04500000000, O00000000000 / C DATA LOG10(1),LOG10(2) / O07746420232, O20476747770 / C C MACHINE CONSTANTS FOR PDP-11 FORTRAN SUPPORTING C 16-BIT INTEGERS (EXPRESSED IN INTEGER AND OCTAL). C C DATA SMALL(1),SMALL(2) / 128, 0 / C DATA SMALL(3),SMALL(4) / 0, 0 / C C DATA LARGE(1),LARGE(2) / 32767, -1 / C DATA LARGE(3),LARGE(4) / -1, -1 / C C DATA RIGHT(1),RIGHT(2) / 9344, 0 / C DATA RIGHT(3),RIGHT(4) / 0, 0 / C C DATA DIVER(1),DIVER(2) / 9472, 0 / C DATA DIVER(3),DIVER(4) / 0, 0 / C C DATA LOG10(1),LOG10(2) / 16282, 8346 / C DATA LOG10(3),LOG10(4) / -31493, -12296 / C C DATA SMALL(1),SMALL(2) / O000200, O000000 / C DATA SMALL(3),SMALL(4) / O000000, O000000 / C C DATA LARGE(1),LARGE(2) / O077777, O177777 / C DATA LARGE(3),LARGE(4) / O177777, O177777 / C C DATA RIGHT(1),RIGHT(2) / O022200, O000000 / C DATA RIGHT(3),RIGHT(4) / O000000, O000000 / C C DATA DIVER(1),DIVER(2) / O022400, O000000 / C DATA DIVER(3),DIVER(4) / O000000, O000000 / C C DATA LOG10(1),LOG10(2) / O037632, O020232 / C DATA LOG10(3),LOG10(4) / O102373, O147770 / C C MACHINE CONSTANTS FOR THE UNIVAC 1100 SERIES. FTN COMPILER C C DATA SMALL(1),SMALL(2) / O000040000000, O000000000000 / C DATA LARGE(1),LARGE(2) / O377777777777, O777777777777 / C DATA RIGHT(1),RIGHT(2) / O170540000000, O000000000000 / C DATA DIVER(1),DIVER(2) / O170640000000, O000000000000 / C DATA LOG10(1),LOG10(2) / O177746420232, O411757177572 / C C MACHINE CONSTANTS FOR THE UNIVAC 1100 SERIES. FOR COMPILER C C DATA SMALL(1), SMALL(2) / O000040000000, O000000000000 / C DATA LARGE(1), LARGE(2) / O377777777777, O777777777777 / C DATA RIGHT(1), RIGHT(2) / O170540000000, O000000000000 / C DATA DIVER(1), DIVER(2) / O170640000000, O000000000000 / C DATA LOG10(1), LOG10(2) / O177746420232, O411757177572/ C C C MACHINE CONSTANTS FOR VAX 11/780 C (EXPRESSED IN INTEGER AND HEXADECIMAL) C ***THE HEX FORMAT BELOW MAY NOT BE SUITABLE FOR UNIX SYSYEMS*** C *** THE INTEGER FORMAT SHOULD BE OK FOR UNIX SYSTEMS*** C C DATA SMALL(1), SMALL(2) / 128, 0 / C DATA LARGE(1), LARGE(2) / -32769, -1 / C DATA RIGHT(1), RIGHT(2) / 9344, 0 / C DATA DIVER(1), DIVER(2) / 9472, 0 / C DATA LOG10(1), LOG10(2) / 546979738, -805665541 / C C DATA SMALL(1), SMALL(2) / Z00000080, Z00000000 / C DATA LARGE(1), LARGE(2) / ZFFFF7FFF, ZFFFFFFFF / C DATA RIGHT(1), RIGHT(2) / Z00002480, Z00000000 / C DATA DIVER(1), DIVER(2) / Z00002500, Z00000000 / C DATA LOG10(1), LOG10(2) / Z209A3F9A, ZCFFA84FB / C C***FIRST EXECUTABLE STATEMENT D1MACH C IF (I .LT. 1 .OR. I .GT. 5) 1 CALL XERROR(25HD1MACH -- I OUT OF BOUNDS,25,1,2) C D1MACH = DMACH(I) RETURN C END REAL FUNCTION R1MACH(I) C***BEGIN PROLOGUE R1MACH C***DATE WRITTEN 790101 (YYMMDD) C***REVISION DATE 820801 (YYMMDD) C***CATEGORY NO. R1 C***KEYWORDS MACHINE CONSTANTS C***AUTHOR FOX, P. A., (BELL LABS) C HALL, A. D., (BELL LABS) C SCHRYER, N. L., (BELL LABS) C***PURPOSE RETURNS SINGLE PRECISION MACHINE DEPENDENT CONSTANTS C***DESCRIPTION C R1MACH CAN BE USED TO OBTAIN MACHINE-DEPENDENT PARAMETERS C FOR THE LOCAL MACHINE ENVIRONMENT. IT IS A FUNCTION C SUBROUTINE WITH ONE (INPUT) ARGUMENT, AND CAN BE CALLED C AS FOLLOWS, FOR EXAMPLE C C A = R1MACH(I) C C WHERE I=1,...,5. THE (OUTPUT) VALUE OF A ABOVE IS C DETERMINED BY THE (INPUT) VALUE OF I. THE RESULTS FOR C VARIOUS VALUES OF I ARE DISCUSSED BELOW. C C SINGLE-PRECISION MACHINE CONSTANTS C R1MACH(1) = B**(EMIN-1), THE SMALLEST POSITIVE MAGNITUDE. C R1MACH(2) = B**EMAX*(1 - B**(-T)), THE LARGEST MAGNITUDE. C R1MACH(3) = B**(-T), THE SMALLEST RELATIVE SPACING. C R1MACH(4) = B**(1-T), THE LARGEST RELATIVE SPACING. C R1MACH(5) = LOG10(B) C***REFERENCES FOX, P.A., HALL, A.D., SCHRYER, N.L, *FRAMEWORK FOR C A PORTABLE LIBRARY*, ACM TRANSACTIONS ON MATHE- C MATICAL SOFTWARE, VOL. 4, NO. 2, JUNE 1978, C PP. 177-188. C***ROUTINES CALLED XERROR C***END PROLOGUE R1MACH C INTEGER SMALL(2) INTEGER LARGE(2) INTEGER RIGHT(2) INTEGER DIVER(2) INTEGER LOG10(2) C REAL RMACH(5) C EQUIVALENCE (RMACH(1),SMALL(1)) EQUIVALENCE (RMACH(2),LARGE(1)) EQUIVALENCE (RMACH(3),RIGHT(1)) EQUIVALENCE (RMACH(4),DIVER(1)) EQUIVALENCE (RMACH(5),LOG10(1)) C C MACHINE CONSTANTS FOR THE BURROUGHS 1700 SYSTEM. C C DATA RMACH(1) / Z400800000 / C DATA RMACH(2) / Z5FFFFFFFF / C DATA RMACH(3) / Z4E9800000 / C DATA RMACH(4) / Z4EA800000 / C DATA RMACH(5) / Z500E730E8 / C C MACHINE CONSTANTS FOR THE BURROUGHS 5700/6700/7700 SYSTEMS. C C DATA RMACH(1) / O1771000000000000 / C DATA RMACH(2) / O0777777777777777 / C DATA RMACH(3) / O1311000000000000 / C DATA RMACH(4) / O1301000000000000 / C DATA RMACH(5) / O1157163034761675 / C C MACHINE CONSTANTS FOR THE CDC 6000/7000 SERIES. C C DATA RMACH(1) / 00014000000000000000B / C DATA RMACH(2) / 37767777777777777777B / C DATA RMACH(3) / 16404000000000000000B / C DATA RMACH(4) / 16414000000000000000B / C DATA RMACH(5) / 17164642023241175720B / C C MACHINE CONSTANTS FOR THE CRAY 1 C C DATA RMACH(1) / 200004000000000000000B / C DATA RMACH(2) / 577777777777777777777B / C DATA RMACH(3) / 377214000000000000000B / C DATA RMACH(4) / 377224000000000000000B / C DATA RMACH(5) / 377774642023241175720B / C C MACHINE CONSTANTS FOR THE DATA GENERAL ECLIPSE S/200 C C NOTE - IT MAY BE APPROPRIATE TO INCLUDE THE FOLLOWING CARD - C STATIC RMACH(5) C C DATA SMALL/20K,0/,LARGE/77777K,177777K/ C DATA RIGHT/35420K,0/,DIVER/36020K,0/ C DATA LOG10/40423K,42023K/ C C MACHINE CONSTANTS FOR THE HARRIS 220 C C DATA SMALL(1),SMALL(2) / [20000000, [00000201 / C DATA LARGE(1),LARGE(2) / [37777777, [00000177 / C DATA RIGHT(1),RIGHT(2) / [20000000, [00000352 / C DATA DIVER(1),DIVER(2) / [20000000, [00000353 / C DATA LOG10(1),LOG10(2) / [23210115, [00000377 / C C MACHINE CONSTANTS FOR THE HONEYWELL 600/6000 SERIES. C C DATA RMACH(1) / O402400000000 / C DATA RMACH(2) / O376777777777 / C DATA RMACH(3) / O714400000000 / C DATA RMACH(4) / O716400000000 / C DATA RMACH(5) / O776464202324 / C C MACHINE CONSTANTS FOR THE HP 2100 C C 3 WORD DOUBLE PRECISION WITH FTN4 C C DATA SMALL(1), SMALL(2) / 40000B, 1 / C DATA LARGE(1), LARGE(2) / 77777B, 177776B / C DATA RIGHT(1), RIGHT(2) / 40000B, 325B / C DATA DIVER(1), DIVER(2) / 40000B, 327B / C DATA LOG10(1), LOG10(2) / 46420B, 46777B / C C MACHINE CONSTANTS FOR THE HP 2100 C 4 WORD DOUBLE PRECISION WITH FTN4 C C DATA SMALL(1), SMALL(2) / 40000B, 1 / C DATA LARGE91), LARGE(2) / 77777B, 177776B / C DATA RIGHT(1), RIGHT(2) / 40000B, 325B / C DATA DIVER(1), DIVER(2) / 40000B, 327B / C DATA LOG10(1), LOG10(2) / 46420B, 46777B / C C MACHINE CONSTANTS FOR THE IBM 360/370 SERIES, C THE XEROX SIGMA 5/7/9, THE SEL SYSTEMS 85/86 AND C THE PERKIN ELMER (INTERDATA) 7/32. C C DATA RMACH(1) / Z00100000 / C DATA RMACH(2) / Z7FFFFFFF / C DATA RMACH(3) / Z3B100000 / C DATA RMACH(4) / Z3C100000 / C DATA RMACH(5) / Z41134413 / C C MACHINE CONSTANTS FOR THE PDP-10 (KA OR KI PROCESSOR). C C DATA RMACH(1) / "000400000000 / C DATA RMACH(2) / "377777777777 / C DATA RMACH(3) / "146400000000 / C DATA RMACH(4) / "147400000000 / C DATA RMACH(5) / "177464202324 / C C MACHINE CONSTANTS FOR PDP-11 FORTRAN SUPPORTING C 32-BIT INTEGERS (EXPRESSED IN INTEGER AND OCTAL). C C DATA SMALL(1) / 8388608 / C DATA LARGE(1) / 2147483647 / C DATA RIGHT(1) / 880803840 / C DATA DIVER(1) / 889192448 / C DATA LOG10(1) / 1067065499 / C C DATA RMACH(1) / O00040000000 / C DATA RMACH(2) / O17777777777 / C DATA RMACH(3) / O06440000000 / C DATA RMACH(4) / O06500000000 / C DATA RMACH(5) / O07746420233 / C C MACHINE CONSTANTS FOR PDP-11 FORTRAN SUPPORTING C 16-BIT INTEGERS (EXPRESSED IN INTEGER AND OCTAL). C C DATA SMALL(1),SMALL(2) / 128, 0 / C DATA LARGE(1),LARGE(2) / 32767, -1 / C DATA RIGHT(1),RIGHT(2) / 13440, 0 / C DATA DIVER(1),DIVER(2) / 13568, 0 / C DATA LOG10(1),LOG10(2) / 16282, 8347 / C C DATA SMALL(1),SMALL(2) / O000200, O000000 / C DATA LARGE(1),LARGE(2) / O077777, O177777 / C DATA RIGHT(1),RIGHT(2) / O032200, O000000 / C DATA DIVER(1),DIVER(2) / O032400, O000000 / C DATA LOG10(1),LOG10(2) / O037632, O020233 / C C MACHINE CONSTANTS FOR THE UNIVAC 1100 SERIES. C C DATA RMACH(1) / O000400000000 / C DATA RMACH(2) / O377777777777 / C DATA RMACH(3) / O146400000000 / C DATA RMACH(4) / O147400000000 / C DATA RMACH(5) / O177464202324 / C C MACHINE CONSTANTS FOR THE VAX 11/780 C (EXPRESSED IN INTEGER AND HEXADECIMAL) C ***THE HEX FORMAT BELOW MAY NOT BE SUITABLE FOR UNIX SYSTEMS*** C *** THE INTEGER FORMAT SHOULD BE OK FOR UNIX SYSTEMS*** C C DATA SMALL(1) / 128 / C DATA LARGE(1) / -32769 / C DATA RIGHT(1) / 13440 / C DATA DIVER(1) / 13568 / C DATA LOG10(1) / 547045274 / C C DATA SMALL(1) / Z00000080 / C DATA LARGE(1) / ZFFFF7FFF / C DATA RIGHT(1) / Z00003480 / C DATA DIVER(1) / Z00003500 / C DATA LOG10(1) / Z209B3F9A / C C***FIRST EXECUTABLE STATEMENT R1MACH C IF (I .LT. 1 .OR. I .GT. 5) 1 CALL XERROR (25HR1MACH -- I OUT OF BOUNDS,25,1,2) C R1MACH = RMACH(I) RETURN C END C.................................................................. MAN 10 C MAN 20 C INTEGRATION OVER A TRIANGLE TEST TRIEX MAN 30 C MAN 40 C THE FORM IN WHICH THE TEST INTEGRALS ARE SUPPLIED TO THE MAN 50 C INTEGRATOR CORRESPONDS TO MAN 60 C I(PA,PB) I(PC*X+PD,PG*X+PH) F(X,Y) DXDY MAN 70 C WITH MAN 80 C F(X,Y) = G(U(X),V(X,Y))*V1(X)/D MAN 90 C WHERE MAN 100 C D = PB-PA MAN 110 C V1(X) = (1-(X-PA)/D)/((PG-PC)*X+PH-PD) MAN 120 C U(X) = (X-PA)/D MAN 130 C V(X,Y) = V1(X)*(Y-PC*X-PD) MAN 140 C HENCE THIS IS THE TRANSFORMED VERSION OF THE INTEGRAL OVER MAN 150 C THE UNIT TRIANGLE GIVEN BY MAN 160 C I(0,1) I(0,1-X) G(X,Y) DXDY MAN 170 C THE FOLLOWING LIST PROVIDES, FOR EACH TEST INTEGRAL, THE MAN 180 C PARAMETERS PA, PB, PC, PD, PG, PH, AND THE VERTICES OF THE MAN 190 C TRIANGLE OVER WHICH IS THUS INTEGRATED: VER(1,K) = ABSCISSAE MAN 200 C AND VER(2,K) = ORDINATES OF THE VERTICES INDEXED K. MAN 210 C MAN 220 C INTEGRAL NUMBER 1 MAN 230 C G(X,Y) = X**(-0.2) MAN 240 C PA=1 PB=5 PC=0 PD=0 PG=0.25 PH=-0.25 MAN 250 C VER(1,1)=1 VER(1,2)=5 VER(1,3)=5 MAN 260 C VER(2,1)=0 VER(2,2)=0 VER(2,3)=1 MAN 270 C MAN 280 C INTEGRAL NUMBER 2 MAN 290 C G(X,Y) = (X+Y)**(-0.2) MAN 300 C PA=0 PB=1 PC=0 PD=0 PG=-1 PH=1 MAN 310 C VER(1,1)=0 VER(1,2)=1 VER(1,3)=0 MAN 320 C VER(2,1)=0 VER(2,2)=0 VER(2,3)=1 MAN 330 C MAN 340 C INTEGRAL NUMBER 3 MAN 350 C G(X,Y) = (1-X-Y)**(-0.2) MAN 360 C PA=-1 PB=3 PC=0.25 PD=-2.75 PD=-1 PH=1 MAN 370 C VER(1,1)=-1 VER(1,2)=3 VER(1,3)=-1 MAN 380 C VER(2,1)=-3 VER(2,2)=-2 VER(2,3)=2 MAN 390 C MAN 400 C INTEGRAL NUMBER 4 MAN 410 C G(X,Y) = (X*Y)**(-0.2) MAN 420 C PA=0 PB=-7 PC=0 PD=0 PG=-3/7 PH=-3 MAN 430 C VER(1,1)=0 VER(1,2)=-7 VER(1,3)=0 MAN 440 C VER(2,1)=0 VER(2,2)=0 VER(2,3)=-3 MAN 450 C MAN 460 C.................................................................. MAN 470 DOUBLE PRECISION ABSERR, C, D, EPSABS, EPSREL, ER, ES, EXA, MAN 480 * EXACT, F, G, H, P, PA, PB, PC, PD, PG, PH, Q, R, RESULT, VER MAN 490 INTEGER IER, INDACC, INUM, KOUNT, LIMEV, NOUT, NUM, N2, N2M1 MAN 500 DOUBLE PRECISION DABS MAN 510 DIMENSION C(4), D(4), EXACT(4), G(4), H(4), P(8), Q(8), R(8), MAN 520 * VER(2,3) MAN 530 COMMON /CNT/ KOUNT, NUM MAN 540 COMMON /PAR/ PA, PB, PC, PD, PG, PH MAN 550 EXTERNAL F MAN 560 C MAN 570 C EXACT VALUES MAN 580 C MAN 590 DATA EXACT(1), EXACT(2), EXACT(3), EXACT(4) /0.6944444444444444D+0MAN 600 * 0,0.5555555555555556D+00,0.6944444444444444D+00, MAN 610 * 0.9481026454955768D+00/ MAN 620 C MAN 630 C (P(2*NUM-1),P(2*NUM)) IS 1ST VERTEX FOR INTEGRAL NUM=1,2,3,4 MAN 640 C MAN 650 DATA P(1), P(2), P(3), P(4), P(5), P(6), P(7), P(8) MAN 660 * /1.0D+00,0.0D+00,0.0D+00,0.0D+00,-1.0D+00,-3.0D+00,0.0D+00, MAN 670 * 0.0D+00/ MAN 680 C MAN 690 C (Q(2*NUM-1),Q(2*NUM)) IS 2ND VERTEX FOR INTEGRAL NUM=1,2,3,4 MAN 700 C MAN 710 DATA Q(1), Q(2), Q(3), Q(4), Q(5), Q(6), Q(7), Q(8) MAN 720 * /5.0D+00,0.0D+00,1.0D+00,0.0D+00,3.0D+00,-2.0D+00,-7.0D+00, MAN 730 * 0.0D+00/ MAN 740 C MAN 750 C (R(2*NUM-1),R(2*NUM)) IS 3RD VERTEX FOR INTEGRAL NUM=1,2,3,4 MAN 760 C MAN 770 DATA R(1), R(2), R(3), R(4), R(5), R(6), R(7), R(8) MAN 780 * /5.0D+00,1.0D+00,0.0D+00,1.0D+00,-1.0D+00,2.0D+00,0.0D+00, MAN 790 * -3.0D+00/ MAN 800 C MAN 810 C PARAMETERS MAN 820 C MAN 830 DATA C(1), C(2), C(3), C(4) /0.0D+00,0.0D+00,0.25D+00,0.0D+00/ MAN 840 DATA D(1), D(2), D(3), D(4) /0.0D+00,0.0D+00,-2.75D+00,0.0D+00/ MAN 850 DATA G(1), G(2), G(3), G(4) /0.25D+00,-1.0D+00,-1.0D+00, MAN 860 * -0.428571428571428571D+00/ MAN 870 DATA H(1), H(2), H(3), H(4) /-0.25D+00,1.0D+00,1.0D+00,-3.0D+00/ MAN 880 C MAN 890 C OUTPUT UNIT MAN 900 C MAN 910 DATA NOUT /21/ MAN 920 C MAN 930 C HEADINGS MAN 940 C MAN 950 WRITE (NOUT,99999) MAN 960 99999 FORMAT (1X, 27HINTEGRATION OVER A TRIANGLE/1X, 27(1H.)/) MAN 970 WRITE (NOUT,99998) MAN 980 99998 FORMAT (6X, 8HRELATIVE, 9X, 8HCOMPUTED, 10X, 8HRELATIVE, 2X, MAN 990 * 6HESTIM., 3X, 6HNUMBER/1X, 3HNR., 2X, 8HACCURACY, 9X, 8HINTEGRAL,MAN 1000 * 10X, 5HERROR, 5X, 8HRELATIVE, 1X, 6HOF F , 4H IER/6X, 7HREQUEST,MAN 1010 * 2HED, 36X, 5HERROR, 4X, 5HEVAL.) MAN 1020 C MAN 1030 C INTEGRATION MAN 1040 C MAN 1050 EPSABS = 0.0D+00 MAN 1060 LIMEV = 50000 MAN 1070 DO 40 INUM=1,4 MAN 1080 EXA = EXACT(INUM) MAN 1090 EPSREL = 1.0D-02 MAN 1100 DO 30 INDACC=1,3 MAN 1110 KOUNT = 0 MAN 1120 NUM = INUM MAN 1130 N2 = NUM + NUM MAN 1140 N2M1 = N2 - 1 MAN 1150 VER(1,1) = P(N2M1) MAN 1160 VER(2,1) = P(N2) MAN 1170 VER(1,2) = Q(N2M1) MAN 1180 VER(2,2) = Q(N2) MAN 1190 VER(1,3) = R(N2M1) MAN 1200 VER(2,3) = R(N2) MAN 1210 PA = P(N2M1) MAN 1220 PB = Q(N2M1) MAN 1230 PC = C(NUM) MAN 1240 PD = D(NUM) MAN 1250 PG = G(NUM) MAN 1260 PH = H(NUM) MAN 1270 CALL TRIEX(F, VER, EPSABS, EPSREL, LIMEV, RESULT, ABSERR, IER)MAN 1280 C ....................................................... MAN 1290 ER = DABS((RESULT-EXA)/EXA) MAN 1300 ES = ABSERR/DABS(RESULT) MAN 1310 IF (INDACC.GT.1) GO TO 10 MAN 1320 WRITE (NOUT,99997) NUM, EPSREL, RESULT, ER, ES, KOUNT, IER MAN 1330 99997 FORMAT (/1X, I2, D11.2, D25.16, 2D10.2, I7, I4) MAN 1340 GO TO 20 MAN 1350 10 WRITE (NOUT,99996) EPSREL, RESULT, ER, ES, KOUNT, IER MAN 1360 99996 FORMAT (3X, D11.2, D25.16, 2D10.2, I7, I4) MAN 1370 20 EPSREL = EPSREL*1.0D-02 MAN 1380 30 CONTINUE MAN 1390 WRITE (NOUT,99995) EXA MAN 1400 99995 FORMAT (7X, 7HEXACT =, D25.16) MAN 1410 40 CONTINUE MAN 1420 STOP MAN 1430 END MAN 1440 C F 10 C TEST INTEGRANDS F 20 C F 30 DOUBLE PRECISION FUNCTION F(X, Y) F 40 DOUBLE PRECISION A, D, DRELPR, D1MACH, P, PA, PB, PC, PD, PG, PH, * X, Y, U, V, V1 INTEGER KOUNT, NUM COMMON /CNT/ KOUNT, NUM COMMON /PAR/ PA, PB, PC, PD, PG, PH C C DRELPR IS THE LARGEST RELATIVE SPACING C C CALL FUNCTION SUBPROGRAM D1MACH FROM P O R T MATHEMATICAL C SUBROUTINE LIBRARY (BELL LABORATORIES) C DRELPR = D1MACH(4) P = -0.2D+00 KOUNT = KOUNT + 1 F = 0.0D+00 D = PB - PA U = (X-PA)/D V1 = (1.0D+00-U)/((PG-PC)*X+PH-PD) V = V1*(Y-PC*X-PD) GO TO (10, 20, 30, 40), NUM C SINGULARITY AT U=0, I.E. X=PA C (VERTEX OF 1ST SAMPLE TRIANGLE) 10 IF (U.GT.0.0D+00) F = U**P GO TO 50 C SINGULARITY AT U=V=0, I.E. X=PA AND Y=PC*X+PD C (VERTEX OF 2ND SAMPLE TRIANGLE) 20 A = U + V IF (A.GT.0.0D+00) F = A**P GO TO 50 C SINGULARITY ALONG U+V=1 C (SIDE OF 3RD SAMPLE TRIANGLE) 30 A = 1.0D+00 - U - V IF (A.GE.DRELPR) F = A**P GO TO 50 C SINGULARITY ALONG U=0 AND V=0 C (TWO SIDES OF 4TH SAMPLE TRIANGLE) 40 A = U*V IF (A.GT.0.0D+00) F = A**P 50 F = F*V1/D RETURN END SUBROUTINE TRIEX(F, VER, EPSABS, EPSREL, LIMEV, RESULT, ABSERR, TRI 10 * IER) C ..................................................................... C 1. TRIEX C DOUBLE INTEGRATION OVER A TRIANGLE C STANDARD FORTRAN SUBROUTINE C C 2. PURPOSE C THE ROUTINE CALCULATES AN APPROXIMATION RESULT TO A C DOUBLE INTEGRAL I = INTEGRAL OF F OVER THE TRIANGLE C WITH VERTICES (VER(1,K),VER(2,K)) FOR K=1,2,3, C HOPEFULLY SATISFYING FOLLOWING CLAIM FOR ACCURACY C ABS(I-RESULT) .LE. MAX(EPSABS,EPSREL*ABS(I)). C C 3. CALLING SEQUENCE C CALL TRIEX (F,VER,EPSABS,EPSREL,LIMEV,RESULT,ABSERR,IER) C C PARAMETERS C F - DOUBLE PRECISION C FUNCTION DEFINING THE INTEGRAND F(X,Y) C THE ACTUAL NAME FOR F NEEDS TO BE DECLARED C E X T E R N A L IN THE DRIVER PROGRAM. C VER - DOUBLE PRECISION C ARRAY OF DIMENSION (2,3), WHERE VER(1,K) AND C VER(2,K) REPRESENT THE ABSCISSAE AND THE C ORDINATES RESPECTIVELY OF THE K-TH VERTEX C EPSABS - DOUBLE PRECISION C ABSOLUTE ACCURACY REQUESTED C EPSREL - DOUBLE PRECISION C RELATIVE ACCURACY REQUESTED C IF EPSABS .LE. 0 AND EPSREL .LE. 0, C THE ROUTINE WILL END WITH IER .NE. 0. C LIMEV - LIMIT ON THE NUMBER OF INTEGRAND EVALUATIONS C IF MORE THAN XXX EVALUATIONS ARE TO BE ALLOWED, C DIMENSION STATEMENTS AND DATA N IN TRIEX C (WHERE N IS THE MAXIMUM NUMBER OF TRIANGLES C ALLOWED) NEED TO BE ADJUSTED ACCORDING TO C XXX .LE. (N-1)*178/3 + 46, C I.E. PRESENT LIMIT N = 1000 IMPLIES C LIMEV .LE. 59320. C RESULT - DOUBLE PRECISION C APPROXIMATION TO THE INTEGRAL I C ABSERR - DOUBLE PRECISION C ESTIMATE OF THE ABSOLUTE ERROR ABS(I-RESULT) C IER - IER = 0 NORMAL AND RELIABLE TERMINATION OF C THE ROUTINE. IT IS ASSUMED THAT THE C REQUESTED ACCURACY HAS BEEN OBTAINED. C - IER .NE. 0 ABNORMAL TERMINATION OF THE ROUTINE C IT IS ASSUMED THAT THE REQUESTED C ACCURACY HAS NOT BEEN OBTAINED. C = 1 MAXIMUM NUMBER OF FUNCTION EVALUATIONS C ALLOWED HAS BEEN REACHED. C = 2 THE PRESENCE OF ROUNDOFF ERROR HAS BEEN C DETECTED, WHICH PREVENTS THE REQUESTED C ACCURACY FROM BEING ACHIEVED. C = 3 IT IS PRESUMED THAT THE TOLERANCE CANNOT C BE ACHIEVED WITHIN THE SPECIFIED NUMBER C OF INTEGRAND EVALUATIONS, AND THAT THE C RESULT RETURNED IS THE BEST WHICH CAN BE C OBTAINED WITHIN THAT LIMIT. C C 4. REMARKS C SUBROUTINES OR FUNCTIONS CALLED C - QUATRI C - DIVIN C - ORDER C - EPSALG C - D1MACH C ..................................................................... C DOUBLE PRECISION ABSERR, ABSER1, AREA, AREA4, BNDEC, DOVFLO, * DRELPR, DRES, DRESC, DRES0, DUNFLO, D1MACH, ELIST, EMACH, EPRN, * EPSABS, EPSREL, ERLARG, ERLAST, ERRBND, ERRMAX, ERRO4, ERRSUM, * F, OFRN, Q, RCOPY, RESEX, RESULT, RESUL1, RESLA, REX, RLAST, * RLIST, STOR, UFRN, VER, X, Y DOUBLE PRECISION DABS, DMAX1, DMIN1 INTEGER I, ICNT, ID, IER, IORD, J, J1, J2, J3, JPOINT, JROFF, * JSLAST, JSTAGE, JVRTX, K, KBIGER, LAST, LIMEV, LIMIT, L4, * MAXERR, MAXER0, N, NRMAX, NTR, NUMRL2, NVRTX, N2 INTEGER MAX0, MIN0 LOGICAL EXTRAP DIMENSION ELIST(1000), IORD(1000), JPOINT(4), JSTAGE(1000), * JVRTX(1000,3), RCOPY(52), REX(1000), RLIST(1000), VER(2,3), * X(3000), Y(3000) C COMMON /MACH/ EMACH, EPRN, UFRN, OFRN C C DIMENSIONS C ---------- C N IS THE MAXIMUM NUMBER OF TRIANGLES ALLOWED. C INCREASING N REQUIRES INCREASING DIMENSIONS IN C TRIEX (AND IN SUBROUTINE DIVIN CORRESPONDINGLY) C DIMENSIONS OF ELIST(N),IORD(N),JSTAGE(N),JVRTX(N,3), C REX(N),RLIST(N),X(3*N),Y(3*N) C SHOULD BE SUCH THAT C LIMEV .LE. (N-1)*178/3 + 46 EXTERNAL F DATA N /1000/ C C RCOPY IS OF DIMENSION AT LEAST (LIMEXP+2) C WITH DATA LIMEXP GIVEN IN SUBROUTINE EPSALG C C C MACHINE DEPENDENT CONSTANTS C --------------------------- C DRELPR IS THE LARGEST RELATIVE SPACING C DUNFLO IS THE SMALLEST POSITIVE MAGNITUDE C DOVFLO IS THE LARGEST MAGNITUDE C DRELPR = D1MACH(4) DUNFLO = D1MACH(1) DOVFLO = D1MACH(2) C C LIST OF MAJOR VARIABLES C ----------------------- C C X - LIST OF ABSCISSAE OF ALL VERTICES C Y - LIST OF ORDINATES OF ALL VERTICES C JVRTX - ARRAY OF NUMBERS OF ALL VERTICES C RLIST - LIST OF DEGREE-11 INTEGRAL APPROXIMATIONS FOR C ALL SUBTRIANGLES C REX - LIST OF DEGREE-9 INTEGRAL APPROXIMATIONS FOR C ALL SUBTRIANGLES C ELIST - LIST OF ERROR ESTIMATES FOR ALL SUBTRIANGLES C JSTAGE - LIST OF SUBDIVISION LEVELS FOR ALL C SUBTRIANGLES C AREA - CURRENT SUM OF DEGREE-11 RESULTS C ERRSUM - CURRENT SUM OF ERROR ESTIMATES C RESEX - APPROXIMATION TO THE TOTAL INTEGRAL COMPUTED C AS A SUM OF DEGREE-9 RESULTS OVER THE SMALLEST C SUBTRIANGLES (HIGHEST LEVEL OF SUBDIVISION) AND C OF DEGREE-11 RESULTS OVER THE LARGER TRIANGLES C ERRBND - ESTIMATE OF THE GLOBAL ABSOLUTE TOLERANCE C BNDEC - GLOBAL TOLERANCE OVER THE SET OF LARGER C TRIANGLES C ERLARG - SUM OF ERROR ESTIMATES OVER THE CURRENT SET OF C LARGE TRIANGLES C NTR - CURRENT NUMBER OF TRIANGLES IN THE PARTITION C OF THE ORIGINAL INTEGRATION REGION C MAXERR - POINTER TO SUBTRIANGLE WITH LARGEST ERROR C ESTIMATE C RCOPY - LIST OF THE ENTRIES TO THE EXTRAPOLATION C PROCEDURE C NUMRL2 - NUMBER OF ELEMENTS IN RCOPY C EXTRAP - LOGICAL VARIABLE DENOTING THAT THE ROUTINE C IS ATTEMPTING TO DECREASE THE VALUE OF ERLARG C THROUGH SUBDIVISION OF LARGE (LOW-LEVEL) C TRIANGLES, BEFORE A NEW EXTRAPOLATION IS C CARRIED OUT C C FIRST APPROXIMATION TO THE INTEGRAL C ----------------------------------- IER = 0 I = 0 K = 1 EPRN = 0.5D+01*DRELPR UFRN = DUNFLO/EPRN STOR = 0.0D+00 DO 10 J=1,3 X(J) = VER(1,J) Y(J) = VER(2,J) JVRTX(1,J) = J 10 CONTINUE CALL QUATRI(F, X(1), Y(1), X(2), Y(2), X(3), Y(3), RESEX, RESULT, * ABSERR, DRESC, I, K, STOR) DRES = DABS(RESEX) ERRBND = DMAX1(EPSABS,EPSREL*DRES,EPRN*DRES) IF (ABSERR.LE.ERRBND) GO TO 170 C C INITIALIZATION C -------------- C RLIST(1) = RESULT RCOPY(1) = RESEX ELIST(1) = ABSERR OFRN = DOVFLO JSTAGE(1) = 1 AREA = RESULT ERRSUM = ABSERR ERRMAX = ABSERR EMACH = DRELPR MAXERR = 1 NTR = 1 NRMAX = 1 NVRTX = 3 NUMRL2 = 2 N2 = 2 JROFF = 0 EXTRAP = .FALSE. C C LIMIT = NUMBER OF TRIANGLES ALLOWED IN FINAL PARTITION OF C ORIGINAL TRIANGLE C LIMIT = MAX0(4,MIN0(N,1+((LIMEV-46)/178)*3)) L4 = (LIMIT-8)/4 C C MAIN DO-LOOP C ------------ C DO 130 LAST=2,LIMIT,3 C C SUBDIVIDE THE TRIANGLE WITH NRMAX-TH LARGEST ERROR ESTIMATE C AND INTEGRATE OVER ITS SUBTRIANGLES C CALL DIVIN(MAXERR, X, Y, JVRTX, NTR, NVRTX, JPOINT) MAXER0 = MAXERR RLAST = RLIST(MAXERR) JSLAST = JSTAGE(MAXERR) ERLAST = ERRMAX AREA4 = 0.0D+00 ERRO4 = 0.0D+00 K = 1 ICNT = 0 STOR = 0.0D+00 DO 20 I=1,4 J = JPOINT(I) J1 = JVRTX(J,1) J2 = JVRTX(J,2) J3 = JVRTX(J,3) CALL QUATRI(F, X(J1), Y(J1), X(J2), Y(J2), X(J3), Y(J3), * REX(J), RLIST(J), ELIST(J), DRESC, I, K, STOR) AREA4 = AREA4 + RLIST(J) ERRO4 = ERRO4 + ELIST(J) JSTAGE(J) = JSLAST + 1 IF (ELIST(J).EQ.DRESC) ICNT = 1 20 CONTINUE C C IMPROVE PREVIOUS INTEGRAL APPROXIMATION AND ERROR C ESTIMATE AND TEST FOR ACCURACY C AREA = AREA + AREA4 - RLAST ERRSUM = ERRSUM + ERRO4 - ERRMAX IF (ERRSUM.LE.ERRBND) GO TO 150 C C TEST FOR ROUNDOFF ERROR AND EVENTUALLY SET ERROR FLAG C IF (ICNT.EQ.1) GO TO 30 Q = 0.0D+00 IF (AREA4.NE.0.0D+00) Q = RLAST/AREA4 IF (0.9999D+00.LT.Q .AND. Q.LT.1.0001D+00 .AND. * ERRMAX.LE.ERRO4) JROFF = JROFF + 1 IF (JROFF.GE.10) IER = 2 C C SET ERROR FLAG IN THE CASE THAT MAXIMUM NUMBER OF C SUBTRIANGLES IS REACHED C 30 IF (LAST.GE.LIMIT-2) IER = 1 IF (IER.NE.0) GO TO 140 C C CALL SUBROUTINE ORDER TO MAINTAIN THE DESCENDING C ORDERING IN THE LIST OF ERROR ESTIMATES, AND TO SELECT C THE SUBTRIANGLE WITH NRMAX-TH LARGEST ERROR ESTIMATE C (TO BE SUBDIVIDED NEXT) C CALL ORDER(LIMIT, LAST, MAXERR, ERRMAX, ELIST, IORD, NRMAX, NTR) IF (LAST.EQ.2) GO TO 120 C C CHECK WHETHER THE TRIANGLES OBTAINED IN THE FOREGOING C SUBDIVISION ARE LARGE, AND ADJUST ERLARG IF SO C ERLARG = ERLARG - ERLAST IF (JSTAGE(MAXER0).LE.NUMRL2) ERLARG = ERLARG + ERRO4 IF (EXTRAP) GO TO 40 C C CHECK WHETHER THE FOLLOWING SUBDIVISION WOULD PROVIDE C LARGE TRIANGLES, AND PROCEED IMMEDIATELY TO THIS C SUBDIVISION IF SO C IF (JSTAGE(MAXERR).LE.NUMRL2) GO TO 130 C C DECREASE THE INTEGRATION ERROR OVER THE COLLECTION OF C LARGE TRIANGLES C EXTRAP = .TRUE. NRMAX = 2 40 BNDEC = DMAX1(EPRN*DRES,DMIN1(0.1D+00*ERRBND,0.1D-02*DRES)) IF (ERLARG.LE.BNDEC) GO TO 60 ID = NRMAX DO 50 J=ID,NTR MAXERR = IORD(NRMAX) ERRMAX = ELIST(MAXERR) IF (JSTAGE(MAXERR).LE.NUMRL2) GO TO 130 NRMAX = NRMAX + 1 50 CONTINUE C C COMPUTE NEXT ENTRY FOR THE EXTRAPOLATION PROCEDURE C 60 RESEX = 0.0D+00 DO 70 J=1,NTR IF (JSTAGE(J).GT.NUMRL2) RESEX = RESEX + REX(J) IF (JSTAGE(J).LE.NUMRL2) RESEX = RESEX + RLIST(J) 70 CONTINUE DRES = DABS(RESEX) ERRBND = DMAX1(EPSABS,EPSREL*DRES,EPRN*DRES) C C EXTRAPOLATE AND TEST FOR ACCURACY C NUMRL2 = NUMRL2 + 1 N2 = N2 + 1 RCOPY(NUMRL2) = RESEX CALL EPSALG(NUMRL2, RCOPY, RESUL1, ABSER1, RESLA) ABSER1 = ABSER1 + ERLARG IF (NUMRL2.EQ.N2 .OR. N2.EQ.3) GO TO 80 IF (ABSER1.GT.ABSERR) GO TO 90 80 RESULT = RESUL1 ABSERR = ABSER1 DRES0 = DABS(RESULT) IF (ABSERR.LE.DMAX1(EPSABS,EPSREL*DRES0,EPRN*DRES0) .AND. * N2.GT.3) GO TO 170 C C IN THE CASE THAT THE EXTRAPOLATION ERROR ESTIMATE IS C CONSIDERABLY SMALLER THAN THE ERROR SUM OVER ALL C SUBTRIANGLES, CHECK IF ENOUGH STEPS ARE LEFT OVER C IN ORDER TO COMPLETE ANOTHER EXTRAPOLATION C STAGE. SET ERROR FLAG FOR ABNORMAL TERMINATION C IF NECESSARY C 90 IF (ABSERR*DABS(RESEX).GE.0.5D-01*ERRSUM*DRES0 .OR. LAST.LE.L4) * GO TO 110 KBIGER = 0 DO 100 J=1,NTR IF (ELIST(J).GT.BNDEC) KBIGER = KBIGER + 1 100 CONTINUE IF (KBIGER.GE.(LIMIT-LAST-2)/3) IER = 3 110 IF (IER.NE.0) GO TO 140 C C START NEW EXTRAPOLATION STAGE, BY INTEGRATING OVER C SUBTRIANGLE WITH LARGEST ERROR ESTIMATE C MAXERR = IORD(1) ERRMAX = ELIST(MAXERR) NRMAX = 1 ERLARG = ERRSUM GO TO 130 C C SECOND EXTRAPOLATION ENTRY C 120 ERLARG = ERRSUM RESEX = REX(1) + REX(2) + REX(3) + REX(4) RCOPY(2) = RESEX DRES = DABS(RESEX) RESLA = RESEX 130 CONTINUE C C IN THE CASE THAT IER .NE. 0, SELECT THE BEST APPROXIMATION C ---------------------------------------------------------- C (RESULT (EXTRAPOLATED), OR AREA (SUM OF INTEGRALS OVER C SUBTRIANGLES) ) C 140 DRES = DABS(AREA) IF (DRES*ABSERR.LE.DRES0*ERRSUM) GO TO 170 IF (DABS(AREA-RESEX)*DRES0.GT.DABS(RESULT-RESEX)*DRES) GO TO 170 150 RESULT = 0.0D+00 DO 160 J=1,NTR RESULT = RESULT + RLIST(J) 160 CONTINUE ABSERR = ERRSUM 170 RETURN END SUBROUTINE DIVIN(MAX, X, Y, IVRTX, NTR, NVRTX, IPOINT) DIV 10 C ..................................................................... C C 1. DIVIN C C SUBDIVISION OF A GIVEN TRIANGLE C STANDARD FORTRAN SUBROUTINE C C 2. PURPOSE C C THE ROUTINE SUBDIVIDES THE GIVEN TRIANGLE WITH INDEX MAX C INTO 4 SIMILAR SUBTRIANGLES, BY HALVING THE SIDES OF THE C GIVEN TRIANGLE. IT CALCULATES THE COORDINATES OF THE NEW C VERTICES, AND DEFINES THE SUBTRIANGLES BY DEFINING POINTERS C TO THESE TRIANGLES AND TO THEIR VERTICES. C C 3. CALLING SEQUENCE C C CALL DIVIN(MAX,X,Y,IVRTX,NTR,NVRTX,IPOINT) C C PARAMETERS C MAX - INTEGER C INDEX OF THE TRIANGLE WHICH MUST BE SUBDIVIDED C C X - DOUBLE PRECISION C ABSCISSAE OF THE VERTICES OF ALL TRIANGLES IN THE C PARTITION OF THE INTEGRATION DOMAIN C C Y - DOUBLE PRECISION C ORDINATES OF THE VERTICES OF ALL TRIANGLES IN THE C PARTITION OF THE INTEGRATION DOMAIN C C IVRTX - INTEGER C IVRTX(J,K) POINTS TO THE K-TH VERTEX (K=1,2,3) C OF THE TRIANGLE WITH INDEX J C C NTR - INTEGER C TOTAL NUMBER OF TRIANGLES IN THE PARTITION C NTR IS INCREASED BY 3 DURING THE COURSE OF THIS C SUBROUTINE C C NVRTX - INTEGER C TOTAL NUMBER OF VERTICES IN THE PARTITION C NVRTX IS INCREASED BY 3 DURING THIS CALL. C C IPOINT - INTEGER C POINTERS TO THE 4 NEW SUBTRIANGLES, SUCH THAT C IPOINT(1)=MAX, IPOINT(2)=NTR+1, IPOINT(3)=NTR+2, C IPOINT(4)=NTR+3 C C 4. REMARKS C C SUBROUTINE OR FUNCTIONS CALLED NONE C C ..................................................................... DOUBLE PRECISION X, Y INTEGER I, I1, I2, I3, I4, I5, I6, IP, IPOINT, IVRTX, MAX, NTR, * NVRTX C DIMENSION IPOINT(4), IVRTX(1000,3), X(3000), Y(3000) C C DIMENSIONS C ---------- C FOR ALL DIMENSION SPECIFICATIONS SEE CALLING ROUTINE TRIEX C C LIST OF MAJOR VARIABLES C ----------------------- C I1,...,I6 - INDICES OF VERTICES OF THE 4 NEW SUBTRIANGLES C C C SET POINTERS TO THE 4 NEW SUBTRIANGLES C DO 10 I=1,3 IP = NTR + I IPOINT(I+1) = IP 10 CONTINUE IPOINT(1) = MAX C C DEFINE VERTICES OF THE SUBTRIANGLES C I1 = IVRTX(MAX,1) I2 = IVRTX(MAX,2) I3 = IVRTX(MAX,3) I4 = NVRTX + 1 I5 = NVRTX + 2 I6 = NVRTX + 3 X(I4) = (X(I2)+X(I3))*0.5D+00 Y(I4) = (Y(I2)+Y(I3))*0.5D+00 X(I5) = (X(I3)+X(I1))*0.5D+00 Y(I5) = (Y(I3)+Y(I1))*0.5D+00 X(I6) = (X(I1)+X(I2))*0.5D+00 Y(I6) = (Y(I1)+Y(I2))*0.5D+00 C C ADJUST NVRTX C NVRTX = I6 C C THE SUBDIVISION IS ORGANIZED AS IS INDICATED IN THE C PICTURES BELOW. THE POINTERS TO THE VERTICES OF THE C SUBTRIANGLES ARE SET IN SUCH A WAY THAT C SUBDIVISION OF THE TRIANGLE WITH INDEX MAX PRODUCES THE C TRIANGLES WITH INDICES MAX, NTR+1, NTR+2, NTR+3 C (DETERMINED BY (I1,I6,I5), (I6,I2,I4), (I5,I4,I3) AND C (I4,I6,I5) RESPECTIVELY). C C I3 C C C I5 I4 C C C I1 I6 I2 C IVRTX(MAX,1) = I1 IVRTX(MAX,2) = I6 IVRTX(MAX,3) = I5 NTR = NTR + 1 IVRTX(NTR,1) = I6 IVRTX(NTR,2) = I2 IVRTX(NTR,3) = I4 NTR = NTR + 1 IVRTX(NTR,1) = I5 IVRTX(NTR,2) = I4 IVRTX(NTR,3) = I3 NTR = NTR + 1 IVRTX(NTR,1) = I4 IVRTX(NTR,2) = I6 IVRTX(NTR,3) = I5 RETURN END SUBROUTINE QUATRI(F, U1, V1, U2, V2, U3, V3, RES9, RES11, EST, QUA 10 * DRESC, I, K, STOR) C ..................................................................... C C 1. QUATRI C C INTEGRATION RULES C STANDARD FORTRAN SUBROUTINE C C 2. PURPOSE C C TO COMPUTE - IF = INTEGRAL OF F OVER THE TRIANGLE C WITH VERTICES (U1,V1),(U2,V2),(U3,V3), AND ERROR C ESTIMATE, C - INTEGRAL OF ABS(F) OVER THIS TRIANGLE C C 3. CALLING SEQUENCE C C CALL QUATRI(F,U1,V1,U2,V2,U3,V3,RES9,RES11,EST,DRESC, C I,K,STOR) C C PARAMETERS C F - DOUBLE PRECISION C FUNCTION SUBPROGRAM DEFINING THE INTEGRAND C F(X,Y) THE ACTUAL NAME FOR F NEEDS TO BE C DECLARED E X T E R N A L IN THE CALLING C PROGRAM. C C U1,U2,U3- DOUBLE PRECISION C ABSCISSAE OF VERTICES C C V1,V2,V3- DOUBLE PRECISION C ORDINATES OF VERTICES C C RES9 - DOUBLE PRECISION C APPROXIMATION TO THE INTEGRAL IF, OBTAINED BY THE C LYNESS AND JESPERSEN RULE OF DEGREE 9, USING C 19 POINTS C C RES11 - DOUBLE PRECISION C APPROXIMATION TO THE INTEGRAL IF, OBTAINED BY THE C LYNESS AND JESPERSEN RULE OF DEGREE 11, C USING 28 POINTS C C EST - DOUBLE PRECISION C ESTIMATE OF THE ABSOLUTE ERROR C C DRESC - DOUBLE PRECISION C APPROXIMATION TO THE INTEGRAL OF ABS(F- IF/DJ), C OBTAINED BY THE RULE OF DEGREE 9, AND USED FOR C THE COMPUTATION OF EST C C I - INTEGER C NUMBER OF TRIANGLES IN CURRENT SUBDIVISION C (1 .LE. I .LE. 4) C C K - INTEGER C INDEX OF POINTS TO BE SHARED WITH 4TH TRIANGLE C OF SUBDIVISION C C STOR - DOUBLE PRECISION C SUM OF FUNCTION VALUES AT THESE POINTS C C 4. REMARKS C C SUBROUTINES OR FUNCTIONS CALLED C - F (USER-SUPPLIED INTEGRAND FUNCTION) C C ..................................................................... C DOUBLE PRECISION DJ, DF0, DRESC, EMACH, EPRN, EST, F, FSTOR, FV, * F0, OFRN, RESAB9, RES9, RES11, STOR, U1, U2, U3, UFRN, V1, V2, * V3, W, W90, W110, X, Y, ZETA1, ZETA2, Z1, Z2, Z3 DOUBLE PRECISION DMAX1, DMIN1 INTEGER I, J, JEND, K, KOUNT, L C DIMENSION FSTOR(3), FV(19), W(15), X(3), Y(3), ZETA1(15), * ZETA2(15) C COMMON /MACH/ EMACH, EPRN, UFRN, OFRN C C FIRST HOMOGENEOUS COORDINATES OF POINTS IN DEGREE-9 C AND DEGREE-11 FORMULA, TAKEN WITH MULTIPLICITY 3 DATA ZETA1(1), ZETA1(2), ZETA1(3), ZETA1(4), ZETA1(5), ZETA1(6), * ZETA1(7), ZETA1(8), ZETA1(9), ZETA1(10), ZETA1(11), ZETA1(12), * ZETA1(13), ZETA1(14), ZETA1(15) /0.2063496160252593D-01, * 0.1258208170141290D+00,0.6235929287619356D+00, * 0.9105409732110941D+00,0.3683841205473626D-01, * 0.7411985987844980D+00,0.9480217181434233D+00, * 0.8114249947041546D+00,0.1072644996557060D-01, * 0.5853132347709715D+00,0.1221843885990187D+00, * 0.4484167758913055D-01,0.6779376548825902D+00,0.0D+00, * 0.8588702812826364D+00/ C SECOND HOMOGENEOUS COORDINATES OF POINTS IN DEGREE-9 C AND DEGREE-11 FORMULA, TAKEN WITH MULTIPLICITY 3 DATA ZETA2(1), ZETA2(2), ZETA2(3), ZETA2(4), ZETA2(5), ZETA2(6), * ZETA2(7), ZETA2(8), ZETA2(9), ZETA2(10), ZETA2(11), ZETA2(12), * ZETA2(13), ZETA2(14), ZETA2(15) /0.4896825191987370D+00, * 0.4370895914929355D+00,0.1882035356190322D+00, * 0.4472951339445297D-01,0.7411985987844980D+00, * 0.3683841205473626D-01,0.2598914092828833D-01, * 0.9428750264792270D-01,0.4946367750172147D+00, * 0.2073433826145142D+00,0.4389078057004907D+00, * 0.6779376548825902D+00,0.4484167758913055D-01, * 0.8588702812826364D+00,0.0D+00/ C WEIGHTS OF MIDPOINT OF TRIANGLE IN DEGREE-9 C RESP. DEGREE-11 FORMULAE DATA W90 /0.9713579628279610D-01/ DATA W110 /0.8797730116222190D-01/ C WEIGHTS IN DEGREE-9 AND DEGREE-11 RULE DATA W(1), W(2), W(3), W(4), W(5), W(6), W(7), W(8), W(9), W(10), * W(11), W(12), W(13), W(14), W(15) /0.3133470022713983D-01, * 0.7782754100477543D-01,0.7964773892720910D-01, * 0.2557767565869810D-01,0.4328353937728940D-01, * 0.4328353937728940D-01,0.8744311553736190D-02, * 0.3808157199393533D-01,0.1885544805613125D-01, * 0.7215969754474100D-01,0.6932913870553720D-01, * 0.4105631542928860D-01,0.4105631542928860D-01, * 0.7362383783300573D-02,0.7362383783300573D-02/ C C LIST OF MAJOR VARIABLES C ---------------------- C DJ - AREA OF THE TRIANGLE C DRESC - APPROXIMATION TO INTEGRAL OF C ABS(F- IF/DJ) OVER THE TRIANGLE C RESAB9 - APPROXIMATION TO INTEGRAL OF C ABS(F) OVER THE TRIANGLE C X - CARTESIAN ABSCISSAE OF THE INTEGRATION C POINTS C Y - CARTESIAN ORDINATES OF THE INTEGRATION C POINTS C FV - FUNCTION VALUES C FSTOR C C COMPUTE DEGREE-9 AND DEGREE-11 RESULTS FOR IF/DJ AND C DEGREE-9 APPROXIMATION FOR ABS(F) C DJ = DABS(U1*V2-U2*V1-U1*V3+V1*U3+U2*V3-V2*U3)*0.5D+00 F0 = F((U1+U2+U3)/3.0D+00,(V1+V2+V3)/3.0D+00) RES9 = F0*W90 RESAB9 = DABS(F0)*W90 FV(1) = F0 KOUNT = 1 RES11 = F0*W110 IF (I.EQ.4) RES11 = RES11 + STOR*W(15) JEND = 15 IF (I.EQ.4) JEND = 13 DO 50 J=1,JEND Z1 = ZETA1(J) Z2 = ZETA2(J) Z3 = 1.0D+00 - Z1 - Z2 X(1) = Z1*U1 + Z2*U2 + Z3*U3 Y(1) = Z1*V1 + Z2*V2 + Z3*V3 X(2) = Z2*U1 + Z3*U2 + Z1*U3 Y(2) = Z2*V1 + Z3*V2 + Z1*V3 X(3) = Z3*U1 + Z1*U2 + Z2*U3 Y(3) = Z3*V1 + Z1*V2 + Z2*V3 IF (J.GT.6) GO TO 20 F0 = 0.0D+00 DF0 = 0.0D+00 DO 10 L=1,3 KOUNT = KOUNT + 1 FV(KOUNT) = F(X(L),Y(L)) F0 = F0 + FV(KOUNT) DF0 = DF0 + DABS(FV(KOUNT)) 10 CONTINUE RES9 = RES9 + F0*W(J) RESAB9 = RESAB9 + DF0*W(J) GO TO 50 20 F0 = 0.0D+00 DO 30 L=1,3 FSTOR(L) = F(X(L),Y(L)) F0 = F0 + FSTOR(L) 30 CONTINUE IF (J.LT.14) GO TO 40 STOR = STOR + FSTOR(K) K = K + 1 IF (K.GT.3) K = 1 40 RES11 = RES11 + F0*W(J) 50 CONTINUE C C COMPUTE DEGREE-9 APPROXIMATION FOR THE INTEGRAL OF C ABS(F-IF/DJ) C DRESC = DABS(FV(1)-RES9)*W90 KOUNT = 2 DO 60 J=1,6 DRESC = DRESC + (DABS(FV(KOUNT)-RES9)+DABS(FV(KOUNT+1)-RES9) * +DABS(FV(KOUNT+2)-RES9))*W(J) KOUNT = KOUNT + 3 60 CONTINUE C C COMPUTE DEGREE-9 AND DEGREE-11 APPROXIMATIONS FOR IF, C AND ERROR ESTIMATE C RES9 = RES9*DJ RES11 = RES11*DJ RESAB9 = RESAB9*DJ DRESC = DRESC*DJ EST = DABS(RES9-RES11) IF (DRESC.NE.0.0D+00) EST = DMAX1(EST,DRESC*DMIN1(1.0D+00, * (20.0D+00*EST/DRESC)**1.5D+00)) IF (RESAB9.GT.UFRN) EST = DMAX1(EPRN*RESAB9,EST) RETURN END SUBROUTINE EPSALG(N, EPSTAB, RESULT, ABSERR, RESLA) EPS 10 C C ................................................................ C C 1. EPSALG C C EPSILON ALGORITHM C STANDARD FORTRAN SUBROUTINE C C 2. PURPOSE C C THE ROUTINE TRANSFORMS A GIVEN SEQUENCE USING THE C EPSILON ALGORITHM OF P. WYNN. C AN ESTIMATE OF THE ABSOLUTE ERROR IS ALSO GIVEN. C THE CONDENSED EPSILON TABLE IS COMPUTED. ONLY THOSE C ELEMENTS NEEDED FOR THE COMPUTATION OF THE NEXT DIAGONAL C ARE PRESERVED. C C 3. CALLING SEQUENCE C C CALL EPSALG (N,EPSTAB,RESULT,ABSERR,RES3LA,NRES) C C PARAMETERS C N - INTEGER C EPSTAB(N) CONTAINS THE NEW ELEMENT IN THE C FIRST COLUMN OF THE EPSILON TABLE. C C EPSTAB - DOUBLE PRECISION C ONE DIMENSIONAL ARRAY CONTAINING THE C ELEMENTS OF THE LAST TWO DIAGONALS OF C THE TRIANGULAR EPSILON TABLE C THE ELEMENTS ARE NUMBERED STARTING AT THE C RIGHT-HAND CORNER OF THE TRIANGLE. C THE DIMENSION SHOULD BE AT LEAST (LIMEXP+2) C (SEE DATA LIMEXP). C C RESULT - DOUBLE PRECISION C RESULTING APPROXIMATION TO THE INTEGRAL C C ABSERR - DOUBLE PRECISION C ESTIMATE OF THE ABSOLUTE ERROR C C RESLA - DOUBLE PRECISION C PREVIOUS RESULT C C 4. REMARKS C C SUBROUTINES OR FUNCTIONS CALLED NONE C C .................................................................. C DOUBLE PRECISION ABSERR, DELTA1, DELTA2, DELTA3, EMACH, EPRN, * EPSINF, EPSTAB, ERROR, ERR1, ERR2, ERR3, E0, E1, E2, E3, E1ABS, * OFRN, RES, RESULT, RESLA, SS, TOL1, TOL2, TOL3, UFRN DOUBLE PRECISION DABS, DMAX1 INTEGER I, IB, IB2, IE, IN, K1, K2, K3, LIMEXP, N, NEWELM, NUM DIMENSION EPSTAB(52) C COMMON /MACH/ EMACH, EPRN, UFRN, OFRN C C LIMEXP IS THE MAXIMUM NUMBER OF ELEMENTS THE EPSILON C TABLE CAN CONTAIN. IF THIS NUMBER IS REACHED, THE UPPER C DIAGONAL OF THE EPSILON TABLE IS DELETED. C ( EPSTAB IS OF DIMENSION (LIMEXP+2) AT LEAST.) C DATA LIMEXP /50/ C C LIST OF MAJOR VARIABLES C ----------------------- C E0 - THE 4 ELEMENTS ON WHICH THE C E1 COMPUTATION OF A NEW ELEMENT IN C E2 THE EPSILON TABLE IS BASED C E3 E0 C E3 E1 NEW C E2 C NEWELM - NUMBER OF ELEMENTS TO BE COMPUTED IN THE NEW C DIAGONAL C ERROR - ERROR = ABS(E0-E1)+ABS(E1-E2)+ABS(E2-NEW) C RESULT - THE ELEMENT IN THE NEW DIAGONAL WITH LEAST C ERROR C EPSTAB(N+2) = EPSTAB(N) NEWELM = (N-1)/2 EPSTAB(N) = OFRN ABSERR = OFRN NUM = N K1 = N DO 40 I=1,NEWELM K2 = K1 - 1 K3 = K1 - 2 RES = EPSTAB(K1+2) E0 = EPSTAB(K3) E1 = EPSTAB(K2) E2 = RES E1ABS = DABS(E1) DELTA2 = E2 - E1 ERR2 = DABS(DELTA2) TOL2 = DMAX1(DABS(E2),E1ABS)*EMACH DELTA3 = E1 - E0 ERR3 = DABS(DELTA3) TOL3 = DMAX1(E1ABS,DABS(E0))*EMACH IF (ERR2.GT.TOL2 .OR. ERR3.GT.TOL3) GO TO 10 C C IF E0, E1 AND E2 ARE EQUAL TO WITHIN MACHINE C ACCURACY, CONVERGENCE IS ASSUMED. C RESULT = E2 C ABSERR = ABS(E1-E0)+ABS(E2-E1) C RESULT = RES ABSERR = ERR2 + ERR3 GO TO 90 10 E3 = EPSTAB(K1) EPSTAB(K1) = E1 DELTA1 = E1 - E3 ERR1 = DABS(DELTA1) TOL1 = DMAX1(E1ABS,DABS(E3))*EMACH C C IF TWO ELEMENTS ARE VERY CLOSE TO EACH OTHER, OMIT C A PART OF THE TABLE BY ADJUSTING THE VALUE OF N. C IF (ERR1.LE.TOL1 .OR. ERR2.LE.TOL2 .OR. ERR3.LE.TOL3) GO TO 20 SS = 0.1D+01/DELTA1 + 0.1D+01/DELTA2 - 0.1D+01/DELTA3 EPSINF = DABS(SS*E1) C C TEST TO DETECT IRREGULAR BEHAVIOUR IN THE TABLE, AND C EVENTUALLY OMIT A PART OF THE TABLE ADJUSTING THE VALUE C OF N. C IF (EPSINF.GT.0.1D-03) GO TO 30 20 N = I + I - 1 GO TO 50 C C COMPUTE A NEW ELEMENT AND EVENTUALLY ADJUST C THE VALUE OF RESULT C 30 RES = E1 + 0.1D+01/SS EPSTAB(K1) = RES K1 = K1 - 2 ERROR = ERR2 + DABS(RES-E2) + ERR3 IF (ERROR.GT.ABSERR) GO TO 40 ABSERR = ERROR RESULT = RES 40 CONTINUE C C SHIFT THE TABLE C 50 IF (N.EQ.LIMEXP) N = 2*(LIMEXP/2) - 1 IB = 1 IF ((NUM/2)*2.EQ.NUM) IB = 2 IE = NEWELM + 1 DO 60 I=1,IE IB2 = IB + 2 EPSTAB(IB) = EPSTAB(IB2) IB = IB2 60 CONTINUE IF (NUM.EQ.N) GO TO 80 IN = NUM - N + 1 DO 70 I=1,N EPSTAB(I) = EPSTAB(IN) IN = IN + 1 70 CONTINUE C C COMPUTE ERROR ESTIMATE C 80 ABSERR = DABS(RESULT-RESLA) RESLA = RESULT 90 ABSERR = DMAX1(ABSERR,EPRN*DABS(RESULT)) RETURN END SUBROUTINE ORDER(LIMIT, LAST, MAXERR, ERMAX, ELIST, IORD, NRMAX, ORD 10 * NAVAIL) C C ..................................................................... C C 1. ORDER C C SORTING ROUTINE C STANDARD FORTRAN SUBROUTINE C C 2. PURPOSE C C THIS ROUTINE MAINTAINS THE DESCENDING ORDERING C IN THE LIST OF THE ERROR ESTIMATES RESULTING FROM THE C SUBDIVISION PROCESS. AT EACH CALL 4 ERROR ESTIMATES C ARE INSERTED USING THE SEQUENTIAL SEARCH METHOD C TOP-DOWN FOR THE LARGEST ESTIMATE, BOTTOM UP C FOR THE OTHER ONES. C C 3. CALLING SEQUENCE C C CALL ORDER(LIMIT,LAST,MAXERR,ERMAX,ELIST,IORD,NRMAX, C NAVAIL) C C PARAMETERS C LIMIT - INTEGER C MAXIMUM NUMBER OF ERROR ESTIMATES C THE LIST CAN CONTAIN C C LAST - INTEGER C NUMBER OF ERROR ESTIMATES CURRENTLY IN C THE LIST ELIST(LAST) CONTAINS THE SMALLEST C ERROR ESTIMATE C C MAXERR - INTEGER C POINTS TO THE NRMAX-TH LARGEST ERROR C ESTIMATE CURRENTLY IN THE LIST C C ERMAX - DOUBLE PRECISION C NRMAX-TH LARGEST ERROR ESTIMATE C ERMAX = ELIST(MAXERR) C C ELIST - DOUBLE PRECISION C ARRAY OF DIMENSION LIMIT CONTAINING C THE ERROR ESTIMATES C C IORD - INTEGER C ARRAY CONTAINING POINTERS TO ELIST SO C THAT IORD(1) POINTS TO THE LARGEST C ERROR ESTIMATE,..., IORD(LAST) TO THE C SMALLEST ERROR ESTIMATE C C NRMAX - INTEGER C MAXERR = IORD(NRMAX) C C NAVAIL - INTEGER C NAVAIL = LAST+2 C C 4. REMARKS C C SUBROUTINES OR FUNCTIONS CALLED NONE C C ..................................................................... DOUBLE PRECISION ELIST, ERMAX, ERRMAX, ERRMIN INTEGER I, IBEG, II, INCR, IND, INDI, INDJ, IORD, IORD4, ISUCC, * J, JBN, K, LAST, LIMIT, MAXERR, MINERR, NAVAIL, NRMAX, NTEST, NUM DIMENSION ELIST(LIMIT), IORD(LIMIT), IORD4(4) C C SORT THE 4 NEW ELIST-ENTRIES INTO DESCENDING ORDERING C IORD4(1) = MAXERR IORD4(2) = LAST IORD4(3) = LAST + 1 IORD4(4) = LAST + 2 DO 20 K=1,3 I = K DO 10 II=1,3 J = I + 1 INDI = IORD4(I) INDJ = IORD4(J) IF (ELIST(INDI).GE.ELIST(INDJ)) GO TO 20 IORD4(I) = INDJ IORD4(J) = INDI I = I - 1 IF (I.EQ.0) GO TO 20 10 CONTINUE 20 CONTINUE C C IF THE LIST CONTAINS ONLY 4 ELEMENTS, RETURN C IF (LAST.GT.2) GO TO 40 DO 30 I=1,4 IORD(I) = IORD4(I) 30 CONTINUE GO TO 170 C C THIS PART OF THE ROUTINE IS ONLY EXECUTED IF, DUE TO A C DIFFICULT INTEGRAND, SUBDIVISION INCREASED THE ERROR C ESTIMATE. IN THE NORMAL CASE THE INSERT PROCEDURE C SHOULD START AFTER THE NRMAX-TH LARGEST ERROR ESTIMATE. C 40 MAXERR = IORD4(1) ERRMAX = ELIST(MAXERR) IF (NRMAX.EQ.1) GO TO 60 II = NRMAX - 1 DO 50 I=1,II ISUCC = IORD(NRMAX-1) IF (ERRMAX.LE.ELIST(ISUCC)) GO TO 60 IORD(NRMAX) = ISUCC NRMAX = NRMAX - 1 50 CONTINUE C C INSERT ERRMAX BY TRAVERSING THE LIST TOP-DOWN C STARTING COMPARISON FROM THE ELEMENT ELIST (IORD(NRMAX+1)) C 60 NAVAIL = LAST + 2 JBN = NAVAIL - 3 IBEG = NRMAX + 1 IF (IBEG.GT.JBN) GO TO 80 DO 70 I=IBEG,JBN ISUCC = IORD(I) IF (ERRMAX.GE.ELIST(ISUCC)) GO TO 100 IORD(I-1) = ISUCC 70 CONTINUE 80 IND = JBN DO 90 I=1,4 IORD(IND) = IORD4(I) IND = IND + 1 90 CONTINUE GO TO 170 C C INSERT THE OTHER ESTIMATES BY TRAVERSING THE LIST BOTTOM-UP C 100 IORD(I-1) = MAXERR K = JBN NUM = JBN - I + 1 INCR = 3 DO 160 II=1,3 NTEST = NUM MINERR = IORD4(INCR+1) ERRMIN = ELIST(MINERR) DO 110 J=1,NTEST ISUCC = IORD(K) IF (ERRMIN.LT.ELIST(ISUCC)) GO TO 150 IND = K + INCR IORD(IND) = ISUCC K = K - 1 NUM = NUM - 1 110 CONTINUE GO TO (140, 130, 120), INCR 120 IORD(I+2) = IORD4(4) 130 IORD(I+1) = IORD4(3) 140 IORD(I) = IORD4(2) GO TO 170 150 IND = K + INCR IORD(IND) = MINERR INCR = INCR - 1 160 CONTINUE C C SET MAXERR AND ERMAX C 170 MAXERR = IORD(NRMAX) ERMAX = ELIST(MAXERR) RETURN END C.................................................................. MAN 10 C MAN 20 C INTEGRATION OVER A TRIANGLE TEST TRIEX MAN 30 C MAN 40 C THE FORM IN WHICH THE TEST INTEGRALS ARE SUPPLIED TO THE MAN 50 C INTEGRATOR CORRESPONDS TO MAN 60 C I(PA,PB) I(PC*X+PD,PG*X+PH) F(X,Y) DXDY MAN 70 C WITH MAN 80 C F(X,Y) = G(U(X),V(X,Y))*V1(X)/D MAN 90 C WHERE MAN 100 C D = PB-PA MAN 110 C V1(X) = (1-(X-PA)/D)/((PG-PC)*X+PH-PD) MAN 120 C U(X) = (X-PA)/D MAN 130 C V(X,Y) = V1(X)*(Y-PC*X-PD) MAN 140 C HENCE THIS IS THE TRANSFORMED VERSION OF THE INTEGRAL OVER MAN 150 C THE UNIT TRIANGLE GIVEN BY MAN 160 C I(0,1) I(0,1-X) G(X,Y) DXDY MAN 170 C THE FOLLOWING LIST PROVIDES, FOR EACH TEST INTEGRAL, THE MAN 180 C PARAMETERS PA, PB, PC, PD, PG, PH, AND THE VERTICES OF THE MAN 190 C TRIANGLE OVER WHICH IS THUS INTEGRATED: VER(1,K) = ABSCISSAE MAN 200 C AND VER(2,K) = ORDINATES OF THE VERTICES INDEXED K. MAN 210 C MAN 220 C INTEGRAL NUMBER 1 MAN 230 C G(X,Y) = X**(-0.2) MAN 240 C PA=1 PB=5 PC=0 PD=0 PG=0.25 PH=-0.25 MAN 250 C VER(1,1)=1 VER(1,2)=5 VER(1,3)=5 MAN 260 C VER(2,1)=0 VER(2,2)=0 VER(2,3)=1 MAN 270 C MAN 280 C INTEGRAL NUMBER 2 MAN 290 C G(X,Y) = (X+Y)**(-0.2) MAN 300 C PA=0 PB=1 PC=0 PD=0 PG=-1 PH=1 MAN 310 C VER(1,1)=0 VER(1,2)=1 VER(1,3)=0 MAN 320 C VER(2,1)=0 VER(2,2)=0 VER(2,3)=1 MAN 330 C MAN 340 C INTEGRAL NUMBER 3 MAN 350 C G(X,Y) = (1-X-Y)**(-0.2) MAN 360 C PA=-1 PB=3 PC=0.25 PD=-2.75 PD=-1 PH=1 MAN 370 C VER(1,1)=-1 VER(1,2)=3 VER(1,3)=-1 MAN 380 C VER(2,1)=-3 VER(2,2)=-2 VER(2,3)=2 MAN 390 C MAN 400 C INTEGRAL NUMBER 4 MAN 410 C G(X,Y) = (X*Y)**(-0.2) MAN 420 C PA=0 PB=-7 PC=0 PD=0 PG=-3/7 PH=-3 MAN 430 C VER(1,1)=0 VER(1,2)=-7 VER(1,3)=0 MAN 440 C VER(2,1)=0 VER(2,2)=0 VER(2,3)=-3 MAN 450 C MAN 460 C.................................................................. MAN 470 REAL ABSERR, C, D, EPSABS, EPSREL, ER, ES, EXA, EXACT, F, G, H, MAN 480 * P, PA, PB, PC, PD, PG, PH, Q, R, RESULT, VER MAN 490 INTEGER IER, INDACC, INUM, KOUNT, LIMEV, NOUT, NUM, N2, N2M1 MAN 500 REAL ABS MAN 510 DIMENSION C(4), D(4), EXACT(4), G(4), H(4), P(8), Q(8), R(8), MAN 520 * VER(2,3) MAN 530 COMMON /CNT/ KOUNT, NUM MAN 540 COMMON /PAR/ PA, PB, PC, PD, PG, PH MAN 550 EXTERNAL F MAN 560 C MAN 570 C EXACT VALUES MAN 580 C MAN 590 DATA EXACT(1), EXACT(2), EXACT(3), EXACT(4) /0.6944444444444444E+0MAN 600 * 0,0.5555555555555556E+00,0.6944444444444444E+00, MAN 610 * 0.9481026454955768E+00/ MAN 620 C MAN 630 C (P(2*NUM-1),P(2*NUM)) IS 1ST VERTEX FOR INTEGRAL NUM=1,2,3,4 MAN 640 C MAN 650 DATA P(1), P(2), P(3), P(4), P(5), P(6), P(7), P(8) MAN 660 * /1.0E+00,0.0E+00,0.0E+00,0.0E+00,-1.0E+00,-3.0E+00,0.0E+00, MAN 670 * 0.0E+00/ MAN 680 C MAN 690 C (Q(2*NUM-1),Q(2*NUM)) IS 2ND VERTEX FOR INTEGRAL NUM=1,2,3,4 MAN 700 C MAN 710 DATA Q(1), Q(2), Q(3), Q(4), Q(5), Q(6), Q(7), Q(8) MAN 720 * /5.0E+00,0.0E+00,1.0E+00,0.0E+00,3.0E+00,-2.0E+00,-7.0E+00, MAN 730 * 0.0E+00/ MAN 740 C MAN 750 C (R(2*NUM-1),R(2*NUM)) IS 3RD VERTEX FOR INTEGRAL NUM=1,2,3,4 MAN 760 C MAN 770 DATA R(1), R(2), R(3), R(4), R(5), R(6), R(7), R(8) MAN 780 * /5.0E+00,1.0E+00,0.0E+00,1.0E+00,-1.0E+00,2.0E+00,0.0E+00, MAN 790 * -3.0E+00/ MAN 800 C MAN 810 C PARAMETERS MAN 820 C MAN 830 DATA C(1), C(2), C(3), C(4) /0.0E+00,0.0E+00,0.25E+00,0.0E+00/ MAN 840 DATA D(1), D(2), D(3), D(4) /0.0E+00,0.0E+00,-2.75E+00,0.0E+00/ MAN 850 DATA G(1), G(2), G(3), G(4) /0.25E+00,-1.0E+00,-1.0E+00, MAN 860 * -0.428571428571428571E+00/ MAN 870 DATA H(1), H(2), H(3), H(4) /-0.25E+00,1.0E+00,1.0E+00,-3.0E+00/ MAN 880 C MAN 890 C OUTPUT UNIT MAN 900 C MAN 910 DATA NOUT /21/ MAN 920 C MAN 930 C HEADINGS MAN 940 C MAN 950 WRITE (NOUT,99999) MAN 960 99999 FORMAT (1X, 27HINTEGRATION OVER A TRIANGLE/1X, 27(1H.)/) MAN 970 WRITE (NOUT,99998) MAN 980 99998 FORMAT (6X, 8HRELATIVE, 9X, 8HCOMPUTED, 10X, 8HRELATIVE, 2X, MAN 990 * 6HESTIM., 3X, 6HNUMBER/1X, 3HNR., 2X, 8HACCURACY, 9X, 8HINTEGRAL,MAN 1000 * 10X, 5HERROR, 5X, 8HRELATIVE, 1X, 6HOF F , 4H IER/6X, 7HREQUEST,MAN 1010 * 2HED, 36X, 5HERROR, 4X, 5HEVAL.) MAN 1020 C MAN 1030 C INTEGRATION MAN 1040 C MAN 1050 EPSABS = 0.0E+00 MAN 1060 LIMEV = 50000 MAN 1070 DO 40 INUM=1,4 MAN 1080 EXA = EXACT(INUM) MAN 1090 EPSREL = 1.0E-01 MAN 1100 DO 30 INDACC=1,3 MAN 1110 KOUNT = 0 MAN 1120 NUM = INUM MAN 1130 N2 = NUM + NUM MAN 1140 N2M1 = N2 - 1 MAN 1150 VER(1,1) = P(N2M1) MAN 1160 VER(2,1) = P(N2) MAN 1170 VER(1,2) = Q(N2M1) MAN 1180 VER(2,2) = Q(N2) MAN 1190 VER(1,3) = R(N2M1) MAN 1200 VER(2,3) = R(N2) MAN 1210 PA = P(N2M1) MAN 1220 PB = Q(N2M1) MAN 1230 PC = C(NUM) MAN 1240 PD = D(NUM) MAN 1250 PG = G(NUM) MAN 1260 PH = H(NUM) MAN 1270 CALL TRIEX(F, VER, EPSABS, EPSREL, LIMEV, RESULT, ABSERR, IER)MAN 1280 C ....................................................... MAN 1290 ER = ABS((RESULT-EXA)/EXA) MAN 1300 ES = ABSERR/ABS(RESULT) MAN 1310 IF (INDACC.GT.1) GO TO 10 MAN 1320 WRITE (NOUT,99997) NUM, EPSREL, RESULT, ER, ES, KOUNT, IER MAN 1330 99997 FORMAT (/1X, I2, E11.2, 6X, E14.7, 5X, 2E10.2, I7, I4) MAN 1340 GO TO 20 MAN 1350 10 WRITE (NOUT,99996) EPSREL, RESULT, ER, ES, KOUNT, IER MAN 1360 99996 FORMAT (3X, E11.2, 6X, E14.7, 5X, 2E10.2, I7, I4) MAN 1370 20 EPSREL = EPSREL*1.0E-02 MAN 1380 30 CONTINUE MAN 1390 WRITE (NOUT,99995) EXA MAN 1400 99995 FORMAT (7X, 7HEXACT =, 6X, E14.7) MAN 1410 40 CONTINUE MAN 1420 STOP MAN 1430 END MAN 1440 C F 10 C TEST INTEGRANDS F 20 C F 30 REAL FUNCTION F(X, Y) F 40 REAL A, D, SRELPR, R1MACH, P, PA, PB, PC, PD, PG, PH, X, Y, U, V, * V1 INTEGER KOUNT, NUM COMMON /CNT/ KOUNT, NUM COMMON /PAR/ PA, PB, PC, PD, PG, PH C C SRELPR IS THE LARGEST RELATIVE SPACING C C CALL FUNCTION SUBPROGRAM R1MACH FROM P O R T MATHEMATICAL C SUBROUTINE LIBRARY (BELL LABORATORIES) C SRELPR = R1MACH(4) P = -0.2E+00 KOUNT = KOUNT + 1 F = 0.0E+00 D = PB - PA U = (X-PA)/D V1 = (1.0E+00-U)/((PG-PC)*X+PH-PD) V = V1*(Y-PC*X-PD) GO TO (10, 20, 30, 40), NUM C SINGULARITY AT U=0, I.E. X=PA C (VERTEX OF 1ST SAMPLE TRIANGLE) 10 IF (U.GT.0.0E+00) F = U**P GO TO 50 C SINGULARITY AT U=V=0, I.E. X=PA AND Y=PC*X+PD C (VERTEX OF 2ND SAMPLE TRIANGLE) 20 A = U + V IF (A.GT.0.0E+00) F = A**P GO TO 50 C SINGULARITY ALONG U+V=1 C (SIDE OF 3RD SAMPLE TRIANGLE) 30 A = 1.0E+00 - U - V IF (A.GT.0.0E+00) F = A**P GO TO 50 C SINGULARITY ALONG U=0 AND V=0 C (TWO SIDES OF 4TH SAMPLE TRIANGLE) 40 A = U*V IF (A.GT.0.0E+00) F = A**P 50 F = F*V1/D RETURN END SUBROUTINE TRIEX(F, VER, EPSABS, EPSREL, LIMEV, RESULT, ABSERR, TRI 10 * IER) C ..................................................................... C 1. TRIEX C DOUBLE INTEGRATION OVER A TRIANGLE C STANDARD FORTRAN SUBROUTINE C C 2. PURPOSE C THE ROUTINE CALCULATES AN APPROXIMATION RESULT TO A C DOUBLE INTEGRAL I = INTEGRAL OF F OVER THE TRIANGLE C WITH VERTICES (VER(1,K),VER(2,K)) FOR K=1,2,3, C HOPEFULLY SATISFYING FOLLOWING CLAIM FOR ACCURACY C ABS(I-RESULT) .LE. MAX(EPSABS,EPSREL*ABS(I)). C C 3. CALLING SEQUENCE C CALL TRIEX (F,VER,EPSABS,EPSREL,LIMEV,RESULT,ABSERR,IER) C C PARAMETERS C F - REAL C FUNCTION DEFINING THE INTEGRAND F(X,Y) C THE ACTUAL NAME FOR F NEEDS TO BE DECLARED C E X T E R N A L IN THE DRIVER PROGRAM. C VER - REAL C ARRAY OF DIMENSION (2,3), WHERE VER(1,K) AND C VER(2,K) REPRESENT THE ABSCISSAE AND THE C ORDINATES RESPECTIVELY OF THE K-TH VERTEX C EPSABS - REAL C ABSOLUTE ACCURACY REQUESTED C EPSREL - REAL C RELATIVE ACCURACY REQUESTED C IF EPSABS .LE. 0 AND EPSREL .LE. 0, C THE ROUTINE WILL END WITH IER .NE. 0. C LIMEV - LIMIT ON THE NUMBER OF INTEGRAND EVALUATIONS C IF MORE THAN XXX EVALUATIONS ARE TO BE ALLOWED, C DIMENSION STATEMENTS AND DATA N IN TRIEX C (WHERE N IS THE MAXIMUM NUMBER OF TRIANGLES C ALLOWED) NEED TO BE ADJUSTED ACCORDING TO C XXX .LE. (N-1)*178/3 + 46, C I.E. PRESENT LIMIT N = 1000 IMPLIES C LIMEV .LE. 59320. C RESULT - REAL C APPROXIMATION TO THE INTEGRAL I C ABSERR - REAL C ESTIMATE OF THE ABSOLUTE ERROR ABS(I-RESULT) C IER - IER = 0 NORMAL AND RELIABLE TERMINATION OF C THE ROUTINE. IT IS ASSUMED THAT THE C REQUESTED ACCURACY HAS BEEN OBTAINED. C - IER .NE. 0 ABNORMAL TERMINATION OF THE ROUTINE C IT IS ASSUMED THAT THE REQUESTED C ACCURACY HAS NOT BEEN OBTAINED. C = 1 MAXIMUM NUMBER OF FUNCTION EVALUATIONS C ALLOWED HAS BEEN REACHED. C = 2 THE PRESENCE OF ROUNDOFF ERROR HAS BEEN C DETECTED, WHICH PREVENTS THE REQUESTED C ACCURACY FROM BEING ACHIEVED. C = 3 IT IS PRESUMED THAT THE TOLERANCE CANNOT C BE ACHIEVED WITHIN THE SPECIFIED NUMBER C OF INTEGRAND EVALUATIONS, AND THAT THE C RESULT RETURNED IS THE BEST WHICH CAN BE C OBTAINED WITHIN THAT LIMIT. C C 4. REMARKS C SUBROUTINES OR FUNCTIONS CALLED C - QUATRI C - DIVIN C - ORDER C - EPSALG C - R1MACH C ..................................................................... C REAL ABSERR, ABSER1, AREA, AREA4, BNDEC, SOVFLO, SRELPR, DRES, * DRESC, DRES0, SUNFLO, R1MACH, ELIST, EMACH, EPRN, EPSABS, * EPSREL, ERLARG, ERLAST, ERRBND, ERRMAX, ERRO4, ERRSUM, F, OFRN, * Q, RCOPY, RESEX, RESULT, RESUL1, RESLA, REX, RLAST, RLIST, STOR, * UFRN, VER, X, Y REAL ABS, AMAX1, AMIN1 INTEGER I, ICNT, ID, IER, IORD, J, J1, J2, J3, JPOINT, JROFF, * JSLAST, JSTAGE, JVRTX, K, KBIGER, LAST, LIMEV, LIMIT, L4, * MAXERR, MAXER0, N, NRMAX, NTR, NUMRL2, NVRTX, N2 INTEGER MAX0, MIN0 LOGICAL EXTRAP DIMENSION ELIST(1000), IORD(1000), JPOINT(4), JSTAGE(1000), * JVRTX(1000,3), RCOPY(52), REX(1000), RLIST(1000), VER(2,3), * X(3000), Y(3000) C COMMON /MACH/ EMACH, EPRN, UFRN, OFRN C C DIMENSIONS C ---------- C N IS THE MAXIMUM NUMBER OF TRIANGLES ALLOWED. C INCREASING N REQUIRES INCREASING DIMENSIONS IN C TRIEX (AND IN SUBROUTINE DIVIN CORRESPONDINGLY) C DIMENSIONS OF ELIST(N),IORD(N),JSTAGE(N),JVRTX(N,3), C REX(N),RLIST(N),X(3*N),Y(3*N) C SHOULD BE SUCH THAT C LIMEV .LE. (N-1)*178/3 + 46 EXTERNAL F DATA N /1000/ C C RCOPY IS OF DIMENSION AT LEAST (LIMEXP+2) C WITH DATA LIMEXP GIVEN IN SUBROUTINE EPSALG C C C MACHINE DEPENDENT CONSTANTS C --------------------------- C SRELPR IS THE LARGEST RELATIVE SPACING C SUNFLO IS THE SMALLEST POSITIVE MAGNITUDE C SOVFLO IS THE LARGEST MAGNITUDE C SRELPR = R1MACH(4) SUNFLO = R1MACH(1) SOVFLO = R1MACH(2) C C LIST OF MAJOR VARIABLES C ----------------------- C C X - LIST OF ABSCISSAE OF ALL VERTICES C Y - LIST OF ORDINATES OF ALL VERTICES C JVRTX - ARRAY OF NUMBERS OF ALL VERTICES C RLIST - LIST OF DEGREE-11 INTEGRAL APPROXIMATIONS FOR C ALL SUBTRIANGLES C REX - LIST OF DEGREE-9 INTEGRAL APPROXIMATIONS FOR C ALL SUBTRIANGLES C ELIST - LIST OF ERROR ESTIMATES FOR ALL SUBTRIANGLES C JSTAGE - LIST OF SUBDIVISION LEVELS FOR ALL C SUBTRIANGLES C AREA - CURRENT SUM OF DEGREE-11 RESULTS C ERRSUM - CURRENT SUM OF ERROR ESTIMATES C RESEX - APPROXIMATION TO THE TOTAL INTEGRAL COMPUTED C AS A SUM OF DEGREE-9 RESULTS OVER THE SMALLEST C SUBTRIANGLES (HIGHEST LEVEL OF SUBDIVISION) AND C OF DEGREE-11 RESULTS OVER THE LARGER TRIANGLES C ERRBND - ESTIMATE OF THE GLOBAL ABSOLUTE TOLERANCE C BNDEC - GLOBAL TOLERANCE OVER THE SET OF LARGER C TRIANGLES C ERLARG - SUM OF ERROR ESTIMATES OVER THE CURRENT SET OF C LARGE TRIANGLES C NTR - CURRENT NUMBER OF TRIANGLES IN THE PARTITION C OF THE ORIGINAL INTEGRATION REGION C MAXERR - POINTER TO SUBTRIANGLE WITH LARGEST ERROR C ESTIMATE C RCOPY - LIST OF THE ENTRIES TO THE EXTRAPOLATION C PROCEDURE C NUMRL2 - NUMBER OF ELEMENTS IN RCOPY C EXTRAP - LOGICAL VARIABLE DENOTING THAT THE ROUTINE C IS ATTEMPTING TO DECREASE THE VALUE OF ERLARG C THROUGH SUBDIVISION OF LARGE (LOW-LEVEL) C TRIANGLES, BEFORE A NEW EXTRAPOLATION IS C CARRIED OUT C C FIRST APPROXIMATION TO THE INTEGRAL C ----------------------------------- IER = 0 I = 0 K = 1 EPRN = 0.5E+01*SRELPR UFRN = SUNFLO/EPRN STOR = 0.0E+00 DO 10 J=1,3 X(J) = VER(1,J) Y(J) = VER(2,J) JVRTX(1,J) = J 10 CONTINUE CALL QUATRI(F, X(1), Y(1), X(2), Y(2), X(3), Y(3), RESEX, RESULT, * ABSERR, DRESC, I, K, STOR) DRES = ABS(RESEX) ERRBND = AMAX1(EPSABS,EPSREL*DRES,EPRN*DRES) IF (ABSERR.LE.ERRBND) GO TO 170 C C INITIALIZATION C -------------- C RLIST(1) = RESULT RCOPY(1) = RESEX ELIST(1) = ABSERR OFRN = SOVFLO JSTAGE(1) = 1 AREA = RESULT ERRSUM = ABSERR ERRMAX = ABSERR EMACH = SRELPR MAXERR = 1 NTR = 1 NRMAX = 1 NVRTX = 3 NUMRL2 = 2 N2 = 2 JROFF = 0 EXTRAP = .FALSE. C C LIMIT = NUMBER OF TRIANGLES ALLOWED IN FINAL PARTITION OF C ORIGINAL TRIANGLE C LIMIT = MAX0(4,MIN0(N,1+((LIMEV-46)/178)*3)) L4 = (LIMIT-8)/4 C C MAIN DO-LOOP C ------------ C DO 130 LAST=2,LIMIT,3 C C SUBDIVIDE THE TRIANGLE WITH NRMAX-TH LARGEST ERROR ESTIMATE C AND INTEGRATE OVER ITS SUBTRIANGLES C CALL DIVIN(MAXERR, X, Y, JVRTX, NTR, NVRTX, JPOINT) MAXER0 = MAXERR RLAST = RLIST(MAXERR) JSLAST = JSTAGE(MAXERR) ERLAST = ERRMAX AREA4 = 0.0E+00 ERRO4 = 0.0E+00 K = 1 ICNT = 0 STOR = 0.0E+00 DO 20 I=1,4 J = JPOINT(I) J1 = JVRTX(J,1) J2 = JVRTX(J,2) J3 = JVRTX(J,3) CALL QUATRI(F, X(J1), Y(J1), X(J2), Y(J2), X(J3), Y(J3), * REX(J), RLIST(J), ELIST(J), DRESC, I, K, STOR) AREA4 = AREA4 + RLIST(J) ERRO4 = ERRO4 + ELIST(J) JSTAGE(J) = JSLAST + 1 IF (ELIST(J).EQ.DRESC) ICNT = 1 20 CONTINUE C C IMPROVE PREVIOUS INTEGRAL APPROXIMATION AND ERROR C ESTIMATE AND TEST FOR ACCURACY C AREA = AREA + AREA4 - RLAST ERRSUM = ERRSUM + ERRO4 - ERRMAX IF (ERRSUM.LE.ERRBND) GO TO 150 C C TEST FOR ROUNDOFF ERROR AND EVENTUALLY SET ERROR FLAG C IF (ICNT.EQ.1) GO TO 30 Q = 0.0E+00 IF (AREA4.NE.0.0E+00) Q = RLAST/AREA4 IF (0.9999E+00.LT.Q .AND. Q.LT.1.0001E+00 .AND. * ERRMAX.LE.ERRO4) JROFF = JROFF + 1 IF (JROFF.GE.10) IER = 2 C C SET ERROR FLAG IN THE CASE THAT MAXIMUM NUMBER OF C SUBTRIANGLES IS REACHED C 30 IF (LAST.GE.LIMIT-2) IER = 1 IF (IER.NE.0) GO TO 140 C C CALL SUBROUTINE ORDER TO MAINTAIN THE DESCENDING C ORDERING IN THE LIST OF ERROR ESTIMATES, AND TO SELECT C THE SUBTRIANGLE WITH NRMAX-TH LARGEST ERROR ESTIMATE C (TO BE SUBDIVIDED NEXT) C CALL ORDER(LIMIT, LAST, MAXERR, ERRMAX, ELIST, IORD, NRMAX, NTR) IF (LAST.EQ.2) GO TO 120 C C CHECK WHETHER THE TRIANGLES OBTAINED IN THE FOREGOING C SUBDIVISION ARE LARGE, AND ADJUST ERLARG IF SO C ERLARG = ERLARG - ERLAST IF (JSTAGE(MAXER0).LE.NUMRL2) ERLARG = ERLARG + ERRO4 IF (EXTRAP) GO TO 40 C C CHECK WHETHER THE FOLLOWING SUBDIVISION WOULD PROVIDE C LARGE TRIANGLES, AND PROCEED IMMEDIATELY TO THIS C SUBDIVISION IF SO C IF (JSTAGE(MAXERR).LE.NUMRL2) GO TO 130 C C DECREASE THE INTEGRATION ERROR OVER THE COLLECTION OF C LARGE TRIANGLES C EXTRAP = .TRUE. NRMAX = 2 40 BNDEC = AMAX1(EPRN*DRES,AMIN1(0.1E+00*ERRBND,0.1E-02*DRES)) IF (ERLARG.LE.BNDEC) GO TO 60 ID = NRMAX DO 50 J=ID,NTR MAXERR = IORD(NRMAX) ERRMAX = ELIST(MAXERR) IF (JSTAGE(MAXERR).LE.NUMRL2) GO TO 130 NRMAX = NRMAX + 1 50 CONTINUE C C COMPUTE NEXT ENTRY FOR THE EXTRAPOLATION PROCEDURE C 60 RESEX = 0.0E+00 DO 70 J=1,NTR IF (JSTAGE(J).GT.NUMRL2) RESEX = RESEX + REX(J) IF (JSTAGE(J).LE.NUMRL2) RESEX = RESEX + RLIST(J) 70 CONTINUE DRES = ABS(RESEX) ERRBND = AMAX1(EPSABS,EPSREL*DRES,EPRN*DRES) C C EXTRAPOLATE AND TEST FOR ACCURACY C NUMRL2 = NUMRL2 + 1 N2 = N2 + 1 RCOPY(NUMRL2) = RESEX CALL EPSALG(NUMRL2, RCOPY, RESUL1, ABSER1, RESLA) ABSER1 = ABSER1 + ERLARG IF (NUMRL2.EQ.N2 .OR. N2.EQ.3) GO TO 80 IF (ABSER1.GT.ABSERR) GO TO 90 80 RESULT = RESUL1 ABSERR = ABSER1 DRES0 = ABS(RESULT) IF (ABSERR.LE.AMAX1(EPSABS,EPSREL*DRES0,EPRN*DRES0) .AND. * N2.GT.3) GO TO 170 C C IN THE CASE THAT THE EXTRAPOLATION ERROR ESTIMATE IS C CONSIDERABLY SMALLER THAN THE ERROR SUM OVER ALL C SUBTRIANGLES, CHECK IF ENOUGH STEPS ARE LEFT OVER C IN ORDER TO COMPLETE ANOTHER EXTRAPOLATION C STAGE. SET ERROR FLAG FOR ABNORMAL TERMINATION C IF NECESSARY C 90 IF (ABSERR*ABS(RESEX).GE.0.5E-01*ERRSUM*DRES0 .OR. LAST.LE.L4) * GO TO 110 KBIGER = 0 DO 100 J=1,NTR IF (ELIST(J).GT.BNDEC) KBIGER = KBIGER + 1 100 CONTINUE IF (KBIGER.GE.(LIMIT-LAST-2)/3) IER = 3 110 IF (IER.NE.0) GO TO 140 C C START NEW EXTRAPOLATION STAGE, BY INTEGRATING OVER C SUBTRIANGLE WITH LARGEST ERROR ESTIMATE C MAXERR = IORD(1) ERRMAX = ELIST(MAXERR) NRMAX = 1 ERLARG = ERRSUM GO TO 130 C C SECOND EXTRAPOLATION ENTRY C 120 ERLARG = ERRSUM RESEX = REX(1) + REX(2) + REX(3) + REX(4) RCOPY(2) = RESEX DRES = ABS(RESEX) RESLA = RESEX 130 CONTINUE C C IN THE CASE THAT IER .NE. 0, SELECT THE BEST APPROXIMATION C ---------------------------------------------------------- C (RESULT (EXTRAPOLATED), OR AREA (SUM OF INTEGRALS OVER C SUBTRIANGLES) ) C 140 DRES = ABS(AREA) IF (DRES*ABSERR.LE.DRES0*ERRSUM) GO TO 170 IF (ABS(AREA-RESEX)*DRES0.GT.ABS(RESULT-RESEX)*DRES) GO TO 170 150 RESULT = 0.0E+00 DO 160 J=1,NTR RESULT = RESULT + RLIST(J) 160 CONTINUE ABSERR = ERRSUM 170 RETURN END SUBROUTINE DIVIN(MAX, X, Y, IVRTX, NTR, NVRTX, IPOINT) DIV 10 C ..................................................................... C C 1. DIVIN C C SUBDIVISION OF A GIVEN TRIANGLE C STANDARD FORTRAN SUBROUTINE C C 2. PURPOSE C C THE ROUTINE SUBDIVIDES THE GIVEN TRIANGLE WITH INDEX MAX C INTO 4 SIMILAR SUBTRIANGLES, BY HALVING THE SIDES OF THE C GIVEN TRIANGLE. IT CALCULATES THE COORDINATES OF THE NEW C VERTICES, AND DEFINES THE SUBTRIANGLES BY DEFINING POINTERS C TO THESE TRIANGLES AND TO THEIR VERTICES. C C 3. CALLING SEQUENCE C C CALL DIVIN(MAX,X,Y,IVRTX,NTR,NVRTX,IPOINT) C C PARAMETERS C MAX - INTEGER C INDEX OF THE TRIANGLE WHICH MUST BE SUBDIVIDED C C X - REAL C ABSCISSAE OF THE VERTICES OF ALL TRIANGLES IN THE C PARTITION OF THE INTEGRATION DOMAIN C C Y - REAL C ORDINATES OF THE VERTICES OF ALL TRIANGLES IN THE C PARTITION OF THE INTEGRATION DOMAIN C C IVRTX - INTEGER C IVRTX(J,K) POINTS TO THE K-TH VERTEX (K=1,2,3) C OF THE TRIANGLE WITH INDEX J C C NTR - INTEGER C TOTAL NUMBER OF TRIANGLES IN THE PARTITION C NTR IS INCREASED BY 3 DURING THE COURSE OF THIS C SUBROUTINE C C NVRTX - INTEGER C TOTAL NUMBER OF VERTICES IN THE PARTITION C NVRTX IS INCREASED BY 3 DURING THIS CALL. C C IPOINT - INTEGER C POINTERS TO THE 4 NEW SUBTRIANGLES, SUCH THAT C IPOINT(1)=MAX, IPOINT(2)=NTR+1, IPOINT(3)=NTR+2, C IPOINT(4)=NTR+3 C C 4. REMARKS C C SUBROUTINE OR FUNCTIONS CALLED NONE C C ..................................................................... REAL X, Y INTEGER I, I1, I2, I3, I4, I5, I6, IP, IPOINT, IVRTX, MAX, NTR, * NVRTX C DIMENSION IPOINT(4), IVRTX(1000,3), X(3000), Y(3000) C C DIMENSIONS C ---------- C FOR ALL DIMENSION SPECIFICATIONS SEE CALLING ROUTINE TRIEX C C LIST OF MAJOR VARIABLES C ----------------------- C I1,...,I6 - INDICES OF VERTICES OF THE 4 NEW SUBTRIANGLES C C C SET POINTERS TO THE 4 NEW SUBTRIANGLES C DO 10 I=1,3 IP = NTR + I IPOINT(I+1) = IP 10 CONTINUE IPOINT(1) = MAX C C DEFINE VERTICES OF THE SUBTRIANGLES C I1 = IVRTX(MAX,1) I2 = IVRTX(MAX,2) I3 = IVRTX(MAX,3) I4 = NVRTX + 1 I5 = NVRTX + 2 I6 = NVRTX + 3 X(I4) = (X(I2)+X(I3))*0.5E+00 Y(I4) = (Y(I2)+Y(I3))*0.5E+00 X(I5) = (X(I3)+X(I1))*0.5E+00 Y(I5) = (Y(I3)+Y(I1))*0.5E+00 X(I6) = (X(I1)+X(I2))*0.5E+00 Y(I6) = (Y(I1)+Y(I2))*0.5E+00 C C ADJUST NVRTX C NVRTX = I6 C C THE SUBDIVISION IS ORGANIZED AS IS INDICATED IN THE C PICTURES BELOW. THE POINTERS TO THE VERTICES OF THE C SUBTRIANGLES ARE SET IN SUCH A WAY THAT C SUBDIVISION OF THE TRIANGLE WITH INDEX MAX PRODUCES THE C TRIANGLES WITH INDICES MAX, NTR+1, NTR+2, NTR+3 C (DETERMINED BY (I1,I6,I5), (I6,I2,I4), (I5,I4,I3) AND C (I4,I6,I5) RESPECTIVELY). C C I3 C C C I5 I4 C C C I1 I6 I2 C IVRTX(MAX,1) = I1 IVRTX(MAX,2) = I6 IVRTX(MAX,3) = I5 NTR = NTR + 1 IVRTX(NTR,1) = I6 IVRTX(NTR,2) = I2 IVRTX(NTR,3) = I4 NTR = NTR + 1 IVRTX(NTR,1) = I5 IVRTX(NTR,2) = I4 IVRTX(NTR,3) = I3 NTR = NTR + 1 IVRTX(NTR,1) = I4 IVRTX(NTR,2) = I6 IVRTX(NTR,3) = I5 RETURN END SUBROUTINE QUATRI(F, U1, V1, U2, V2, U3, V3, RES9, RES11, EST, QUA 10 * DRESC, I, K, STOR) C ..................................................................... C C 1. QUATRI C C INTEGRATION RULES C STANDARD FORTRAN SUBROUTINE C C 2. PURPOSE C C TO COMPUTE - IF = INTEGRAL OF F OVER THE TRIANGLE C WITH VERTICES (U1,V1),(U2,V2),(U3,V3), AND ERROR C ESTIMATE, C - INTEGRAL OF ABS(F) OVER THIS TRIANGLE C C 3. CALLING SEQUENCE C C CALL QUATRI(F,U1,V1,U2,V2,U3,V3,RES9,RES11,EST,DRESC, C I,K,STOR) C C PARAMETERS C F - REAL C FUNCTION SUBPROGRAM DEFINING THE INTEGRAND C F(X,Y) THE ACTUAL NAME FOR F NEEDS TO BE C DECLARED E X T E R N A L IN THE CALLING C PROGRAM. C C U1,U2,U3- REAL C ABSCISSAE OF VERTICES C C V1,V2,V3- REAL C ORDINATES OF VERTICES C C RES9 - REAL C APPROXIMATION TO THE INTEGRAL IF, OBTAINED BY THE C LYNESS AND JESPERSEN RULE OF DEGREE 9, USING C 19 POINTS C C RES11 - REAL C APPROXIMATION TO THE INTEGRAL IF, OBTAINED BY THE C LYNESS AND JESPERSEN RULE OF DEGREE 11, C USING 28 POINTS C C EST - REAL C ESTIMATE OF THE ABSOLUTE ERROR C C DRESC - REAL C APPROXIMATION TO THE INTEGRAL OF ABS(F- IF/DJ), C OBTAINED BY THE RULE OF DEGREE 9, AND USED FOR C THE COMPUTATION OF EST C C I - INTEGER C NUMBER OF TRIANGLES IN CURRENT SUBDIVISION C (1 .LE. I .LE. 4) C C K - INTEGER C INDEX OF POINTS TO BE SHARED WITH 4TH TRIANGLE C OF SUBDIVISION C C STOR - REAL C SUM OF FUNCTION VALUES AT THESE POINTS C C 4. REMARKS C C SUBROUTINES OR FUNCTIONS CALLED C - F (USER-SUPPLIED INTEGRAND FUNCTION) C C ..................................................................... C REAL DJ, DF0, DRESC, EMACH, EPRN, EST, F, FSTOR, FV, F0, OFRN, * RESAB9, RES9, RES11, STOR, U1, U2, U3, UFRN, V1, V2, V3, W, W90, * W110, X, Y, ZETA1, ZETA2, Z1, Z2, Z3 REAL AMAX1, AMIN1 INTEGER I, J, JEND, K, KOUNT, L C DIMENSION FSTOR(3), FV(19), W(15), X(3), Y(3), ZETA1(15), * ZETA2(15) C COMMON /MACH/ EMACH, EPRN, UFRN, OFRN C C FIRST HOMOGENEOUS COORDINATES OF POINTS IN DEGREE-9 C AND DEGREE-11 FORMULA, TAKEN WITH MULTIPLICITY 3 DATA ZETA1(1), ZETA1(2), ZETA1(3), ZETA1(4), ZETA1(5), ZETA1(6), * ZETA1(7), ZETA1(8), ZETA1(9), ZETA1(10), ZETA1(11), ZETA1(12), * ZETA1(13), ZETA1(14), ZETA1(15) /0.2063496160252593E-01, * 0.1258208170141290E+00,0.6235929287619356E+00, * 0.9105409732110941E+00,0.3683841205473626E-01, * 0.7411985987844980E+00,0.9480217181434233E+00, * 0.8114249947041546E+00,0.1072644996557060E-01, * 0.5853132347709715E+00,0.1221843885990187E+00, * 0.4484167758913055E-01,0.6779376548825902E+00,0.0E+00, * 0.8588702812826364E+00/ C SECOND HOMOGENEOUS COORDINATES OF POINTS IN DEGREE-9 C AND DEGREE-11 FORMULA, TAKEN WITH MULTIPLICITY 3 DATA ZETA2(1), ZETA2(2), ZETA2(3), ZETA2(4), ZETA2(5), ZETA2(6), * ZETA2(7), ZETA2(8), ZETA2(9), ZETA2(10), ZETA2(11), ZETA2(12), * ZETA2(13), ZETA2(14), ZETA2(15) /0.4896825191987370E+00, * 0.4370895914929355E+00,0.1882035356190322E+00, * 0.4472951339445297E-01,0.7411985987844980E+00, * 0.3683841205473626E-01,0.2598914092828833E-01, * 0.9428750264792270E-01,0.4946367750172147E+00, * 0.2073433826145142E+00,0.4389078057004907E+00, * 0.6779376548825902E+00,0.4484167758913055E-01, * 0.8588702812826364E+00,0.0E+00/ C WEIGHTS OF MIDPOINT OF TRIANGLE IN DEGREE-9 C RESP. DEGREE-11 FORMULAE DATA W90 /0.9713579628279610E-01/ DATA W110 /0.8797730116222190E-01/ C WEIGHTS IN DEGREE-9 AND DEGREE-11 RULE DATA W(1), W(2), W(3), W(4), W(5), W(6), W(7), W(8), W(9), W(10), * W(11), W(12), W(13), W(14), W(15) /0.3133470022713983E-01, * 0.7782754100477543E-01,0.7964773892720910E-01, * 0.2557767565869810E-01,0.4328353937728940E-01, * 0.4328353937728940E-01,0.8744311553736190E-02, * 0.3808157199393533E-01,0.1885544805613125E-01, * 0.7215969754474100E-01,0.6932913870553720E-01, * 0.4105631542928860E-01,0.4105631542928860E-01, * 0.7362383783300573E-02,0.7362383783300573E-02/ C C LIST OF MAJOR VARIABLES C ---------------------- C DJ - AREA OF THE TRIANGLE C DRESC - APPROXIMATION TO INTEGRAL OF C ABS(F- IF/DJ) OVER THE TRIANGLE C RESAB9 - APPROXIMATION TO INTEGRAL OF C ABS(F) OVER THE TRIANGLE C X - CARTESIAN ABSCISSAE OF THE INTEGRATION C POINTS C Y - CARTESIAN ORDINATES OF THE INTEGRATION C POINTS C FV - FUNCTION VALUES C FSTOR C C COMPUTE DEGREE-9 AND DEGREE-11 RESULTS FOR IF/DJ AND C DEGREE-9 APPROXIMATION FOR ABS(F) C DJ = ABS(U1*V2-U2*V1-U1*V3+V1*U3+U2*V3-V2*U3)*0.5E+00 F0 = F((U1+U2+U3)/3.0E+00,(V1+V2+V3)/3.0E+00) RES9 = F0*W90 RESAB9 = ABS(F0)*W90 FV(1) = F0 KOUNT = 1 RES11 = F0*W110 IF (I.EQ.4) RES11 = RES11 + STOR*W(15) JEND = 15 IF (I.EQ.4) JEND = 13 DO 50 J=1,JEND Z1 = ZETA1(J) Z2 = ZETA2(J) Z3 = 1.0E+00 - Z1 - Z2 X(1) = Z1*U1 + Z2*U2 + Z3*U3 Y(1) = Z1*V1 + Z2*V2 + Z3*V3 X(2) = Z2*U1 + Z3*U2 + Z1*U3 Y(2) = Z2*V1 + Z3*V2 + Z1*V3 X(3) = Z3*U1 + Z1*U2 + Z2*U3 Y(3) = Z3*V1 + Z1*V2 + Z2*V3 IF (J.GT.6) GO TO 20 F0 = 0.0E+00 DF0 = 0.0E+00 DO 10 L=1,3 KOUNT = KOUNT + 1 FV(KOUNT) = F(X(L),Y(L)) F0 = F0 + FV(KOUNT) DF0 = DF0 + ABS(FV(KOUNT)) 10 CONTINUE RES9 = RES9 + F0*W(J) RESAB9 = RESAB9 + DF0*W(J) GO TO 50 20 F0 = 0.0E+00 DO 30 L=1,3 FSTOR(L) = F(X(L),Y(L)) F0 = F0 + FSTOR(L) 30 CONTINUE IF (J.LT.14) GO TO 40 STOR = STOR + FSTOR(K) K = K + 1 IF (K.GT.3) K = 1 40 RES11 = RES11 + F0*W(J) 50 CONTINUE C C COMPUTE DEGREE-9 APPROXIMATION FOR THE INTEGRAL OF C ABS(F-IF/DJ) C DRESC = ABS(FV(1)-RES9)*W90 KOUNT = 2 DO 60 J=1,6 DRESC = DRESC + (ABS(FV(KOUNT)-RES9)+ABS(FV(KOUNT+1)-RES9) * +ABS(FV(KOUNT+2)-RES9))*W(J) KOUNT = KOUNT + 3 60 CONTINUE C C COMPUTE DEGREE-9 AND DEGREE-11 APPROXIMATIONS FOR IF, C AND ERROR ESTIMATE C RES9 = RES9*DJ RES11 = RES11*DJ RESAB9 = RESAB9*DJ DRESC = DRESC*DJ EST = ABS(RES9-RES11) IF (DRESC.NE.0.0E+00) EST = AMAX1(EST,DRESC*AMIN1(1.0E+00, * (20.0E+00*EST/DRESC)**1.5E+00)) IF (RESAB9.GT.UFRN) EST = AMAX1(EPRN*RESAB9,EST) RETURN END SUBROUTINE EPSALG(N, EPSTAB, RESULT, ABSERR, RESLA) EPS 10 C C ................................................................ C C 1. EPSALG C C EPSILON ALGORITHM C STANDARD FORTRAN SUBROUTINE C C 2. PURPOSE C C THE ROUTINE TRANSFORMS A GIVEN SEQUENCE USING THE C EPSILON ALGORITHM OF P. WYNN. C AN ESTIMATE OF THE ABSOLUTE ERROR IS ALSO GIVEN. C THE CONDENSED EPSILON TABLE IS COMPUTED. ONLY THOSE C ELEMENTS NEEDED FOR THE COMPUTATION OF THE NEXT DIAGONAL C ARE PRESERVED. C C 3. CALLING SEQUENCE C C CALL EPSALG (N,EPSTAB,RESULT,ABSERR,RES3LA,NRES) C C PARAMETERS C N - INTEGER C EPSTAB(N) CONTAINS THE NEW ELEMENT IN THE C FIRST COLUMN OF THE EPSILON TABLE. C C EPSTAB - REAL C ONE DIMENSIONAL ARRAY CONTAINING THE C ELEMENTS OF THE LAST TWO DIAGONALS OF C THE TRIANGULAR EPSILON TABLE C THE ELEMENTS ARE NUMBERED STARTING AT THE C RIGHT-HAND CORNER OF THE TRIANGLE. C THE DIMENSION SHOULD BE AT LEAST (LIMEXP+2) C (SEE DATA LIMEXP). C C RESULT - REAL C RESULTING APPROXIMATION TO THE INTEGRAL C C ABSERR - REAL C ESTIMATE OF THE ABSOLUTE ERROR C C RESLA - REAL C PREVIOUS RESULT C C 4. REMARKS C C SUBROUTINES OR FUNCTIONS CALLED NONE C C .................................................................. C REAL ABSERR, DELTA1, DELTA2, DELTA3, EMACH, EPRN, EPSINF, EPSTAB, * ERROR, ERR1, ERR2, ERR3, E0, E1, E2, E3, E1ABS, OFRN, RES, * RESULT, RESLA, SS, TOL1, TOL2, TOL3, UFRN REAL ABS, AMAX1 INTEGER I, IB, IB2, IE, IN, K1, K2, K3, LIMEXP, N, NEWELM, NUM DIMENSION EPSTAB(52) C COMMON /MACH/ EMACH, EPRN, UFRN, OFRN C C LIMEXP IS THE MAXIMUM NUMBER OF ELEMENTS THE EPSILON C TABLE CAN CONTAIN. IF THIS NUMBER IS REACHED, THE UPPER C DIAGONAL OF THE EPSILON TABLE IS DELETED. C ( EPSTAB IS OF DIMENSION (LIMEXP+2) AT LEAST.) C DATA LIMEXP /50/ C C LIST OF MAJOR VARIABLES C ----------------------- C E0 - THE 4 ELEMENTS ON WHICH THE C E1 COMPUTATION OF A NEW ELEMENT IN C E2 THE EPSILON TABLE IS BASED C E3 E0 C E3 E1 NEW C E2 C NEWELM - NUMBER OF ELEMENTS TO BE COMPUTED IN THE NEW C DIAGONAL C ERROR - ERROR = ABS(E0-E1)+ABS(E1-E2)+ABS(E2-NEW) C RESULT - THE ELEMENT IN THE NEW DIAGONAL WITH LEAST C ERROR C EPSTAB(N+2) = EPSTAB(N) NEWELM = (N-1)/2 EPSTAB(N) = OFRN ABSERR = OFRN NUM = N K1 = N DO 40 I=1,NEWELM K2 = K1 - 1 K3 = K1 - 2 RES = EPSTAB(K1+2) E0 = EPSTAB(K3) E1 = EPSTAB(K2) E2 = RES E1ABS = ABS(E1) DELTA2 = E2 - E1 ERR2 = ABS(DELTA2) TOL2 = AMAX1(ABS(E2),E1ABS)*EMACH DELTA3 = E1 - E0 ERR3 = ABS(DELTA3) TOL3 = AMAX1(E1ABS,ABS(E0))*EMACH IF (ERR2.GT.TOL2 .OR. ERR3.GT.TOL3) GO TO 10 C C IF E0, E1 AND E2 ARE EQUAL TO WITHIN MACHINE C ACCURACY, CONVERGENCE IS ASSUMED. C RESULT = E2 C ABSERR = ABS(E1-E0)+ABS(E2-E1) C RESULT = RES ABSERR = ERR2 + ERR3 GO TO 90 10 E3 = EPSTAB(K1) EPSTAB(K1) = E1 DELTA1 = E1 - E3 ERR1 = ABS(DELTA1) TOL1 = AMAX1(E1ABS,ABS(E3))*EMACH C C IF TWO ELEMENTS ARE VERY CLOSE TO EACH OTHER, OMIT C A PART OF THE TABLE BY ADJUSTING THE VALUE OF N. C IF (ERR1.LE.TOL1 .OR. ERR2.LE.TOL2 .OR. ERR3.LE.TOL3) GO TO 20 SS = 0.1E+01/DELTA1 + 0.1E+01/DELTA2 - 0.1E+01/DELTA3 EPSINF = ABS(SS*E1) C C TEST TO DETECT IRREGULAR BEHAVIOUR IN THE TABLE, AND C EVENTUALLY OMIT A PART OF THE TABLE ADJUSTING THE VALUE C OF N. C IF (EPSINF.GT.0.1E-03) GO TO 30 20 N = I + I - 1 GO TO 50 C C COMPUTE A NEW ELEMENT AND EVENTUALLY ADJUST C THE VALUE OF RESULT C 30 RES = E1 + 0.1E+01/SS EPSTAB(K1) = RES K1 = K1 - 2 ERROR = ERR2 + ABS(RES-E2) + ERR3 IF (ERROR.GT.ABSERR) GO TO 40 ABSERR = ERROR RESULT = RES 40 CONTINUE C C SHIFT THE TABLE C 50 IF (N.EQ.LIMEXP) N = 2*(LIMEXP/2) - 1 IB = 1 IF ((NUM/2)*2.EQ.NUM) IB = 2 IE = NEWELM + 1 DO 60 I=1,IE IB2 = IB + 2 EPSTAB(IB) = EPSTAB(IB2) IB = IB2 60 CONTINUE IF (NUM.EQ.N) GO TO 80 IN = NUM - N + 1 DO 70 I=1,N EPSTAB(I) = EPSTAB(IN) IN = IN + 1 70 CONTINUE C C COMPUTE ERROR ESTIMATE C 80 ABSERR = ABS(RESULT-RESLA) RESLA = RESULT 90 ABSERR = AMAX1(ABSERR,EPRN*ABS(RESULT)) RETURN END SUBROUTINE ORDER(LIMIT, LAST, MAXERR, ERMAX, ELIST, IORD, NRMAX, ORD 10 * NAVAIL) C C ..................................................................... C C 1. ORDER C C SORTING ROUTINE C STANDARD FORTRAN SUBROUTINE C C 2. PURPOSE C C THIS ROUTINE MAINTAINS THE DESCENDING ORDERING C IN THE LIST OF THE ERROR ESTIMATES RESULTING FROM THE C SUBDIVISION PROCESS. AT EACH CALL 4 ERROR ESTIMATES C ARE INSERTED USING THE SEQUENTIAL SEARCH METHOD C TOP-DOWN FOR THE LARGEST ESTIMATE, BOTTOM UP C FOR THE OTHER ONES. C C 3. CALLING SEQUENCE C C CALL ORDER(LIMIT,LAST,MAXERR,ERMAX,ELIST,IORD,NRMAX, C NAVAIL) C C PARAMETERS C LIMIT - INTEGER C MAXIMUM NUMBER OF ERROR ESTIMATES C THE LIST CAN CONTAIN C C LAST - INTEGER C NUMBER OF ERROR ESTIMATES CURRENTLY IN C THE LIST ELIST(LAST) CONTAINS THE SMALLEST C ERROR ESTIMATE C C MAXERR - INTEGER C POINTS TO THE NRMAX-TH LARGEST ERROR C ESTIMATE CURRENTLY IN THE LIST C C ERMAX - REAL C NRMAX-TH LARGEST ERROR ESTIMATE C ERMAX = ELIST(MAXERR) C C ELIST - REAL C ARRAY OF DIMENSION LIMIT CONTAINING C THE ERROR ESTIMATES C C IORD - INTEGER C ARRAY CONTAINING POINTERS TO ELIST SO C THAT IORD(1) POINTS TO THE LARGEST C ERROR ESTIMATE,..., IORD(LAST) TO THE C SMALLEST ERROR ESTIMATE C C NRMAX - INTEGER C MAXERR = IORD(NRMAX) C C NAVAIL - INTEGER C NAVAIL = LAST+2 C C 4. REMARKS C C SUBROUTINES OR FUNCTIONS CALLED NONE C C ..................................................................... REAL ELIST, ERMAX, ERRMAX, ERRMIN INTEGER I, IBEG, II, INCR, IND, INDI, INDJ, IORD, IORD4, ISUCC, * J, JBN, K, LAST, LIMIT, MAXERR, MINERR, NAVAIL, NRMAX, NTEST, NUM DIMENSION ELIST(LIMIT), IORD(LIMIT), IORD4(4) C C SORT THE 4 NEW ELIST-ENTRIES INTO DESCENDING ORDERING C IORD4(1) = MAXERR IORD4(2) = LAST IORD4(3) = LAST + 1 IORD4(4) = LAST + 2 DO 20 K=1,3 I = K DO 10 II=1,3 J = I + 1 INDI = IORD4(I) INDJ = IORD4(J) IF (ELIST(INDI).GE.ELIST(INDJ)) GO TO 20 IORD4(I) = INDJ IORD4(J) = INDI I = I - 1 IF (I.EQ.0) GO TO 20 10 CONTINUE 20 CONTINUE C C IF THE LIST CONTAINS ONLY 4 ELEMENTS, RETURN C IF (LAST.GT.2) GO TO 40 DO 30 I=1,4 IORD(I) = IORD4(I) 30 CONTINUE GO TO 170 C C THIS PART OF THE ROUTINE IS ONLY EXECUTED IF, DUE TO A C DIFFICULT INTEGRAND, SUBDIVISION INCREASED THE ERROR C ESTIMATE. IN THE NORMAL CASE THE INSERT PROCEDURE C SHOULD START AFTER THE NRMAX-TH LARGEST ERROR ESTIMATE. C 40 MAXERR = IORD4(1) ERRMAX = ELIST(MAXERR) IF (NRMAX.EQ.1) GO TO 60 II = NRMAX - 1 DO 50 I=1,II ISUCC = IORD(NRMAX-1) IF (ERRMAX.LE.ELIST(ISUCC)) GO TO 60 IORD(NRMAX) = ISUCC NRMAX = NRMAX - 1 50 CONTINUE C C INSERT ERRMAX BY TRAVERSING THE LIST TOP-DOWN C STARTING COMPARISON FROM THE ELEMENT ELIST (IORD(NRMAX+1)) C 60 NAVAIL = LAST + 2 JBN = NAVAIL - 3 IBEG = NRMAX + 1 IF (IBEG.GT.JBN) GO TO 80 DO 70 I=IBEG,JBN ISUCC = IORD(I) IF (ERRMAX.GE.ELIST(ISUCC)) GO TO 100 IORD(I-1) = ISUCC 70 CONTINUE 80 IND = JBN DO 90 I=1,4 IORD(IND) = IORD4(I) IND = IND + 1 90 CONTINUE GO TO 170 C C INSERT THE OTHER ESTIMATES BY TRAVERSING THE LIST BOTTOM-UP C 100 IORD(I-1) = MAXERR K = JBN NUM = JBN - I + 1 INCR = 3 DO 160 II=1,3 NTEST = NUM MINERR = IORD4(INCR+1) ERRMIN = ELIST(MINERR) DO 110 J=1,NTEST ISUCC = IORD(K) IF (ERRMIN.LT.ELIST(ISUCC)) GO TO 150 IND = K + INCR IORD(IND) = ISUCC K = K - 1 NUM = NUM - 1 110 CONTINUE GO TO (140, 130, 120), INCR 120 IORD(I+2) = IORD4(4) 130 IORD(I+1) = IORD4(3) 140 IORD(I) = IORD4(2) GO TO 170 150 IND = K + INCR IORD(IND) = MINERR INCR = INCR - 1 160 CONTINUE C C SET MAXERR AND ERMAX C 170 MAXERR = IORD(NRMAX) ERMAX = ELIST(MAXERR) RETURN END