← Back to branch summary

~maddevelopers/mg5amcnlo/2.9.4

« back to all changes in this revision

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

pass to v2.0.0

Show diffs side-by-side

added added

removed removed

Lines of Context:
vendor/CutTools/src/qcdloop/ffcel4.f
 
1
*###[ ffcel4:
 
2
        subroutine ffcel4(del4,xpi,piDpj,ns,ier)
 
3
***#[*comment:***********************************************************
 
4
*                                                                       *
 
5
*       Calculate del4(piDpj) = det(si.sj)      with                    *
 
6
*       the momenta as follows:                                         *
 
7
*       p(1-4) = s(i)                                                   *
 
8
*       p(4-10) = 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: del4            (real)  det(si.sj)                      *
 
16
*                                                                       *
 
17
***#]*comment:*********************************************************** 
 
18
*  #[ declarations:
 
19
        implicit none
 
20
*
 
21
*       arguments:
 
22
*
 
23
        integer ns,ier
 
24
        DOUBLE COMPLEX del4,xpi(10),piDpj(10,10)
 
25
*
 
26
*       local variables:
 
27
*
 
28
        integer mem,nperm
 
29
        parameter(mem=10,nperm=125)
 
30
        integer i,jj(8),iperm(4,nperm),imem,jmem,memarr(mem,4),memind,
 
31
     +          inow,jnow,icount
 
32
        DOUBLE PRECISION xmax,xmaxp,absc,rloss
 
33
        DOUBLE COMPLEX s(24),del4p,c
 
34
        save iperm,memind,memarr,inow,jnow
 
35
*
 
36
*       common blocks:
 
37
*
 
38
        include 'ff.h'
 
39
*
 
40
*       statement functions:
 
41
*
 
42
        absc(c) = abs(DBLE(c)) + abs(DIMAG(c))
 
43
*  #] declarations: 
 
44
*  #[ data:
 
45
        data memind /0/
 
46
        data memarr /mem*0,mem*0,mem*1,mem*1/
 
47
        data inow /1/
 
48
        data jnow /1/
 
49
*
 
50
*       these are all permutations that give a non-zero result with the
 
51
*       correct sign.  This list was generated with getperm4.
 
52
*       (note: this used to be well-ordened, but then it had more than
 
53
*       19 continuation lines)
 
54
*
 
55
        data iperm/
 
56
     +  1,2,3,4,1,2,3,7,1,2,8,3,1,2,3,10,1,2,6,4,1,2,4,7,1,2,4,9,1,2,6,7
 
57
     +  ,1,2,8,6,1,2,6,10,1,2,7,8,1,2,7,9,1,2,10,7,1,2,9,8,1,2,10,9,1,3,
 
58
     +  4,5,1,3,6,4,1,3,10,4,1,3,7,5,1,3,5,8,1,3,10,5,1,3,6,7,1,3,8,6,1,
 
59
     +  3,6,10,1,3,10,7,1,3,8,10,1,4,5,6,1,4,7,5,1,4,9,5,1,4,6,7,1,4,6,9
 
60
     +  ,1,4,6,10,1,4,10,7,1,4,10,9,1,5,6,7,1,5,8,6,1,5,6,10,1,5,7,8,1,5
 
61
     +  ,7,9,1,5,10,7,1,5,9,8,1,5,10,9,1,6,8,7,1,6,9,7,1,6,8,9,1,6,8,10,
 
62
     +  1,6,9,10,1,7,10,8,1,7,10,9,1,8,9,10,2,3,4,5,2,3,8,4,2,3,9,4,2,3,
 
63
     +  7,5,2,3,5,8,2,3,10,5,2,3,8,7,2,3,9,7,2,3,8,9,2,3,8,10,2,3,9,10,2
 
64
     +  ,4,5,6,2,4,7,5,2,4,9,5,2,4,6,8,2,4,6,9,2,4,8,7,2,4,9,7,2,4,8,9,2
 
65
     +  ,5,6,7,2,5,8,6,2,5,6,10,2,5,7,8,2,5,7,9,2,5,10,7,2,5,9,8,2,5,10,
 
66
     +  9,2,6,8,7,2,6,9,7,2,6,8,9,2,6,8,10,2,6,9,10,2,7,10,8,2,7,10,9,2,
 
67
     +  8,9,10,3,4,5,6,3,4,8,5,3,4,9,5,3,4,5,10,3,4,6,8,3,4,6,9,3,4,10,8
 
68
     +  ,3,4,10,9,3,5,6,7,3,5,8,6,3,5,6,10,3,5,7,8,3,5,7,9,3,5,10,7,3,5,
 
69
     +  9,8,3,5,10,9,3,6,8,7,3,6,9,7,3,6,8,9,3,6,8,10,3,6,9,10,3,7,10,8,
 
70
     +  3,7,10,9,3,8,9,10,4,5,6,7,4,5,8,6,4,5,6,10,4,5,7,8,4,5,7,9,4,5,1
 
71
     +  0,7,4,5,9,8,4,5,10,9,4,6,8,7,4,6,9,7,4,6,8,9,4,6,8,10,4,6,9,10,4
 
72
     +  ,7,10,8,4,7,10,9,4,8,9,10/
 
73
*  #] data: 
 
74
*  #[ get starting point from memory:
 
75
*
 
76
*       see if we know were to start, if not: go on as last time
 
77
*
 
78
        do 5 i=1,mem
 
79
            if ( id .eq. memarr(i,1) .and. idsub .eq. memarr(i,2) ) then
 
80
                inow = memarr(i,3)
 
81
                jnow = memarr(i,4)
 
82
                if ( lwrite ) print *,'ffcel4: from memory: ',id,idsub,
 
83
     +                  inow,jnow
 
84
                goto 6
 
85
            endif
 
86
    5   continue
 
87
    6   continue
 
88
*  #] get starting point from memory: 
 
89
*  #[ calculations:
 
90
        imem = inow
 
91
        jmem = jnow
 
92
        del4 = 0
 
93
        xmax = 0
 
94
        icount = 0
 
95
*
 
96
   10   continue
 
97
 
 
98
        jj(1) = iperm(1,inow)
 
99
        jj(3) = iperm(2,inow)
 
100
        jj(5) = iperm(3,inow)
 
101
        jj(7) = iperm(4,inow)
 
102
 
 
103
        jj(2) = iperm(1,jnow)
 
104
        jj(4) = iperm(2,jnow)
 
105
        jj(6) = iperm(3,jnow)
 
106
        jj(8) = iperm(4,jnow)
 
107
 
 
108
        s( 1) = +piDpj(jj(1),jj(2))*piDpj(jj(3),jj(4))*
 
109
     +          piDpj(jj(5),jj(6))*piDpj(jj(7),jj(8))
 
110
        s( 2) = +piDpj(jj(1),jj(4))*piDpj(jj(3),jj(6))*
 
111
     +          piDpj(jj(5),jj(2))*piDpj(jj(7),jj(8))
 
112
        s( 3) = +piDpj(jj(1),jj(6))*piDpj(jj(3),jj(2))*
 
113
     +          piDpj(jj(5),jj(4))*piDpj(jj(7),jj(8))
 
114
        s( 4) = -piDpj(jj(1),jj(2))*piDpj(jj(3),jj(6))*
 
115
     +          piDpj(jj(5),jj(4))*piDpj(jj(7),jj(8))
 
116
        s( 5) = -piDpj(jj(1),jj(6))*piDpj(jj(3),jj(4))*
 
117
     +          piDpj(jj(5),jj(2))*piDpj(jj(7),jj(8))
 
118
        s( 6) = -piDpj(jj(1),jj(4))*piDpj(jj(3),jj(2))*
 
119
     +          piDpj(jj(5),jj(6))*piDpj(jj(7),jj(8))
 
120
 
 
121
        s( 7) = -piDpj(jj(1),jj(2))*piDpj(jj(3),jj(4))*
 
122
     +          piDpj(jj(7),jj(6))*piDpj(jj(5),jj(8))
 
123
        s( 8) = -piDpj(jj(1),jj(4))*piDpj(jj(3),jj(6))*
 
124
     +          piDpj(jj(7),jj(2))*piDpj(jj(5),jj(8))
 
125
        s( 9) = -piDpj(jj(1),jj(6))*piDpj(jj(3),jj(2))*
 
126
     +          piDpj(jj(7),jj(4))*piDpj(jj(5),jj(8))
 
127
        s(10) = +piDpj(jj(1),jj(2))*piDpj(jj(3),jj(6))*
 
128
     +          piDpj(jj(7),jj(4))*piDpj(jj(5),jj(8))
 
129
        s(11) = +piDpj(jj(1),jj(6))*piDpj(jj(3),jj(4))*
 
130
     +          piDpj(jj(7),jj(2))*piDpj(jj(5),jj(8))
 
131
        s(12) = +piDpj(jj(1),jj(4))*piDpj(jj(3),jj(2))*
 
132
     +          piDpj(jj(7),jj(6))*piDpj(jj(5),jj(8))
 
133
 
 
134
        s(13) = -piDpj(jj(1),jj(2))*piDpj(jj(7),jj(4))*
 
135
     +          piDpj(jj(5),jj(6))*piDpj(jj(3),jj(8))
 
136
        s(14) = -piDpj(jj(1),jj(4))*piDpj(jj(7),jj(6))*
 
137
     +          piDpj(jj(5),jj(2))*piDpj(jj(3),jj(8))
 
138
        s(15) = -piDpj(jj(1),jj(6))*piDpj(jj(7),jj(2))*
 
139
     +          piDpj(jj(5),jj(4))*piDpj(jj(3),jj(8))
 
140
        s(16) = +piDpj(jj(1),jj(2))*piDpj(jj(7),jj(6))*
 
141
     +          piDpj(jj(5),jj(4))*piDpj(jj(3),jj(8))
 
142
        s(17) = +piDpj(jj(1),jj(6))*piDpj(jj(7),jj(4))*
 
143
     +          piDpj(jj(5),jj(2))*piDpj(jj(3),jj(8))
 
144
        s(18) = +piDpj(jj(1),jj(4))*piDpj(jj(7),jj(2))*
 
145
     +          piDpj(jj(5),jj(6))*piDpj(jj(3),jj(8))
 
146
 
 
147
        s(19) = -piDpj(jj(7),jj(2))*piDpj(jj(3),jj(4))*
 
148
     +          piDpj(jj(5),jj(6))*piDpj(jj(1),jj(8))
 
149
        s(20) = -piDpj(jj(7),jj(4))*piDpj(jj(3),jj(6))*
 
150
     +          piDpj(jj(5),jj(2))*piDpj(jj(1),jj(8))
 
151
        s(21) = -piDpj(jj(7),jj(6))*piDpj(jj(3),jj(2))*
 
152
     +          piDpj(jj(5),jj(4))*piDpj(jj(1),jj(8))
 
153
        s(22) = +piDpj(jj(7),jj(2))*piDpj(jj(3),jj(6))*
 
154
     +          piDpj(jj(5),jj(4))*piDpj(jj(1),jj(8))
 
155
        s(23) = +piDpj(jj(7),jj(6))*piDpj(jj(3),jj(4))*
 
156
     +          piDpj(jj(5),jj(2))*piDpj(jj(1),jj(8))
 
157
        s(24) = +piDpj(jj(7),jj(4))*piDpj(jj(3),jj(2))*
 
158
     +          piDpj(jj(5),jj(6))*piDpj(jj(1),jj(8))
 
159
 
 
160
        del4p = 0
 
161
        xmaxp = 0
 
162
        do 20 i=1,24
 
163
            del4p = del4p + s(i)
 
164
            xmaxp = max(xmaxp,absc(s(i)))
 
165
   20   continue
 
166
        if ( absc(del4p) .lt. xloss*xmaxp ) then
 
167
            if ( lwrite ) print *,'del4+',icount,' = ',del4p,xmaxp,inow,
 
168
     +          jnow
 
169
            if ( inow .eq. imem .or. xmaxp .lt. xmax ) then
 
170
                del4 = del4p
 
171
                xmax = xmaxp
 
172
            endif
 
173
*           as the list is ordered we may have more luck stepping
 
174
*           through with large steps
 
175
            inow = inow + 43
 
176
            jnow = jnow + 49
 
177
            if ( inow .gt. nperm ) inow = inow - nperm
 
178
            if ( jnow .gt. nperm ) jnow = jnow - nperm
 
179
            icount = icount + 1
 
180
            if ( icount.gt.15 .or. inow.eq.imem .or. jnow.eq.jmem
 
181
     +                  ) then
 
182
                if ( lwarn ) call ffwarn(143,ier,absc(del4),xmax)
 
183
                goto 800
 
184
            endif
 
185
            goto 10
 
186
        endif
 
187
        if ( inow.ne.imem) then
 
188
            if ( lwrite ) print *,'del4+',icount,' = ',del4p,xmaxp,inow,
 
189
     +          jnow
 
190
        endif
 
191
        del4 = del4p
 
192
        xmax = xmaxp
 
193
*  #] calculations: 
 
194
*  #[ into memory:
 
195
        if ( lwrite ) print *,'ffcel4: into memory: ',id,idsub,inow,jnow
 
196
        memind = memind + 1
 
197
        if ( memind .gt. mem ) memind = 1
 
198
        memarr(memind,1) = id
 
199
        memarr(memind,2) = idsub
 
200
        memarr(memind,3) = inow
 
201
        memarr(memind,4) = jnow
 
202
  800   continue
 
203
*  #] into memory: 
 
204
*  #[ check output:
 
205
        if ( ltest ) then
 
206
 
 
207
            s( 1) = +piDpj(1,1)*piDpj(2,2)*piDpj(3,3)*piDpj(4,4)
 
208
            s( 2) = +piDpj(1,2)*piDpj(2,3)*piDpj(3,1)*piDpj(4,4)
 
209
            s( 3) = +piDpj(1,3)*piDpj(2,1)*piDpj(3,2)*piDpj(4,4)
 
210
            s( 4) = -piDpj(1,1)*piDpj(2,3)*piDpj(3,2)*piDpj(4,4)
 
211
            s( 5) = -piDpj(1,3)*piDpj(2,2)*piDpj(3,1)*piDpj(4,4)
 
212
            s( 6) = -piDpj(1,2)*piDpj(2,1)*piDpj(3,3)*piDpj(4,4)
 
213
 
 
214
            s( 7) = -piDpj(1,1)*piDpj(2,2)*piDpj(4,3)*piDpj(3,4)
 
215
            s( 8) = -piDpj(1,2)*piDpj(2,3)*piDpj(4,1)*piDpj(3,4)
 
216
            s( 9) = -piDpj(1,3)*piDpj(2,1)*piDpj(4,2)*piDpj(3,4)
 
217
            s(10) = +piDpj(1,1)*piDpj(2,3)*piDpj(4,2)*piDpj(3,4)
 
218
            s(11) = +piDpj(1,3)*piDpj(2,2)*piDpj(4,1)*piDpj(3,4)
 
219
            s(12) = +piDpj(1,2)*piDpj(2,1)*piDpj(4,3)*piDpj(3,4)
 
220
 
 
221
            s(13) = -piDpj(1,1)*piDpj(4,2)*piDpj(3,3)*piDpj(2,4)
 
222
            s(14) = -piDpj(1,2)*piDpj(4,3)*piDpj(3,1)*piDpj(2,4)
 
223
            s(15) = -piDpj(1,3)*piDpj(4,1)*piDpj(3,2)*piDpj(2,4)
 
224
            s(16) = +piDpj(1,1)*piDpj(4,3)*piDpj(3,2)*piDpj(2,4)
 
225
            s(17) = +piDpj(1,3)*piDpj(4,2)*piDpj(3,1)*piDpj(2,4)
 
226
            s(18) = +piDpj(1,2)*piDpj(4,1)*piDpj(3,3)*piDpj(2,4)
 
227
 
 
228
            s(19) = -piDpj(4,1)*piDpj(2,2)*piDpj(3,3)*piDpj(1,4)
 
229
            s(20) = -piDpj(4,2)*piDpj(2,3)*piDpj(3,1)*piDpj(1,4)
 
230
            s(21) = -piDpj(4,3)*piDpj(2,1)*piDpj(3,2)*piDpj(1,4)
 
231
            s(22) = +piDpj(4,1)*piDpj(2,3)*piDpj(3,2)*piDpj(1,4)
 
232
            s(23) = +piDpj(4,3)*piDpj(2,2)*piDpj(3,1)*piDpj(1,4)
 
233
            s(24) = +piDpj(4,2)*piDpj(2,1)*piDpj(3,3)*piDpj(1,4)
 
234
 
 
235
            del4p = 0
 
236
            xmaxp = 0
 
237
            do 820 i=1,24
 
238
                del4p = del4p + s(i)
 
239
                xmaxp = max(xmaxp,absc(s(i)))
 
240
  820       continue
 
241
            rloss = xloss*DBLE(10)**(-mod(ier,50)-1)
 
242
            if ( rloss*absc(del4p-del4) .gt. precc*xmaxp ) then
 
243
                print *,'ffcel4: error: result does not agree with',
 
244
     +                  ' normal case'
 
245
                print *,'result: ',del4,xmax
 
246
                print *,'normal: ',del4p,xmaxp
 
247
                print *,'diff.:  ',del4-del4p,ier
 
248
            endif
 
249
        endif
 
250
*  #] check output:
 
251
*###] ffcel4:
 
252
        end
 
253
*###[ ffcl3p:
 
254
        subroutine ffcl3p(dl3p,piDpj,ns,ii,jj,ier)
 
255
***#[*comment:***********************************************************
 
256
*       calculate in a numerically stable way                           *
 
257
*                                                                       *
 
258
*            p1  p2  p3                                                 *
 
259
*       delta                                                           *
 
260
*            p1' p2' p3'                                                *
 
261
*                                                                       *
 
262
*       with pn = xpi(ii(n)), p4 = -p1-p2-p3, p5 = -p1-p2, p6 = p2+p3   *
 
263
*       with pn'= xpi(jj(n)), p4'= etc.       (when ns=15 p5=p1+p2)     *
 
264
*                                                                       *
 
265
*       Input:  piDpj   complex(ns,ns)  dotpruducts                     *
 
266
*               ns      integer         either 10 or 15                 *
 
267
*               ii,jj   integer(6)      location of pi in piDpj         *
 
268
*               ier     integer         number of digits lost so far    *
 
269
*       Output: dl3p    complex         see above                       *
 
270
*               ier     integer         number of digits lost so far    *
 
271
*                                                                       *
 
272
***#]*comment:*********************************************************** 
 
273
*  #[ declarations:
 
274
        implicit none
 
275
*
 
276
*       arguments:
 
277
*
 
278
        integer ns,ii(6),jj(6),ier
 
279
        DOUBLE COMPLEX dl3p,piDpj(ns,ns)
 
280
*
 
281
*       local variables
 
282
*
 
283
        integer i,j,k,l,iperm(3,16),ii1,ii2,ii3,jj1,jj2,jj3,nl
 
284
        DOUBLE PRECISION xmax,smax,absc
 
285
        DOUBLE COMPLEX s(6),som,xheck,c
 
286
*
 
287
*       common blocks
 
288
*
 
289
        include 'ff.h'
 
290
*
 
291
*       statement functions:
 
292
*
 
293
        absc(c) = abs(DBLE(c)) + abs(DIMAG(c))
 
294
*
 
295
*       data
 
296
*
 
297
        data iperm /1,2,3, 2,4,3, 3,4,1, 4,2,1,
 
298
     +              1,2,6, 6,4,3, 3,1,6, 2,4,6,
 
299
     +              2,5,3, 5,4,1, 1,3,5, 2,4,5,
 
300
     +              1,6,5, 2,5,6, 3,6,5, 4,5,6/
 
301
*  #] declarations: 
 
302
*  #[ check input:
 
303
        if ( lwrite ) then
 
304
            print *,'ffcl3p: indices are'
 
305
            print *,ii
 
306
            print *,jj
 
307
        endif
 
308
        if ( ltest ) then
 
309
            if ( ns .ne. 10 .and. ns .ne. 15 ) print *,'ffcl3p: error:',
 
310
     +          ' only tested for ns=10,15'
 
311
            do 10 i=1,ns
 
312
                xheck = +piDpj(i,ii(1))+piDpj(i,ii(2))
 
313
     +                  +piDpj(i,ii(3))+piDpj(i,ii(4))
 
314
                xmax = max(absc(piDpj(i,ii(1))),absc(piDpj(i,ii(2))),
 
315
     +                     absc(piDpj(i,ii(3))),absc(piDpj(i,ii(4))))
 
316
                if ( xloss*absc(xheck) .gt. precc*xmax ) print *,
 
317
     +                  'ffcl3p: error: momenta i1234 do not add to 0:',
 
318
     +                  piDpj(i,ii(1)),piDpj(i,ii(2)),piDpj(i,ii(3)),
 
319
     +                  piDpj(i,ii(4)),xheck,i
 
320
                xheck = piDpj(i,ii(6))-piDpj(i,ii(2))-piDpj(i,ii(3))
 
321
                xmax = max(absc(piDpj(i,ii(6))),absc(piDpj(i,ii(2))),
 
322
     +                     absc(piDpj(i,ii(3))))
 
323
                if ( xloss*absc(xheck) .gt. precc*xmax ) print *,
 
324
     +                  'ffcl3p: error: momenta i623 do not add to 0:',
 
325
     +                  piDpj(i,ii(6)),piDpj(i,ii(2)),piDpj(i,ii(3)),
 
326
     +                  xheck,i
 
327
                if ( ns .eq. 10 ) then
 
328
                    xheck = piDpj(i,ii(5))+piDpj(i,ii(1))+piDpj(i,ii(2))
 
329
                else
 
330
                    xheck = piDpj(i,ii(5))-piDpj(i,ii(1))-piDpj(i,ii(2))
 
331
                endif
 
332
                xmax = max(absc(piDpj(i,ii(5))),absc(piDpj(i,ii(1))),
 
333
     +                     absc(piDpj(i,ii(2))))
 
334
                if ( xloss*absc(xheck) .gt. precc*xmax ) print *,
 
335
     +                  'ffcl3p: error: momenta i512 do not add to 0:',
 
336
     +                  piDpj(i,ii(5)),piDpj(i,ii(1)),piDpj(i,ii(2)),
 
337
     +                  xheck,i
 
338
                xheck = +piDpj(i,jj(1))+piDpj(i,jj(2))
 
339
     +                  +piDpj(i,jj(3))+piDpj(i,jj(4))
 
340
                xmax = max(absc(piDpj(i,jj(1))),absc(piDpj(i,jj(2))),
 
341
     +                     absc(piDpj(i,jj(3))),absc(piDpj(i,jj(4))))
 
342
                if ( xloss*absc(xheck) .gt. precc*xmax ) print *,
 
343
     +                  'ffcl3p: error: momenta j1234 do not add to 0:',
 
344
     +                  piDpj(i,jj(1)),piDpj(i,jj(2)),piDpj(i,jj(3)),
 
345
     +                  piDpj(i,jj(4)),xheck,i
 
346
                xheck = piDpj(i,jj(6))-piDpj(i,jj(2))-piDpj(i,jj(3))
 
347
                xmax = max(absc(piDpj(i,jj(6))),absc(piDpj(i,jj(2))),
 
348
     +                     absc(piDpj(i,jj(3))))
 
349
                if ( xloss*absc(xheck) .gt. precc*xmax ) print *,
 
350
     +                  'ffcl3p: error: momenta j623 do not add to 0:',
 
351
     +                  piDpj(i,jj(6)),piDpj(i,jj(2)),piDpj(i,jj(3)),
 
352
     +                  xheck,i
 
353
                if ( ns .eq. 10 ) then
 
354
                    xheck = piDpj(i,jj(5))+piDpj(i,jj(1))+piDpj(i,jj(2))
 
355
                else
 
356
                    xheck = piDpj(i,jj(5))-piDpj(i,jj(1))-piDpj(i,jj(2))
 
357
                endif
 
358
                xmax = max(absc(piDpj(i,jj(5))),absc(piDpj(i,jj(1))),
 
359
     +                     absc(piDpj(i,jj(2))))
 
360
                if ( xloss*absc(xheck) .gt. precc*xmax ) print *,
 
361
     +                  'ffcl3p: error: momenta j512 do not add to 0:',
 
362
     +                  piDpj(i,jj(5)),piDpj(i,jj(1)),piDpj(i,jj(2)),
 
363
     +                  xheck,i
 
364
   10       continue
 
365
        endif
 
366
*  #] check input: 
 
367
*  #[ calculations:
 
368
        if ( ii(1).eq.jj(1) .and. ii(2).eq.jj(2) .and. ii(3).eq.jj(3) )
 
369
     +          then
 
370
*
 
371
*           symmetric - fewer possibilities
 
372
*
 
373
            nl = 1
 
374
        else
 
375
            nl = 16
 
376
        endif
 
377
*
 
378
*       try all (1,16)*16 permutations
 
379
*
 
380
        xmax = 0
 
381
        do 101 l=1,nl
 
382
        do 100 i=1,16
 
383
            ii1 = ii(iperm(1,i))
 
384
            ii2 = ii(iperm(2,i))
 
385
            ii3 = ii(iperm(3,i))
 
386
            j = i+l-1
 
387
            if ( j .gt. 16 ) j=j-16
 
388
            jj1 = jj(iperm(1,j))
 
389
            jj2 = jj(iperm(2,j))
 
390
            jj3 = jj(iperm(3,j))
 
391
            s(1) = +piDpj(ii1,jj1)*piDpj(ii2,jj2)*piDpj(ii3,jj3)
 
392
            s(2) = +piDpj(ii2,jj1)*piDpj(ii3,jj2)*piDpj(ii1,jj3)
 
393
            s(3) = +piDpj(ii3,jj1)*piDpj(ii1,jj2)*piDpj(ii2,jj3)
 
394
            s(4) = -piDpj(ii1,jj1)*piDpj(ii3,jj2)*piDpj(ii2,jj3)
 
395
            s(5) = -piDpj(ii3,jj1)*piDpj(ii2,jj2)*piDpj(ii1,jj3)
 
396
            s(6) = -piDpj(ii2,jj1)*piDpj(ii1,jj2)*piDpj(ii3,jj3)
 
397
            som = 0
 
398
            smax = 0
 
399
            do 80 k=1,6
 
400
                som = som + s(k)
 
401
                smax = max(smax,absc(som))
 
402
   80       continue
 
403
            if ( ns .eq. 15 .and. (i.gt.8 .neqv. j.gt.8) )
 
404
     +          som = -som
 
405
            if ( i .eq. 1 .or. smax .lt. xmax ) then
 
406
                dl3p = som
 
407
                xmax = smax
 
408
            endif
 
409
            if ( lwrite ) then
 
410
                print *,'dl3p = +',i-1+16*(l-1),' = ',som,smax
 
411
            endif
 
412
            if ( absc(dl3p) .ge. xloss*smax ) goto 110
 
413
  100   continue
 
414
  101   continue
 
415
        if ( lwarn ) call ffwarn(138,ier,absc(dl3p),xmax)
 
416
  110   continue
 
417
*  #] calculations: 
 
418
*###] ffcl3p: 
 
419
        end