← Back to branch summary

~alifson/chiralityflow/trunk

« back to all changes in this revision

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

  • Committer: andrew.lifson at lu
  • Date: 2021-09-01 15:34:39 UTC
  • Revision ID: andrew.lifson@thep.lu.se-20210901153439-7fasjhav4cp4m88r
testing a new repository of a madgraph folder

Show diffs side-by-side

added added

removed removed

Lines of Context:
vendor/CutTools/src/qcdloop/ffdel3.f
 
1
*###[ ffdel3:
 
2
        subroutine ffdel3(del3,xpi,piDpj,ns,ier)
 
3
***#[*comment:***********************************************************
 
4
*                                                                       *
 
5
*       Calculate del3(piDpj) = det(si.sj)      with                    *
 
6
*       the momenta as follows:                                         *
 
7
*       p(1-3) = s(i)                                                   *
 
8
*       p(4-6) = p(i)                                                   *
 
9
*                                                                       *
 
10
*       Input:  xpi(ns)         (real)  m^2(i),i=1,3; p^2(i-3),i=4,10   *
 
11
*               piDpj(ns,ns)    (real)                                  *
 
12
*               ns              (integer)                               *
 
13
*               ier             (integer)                               *
 
14
*                                                                       *
 
15
*       Output: del3            (real)  det(si.sj)                      *
 
16
*                                                                       *
 
17
***#]*comment:***********************************************************
 
18
*  #[ declarations:
 
19
        implicit none
 
20
*
 
21
*       arguments:
 
22
*
 
23
        integer ns,ier
 
24
        DOUBLE PRECISION del3,xpi(6),piDpj(6,6)
 
25
*
 
26
*       local variables:
 
27
*
 
28
        integer mem,nperm
 
29
        parameter(mem=10,nperm=16)
 
30
        integer i,jj(6),iperm(3,nperm),imem,memarr(mem,3),memind,inow
 
31
        DOUBLE PRECISION s(6),xmax,del3p,xmaxp,rloss
 
32
        save iperm,memind,memarr,inow
 
33
*
 
34
*       common blocks:
 
35
*
 
36
        include 'ff.h'
 
37
*  #] declarations:
 
38
*  #[ data:
 
39
        data memind /0/
 
40
        data memarr /mem*0,mem*0,mem*1/
 
41
        data inow /1/
 
42
*
 
43
*       these are all permutations that give a non-zero result with the
 
44
*       correct sign.  This list was generated with getperm3.
 
45
*
 
46
        data iperm/
 
47
     +          1,2,3,  1,2,5,  1,6,2,  1,4,3,
 
48
     +          1,3,5,  1,4,5,  1,6,4,  1,5,6,
 
49
     +          2,4,3,  2,3,6,  2,4,5,  2,6,4,
 
50
     +          2,5,6,  3,4,5,  3,6,4,  3,5,6/
 
51
*  #] data:
 
52
*  #[ starting point in memory?:
 
53
*
 
54
*       see if we know were to start, if not: go on as last time
 
55
*
 
56
        do 5 i=1,mem
 
57
            if ( id .eq. memarr(i,1) .and. idsub .eq. memarr(i,2) ) then
 
58
                inow = memarr(i,3)
 
59
                goto 6
 
60
            endif
 
61
    5   continue
 
62
    6   continue
 
63
*  #] starting point in memory?:
 
64
*  #[ calculations:
 
65
        imem = inow
 
66
        del3 = 0
 
67
        xmax = 0
 
68
 
 
69
   10   continue
 
70
 
 
71
        jj(1) = iperm(1,inow)
 
72
        jj(3) = iperm(2,inow)
 
73
        jj(5) = iperm(3,inow)
 
74
 
 
75
        jj(2) = iperm(1,inow)
 
76
        jj(4) = iperm(2,inow)
 
77
        jj(6) = iperm(3,inow)
 
78
 
 
79
        s(1) = +piDpj(jj(1),jj(2))*piDpj(jj(3),jj(4))*piDpj(jj(5),jj(6))
 
80
        s(2) = +piDpj(jj(1),jj(4))*piDpj(jj(3),jj(6))*piDpj(jj(5),jj(2))
 
81
        s(3) = +piDpj(jj(1),jj(6))*piDpj(jj(3),jj(2))*piDpj(jj(5),jj(4))
 
82
        s(4) = -piDpj(jj(1),jj(2))*piDpj(jj(3),jj(6))*piDpj(jj(5),jj(4))
 
83
        s(5) = -piDpj(jj(1),jj(6))*piDpj(jj(3),jj(4))*piDpj(jj(5),jj(2))
 
84
        s(6) = -piDpj(jj(1),jj(4))*piDpj(jj(3),jj(2))*piDpj(jj(5),jj(6))
 
85
 
 
86
        del3p = 0
 
87
        xmaxp = 0
 
88
        do 20 i=1,6
 
89
            del3p = del3p + s(i)
 
90
            xmaxp = max(xmaxp,abs(s(i)))
 
91
   20   continue
 
92
        if ( abs(del3p) .lt. xloss*xmaxp ) then
 
93
            if ( lwrite ) print *,'del3+',inow,' = ',del3p,xmaxp
 
94
            if ( inow .eq. imem .or. xmaxp .lt. xmax ) then
 
95
                del3 = del3p
 
96
                xmax = xmaxp
 
97
            endif
 
98
            inow = inow + 1
 
99
            if ( inow .gt. nperm ) inow = 1
 
100
            if ( inow .eq. imem ) then
 
101
                if ( lwarn ) call ffwarn(73,ier,del3,xmax)
 
102
                goto 800
 
103
            endif
 
104
            goto 10
 
105
        endif
 
106
        if ( inow .ne. imem ) then
 
107
            if ( lwrite ) print *,'del3+',inow,' = ',del3p,xmaxp
 
108
        endif
 
109
        del3 = del3p
 
110
        xmax = xmaxp
 
111
*  #] calculations:
 
112
*  #[ into memory:
 
113
  800   continue
 
114
        memind = memind + 1
 
115
        if ( memind .gt. mem ) memind = 1
 
116
        memarr(memind,1) = id
 
117
        memarr(memind,2) = idsub
 
118
        memarr(memind,3) = inow
 
119
*  #] into memory:
 
120
*  #[ check output:
 
121
        if ( ltest ) then
 
122
 
 
123
            s(1) = +piDpj(1,1)*piDpj(2,2)*piDpj(3,3)
 
124
            s(2) = +piDpj(1,2)*piDpj(2,3)*piDpj(3,1)
 
125
            s(3) = +piDpj(1,3)*piDpj(2,1)*piDpj(3,2)
 
126
            s(4) = -piDpj(1,1)*piDpj(2,3)*piDpj(3,2)
 
127
            s(5) = -piDpj(1,3)*piDpj(2,2)*piDpj(3,1)
 
128
            s(6) = -piDpj(1,2)*piDpj(2,1)*piDpj(3,3)
 
129
 
 
130
            del3p = 0
 
131
            xmaxp = 0
 
132
            do 820 i=1,6
 
133
                del3p = del3p + s(i)
 
134
                xmaxp = max(xmaxp,abs(s(i)))
 
135
  820       continue
 
136
            rloss = xloss*DBLE(10)**(-mod(ier,50))
 
137
            if ( rloss*abs(del3p-del3) .gt. precx*xmaxp ) then
 
138
                print *,'ffdel3: error: result does not agree with',
 
139
     +                  ' normal case'
 
140
                print *,'result: ',del3,xmax
 
141
                print *,'normal: ',del3p,xmaxp
 
142
                print *,'diff.:  ',del3-del3p
 
143
            endif
 
144
        endif
 
145
*  #] check output:
 
146
*###] ffdel3:
 
147
        end
 
148
*(##[ ffdl3s:
 
149
        subroutine ffdl3s(dl3s,xpi,piDpj,ii,ns,ier)
 
150
***#[*comment:***********************************************************
 
151
*                                                                       *
 
152
*       Calculate dl3s(piDpj) = det(si.sj)      with                    *
 
153
*       the momenta indicated by the indices ii(1-6,1), ii(1-6,2)       *
 
154
*       as follows:                                                     *
 
155
*       p(|ii(1,)|-|ii(3,)|) = s(i)                                     *
 
156
*       p(|ii(4,)|-|ii(6,)|) = p(i) = sgn(ii())*(s(i+1) - s(i))         *
 
157
*                                                                       *
 
158
*       At this moment (26-apr-1990) only the diagonal is tried         *                                                                       
 
159
*                                                                       *                                                                       
 
160
*       Input:  xpi(ns)         (real)  m^2(i),i=1,3; p^2(i-3),i=4,10   *
 
161
*               piDpj(ns,ns)    (real)                                  *
 
162
*               ii(6,2)         (integer)       see above               *
 
163
*               ns              (integer)                               *
 
164
*               ier             (integer)                               *
 
165
*                                                                       *
 
166
*       Output: dl3s            (real)  det(si.sj)                      *
 
167
*                                                                       *
 
168
***#]*comment:***********************************************************
 
169
*  #[ declarations:
 
170
        implicit none
 
171
*
 
172
*       arguments:
 
173
*
 
174
        integer ii(6,2),ns,ier
 
175
        DOUBLE PRECISION dl3s,xpi(ns),piDpj(ns,ns)
 
176
*
 
177
*       local variables:
 
178
*
 
179
        integer mem,nperm
 
180
        parameter(mem=10,nperm=16)
 
181
        integer i,j,jj(6),jsgn,iperm(3,nperm),imem,memarr(mem,3),
 
182
     +          memind,inow
 
183
        DOUBLE PRECISION s(6),xmax,dl3sp,xmaxp,xlosn,xhck,rloss
 
184
        save iperm,memind,memarr,inow
 
185
*
 
186
*       common blocks:
 
187
*
 
188
        include 'ff.h'
 
189
*  #] declarations:
 
190
*  #[ data:
 
191
        data memind /0/
 
192
        data memarr /mem*0,mem*0,mem*1/
 
193
        data inow /1/
 
194
*
 
195
*       these are all permutations that give a non-zero result with the
 
196
*       correct sign.  This list was generated with getperm3.
 
197
*
 
198
        data iperm/
 
199
     +          1,2,3,  1,2,5,  1,6,2,  1,4,3,
 
200
     +          1,3,5,  1,4,5,  1,6,4,  1,5,6,
 
201
     +          2,4,3,  2,3,6,  2,4,5,  2,6,4,
 
202
     +          2,5,6,  3,4,5,  3,6,4,  3,5,6/
 
203
*  #] data:
 
204
*  #[ test input:
 
205
        if ( ltest ) then
 
206
*           print *,'ffdl3s: input: ii(,1) = ',(ii(i,1),i=1,6)
 
207
*           print *,'               ii(,2) = ',(ii(i,2),i=1,6)
 
208
            xlosn = xloss*DBLE(10)**(-mod(ier,50))
 
209
            do 3 j=1,2
 
210
            do 1 i=1,6
 
211
                if ( abs(ii(i,j)) .gt. ns ) print *,'ffdl3s: error: ',
 
212
     +                  '|ii(i,j)| > ns: ',ii(i,j),ns
 
213
                if ( abs(ii(i,j)) .eq. 0 ) print *,'ffdl3s: error: ',
 
214
     +                  '|ii(i,j)| = 0: ',ii(i,j)
 
215
    1       continue
 
216
            do 2 i=1,6
 
217
 
 
218
                xhck = piDpj(abs(ii(i,j)),ii(1,j))
 
219
     +               - piDpj(abs(ii(i,j)),ii(2,j))
 
220
     +               + sign(1,ii(4,j))*piDpj(abs(ii(i,j)),abs(ii(4,j)))
 
221
                xmax = max(abs(piDpj(abs(ii(i,j)),ii(1,j))),
 
222
     +                     abs(piDpj(abs(ii(i,j)),ii(2,j))))
 
223
                if ( xlosn*abs(xhck) .gt. precx*xmax ) print *,'ffdl3s:'
 
224
     +            ,' error: dotproducts 124 with ',i,' do not add to 0:'
 
225
     +            ,piDpj(abs(ii(i,j)),ii(1,j)),
 
226
     +            piDpj(abs(ii(i,j)),ii(2,j)),
 
227
     +            piDpj(abs(ii(i,j)),abs(ii(4,j))),xhck
 
228
 
 
229
                xhck = piDpj(abs(ii(i,j)),ii(2,j))
 
230
     +               - piDpj(abs(ii(i,j)),ii(3,j))
 
231
     +               + sign(1,ii(5,j))*piDpj(abs(ii(i,j)),abs(ii(5,j))) 
 
232
                xmax = max(abs(piDpj(abs(ii(i,j)),ii(2,j))),
 
233
     +                  abs(piDpj(abs(ii(i,j)),ii(3,j))))
 
234
                if ( xlosn*abs(xhck) .gt. precx*xmax ) print *,'ffdl3s:'
 
235
     +            ,' error: dotproducts 235 with ',i,' do not add to 0:'
 
236
     +            ,piDpj(abs(ii(i,j)),ii(2,j)),
 
237
     +            piDpj(abs(ii(i,j)),ii(3,j)),
 
238
     +            piDpj(abs(ii(i,j)),abs(ii(5,j))),xhck
 
239
 
 
240
                xhck = piDpj(abs(ii(i,j)),ii(3,j)) 
 
241
     +               - piDpj(abs(ii(i,j)),ii(1,j))
 
242
     +               + sign(1,ii(6,j))*piDpj(abs(ii(i,j)),abs(ii(6,j))) 
 
243
                xmax = max(abs(piDpj(abs(ii(i,j)),ii(3,j))),
 
244
     +                  abs(piDpj(abs(ii(i,j)),ii(1,j))))
 
245
                if ( xlosn*abs(xhck) .gt. precx*xmax ) print *,'ffdl3s:'
 
246
     +            ,' error: dotproducts 316 with ',i,' do not add to 0:'
 
247
     +            ,piDpj(abs(ii(i,j)),ii(3,j)),
 
248
     +            piDpj(abs(ii(i,j)),ii(1,j)),
 
249
     +            piDpj(abs(ii(i,j)),abs(ii(6,j))),xhck
 
250
 
 
251
                xhck = sign(1,ii(4,j))*piDpj(abs(ii(i,j)),abs(ii(4,j))) 
 
252
     +               + sign(1,ii(5,j))*piDpj(abs(ii(i,j)),abs(ii(5,j)))
 
253
     +               + sign(1,ii(6,j))*piDpj(abs(ii(i,j)),abs(ii(6,j)))
 
254
                xmax = max(abs(piDpj(abs(ii(i,j)),abs(ii(4,j)))),
 
255
     +                     abs(piDpj(abs(ii(i,j)),abs(ii(5,j)))))
 
256
                if ( xlosn*abs(xhck) .gt. precx*xmax ) print *,'ffdl3s:'
 
257
     +            ,' error: dotproducts 456 with ',i,' do not add to 0:'
 
258
     +            ,piDpj(abs(ii(i,j)),abs(ii(4,j))),
 
259
     +            piDpj(abs(ii(i,j)),abs(ii(5,j))),
 
260
     +            piDpj(abs(ii(i,j)),abs(ii(6,j))),xhck
 
261
 
 
262
    2       continue
 
263
    3       continue
 
264
        endif
 
265
*  #] test input:
 
266
*  #[ starting point in memory?:
 
267
*
 
268
*       see if we know were to start, if not: go on as last time
 
269
*
 
270
        do 5 i=1,mem
 
271
            if ( id .eq. memarr(i,1) .and. idsub .eq. memarr(i,2) ) then
 
272
                inow = memarr(i,3)
 
273
                goto 6
 
274
            endif
 
275
    5   continue
 
276
    6   continue
 
277
*  #] starting point in memory?:
 
278
*  #[ calculations:
 
279
        imem = inow
 
280
        dl3s = 0
 
281
        xmax = 0
 
282
 
 
283
   10   continue
 
284
 
 
285
        jj(1) = abs(ii(iperm(1,inow),1))
 
286
        jj(3) = abs(ii(iperm(2,inow),1))
 
287
        jj(5) = abs(ii(iperm(3,inow),1))
 
288
 
 
289
        jj(2) = abs(ii(iperm(1,inow),2))
 
290
        jj(4) = abs(ii(iperm(2,inow),2))
 
291
        jj(6) = abs(ii(iperm(3,inow),2))
 
292
 
 
293
        jsgn =  sign(1,ii(iperm(1,inow),1))
 
294
     +          *sign(1,ii(iperm(2,inow),1))
 
295
     +          *sign(1,ii(iperm(3,inow),1))
 
296
     +          *sign(1,ii(iperm(1,inow),2))
 
297
     +          *sign(1,ii(iperm(2,inow),2))
 
298
     +          *sign(1,ii(iperm(3,inow),2))
 
299
 
 
300
        s(1) = +piDpj(jj(1),jj(2))*piDpj(jj(3),jj(4))*piDpj(jj(5),jj(6))
 
301
        s(2) = +piDpj(jj(1),jj(4))*piDpj(jj(3),jj(6))*piDpj(jj(5),jj(2))
 
302
        s(3) = +piDpj(jj(1),jj(6))*piDpj(jj(3),jj(2))*piDpj(jj(5),jj(4))
 
303
        s(4) = -piDpj(jj(1),jj(2))*piDpj(jj(3),jj(6))*piDpj(jj(5),jj(4))
 
304
        s(5) = -piDpj(jj(1),jj(6))*piDpj(jj(3),jj(4))*piDpj(jj(5),jj(2))
 
305
        s(6) = -piDpj(jj(1),jj(4))*piDpj(jj(3),jj(2))*piDpj(jj(5),jj(6))
 
306
 
 
307
        dl3sp = 0
 
308
        xmaxp = 0
 
309
        do 20 i=1,6
 
310
            dl3sp = dl3sp + s(i)
 
311
            xmaxp = max(xmaxp,abs(s(i)))
 
312
   20   continue
 
313
        if ( abs(dl3sp) .lt. xloss*xmaxp ) then
 
314
            if ( lwrite ) print *,'dl3s+',inow,' = ',dl3sp,xmaxp
 
315
            if ( inow .eq. imem .or. xmaxp .lt. xmax ) then
 
316
                dl3s = jsgn*dl3sp
 
317
                xmax = xmaxp
 
318
            endif
 
319
            inow = inow + 1
 
320
            if ( inow .gt. nperm ) inow = 1
 
321
            if ( inow .eq. imem ) then
 
322
                if ( lwarn ) call ffwarn(85,ier,dl3s,xmax)
 
323
                goto 800
 
324
            endif
 
325
            goto 10
 
326
        endif
 
327
        if ( inow .ne. imem ) then
 
328
            if ( lwrite ) print *,'dl3s+',inow,' = ',dl3sp,xmaxp
 
329
        endif
 
330
        dl3s = jsgn*dl3sp
 
331
        xmax = xmaxp
 
332
*  #] calculations:
 
333
*  #[ into memory:
 
334
  800   continue
 
335
        memind = memind + 1
 
336
        if ( memind .gt. mem ) memind = 1
 
337
        memarr(memind,1) = id
 
338
        memarr(memind,2) = idsub
 
339
        memarr(memind,3) = inow
 
340
*  #] into memory:
 
341
*  #[ check output:
 
342
        if ( ltest ) then
 
343
 
 
344
            s(1) = +piDpj(ii(1,1),ii(1,2))*piDpj(ii(2,1),ii(2,2))*
 
345
     +                                          piDpj(ii(3,1),ii(3,2))
 
346
            s(2) = +piDpj(ii(1,1),ii(2,2))*piDpj(ii(2,1),ii(3,2))*
 
347
     +                                          piDpj(ii(3,1),ii(1,2))
 
348
            s(3) = +piDpj(ii(1,1),ii(3,2))*piDpj(ii(3,1),ii(2,2))*
 
349
     +                                          piDpj(ii(2,1),ii(1,2))
 
350
            s(4) = -piDpj(ii(1,1),ii(1,2))*piDpj(ii(2,1),ii(3,2))*
 
351
     +                                          piDpj(ii(3,1),ii(2,2))
 
352
            s(5) = -piDpj(ii(1,1),ii(3,2))*piDpj(ii(2,1),ii(2,2))*
 
353
     +                                          piDpj(ii(3,1),ii(1,2))
 
354
            s(6) = -piDpj(ii(1,1),ii(2,2))*piDpj(ii(2,1),ii(1,2))*
 
355
     +                                          piDpj(ii(3,1),ii(3,2))
 
356
 
 
357
            dl3sp = 0
 
358
            xmaxp = 0
 
359
            do 820 i=1,6
 
360
                dl3sp = dl3sp + s(i)
 
361
                xmaxp = max(xmaxp,abs(s(i)))
 
362
  820       continue
 
363
            rloss = xloss*DBLE(10)**(-mod(ier,50))
 
364
            if ( rloss*abs(dl3sp-dl3s) .gt. precx*xmaxp ) then
 
365
                print *,'ffdl3s: error: result does not agree with',
 
366
     +                  ' normal case'
 
367
                print *,'result: ',dl3s,xmax
 
368
                print *,'normal: ',dl3sp,xmaxp
 
369
                print *,'diff.:  ',dl3s-dl3sp
 
370
            endif
 
371
        endif
 
372
*  #] check output:
 
373
*)##] ffdl3s:
 
374
        end