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