1C ALGORITHM 584, COLLECTED ALGORITHMS FROM ACM. 2C ALGORITHM APPEARED IN ACM-TRANS. MATH. SOFTWARE, VOL.8, NO. 2, 3C JUN., 1982, P. 210. Author: D.P. Laurie 4C PROGRAM KONYN(OUTPUT,TAPE6=OUTPUT) 5 SUBROUTINE CUBTRI(F, T, EPS, MCALLS, ANS, ERR, NCALLS, W, NW, 6 * IDATA, RDATA, IER) 7c 8c changed on June 6th 2011 in order to allow parallelization: 9c the common statements were commented with cparallel, the data 10C statements were converted. 11C 12C ADAPTIVE CUBATURE OVER A TRIANGLE 13C 14C PARAMETERS 15C F - USER SUPPLIED EXTERNAL FUNCTION OF THE FORM 16C F(X,Y,IDATA,RDATA) 17C WHERE X AND Y ARE THE CARTESIAN COORDINATES OF A 18C POINT IN THE PLANE, AND IDATA AND RDATA ARE INTEGER 19C AND REAL*8 VECTORS IN WHICH DATA MAY BE PASSED. 20C T - ARRAY OF DIMENSION (2,3) WHERE T(1,J) AND T(2,J) 21C ARE THE X AND Y COORDINATES OF THE J-TH VERTEX OF 22C THE GIVEN TRIANGLE (INPUT) 23C EPS - REQUIRED TOLERANCE (INPUT). IF THE COMPUTED 24C INTEGRAL IS BETWEEN-1 AND 1, AN ABSOLUTE ERROR 25C TEST IS USED, ELSE A RELATIVE ERROR TEST IS USED. 26C MCALLS- MAXIMUM PERMITTED NUMBER OF CALLS TO F (INPUT) 27C ANS - ESTIMATE FOR THE VALUE OF THE INTEGRAL OF F OVER 28C THE GIVEN TRIANGLE (OUTPUT) 29C ERR - ESTIMATED ABSOLUTE ERROR IN ANS (OUTPUT) 30C NCALLS- ACTUAL NUMBER OF CALLS TO F (OUTPUT). THIS 31C PARAMETER MUST BE INITIALIZED TO 0 ON THE FIRST 32C CALL TO CUBTRI FOR A GIVEN INTEGRAL (INPUT) 33C W - WORK SPACE. MAY NOT BE DESTROYED BETWEEN CALLS TO 34C CUBTRI IF RESTARTING IS INTENDED 35C NW - LENGTH OF WORK SPACE (INPUT). 36C IF NW .GE. 3*(19+3*MCALLS)/38, TERMINATION DUE TO 37C FULL WORK SPACE WILL NOT OCCUR. 38C IER - TERMINATION INDICATOR (OUTPUT) 39C IER=0 NORMAL TERMINATION, TOLERANCE SATISFIED 40C IER=1 MAXIMUM NUMBER OF CALLS REACHED 41C IER=2 WORK SPACE FULL 42C IER=3 FURTHER SUBDIVISION OF TRIANGLES IMPOSSIBLE 43C IER=4 NO FURTHER IMPROVEMENT IN ACCURACY IS 44C POSSIBLE DUE TO ROUNDING ERRORS IN FUNCTION 45C VALUES 46C IER=5 NO FURTHER IMPROVEMENT IN ACCURACY IS 47C POSSIBLE BECAUSE SUBDIVISION DOES NOT 48C CHANGE THE ESTIMATED INTEGRAL. MACHINE 49C ACCURACY HAS PROBABLY BEEN REACHED BUT 50C THE ERROR ESTIMATE IS NOT SHARP ENOUGH. 51C 52C CUBTRI IS DESIGNED TO BE CALLED REPEATEDLY WITHOUT WASTING 53C EARLIER WORK. THE PARAMETER NCALLS IS USED TO INDICATE TO 54C CUBTRI AT WHAT POINT TO RESTART, AND MUST BE RE-INITIALIZED 55C TO 0 WHEN A NEW INTEGRAL IS TO BE COMPUTED. AT LEAST ONE OF 56C THE PARAMETERS EPS, MCALLS AND NW MUST BE CHANGED BETWEEN 57C CALLS TO CUBTRI, ACCORDING TO THE RETURNED VALUE OF IER. NONE 58C OF THE OTHER PARAMETERS MAY BE CHANGED IF RESTARTING IS DONE. 59C IF IER=3 IS ENCOUNTERED, THERE PROBABLY IS A SINGULARITY 60C SOMEWHERE IN THE REGION. THE ERROR MESSAGE INDICATES THAT 61C FURTHER SUBDIVISION IS IMPOSSIBLE BECAUSE THE VERTICES OF THE 62C SMALLER TRIANGLES PRODUCED WILL BEGIN TO COALESCE TO THE 63C PRECISION OF THE COMPUTER. THIS SITUATION CAN USUALLY BE 64C RELIEVED BY SPECIFYING THE REGION IN SUCH A WAY THAT THE 65C SINGULARITY IS LOCATED AT THE THIRD VERTEX OF THE TRIANGLE. 66C IF IER=4 IS ENCOUNTERED, THE VALUE OF THE INTEGRAL CANNOT BE 67C IMPROVED ANY FURTHER. THE ONLY EXCEPTION TO THIS OCCURS WHEN A 68C FUNCTION WITH HIGHLY IRREGULAR BEHAVIOUR IS INTEGRATED (E.G. 69C FUNCTIONS WITH JUMP DISCONTINUITIES OR VERY HIGHLY OSCILLATORY 70C FUNCTIONS). IN SUCH A CASE THE USER CAN DISABLE THE ROUNDING 71C ERROR TEST BY REMOVING THE IF STATEMENT SHORTLY AFTER STATEMENT 72C NUMBER 70. 73C 74 implicit none 75 EXTERNAL F,rnderr 76 INTEGER IDATA(1), IER, MCALLS, NCALLS, NW,jkp,i,j,k,l,maxc,maxk, 77 & mw,nfe 78 REAL*8 ALFA, ANS, ANSKP, AREA, EPS, ERR, ERRMAX, H, Q1, Q2, R1,R2, 79 * RDATA(1), D(2,4), S(4), T(2,3), VEC(2,3), W(6,NW), X(2),zero, 80 & point5,one,rnderr 81C ACTUAL DIMENSION OF W IS (6,NW/6) 82C 83 REAL*8 TANS, TERR, DZERO 84cparallel COMMON /CUBSTA/ TANS, TERR 85C THIS COMMON IS REQUIRED TO PRESERVE TANS AND TERR BETWEEN CALLS 86C AND TO SAVE VARIABLES IN FUNCTION RNDERR 87 nfe=19 88 s=(/1.d0,1.d0,1.d0,-1.d0/) 89 d=reshape((/0.0d0,0.0d0,0.0d0,1.0d0,1.0d0,0.0d0,1.0d0,1.0d0/), 90 & (/2,4/)) 91 zero=0.d0 92 one=1.d0 93 dzero=0.d0 94 point5=.5d0 95cparallel DATA NFE /19/, S(1), S(2), S(3), S(4) /3*1E0,-1E0/, D(1,1), 96cparallel * D(2,1) /0.0,0.0/, D(1,2), D(2,2) /0.0,1.0/, D(1,3), D(2,3) 97cparallel * /1.0,0.0/, D(1,4), D(2,4) /1.0,1.0/ 98C NFE IS THE NUMBER OF FUNCTION EVALUATIONS PER CALL TO CUBRUL. 99cparallel DATA ZERO /0.E0/, ONE /1.E0/, DZERO /0.D0/, POINT5 /.5E0/ 100C 101C CALCULATE DIRECTION VECTORS, AREA AND MAXIMUM NUMBER 102C OF SUBDIVISIONS THAT MAY BE PERFORMED 103 DO 20 I=1,2 104 VEC(I,3) = T(I,3) 105 DO 10 J=1,2 106 VEC(I,J) = T(I,J) - T(I,3) 107 10 CONTINUE 108 20 CONTINUE 109 MAXC = (MCALLS/NFE+3)/4 110 IER = 1 111 MAXK = MIN0(MAXC,(NW/6+2)/3) 112 IF (MAXC.GT.MAXK) IER = 2 113 AREA = ABS(VEC(1,1)*VEC(2,2)-VEC(1,2)*VEC(2,1))*POINT5 114 K = (NCALLS/NFE+3)/4 115 MW = 3*(K-1) + 1 116 IF (NCALLS.GT.0) GO TO 30 117C 118C TEST FOR TRIVIAL CASES 119 TANS = DZERO 120 TERR = DZERO 121 IF (AREA.EQ.ZERO) GO TO 90 122 IF (MCALLS.LT.NFE) GO TO 100 123 IF (NW.LT.6) GO TO 110 124C 125C INITIALIZE DATA LIST 126 K = 1 127 MW = 1 128 W(1,1) = ZERO 129 W(2,1) = ZERO 130 W(3,1) = ONE 131 CALL CUBRUL(F, VEC, W(1,1), IDATA, RDATA) 132 TANS = W(5,1) 133 TERR = W(6,1) 134 NCALLS = NFE 135C 136C TEST TERMINATION CRITERIA 137 30 ANS = TANS 138 ERR = TERR 139 IF (ERR.LT.DMAX1(ONE,ABS(ANS))*EPS) GO TO 90 140 IF (K.EQ.MAXK) GO TO 120 141C 142C FIND TRIANGLE WITH LARGEST ERROR 143 ERRMAX = ZERO 144 DO 40 I=1,MW 145 IF (W(6,I).LE.ERRMAX) GO TO 40 146 ERRMAX = W(6,I) 147 J = I 148 40 CONTINUE 149C 150C SUBDIVIDE TRIANGLE INTO FOUR SUBTRIANGLES AND UPDATE DATA LIST 151 DO 50 I=1,2 152 X(I) = W(I,J) 153 50 CONTINUE 154 H = W(3,J)*POINT5 155 IF (RNDERR(X(1),H,X(1),H).NE.ZERO) GO TO 130 156 IF (RNDERR(X(2),H,X(2),H).NE.ZERO) GO TO 130 157 ANSKP = (TANS) 158 TANS = TANS - (W(5,J)) 159 TERR = TERR - (W(6,J)) 160 R1 = W(4,J) 161 R2 = W(5,J) 162 JKP = J 163 Q1 = ZERO 164 Q2 = ZERO 165 DO 70 I=1,4 166 DO 60 L=1,2 167 W(L,J) = X(L) + H*D(L,I) 168 60 CONTINUE 169 W(3,J) = H*S(I) 170 CALL CUBRUL(F, VEC, W(1,J), IDATA, RDATA) 171 Q2 = Q2 + W(5,J) 172 Q1 = Q1 + W(4,J) 173 J = MW + I 174 70 CONTINUE 175 ALFA = 1E15 176 IF (Q2.NE.R2) ALFA = ABS((Q1-R1)/(Q2-R2)-ONE) 177 J = JKP 178 DO 80 I=1,4 179 W(6,J) = W(6,J)/ALFA 180 TANS = TANS + W(5,J) 181 TERR = TERR + W(6,J) 182 J = MW + I 183 80 CONTINUE 184 MW = MW + 3 185 NCALLS = NCALLS + 4*NFE 186 K = K + 1 187C 188C IF ANSWER IS UNCHANGED, IT CANNOT BE IMPROVED 189 IF (ANSKP.EQ.(TANS)) GO TO 150 190C 191C REMOVE THIS IF STATEMENT TO DISABLE ROUNDING ERROR TEST 192 IF (K.GT.3 .AND. ABS(Q2-R2).GT.ABS(Q1-R1)) GO TO 140 193 GO TO 30 194C 195C EXITS FROM SUBROUTINE 196 90 IER = 0 197 GO TO 120 198 100 IER = 1 199 GO TO 120 200 110 IER = 2 201 120 ANS = TANS 202 ERR = TERR 203 RETURN 204 130 IER = 3 205 GO TO 120 206 140 IER = 4 207 GO TO 120 208 150 IER = 5 209 GO TO 120 210 END 211 real*8 FUNCTION RNDERR(X, A, Y, B) 212C THIS FUNCTION COMPUTES THE ROUNDING ERROR COMMITTED WHEN THE 213C SUM X+A IS FORMED. IN THE CALLING PROGRAM, Y MUST BE THE SAME 214C AS X AND B MUST BE THE SAME AS A. THEY ARE DECLARED AS 215C DISTINCT VARIABLES IN THIS FUNCTION, AND THE INTERMEDIATE 216C VARIABLES S AND T ARE PUT INTO COMMON, IN ORDER TO DEFEND 217C AGAINST THE WELL-MEANING ACTIONS OF SOME OFFICIOUS OPTIMIZING 218C FORTRAN COMPILERS. 219 implicit none 220 real*8 x,a,y,b,s,t 221cparallel COMMON /CUBATB/ S, T 222 S = X + A 223 T = S - Y 224 RNDERR = T - B 225 RETURN 226 END 227 SUBROUTINE CUBRUL(F, VEC, P, IDATA, RDATA) 228C 229C BASIC CUBATURE RULE PAIR OVER A TRIANGLE 230C 231C PARAMETERS 232C F - EXTERNAL FUNCTION - SEE COMMENTS TO CUBTRI 233C VEC- MATRIX OF BASE VECTORS AND ORIGIN (INPUT) 234C P - TRIANGLE DESCRIPTION VECTOR OF DIMENSION 6 235C P(1) - TRANSFORMED X COORDINATE OF ORIGIN VERTEX(INPUT) 236C P(2) - TRANSFORMED Y COORDINATE OF ORIGIN VERTEX(INPUT) 237C P(3) - DISTANCE OF OTHER VERTICES IN THE DIRECTIONS 238C OF THE BASE VECTORS (INPUT) 239C P(4) - LESS ACCURATE ESTIMATED INTEGRAL (OUTPUT) 240C P(5) - MORE ACCURATE ESTIMATED INTEGRAL (OUTPUT) 241C P(6) - ABS(P(5)-P(4)) (OUTPUT) 242C 243C CUBRUL EVALUATES A LINEAR COMBINATION OF BASIC INTEGRATION 244C RULES HAVING D3 SYMMETRY. THE AREAL*8 COORDINATES PERTAINING TO 245C THE J-TH RULE ARE STORED IN W(I,J),I=1,2,3. THE CORRESPONDING 246C WEIGHTS ARE W(4,J) AND W(5,J), WITH W(5,J) BELONGING TO THE 247C MORE ACCURATE FORMULA. IF W(1,J).EQ.W(2,J), THE INTEGRATION 248C POINT IS THE CENTROID, ELSE IF W(2,J).EQ.W(3,J), THE EVALUATION 249C POINTS ARE ON THE MEDIANS. IN BOTH CASES ADVANTAGE IS TAKEN OF 250C SYMMETRY TO AVOID REPEATING FUNCTION EVALUATIONS. 251C 252C THE FOLLOWING REAL*8 VARIABLES ARE USED TO AVOID 253C UNNECESSARY ROUNDING ERRORS IN FLOATING POINT ADDITION. 254C THEY MAY BE DECLARED SINGLE PRECISION IF REAL*8 IS 255C NOT AVAILABLE AND FULL ACCURACY IS NOT NEEDED. 256 implicit none 257 REAL*8 A1, A2, S, SN, DZERO, DONE, DTHREE, DSIX,f, 258 & point5,x,y 259 REAL*8 AREA, ORIGIN(2), P(6), RDATA(1), TVEC(2,3), VEC(2,3),W(5,6) 260 INTEGER IDATA(1),nquad,i,j,k 261C 262C W CONTAINS POINTS AND WEIGHTS OF THE INTEGRATION FORMULAE 263C NQUAD - NUMBER OF BASIC RULES USED 264C 265C THIS PARTICULAR RULE IS THE 19 POINT EXTENSION (DEGREE 8) OF 266C THE FAMILIAR 7 POINT RULE (DEGREE 5). 267C 268C SIGMA=SQRT(7) 269C PHI=SQRT(15) 270C W(1,1),W(2,1),W(3,1) = 1/3 271C W(4,1) = 9/40 272C W(5,1) = 7137/62720 - 45*SIGMA/1568 273C W(1,2) = 3/7 + 2*PHI/21 274C W(2,2),W(3,2) = 2/7 - PHI/21 275C W(4,2) = 31/80 - PHI/400 276C W(5,2) = - 9301697/4695040 - 13517313*PHI/23475200 277C + 764885*SIGMA/939008 + 198763*PHI*SIGMA/939008 278C W(*,3) = W(*,2) WITH PHI REPLACED BY -PHI 279C W(1,5) = 4/9 + PHI/9 + SIGMA/9 - SIGMA*PHI/45 280C W(2,5),W(3,5) = 5/18 - PHI/18 - SIGMA/18 + SIGMA*PHI/90 281C W(4,5) = 0 282C W(5,5) = 102791225/59157504 + 23876225*PHI/59157504 283C - 34500875*SIGMA/59157504 - 9914825*PHI*SIGMA/59157504 284C W(*,4) = W(*,5) WITH PHI REPLACED BY -PHI 285C W(1,6) = 4/9 + SIGMA/9 286C W(2,6) = W(2,4) 287C W(3,6) = W(2,5) 288C W(4,6) = 0 289C W(5,6) = 11075/8064 - 125*SIGMA/288 290C 291 nquad=6 292 w=reshape((/.3333333333333333333333333d0, 293 & .3333333333333333333333333d0,.3333333333333333333333333d0,.225d0, 294 & .3786109120031468330830822d-1,.7974269853530873223980253d0, 295 & .1012865073234563388009874d0,.1012865073234563388009874d0, 296 & .3778175416344814577870518d0, 297 & .1128612762395489164329420d0,.5971587178976982045911758d-1, 298 & .4701420641051150897704412d0,.4701420641051150897704412d0, 299 & .3971824583655185422129482d0, 300 & .2350720567323520126663380d0,.5357953464498992646629509d0, 301 & .2321023267750503676685246d0,.2321023267750503676685246d0, 302 & 0.d0,.3488144389708976891842461d0, 303 & .9410382782311208665596304d0, 304 & .2948086088443956672018481d-1,.2948086088443956672018481d-1, 305 & 0.d0,.4033280212549620569433320d-1,.7384168123405100656112906d0, 306 & .2321023267750503676685246d0,.2948086088443956672018481d-1, 307 & 0.d0,.2250583347313904927138324d0/),(/5,6/)) 308cparallel 309cparallel DATA NQUAD /6/, W(1,1), W(2,1), W(3,1)/3*.33333333333333333333333 310cparallel * 33E0/, W(4,1), W(5,1) /.225E0,.3786109120031468330830822E-1/, 311cparallel * W(1,2), W(2,2), W(3,2) /.7974269853530873223980253E0,2* 312cparallel * .1012865073234563388009874E0/, W(4,2), W(5,2) 313cparallel * /.3778175416344814577870518E0,.1128612762395489164329420E0/, 314cparallel * W(1,3), W(2,3), W(3,3) /.5971587178976982045911758E-1,2* 315cparallel * .4701420641051150897704412E0/, W(4,3), W(5,3) 316cparallel * /.3971824583655185422129482E0,.2350720567323520126663380E0/ 317cparallel DATA W(1,4), W(2,4), W(3,4) /.5357953464498992646629509E0,2* 318cparallel * .2321023267750503676685246E0/, W(4,4), W(5,4) 319cparallel * /0.E0,.3488144389708976891842461E0/, W(1,5), W(2,5), W(3,5) 320cparallel * /.9410382782311208665596304E0,2*.2948086088443956672018481E-1/, 321cparallel * W(4,5), W(5,5) /0.E0,.4033280212549620569433320E-1/, W(1,6), 322cparallel * W(2,6), W(3,6) /.7384168123405100656112906E0, 323cparallel * .2321023267750503676685246E0,.2948086088443956672018481E-1/, 324cparallel * W(4,6), W(5,6) /0.E0,.2250583347313904927138324E0/ 325cparallel 326C 327 dzero=0.d0 328 done=1.d0 329 dthree=3.d0 330 dsix=6.d0 331 point5=.5d0 332cparallel DATA DZERO /0.D0/, DONE /1.D0/, DTHREE /3.D0/, DSIX /6.D0/, 333cparallel * POINT5 /.5E0/ 334C 335C SCALE BASE VECTORS AND OBTAIN AREA 336 DO 20 I=1,2 337 ORIGIN(I) = VEC(I,3) + P(1)*VEC(I,1) + P(2)*VEC(I,2) 338 DO 10 J=1,2 339 TVEC(I,J) = P(3)*VEC(I,J) 340 10 CONTINUE 341 20 CONTINUE 342 AREA = POINT5*ABS(TVEC(1,1)*TVEC(2,2)-TVEC(1,2)*TVEC(2,1)) 343 A1 = DZERO 344 A2 = DZERO 345C 346C COMPUTE ESTIMATES FOR INTEGRAL AND ERROR 347 DO 40 K=1,NQUAD 348 X = ORIGIN(1) + W(1,K)*TVEC(1,1) + W(2,K)*TVEC(1,2) 349 Y = ORIGIN(2) + W(1,K)*TVEC(2,1) + W(2,K)*TVEC(2,2) 350 S = (F(X,Y,IDATA,RDATA)) 351 SN = DONE 352 IF (W(1,K).EQ.W(2,K)) GO TO 30 353 X = ORIGIN(1) + W(2,K)*TVEC(1,1) + W(1,K)*TVEC(1,2) 354 Y = ORIGIN(2) + W(2,K)*TVEC(2,1) + W(1,K)*TVEC(2,2) 355 S = S + (F(X,Y,IDATA,RDATA)) 356 X = ORIGIN(1) + W(2,K)*TVEC(1,1) + W(3,K)*TVEC(1,2) 357 Y = ORIGIN(2) + W(2,K)*TVEC(2,1) + W(3,K)*TVEC(2,2) 358 S = S + (F(X,Y,IDATA,RDATA)) 359 SN = DTHREE 360 IF (W(2,K).EQ.W(3,K)) GO TO 30 361 X = ORIGIN(1) + W(1,K)*TVEC(1,1) + W(3,K)*TVEC(1,2) 362 Y = ORIGIN(2) + W(1,K)*TVEC(2,1) + W(3,K)*TVEC(2,2) 363 S = S + (F(X,Y,IDATA,RDATA)) 364 X = ORIGIN(1) + W(3,K)*TVEC(1,1) + W(1,K)*TVEC(1,2) 365 Y = ORIGIN(2) + W(3,K)*TVEC(2,1) + W(1,K)*TVEC(2,2) 366 S = S + (F(X,Y,IDATA,RDATA)) 367 X = ORIGIN(1) + W(3,K)*TVEC(1,1) + W(2,K)*TVEC(1,2) 368 Y = ORIGIN(2) + W(3,K)*TVEC(2,1) + W(2,K)*TVEC(2,2) 369 S = S + (F(X,Y,IDATA,RDATA)) 370 SN = DSIX 371 30 S = S/SN 372 A1 = A1 + W(4,K)*S 373 A2 = A2 + W(5,K)*S 374 40 CONTINUE 375 P(4) = (A1)*AREA 376 P(5) = (A2)*AREA 377 P(6) = ABS(P(5)-P(4)) 378 RETURN 379 END 380