C ALGORITHM 706 , COLLECTED ALGORITHMS FROM ACM. C THIS WORK PUBLISHED IN TRANSACTIONS ON MATHEMATICAL SOFTWARE, C VOL. 18, NO. 3, SEPTEMBER, 1992, PP. 329-342. C C Remark: TOMS Vol. 24, No. 3, September, 1998, pp. 334-335. C C This file summarizes the DOUBLE PRECISION routines C in this directory. SINGLE PRECISION names are obtained C by replacing the leading D with S in all names. C The machine constant routines D1MACH and R1MACH C are also included. C C DTEST1 Simple test driver for DCUTRI. C Sample output from a SUN SPARC station 1 is included. C DTEST2 More detailed test driver for DCUTRI. C Sample output from a SUN SPARC station 1 is included. C DCUTRI Main Integrator.DCUTRI calls DCHTRI and DADTRI. C DCHTRI Checks the input to DCUTRI. C DADTRI The adaptive integration routine. C DADTRI calls DTRTRI and DRLTRI. C DTRTRI Maintaines the heap of subtriangles. C DRLTRI Computes estimates of integral and error over C subtriangles. dtest1: dtest1.o dcutri.o dadtri.o dchtri.o \ dtrtri.o drltri.o d1mach.o f77 -o dtest1 dtest1.o dcutri.o dadtri.o dchtri.o \ dtrtri.o drltri.o d1mach.o dtest2: dtest2.o dcutri.o dadtri.o dchtri.o \ dtrtri.o drltri.o d1mach.o f77 -o dtest2 dtest2.o dcutri.o dadtri.o dchtri.o \ dtrtri.o drltri.o d1mach.o dtest1.o: dtest1.f f77 -ansi -c dtest1.f dtest2.o: dtest2.f f77 -ansi -c dtest2.f dcutri.o: dcutri.f f77 -ansi -c dcutri.f dadtri.o: dadtri.f f77 -ansi -c dadtri.f dchtri.o: dchtri.f f77 -ansi -c dchtri.f dtrtri.o: dtrtri.f f77 -ansi -c dtrtri.f drltri.o: drltri.f f77 -ansi -c drltri.f d1mach.o: d1mach.f f77 -ansi -c d1mach.f stest1: stest1.o scutri.o sadtri.o schtri.o \ strtri.o srltri.o r1mach.o f77 -o stest1 stest1.o scutri.o sadtri.o schtri.o \ strtri.o srltri.o r1mach.o stest2: stest2.o scutri.o sadtri.o schtri.o \ strtri.o srltri.o r1mach.o f77 -o stest2 stest2.o scutri.o sadtri.o schtri.o \ strtri.o srltri.o r1mach.o stest1.o: stest1.f f77 -ansi -c stest1.f stest2.o: stest2.f f77 -ansi -c stest2.f scutri.o: scutri.f f77 -ansi -c scutri.f sadtri.o: sadtri.f f77 -ansi -c sadtri.f schtri.o: schtri.f f77 -ansi -c schtri.f strtri.o: strtri.f f77 -ansi -c strtri.f srltri.o: srltri.f f77 -ansi -c srltri.f r1mach.o: r1mach.f f77 -ansi -c r1mach.f DOUBLE PRECISION FUNCTION D1MACH(I) C C DOUBLE-PRECISION MACHINE CONSTANTS C C D1MACH( 1) = B**(EMIN-1), THE SMALLEST POSITIVE MAGNITUDE. C C D1MACH( 2) = B**EMAX*(1 - B**(-T)), THE LARGEST MAGNITUDE. C C D1MACH( 3) = B**(-T), THE SMALLEST RELATIVE SPACING. C C D1MACH( 4) = B**(1-T), THE LARGEST RELATIVE SPACING. C C D1MACH( 5) = LOG10(B) C C TO ALTER THIS FUNCTION FOR A PARTICULAR ENVIRONMENT, C THE DESIRED SET OF DATA STATEMENTS SHOULD BE ACTIVATED BY C REMOVING THE C FROM COLUMN 1. C ON RARE MACHINES A STATIC STATEMENT MAY NEED TO BE ADDED. C (BUT PROBABLY MORE SYSTEMS PROHIBIT IT THAN REQUIRE IT.) C C FOR IEEE-ARITHMETIC MACHINES (BINARY STANDARD), ONE OF THE FIRST C TWO SETS OF CONSTANTS BELOW SHOULD BE APPROPRIATE. C C WHERE POSSIBLE, DECIMAL, OCTAL OR HEXADECIMAL CONSTANTS ARE USED C TO SPECIFY THE CONSTANTS EXACTLY. SOMETIMES THIS REQUIRES USING C EQUIVALENT INTEGER ARRAYS. IF YOUR COMPILER USES HALF-WORD C INTEGERS BY DEFAULT (SOMETIMES CALLED INTEGER*2), YOU MAY NEED TO C CHANGE INTEGER TO INTEGER*4 OR OTHERWISE INSTRUCT YOUR COMPILER C TO USE FULL-WORD INTEGERS IN THE NEXT 5 DECLARATIONS. C INTEGER SMALL(4) INTEGER LARGE(4) INTEGER RIGHT(4) INTEGER DIVER(4) INTEGER LOG10(4) INTEGER SC 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 IEEE ARITHMETIC MACHINES, SUCH AS THE AT&T C 3B SERIES AND MOTOROLA 68000 BASED MACHINES (E.G. SUN 3 AND AT&T C PC 7300), IN WHICH THE MOST SIGNIFICANT BYTE IS STORED FIRST. C DATA SMALL(1),SMALL(2) / 1048576, 0 / DATA LARGE(1),LARGE(2) / 2146435071, -1 / DATA RIGHT(1),RIGHT(2) / 1017118720, 0 / DATA DIVER(1),DIVER(2) / 1018167296, 0 / DATA LOG10(1),LOG10(2) / 1070810131, 1352628735 /, SC/987/ C C MACHINE CONSTANTS FOR IEEE ARITHMETIC MACHINES AND 8087-BASED C MICROS, SUCH AS THE IBM PC AND AT&T 6300, IN WHICH THE LEAST C SIGNIFICANT BYTE IS STORED FIRST. C C DATA SMALL(1),SMALL(2) / 0, 1048576 / C DATA LARGE(1),LARGE(2) / -1, 2146435071 / C DATA RIGHT(1),RIGHT(2) / 0, 1017118720 / C DATA DIVER(1),DIVER(2) / 0, 1018167296 / C DATA LOG10(1),LOG10(2) / 1352628735, 1070810131 /, SC/987/ C C MACHINE CONSTANTS FOR AMDAHL MACHINES. C C DATA SMALL(1),SMALL(2) / 1048576, 0 / C DATA LARGE(1),LARGE(2) / 2147483647, -1 / C DATA RIGHT(1),RIGHT(2) / 856686592, 0 / C DATA DIVER(1),DIVER(2) / 873463808, 0 / C DATA LOG10(1),LOG10(2) / 1091781651, 1352628735 /, SC/987/ 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 /, SC/987/ 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 /, SC/987/ 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 /, SC/987/ C C MACHINE CONSTANTS FOR FTN4 ON THE CDC 6000/7000 SERIES. C C DATA SMALL(1) / 00564000000000000000B / C DATA SMALL(2) / 00000000000000000000B / C C DATA LARGE(1) / 37757777777777777777B / C DATA LARGE(2) / 37157777777777777774B / C C DATA RIGHT(1) / 15624000000000000000B / C DATA RIGHT(2) / 00000000000000000000B / C C DATA DIVER(1) / 15634000000000000000B / C DATA DIVER(2) / 00000000000000000000B / C C DATA LOG10(1) / 17164642023241175717B / C DATA LOG10(2) / 16367571421742254654B /, SC/987/ C C MACHINE CONSTANTS FOR FTN5 ON THE CDC 6000/7000 SERIES. C C DATA SMALL(1) / O"00564000000000000000" / C DATA SMALL(2) / O"00000000000000000000" / C C DATA LARGE(1) / O"37757777777777777777" / C DATA LARGE(2) / O"37157777777777777774" / C C DATA RIGHT(1) / O"15624000000000000000" / C DATA RIGHT(2) / O"00000000000000000000" / C C DATA DIVER(1) / O"15634000000000000000" / C DATA DIVER(2) / O"00000000000000000000" / C C DATA LOG10(1) / O"17164642023241175717" / C DATA LOG10(2) / O"16367571421742254654" /, SC/987/ C C MACHINE CONSTANTS FOR CONVEX C-1 C C DATA SMALL(1),SMALL(2) / '00100000'X, '00000000'X / C DATA LARGE(1),LARGE(2) / '7FFFFFFF'X, 'FFFFFFFF'X / C DATA RIGHT(1),RIGHT(2) / '3CC00000'X, '00000000'X / C DATA DIVER(1),DIVER(2) / '3CD00000'X, '00000000'X / C DATA LOG10(1),LOG10(2) / '3FF34413'X, '509F79FF'X /, SC/987/ C C MACHINE CONSTANTS FOR THE CRAY 1, XMP, 2, AND 3. C C DATA SMALL(1) / 201354000000000000000B / C DATA SMALL(2) / 000000000000000000000B / C C DATA LARGE(1) / 577767777777777777777B / C DATA LARGE(2) / 000007777777777777776B / C C DATA RIGHT(1) / 376434000000000000000B / C DATA RIGHT(2) / 000000000000000000000B / C C DATA DIVER(1) / 376444000000000000000B / C DATA DIVER(2) / 000000000000000000000B / C C DATA LOG10(1) / 377774642023241175717B / C DATA LOG10(2) / 000007571421742254654B /, SC/987/ C C MACHINE CONSTANTS FOR THE DATA GENERAL ECLIPSE S/200 C C NOTE - IT MAY BE APPROPRIATE TO INCLUDE THE FOLLOWING LINE - 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/, SC/987/ C C MACHINE CONSTANTS FOR THE HARRIS SLASH 6 AND SLASH 7 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 /, SC/987/ C C MACHINE CONSTANTS FOR THE HONEYWELL DPS 8/70 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 /, SC/987/ C C MACHINE CONSTANTS FOR THE IBM 360/370 SERIES, C THE XEROX SIGMA 5/7/9 AND THE SEL SYSTEMS 85/86. 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 /, SC/987/ C C MACHINE CONSTANTS FOR THE INTERDATA 8/32 C WITH THE UNIX SYSTEM FORTRAN 77 COMPILER. C C FOR THE INTERDATA FORTRAN VII COMPILER REPLACE C THE Z'S SPECIFYING HEX CONSTANTS WITH Y'S. C C DATA SMALL(1),SMALL(2) / Z'00100000', Z'00000000' / C DATA LARGE(1),LARGE(2) / Z'7EFFFFFF', Z'FFFFFFFF' / C DATA RIGHT(1),RIGHT(2) / Z'33100000', Z'00000000' / C DATA DIVER(1),DIVER(2) / Z'34100000', Z'00000000' / C DATA LOG10(1),LOG10(2) / Z'41134413', Z'509F79FF' /, SC/987/ 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 /, SC/987/ 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, "047674776746 /, SC/987/ C C MACHINE CONSTANTS FOR PDP-11 FORTRANS 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 /, SC/987/ 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 /, SC/987/ C C MACHINE CONSTANTS FOR PDP-11 FORTRANS 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 /, SC/987/ 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 /, SC/987/ C C MACHINE CONSTANTS FOR THE PRIME 50 SERIES SYSTEMS C WITH 32-BIT INTEGERS AND 64V MODE INSTRUCTIONS, C SUPPLIED BY IGOR BRAY. C C DATA SMALL(1),SMALL(2) / :10000000000, :00000100001 / C DATA LARGE(1),LARGE(2) / :17777777777, :37777677775 / C DATA RIGHT(1),RIGHT(2) / :10000000000, :00000000122 / C DATA DIVER(1),DIVER(2) / :10000000000, :00000000123 / C DATA LOG10(1),LOG10(2) / :11504046501, :07674600177 /, SC/987/ C C MACHINE CONSTANTS FOR THE SEQUENT BALANCE 8000 C C DATA SMALL(1),SMALL(2) / $00000000, $00100000 / C DATA LARGE(1),LARGE(2) / $FFFFFFFF, $7FEFFFFF / C DATA RIGHT(1),RIGHT(2) / $00000000, $3CA00000 / C DATA DIVER(1),DIVER(2) / $00000000, $3CB00000 / C DATA LOG10(1),LOG10(2) / $509F79FF, $3FD34413 /, SC/987/ C C MACHINE CONSTANTS FOR THE UNIVAC 1100 SERIES. 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 /, SC/987/ C C MACHINE CONSTANTS FOR THE VAX UNIX F77 COMPILER 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, -805796613 /, SC/987/ C C MACHINE CONSTANTS FOR THE VAX-11 WITH C FORTRAN IV-PLUS COMPILER 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, ZCFF884FB /, SC/987/ C C MACHINE CONSTANTS FOR VAX/VMS VERSION 2.2 C C DATA SMALL(1),SMALL(2) / '80'X, '0'X / C DATA LARGE(1),LARGE(2) / 'FFFF7FFF'X, 'FFFFFFFF'X / C DATA RIGHT(1),RIGHT(2) / '2480'X, '0'X / C DATA DIVER(1),DIVER(2) / '2500'X, '0'X / C DATA LOG10(1),LOG10(2) / '209A3F9A'X, 'CFF884FB'X /, SC/987/ C C *** ISSUE STOP 779 IF ALL DATA STATEMENTS ARE COMMENTED... IF (SC .NE. 987) STOP 779 IF (I .LT. 1 .OR. I .GT. 5) GOTO 999 D1MACH = DMACH(I) RETURN 999 WRITE(*,1999) I 1999 FORMAT(' D1MACH - I OUT OF BOUNDS',I10) STOP END SUBROUTINE DCHTRI(NUMFUN,MDIV,VER,NUMTRI,MINPTS,MAXPTS,EPSABS, + EPSREL,LENVER,NW,RESTAR,MAXSUB,MINSUB,IFAIL) C***BEGIN PROLOGUE DCHTRI C***REFER TO DCUTRI C***PURPOSE DCHTRI checks the validity of the C input parameters to DCUTRI. C***DESCRIPTION C DCHTRI computes MAXSUB, MINSUB and IFAIL as C functions of the input parameters to DCUTRI, C and checks the validity of the input parameters to DCUTRI. C C ON ENTRY C C NUMFUN Integer. C Number of components of the integral. C MDIV Integer. C MDIV is the number of triangles that are divided in C each subdivision step in DADTRI. C MDIV is chosen default to 1. C For efficient execution on parallel computers C with NPROC processors MDIV should be set equal to C the smallest integer such that MOD(4*MDIV,NPROC) = 0. C VER Real array of dimension (2,3,LENVER). C VER(1,K,L) and VER(2,K,L) are the x and y coordinates C respectively of vertex K in triangle L. C On entry VER(*,*,L) must contain the vertices of the C NUMTRI user specified triangles for L = 1,2,...,NUMTRI. C NUMTRI Integer. C The number of triangles. C MINPTS Integer. C Minimum number of function evaluations. C MAXPTS Integer. C Maximum number of function evaluations. C The number of function evaluations over each subregion C is 37. C MAXPTS >= 37*NUMTRI and MAXPTS >= MINPTS C C EPSABS Real. C Requested absolute error. C EPSREL Real. C Requested relative error. C LENVER Integer. C Defines the length of the array VER. C C We let C MAXSUB denote the maximum allowed number of subregions C for the given values of MAXPTS. C MAXSUB = 3*((MAXPTS-37*NUMTRI)/(4*37)) + NUMTRI C LENVER should then be greater or equal to MAXSUB. C C NW Integer. C Defines the length of the working array WORK. C C We let C MAXSUB denote the maximum allowed number of subregions C for the given values of MAXPTS. C MAXSUB = 3*((MAXPTS-37*NUMTRI)/(4*37)) + NUMTRI C NW should then be greater or equal to C MAXSUB*(2*NUMFUN+1) + MAX(32*MDIV,8*NUMTRI)*NUMFUN + 1 C C RESTAR Integer. C If RESTAR = 0, this is the first attempt to compute C the integral. C If RESTAR = 1, C then we restart a previous attempt. C In this case the only parameters for DCUTRI that may C be changed (with respect to the previous call of DCUTRI) C are MINPTS, MAXPTS, EPSABS, EPSREL and RESTAR. C C C ON RETURN C C MAXSUB Integer. C The maximum allowed number of subregions for the C given values of MAXPTS. C MINSUB Integer. C The minimum allowed number of subregions for the given C values of MINPTS. C IFAIL Integer. C IFAIL = 0 for normal exit. C IFAIL = 2 if NUMFUN is less than 1. C IFAIL = 3 if area of one of the initially given C triangles is zero. C IFAIL = 4 if MAXPTS is less than 37*NUMTRI. C IFAIL = 5 if MAXPTS is less than MINPTS. C IFAIL = 6 if EPSABS <= 0 and EPSREL <= 0. C IFAIL = 7 if LENVER is less than MAXSUB. C IFAIL = 8 if NW is less than MAXSUB*(2*NUMFUN+1) + C NUMFUN*MAX(32*MDIV,8*NUMTRI) + 1. C IFAIL = 9 if illegal RESTAR. C C***ROUTINES CALLED-NONE C***END PROLOGUE DCHTRI C C Global variables. C INTEGER NUMFUN,MDIV,NUMTRI,MINPTS,MAXPTS,NW,MAXSUB,MINSUB,RESTAR, + LENVER,IFAIL DOUBLE PRECISION VER(2,3,NUMTRI),EPSABS,EPSREL C C Local variables. C INTEGER LIMIT,J DOUBLE PRECISION AREA C C***FIRST EXECUTABLE STATEMENT DCHTRI C IFAIL = 0 C C We compute MAXSUB and MINSUB. C MAXSUB = 3* ((MAXPTS-37*NUMTRI)/ (4*37)) + NUMTRI MINSUB = 3* ((MINPTS-37*NUMTRI)/ (4*37)) + NUMTRI IF (MOD((MINPTS-37*NUMTRI),4*37).GT.0) THEN MINSUB = MINSUB + 3 END IF MINSUB = MAX(NUMTRI,MINSUB) C C Check on positive NUMFUN. C IF (NUMFUN.LT.1) THEN IFAIL = 2 RETURN END IF C C Check on legal area of the region of integration. C DO 10 J = 1,NUMTRI AREA = ABS(VER(1,1,J)*VER(2,2,J)-VER(1,2,J)*VER(2,1,J)- + VER(1,1,J)*VER(2,3,J)+VER(1,3,J)*VER(2,1,J)+ + VER(1,2,J)*VER(2,3,J)-VER(1,3,J)*VER(2,2,J))*0.5 IF (AREA.EQ.0) THEN IFAIL = 3 RETURN END IF 10 CONTINUE C C Check on MAXPTS >= 37*NUMTRI. C IF (MAXPTS.LT.37*NUMTRI) THEN IFAIL = 4 RETURN END IF C C Check on MAXPTS >= MINPTS. C IF (MAXPTS.LT.MINPTS) THEN IFAIL = 5 RETURN END IF C C Check on legal accuracy requests. C IF (EPSABS.LE.0 .AND. EPSREL.LE.0) THEN IFAIL = 6 RETURN END IF C C Check on big enough LENVER. C IF (LENVER.LT.MAXSUB) THEN IFAIL = 7 RETURN END IF C C Check on big enough NW. C LIMIT = MAXSUB* (2*NUMFUN+1) + MAX(32*MDIV,8*NUMTRI)*NUMFUN + 1 IF (NW.LT.LIMIT) THEN IFAIL = 8 RETURN END IF C C Check on legal RESTAR. C IF (RESTAR.NE.0 .AND. RESTAR.NE.1) THEN IFAIL = 9 RETURN END IF RETURN C C***END DCHTRI C END SUBROUTINE DCUTRI(FUNSUB,NUMFUN,VER,NUMTRI,MINPTS,MAXPTS,EPSABS, + EPSREL,LENVER,NW,RESTAR,RESULT,ABSERR,NEVAL, + IFAIL,WORK,IWORK) C***BEGIN PROLOGUE DCUTRI C***DATE WRITTEN 891023 (YYMMDD) C***REVISION DATE 910214 (YYMMDD) C***CATEGORY NO. H2B2A1/H2B2A2 C***KEYWORDS QUADRATURE, TWO DIMENSIONAL, ADAPTIVE, CUBATURE C***AUTHOR C C Jarle Berntsen, Institute of Marine Research, C Nordnesparken 2, Postboks 1870, C N-5024 Bergen-Nordnes, Norway C Phone.. 47-5-238461 C Email.. jarleb@imr.no C C Terje O. Espelid, Department of Informatics, C University of Bergen, Thormohlens gt. 55, C N-5008 Bergen, Norway C Phone.. 47-5-544180 C Email.. terje@eik.ii.uib.no C C***PURPOSE To compute the two-dimensional integral over a region C consisting of NUMTRI triangles. C***DESCRIPTION C C The routine calculates an approximation to a given C vector of definite integrals C C I I (F ,F ,...,F ) DX(2)DX(1), C 1 2 NUMFUN C C where the region of integration is a selection of C NUMTRI triangles and C where F = F (X(1),X(2)), J = 1,2,...,NUMFUN. C J J C C C A globally adaptive strategy is applied in order to C compute approximations RESULT(K) C hopefully satisfying for each component of I the following C claim for accuracy: C ABS(I(K)-RESULT(K)).LE.MAX(EPSABS,EPSREL*ABS(I(K))) C C DCUTRI is a driver for the integration routine C DADTRI. C C DADTRI repeatedly C subdivides the triangles with greatest estimated errors C and estimates the integrals and the C errors over the new subtriangles C until either the error request C is met or MAXPTS function evaluations have been used. C C A 37 point integration rule C with all evaluation points inside the triangle C is applied. The rule has polynomial degree 13. C C If the values of the input parameters EPSABS C or EPSREL are selected great enough, C an integration rule is applied over each triangle and C the results are added up to give the approximations C RESULT(K). No further subdivision C of the triangles will then be applied. C C When DCUTRI computes estimates to a vector of C integrals, all components of the vector are given C the same treatment. That is, I(F ) and I(F ) for C J K C J not equal to K, are estimated with the same C subdivision of the region of integration. C For integrals with enough similarity, we may save C time by applying DCUTRI to all integrands in one call. C For integrals that varies continuously as functions of C some parameter, the estimates produced by DCUTRI will C also vary continuously when the same subdivision is C applied to all components. This will generally not be C the case when the different components are given C separate treatment. C C On the other hand this feature should be used with C caution when the different components of the integrals C require clearly different subdivisions. C C ON ENTRY C C FUNSUB Externally declared subroutine for computing C all components of the integrand at the given C evaluation point. C It must have parameters (X,NUMFUN,FUNVLS) C Input parameters: C X(1) The x-coordinate of the evaluation point. C X(2) The y-coordinate of the evaluation point. C NUMFUN Integer that defines the number of C components of I. C Output parameter: C FUNVLS Real array of dimension NUMFUN C that defines NUMFUN components of the integrand. C C NUMFUN Integer. C Number of components of the integral. C VER Real array of dimension (2,3,LENVER). C VER(1,K,L) and VER(2,K,L) are the x and y coordinates C respectively of vertex K in triangle L. C On entry VER(*,*,L) must contain the vertices of the C NUMTRI user specified triangles for L = 1,2,...,NUMTRI. C VER may be changed on exit. C C NUMTRI Integer. C The number of triangles. C MINPTS Integer. C Minimum number of function evaluations. C MAXPTS Integer. C Maximum number of function evaluations. C The number of function evaluations over each subregion C is 37. C C MAXPTS >= 37*NUMTRI and MAXPTS >= MINPTS C C EPSABS Real. C Requested absolute error. C EPSREL Real. C Requested relative error. C LENVER Integer. C Defines the length of the array VER. C C We let C MAXSUB denote the maximum allowed number of subregions C for the given value of MAXPTS. C MAXSUB = 3*((MAXPTS-37*NUMTRI)/(4*37)) + NUMTRI C LENVER should be greater or equal to MAXSUB. C C NW Integer. C Defines the length of the working array WORK. C C If LENVER is chosen correctly as a function of MAXPTS C and NUMTRI and NW is selected to be greater or equal to C LENVER*(2*NUMFUN+1) + MAX(32*MDIV,8*NUMTRI)*NUMFUN + 1 , C then the size of the working storage will be great enough. C MDIV is the number of triangles that are divided in C each subdivision step. C MDIV is default set internally in DCUTRI equal to 1. C For efficient execution on parallel computers C with NPROC processors MDIV should be set equal to C the smallest integer such that MOD(4*MDIV,NPROC) = 0. C C RESTAR Integer. C If RESTAR = 0, this is the first attempt to compute C the integral. C If RESTAR = 1, C then we restart a previous attempt. C In this case the only parameters for DCUTRI that may C be changed (with respect to the previous call of DCUTRI) C are MINPTS, MAXPTS, EPSABS, EPSREL and RESTAR. C Note: If MAXPTS was too small in the previous call, C then we get a second chance to continue the computations C with MAXPTS increased. C LENVER can not be changed, therefore the new value of MAXPTS C should not be chosen greater than the value of LENVER C allows us to do. C WORK Real array of dimension NW. C Used as working storage. C IWORK Integer array of dimension LENVER + MDIV. C Used as working storage. C C ON RETURN C C RESULT Real array of dimension NUMFUN. C Approximations to all components of the integral. C ABSERR Real array of dimension NUMFUN. C Estimates of absolute errors. C NEVAL Integer. C Number of function evaluations used by DCUTRI. C IFAIL Integer. C IFAIL = 0 for normal exit. C C ABSERR(K) <= EPSABS or C ABSERR(K) <= ABS(RESULT(K))*EPSREL with MAXPTS or less C function evaluations for all values of K, C 1 <= K <= NUMFUN . C C IFAIL = 1 if MAXPTS was too small for DCUTRI C to obtain the required accuracy. In this case DCUTRI C returns values of RESULT with estimated absolute C errors ABSERR. C IFAIL = 2 if NUMFUN is less than 1. C IFAIL = 3 if area of one of the initially given C triangles is zero. C IFAIL = 4 if MAXPTS is less than 37*NUMTRI. C IFAIL = 5 if MAXPTS is less than MINPTS. C IFAIL = 6 if EPSABS <= 0 and EPSREL <= 0. C IFAIL = 7 if LENVER is less than MAXSUB. C IFAIL = 8 if NW is less than MAXSUB*(2*NUMFUN+1) + C NUMFUN*MAX(32*MDIV,8*NUMTRI) + 1. C IFAIL = 9 if unlegal RESTAR. C VER Real array of dimension (2,3,LENVER). C On exit VER C contains the vertices of the triangles used to produce C the approximations to the integrals. C WORK Real array of dimension NW. C Used as working storage. C WORK(NW) = NSUB, the number of subregions in the data C structure. C Let WRKSUB=(NW-1-NUMFUN*MAX(32*MDIV,8*NUMTRI))/ C (2*NUMFUN+1). C WORK(1),...,WORK(NUMFUN*WRKSUB) contain C the estimated components of the integrals over the C subregions. C WORK(NUMFUN*WRKSUB+1),...,WORK(2*NUMFUN*WRKSUB) contain C the estimated errors over the subregions. C WORK(2*NUMFUN*WRKSUB+1),...,WORK(2*NUMFUN* C WRKSUB+WRKSUB) contain the greatest errors C in each subregion. C WORK((2*NUMFUN+1)*WRKSUB),...,WORK(NW - 1) is used as C temporary storage in DADTRI. C IWORK Integer array of dimension LENVER + MDIV. C Used as working storage. C C SAMPLE PROGRAM C The following program will approximate the integral of C exp(x*x+y*y) C over the triangle (0.,0.),(2.,0.),(0.,2.) and print out the C values of the estimated integral, the estimated error, the number C of function evaluations, and IFAIL. C PROGRAM DTEST1 C INTEGER NUMFUN,NW,MDIV,MINPTS,LENVER,NUMTRI C PARAMETER (NUMFUN=1,MDIV=1,NUMTRI=1,MINPTS=37) C PARAMETER ( LENVER = 100 ) C C We need to choose LENVER such that: C LENVER >= 3*((MAXPTS-37*NUMTRI)/(4*37)) + NUMTRI C This simply means that we have enough space to achieve MAXPTS C function evaluations. By choosing LENVER bigger than C this lower bound we may later change MAXPTS and RESTAR and restart C the computations from the point where the last run stopped. C The reason for stopping may have been that MAXPTS was too small C to achieve the requested error. C With our choice LENVER = 100 any value of MAXPTS C between 37 and 5068 will be accepted by the code. Choosing C MAXPTS = 2000 allows us to restart with a greater value C if necessary. C C PARAMETER (NW = LENVER*(2*NUMFUN+1) + C + MAX(32*MDIV,8*NUMTRI)*NUMFUN + 1) C DOUBLE PRECISION VER(2,3,LENVER),RESULT(NUMFUN), C + ABSERR(NUMFUN),WORK(NW),EPSABS,EPSREL C INTEGER RESTAR,NEVAL,IFAIL,MDIV,IWORK(LENVER+MDIV),MAXPTS C EXTERNAL F C VER(1,1,1) = 0.D0 C VER(1,2,1) = 2.D0 C VER(1,3,1) = 0.D0 C VER(2,1,1) = 0.D0 C VER(2,2,1) = 0.D0 C VER(2,3,1) = 2.D0 C EPSABS = 0.D0 C EPSREL = 1.D-5 C RESTAR = 0 C MAXPTS = 2000 C CALL DCUTRI(F,NUMFUN,VER,NUMTRI,MINPTS,MAXPTS,EPSABS, C + EPSREL,LENVER,NW,RESTAR,RESULT,ABSERR,NEVAL, C + IFAIL,WORK,IWORK) C WRITE(*,*)'RESULT=',RESULT(1) C WRITE(*,*)'ERROR =',ABSERR(1) C WRITE(*,*)'NEVAL =',NEVAL C WRITE(*,*)'IFAIL =',IFAIL C END C SUBROUTINE F(X,NUMFUN,FUNVLS) C INTEGER NUMFUN C DOUBLE PRECISION X(2),FUNVLS(NUMFUN) C FUNVLS(1) = EXP(X(1)*X(1)+X(2)*X(2)) C RETURN C END C C Output produced on a SUN SPARC station 1. C C RESULT= 11.181284417019 C ERROR = 4.7009048669038D-06 C NEVAL = 185 C IFAIL = 0 C C C***LONG DESCRIPTION C C The information for each triangle is contained in the C data structures VER, WORK and IWORK. C VER contains the coordinates of the triangles. C When passed on to DADTRI, WORK is split into four C arrays VALUES, ERRORS, GREATE and WORK2. C VALUES contains the estimated values of the integrals. C ERRORS contains the estimated errors. C GREATE contains the greatest estimated error for each triangle. C The data structures for the triangles are in DADTRI organized C as a heap, and the size of GREATE(I) defines the position of C triangle I in the heap. The heap is maintained by the program C DTRTRI and we use a partially ordered list of pointers, saved C in IWORK. When passed on to DADTRI, IWORK is split into two C arrays LIST and VACANT. LIST is a partially ordered list C such that GREATE(LIST(1)) is the maximum error estimate for C all sub-triangles in our datastructure. VACANT is a work space C needed in the updating process of the list. C C The subroutine DADTRI is written for efficient execution on shared C memory parallel computer. On a computer with NPROC processors we will C in each subdivision step divide MDIV triangles, where MDIV is C chosen such that MOD(4*MDIV,NPROC) = 0, in totally 4*MDIV new triangles. C Each processor will then compute estimates of the integrals and errors C over 4*MDIV/NPROC triangles in each subdivision step. C The subroutine for estimating the integral and the error over C each subregion, DRLTRI, uses WORK2 as a work array. C We must make sure that each processor writes its results to C separate parts of the memory, and therefore the sizes of WORK and C WORK2 are functions of MDIV. C In order to achieve parallel processing of triangles, compiler C directives should be placed in front of the DO 20 and the DO 200 C loops in DADTRI on machines like Alliant and CRAY. C C***REFERENCES C J.Berntsen and T.O.Espelid, An Algorithm for Adaptive C Cubature Over a Collection of Triangles, to be published. C C***ROUTINES CALLED DCHTRI,DADTRI C***END PROLOGUE DCUTRI C C Parameters C C MDIV Integer. C MDIV is the number of triangles that are divided in C each subdivision step in DADTRI. C MDIV is chosen default to 1. C For efficient execution on parallel computers C with NPROC processors MDIV should be set equal to C the smallest integer such that MOD(4*MDIV,NPROC) = 0. INTEGER MDIV PARAMETER (MDIV=1) C C Global variables. C EXTERNAL FUNSUB INTEGER NUMFUN,NUMTRI,MINPTS,MAXPTS,LENVER,NW,RESTAR,NEVAL,IFAIL, + IWORK(LENVER+MDIV) DOUBLE PRECISION VER(2,3,LENVER),EPSABS,EPSREL,RESULT(NUMFUN), + ABSERR(NUMFUN),WORK(NW) C C Local variables. C C MAXSUB Integer. C The maximum allowed number of subregions C for the given values of MAXPTS. C MINSUB Integer. C The minimum allowed number of subregions for the given C values of MINPTS. C WRKSUB Integer. C Determines the length of the main work arrays. C INTEGER MAXSUB,MINSUB,NSUB,LENW INTEGER WRKSUB,I1,I2,I3,I4 C C***FIRST EXECUTABLE STATEMENT DCUTRI C C Compute MAXSUB and MINSUB, C and check the input parameters. C CALL DCHTRI(NUMFUN,MDIV,VER,NUMTRI,MINPTS,MAXPTS,EPSABS,EPSREL, + LENVER,NW,RESTAR,MAXSUB,MINSUB,IFAIL) WRKSUB = (NW-1-NUMFUN*MAX(32*MDIV,8*NUMTRI))/ (2*NUMFUN+1) IF (IFAIL.NE.0) THEN RETURN END IF C C Split up the work space. C I1 = 1 I2 = I1 + WRKSUB*NUMFUN I3 = I2 + WRKSUB*NUMFUN I4 = I3 + WRKSUB C C On restart runs the number of subregions from the C previous call is assigned to NSUB. C IF (RESTAR.EQ.1) THEN NSUB = WORK(NW) END IF C C Compute the size of the temporary work space needed in DADTRI. C LENW = MAX(32*MDIV,8*NUMTRI)*NUMFUN CALL DADTRI(NUMFUN,MDIV,VER,NUMTRI,MINSUB,MAXSUB,FUNSUB,EPSABS, + EPSREL,LENVER,RESTAR,LENW,RESULT,ABSERR,NEVAL,NSUB, + IFAIL,WORK(I1),WORK(I2),WORK(I3),WORK(I4), + IWORK(1),IWORK(1+LENVER)) WORK(NW) = NSUB RETURN C C***END DCUTRI C END SUBROUTINE DRLTRI(VER,NUMFUN,FUNSUB,NULL,BASVAL,RGNERR,GREATE) C***BEGIN PROLOGUE DRLTRI C***REFER TO DCUTRI C***PURPOSE To compute basic integration rule values and C corresponding error estimates. C***DESCRIPTION DRLTRI computes basic integration rule values C for a vector of integrands over a triangular region. C DRLTRI also computes estimates for the errors by C using several null rule approximations. C ON ENTRY C C VER Real array of dimension (2,3). C The coordinates of the vertices of the triangle. C NUMFUN Integer. C Number of components of the vector integrand. C FUNSUB Externally declared subroutine for computing C all components of the integrand at the given C evaluation point. C It must have parameters (X,NUMFUN,FUNVLS) C Input parameters: C X(1) The x-coordinate of the evaluation point. C X(2) The y-coordinate of the evaluation point. C NUMFUN Integer that defines the number of C components of I. C Output parameter: C FUNVLS Real array of dimension NUMFUN C that defines NUMFUN components of the integrand. C C NULL Real array of dimension (NUMFUN,8). C A work array. C C ON RETURN C C BASVAL Real array of dimension NUMFUN. C The values for the basic rule for each component C of the integrand. C RGNERR Real array of dimension NUMFUN. C The error estimates for each component of the integrand. C GREATE Real. C The greatest error component of the integrand. C C***REFERENCES Berntsen,J. and Espelid,T.O., Degree 13 Symmetric C Quadrature Rules for the Triangle, Reports in C Informatics 44, Dept. of Informatics, University of C Bergen, 1990. C***ROUTINES CALLED D1MACH,FUNSUB C***END PROLOGUE DRLTRI C C Global variables. C EXTERNAL FUNSUB INTEGER NUMFUN DOUBLE PRECISION VER(2,3),BASVAL(NUMFUN),RGNERR(NUMFUN), + NULL(NUMFUN,8),GREATE,NOISE,D1MACH,TRES C C Local variables. C C WTLENG The number of weights of the integration rule. C G Real array of dimension (2,WTLENG). C The homogeneous coordinates for the generators of C the evaluation points. C The integration rule is using symmetric evaluation C points and has the structure (1,6,3). That is, C 1 point of multiplicity 1, C 6 sets of points of multiplicity 3 and C 3 sets of points of multiplicity 6. C This gives totally 37 evaluation points. C In order to reduce the number of loops in DRLTRI, C the 3 loops for the sets of multiplicity 6 are split C into 6 loops and added to the loops for the sets of C multiplicity 3. C The number of weights we have to give with C this splitting is WTLENG = 13. C C W Real array of dimension (9,WTLENG). C The weights of the basic rule and the null rules. C W(1,1),...,W(1,WTLENG) are weights for the basic rule. C W(I,1),...,W(I,WTLENG) for I>1 are null rule weights. C INTEGER WTLENG PARAMETER (WTLENG=13) DOUBLE PRECISION AREA,X(2,3),Z1,Z2,Z3,R1,R2,R3,R,DEG7,DEG5,DEG3, + DEG1,G(2,13),W(9,13) INTEGER I,J,K,L C C The abscissas are given in homogeneous coordinates. C DATA (G(1,I),I=1,13)/0.333333333333333333333333333333D0, + 0.950275662924105565450352089520D0, + 0.171614914923835347556304795551D0, + 0.539412243677190440263092985511D0, + 0.772160036676532561750285570113D0, + 0.009085399949835353883572964740D0, + 0.062277290305886993497083640527D0, + 0.022076289653624405142446876931D0, + 0.018620522802520968955913511549D0, + 0.096506481292159228736516560903D0, + 0.851306504174348550389457672223D0, + 0.689441970728591295496647976487D0, + 0.635867859433872768286976979827D0/ DATA (G(2,I),I=1,13)/0.333333333333333333333333333333D0, + 0.024862168537947217274823955239D0, + 0.414192542538082326221847602214D0, + 0.230293878161404779868453507244D0, + 0.113919981661733719124857214943D0, + 0.495457300025082323058213517632D0, + 0.468861354847056503251458179727D0, + 0.851306504174348550389457672223D0, + 0.689441970728591295496647976487D0, + 0.635867859433872768286976979827D0, + 0.022076289653624405142446876931D0, + 0.018620522802520968955913511549D0, + 0.096506481292159228736516560903D0/ C C Weights of the degree 13 quadrature rule. C DATA (W(1,I),I=1,13)/0.051739766065744133555179145422D0, + 0.008007799555564801597804123460D0, + 0.046868898981821644823226732071D0, + 0.046590940183976487960361770070D0, + 0.031016943313796381407646220131D0, + 0.010791612736631273623178240136D0, + 0.032195534242431618819414482205D0, + 0.015445834210701583817692900053D0, + 0.017822989923178661888748319485D0, + 0.037038683681384627918546472190D0, + 0.015445834210701583817692900053D0, + 0.017822989923178661888748319485D0, + 0.037038683681384627918546472190D0/ C C Weights of the first null rule of degree 7. C DATA (W(2,I),I=1,13)/-0.077738051051462052051304462750D0, + 0.001640389740236881582083124927D0, + 0.078124083459915167386776552733D0, + -0.030706528522391137165581298102D0, + 0.010246307817678312345028512621D0, + 0.012586300774453821540476193059D0, + -0.043630506151410607808929481439D0, + -0.004567055157220063810223671248D0, + 0.003393373439889186878847613140D0, + 0.000000000000000000000000000000D0, + -0.004567055157220063810223671248D0, + 0.003393373439889186878847613140D0, + 0.000000000000000000000000000000D0/ C C Weights of the second null rule of degree 7. C DATA (W(3,I),I=1,13)/-0.064293709240668260928898888457D0, + 0.003134665264639380635175608661D0, + 0.007822550509742830478456728602D0, + 0.048653051907689492781049400973D0, + 0.032883327334384971735434067029D0, + -0.017019508374229390108580829589D0, + 0.025973557893399824586684707198D0, + -0.010716753326806275930657622320D0, + 0.018315629578968063765722278290D0, + -0.047607080313197299401024682666D0, + -0.010716753326806275930657622320D0, + 0.018315629578968063765722278290D0, + -0.047607080313197299401024682666D0/ C C Weights of the first degree 5 null rule. C C DATA (W(4,I),I=1,13)/0.021363205584741860993131879186D0, + 0.022716410154120323440432428315D0, + -0.026366191271182090678117381002D0, + 0.029627021479068212693155637482D0, + 0.004782834546596399307634111034D0, + 0.004178667433984132052378990240D0, + -0.065398996748953861618846710897D0, + -0.033589813176131630980793760168D0, + 0.033018320112481615757912576257D0, + 0.012241086002709814125707333127D0, + -0.033589813176131630980793760168D0, + 0.033018320112481615757912576257D0, + 0.012241086002709814125707333127D0/ C C Weights of the second degree 5 null rule. C DATA (W(5,I),I=1,13)/-0.046058756832790538620830792345D0, + 0.005284159186732627192774759959D0, + 0.009325799301158899112648198129D0, + -0.006101110360950124560783393745D0, + -0.056223328794664871336486737231D0, + -0.062516479198185693171971930698D0, + 0.022428226812039547178810743269D0, + -0.000026014926110604563130107142D0, + 0.032882099937471182365626663487D0, + 0.018721740987705986426812755881D0, + -0.000026014926110604563130107142D0, + 0.032882099937471182365626663487D0, + 0.018721740987705986426812755881D0/ C C Weights of first degree 3 null rule. C DATA (W(6,I),I=1,13)/0.080867117677405246540283712799D0, + -0.033915806661511608094988607349D0, + 0.014813362053697845461526433401D0, + 0.001442315416337389214102507204D0, + -0.024309696484708683486455879210D0, + -0.005135085639122398522835391664D0, + -0.034649417896235909885490654650D0, + 0.035748423431577326597742956780D0, + 0.024548155266816447583155562333D0, + -0.032897267038856299280541675107D0, + 0.035748423431577326597742956780D0, + 0.024548155266816447583155562333D0, + -0.032897267038856299280541675107D0/ C C Weights of second degree 3 null rule. C DATA (W(7,I),I=1,13)/-0.038457863913548248582247346193D0, + -0.055143631258696406147982448269D0, + -0.021536994314510083845999131455D0, + 0.001547467894857008228010564582D0, + 0.057409361764652373776043522086D0, + -0.040636938884669694118908764512D0, + -0.020801144746964801777584428369D0, + 0.019490770404993674256256421103D0, + 0.002606109985826399625043764771D0, + 0.023893703367437102825618048130D0, + 0.019490770404993674256256421103D0, + 0.002606109985826399625043764771D0, + 0.023893703367437102825618048130D0/ C C Weights of first degree 1 null rule. C DATA (W(8,I),I=1,13)/0.074839568911184074117081012527D0, + -0.004270103034833742737299816615D0, + 0.049352639555084484177095781183D0, + 0.048832124609719176627453278550D0, + 0.001042698696559292759051590242D0, + -0.044445273029113458906055765365D0, + -0.004670751812662861209726508477D0, + -0.015613390485814379318605247424D0, + -0.030581651696100000521074498679D0, + 0.010801113204340588798240297593D0, + -0.015613390485814379318605247424D0, + -0.030581651696100000521074498679D0, + 0.010801113204340588798240297593D0/ C C Weights of second degree 1 null rule. C DATA (W(9,I),I=1,13)/0.009373028261842556370231264134D0, + -0.074249368848508554545399978725D0, + 0.014709707700258308001897299938D0, + 0.009538502545163567494354463302D0, + -0.014268362488069444905870465047D0, + 0.040126396495352694403045023109D0, + 0.028737181842214741174950928350D0, + -0.031618075834734607275229608099D0, + 0.016879961075872039084307382161D0, + 0.010878914758683152984395046434D0, + -0.031618075834734607275229608099D0, + 0.016879961075872039084307382161D0, + 0.010878914758683152984395046434D0/ C C***FIRST EXECUTABLE STATEMENT DRLTRI C TRES = 50*D1MACH(4) C C Compute area of triangle. C AREA = ABS(VER(1,1)*VER(2,2)-VER(1,2)*VER(2,1)-VER(1,1)*VER(2,3)+ + VER(1,3)*VER(2,1)+VER(1,2)*VER(2,3)-VER(1,3)*VER(2,2))/2 C C Compute contributions from the center of the triangle. C X(1,1) = (VER(1,1)+VER(1,2)+VER(1,3))/3 X(2,1) = (VER(2,1)+VER(2,2)+VER(2,3))/3 CALL FUNSUB(X,NUMFUN,RGNERR) DO 20 J = 1,NUMFUN BASVAL(J) = W(1,1)*RGNERR(J) DO 10 K = 1,8 NULL(J,K) = W(K+1,1)*RGNERR(J) 10 CONTINUE 20 CONTINUE C C Compute the contributions from the points with C multiplicity 3. C DO 100 I = 2,WTLENG Z1 = G(1,I) Z2 = G(2,I) Z3 = 1 - Z1 - Z2 X(1,1) = Z1*VER(1,1) + Z2*VER(1,2) + Z3*VER(1,3) X(2,1) = Z1*VER(2,1) + Z2*VER(2,2) + Z3*VER(2,3) X(1,2) = Z2*VER(1,1) + Z3*VER(1,2) + Z1*VER(1,3) X(2,2) = Z2*VER(2,1) + Z3*VER(2,2) + Z1*VER(2,3) X(1,3) = Z3*VER(1,1) + Z1*VER(1,2) + Z2*VER(1,3) X(2,3) = Z3*VER(2,1) + Z1*VER(2,2) + Z2*VER(2,3) DO 90 L = 1,3 CALL FUNSUB(X(1,L),NUMFUN,RGNERR) DO 80 J = 1,NUMFUN BASVAL(J) = BASVAL(J) + W(1,I)*RGNERR(J) DO 70 K = 1,8 NULL(J,K) = NULL(J,K) + W(K+1,I)*RGNERR(J) 70 CONTINUE 80 CONTINUE 90 CONTINUE 100 CONTINUE C C Compute the errors. C GREATE = 0 DO 130 J = 1,NUMFUN DEG7 = AREA*SQRT(NULL(J,1)**2+NULL(J,2)**2) DEG5 = AREA*SQRT(NULL(J,3)**2+NULL(J,4)**2) DEG3 = AREA*SQRT(NULL(J,5)**2+NULL(J,6)**2) DEG1 = AREA*SQRT(NULL(J,7)**2+NULL(J,8)**2) IF (DEG5.NE.0) THEN R1 = DEG7/DEG5 ELSE R1 = 1 END IF IF (DEG3.NE.0) THEN R2 = DEG5/DEG3 ELSE R2 = 1 END IF IF (DEG1.NE.0) THEN R3 = DEG3/DEG1 ELSE R3 = 1 END IF R = MAX(R1,R2,R3) IF (R.GE.1) THEN RGNERR(J) = 10*MAX(DEG1,DEG3,DEG5,DEG7) ELSE IF (R.GE.0.5D0) THEN RGNERR(J) = 10*R*DEG7 ELSE RGNERR(J) = 40* (R**3)*DEG7 END IF BASVAL(J) = AREA*BASVAL(J) NOISE = ABS(BASVAL(J))*TRES C The following 6 statements added as the result of authors remark C C The following statement is included to set the error to the noise C level if the two error estimates, assumed to be the best, are both C below or on the noise level C IF ((DEG7.LE.NOISE).AND.(DEG5.LE.NOISE)) RGNERR(J)=NOISE RGNERR(J) = MAX(NOISE,RGNERR(J)) IF (RGNERR(J).GT.GREATE) THEN GREATE = RGNERR(J) END IF 130 CONTINUE C C***END DRLTRI C END C SAMPLE PROGRAM C The following program will approximate the integral of C exp(x*x+y*y) C over the triangle (0.,0.),(2.,0.),(0.,2.) and print out the C values of the estimated integral, the estimated error, the number C of function evaluations, and IFAIL. PROGRAM DTEST1 INTEGER NUMFUN,NW,MDIV,MINPTS,LENVER,NUMTRI PARAMETER (NUMFUN=1,MDIV=1,NUMTRI=1,MINPTS=37) PARAMETER ( LENVER = 100 ) C C We need to choose LENVER such that: C LENVER >= 3*((MAXPTS-37*NUMTRI)/(4*37)) + NUMTRI C This simply means that we have enough space to achieve MAXPTS C function evaluations. By choosing LENVER bigger than C this lower bound we may later change MAXPTS and RESTAR and restart C the computations from the point where the last run stopped. C The reason for stopping may have been that MAXPTS was too small C to achieve the requested error. C With our choice LENVER = 100 any value of MAXPTS C between 37 and 5068 will be accepted by the code. Choosing C MAXPTS = 2000 allows us to restart with a greater value C if necessary. C PARAMETER (NW = LENVER*(2*NUMFUN+1) + + MAX(32*MDIV,8*NUMTRI)*NUMFUN + 1) DOUBLE PRECISION VER(2,3,LENVER),RESULT(NUMFUN), + ABSERR(NUMFUN),WORK(NW),EPSABS,EPSREL INTEGER RESTAR,NEVAL,IFAIL,IWORK(LENVER+MDIV),MAXPTS EXTERNAL F VER(1,1,1) = 0.D0 VER(1,2,1) = 2.D0 VER(1,3,1) = 0.D0 VER(2,1,1) = 0.D0 VER(2,2,1) = 0.D0 VER(2,3,1) = 2.D0 EPSABS = 0.D0 EPSREL = 1.D-5 RESTAR = 0 MAXPTS = 2000 CALL DCUTRI(F,NUMFUN,VER,NUMTRI,MINPTS,MAXPTS,EPSABS, + EPSREL,LENVER,NW,RESTAR,RESULT,ABSERR,NEVAL, + IFAIL,WORK,IWORK) WRITE(*,*)'RESULT=',RESULT(1) WRITE(*,*)'ERROR =',ABSERR(1) WRITE(*,*)'NEVAL =',NEVAL WRITE(*,*)'IFAIL =',IFAIL END SUBROUTINE F(X,NUMFUN,FUNVLS) INTEGER NUMFUN DOUBLE PRECISION X(2),FUNVLS(NUMFUN) FUNVLS(1) = EXP(X(1)*X(1)+X(2)*X(2)) RETURN END C C Output produced on a SUN SPARC station 1. C C RESULT= 11.181284417019 C ERROR = 4.7009048669038D-06 C NEVAL = 185 C IFAIL = 0 C PROGRAM DTEST2 C C DTEST2 tests some of the features of DCUTRI. C DTEST2 checks that DCUTRI integrates to machine C precision all monomials of degree less or equal to 13. C DTEST2 checks that DCUTRI integrates correctly a vector C of integrands. C DTEST2 checks that the restart feature of DCUTRI works. C EXTERNAL FTESTP,FTESTQ,FTESTO C C DCUTRI TEST PROGRAM C INTEGER NW,IWORK(501) PARAMETER (NW=50000) DOUBLE PRECISION VERTCE(2,3,500),WRKSTR(NW),ABSERR, + VERTCQ(2,3,500) VERTCE(1,1,1) = 0 VERTCE(2,1,1) = 0 VERTCE(1,2,1) = 1 VERTCE(2,2,1) = 0 VERTCE(1,3,1) = 0 VERTCE(2,3,1) = 1 VERTCE(1,1,2) = 1 VERTCE(2,1,2) = 1 VERTCE(1,2,2) = 1 VERTCE(2,2,2) = 0 VERTCE(1,3,2) = 0 VERTCE(2,3,2) = 1 VERTCQ(1,1,1) = -1 VERTCQ(2,1,1) = 0 VERTCQ(1,2,1) = 0.5D0 VERTCQ(2,2,1) = -SQRT(3.D0)/2 VERTCQ(1,3,1) = 0.5D0 VERTCQ(2,3,1) = SQRT(3.D0)/2 ABSERR = 0.6 C C TEST FOR INTEGRATING MONOMIALS C Selected monomials, degrees 0-13 C C We first check the 105 monomials in cartesian C coordinates which DCUTRI should integrate to machine C precision. C CALL ATEST(VERTCE,1,500,1000,105,FTESTP,ABSERR,NW,0,WRKSTR, + IWORK) C C Then the same for the 21 monomials in polar coordinates. C CALL ATEST(VERTCQ,1,500,1000,21,FTESTQ,ABSERR,NW,0,WRKSTR, + IWORK) C C Test restart on a peak problem. C VERTCE(1,1,1) = 0 VERTCE(2,1,1) = 0 VERTCE(1,2,1) = 1 VERTCE(2,2,1) = 0 VERTCE(1,3,1) = 0 VERTCE(2,3,1) = 1 VERTCE(1,1,2) = 1 VERTCE(2,1,2) = 1 VERTCE(1,2,2) = 1 VERTCE(2,2,2) = 0 VERTCE(1,3,2) = 0 VERTCE(2,3,2) = 1 ABSERR = 0.1 CALL ATEST(VERTCE,2,500,74,1,FTESTO,ABSERR,NW,0,WRKSTR, + IWORK) ABSERR = 1D-3 CALL ATEST(VERTCE,2,500,1000,1,FTESTO,ABSERR,NW,1,WRKSTR, + IWORK) CALL ATEST(VERTCE,2,500,2000,1,FTESTO,ABSERR,NW,1,WRKSTR, + IWORK) END SUBROUTINE ATEST(VERTCE,NUMTRI,LENVER,MAXCLS,NFUN,TSTSUB,ABSERR, + LENWRK,IREST,WRKSTR,IWORK) EXTERNAL TSTSUB INTEGER LENWRK,IREST,NEVAL,NUMTRI,LENVER,IWORK(LENVER+1) DOUBLE PRECISION VERTCE(2,3,LENVER),ABSEST(105),FINEST(105), + WRKSTR(LENWRK),ABSERR,REL SAVE NEVAL,ABSEST,FINEST INTEGER N,MAXCLS,NFUN,IFAIL REL = 0 CALL DCUTRI(TSTSUB,NFUN,VERTCE,NUMTRI,0,MAXCLS,ABSERR,REL,LENVER, + LENWRK,IREST,FINEST,ABSEST,NEVAL,IFAIL,WRKSTR,IWORK) WRITE (*,99999) 99999 FORMAT (' '/5X,'DCUTRI TEST ') WRITE (*,99998) NEVAL,IFAIL 99998 FORMAT (5X,'SUBROUTINE CALLS = ',I6,', IFAIL = ',I2) WRITE (*,99997) 99997 FORMAT (7X,'N',3X,'ABSOLUTE ERROR',4X,'INTEGRAL') DO 10 N = 1,NFUN WRITE (*,99996) N,ABS(FINEST(N)-1),FINEST(N) 99996 FORMAT (6X,I3,E14.2,F21.16) 10 CONTINUE END SUBROUTINE FTESTP(Z,NFUN,F) C C Selected monomials, degree 0-13 C INTEGER ISAVE,NFUN,M,K,L,INDEX DOUBLE PRECISION Z(2),F(NFUN),EX(120),EXACT SAVE ISAVE,EX DATA ISAVE/0/ IF (ISAVE.EQ.0) THEN DO 100 M = 0,13 DO 90 K = 0,M L = M - K INDEX = (M* (M+1))/2 + K + 1 EX(INDEX) = EXACT(K,L) 90 CONTINUE 100 CONTINUE ISAVE = 1 END IF DO 200 M = 0,13 DO 190 K = 0,M L = M - K INDEX = (M* (M+1))/2 + K + 1 F(INDEX) = Z(1)**K*Z(2)**L/EX(INDEX) 190 CONTINUE 200 CONTINUE END SUBROUTINE FTESTQ(Z,NFUN,F) C C Selected monomials in polar coordinates, degree 0-13 C INTEGER NFUN,J,K,I,L DOUBLE PRECISION Z(2),F(NFUN),M(0:15,0:15),R,ALPHA,PI PI = 3.1415926535897932384D0 M(0,0) = 1.000000000000000000000000000000D0 M(2,0) = 0.250000000000000000000000000000D0 M(4,0) = 0.099999999999999999999999999997D0 M(6,0) = 0.051785714285714285714285714284D0 M(8,0) = 0.031428571428571428571428571427D0 M(10,0) = 0.021103896103896103896103896103D0 M(12,0) = 0.015163408020550877693734836591D0 M(14,0) = 0.011429195804195804195804195802D0 M(3,3) = -0.099999999999999999999999999997D0 M(5,3) = -0.057142857142857142857142857141D0 M(7,3) = -0.035714285714285714285714285713D0 M(9,3) = -0.024025974025974025974025974024D0 M(11,3) = -0.017132867132867132867132867132D0 M(13,3) = -0.012787212787212787212787212785D0 M(15,3) = -0.009891579009226068049597461361D0 M(6,6) = 0.035714285714285714285714285713D0 M(8,6) = 0.024999999999999999999999999997D0 M(10,6) = 0.018181818181818181818181818181D0 M(12,6) = 0.013686313686313686313686313685D0 M(14,6) = 0.010614385614385614385614385614D0 M(9,9) = -0.018181818181818181818181818181D0 M(11,9) = -0.013986013986013986013986013985D0 M(13,9) = -0.010989010989010989010989010988D0 M(15,9) = -0.008807369101486748545572074984D0 M(12,12) = 0.010989010989010989010989010988D0 M(14,12) = 0.008928571428571428571428571429D0 M(15,15) = -0.007352941176470588235294117647D0 I = 0 R = SQRT(ABS(Z(1))**2+ABS(Z(2))**2) IF (R.EQ.0) THEN ALPHA = 0 ELSE ALPHA = ACOS(Z(1)/R) END IF IF (Z(2).LT.0) THEN ALPHA = 2*PI - ALPHA END IF DO 220 J = 0,13 DO 210 K = 0,J/3 IF (MOD(J+3*K,2).EQ.0) THEN I = I + 1 IF (R.EQ.0) THEN IF (J.EQ.0) THEN F(I) = 4/ (3*SQRT(3.D0)) ELSE F(I) = 0 END IF ELSE IF (J.EQ.0) THEN F(I) = 1 ELSE F(I) = 1 DO 205 L = 1,J F(I) = F(I)*R 205 CONTINUE END IF F(I) = F(I)*COS(3*K*ALPHA) F(I) = 4*F(I)/ (3*SQRT(3.D0)*M(J,3*K)) END IF END IF 210 CONTINUE 220 CONTINUE END SUBROUTINE FTESTO(Z,NFUN,F) C C Product Peak C INTEGER NFUN DOUBLE PRECISION Z(2),F(NFUN) DOUBLE PRECISION ALPHA1,ALPHA2,BETA1,BETA2,TINT F(1) = 1.D0 ALPHA1 = 10 ALPHA2 = 10 BETA1 = 0.2D0 BETA2 = 0.3D0 TINT = 1 TINT = TINT* (ATAN((1-BETA1)*ALPHA1)-ATAN( - BETA1*ALPHA1))* + ALPHA1 TINT = TINT* (ATAN((1-BETA2)*ALPHA2)-ATAN( - BETA2*ALPHA2))* + ALPHA2 F(1) = F(1)/ (1.D0/ALPHA1**2+ (Z(1)-BETA1)**2) F(1) = F(1)/ (1.D0/ALPHA2**2+ (Z(2)-BETA2)**2) F(1) = F(1)/TINT END DOUBLE PRECISION FUNCTION EXACT(K,L) INTEGER K,L,I DOUBLE PRECISION BIN,BINOM EXACT = 0 DO 100 I = 0,L + 1 BIN = BINOM(L+1,I) EXACT = EXACT + (-1)** (I)*BIN/DBLE(I+K+1) 100 CONTINUE EXACT = EXACT/DBLE(L+1) RETURN END DOUBLE PRECISION FUNCTION BINOM(N,R) INTEGER N,R,I BINOM = 1 IF (R.EQ.0 .OR. R.EQ.N) THEN GO TO 999 END IF DO 10 I = N,N - R + 1,-1 BINOM = BINOM*I 10 CONTINUE DO 20 I = 1,R BINOM = BINOM/I 20 CONTINUE 999 RETURN END C C Output produced on a SUN SPARC station 1. C C C DCUTRI TEST C SUBROUTINE CALLS = 925, IFAIL = 0 C N ABSOLUTE ERROR INTEGRAL C 1 0.11E-15 0.9999999999999999 C 2 0.00E+00 1.0000000000000000 C 3 0.11E-15 0.9999999999999999 C 4 0.00E+00 1.0000000000000000 C 5 0.22E-15 0.9999999999999998 C 6 0.22E-15 1.0000000000000002 C 7 0.00E+00 1.0000000000000000 C 8 0.00E+00 1.0000000000000000 C 9 0.22E-15 1.0000000000000002 C 10 0.22E-15 1.0000000000000002 C 11 0.11E-14 0.9999999999999989 C 12 0.56E-15 0.9999999999999994 C 13 0.18E-14 1.0000000000000018 C 14 0.20E-14 1.0000000000000020 C 15 0.56E-15 0.9999999999999994 C 16 0.00E+00 1.0000000000000000 C 17 0.11E-14 1.0000000000000011 C 18 0.93E-14 1.0000000000000093 C 19 0.42E-14 0.9999999999999958 C 20 0.23E-14 0.9999999999999977 C 21 0.22E-15 1.0000000000000002 C 22 0.22E-15 1.0000000000000002 C 23 0.27E-14 0.9999999999999973 C 24 0.22E-13 1.0000000000000222 C 25 0.36E-14 1.0000000000000036 C 26 0.36E-14 1.0000000000000036 C 27 0.89E-15 0.9999999999999991 C 28 0.44E-15 1.0000000000000004 C 29 0.00E+00 1.0000000000000000 C 30 0.32E-13 1.0000000000000315 C 31 0.73E-13 1.0000000000000735 C 32 0.25E-13 0.9999999999999755 C 33 0.33E-14 1.0000000000000033 C 34 0.24E-13 0.9999999999999758 C 35 0.38E-14 1.0000000000000038 C 36 0.33E-15 0.9999999999999997 C 37 0.73E-14 1.0000000000000073 C 38 0.69E-13 0.9999999999999311 C 39 0.66E-13 0.9999999999999339 C 40 0.20E-13 0.9999999999999800 C 41 0.10E-12 0.9999999999998987 C 42 0.32E-13 0.9999999999999685 C 43 0.27E-13 1.0000000000000269 C 44 0.64E-14 0.9999999999999936 C 45 0.11E-14 1.0000000000000011 C 46 0.22E-15 1.0000000000000002 C 47 0.76E-13 1.0000000000000759 C 48 0.61E-12 1.0000000000006117 C 49 0.14E-12 1.0000000000001399 C 50 0.21E-13 1.0000000000000209 C 51 0.29E-12 1.0000000000002902 C 52 0.15E-12 1.0000000000001492 C 53 0.67E-14 0.9999999999999933 C 54 0.73E-14 1.0000000000000073 C 55 0.33E-15 0.9999999999999997 C 56 0.71E-13 1.0000000000000711 C 57 0.39E-12 1.0000000000003877 C 58 0.29E-12 0.9999999999997095 C 59 0.11E-11 1.0000000000010651 C 60 0.10E-11 1.0000000000010376 C 61 0.29E-12 1.0000000000002902 C 62 0.14E-12 1.0000000000001363 C 63 0.48E-13 1.0000000000000484 C 64 0.21E-13 1.0000000000000211 C 65 0.29E-14 1.0000000000000029 C 66 0.89E-15 0.9999999999999991 C 67 0.00E+00 1.0000000000000000 C 68 0.73E-12 0.9999999999992676 C 69 0.92E-12 0.9999999999990812 C 70 0.16E-11 1.0000000000016056 C 71 0.47E-11 0.9999999999952944 C 72 0.15E-11 1.0000000000015101 C 73 0.10E-11 1.0000000000010338 C 74 0.20E-12 0.9999999999997958 C 75 0.30E-12 1.0000000000002958 C 76 0.61E-13 0.9999999999999387 C 77 0.13E-13 0.9999999999999866 C 78 0.13E-14 1.0000000000000013 C 79 0.57E-13 1.0000000000000566 C 80 0.73E-12 0.9999999999992712 C 81 0.88E-11 0.9999999999911836 C 82 0.51E-11 0.9999999999949216 C 83 0.12E-10 0.9999999999877534 C 84 0.13E-12 0.9999999999998671 C 85 0.24E-11 1.0000000000023679 C 86 0.16E-11 0.9999999999983668 C 87 0.19E-11 1.0000000000018676 C 88 0.53E-12 1.0000000000005338 C 89 0.77E-13 0.9999999999999228 C 90 0.19E-13 1.0000000000000191 C 91 0.13E-14 0.9999999999999987 C 92 0.27E-12 1.0000000000002667 C 93 0.10E-12 1.0000000000000999 C 94 0.38E-11 1.0000000000037732 C 95 0.15E-10 1.0000000000152343 C 96 0.18E-10 1.0000000000178288 C 97 0.27E-10 1.0000000000267895 C 98 0.16E-10 1.0000000000157454 C 99 0.86E-11 0.9999999999913655 C 100 0.13E-10 1.0000000000132447 C 101 0.15E-11 0.9999999999984502 C 102 0.92E-12 0.9999999999990755 C 103 0.23E-12 1.0000000000002311 C 104 0.15E-13 0.9999999999999850 C 105 0.67E-15 1.0000000000000007 C C DCUTRI TEST C SUBROUTINE CALLS = 925, IFAIL = 1 C N ABSOLUTE ERROR INTEGRAL C 1 0.11E-15 0.9999999999999999 C 2 0.11E-15 0.9999999999999999 C 3 0.00E+00 1.0000000000000000 C 4 0.00E+00 1.0000000000000000 C 5 0.00E+00 1.0000000000000000 C 6 0.22E-15 0.9999999999999998 C 7 0.22E-15 1.0000000000000002 C 8 0.00E+00 1.0000000000000000 C 9 0.22E-15 0.9999999999999998 C 10 0.22E-15 0.9999999999999998 C 11 0.11E-15 0.9999999999999999 C 12 0.00E+00 1.0000000000000000 C 13 0.00E+00 1.0000000000000000 C 14 0.00E+00 1.0000000000000000 C 15 0.33E-15 0.9999999999999997 C 16 0.22E-15 0.9999999999999998 C 17 0.56E-15 0.9999999999999994 C 18 0.33E-15 0.9999999999999997 C 19 0.11E-15 0.9999999999999999 C 20 0.22E-15 0.9999999999999998 C 21 0.44E-15 0.9999999999999996 C C DCUTRI TEST C SUBROUTINE CALLS = 74, IFAIL = 1 C N ABSOLUTE ERROR INTEGRAL C 1 0.13E-01 0.9869417152262921 C C DCUTRI TEST C SUBROUTINE CALLS = 962, IFAIL = 1 C N ABSOLUTE ERROR INTEGRAL C 1 0.29E-02 0.9971072795998221 C C DCUTRI TEST C SUBROUTINE CALLS = 1702, IFAIL = 0 C N ABSOLUTE ERROR INTEGRAL C 1 0.14E-05 1.0000013854850509 SUBROUTINE DTRTRI(DVFLAG,SBRGNS,GREATE,LIST,NEW) C***BEGIN PROLOGUE DTRTRI C***REFER TO DCUTRI C***PURPOSE DTRTRI maintains a heap of subregions. C***DESCRIPTION DTRTRI maintains a heap of subregions. C The subregions are stored in a partially sorted C binary tree, ordered according to the size of the C greatest error estimates of each subregion(GREATE). C The subregion with greatest error estimate is in the C first position of the heap. C C PARAMETERS C C DVFLAG Integer. C If DVFLAG = 1, we remove the subregion with C greatest error from the heap. C If DVFLAG = 2, we insert a new subregion in the heap. C SBRGNS Integer. C Number of subregions in the heap. C GREATE Real array of dimension SBRGNS. C Used to store the greatest estimated errors in C all subregions. C LIST Integer array of dimension SBRGNS. C Used as a partially ordered list of pointers to the C different subregions. This list is a heap where the C element on top of the list is the subregion with the C greatest error estimate. C NEW Integer. C Index to the new region to be inserted in the heap. C***ROUTINES CALLED-NONE C***END PROLOGUE DTRTRI C C Global variables. C INTEGER DVFLAG,NEW,SBRGNS,LIST(*) DOUBLE PRECISION GREATE(*) C C Local variables. C C GREAT is used as intermediate storage for the greatest error of a C subregion. C SUBRGN Position of child/parent subregion in the heap. C SUBTMP Position of parent/child subregion in the heap. INTEGER SUBRGN,SUBTMP DOUBLE PRECISION GREAT C C***FIRST EXECUTABLE STATEMENT DTRTRI C C If DVFLAG = 1, we will reduce the partial ordered list by the C element with greatest estimated error. Thus the element in C in the heap with index LIST(1) is vacant and can be used later. C Reducing the heap by one element implies that the last element C should be re-positioned. C IF (DVFLAG.EQ.1) THEN GREAT = GREATE(LIST(SBRGNS)) SBRGNS = SBRGNS - 1 SUBRGN = 1 20 SUBTMP = 2*SUBRGN IF (SUBTMP.LE.SBRGNS) THEN IF (SUBTMP.NE.SBRGNS) THEN C C Find max. of left and right child. C IF (GREATE(LIST(SUBTMP)).LT. + GREATE(LIST(SUBTMP+1))) THEN SUBTMP = SUBTMP + 1 END IF END IF C C Compare max.child with parent. C If parent is max., then done. C IF (GREAT.LT.GREATE(LIST(SUBTMP))) THEN C C Move the pointer at position SUBTMP up the heap. C LIST(SUBRGN) = LIST(SUBTMP) SUBRGN = SUBTMP GO TO 20 END IF END IF C C Update the pointer. C IF (SBRGNS.GT.0) THEN LIST(SUBRGN) = LIST(SBRGNS+1) END IF ELSE IF (DVFLAG.EQ.2) THEN C C If DVFLAG = 2, find the position for the NEW region in the heap. C GREAT = GREATE(NEW) SUBRGN = SBRGNS 40 SUBTMP = SUBRGN/2 IF (SUBTMP.GE.1) THEN C C Compare max.child with parent. C If parent is max, then done. C IF (GREAT.GT.GREATE(LIST(SUBTMP))) THEN C C Move the pointer at position SUBTMP down the heap. C LIST(SUBRGN) = LIST(SUBTMP) SUBRGN = SUBTMP GO TO 40 END IF END IF C C Set the pointer to the new region in the heap. C LIST(SUBRGN) = NEW END IF C C***END DTRTRI C RETURN END SUBROUTINE DADTRI(NUMFUN,MDIV,VER,NUMTRI,MINSUB,MAXSUB,FUNSUB, + EPSABS,EPSREL,LENVER,RESTAR,LENW,RESULT,ABSERR, + NEVAL,NSUB,IFAIL,VALUES,ERRORS,GREATE,WORK2, + LIST,VACANT) C***BEGIN PROLOGUE DADTRI C***REFER TO DCUTRI C***DESCRIPTION Computation of integrals over triangular C regions. C C DADTRI repeatedly C subdivides the triangles with greatest estimated errors C and estimates the integrals and the C errors over the new subtriangles C until either the error request C is met or MAXPTS function evaluations have been used. C C A 37 point integration rule C with all evaluation points inside the triangle C is applied. The rule has polynomial degree 13. C C ON ENTRY C C NUMFUN Integer. C Number of components of the integral. C MDIV Integer. C MDIV is the number of triangles that are divided in C each subdivision step in DADTRI. C MDIV is chosen by default to be 1. C For efficient execution on parallel computers C with NPROC processors MDIV should be set equal to C the smallest integer such that MOD(4*MDIV,NPROC) = 0. C VER Real array of dimension (2,3,LENVER). C VER(1,K,L) and VER(2,K,L) are the x and y coordinates C respectively of vertex K in triangle L. C On entry VER(*,*,L) must contain the vertices of the C NUMTRI user specified triangles for L = 1,2,...,NUMTRI. C NUMTRI Integer. C The number of triangles. C MINSUB Integer. C The minimum allowed number of subregions. C MAXSUB Integer. C The maximum allowed number of subregions. C FUNSUB Externally declared subroutine for computing C all components of the integrand at the given C evaluation point. C It must have parameters (X,NUMFUN,FUNVLS) C Input parameters: C X(1) The x-coordinate of the evaluation point. C X(2) The y-coordinate of the evaluation point. C NUMFUN Integer that defines the number of C components of I. C Output parameter: C FUNVLS Real array of dimension NUMFUN C that defines NUMFUN components of the integrand. C C EPSABS Real. C Requested absolute error. C EPSREL Real. C Requested relative error. C C LENVER Integer. C Defines the length of the array VER. C C We let C MAXSUB denote the maximum allowed number of subregions C for the given values of MAXPTS. C MAXSUB = 3*((MAXPTS-37*NUMTRI)/(4*37)) + NUMTRI C LENVER should then be greater or equal to MAXSUB. C C RESTAR Integer. C If RESTAR = 0, this is the first attempt to compute C the integral. C If RESTAR = 1, C then we restart a previous attempt. C In this case the only parameters for DCUTRI that may C be changed (with respect to the previous call of DCUTRI) C are MINPTS, MAXPTS, EPSABS, EPSREL and RESTAR. C LENW Integer. C Length of the workspace WORK2. C ON RETURN C C RESULT Real array of dimension NUMFUN. C Approximations to all components of the integral. C ABSERR Real array of dimension NUMFUN. C Estimates of absolute errors. C NEVAL Integer. C Number of function evaluations used by DCUTRI. C NSUB Integer. C The number of triangles in the data structure. C IFAIL Integer. C IFAIL = 0 for normal exit. C C ABSERR(K) <= EPSABS or C ABSERR(K) <= ABS(RESULT(K))*EPSREL with MAXPTS or less C function evaluations for all values of K, C 1 <= K <= NUMFUN . C C IFAIL = 1 if MAXSUB was too small for DADTRI C to obtain the required accuracy. In this case DADTRI C returns values of RESULT with estimated absolute C errors ABSERR. C VALUES Real array of dimension (NUMFUN,MAXSUB). C The estimated components of the integrals over the C subregions. C ERRORS Real array of dimension (NUMFUN,MAXSUB). C The estimated errors over the subregions. C GREATE Real array of dimension MAXSUB. C The greatest errors in each subregion. C WORK2 Real array of dimension LENW. C Work array used in DTRTRI and DRLTRI. C LIST Integer array used in DTRTRI of dimension LENVER. C Is a partially sorted list, where LIST(1) is the top C element in a heap of subregions. C VACANT Integer array of dimension MDIV. C Used as intermediate storage of pointers. C C***ROUTINES CALLED DTRTRI,DRLTRI C***END PROLOGUE DADTRI C C Global variables. C EXTERNAL FUNSUB INTEGER NUMFUN,MDIV,NUMTRI,MINSUB,MAXSUB,LENW,RESTAR,LENVER,NEVAL, + NSUB,IFAIL,LIST(LENVER),VACANT(MDIV) DOUBLE PRECISION VER(2,3,LENVER),EPSABS,EPSREL,RESULT(NUMFUN), + ABSERR(NUMFUN),VALUES(NUMFUN,MAXSUB), + ERRORS(NUMFUN,MAXSUB),GREATE(MAXSUB),WORK2(LENW) C C Local variables. C C SBRGNS is the number of stored subregions. C NDIV The number of subregions to be divided in each main step. C POINTR Pointer to the position in the data structure where C the new subregions are to be stored. C TOP is a pointer to the top element in the heap of subregions. C INTEGER I,J,K,L INTEGER SBRGNS,I1,I2,I3,I4,SIZE INTEGER L1 INTEGER NDIV,POINTR,INDEX,TOP DOUBLE PRECISION VEROLD(2,3) C C***FIRST EXECUTABLE STATEMENT DADTRI C IF (RESTAR.EQ.1) THEN SBRGNS = NSUB GO TO 110 END IF C C Initiate SBRGNS, RESULT, ABSERR and NEVAL. C SBRGNS = 0 DO 10 J = 1,NUMFUN RESULT(J) = 0 ABSERR(J) = 0 10 CONTINUE NEVAL = 0 C C Apply DRLTRI over the NUMTRI triangles. C This loop may be run in parallel. C DO 20 I = 1,NUMTRI L1 = 1 + (I-1)*8*NUMFUN CALL DRLTRI(VER(1,1,I),NUMFUN,FUNSUB,WORK2(L1),VALUES(1,I), + ERRORS(1,I),GREATE(I)) NEVAL = NEVAL + 37 SBRGNS = SBRGNS + 1 20 CONTINUE C C Add the computed values to RESULT and ABSERR. C DO 40 I = 1,NUMTRI DO 30 J = 1,NUMFUN RESULT(J) = RESULT(J) + VALUES(J,I) ABSERR(J) = ABSERR(J) + ERRORS(J,I) 30 CONTINUE 40 CONTINUE C C Store results in heap. C DO 50 I = 1,NUMTRI INDEX = I CALL DTRTRI(2,INDEX,GREATE,LIST,I) 50 CONTINUE C C We check for termination. C IF (SBRGNS.LT.MINSUB) THEN GO TO 110 END IF DO 60 J = 1,NUMFUN IF (ABSERR(J).GT.EPSREL*ABS(RESULT(J)) .AND. + ABSERR(J).GT.EPSABS) THEN GO TO 110 END IF 60 CONTINUE IFAIL = 0 GO TO 499 C C***End initiation. C C***Begin loop while the error is too great, C and SBRGNS+3 is less than MAXSUB. C 110 IF (SBRGNS+3.LE.MAXSUB) THEN C C If we are allowed to divide further, C prepare to apply basic rule over the triangles produced C by dividing the C NDIV subtriangles with greatest errors. C If MAXSUB and SBRGNS are great enough, NDIV = MDIV. C NDIV = MAXSUB - SBRGNS NDIV = MIN(NDIV,MDIV,SBRGNS) C C Divide the NDIV subtriangles in four new subtriangles, and compute C integral and error over each. C When we pick a triangle to divide in four, then one of the new C sub-triangles is stored in the original triangle's position in the C datastructure. Thus we get 3*NDIV new elements in the datastructure C after the subdivision. The size of the datastructure before the C subdivision is stored in the variable SIZE, while SBRGNS is the C size of the heap at any time. C Hence. C SIZE = SBRGNS DO 150 I = 1,NDIV POINTR = SIZE + 3* (NDIV+1-I) C C Adjust RESULT and ABSERR. TOP is a pointer to the top of the heap. C TOP = LIST(1) VACANT(I) = TOP DO 115 J = 1,NUMFUN RESULT(J) = RESULT(J) - VALUES(J,TOP) ABSERR(J) = ABSERR(J) - ERRORS(J,TOP) 115 CONTINUE C C Save the vertices. C DO 130 L = 1,2 DO 135 K = 1,3 VEROLD(L,K) = VER(L,K,TOP) 135 CONTINUE 130 CONTINUE C C Adjust the heap. C CALL DTRTRI(1,SBRGNS,GREATE,LIST,K) C C Compute the four new triangles. C I1 = TOP I2 = POINTR - 2 I3 = POINTR - 1 I4 = POINTR VER(1,1,I1) = VEROLD(1,1) VER(2,1,I1) = VEROLD(2,1) VER(1,2,I1) = (VEROLD(1,1)+VEROLD(1,2))/2 VER(2,2,I1) = (VEROLD(2,1)+VEROLD(2,2))/2 VER(1,3,I1) = (VEROLD(1,1)+VEROLD(1,3))/2 VER(2,3,I1) = (VEROLD(2,1)+VEROLD(2,3))/2 VER(1,1,I2) = VER(1,2,I1) VER(2,1,I2) = VER(2,2,I1) VER(1,2,I2) = VEROLD(1,2) VER(2,2,I2) = VEROLD(2,2) VER(1,3,I2) = (VEROLD(1,2)+VEROLD(1,3))/2 VER(2,3,I2) = (VEROLD(2,2)+VEROLD(2,3))/2 VER(1,1,I3) = VER(1,3,I1) VER(2,1,I3) = VER(2,3,I1) VER(1,2,I3) = VER(1,3,I2) VER(2,2,I3) = VER(2,3,I2) VER(1,3,I3) = VEROLD(1,3) VER(2,3,I3) = VEROLD(2,3) VER(1,1,I4) = VER(1,3,I2) VER(2,1,I4) = VER(2,3,I2) VER(1,2,I4) = VER(1,3,I1) VER(2,2,I4) = VER(2,3,I1) VER(1,3,I4) = VER(1,2,I1) VER(2,3,I4) = VER(2,2,I1) 150 CONTINUE C C Apply basic rule over 4*NDIV triangles. C This loop may be run in parallel. C DO 200 I = 1,4*NDIV IF (I.LE.NDIV) THEN INDEX = VACANT(I) ELSE INDEX = SBRGNS + I END IF L1 = 1 + (I-1)*8*NUMFUN CALL DRLTRI(VER(1,1,INDEX),NUMFUN,FUNSUB,WORK2(L1), + VALUES(1,INDEX),ERRORS(1,INDEX),GREATE(INDEX)) 200 CONTINUE NEVAL = NEVAL + 4*NDIV*37 C C Add new contributions to RESULT and ABSERR. C DO 220 I = 1,4*NDIV IF (I.LE.NDIV) THEN INDEX = VACANT(I) ELSE INDEX = SBRGNS + I END IF DO 210 J = 1,NUMFUN RESULT(J) = RESULT(J) + VALUES(J,INDEX) ABSERR(J) = ABSERR(J) + ERRORS(J,INDEX) 210 CONTINUE 220 CONTINUE C C Store results in heap. C DO 250 I = 1,4*NDIV IF (I.LE.NDIV) THEN INDEX = VACANT(I) ELSE INDEX = SBRGNS + I END IF J = SBRGNS + I CALL DTRTRI(2,J,GREATE,LIST,INDEX) 250 CONTINUE SBRGNS = SBRGNS + 4*NDIV C C Check for termination. C IF (SBRGNS.LT.MINSUB) THEN GO TO 110 END IF DO 255 J = 1,NUMFUN IF (ABSERR(J).GT.EPSREL*ABS(RESULT(J)) .AND. + ABSERR(J).GT.EPSABS) THEN GO TO 110 END IF 255 CONTINUE IFAIL = 0 C C Else we did not succeed with the C given value of MAXSUB. C ELSE IFAIL = 1 END IF C C Compute more accurate values of RESULT and ABSERR. C 499 CONTINUE DO 500 J = 1,NUMFUN RESULT(J) = 0 ABSERR(J) = 0 500 CONTINUE DO 510 I = 1,SBRGNS DO 505 J = 1,NUMFUN RESULT(J) = RESULT(J) + VALUES(J,I) ABSERR(J) = ABSERR(J) + ERRORS(J,I) 505 CONTINUE 510 CONTINUE C NSUB = SBRGNS RETURN C C***END DADTRI C END REAL FUNCTION R1MACH(I) C C SINGLE-PRECISION MACHINE CONSTANTS C C R1MACH(1) = B**(EMIN-1), THE SMALLEST POSITIVE MAGNITUDE. C C R1MACH(2) = B**EMAX*(1 - B**(-T)), THE LARGEST MAGNITUDE. C C R1MACH(3) = B**(-T), THE SMALLEST RELATIVE SPACING. C C R1MACH(4) = B**(1-T), THE LARGEST RELATIVE SPACING. C C R1MACH(5) = LOG10(B) C C TO ALTER THIS FUNCTION FOR A PARTICULAR ENVIRONMENT, C THE DESIRED SET OF DATA STATEMENTS SHOULD BE ACTIVATED BY C REMOVING THE C FROM COLUMN 1. C ON RARE MACHINES A STATIC STATEMENT MAY NEED TO BE ADDED. C (BUT PROBABLY MORE SYSTEMS PROHIBIT IT THAN REQUIRE IT.) C C FOR IEEE-ARITHMETIC MACHINES (BINARY STANDARD), THE FIRST C SET OF CONSTANTS BELOW SHOULD BE APPROPRIATE. C C WHERE POSSIBLE, DECIMAL, OCTAL OR HEXADECIMAL CONSTANTS ARE USED C TO SPECIFY THE CONSTANTS EXACTLY. SOMETIMES THIS REQUIRES USING C EQUIVALENT INTEGER ARRAYS. IF YOUR COMPILER USES HALF-WORD C INTEGERS BY DEFAULT (SOMETIMES CALLED INTEGER*2), YOU MAY NEED TO C CHANGE INTEGER TO INTEGER*4 OR OTHERWISE INSTRUCT YOUR COMPILER C TO USE FULL-WORD INTEGERS IN THE NEXT 5 DECLARATIONS. C INTEGER SMALL(2) INTEGER LARGE(2) INTEGER RIGHT(2) INTEGER DIVER(2) INTEGER LOG10(2) INTEGER SC 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 IEEE ARITHMETIC MACHINES, SUCH AS THE AT&T C 3B SERIES, MOTOROLA 68000 BASED MACHINES (E.G. SUN 3 AND AT&T C PC 7300), AND 8087 BASED MICROS (E.G. IBM PC AND AT&T 6300). C DATA SMALL(1) / 8388608 / DATA LARGE(1) / 2139095039 / DATA RIGHT(1) / 864026624 / DATA DIVER(1) / 872415232 / DATA LOG10(1) / 1050288283 /, SC/987/ C C MACHINE CONSTANTS FOR AMDAHL MACHINES. C C DATA SMALL(1) / 1048576 / C DATA LARGE(1) / 2147483647 / C DATA RIGHT(1) / 990904320 / C DATA DIVER(1) / 1007681536 / C DATA LOG10(1) / 1091781651 /, SC/987/ 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 /, SC/987/ 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 /, SC/987/ C C MACHINE CONSTANTS FOR FTN4 ON THE CDC 6000/7000 SERIES. C C DATA RMACH(1) / 00564000000000000000B / C DATA RMACH(2) / 37767777777777777776B / C DATA RMACH(3) / 16414000000000000000B / C DATA RMACH(4) / 16424000000000000000B / C DATA RMACH(5) / 17164642023241175720B /, SC/987/ C C MACHINE CONSTANTS FOR FTN5 ON THE CDC 6000/7000 SERIES. C C DATA RMACH(1) / O"00564000000000000000" / C DATA RMACH(2) / O"37767777777777777776" / C DATA RMACH(3) / O"16414000000000000000" / C DATA RMACH(4) / O"16424000000000000000" / C DATA RMACH(5) / O"17164642023241175720" /, SC/987/ C C MACHINE CONSTANTS FOR CONVEX C-1. C C DATA RMACH(1) / '00800000'X / C DATA RMACH(2) / '7FFFFFFF'X / C DATA RMACH(3) / '34800000'X / C DATA RMACH(4) / '35000000'X / C DATA RMACH(5) / '3F9A209B'X /, SC/987/ C C MACHINE CONSTANTS FOR THE CRAY 1, XMP, 2, AND 3. C C DATA RMACH(1) / 200034000000000000000B / C DATA RMACH(2) / 577767777777777777776B / C DATA RMACH(3) / 377224000000000000000B / C DATA RMACH(4) / 377234000000000000000B / C DATA RMACH(5) / 377774642023241175720B /, SC/987/ C C MACHINE CONSTANTS FOR THE DATA GENERAL ECLIPSE S/200. C C NOTE - IT MAY BE APPROPRIATE TO INCLUDE THE FOLLOWING LINE - 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/, SC/987/ C C MACHINE CONSTANTS FOR THE HARRIS SLASH 6 AND SLASH 7. 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 /, SC/987/ C C MACHINE CONSTANTS FOR THE HONEYWELL DPS 8/70 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 /, SC/987/ C C MACHINE CONSTANTS FOR THE IBM 360/370 SERIES, C THE XEROX SIGMA 5/7/9 AND THE SEL SYSTEMS 85/86. 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 /, SC/987/ C C MACHINE CONSTANTS FOR THE INTERDATA 8/32 C WITH THE UNIX SYSTEM FORTRAN 77 COMPILER. C C FOR THE INTERDATA FORTRAN VII COMPILER REPLACE C THE Z'S SPECIFYING HEX CONSTANTS WITH Y'S. C C DATA RMACH(1) / Z'00100000' / C DATA RMACH(2) / Z'7EFFFFFF' / C DATA RMACH(3) / Z'3B100000' / C DATA RMACH(4) / Z'3C100000' / C DATA RMACH(5) / Z'41134413' /, SC/987/ 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 /, SC/987/ C C MACHINE CONSTANTS FOR PDP-11 FORTRANS 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 /, SC/987/ 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 /, SC/987/ C C MACHINE CONSTANTS FOR PDP-11 FORTRANS 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 /, SC/987/ 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 /, SC/987/ C C MACHINE CONSTANTS FOR THE SEQUENT BALANCE 8000. C C DATA SMALL(1) / $00800000 / C DATA LARGE(1) / $7F7FFFFF / C DATA RIGHT(1) / $33800000 / C DATA DIVER(1) / $34000000 / C DATA LOG10(1) / $3E9A209B /, SC/987/ 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 /, SC/987/ C C MACHINE CONSTANTS FOR THE VAX UNIX F77 COMPILER. 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 /, SC/987/ C C MACHINE CONSTANTS FOR THE VAX-11 WITH C FORTRAN IV-PLUS COMPILER. C C DATA RMACH(1) / Z00000080 / C DATA RMACH(2) / ZFFFF7FFF / C DATA RMACH(3) / Z00003480 / C DATA RMACH(4) / Z00003500 / C DATA RMACH(5) / Z209B3F9A /, SC/987/ C C MACHINE CONSTANTS FOR VAX/VMS VERSION 2.2. C C DATA RMACH(1) / '80'X / C DATA RMACH(2) / 'FFFF7FFF'X / C DATA RMACH(3) / '3480'X / C DATA RMACH(4) / '3500'X / C DATA RMACH(5) / '209B3F9A'X /, SC/987/ C C *** ISSUE STOP 778 IF ALL DATA STATEMENTS ARE COMMENTED... IF (SC .NE. 987) STOP 778 IF (I .LT. 1 .OR. I .GT. 5) GOTO 999 R1MACH = RMACH(I) RETURN 999 WRITE(*,1999) I 1999 FORMAT(' R1MACH - I OUT OF BOUNDS',I10) STOP END SUBROUTINE SADTRI(NUMFUN,MDIV,VER,NUMTRI,MINSUB,MAXSUB,FUNSUB, + EPSABS,EPSREL,LENVER,RESTAR,LENW,RESULT,ABSERR, + NEVAL,NSUB,IFAIL,VALUES,ERRORS,GREATE,WORK2, + LIST,VACANT) C***BEGIN PROLOGUE SADTRI C***REFER TO SCUTRI C***DESCRIPTION Computation of integrals over triangular C regions. C C SADTRI repeatedly C subdivides the triangles with greatest estimated errors C and estimates the integrals and the C errors over the new subtriangles C until either the error request C is met or MAXPTS function evaluations have been used. C C A 37 point integration rule C with all evaluation points inside the triangle C is applied. The rule has polynomial degree 13. C C ON ENTRY C C NUMFUN Integer. C Number of components of the integral. C MDIV Integer. C MDIV is the number of triangles that are divided in C each subdivision step in SADTRI. C MDIV is chosen by default to be 1. C For efficient execution on parallel computers C with NPROC processors MDIV should be set equal to C the smallest integer such that MOD(4*MDIV,NPROC) = 0. C VER Real array of dimension (2,3,LENVER). C VER(1,K,L) and VER(2,K,L) are the x and y coordinates C respectively of vertex K in triangle L. C On entry VER(*,*,L) must contain the vertices of the C NUMTRI user specified triangles for L = 1,2,...,NUMTRI. C NUMTRI Integer. C The number of triangles. C MINSUB Integer. C The minimum allowed number of subregions. C MAXSUB Integer. C The maximum allowed number of subregions. C FUNSUB Externally declared subroutine for computing C all components of the integrand at the given C evaluation point. C It must have parameters (X,NUMFUN,FUNVLS) C Input parameters: C X(1) The x-coordinate of the evaluation point. C X(2) The y-coordinate of the evaluation point. C NUMFUN Integer that defines the number of C components of I. C Output parameter: C FUNVLS Real array of dimension NUMFUN C that defines NUMFUN components of the integrand. C C EPSABS Real. C Requested absolute error. C EPSREL Real. C Requested relative error. C C LENVER Integer. C Defines the length of the array VER. C C We let C MAXSUB denote the maximum allowed number of subregions C for the given values of MAXPTS. C MAXSUB = 3*((MAXPTS-37*NUMTRI)/(4*37)) + NUMTRI C LENVER should then be greater or equal to MAXSUB. C C RESTAR Integer. C If RESTAR = 0, this is the first attempt to compute C the integral. C If RESTAR = 1, C then we restart a previous attempt. C In this case the only parameters for SCUTRI that may C be changed (with respect to the previous call of SCUTRI) C are MINPTS, MAXPTS, EPSABS, EPSREL and RESTAR. C LENW Integer. C Length of the workspace WORK2. C ON RETURN C C RESULT Real array of dimension NUMFUN. C Approximations to all components of the integral. C ABSERR Real array of dimension NUMFUN. C Estimates of absolute errors. C NEVAL Integer. C Number of function evaluations used by SCUTRI. C NSUB Integer. C The number of triangles in the data structure. C IFAIL Integer. C IFAIL = 0 for normal exit. C C ABSERR(K) <= EPSABS or C ABSERR(K) <= ABS(RESULT(K))*EPSREL with MAXPTS or less C function evaluations for all values of K, C 1 <= K <= NUMFUN . C C IFAIL = 1 if MAXSUB was too small for SADTRI C to obtain the required accuracy. In this case SADTRI C returns values of RESULT with estimated absolute C errors ABSERR. C VALUES Real array of dimension (NUMFUN,MAXSUB). C The estimated components of the integrals over the C subregions. C ERRORS Real array of dimension (NUMFUN,MAXSUB). C The estimated errors over the subregions. C GREATE Real array of dimension MAXSUB. C The greatest errors in each subregion. C WORK2 Real array of dimension LENW. C Work array used in STRTRI and SRLTRI. C LIST Integer array used in STRTRI of dimension LENVER. C Is a partially sorted list, where LIST(1) is the top C element in a heap of subregions. C VACANT Integer array of dimension MDIV. C Used as intermediate storage of pointers. C C***ROUTINES CALLED STRTRI,SRLTRI C***END PROLOGUE SADTRI C C Global variables. C EXTERNAL FUNSUB INTEGER NUMFUN,MDIV,NUMTRI,MINSUB,MAXSUB,LENW,RESTAR,LENVER,NEVAL, + NSUB,IFAIL,LIST(LENVER),VACANT(MDIV) REAL VER(2,3,LENVER),EPSABS,EPSREL,RESULT(NUMFUN), + ABSERR(NUMFUN),VALUES(NUMFUN,MAXSUB), + ERRORS(NUMFUN,MAXSUB),GREATE(MAXSUB),WORK2(LENW) C C Local variables. C C SBRGNS is the number of stored subregions. C NDIV The number of subregions to be divided in each main step. C POINTR Pointer to the position in the data structure where C the new subregions are to be stored. C TOP is a pointer to the top element in the heap of subregions. C INTEGER I,J,K,L INTEGER SBRGNS,I1,I2,I3,I4,SIZE INTEGER L1 INTEGER NDIV,POINTR,INDEX,TOP REAL VEROLD(2,3) C C***FIRST EXECUTABLE STATEMENT SADTRI C IF (RESTAR.EQ.1) THEN SBRGNS = NSUB GO TO 110 END IF C C Initiate SBRGNS, RESULT, ABSERR and NEVAL. C SBRGNS = 0 DO 10 J = 1,NUMFUN RESULT(J) = 0 ABSERR(J) = 0 10 CONTINUE NEVAL = 0 C C Apply SRLTRI over the NUMTRI triangles. C This loop may be run in parallel. C DO 20 I = 1,NUMTRI L1 = 1 + (I-1)*8*NUMFUN CALL SRLTRI(VER(1,1,I),NUMFUN,FUNSUB,WORK2(L1),VALUES(1,I), + ERRORS(1,I),GREATE(I)) NEVAL = NEVAL + 37 SBRGNS = SBRGNS + 1 20 CONTINUE C C Add the computed values to RESULT and ABSERR. C DO 40 I = 1,NUMTRI DO 30 J = 1,NUMFUN RESULT(J) = RESULT(J) + VALUES(J,I) ABSERR(J) = ABSERR(J) + ERRORS(J,I) 30 CONTINUE 40 CONTINUE C C Store results in heap. C DO 50 I = 1,NUMTRI INDEX = I CALL STRTRI(2,INDEX,GREATE,LIST,I) 50 CONTINUE C C We check for termination. C IF (SBRGNS.LT.MINSUB) THEN GO TO 110 END IF DO 60 J = 1,NUMFUN IF (ABSERR(J).GT.EPSREL*ABS(RESULT(J)) .AND. + ABSERR(J).GT.EPSABS) THEN GO TO 110 END IF 60 CONTINUE IFAIL = 0 GO TO 499 C C***End initiation. C C***Begin loop while the error is too great, C and SBRGNS+3 is less than MAXSUB. C 110 IF (SBRGNS+3.LE.MAXSUB) THEN C C If we are allowed to divide further, C prepare to apply basic rule over the triangles produced C by dividing the C NDIV subtriangles with greatest errors. C If MAXSUB and SBRGNS are great enough, NDIV = MDIV. C NDIV = MAXSUB - SBRGNS NDIV = MIN(NDIV,MDIV,SBRGNS) C C Divide the NDIV subtriangles in four new subtriangles, and compute C integral and error over each. C When we pick a triangle to divide in four, then one of the new C sub-triangles is stored in the original triangle's position in the C datastructure. Thus we get 3*NDIV new elements in the datastructure C after the subdivision. The size of the datastructure before the C subdivision is stored in the variable SIZE, while SBRGNS is the C size of the heap at any time. C Hence. C SIZE = SBRGNS DO 150 I = 1,NDIV POINTR = SIZE + 3* (NDIV+1-I) C C Adjust RESULT and ABSERR. TOP is a pointer to the top of the heap. C TOP = LIST(1) VACANT(I) = TOP DO 115 J = 1,NUMFUN RESULT(J) = RESULT(J) - VALUES(J,TOP) ABSERR(J) = ABSERR(J) - ERRORS(J,TOP) 115 CONTINUE C C Save the vertices. C DO 130 L = 1,2 DO 135 K = 1,3 VEROLD(L,K) = VER(L,K,TOP) 135 CONTINUE 130 CONTINUE C C Adjust the heap. C CALL STRTRI(1,SBRGNS,GREATE,LIST,K) C C Compute the four new triangles. C I1 = TOP I2 = POINTR - 2 I3 = POINTR - 1 I4 = POINTR VER(1,1,I1) = VEROLD(1,1) VER(2,1,I1) = VEROLD(2,1) VER(1,2,I1) = (VEROLD(1,1)+VEROLD(1,2))/2 VER(2,2,I1) = (VEROLD(2,1)+VEROLD(2,2))/2 VER(1,3,I1) = (VEROLD(1,1)+VEROLD(1,3))/2 VER(2,3,I1) = (VEROLD(2,1)+VEROLD(2,3))/2 VER(1,1,I2) = VER(1,2,I1) VER(2,1,I2) = VER(2,2,I1) VER(1,2,I2) = VEROLD(1,2) VER(2,2,I2) = VEROLD(2,2) VER(1,3,I2) = (VEROLD(1,2)+VEROLD(1,3))/2 VER(2,3,I2) = (VEROLD(2,2)+VEROLD(2,3))/2 VER(1,1,I3) = VER(1,3,I1) VER(2,1,I3) = VER(2,3,I1) VER(1,2,I3) = VER(1,3,I2) VER(2,2,I3) = VER(2,3,I2) VER(1,3,I3) = VEROLD(1,3) VER(2,3,I3) = VEROLD(2,3) VER(1,1,I4) = VER(1,3,I2) VER(2,1,I4) = VER(2,3,I2) VER(1,2,I4) = VER(1,3,I1) VER(2,2,I4) = VER(2,3,I1) VER(1,3,I4) = VER(1,2,I1) VER(2,3,I4) = VER(2,2,I1) 150 CONTINUE C C Apply basic rule over 4*NDIV triangles. C This loop may be run in parallel. C DO 200 I = 1,4*NDIV IF (I.LE.NDIV) THEN INDEX = VACANT(I) ELSE INDEX = SBRGNS + I END IF L1 = 1 + (I-1)*8*NUMFUN CALL SRLTRI(VER(1,1,INDEX),NUMFUN,FUNSUB,WORK2(L1), + VALUES(1,INDEX),ERRORS(1,INDEX),GREATE(INDEX)) 200 CONTINUE NEVAL = NEVAL + 4*NDIV*37 C C Add new contributions to RESULT and ABSERR. C DO 220 I = 1,4*NDIV IF (I.LE.NDIV) THEN INDEX = VACANT(I) ELSE INDEX = SBRGNS + I END IF DO 210 J = 1,NUMFUN RESULT(J) = RESULT(J) + VALUES(J,INDEX) ABSERR(J) = ABSERR(J) + ERRORS(J,INDEX) 210 CONTINUE 220 CONTINUE C C Store results in heap. C DO 250 I = 1,4*NDIV IF (I.LE.NDIV) THEN INDEX = VACANT(I) ELSE INDEX = SBRGNS + I END IF J = SBRGNS + I CALL STRTRI(2,J,GREATE,LIST,INDEX) 250 CONTINUE SBRGNS = SBRGNS + 4*NDIV C C Check for termination. C IF (SBRGNS.LT.MINSUB) THEN GO TO 110 END IF DO 255 J = 1,NUMFUN IF (ABSERR(J).GT.EPSREL*ABS(RESULT(J)) .AND. + ABSERR(J).GT.EPSABS) THEN GO TO 110 END IF 255 CONTINUE IFAIL = 0 C C Else we did not succeed with the C given value of MAXSUB. C ELSE IFAIL = 1 END IF C C Compute more accurate values of RESULT and ABSERR. C 499 CONTINUE DO 500 J = 1,NUMFUN RESULT(J) = 0 ABSERR(J) = 0 500 CONTINUE DO 510 I = 1,SBRGNS DO 505 J = 1,NUMFUN RESULT(J) = RESULT(J) + VALUES(J,I) ABSERR(J) = ABSERR(J) + ERRORS(J,I) 505 CONTINUE 510 CONTINUE C NSUB = SBRGNS RETURN C C***END SADTRI C END SUBROUTINE SCHTRI(NUMFUN,MDIV,VER,NUMTRI,MINPTS,MAXPTS,EPSABS, + EPSREL,LENVER,NW,RESTAR,MAXSUB,MINSUB,IFAIL) C***BEGIN PROLOGUE SCHTRI C***REFER TO SCUTRI C***PURPOSE SCHTRI checks the validity of the C input parameters to SCUTRI. C***DESCRIPTION C SCHTRI computes MAXSUB, MINSUB and IFAIL as C functions of the input parameters to SCUTRI, C and checks the validity of the input parameters to SCUTRI. C C ON ENTRY C C NUMFUN Integer. C Number of components of the integral. C MDIV Integer. C MDIV is the number of triangles that are divided in C each subdivision step in SADTRI. C MDIV is chosen default to 1. C For efficient execution on parallel computers C with NPROC processors MDIV should be set equal to C the smallest integer such that MOD(4*MDIV,NPROC) = 0. C VER Real array of dimension (2,3,LENVER). C VER(1,K,L) and VER(2,K,L) are the x and y coordinates C respectively of vertex K in triangle L. C On entry VER(*,*,L) must contain the vertices of the C NUMTRI user specified triangles for L = 1,2,...,NUMTRI. C NUMTRI Integer. C The number of triangles. C MINPTS Integer. C Minimum number of function evaluations. C MAXPTS Integer. C Maximum number of function evaluations. C The number of function evaluations over each subregion C is 37. C MAXPTS >= 37*NUMTRI and MAXPTS >= MINPTS C C EPSABS Real. C Requested absolute error. C EPSREL Real. C Requested relative error. C LENVER Integer. C Defines the length of the array VER. C C We let C MAXSUB denote the maximum allowed number of subregions C for the given values of MAXPTS. C MAXSUB = 3*((MAXPTS-37*NUMTRI)/(4*37)) + NUMTRI C LENVER should then be greater or equal to MAXSUB. C C NW Integer. C Defines the length of the working array WORK. C C We let C MAXSUB denote the maximum allowed number of subregions C for the given values of MAXPTS. C MAXSUB = 3*((MAXPTS-37*NUMTRI)/(4*37)) + NUMTRI C NW should then be greater or equal to C MAXSUB*(2*NUMFUN+1) + MAX(32*MDIV,8*NUMTRI)*NUMFUN + 1 C C RESTAR Integer. C If RESTAR = 0, this is the first attempt to compute C the integral. C If RESTAR = 1, C then we restart a previous attempt. C In this case the only parameters for SCUTRI that may C be changed (with respect to the previous call of SCUTRI) C are MINPTS, MAXPTS, EPSABS, EPSREL and RESTAR. C C C ON RETURN C C MAXSUB Integer. C The maximum allowed number of subregions for the C given values of MAXPTS. C MINSUB Integer. C The minimum allowed number of subregions for the given C values of MINPTS. C IFAIL Integer. C IFAIL = 0 for normal exit. C IFAIL = 2 if NUMFUN is less than 1. C IFAIL = 3 if area of one of the initially given C triangles is zero. C IFAIL = 4 if MAXPTS is less than 37*NUMTRI. C IFAIL = 5 if MAXPTS is less than MINPTS. C IFAIL = 6 if EPSABS <= 0 and EPSREL <= 0. C IFAIL = 7 if LENVER is less than MAXSUB. C IFAIL = 8 if NW is less than MAXSUB*(2*NUMFUN+1) + C NUMFUN*MAX(32*MDIV,8*NUMTRI) + 1. C IFAIL = 9 if illegal RESTAR. C C***ROUTINES CALLED-NONE C***END PROLOGUE SCHTRI C C Global variables. C INTEGER NUMFUN,MDIV,NUMTRI,MINPTS,MAXPTS,NW,MAXSUB,MINSUB,RESTAR, + LENVER,IFAIL REAL VER(2,3,NUMTRI),EPSABS,EPSREL C C Local variables. C INTEGER LIMIT,J REAL AREA C C***FIRST EXECUTABLE STATEMENT SCHTRI C IFAIL = 0 C C We compute MAXSUB and MINSUB. C MAXSUB = 3* ((MAXPTS-37*NUMTRI)/ (4*37)) + NUMTRI MINSUB = 3* ((MINPTS-37*NUMTRI)/ (4*37)) + NUMTRI IF (MOD((MINPTS-37*NUMTRI),4*37).GT.0) THEN MINSUB = MINSUB + 3 END IF MINSUB = MAX(NUMTRI,MINSUB) C C Check on positive NUMFUN. C IF (NUMFUN.LT.1) THEN IFAIL = 2 RETURN END IF C C Check on legal area of the region of integration. C DO 10 J = 1,NUMTRI AREA = ABS(VER(1,1,J)*VER(2,2,J)-VER(1,2,J)*VER(2,1,J)- + VER(1,1,J)*VER(2,3,J)+VER(1,3,J)*VER(2,1,J)+ + VER(1,2,J)*VER(2,3,J)-VER(1,3,J)*VER(2,2,J))*0.5 IF (AREA.EQ.0) THEN IFAIL = 3 RETURN END IF 10 CONTINUE C C Check on MAXPTS >= 37*NUMTRI. C IF (MAXPTS.LT.37*NUMTRI) THEN IFAIL = 4 RETURN END IF C C Check on MAXPTS >= MINPTS. C IF (MAXPTS.LT.MINPTS) THEN IFAIL = 5 RETURN END IF C C Check on legal accuracy requests. C IF (EPSABS.LE.0 .AND. EPSREL.LE.0) THEN IFAIL = 6 RETURN END IF C C Check on big enough LENVER. C IF (LENVER.LT.MAXSUB) THEN IFAIL = 7 RETURN END IF C C Check on big enough NW. C LIMIT = MAXSUB* (2*NUMFUN+1) + MAX(32*MDIV,8*NUMTRI)*NUMFUN + 1 IF (NW.LT.LIMIT) THEN IFAIL = 8 RETURN END IF C C Check on legal RESTAR. C IF (RESTAR.NE.0 .AND. RESTAR.NE.1) THEN IFAIL = 9 RETURN END IF RETURN C C***END SCHTRI C END SUBROUTINE SCUTRI(FUNSUB,NUMFUN,VER,NUMTRI,MINPTS,MAXPTS,EPSABS, + EPSREL,LENVER,NW,RESTAR,RESULT,ABSERR,NEVAL, + IFAIL,WORK,IWORK) C***BEGIN PROLOGUE SCUTRI C***DATE WRITTEN 891023 (YYMMDD) C***REVISION DATE 910214 (YYMMDD) C***CATEGORY NO. H2B2A1/H2B2A2 C***KEYWORDS QUADRATURE, TWO DIMENSIONAL, ADAPTIVE, CUBATURE C***AUTHOR C C Jarle Berntsen, Institute of Marine Research, C Nordnesparken 2, Postboks 1870, C N-5024 Bergen-Nordnes, Norway C Phone.. 47-5-238461 C Email.. jarleb@imr.no C C Terje O. Espelid, Department of Informatics, C University of Bergen, Thormohlens gt. 55, C N-5008 Bergen, Norway C Phone.. 47-5-544180 C Email.. terje@eik.ii.uib.no C C***PURPOSE To compute the two-dimensional integral over a region C consisting of NUMTRI triangles. C***DESCRIPTION C C The routine calculates an approximation to a given C vector of definite integrals C C I I (F ,F ,...,F ) DX(2)DX(1), C 1 2 NUMFUN C C where the region of integration is a selection of C NUMTRI triangles and C where F = F (X(1),X(2)), J = 1,2,...,NUMFUN. C J J C C C A globally adaptive strategy is applied in order to C compute approximations RESULT(K) C hopefully satisfying for each component of I the following C claim for accuracy: C ABS(I(K)-RESULT(K)).LE.MAX(EPSABS,EPSREL*ABS(I(K))) C C SCUTRI is a driver for the integration routine C SADTRI. C C SADTRI repeatedly C subdivides the triangles with greatest estimated errors C and estimates the integrals and the C errors over the new subtriangles C until either the error request C is met or MAXPTS function evaluations have been used. C C A 37 point integration rule C with all evaluation points inside the triangle C is applied. The rule has polynomial degree 13. C C If the values of the input parameters EPSABS C or EPSREL are selected great enough, C an integration rule is applied over each triangle and C the results are added up to give the approximations C RESULT(K). No further subdivision C of the triangles will then be applied. C C When SCUTRI computes estimates to a vector of C integrals, all components of the vector are given C the same treatment. That is, I(F ) and I(F ) for C J K C J not equal to K, are estimated with the same C subdivision of the region of integration. C For integrals with enough similarity, we may save C time by applying SCUTRI to all integrands in one call. C For integrals that varies continuously as functions of C some parameter, the estimates produced by SCUTRI will C also vary continuously when the same subdivision is C applied to all components. This will generally not be C the case when the different components are given C separate treatment. C C On the other hand this feature should be used with C caution when the different components of the integrals C require clearly different subdivisions. C C ON ENTRY C C FUNSUB Externally declared subroutine for computing C all components of the integrand at the given C evaluation point. C It must have parameters (X,NUMFUN,FUNVLS) C Input parameters: C X(1) The x-coordinate of the evaluation point. C X(2) The y-coordinate of the evaluation point. C NUMFUN Integer that defines the number of C components of I. C Output parameter: C FUNVLS Real array of dimension NUMFUN C that defines NUMFUN components of the integrand. C C NUMFUN Integer. C Number of components of the integral. C VER Real array of dimension (2,3,LENVER). C VER(1,K,L) and VER(2,K,L) are the x and y coordinates C respectively of vertex K in triangle L. C On entry VER(*,*,L) must contain the vertices of the C NUMTRI user specified triangles for L = 1,2,...,NUMTRI. C VER may be changed on exit. C C NUMTRI Integer. C The number of triangles. C MINPTS Integer. C Minimum number of function evaluations. C MAXPTS Integer. C Maximum number of function evaluations. C The number of function evaluations over each subregion C is 37. C C MAXPTS >= 37*NUMTRI and MAXPTS >= MINPTS C C EPSABS Real. C Requested absolute error. C EPSREL Real. C Requested relative error. C LENVER Integer. C Defines the length of the array VER. C C We let C MAXSUB denote the maximum allowed number of subregions C for the given value of MAXPTS. C MAXSUB = 3*((MAXPTS-37*NUMTRI)/(4*37)) + NUMTRI C LENVER should be greater or equal to MAXSUB. C C NW Integer. C Defines the length of the working array WORK. C C If LENVER is chosen correctly as a function of MAXPTS C and NUMTRI and NW is selected to be greater or equal to C LENVER*(2*NUMFUN+1) + MAX(32*MDIV,8*NUMTRI)*NUMFUN + 1 , C then the size of the working storage will be great enough. C MDIV is the number of triangles that are divided in C each subdivision step. C MDIV is default set internally in SCUTRI equal to 1. C For efficient execution on parallel computers C with NPROC processors MDIV should be set equal to C the smallest integer such that MOD(4*MDIV,NPROC) = 0. C C RESTAR Integer. C If RESTAR = 0, this is the first attempt to compute C the integral. C If RESTAR = 1, C then we restart a previous attempt. C In this case the only parameters for SCUTRI that may C be changed (with respect to the previous call of SCUTRI) C are MINPTS, MAXPTS, EPSABS, EPSREL and RESTAR. C Note: If MAXPTS was too small in the previous call, C then we get a second chance to continue the computations C with MAXPTS increased. C LENVER can not be changed, therefore the new value of MAXPTS C should not be chosen greater than the value of LENVER C allows us to do. C WORK Real array of dimension NW. C Used as working storage. C IWORK Integer array of dimension LENVER + MDIV. C Used as working storage. C C ON RETURN C C RESULT Real array of dimension NUMFUN. C Approximations to all components of the integral. C ABSERR Real array of dimension NUMFUN. C Estimates of absolute errors. C NEVAL Integer. C Number of function evaluations used by SCUTRI. C IFAIL Integer. C IFAIL = 0 for normal exit. C C ABSERR(K) <= EPSABS or C ABSERR(K) <= ABS(RESULT(K))*EPSREL with MAXPTS or less C function evaluations for all values of K, C 1 <= K <= NUMFUN . C C IFAIL = 1 if MAXPTS was too small for SCUTRI C to obtain the required accuracy. In this case SCUTRI C returns values of RESULT with estimated absolute C errors ABSERR. C IFAIL = 2 if NUMFUN is less than 1. C IFAIL = 3 if area of one of the initially given C triangles is zero. C IFAIL = 4 if MAXPTS is less than 37*NUMTRI. C IFAIL = 5 if MAXPTS is less than MINPTS. C IFAIL = 6 if EPSABS <= 0 and EPSREL <= 0. C IFAIL = 7 if LENVER is less than MAXSUB. C IFAIL = 8 if NW is less than MAXSUB*(2*NUMFUN+1) + C NUMFUN*MAX(32*MDIV,8*NUMTRI) + 1. C IFAIL = 9 if unlegal RESTAR. C VER Real array of dimension (2,3,LENVER). C On exit VER C contains the vertices of the triangles used to produce C the approximations to the integrals. C WORK Real array of dimension NW. C Used as working storage. C WORK(NW) = NSUB, the number of subregions in the data C structure. C Let WRKSUB=(NW-1-NUMFUN*MAX(32*MDIV,8*NUMTRI))/ C (2*NUMFUN+1). C WORK(1),...,WORK(NUMFUN*WRKSUB) contain C the estimated components of the integrals over the C subregions. C WORK(NUMFUN*WRKSUB+1),...,WORK(2*NUMFUN*WRKSUB) contain C the estimated errors over the subregions. C WORK(2*NUMFUN*WRKSUB+1),...,WORK(2*NUMFUN* C WRKSUB+WRKSUB) contain the greatest errors C in each subregion. C WORK((2*NUMFUN+1)*WRKSUB),...,WORK(NW - 1) is used as C temporary storage in SADTRI. C IWORK Integer array of dimension LENVER + MDIV. C Used as working storage. C C SAMPLE PROGRAM C The following program will approximate the integral of C exp(x*x+y*y) C over the triangle (0.,0.),(2.,0.),(0.,2.) and print out the C values of the estimated integral, the estimated error, the number C of function evaluations, and IFAIL. C PROGRAM STEST1 C INTEGER NUMFUN,NW,MDIV,MINPTS,LENVER,NUMTRI C PARAMETER (NUMFUN=1,MDIV=1,NUMTRI=1,MINPTS=37) C PARAMETER ( LENVER = 100 ) C C We need to choose LENVER such that: C LENVER >= 3*((MAXPTS-37*NUMTRI)/(4*37)) + NUMTRI C This simply means that we have enough space to achieve MAXPTS C function evaluations. By choosing LENVER bigger than C this lower bound we may later change MAXPTS and RESTAR and restart C the computations from the point where the last run stopped. C The reason for stopping may have been that MAXPTS was too small C to achieve the requested error. C With our choice LENVER = 100 any value of MAXPTS C between 37 and 5068 will be accepted by the code. Choosing C MAXPTS = 2000 allows us to restart with a greater value C if necessary. C C PARAMETER (NW = LENVER*(2*NUMFUN+1) + C + MAX(32*MDIV,8*NUMTRI)*NUMFUN + 1) C REAL VER(2,3,LENVER),RESULT(NUMFUN), C + ABSERR(NUMFUN),WORK(NW),EPSABS,EPSREL C INTEGER RESTAR,NEVAL,IFAIL,IWORK(LENVER+MDIV),MAXPTS C EXTERNAL F C VER(1,1,1) = 0.E0 C VER(1,2,1) = 2.E0 C VER(1,3,1) = 0.E0 C VER(2,1,1) = 0.E0 C VER(2,2,1) = 0.E0 C VER(2,3,1) = 2.E0 C EPSABS = 0.E0 C EPSREL = 1.E-5 C RESTAR = 0 C MAXPTS = 2000 C CALL SCUTRI(F,NUMFUN,VER,NUMTRI,MINPTS,MAXPTS,EPSABS, C + EPSREL,LENVER,NW,RESTAR,RESULT,ABSERR,NEVAL, C + IFAIL,WORK,IWORK) C WRITE(*,*)'RESULT=',RESULT(1) C WRITE(*,*)'ERROR =',ABSERR(1) C WRITE(*,*)'NEVAL =',NEVAL C WRITE(*,*)'IFAIL =',IFAIL C END C SUBROUTINE F(X,NUMFUN,FUNVLS) C INTEGER NUMFUN C REAL X(2),FUNVLS(NUMFUN) C FUNVLS(1) = EXP(X(1)*X(1)+X(2)*X(2)) C RETURN C END C C Output produced on a SUN SPARC station 1. C C RESULT= 11.1813 C ERROR = 4.71008E-06 C NEVAL = 185 C IFAIL = 0 C C C***LONG DESCRIPTION C C The information for each triangle is contained in the C data structures VER, WORK and IWORK. C VER contains the coordinates of the triangles. C When passed on to SADTRI, WORK is split into four C arrays VALUES, ERRORS, GREATE and WORK2. C VALUES contains the estimated values of the integrals. C ERRORS contains the estimated errors. C GREATE contains the greatest estimated error for each triangle. C The data structures for the triangles are in SADTRI organized C as a heap, and the size of GREATE(I) defines the position of C triangle I in the heap. The heap is maintained by the program C DTRTRI and we use a partially ordered list of pointers, saved C in IWORK. When passed on to SADTRI, IWORK is split into two C arrays LIST and VACANT. LIST is a partially ordered list C such that GREATE(LIST(1)) is the maximum error estimate for C all sub-triangles in our datastructure. VACANT is a work space C needed in the updating process of the list. C C The subroutine SADTRI is written for efficient execution on shared C memory parallel computer. On a computer with NPROC processors we will C in each subdivision step divide MDIV triangles, where MDIV is C chosen such that MOD(4*MDIV,NPROC) = 0, in totally 4*MDIV new triangles. C Each processor will then compute estimates of the integrals and errors C over 4*MDIV/NPROC triangles in each subdivision step. C The subroutine for estimating the integral and the error over C each subregion, DRLTRI, uses WORK2 as a work array. C We must make sure that each processor writes its results to C separate parts of the memory, and therefore the sizes of WORK and C WORK2 are functions of MDIV. C In order to achieve parallel processing of triangles, compiler C directives should be placed in front of the DO 20 and the DO 200 C loops in SADTRI on machines like Alliant and CRAY. C C***REFERENCES C J.Berntsen and T.O.Espelid, An Algorithm for Adaptive C Cubature Over a Collection of Triangles, to be published. C C***ROUTINES CALLED SCHTRI,SADTRI C***END PROLOGUE SCUTRI C C Parameters C C MDIV Integer. C MDIV is the number of triangles that are divided in C each subdivision step in SADTRI. C MDIV is chosen default to 1. C For efficient execution on parallel computers C with NPROC processors MDIV should be set equal to C the smallest integer such that MOD(4*MDIV,NPROC) = 0. INTEGER MDIV PARAMETER (MDIV=1) C C Global variables. C EXTERNAL FUNSUB INTEGER NUMFUN,NUMTRI,MINPTS,MAXPTS,LENVER,NW,RESTAR,NEVAL,IFAIL, + IWORK(LENVER+MDIV) REAL VER(2,3,LENVER),EPSABS,EPSREL,RESULT(NUMFUN), + ABSERR(NUMFUN),WORK(NW) C C Local variables. C C MAXSUB Integer. C The maximum allowed number of subregions C for the given values of MAXPTS. C MINSUB Integer. C The minimum allowed number of subregions for the given C values of MINPTS. C WRKSUB Integer. C Determines the length of the main work arrays. C INTEGER MAXSUB,MINSUB,NSUB,LENW INTEGER WRKSUB,I1,I2,I3,I4 C C***FIRST EXECUTABLE STATEMENT SCUTRI C C Compute MAXSUB and MINSUB, C and check the input parameters. C CALL SCHTRI(NUMFUN,MDIV,VER,NUMTRI,MINPTS,MAXPTS,EPSABS,EPSREL, + LENVER,NW,RESTAR,MAXSUB,MINSUB,IFAIL) WRKSUB = (NW-1-NUMFUN*MAX(32*MDIV,8*NUMTRI))/ (2*NUMFUN+1) IF (IFAIL.NE.0) THEN RETURN END IF C C Split up the work space. C I1 = 1 I2 = I1 + WRKSUB*NUMFUN I3 = I2 + WRKSUB*NUMFUN I4 = I3 + WRKSUB C C On restart runs the number of subregions from the C previous call is assigned to NSUB. C IF (RESTAR.EQ.1) THEN NSUB = WORK(NW) END IF C C Compute the size of the temporary work space needed in SADTRI. C LENW = MAX(32*MDIV,8*NUMTRI)*NUMFUN CALL SADTRI(NUMFUN,MDIV,VER,NUMTRI,MINSUB,MAXSUB,FUNSUB,EPSABS, + EPSREL,LENVER,RESTAR,LENW,RESULT,ABSERR,NEVAL,NSUB, + IFAIL,WORK(I1),WORK(I2),WORK(I3),WORK(I4), + IWORK(1),IWORK(1+LENVER)) WORK(NW) = NSUB RETURN C C***END SCUTRI C END SUBROUTINE SRLTRI(VER,NUMFUN,FUNSUB,NULL,BASVAL,RGNERR,GREATE) C***BEGIN PROLOGUE SRLTRI C***REFER TO SCUTRI C***PURPOSE To compute basic integration rule values and C corresponding error estimates. C***DESCRIPTION SRLTRI computes basic integration rule values C for a vector of integrands over a triangular region. C SRLTRI also computes estimates for the errors by C using several null rule approximations. C ON ENTRY C C VER Real array of dimension (2,3). C The coordinates of the vertices of the triangle. C NUMFUN Integer. C Number of components of the vector integrand. C FUNSUB Externally declared subroutine for computing C all components of the integrand at the given C evaluation point. C It must have parameters (X,NUMFUN,FUNVLS) C Input parameters: C X(1) The x-coordinate of the evaluation point. C X(2) The y-coordinate of the evaluation point. C NUMFUN Integer that defines the number of C components of I. C Output parameter: C FUNVLS Real array of dimension NUMFUN C that defines NUMFUN components of the integrand. C C NULL Real array of dimension (NUMFUN,8). C A work array. C C ON RETURN C C BASVAL Real array of dimension NUMFUN. C The values for the basic rule for each component C of the integrand. C RGNERR Real array of dimension NUMFUN. C The error estimates for each component of the integrand. C GREATE Real. C The greatest error component of the integrand. C C***REFERENCES Berntsen,J. and Espelid,T.O., Degree 13 Symmetric C Quadrature Rules for the Triangle, Reports in C Informatics 44, Dept. of Informatics, University of C Bergen, 1990. C***ROUTINES CALLED R1MACH,FUNSUB C***END PROLOGUE SRLTRI C C Global variables. C EXTERNAL FUNSUB INTEGER NUMFUN REAL VER(2,3),BASVAL(NUMFUN),RGNERR(NUMFUN), + NULL(NUMFUN,8),GREATE,NOISE,R1MACH,TRES C C Local variables. C C WTLENG The number of weights of the integration rule. C G Real array of dimension (2,WTLENG). C The homogeneous coordinates for the generators of C the evaluation points. C The integration rule is using symmetric evaluation C points and has the structure (1,6,3). That is, C 1 point of multiplicity 1, C 6 sets of points of multiplicity 3 and C 3 sets of points of multiplicity 6. C This gives totally 37 evaluation points. C In order to reduce the number of loops in SRLTRI, C the 3 loops for the sets of multiplicity 6 are split C into 6 loops and added to the loops for the sets of C multiplicity 3. C The number of weights we have to give with C this splitting is WTLENG = 13. C C W Real array of dimension (9,WTLENG). C The weights of the basic rule and the null rules. C W(1,1),...,W(1,WTLENG) are weights for the basic rule. C W(I,1),...,W(I,WTLENG) for I>1 are null rule weights. C INTEGER WTLENG PARAMETER (WTLENG=13) REAL AREA,X(2,3),Z1,Z2,Z3,R1,R2,R3,R,DEG7,DEG5,DEG3, + DEG1,G(2,13),W(9,13) INTEGER I,J,K,L C C The abscissas are given in homogeneous coordinates. C DATA (G(1,I),I=1,13)/0.333333333333333333333333333333E0, + 0.950275662924105565450352089520E0, + 0.171614914923835347556304795551E0, + 0.539412243677190440263092985511E0, + 0.772160036676532561750285570113E0, + 0.009085399949835353883572964740E0, + 0.062277290305886993497083640527E0, + 0.022076289653624405142446876931E0, + 0.018620522802520968955913511549E0, + 0.096506481292159228736516560903E0, + 0.851306504174348550389457672223E0, + 0.689441970728591295496647976487E0, + 0.635867859433872768286976979827E0/ DATA (G(2,I),I=1,13)/0.333333333333333333333333333333E0, + 0.024862168537947217274823955239E0, + 0.414192542538082326221847602214E0, + 0.230293878161404779868453507244E0, + 0.113919981661733719124857214943E0, + 0.495457300025082323058213517632E0, + 0.468861354847056503251458179727E0, + 0.851306504174348550389457672223E0, + 0.689441970728591295496647976487E0, + 0.635867859433872768286976979827E0, + 0.022076289653624405142446876931E0, + 0.018620522802520968955913511549E0, + 0.096506481292159228736516560903E0/ C C Weights of the degree 13 quadrature rule. C DATA (W(1,I),I=1,13)/0.051739766065744133555179145422E0, + 0.008007799555564801597804123460E0, + 0.046868898981821644823226732071E0, + 0.046590940183976487960361770070E0, + 0.031016943313796381407646220131E0, + 0.010791612736631273623178240136E0, + 0.032195534242431618819414482205E0, + 0.015445834210701583817692900053E0, + 0.017822989923178661888748319485E0, + 0.037038683681384627918546472190E0, + 0.015445834210701583817692900053E0, + 0.017822989923178661888748319485E0, + 0.037038683681384627918546472190E0/ C C Weights of the first null rule of degree 7. C DATA (W(2,I),I=1,13)/-0.077738051051462052051304462750E0, + 0.001640389740236881582083124927E0, + 0.078124083459915167386776552733E0, + -0.030706528522391137165581298102E0, + 0.010246307817678312345028512621E0, + 0.012586300774453821540476193059E0, + -0.043630506151410607808929481439E0, + -0.004567055157220063810223671248E0, + 0.003393373439889186878847613140E0, + 0.000000000000000000000000000000E0, + -0.004567055157220063810223671248E0, + 0.003393373439889186878847613140E0, + 0.000000000000000000000000000000E0/ C C Weights of the second null rule of degree 7. C DATA (W(3,I),I=1,13)/-0.064293709240668260928898888457E0, + 0.003134665264639380635175608661E0, + 0.007822550509742830478456728602E0, + 0.048653051907689492781049400973E0, + 0.032883327334384971735434067029E0, + -0.017019508374229390108580829589E0, + 0.025973557893399824586684707198E0, + -0.010716753326806275930657622320E0, + 0.018315629578968063765722278290E0, + -0.047607080313197299401024682666E0, + -0.010716753326806275930657622320E0, + 0.018315629578968063765722278290E0, + -0.047607080313197299401024682666E0/ C C Weights of the first degree 5 null rule. C C DATA (W(4,I),I=1,13)/0.021363205584741860993131879186E0, + 0.022716410154120323440432428315E0, + -0.026366191271182090678117381002E0, + 0.029627021479068212693155637482E0, + 0.004782834546596399307634111034E0, + 0.004178667433984132052378990240E0, + -0.065398996748953861618846710897E0, + -0.033589813176131630980793760168E0, + 0.033018320112481615757912576257E0, + 0.012241086002709814125707333127E0, + -0.033589813176131630980793760168E0, + 0.033018320112481615757912576257E0, + 0.012241086002709814125707333127E0/ C C Weights of the second degree 5 null rule. C DATA (W(5,I),I=1,13)/-0.046058756832790538620830792345E0, + 0.005284159186732627192774759959E0, + 0.009325799301158899112648198129E0, + -0.006101110360950124560783393745E0, + -0.056223328794664871336486737231E0, + -0.062516479198185693171971930698E0, + 0.022428226812039547178810743269E0, + -0.000026014926110604563130107142E0, + 0.032882099937471182365626663487E0, + 0.018721740987705986426812755881E0, + -0.000026014926110604563130107142E0, + 0.032882099937471182365626663487E0, + 0.018721740987705986426812755881E0/ C C Weights of first degree 3 null rule. C DATA (W(6,I),I=1,13)/0.080867117677405246540283712799E0, + -0.033915806661511608094988607349E0, + 0.014813362053697845461526433401E0, + 0.001442315416337389214102507204E0, + -0.024309696484708683486455879210E0, + -0.005135085639122398522835391664E0, + -0.034649417896235909885490654650E0, + 0.035748423431577326597742956780E0, + 0.024548155266816447583155562333E0, + -0.032897267038856299280541675107E0, + 0.035748423431577326597742956780E0, + 0.024548155266816447583155562333E0, + -0.032897267038856299280541675107E0/ C C Weights of second degree 3 null rule. C DATA (W(7,I),I=1,13)/-0.038457863913548248582247346193E0, + -0.055143631258696406147982448269E0, + -0.021536994314510083845999131455E0, + 0.001547467894857008228010564582E0, + 0.057409361764652373776043522086E0, + -0.040636938884669694118908764512E0, + -0.020801144746964801777584428369E0, + 0.019490770404993674256256421103E0, + 0.002606109985826399625043764771E0, + 0.023893703367437102825618048130E0, + 0.019490770404993674256256421103E0, + 0.002606109985826399625043764771E0, + 0.023893703367437102825618048130E0/ C C Weights of first degree 1 null rule. C DATA (W(8,I),I=1,13)/0.074839568911184074117081012527E0, + -0.004270103034833742737299816615E0, + 0.049352639555084484177095781183E0, + 0.048832124609719176627453278550E0, + 0.001042698696559292759051590242E0, + -0.044445273029113458906055765365E0, + -0.004670751812662861209726508477E0, + -0.015613390485814379318605247424E0, + -0.030581651696100000521074498679E0, + 0.010801113204340588798240297593E0, + -0.015613390485814379318605247424E0, + -0.030581651696100000521074498679E0, + 0.010801113204340588798240297593E0/ C C Weights of second degree 1 null rule. C DATA (W(9,I),I=1,13)/0.009373028261842556370231264134E0, + -0.074249368848508554545399978725E0, + 0.014709707700258308001897299938E0, + 0.009538502545163567494354463302E0, + -0.014268362488069444905870465047E0, + 0.040126396495352694403045023109E0, + 0.028737181842214741174950928350E0, + -0.031618075834734607275229608099E0, + 0.016879961075872039084307382161E0, + 0.010878914758683152984395046434E0, + -0.031618075834734607275229608099E0, + 0.016879961075872039084307382161E0, + 0.010878914758683152984395046434E0/ C C***FIRST EXECUTABLE STATEMENT SRLTRI C TRES = 50*R1MACH(4) C C Compute area of triangle. C AREA = ABS(VER(1,1)*VER(2,2)-VER(1,2)*VER(2,1)-VER(1,1)*VER(2,3)+ + VER(1,3)*VER(2,1)+VER(1,2)*VER(2,3)-VER(1,3)*VER(2,2))/2 C C Compute contributions from the center of the triangle. C X(1,1) = (VER(1,1)+VER(1,2)+VER(1,3))/3 X(2,1) = (VER(2,1)+VER(2,2)+VER(2,3))/3 CALL FUNSUB(X,NUMFUN,RGNERR) DO 20 J = 1,NUMFUN BASVAL(J) = W(1,1)*RGNERR(J) DO 10 K = 1,8 NULL(J,K) = W(K+1,1)*RGNERR(J) 10 CONTINUE 20 CONTINUE C C Compute the contributions from the points with C multiplicity 3. C DO 100 I = 2,WTLENG Z1 = G(1,I) Z2 = G(2,I) Z3 = 1 - Z1 - Z2 X(1,1) = Z1*VER(1,1) + Z2*VER(1,2) + Z3*VER(1,3) X(2,1) = Z1*VER(2,1) + Z2*VER(2,2) + Z3*VER(2,3) X(1,2) = Z2*VER(1,1) + Z3*VER(1,2) + Z1*VER(1,3) X(2,2) = Z2*VER(2,1) + Z3*VER(2,2) + Z1*VER(2,3) X(1,3) = Z3*VER(1,1) + Z1*VER(1,2) + Z2*VER(1,3) X(2,3) = Z3*VER(2,1) + Z1*VER(2,2) + Z2*VER(2,3) DO 90 L = 1,3 CALL FUNSUB(X(1,L),NUMFUN,RGNERR) DO 80 J = 1,NUMFUN BASVAL(J) = BASVAL(J) + W(1,I)*RGNERR(J) DO 70 K = 1,8 NULL(J,K) = NULL(J,K) + W(K+1,I)*RGNERR(J) 70 CONTINUE 80 CONTINUE 90 CONTINUE 100 CONTINUE C C Compute the errors. C GREATE = 0 DO 130 J = 1,NUMFUN DEG7 = AREA*SQRT(NULL(J,1)**2+NULL(J,2)**2) DEG5 = AREA*SQRT(NULL(J,3)**2+NULL(J,4)**2) DEG3 = AREA*SQRT(NULL(J,5)**2+NULL(J,6)**2) DEG1 = AREA*SQRT(NULL(J,7)**2+NULL(J,8)**2) IF (DEG5.NE.0) THEN R1 = DEG7/DEG5 ELSE R1 = 1 END IF IF (DEG3.NE.0) THEN R2 = DEG5/DEG3 ELSE R2 = 1 END IF IF (DEG1.NE.0) THEN R3 = DEG3/DEG1 ELSE R3 = 1 END IF R = MAX(R1,R2,R3) IF (R.GE.1) THEN RGNERR(J) = 10*MAX(DEG1,DEG3,DEG5,DEG7) ELSE IF (R.GE.0.5E0) THEN RGNERR(J) = 10*R*DEG7 ELSE RGNERR(J) = 40* (R**3)*DEG7 END IF BASVAL(J) = AREA*BASVAL(J) NOISE = ABS(BASVAL(J))*TRES C C The following statement is included to set the error to the noise C level if the two error estimates, assumed to be the best, are both C below or on the noise level C IF ((DEG7.LE.NOISE).AND.(DEG5.LE.NOISE)) RGNERR(J)=NOISE RGNERR(J) = MAX(NOISE,RGNERR(J)) IF (RGNERR(J).GT.GREATE) THEN GREATE = RGNERR(J) END IF 130 CONTINUE C C***END SRLTRI C END C SAMPLE PROGRAM C The following program will approximate the integral of C exp(x*x+y*y) C over the triangle (0.,0.),(2.,0.),(0.,2.) and print out the C values of the estimated integral, the estimated error, the number C of function evaluations, and IFAIL. PROGRAM STEST1 INTEGER NUMFUN,NW,MDIV,MINPTS,LENVER,NUMTRI PARAMETER (NUMFUN=1,MDIV=1,NUMTRI=1,MINPTS=37) PARAMETER ( LENVER = 100 ) C C We need to choose LENVER such that: C LENVER >= 3*((MAXPTS-37*NUMTRI)/(4*37)) + NUMTRI C This simply means that we have enough space to achieve MAXPTS C function evaluations. By choosing LENVER bigger than C this lower bound we may later change MAXPTS and RESTAR and restart C the computations from the point where the last run stopped. C The reason for stopping may have been that MAXPTS was too small C to achieve the requested error. C With our choice LENVER = 100 any value of MAXPTS C between 37 and 5068 will be accepted by the code. Choosing C MAXPTS = 2000 allows us to restart with a greater value C if necessary. C PARAMETER (NW = LENVER*(2*NUMFUN+1) + + MAX(32*MDIV,8*NUMTRI)*NUMFUN + 1) REAL VER(2,3,LENVER),RESULT(NUMFUN), + ABSERR(NUMFUN),WORK(NW),EPSABS,EPSREL INTEGER RESTAR,NEVAL,IFAIL,IWORK(LENVER+MDIV),MAXPTS EXTERNAL F VER(1,1,1) = 0.E0 VER(1,2,1) = 2.E0 VER(1,3,1) = 0.E0 VER(2,1,1) = 0.E0 VER(2,2,1) = 0.E0 VER(2,3,1) = 2.E0 EPSABS = 0.E0 EPSREL = 1.E-5 RESTAR = 0 MAXPTS = 2000 CALL SCUTRI(F,NUMFUN,VER,NUMTRI,MINPTS,MAXPTS,EPSABS, + EPSREL,LENVER,NW,RESTAR,RESULT,ABSERR,NEVAL, + IFAIL,WORK,IWORK) WRITE(*,*)'RESULT=',RESULT(1) WRITE(*,*)'ERROR =',ABSERR(1) WRITE(*,*)'NEVAL =',NEVAL WRITE(*,*)'IFAIL =',IFAIL END SUBROUTINE F(X,NUMFUN,FUNVLS) INTEGER NUMFUN REAL X(2),FUNVLS(NUMFUN) FUNVLS(1) = EXP(X(1)*X(1)+X(2)*X(2)) RETURN END C C Output produced on a SUN SPARC station 1. C C RESULT= 11.1813 C ERROR = 4.71008E-06 C NEVAL = 185 C IFAIL = 0 PROGRAM STEST2 C C STEST2 tests some of the features of SCUTRI. C STEST2 checks that SCUTRI integrates to machine C precision all monomials of degree less or equal to 13. C STEST2 checks that SCUTRI integrates correctly a vector C of integrands. C STEST2 checks that the restart feature of SCUTRI works. C EXTERNAL FTESTP,FTESTQ,FTESTO C C SCUTRI TEST PROGRAM C INTEGER NW,IWORK(501) PARAMETER (NW=50000) REAL VERTCE(2,3,500),WRKSTR(NW),ABSERR, + VERTCQ(2,3,500) VERTCE(1,1,1) = 0 VERTCE(2,1,1) = 0 VERTCE(1,2,1) = 1 VERTCE(2,2,1) = 0 VERTCE(1,3,1) = 0 VERTCE(2,3,1) = 1 VERTCE(1,1,2) = 1 VERTCE(2,1,2) = 1 VERTCE(1,2,2) = 1 VERTCE(2,2,2) = 0 VERTCE(1,3,2) = 0 VERTCE(2,3,2) = 1 VERTCQ(1,1,1) = -1 VERTCQ(2,1,1) = 0 VERTCQ(1,2,1) = 0.5E0 VERTCQ(2,2,1) = -SQRT(3.E0)/2 VERTCQ(1,3,1) = 0.5E0 VERTCQ(2,3,1) = SQRT(3.E0)/2 ABSERR = 0.6 C C TEST FOR INTEGRATING MONOMIALS C Selected monomials, degrees 0-13 C C We first check the 105 monomials in cartesian C coordinates which SCUTRI should integrate to machine C precision. C CALL ATEST(VERTCE,1,500,1000,105,FTESTP,ABSERR,NW,0,WRKSTR, + IWORK) C C Then the same for the 21 monomials in polar coordinates. C CALL ATEST(VERTCQ,1,500,1000,21,FTESTQ,ABSERR,NW,0,WRKSTR, + IWORK) C C Test restart on a peak problem. C VERTCE(1,1,1) = 0 VERTCE(2,1,1) = 0 VERTCE(1,2,1) = 1 VERTCE(2,2,1) = 0 VERTCE(1,3,1) = 0 VERTCE(2,3,1) = 1 VERTCE(1,1,2) = 1 VERTCE(2,1,2) = 1 VERTCE(1,2,2) = 1 VERTCE(2,2,2) = 0 VERTCE(1,3,2) = 0 VERTCE(2,3,2) = 1 ABSERR = 0.1 CALL ATEST(VERTCE,2,500,74,1,FTESTO,ABSERR,NW,0,WRKSTR, + IWORK) ABSERR = 1D-3 CALL ATEST(VERTCE,2,500,1000,1,FTESTO,ABSERR,NW,1,WRKSTR, + IWORK) CALL ATEST(VERTCE,2,500,2000,1,FTESTO,ABSERR,NW,1,WRKSTR, + IWORK) END SUBROUTINE ATEST(VERTCE,NUMTRI,LENVER,MAXCLS,NFUN,TSTSUB,ABSERR, + LENWRK,IREST,WRKSTR,IWORK) EXTERNAL TSTSUB INTEGER LENWRK,IREST,NEVAL,NUMTRI,LENVER,IWORK(LENVER+1) REAL VERTCE(2,3,LENVER),ABSEST(105),FINEST(105), + WRKSTR(LENWRK),ABSERR,REL SAVE NEVAL,ABSEST,FINEST INTEGER N,MAXCLS,NFUN,IFAIL REL = 0 CALL SCUTRI(TSTSUB,NFUN,VERTCE,NUMTRI,0,MAXCLS,ABSERR,REL,LENVER, + LENWRK,IREST,FINEST,ABSEST,NEVAL,IFAIL,WRKSTR,IWORK) WRITE (*,99999) 99999 FORMAT (' '/5X,'SCUTRI TEST ') WRITE (*,99998) NEVAL,IFAIL 99998 FORMAT (5X,'SUBROUTINE CALLS = ',I6,', IFAIL = ',I2) WRITE (*,99997) 99997 FORMAT (7X,'N',3X,'ABSOLUTE ERROR',4X,'INTEGRAL') DO 10 N = 1,NFUN WRITE (*,99996) N,ABS(FINEST(N)-1),FINEST(N) 99996 FORMAT (6X,I3,E14.2,F21.16) 10 CONTINUE END SUBROUTINE FTESTP(Z,NFUN,F) C C Selected monomials, degree 0-13 C INTEGER ISAVE,NFUN,M,K,L,INDEX REAL Z(2),F(NFUN),EX(120),EXACT SAVE ISAVE,EX DATA ISAVE/0/ IF (ISAVE.EQ.0) THEN DO 100 M = 0,13 DO 90 K = 0,M L = M - K INDEX = (M* (M+1))/2 + K + 1 EX(INDEX) = EXACT(K,L) 90 CONTINUE 100 CONTINUE ISAVE = 1 END IF DO 200 M = 0,13 DO 190 K = 0,M L = M - K INDEX = (M* (M+1))/2 + K + 1 F(INDEX) = Z(1)**K*Z(2)**L/EX(INDEX) 190 CONTINUE 200 CONTINUE END SUBROUTINE FTESTQ(Z,NFUN,F) C C Selected monomials in polar coordinates, degree 0-13 C INTEGER NFUN,J,K,I,L REAL Z(2),F(NFUN),M(0:15,0:15),R,ALPHA,PI PI = 3.1415926535897932384E0 M(0,0) = 1.000000000000000000000000000000E0 M(2,0) = 0.250000000000000000000000000000E0 M(4,0) = 0.099999999999999999999999999997E0 M(6,0) = 0.051785714285714285714285714284E0 M(8,0) = 0.031428571428571428571428571427E0 M(10,0) = 0.021103896103896103896103896103E0 M(12,0) = 0.015163408020550877693734836591E0 M(14,0) = 0.011429195804195804195804195802E0 M(3,3) = -0.099999999999999999999999999997E0 M(5,3) = -0.057142857142857142857142857141E0 M(7,3) = -0.035714285714285714285714285713E0 M(9,3) = -0.024025974025974025974025974024E0 M(11,3) = -0.017132867132867132867132867132E0 M(13,3) = -0.012787212787212787212787212785E0 M(15,3) = -0.009891579009226068049597461361E0 M(6,6) = 0.035714285714285714285714285713E0 M(8,6) = 0.024999999999999999999999999997E0 M(10,6) = 0.018181818181818181818181818181E0 M(12,6) = 0.013686313686313686313686313685E0 M(14,6) = 0.010614385614385614385614385614E0 M(9,9) = -0.018181818181818181818181818181E0 M(11,9) = -0.013986013986013986013986013985E0 M(13,9) = -0.010989010989010989010989010988E0 M(15,9) = -0.008807369101486748545572074984E0 M(12,12) = 0.010989010989010989010989010988E0 M(14,12) = 0.008928571428571428571428571429E0 M(15,15) = -0.007352941176470588235294117647E0 I = 0 R = SQRT(ABS(Z(1))**2+ABS(Z(2))**2) IF (R.EQ.0) THEN ALPHA = 0 ELSE ALPHA = ACOS(Z(1)/R) END IF IF (Z(2).LT.0) THEN ALPHA = 2*PI - ALPHA END IF DO 220 J = 0,13 DO 210 K = 0,J/3 IF (MOD(J+3*K,2).EQ.0) THEN I = I + 1 IF (R.EQ.0) THEN IF (J.EQ.0) THEN F(I) = 4/ (3*SQRT(3.E0)) ELSE F(I) = 0 END IF ELSE IF (J.EQ.0) THEN F(I) = 1 ELSE F(I) = 1 DO 205 L = 1,J F(I) = F(I)*R 205 CONTINUE END IF F(I) = F(I)*COS(3*K*ALPHA) F(I) = 4*F(I)/ (3*SQRT(3.E0)*M(J,3*K)) END IF END IF 210 CONTINUE 220 CONTINUE END SUBROUTINE FTESTO(Z,NFUN,F) C C Product Peak C INTEGER NFUN REAL Z(2),F(NFUN) REAL ALPHA1,ALPHA2,BETA1,BETA2,TINT F(1) = 1.E0 ALPHA1 = 10 ALPHA2 = 10 BETA1 = 0.2E0 BETA2 = 0.3E0 TINT = 1 TINT = TINT* (ATAN((1-BETA1)*ALPHA1)-ATAN( - BETA1*ALPHA1))* + ALPHA1 TINT = TINT* (ATAN((1-BETA2)*ALPHA2)-ATAN( - BETA2*ALPHA2))* + ALPHA2 F(1) = F(1)/ (1.E0/ALPHA1**2+ (Z(1)-BETA1)**2) F(1) = F(1)/ (1.E0/ALPHA2**2+ (Z(2)-BETA2)**2) F(1) = F(1)/TINT END REAL FUNCTION EXACT(K,L) INTEGER K,L,I REAL BIN,BINOM EXACT = 0 DO 100 I = 0,L + 1 BIN = BINOM(L+1,I) EXACT = EXACT + (-1)** (I)*BIN/DBLE(I+K+1) 100 CONTINUE EXACT = EXACT/DBLE(L+1) RETURN END REAL FUNCTION BINOM(N,R) INTEGER N,R,I BINOM = 1 IF (R.EQ.0 .OR. R.EQ.N) THEN GO TO 999 END IF DO 10 I = N,N - R + 1,-1 BINOM = BINOM*I 10 CONTINUE DO 20 I = 1,R BINOM = BINOM/I 20 CONTINUE 999 RETURN END C C Output produced on a SUN SPARC station 1. C C C SCUTRI TEST C SUBROUTINE CALLS = 925, IFAIL = 0 C N ABSOLUTE ERROR INTEGRAL C 1 0.60E-07 0.9999999403953552 C 2 0.00E+00 1.0000000000000000 C 3 0.00E+00 1.0000000000000000 C 4 0.00E+00 1.0000000000000000 C 5 0.00E+00 1.0000000000000000 C 6 0.00E+00 1.0000000000000000 C 7 0.00E+00 1.0000000000000000 C 8 0.60E-07 0.9999999403953552 C 9 0.24E-06 0.9999997615814209 C 10 0.00E+00 1.0000000000000000 C 11 0.12E-06 0.9999998807907104 C 12 0.72E-06 0.9999992847442627 C 13 0.12E-06 0.9999998807907104 C 14 0.24E-06 1.0000002384185791 C 15 0.60E-07 0.9999999403953552 C 16 0.00E+00 1.0000000000000000 C 17 0.14E-05 0.9999985694885254 C 18 0.22E-05 0.9999977946281433 C 19 0.77E-06 0.9999992251396179 C 20 0.72E-06 0.9999992847442627 C 21 0.12E-06 0.9999998807907104 C 22 0.12E-06 0.9999998807907104 C 23 0.00E+00 1.0000000000000000 C 24 0.36E-05 0.9999963641166687 C 25 0.42E-05 1.0000041723251343 C 26 0.21E-05 1.0000021457672119 C 27 0.12E-05 0.9999988079071045 C 28 0.42E-06 0.9999995827674866 C 29 0.22E-05 0.9999977946281433 C 30 0.17E-04 1.0000172853469849 C 31 0.21E-05 1.0000021457672119 C 32 0.23E-04 0.9999770522117615 C 33 0.39E-05 0.9999960660934448 C 34 0.77E-05 0.9999923110008240 C 35 0.15E-05 0.9999984502792358 C 36 0.00E+00 1.0000000000000000 C 37 0.19E-05 1.0000019073486328 C 38 0.26E-04 1.0000257492065430 C 39 0.16E-04 0.9999836087226868 C 40 0.48E-05 0.9999952316284180 C 41 0.26E-04 1.0000263452529907 C 42 0.11E-04 1.0000113248825073 C 43 0.75E-05 0.9999925494194031 C 44 0.60E-06 1.0000005960464478 C 45 0.18E-06 0.9999998211860657 C 46 0.60E-07 0.9999999403953552 C 47 0.21E-04 1.0000211000442505 C 48 0.49E-04 1.0000492334365845 C 49 0.96E-04 1.0000958442687988 C 50 0.48E-05 1.0000047683715820 C 51 0.26E-04 0.9999737739562988 C 52 0.39E-04 1.0000392198562622 C 53 0.26E-05 0.9999974370002747 C 54 0.12E-05 1.0000011920928955 C 55 0.18E-06 0.9999998211860657 C 56 0.13E-04 1.0000134706497192 C 57 0.18E-03 0.9998213052749634 C 58 0.40E-03 1.0004028081893921 C 59 0.31E-03 1.0003107786178589 C 60 0.44E-03 1.0004438161849976 C 61 0.13E-03 0.9998712539672852 C 62 0.91E-04 1.0000914335250854 C 63 0.50E-04 1.0000503063201904 C 64 0.58E-05 1.0000058412551880 C 65 0.23E-05 0.9999976754188538 C 66 0.48E-06 0.9999995231628418 C 67 0.00E+00 1.0000000000000000 C 68 0.25E-03 1.0002484321594238 C 69 0.14E-02 1.0013911724090576 C 70 0.21E-03 1.0002107620239258 C 71 0.10E-03 0.9998965263366699 C 72 0.11E-02 0.9988908171653748 C 73 0.27E-03 0.9997336268424988 C 74 0.27E-03 1.0002706050872803 C 75 0.18E-03 1.0001844167709351 C 76 0.22E-05 0.9999977946281433 C 77 0.44E-05 0.9999955892562866 C 78 0.36E-06 0.9999996423721313 C 79 0.46E-04 0.9999543428421021 C 80 0.10E-03 0.9999003410339355 C 81 0.13E-02 1.0012675523757935 C 82 0.20E-02 1.0019829273223877 C 83 0.42E-02 0.9957622885704041 C 84 0.45E-02 1.0045202970504761 C 85 0.20E-02 1.0020430088043213 C 86 0.40E-03 1.0004023313522339 C 87 0.18E-03 0.9998205900192261 C 88 0.16E-03 0.9998353123664856 C 89 0.15E-04 0.9999846220016479 C 90 0.46E-05 0.9999954104423523 C 91 0.66E-06 0.9999993443489075 C 92 0.15E-03 0.9998496174812317 C 93 0.99E-03 0.9990144968032837 C 94 0.47E-03 1.0004729032516479 C 95 0.50E-02 0.9950026273727417 C 96 0.73E-02 1.0072541236877441 C 97 0.10E-01 1.0101681947708130 C 98 0.35E-04 1.0000345706939697 C 99 0.16E-02 0.9983593225479126 C 100 0.21E-02 1.0020525455474854 C 101 0.23E-03 0.9997737407684326 C 102 0.18E-03 0.9998184442520142 C 103 0.60E-04 0.9999404549598694 C 104 0.54E-05 0.9999946355819702 C 105 0.77E-06 0.9999992251396179 C C SCUTRI TEST C SUBROUTINE CALLS = 925, IFAIL = 0 C N ABSOLUTE ERROR INTEGRAL C 1 0.24E-06 1.0000002384185791 C 2 0.00E+00 1.0000000000000000 C 3 0.00E+00 1.0000000000000000 C 4 0.00E+00 1.0000000000000000 C 5 0.00E+00 1.0000000000000000 C 6 0.12E-06 1.0000001192092896 C 7 0.12E-06 0.9999998807907104 C 8 0.00E+00 1.0000000000000000 C 9 0.12E-06 1.0000001192092896 C 10 0.00E+00 1.0000000000000000 C 11 0.00E+00 1.0000000000000000 C 12 0.12E-06 0.9999998807907104 C 13 0.24E-06 1.0000002384185791 C 14 0.00E+00 1.0000000000000000 C 15 0.00E+00 1.0000000000000000 C 16 0.00E+00 1.0000000000000000 C 17 0.24E-06 1.0000002384185791 C 18 0.12E-06 1.0000001192092896 C 19 0.24E-06 0.9999997615814209 C 20 0.12E-06 1.0000001192092896 C 21 0.18E-06 0.9999998211860657 C C SCUTRI TEST C SUBROUTINE CALLS = 74, IFAIL = 1 C N ABSOLUTE ERROR INTEGRAL C 1 0.13E-01 0.9869416356086731 C C SCUTRI TEST C SUBROUTINE CALLS = 962, IFAIL = 1 C N ABSOLUTE ERROR INTEGRAL C 1 0.29E-02 0.9971073269844055 C C SCUTRI TEST C SUBROUTINE CALLS = 1702, IFAIL = 0 C N ABSOLUTE ERROR INTEGRAL C 1 0.14E-05 1.0000014305114746 SUBROUTINE STRTRI(DVFLAG,SBRGNS,GREATE,LIST,NEW) C***BEGIN PROLOGUE STRTRI C***REFER TO SCUTRI C***PURPOSE STRTRI maintains a heap of subregions. C***DESCRIPTION STRTRI maintains a heap of subregions. C The subregions are stored in a partially sorted C binary tree, ordered according to the size of the C greatest error estimates of each subregion(GREATE). C The subregion with greatest error estimate is in the C first position of the heap. C C PARAMETERS C C DVFLAG Integer. C If DVFLAG = 1, we remove the subregion with C greatest error from the heap. C If DVFLAG = 2, we insert a new subregion in the heap. C SBRGNS Integer. C Number of subregions in the heap. C GREATE Real array of dimension SBRGNS. C Used to store the greatest estimated errors in C all subregions. C LIST Integer array of dimension SBRGNS. C Used as a partially ordered list of pointers to the C different subregions. This list is a heap where the C element on top of the list is the subregion with the C greatest error estimate. C NEW Integer. C Index to the new region to be inserted in the heap. C***ROUTINES CALLED-NONE C***END PROLOGUE STRTRI C C Global variables. C INTEGER DVFLAG,NEW,SBRGNS,LIST(*) REAL GREATE(*) C C Local variables. C C GREAT is used as intermediate storage for the greatest error of a C subregion. C SUBRGN Position of child/parent subregion in the heap. C SUBTMP Position of parent/child subregion in the heap. INTEGER SUBRGN,SUBTMP REAL GREAT C C***FIRST EXECUTABLE STATEMENT STRTRI C C If DVFLAG = 1, we will reduce the partial ordered list by the C element with greatest estimated error. Thus the element in C in the heap with index LIST(1) is vacant and can be used later. C Reducing the heap by one element implies that the last element C should be re-positioned. C IF (DVFLAG.EQ.1) THEN GREAT = GREATE(LIST(SBRGNS)) SBRGNS = SBRGNS - 1 SUBRGN = 1 20 SUBTMP = 2*SUBRGN IF (SUBTMP.LE.SBRGNS) THEN IF (SUBTMP.NE.SBRGNS) THEN C C Find max. of left and right child. C IF (GREATE(LIST(SUBTMP)).LT. + GREATE(LIST(SUBTMP+1))) THEN SUBTMP = SUBTMP + 1 END IF END IF C C Compare max.child with parent. C If parent is max., then done. C IF (GREAT.LT.GREATE(LIST(SUBTMP))) THEN C C Move the pointer at position SUBTMP up the heap. C LIST(SUBRGN) = LIST(SUBTMP) SUBRGN = SUBTMP GO TO 20 END IF END IF C C Update the pointer. C IF (SBRGNS.GT.0) THEN LIST(SUBRGN) = LIST(SBRGNS+1) END IF ELSE IF (DVFLAG.EQ.2) THEN C C If DVFLAG = 2, find the position for the NEW region in the heap. C GREAT = GREATE(NEW) SUBRGN = SBRGNS 40 SUBTMP = SUBRGN/2 IF (SUBTMP.GE.1) THEN C C Compare max.child with parent. C If parent is max, then done. C IF (GREAT.GT.GREATE(LIST(SUBTMP))) THEN C C Move the pointer at position SUBTMP down the heap. C LIST(SUBRGN) = LIST(SUBTMP) SUBRGN = SUBTMP GO TO 40 END IF END IF C C Set the pointer to the new region in the heap. C LIST(SUBRGN) = NEW END IF C C***END STRTRI C RETURN END