← Back to branch summary

~maddevelopers/mg5amcnlo/2.9.4

« back to all changes in this revision

Viewing changes to vendor/CutTools/src/qcdloop/aaxex.f

pass to v2.0.0

Show diffs side-by-side

added added

removed removed

Lines of Context:
vendor/CutTools/src/qcdloop/aaxex.f
 
1
 
 
2
* file aaxex  23-sep-1990
 
3
 
 
4
*###[ aaxex:
 
5
        subroutine aaxex(ccxi,cdxi,cexi,d0,xmm,xpi,level,ier)
 
6
***#[*comment:***********************************************************
 
7
*                                                                       *
 
8
*       Calculation of four point tensor integrals.  Just a wrapper     *
 
9
*       for ffxdx nowadays, see there for the real description.         *
 
10
*                                                                       *
 
11
*       Input:  xpi(20)   real     the same as in FF                    *
 
12
*               level     integer  rank of tensor integral              *
 
13
*       Output: ccxi(30)  complex  cc0(1),cc1( ),[cc2( ),cc3( ) ]  i,j  *
 
14
*               cdxi(55)  complex  cd0(1),cd1(3),cd2(7),[cd3(13)]       *
 
15
*                                       i=1,2,3,4,5                     *
 
16
*               cexi(35)  complex  ce0(1),ce1(4),ce2(10),ce3(20)        *
 
17
*               ier       integer  FF error flag                        *
 
18
*                                                                       *
 
19
***#]*comment:***********************************************************
 
20
*  #[ declarations:
 
21
        implicit none
 
22
*
 
23
*       arguments
 
24
*
 
25
        integer ier,level
 
26
        DOUBLE PRECISION xpi(20),d0,xmm
 
27
        DOUBLE COMPLEX ccxi(30),cdxi(55),cexi(35)
 
28
*
 
29
*       local variables
 
30
*
 
31
        DOUBLE COMPLEX caxi(5),cbxi(30)
 
32
        DOUBLE PRECISION maxi(5),mbxi(30),mcxi(30),mdxi(55),mexi(35)
 
33
*
 
34
*  #] declarations:
 
35
*  #[ call ffxex:
 
36
*
 
37
        call ffxex(caxi,maxi,cbxi,mbxi,ccxi,mcxi,cdxi,mdxi,cexi,mexi,
 
38
     +          d0,xmm,xpi,level,ier)
 
39
*
 
40
*  #] call ffxex:
 
41
*###] aaxex:
 
42
        end
 
43
*###[ ffxex:
 
44
        subroutine ffxex(caxi,maxi,cbxi,mbxi,ccxi,mcxi,cdxi,mdxi,cexi,
 
45
     +          mexi,d0,xmm,xpi,level,ier)
 
46
***#[*comment:***********************************************************
 
47
*                                                                       *
 
48
*       Calculation of five-point formfactors, notation is defined in   *
 
49
*       adapt.prc and covdec5.prc (WW), leaving out the d_ terms        *
 
50
*       adapted by GJvO 4-apr-1995                                      *
 
51
*                                                                       *
 
52
*       cexi(1)         E0                                              *
 
53
*       cexi(2-5)       E1i, coeff of pi                                *
 
54
*       cexi(6-15)      E2i, coeff of p1p1,p2p2,p3p3,p4p4,(p1p2),       *
 
55
*                       (p1p3),(p1p4),(p2p3),(p2p4),(p3p4)              *
 
56
*       cexi(16-35)     E2i, coeff of p1p1p1,p2p2p2,p3p3p3,p4p4p4,      *
 
57
*                       (p1p1p2),(p1p1p3),(p1p1p4),(p2p2p1),(p2p2p3),   *
 
58
*                       (p2p2p4),(p3p3p1),(p3p3p2),(p3p3p4),(p4p4p1),   *
 
59
*                       (p4p4p2),(p4p4p3),(p1p2p3),(p1p2p4),(p1p3p4),   *
 
60
*                       (p2p3p4)                                        *
 
61
*       cdxi(55)        5*(D0(1),D1i(3),D2i(7)), s_i missing            *
 
62
*       ccxi(30)        15*(C0(1),C1i(2))  ??????                       *
 
63
*       cbxi(30)        not used yet                                    *
 
64
*       caxi(5)         not used yet                                    *
 
65
*       m[abcde]xi()    if c[abcde]xi() were written as a sum of stable *
 
66
*                       terms, the largest term in that sum, i.e., the  *
 
67
*                       accuracy of c[abcde]xi() is precc*m[abcde]xi    *
 
68
*                                                                       *
 
69
*       Input:  xpi(20)   real     the same as in FF                    *
 
70
*               level     integer  rank of tensor integral              *
 
71
*       Output: ccxi(30)  complex  cc0(1),cc1( ),[cc2( ),cc3( ) ]  i,j  *
 
72
*               cdxi(55)  complex  cd0(1),cd1(3),cd2(7),[cd3(13)]       *
 
73
*                                       i=1,2,3,4,5                     *
 
74
*               cexi(35)  complex  ce0(1),ce1(5),ce2(10),ce3(20)        *
 
75
*               ier       integer  FF error flag                        *
 
76
*                                                                       *
 
77
***#]*comment:***********************************************************
 
78
*  #[ declarations:
 
79
        implicit none
 
80
*
 
81
*       arguments
 
82
*
 
83
        integer ier,level
 
84
        DOUBLE PRECISION xpi(20),d0,xmm
 
85
        DOUBLE COMPLEX caxi(5),cbxi(30),ccxi(30),cdxi(55),cexi(35)
 
86
        DOUBLE PRECISION maxi(5),mbxi(30),mcxi(30),mdxi(55),mexi(35)
 
87
*
 
88
*       local variables
 
89
*
 
90
        integer i,j,dl,iplace(11,5),ier0,ier1
 
91
        DOUBLE PRECISION xpj(13),absc,big
 
92
        DOUBLE COMPLEX cd0i(5),cdxj(120),ccxj(140),cbxj(60),caxj(20),cc
 
93
        DOUBLE PRECISION mdxj(120),mcxj(140),mbxj(60),maxj(20)
 
94
        save iplace
 
95
*
 
96
*       common blocks
 
97
*
 
98
        include 'ff.h'
 
99
        include 'aa.h'
 
100
*
 
101
*       statement functions
 
102
*
 
103
        absc(cc) = abs(DBLE(cc)) + abs(DIMAG(cc))
 
104
*
 
105
*       data
 
106
*
 
107
        data iplace /
 
108
     +          2,3,4,5, 07,08,09,15, 12,13, 17,
 
109
     +          1,3,4,5, 11,08,09,10, 14,13, 18,
 
110
     +          1,2,4,5, 06,12,09,10, 14,15, 19,
 
111
     +          1,2,3,5, 06,07,13,10, 11,15, 20,
 
112
     +          1,2,3,4, 06,07,08,14, 11,12, 16/
 
113
*
 
114
*  #] declarations:
 
115
*  #[ init:
 
116
*
 
117
*       initialize to something ridiculous so that one immediately 
 
118
*       notices when it is accidentally used...
 
119
*
 
120
        big = 1/(1d20*xclogm)
 
121
        do i=1,5
 
122
            caxi(i) = big
 
123
        enddo
 
124
        do i=1,30
 
125
            cbxi(i) = big
 
126
        enddo
 
127
        do i=1,30
 
128
            ccxi(i) = big
 
129
        enddo
 
130
        do i=1,55
 
131
            cdxi(i) = big
 
132
        enddo
 
133
        do i=1,35
 
134
            cexi(i) = big
 
135
        enddo
 
136
*
 
137
*  #] init:
 
138
*  #[ level 0: E0, and kinematical quantities for 5 point PV-red
 
139
*       E0-function (ff)
 
140
*
 
141
        ldot = .TRUE.
 
142
        ier1 = ier
 
143
        call ffxe0(cexi(1),cd0i,xpi,ier1)
 
144
        mexi(1) = absc(cexi(1))*DBLE(10)**mod(ier1,50)
 
145
        do i=1,5
 
146
            cdxi(1+11*(i-1)) = cd0i(i)
 
147
        enddo
 
148
*
 
149
        if ( awrite ) then
 
150
            print *,'    '
 
151
            print *,'aaxex : level 0, imported from ff '
 
152
            print *,'E0    = ',cexi(1),mexi(1)
 
153
            print *,'D0(1) = ',cd0i(1)
 
154
            print *,'D0(2) = ',cd0i(2)
 
155
            print *,'D0(3) = ',cd0i(3)
 
156
            print *,'D0(4) = ',cd0i(4)
 
157
            print *,'D0(5) = ',cd0i(5)
 
158
            print *,'xpi used:'
 
159
            do i =1,15
 
160
                print *,i,xpi(i)
 
161
            enddo
 
162
            print *,'imported stuff via ff.h:'
 
163
            print *,'    kin determinat = ',fdel4
 
164
            print *,'dotpr(1,1)= ',fpij5(6,6)
 
165
            print *,'dotpr(2,2)= ',fpij5(7,7)
 
166
            print *,'dotpr(1,2)= ',fpij5(6,7)
 
167
        endif
 
168
*
 
169
        if (level .eq. 0) return
 
170
*
 
171
*  #] level 0: E0, and kinematical quantities for 5 point PV-red
 
172
*  #[ need D-functions till d-level=(level-1):
 
173
        dl=level-1
 
174
*
 
175
*       go trough the 5 different cancellation patterns
 
176
*
 
177
        if ( awrite ) then
 
178
           print *,'     '
 
179
           print *,'------>underlying D-functions up to level:',dl
 
180
        endif
 
181
        xpj(12) = 0
 
182
        xpj(13) = 0
 
183
        do j=1,5
 
184
*           D(j) is the D0 with leg j missing.
 
185
            do i=1,11
 
186
                xpj(i) = xpi(iplace(i,j))
 
187
            enddo
 
188
*           note that we recompute the D0 (or get it from cache)
 
189
            ier0 = ier
 
190
            call ffxdx(caxj(1+4*(j-1)),maxj(1+4*(j-1)),cbxj(1+12*(j-1))
 
191
     +          ,mbxj(1+12*(j-1)),ccxj(1+28*(j-1)),mcxj(1+28*(j-1)),
 
192
     +          cdxj(1+24*(j-1)),mdxj(1+24*(j-1)),d0,xmm,xpj,dl,ier0)
 
193
            ier1 = max(ier1,ier0)
 
194
        enddo
 
195
        ier = ier1
 
196
        if ( awrite ) then
 
197
            print *,' '
 
198
            print *,'ier = ',ier
 
199
            print *,'---->end of D-function output--------------------'
 
200
        endif
 
201
        if ( atest ) then
 
202
*           (although these should come from cache) (but don't yet!)
 
203
            do i=1,5
 
204
                if ( xloss*10d0**(-mod(ier,50))*absc(cd0i(i)-
 
205
     +                  cdxj(1+24*(i-1))) .gt. precx*absc(cd0i(i)) )
 
206
     +                  then
 
207
                    print *,'aaxex: error: D0 does not agree with ',
 
208
     +                  'recomputed: ',cd0i(i),cdxj(1+24*(i-1)),ier
 
209
                endif
 
210
            enddo
 
211
        endif
 
212
*
 
213
*  #] need D-functions till d-level=(level-1):
 
214
*  #[ break to let ffzez tie in:
 
215
*
 
216
*       convert ??xj to ??xi
 
217
*
 
218
        call ffeji(cdxi,mcxi,ccxi,mcxi,cbxi,mbxi,caxi,maxi,
 
219
     +          cdxj,mdxj,ccxj,mcxj,cbxj,mbxj,caxj,maxj,level)
 
220
*
 
221
*       and call the real routine for the rest
 
222
*
 
223
        call ffxexp(caxi,maxi,cbxi,mbxi,ccxi,mcxi,cdxi,mdxi,cexi,mexi,
 
224
     +          xpi,level,ier)
 
225
*
 
226
*  #] break to let ffzez tie in:
 
227
        end
 
228
        subroutine ffxexp(caxi,maxi,cbxi,mbxi,ccxi,mcxi,cdxi,mdxi,cexi,
 
229
     +          mexi,xpi,level,ier)
 
230
*  #[ declarations:
 
231
        implicit none
 
232
*
 
233
*       arguments
 
234
*
 
235
        integer ier,level
 
236
        DOUBLE PRECISION xpi(20),d0,xmm
 
237
        DOUBLE COMPLEX caxi(5),cbxi(30),ccxi(30),cdxi(55),cexi(35)
 
238
        DOUBLE PRECISION maxi(5),mbxi(30),mcxi(30),mdxi(55),mexi(35)
 
239
*
 
240
*       local variables
 
241
*
 
242
        integer i,j,ier0,ier1,ij2k(4,4),m2ijk(3,20)
 
243
        DOUBLE PRECISION xi5(10),f1,f2,f3,f4,absc
 
244
        DOUBLE COMPLEX R(70),cd0i(5),cd1ij(3,5),ce2ij(4,4),ce3ijk(4,4,4)
 
245
     +          ,cd2ijk(3,3,5),cd2i(5),cxy(70),cc,rg(4),cexj(39)
 
246
        DOUBLE PRECISION mdxj(120),mcxj(140),mbxj(60),maxj(20),
 
247
     +          del3ij(5,5),del4i(4)
 
248
        save ij2k,m2ijk
 
249
*
 
250
*       common blocks
 
251
*
 
252
        include 'ff.h'
 
253
        include 'aa.h'
 
254
*
 
255
*       statement functions
 
256
*
 
257
        absc(cc) = abs(DBLE(cc)) + abs(DIMAG(cc))
 
258
*
 
259
*       data
 
260
*
 
261
        data ij2k / 6,10,11,12,
 
262
     +             10, 7,13,14,
 
263
     +             11,13, 8,15,
 
264
     +             12,14,15, 9/
 
265
        data m2ijk /1,1,1,  2,2,2,
 
266
     +              3,3,3,  4,4,4,
 
267
     +              1,1,2,  1,1,3,
 
268
     +              1,1,4,  2,2,1,
 
269
     +              2,2,3,  2,2,4,
 
270
     +              3,3,1,  3,3,2,
 
271
     +              3,3,4,  4,4,1,
 
272
     +              4,4,2,  4,4,3,
 
273
     +              1,2,3,  1,2,4,
 
274
     +              1,3,4,  2,3,4/
 
275
*
 
276
*  #] declarations:
 
277
*  #[ kinematical quatities for 5pv-red:
 
278
*
 
279
*       inverse kinematical matrix xi5  (4X4)
 
280
*
 
281
        call aaxi5(xi5,ier)
 
282
*
 
283
*       AAs f-functions:
 
284
*
 
285
        f1 = 2*fpij5(1,6)
 
286
        f2 = 2*fpij5(1,7)
 
287
        f3 = 2*fpij5(1,8)
 
288
        f4 = 2*fpij5(1,9)
 
289
*
 
290
*  #] kinematical quatities for 5pv-red:
 
291
*  #[ level 1: E11,E12,E13,E14,D0(I)
 
292
        do i=1,5
 
293
            cd0i(i) = cdxi(1+11*(i-1))
 
294
        enddo
 
295
        call ffxe1(cexi(2),cexi(1),del3ij,del4i,cd0i,xpi,fpij5,fdel4,
 
296
     +          ier)
 
297
        if ( atest ) then
 
298
*
 
299
*       PV-reduction
 
300
*
 
301
            R(1) = (f1*cexi(1) + cd0i(2) - cd0i(1))/2
 
302
            R(2) = (f2*cexi(1) + cd0i(3) - cd0i(2))/2
 
303
            R(3) = (f3*cexi(1) + cd0i(4) - cd0i(3))/2
 
304
            R(4) = (f4*cexi(1) + cd0i(5) - cd0i(4))/2
 
305
            cxy(2)=xi5(1)*R(1)+xi5(5)*R(2)+xi5(6) *R(3)+xi5(7) *R(4)
 
306
            cxy(3)=xi5(5)*R(1)+xi5(2)*R(2)+xi5(8) *R(3)+xi5(9) *R(4)
 
307
            cxy(4)=xi5(6)*R(1)+xi5(8)*R(2)+xi5(3) *R(3)+xi5(10)*R(4)
 
308
            cxy(5)=xi5(7)*R(1)+xi5(9)*R(2)+xi5(10)*R(3)+xi5(4) *R(4)
 
309
            do i=2,5
 
310
                if ( xloss*10d0**(-mod(ier,50))*absc(cexi(i)-cxy(i))
 
311
     +                  .gt. precc*absc(cxy(i)) ) then
 
312
                    print *,'aaxex: error: E1 from ffxe1 is not correct'
 
313
     +                  ,i,cexi(i),cxy(i),ier
 
314
                endif
 
315
            enddo
 
316
        endif
 
317
        if ( awrite ) then
 
318
            print *,'     '
 
319
            print *,'aaxex: level 1: id,nevent ',id,nevent
 
320
            print *,'E11=',cexi(2),ier
 
321
            print *,'E12=',cexi(3),ier
 
322
            print *,'E13=',cexi(4),ier
 
323
            print *,'E14=',cexi(5),ier
 
324
        endif
 
325
 
 
326
        if ( level.le.1 ) return
 
327
*
 
328
*  #] level 1:
 
329
*  #[ level 2: E21,E22,E23,E24,E25,E26,E27,E28,E29,E210
 
330
*
 
331
*       D11(I),D12(I),D13(I) need 5 diff D1-functions
 
332
*
 
333
        do i=1,5
 
334
            j = 2 + (i-1)*11
 
335
            cd1ij(1,i)=cdxi(j)
 
336
            cd1ij(2,i)=cdxi(j+1)
 
337
            cd1ij(3,i)=cdxi(j+2)
 
338
        enddo
 
339
*
 
340
*       GJ reduction:
 
341
*
 
342
        call ffxe2(ce2ij,cexi(2),cexi(1),cd1ij,cd0i,xpi,fpij5,
 
343
     +          del3ij,del4i,fdel4,ier)
 
344
        if ( atest ) then
 
345
*
 
346
*       PV-reduction
 
347
*
 
348
        R(11) = (f1*cexi(2) + cd1ij(1,2) + cd0i(1) )/2
 
349
        R(12) = (f2*cexi(2) + cd1ij(1,3) - cd1ij(1,2))/2
 
350
        R(13) = (f3*cexi(2) + cd1ij(1,4) - cd1ij(1,3))/2
 
351
        R(14) = (f4*cexi(2) + cd1ij(1,5) - cd1ij(1,4))/2
 
352
*
 
353
        R(15) = (f1*cexi(3) + cd1ij(1,2) - cd1ij(1,1))/2
 
354
        R(16) = (f2*cexi(3) + cd1ij(2,3) - cd1ij(1,2))/2
 
355
        R(17) = (f3*cexi(3) + cd1ij(2,4) - cd1ij(2,3))/2
 
356
        R(18) = (f4*cexi(3) + cd1ij(2,5) - cd1ij(2,4))/2
 
357
*
 
358
        R(19) = (f1*cexi(4) + cd1ij(2,2) - cd1ij(2,1))/2
 
359
        R(20) = (f2*cexi(4) + cd1ij(2,3) - cd1ij(2,2))/2
 
360
        R(21) = (f3*cexi(4) + cd1ij(3,4) - cd1ij(2,3))/2
 
361
        R(22) = (f4*cexi(4) + cd1ij(3,5) - cd1ij(3,4))/2
 
362
*
 
363
        R(23) = (f1*cexi(5) + cd1ij(3,2) - cd1ij(3,1))/2
 
364
        R(24) = (f2*cexi(5) + cd1ij(3,3) - cd1ij(3,2))/2
 
365
        R(25) = (f3*cexi(5) + cd1ij(3,4) - cd1ij(3,3))/2
 
366
        R(26) = (f4*cexi(5)            - cd1ij(3,4))/2
 
367
*
 
368
        cexi(6) = xi5(1)*R(11)+xi5(5)*R(12)+xi5(6) *R(13)+xi5(7) *R(14)
 
369
        cexi(7) = xi5(5)*R(15)+xi5(2)*R(16)+xi5(8) *R(17)+xi5(9) *R(18)
 
370
        cexi(8) = xi5(6)*R(19)+xi5(8)*R(20)+xi5(3) *R(21)+xi5(10)*R(22)
 
371
        cexi(9) = xi5(7)*R(23)+xi5(9)*R(24)+xi5(10)*R(25)+xi5(4) *R(26)
 
372
        cexi(10)= xi5(5)*R(11)+xi5(2)*R(12)+xi5(8) *R(13)+xi5(9) *R(14)
 
373
        cexi(11)= xi5(6)*R(11)+xi5(8)*R(12)+xi5(3) *R(13)+xi5(10)*R(14)
 
374
        cexi(12)= xi5(7)*R(11)+xi5(9)*R(12)+xi5(10)*R(13)+xi5(4) *R(14)
 
375
        cexi(13)= xi5(6)*R(15)+xi5(8)*R(16)+xi5(3) *R(17)+xi5(10)*R(18)
 
376
        cexi(14)= xi5(7)*R(15)+xi5(9)*R(16)+xi5(10)*R(17)+xi5(4) *R(18)
 
377
        cexi(15)= xi5(7)*R(19)+xi5(9)*R(20)+xi5(10)*R(21)+xi5(4) *R(22)
 
378
*
 
379
        if ( atest ) then
 
380
*
 
381
*           redundancy check
 
382
*
 
383
        cxy(10) = xi5(1)*R(15)+xi5(5)*R(16)+xi5(6)*R(17)+xi5(7) *R(18)
 
384
        cxy(11) = xi5(1)*R(19)+xi5(5)*R(20)+xi5(6)*R(21)+xi5(7) *R(22)
 
385
        cxy(12) = xi5(1)*R(23)+xi5(5)*R(24)+xi5(6)*R(25)+xi5(7) *R(26)
 
386
        cxy(13) = xi5(5)*R(19)+xi5(2)*R(20)+xi5(8)*R(21)+xi5(9) *R(22)
 
387
        cxy(14) = xi5(5)*R(23)+xi5(2)*R(24)+xi5(8)*R(25)+xi5(9) *R(26)
 
388
        cxy(15) = xi5(6)*R(23)+xi5(8)*R(24)+xi5(3)*R(25)+xi5(10)*R(26)
 
389
            do i=10,15
 
390
                if ( absc(cxy(i)-cexi(i))  .gt. .1d-3  ) then
 
391
                    print *,'redundancy check at level 2 failed'
 
392
                endif
 
393
            enddo
 
394
*###] :
 
395
        endif
 
396
*
 
397
*       check against GJ
 
398
*
 
399
        do i=1,4
 
400
            do j=1,4
 
401
                if ( xloss*10d0**(-mod(ier,50))*absc(ce2ij(i,j) -
 
402
     +                  cexi(ij2k(i,j))) .gt. precc*absc(ce2ij(i,j)) )
 
403
     +                  then
 
404
                    print *,'ffxdxp: error: GJ does not agree with PV:',
 
405
     +                  i,j,ce2ij(i,j),cexi(ij2k(i,j)),ce2ij(i,j) -
 
406
     +                  cexi(ij2k(i,j)),ier
 
407
                endif
 
408
            enddo
 
409
        enddo
 
410
        endif
 
411
*
 
412
*       copy to AAs arrays
 
413
*
 
414
        do i=1,4
 
415
            do j=1,4
 
416
                cexi(ij2k(j,i)) = ce2ij(j,i)
 
417
            enddo
 
418
        enddo
 
419
*
 
420
        if ( awrite ) then
 
421
            print *,'     '
 
422
            print *,'aaxex: level 2: id,nevent ',id,nevent
 
423
            print *,'E21 =',cexi(6),ier
 
424
            print *,'E22 =',cexi(7),ier
 
425
            print *,'E23 =',cexi(8),ier
 
426
            print *,'E24 =',cexi(9),ier
 
427
            print *,'E25 =',cexi(10),ier
 
428
            print *,'E26 =',cexi(11),ier
 
429
            print *,'E27 =',cexi(12),ier
 
430
            print *,'E28 =',cexi(13),ier
 
431
            print *,'E28 =',cexi(14),ier
 
432
            print *,'E210=',cexi(15),ier
 
433
        endif
 
434
*
 
435
        if (level .le. 2) return
 
436
*
 
437
*  #] level 2:
 
438
*###[ : level 3 :
 
439
*
 
440
*       first recast the D2's to a more useful form
 
441
*
 
442
        do 15 i=1,5
 
443
            j = 5 +(i-1)*11
 
444
            cd2ijk(1,1,i) = cdxi(j)
 
445
            cd2ijk(2,2,i) = cdxi(j+1)
 
446
            cd2ijk(3,3,i) = cdxi(j+2)
 
447
            cd2ijk(1,2,i) = cdxi(j+3)
 
448
            cd2ijk(2,1,i) = cdxi(j+3)
 
449
            cd2ijk(1,3,i) = cdxi(j+4)
 
450
            cd2ijk(3,1,i) = cdxi(j+4)
 
451
            cd2ijk(2,3,i) = cdxi(j+5)
 
452
            cd2ijk(3,2,i) = cdxi(j+5)
 
453
            cd2i(i) = cdxi(j+6)
 
454
 15     continue
 
455
*
 
456
*       FF function
 
457
*
 
458
        call ffxe3(ce3ijk,ce2ij,cexi(2),cexi(1), cd2ijk,cd2i,cd1ij,cd0i,
 
459
     +          xpi,fpij5, del3ij,del4i,fdel4, ier)
 
460
*
 
461
*       copy FF structs to AA
 
462
*
 
463
        do i=16,35
 
464
            cexi(i) = ce3ijk(m2ijk(1,i),m2ijk(2,i),m2ijk(3,i))
 
465
        enddo
 
466
        if ( atest ) then
 
467
*
 
468
*
 
469
*               PV-reduction
 
470
*       g-terms
 
471
        rg(1)=1/2d0*( f1*cexi(6)+f2*cexi(10)+f3*cexi(11)+f4*cexi(12) )
 
472
        rg(2)=1/2d0*( f1*cexi(10)+f2*cexi(7)+f3*cexi(13)+f4*cexi(14) )
 
473
        rg(3)=1/2d0*( f1*cexi(11)+f2*cexi(13)+f3*cexi(8)+f4*cexi(15) )
 
474
        rg(4)=1/2d0*( f1*cexi(12)+f2*cexi(14)+f3*cexi(15)+f4*cexi(9) )
 
475
*
 
476
        cexj(36)=xpi(1)*cexj(2)-1/2d0*cd0i(1) -rg(1)
 
477
        cexj(37)=xpi(1)*cexj(3)+1/2d0*cd1ij(1,1)-rg(2)
 
478
        cexj(38)=xpi(1)*cexj(4)+1/2d0*cd1ij(2,1)-rg(3)
 
479
        cexj(39)=xpi(1)*cexj(5)+1/2d0*cd1ij(3,1)-rg(4)
 
480
*
 
481
*       terms ~pipi
 
482
*       1)
 
483
        R(31)=1/2d0*(f1*cexj(6)+cd2ijk(1,1,2)-cd0i(1)) -2*cexj(36)
 
484
        R(32)=1/2d0*(f2*cexj(6)+cd2ijk(1,1,3)-cd2ijk(1,1,2))
 
485
        R(33)=1/2d0*(f3*cexj(6)+cd2ijk(1,1,4)-cd2ijk(1,1,3))
 
486
        R(34)=1/2d0*(f4*cexj(6)+cd2ijk(1,1,5)-cd2ijk(1,1,4))
 
487
*       2)
 
488
        R(35)=1/2d0*(f1*cexj(7)+cd2ijk(1,1,2)-cd2ijk(1,1,1))
 
489
        R(36)=1/2d0*(f2*cexj(7)+cd2ijk(2,2,3)-cd2ijk(1,1,2)) -2*cexj(37)
 
490
        R(37)=1/2d0*(f3*cexj(7)+cd2ijk(2,2,4)-cd2ijk(2,2,3))
 
491
        R(38)=1/2d0*(f4*cexj(7)+cd2ijk(2,2,5)-cd2ijk(2,2,4))
 
492
*       3)
 
493
        R(39)=1/2d0*(f1*cexj(8)+cd2ijk(2,2,2)-cd2ijk(2,2,1))
 
494
        R(40)=1/2d0*(f2*cexj(8)+cd2ijk(2,2,3)-cd2ijk(2,2,2))
 
495
        R(41)=1/2d0*(f3*cexj(8)+cd2ijk(3,3,4)-cd2ijk(2,2,3)) -2*cexj(38)
 
496
        R(42)=1/2d0*(f4*cexj(8)+cd2ijk(3,3,5)-cd2ijk(3,3,4))
 
497
*       4)
 
498
        R(43)=1/2d0*(f1*cexj(9)+cd2ijk(3,3,2)-cd2ijk(3,3,1))
 
499
        R(44)=1/2d0*(f2*cexj(9)+cd2ijk(3,3,3)-cd2ijk(3,3,2))
 
500
        R(45)=1/2d0*(f3*cexj(9)+cd2ijk(3,3,4)-cd2ijk(3,3,3))
 
501
        R(46)=1/2d0*(f4*cexj(9)         -cd2ijk(3,3,4) ) -2*cexj(39)
 
502
*
 
503
*       terms ~p1pi
 
504
*       1)
 
505
        R(47)=1/2d0*(f1*cexj(10)+cd2ijk(1,1,2)+cd1ij(1,1)) -cexj(37)
 
506
        R(48)=1/2d0*(f2*cexj(10)+cd2ijk(1,2,3)-cd2ijk(1,1,2)) -cexj(36)
 
507
        R(49)=1/2d0*(f3*cexj(10)+cd2ijk(1,2,4)-cd2ijk(1,2,3))
 
508
        R(50)=1/2d0*(f4*cexj(10)+cd2ijk(1,2,5)-cd2ijk(1,2,4))
 
509
*       2)
 
510
        R(51)=1/2d0*(f1*cexj(11)+cd2ijk(1,2,2)+cd1ij(2,1) ) -cexj(38)
 
511
        R(52)=1/2d0*(f2*cexj(11)+cd2ijk(1,2,3)-cd2ijk(1,2,2))
 
512
        R(53)=1/2d0*(f3*cexj(11)+cd2ijk(1,3,4)-cd2ijk(1,2,3)) -cexj(36)
 
513
        R(54)=1/2d0*(f4*cexj(11)+cd2ijk(1,3,5)-cd2ijk(1,3,4))
 
514
*       3)
 
515
        R(55)=1/2d0*(f1*cexj(12)+cd2ijk(1,3,2)+cd1ij(3,1) ) -cexj(39)
 
516
        R(56)=1/2d0*(f2*cexj(12)+cd2ijk(1,3,3)-cd2ijk(1,3,2))
 
517
        R(57)=1/2d0*(f3*cexj(12)+cd2ijk(1,3,4)-cd2ijk(1,3,3))
 
518
        R(58)=1/2d0*(f4*cexj(12)         -cd2ijk(1,3,4) ) -cexj(36)
 
519
*
 
520
*       terms ~p2pi
 
521
*       1)
 
522
        R(59)=1/2d0*(f1*cexj(13)+cd2ijk(1,2,2)-cd2ijk(1,2,1))
 
523
        R(60)=1/2d0*(f2*cexj(13)+cd2ijk(2,2,3)-cd2ijk(1,2,2)) -cexj(38)
 
524
        R(61)=1/2d0*(f3*cexj(13)+cd2ijk(2,3,4)-cd2ijk(2,2,3)) -cexj(37)
 
525
        R(62)=1/2d0*(f4*cexj(13)+cd2ijk(2,3,5)-cd2ijk(2,3,4))
 
526
*       2)
 
527
        R(63)=1/2d0*(f1*cexj(14)+cd2ijk(1,3,2)-cd2ijk(1,3,1))
 
528
        R(64)=1/2d0*(f2*cexj(14)+cd2ijk(2,3,3)-cd2ijk(1,3,2)) -cexj(39)
 
529
        R(65)=1/2d0*(f3*cexj(14)+cd2ijk(2,3,4)-cd2ijk(2,3,3))
 
530
        R(66)=1/2d0*(f4*cexj(14)         -cd2ijk(2,3,4) ) -cexj(37)
 
531
*
 
532
*       terms ~p3pi
 
533
*       1)
 
534
        R(67)=1/2d0*(f1*cexj(15)+cd2ijk(2,3,2)-cd2ijk(2,3,1))
 
535
        R(68)=1/2d0*(f2*cexj(15)+cd2ijk(2,3,3)-cd2ijk(2,3,2))
 
536
        R(69)=1/2d0*(f3*cexj(15)+cd2ijk(3,3,4)-cd2ijk(2,3,3)) -cexj(39)
 
537
        R(70)=1/2d0*(f4*cexj(15)         -cd2ijk(3,3,4) ) -cexj(38)
 
538
 
 
539
        cexj(16)=xi5(1)*R(31)+xi5(5)*R(32)+xi5(6)*R(33)+xi5(7) *R(34)
 
540
        cexj(17)=xi5(5)*R(35)+xi5(2)*R(36)+xi5(8)*R(37)+xi5(9) *R(38)
 
541
        cexj(18)=xi5(6)*R(39)+xi5(8)*R(40)+xi5(3)*R(41)+xi5(10)*R(42)
 
542
        cexj(19)=xi5(7)*R(43)+xi5(9)*R(44)+xi5(10)*R(45)+xi5(4)*R(46)
 
543
        cexj(20)=xi5(5)*R(31)+xi5(2)*R(32)+xi5(8)*R(33)+xi5(9) *R(34)
 
544
        cexj(21)=xi5(6)*R(31)+xi5(8)*R(32)+xi5(3)*R(33)+xi5(10)*R(34)
 
545
        cexj(22)=xi5(7)*R(31)+xi5(9)*R(32)+xi5(10)*R(33)+xi5(4)*R(34)
 
546
        cexj(23)=xi5(1)*R(35)+xi5(5)*R(36)+xi5(6)*R(37)+xi5(7) *R(38)
 
547
        cexj(24)=xi5(6)*R(35)+xi5(8)*R(36)+xi5(3)*R(37)+xi5(10)*R(38)
 
548
        cexj(25)=xi5(7)*R(35)+xi5(9)*R(36)+xi5(10)*R(37)+xi5(4)*R(38)
 
549
        cexj(26)=xi5(1)*R(39)+xi5(5)*R(40)+xi5(6)*R(41)+xi5(7) *R(42)
 
550
        cexj(27)=xi5(5)*R(39)+xi5(2)*R(40)+xi5(8)*R(41)+xi5(9) *R(42)
 
551
        cexj(28)=xi5(7)*R(39)+xi5(9)*R(40)+xi5(10)*R(41)+xi5(4)*R(42)
 
552
        cexj(29)=xi5(1)*R(43)+xi5(5)*R(44)+xi5(6)*R(45)+xi5(7) *R(46)
 
553
        cexj(30)=xi5(5)*R(43)+xi5(2)*R(44)+xi5(8)*R(45)+xi5(9) *R(46)
 
554
        cexj(31)=xi5(6)*R(43)+xi5(8)*R(44)+xi5(3)*R(45)+xi5(10)*R(46)
 
555
        cexj(32)=xi5(6)*R(47)+xi5(8)*R(48)+xi5(3)*R(49)+xi5(10)*R(50)
 
556
        cexj(33)=xi5(7)*R(47)+xi5(9)*R(48)+xi5(10)*R(49)+xi5(4)*R(50)
 
557
        cexj(34)=xi5(7)*R(51)+xi5(9)*R(52)+xi5(10)*R(53)+xi5(4)*R(54)
 
558
        cexj(35)=xi5(7)*R(59)+xi5(9)*R(60)+xi5(10)*R(61)+xi5(4)*R(62)
 
559
*
 
560
 
 
561
*###[ : redundancy check
 
562
        cxy(20)=xi5(1)*R(47)+xi5(5)*R(48)+xi5(6)*R(49)+xi5(7) *R(50)
 
563
        cxy(21)=xi5(1)*R(51)+xi5(5)*R(52)+xi5(6)*R(53)+xi5(7) *R(54)
 
564
        cxy(22)=xi5(1)*R(55)+xi5(5)*R(56)+xi5(6)*R(57)+xi5(7) *R(58)
 
565
        cxy(23)=xi5(5)*R(47)+xi5(2)*R(48)+xi5(8)*R(49)+xi5(9) *R(50)
 
566
        cxy(24)=xi5(5)*R(59)+xi5(2)*R(60)+xi5(8)*R(61)+xi5(9) *R(62)
 
567
        cxy(25)=xi5(5)*R(63)+xi5(2)*R(64)+xi5(8)*R(65)+xi5(9) *R(66)
 
568
        cxy(26)=xi5(6)*R(51)+xi5(8)*R(52)+xi5(3)*R(53)+xi5(10)*R(54)
 
569
        cxy(27)=xi5(6)*R(59)+xi5(8)*R(60)+xi5(3)*R(61)+xi5(10)*R(62)
 
570
        cxy(28)=xi5(6)*R(67)+xi5(8)*R(68)+xi5(3)*R(69)+xi5(10)*R(70)
 
571
        cxy(29)=xi5(7)*R(55)+xi5(9)*R(56)+xi5(10)*R(57)+xi5(4)*R(58)
 
572
        cxy(30)=xi5(7)*R(63)+xi5(9)*R(64)+xi5(10)*R(65)+xi5(4)*R(66)
 
573
        cxy(31)=xi5(7)*R(67)+xi5(9)*R(68)+xi5(10)*R(69)+xi5(4)*R(70)
 
574
        cxy(32)=xi5(5)*R(51)+xi5(2)*R(52)+xi5(8)*R(53)+xi5(9) *R(54)
 
575
        cxy(32)=xi5(1)*R(59)+xi5(5)*R(60)+xi5(6)*R(61)+xi5(7) *R(62)
 
576
        cxy(33)=xi5(5)*R(55)+xi5(2)*R(56)+xi5(8)*R(57)+xi5(9) *R(58)
 
577
        cxy(33)=xi5(1)*R(63)+xi5(5)*R(64)+xi5(6)*R(65)+xi5(7) *R(66)
 
578
        cxy(34)=xi5(6)*R(55)+xi5(8)*R(56)+xi5(3)*R(57)+xi5(10)*R(58)
 
579
        cxy(34)=xi5(1)*R(67)+xi5(5)*R(68)+xi5(6)*R(69)+xi5(7) *R(70)
 
580
        cxy(35)=xi5(6)*R(63)+xi5(8)*R(64)+xi5(3)*R(65)+xi5(10)*R(66)
 
581
        cxy(35)=xi5(5)*R(67)+xi5(2)*R(68)+xi5(8)*R(69)+xi5(9) *R(70)
 
582
        do 16 i=20,35
 
583
           if ( absc(cxy(i)-cexj(i)) .gt. .1d-3  ) then
 
584
              print *,'ffxdxp: redundancy check at level 3 failed ',
 
585
     +          i,cxy(i),cexj(i),cxy(i)-cexj(i),ier
 
586
              return
 
587
           endif
 
588
 16     continue
 
589
*       
 
590
*       check against FF
 
591
*       
 
592
        do 17 i=16,35
 
593
            if ( xloss*10d0**(-mod(ier,50))*absc(cexi(i)-cexj(i)) .gt. 
 
594
     +                  precc*absc(cexi(i)) ) then
 
595
                print *,'ffxexp: error: FF disagrees with PV: ',i,
 
596
     +                  cexi(i),cexj(i),cexi(i)-cexj(i),ier
 
597
            endif
 
598
   17   continue
 
599
*###] :
 
600
        endif
 
601
        if ( awrite ) then
 
602
*###[ : print level 3
 
603
            print *,'     '
 
604
            print *,'aaxex : level 3 '
 
605
            print *,'E31 =',cexi(16)
 
606
            print *,'E32 =',cexi(17)
 
607
            print *,'E33 =',cexi(18)
 
608
            print *,'E34 =',cexi(19)
 
609
            print *,'E35 =',cexi(20)
 
610
            print *,'E36 =',cexi(21)
 
611
            print *,'E37 =',cexi(22)
 
612
            print *,'E38 =',cexi(23)
 
613
            print *,'E39 =',cexi(24)
 
614
            print *,'E310=',cexi(25)
 
615
            print *,'E311=',cexi(26)
 
616
            print *,'E312=',cexi(27)
 
617
            print *,'E313=',cexi(28)
 
618
            print *,'E314=',cexi(29)
 
619
            print *,'E315=',cexi(30)
 
620
            print *,'E316=',cexi(31)
 
621
            print *,'E317=',cexi(32)
 
622
            print *,'E318=',cexi(33)
 
623
            print *,'E319=',cexi(34)
 
624
            print *,'E320=',cexi(35)
 
625
*###] :
 
626
        endif
 
627
*
 
628
        if (level .eq. 3) return
 
629
*
 
630
*###] : level 3 :
 
631
        if (level .gt. 3) then
 
632
            print *,'higher than third rank not yet implemented'
 
633
            stop
 
634
        endif
 
635
 
 
636
*###] ffxex:
 
637
        end
 
638
*###[ ffeji:
 
639
        subroutine ffeji(cdxi,mdxi,ccxi,mcxi,cbxi,mbxi,caxi,maxi,
 
640
     +          cdxj,mdxj,ccxj,mcxj,cbxj,mbxj,caxj,maxj,level)
 
641
***#[*comment:***********************************************************
 
642
*                                                                       *
 
643
*       Copy the fourpoint arrays [cm][dcba]xj to the five point arrays *
 
644
*       [cm][dcba]xi.                                                   *
 
645
*                                                                       *
 
646
***#]*comment:***********************************************************
 
647
*  #[ declarations:
 
648
*
 
649
        integer level
 
650
        DOUBLE PRECISION mdxi(55),mcxi(30),mbxi(30),maxi(5),
 
651
     +          mdxj(120),mcxj(140),mbxj(60),maxj(20)
 
652
        DOUBLE COMPLEX cdxi(55),ccxi(30),cbxi(30),caxi(5),
 
653
     +          cdxj(120),ccxj(140),cbxj(60),caxj(20)
 
654
*
 
655
        integer i,j,k
 
656
*
 
657
        include 'aa.h'
 
658
*
 
659
*  #] declarations:
 
660
*  #[ copy necessary parts to E0 arrays:
 
661
*
 
662
*       1)D-output: reduce the array cdxj(5*24) to cdxi(5*11)
 
663
*                   d's are calculated only to (level-1)
 
664
        do j=1,5
 
665
            do i=1,11
 
666
                cdxi(i+(j-1)*11) = cdxj(i+(j-1)*24)
 
667
            enddo
 
668
        enddo
 
669
*
 
670
*       2)C-output: reduce the array ccxj(20*7) to ccxi(10*3)
 
671
*                   c's are calculated only to (level-2)
 
672
*
 
673
        do j=1,4
 
674
            do i=1,3
 
675
                ccxi(0+i+(j-1)*3) = ccxj(0+i+(j-1)*7)
 
676
            enddo
 
677
        enddo
 
678
        do j=1,3
 
679
            do i=1,3
 
680
                ccxi(12+i+(j-1)*3) = ccxj(35+i+(j-1)*7)
 
681
            enddo
 
682
        enddo
 
683
        do j=1,2
 
684
            do i=1,3
 
685
                ccxi(21+i+(j-1)*3) = ccxj(70+i+(j-1)*7)
 
686
            enddo
 
687
        enddo
 
688
        ccxi(28) = ccxj(106)
 
689
        ccxi(29) = ccxj(107)
 
690
        ccxi(30) = ccxj(108)
 
691
*
 
692
*       check the symmetry in B0(i,j)
 
693
*       if ( atest ) then
 
694
*           do 13 i=1,4
 
695
*               j=4+i
 
696
*               k=8+i
 
697
*               if (    ( cbxj(i)     - cbxj(i+1*12) ) .ne. 0. .or.
 
698
*     +                 ( cbxj(j)     - cbxj(i+2*12) ) .ne. 0. .or.
 
699
*     +                 ( cbxj(k)     - cbxj(i+3*12) ) .ne. 0. .or.
 
700
*     +                 ( cbxj(j+1*12)- cbxj(j+2*12) ) .ne. 0. .or.
 
701
*     +                 ( cbxj(k+1*12)- cbxj(j+3*12) ) .ne. 0. .or.
 
702
*     +                 ( cbxj(k+2*12)- cbxj(k+3*12) ) .ne. 0. ) then
 
703
*                   print *,'error in B0-calculations in aaxcx.for'
 
704
*               endif
 
705
* 13        continue
 
706
*       endif
 
707
*
 
708
*  #] copy necessary parts to E0 arrays:
 
709
*###] ffeji:
 
710
        end