~maddevelopers/mg5amcnlo/3.0.1

« back to all changes in this revision

Viewing changes to vendor/IREGI/src/matrices.f90

  • Committer: Marco Zaro
  • Date: 2014-01-27 16:54:10 UTC
  • mfrom: (78.124.55 MG5_aMC_2.1)
  • Revision ID: marco.zaro@gmail.com-20140127165410-5lma8c2hzbzm426j
merged with lp:~maddevelopers/madgraph5/MG5_aMC_2.1 r 267

Show diffs side-by-side

added added

removed removed

Lines of Context:
 
1
MODULE matrices
 
2
USE global
 
3
USE funlib
 
4
USE kinematics
 
5
USE matrix_base
 
6
USE cmatrix_base
 
7
USE binary_tree
 
8
IMPLICIT NONE
 
9
CONTAINS
 
10
  SUBROUTINE CaleyMatrix(NLOOPLINE,PCL,M2L,matrix,det)
 
11
    IMPLICIT NONE
 
12
    INTEGER,INTENT(IN)::NLOOPLINE
 
13
    REAL(KIND(1d0)),DIMENSION(1:NLOOPLINE,0:3),INTENT(IN)::PCL
 
14
    REAL(KIND(1d0)),DIMENSION(1:NLOOPLINE),INTENT(IN)::M2L
 
15
    REAL(KIND(1d0)),DIMENSION(0:3)::rij
 
16
    REAL(KIND(1d0)),DIMENSION(NLOOPLINE,NLOOPLINE),INTENT(OUT)::matrix
 
17
    REAL(KIND(1d0)),INTENT(OUT)::det
 
18
    INTEGER::i,j
 
19
    DO i=1,NLOOPLINE
 
20
       DO j=i,NLOOPLINE
 
21
          rij(0:3)=PCL(i,0:3)-PCL(j,0:3)
 
22
          matrix(i,j)=scalarprod(rij(0:3),rij(0:3))-M2L(i)-M2L(j)
 
23
          IF(i.LT.j)THEN
 
24
             matrix(j,i)=matrix(i,j)
 
25
          ENDIF
 
26
       ENDDO
 
27
    ENDDO
 
28
    det=MNXNDET(NLOOPLINE,matrix)
 
29
    RETURN
 
30
  END SUBROUTINE CaleyMatrix
 
31
 
 
32
  SUBROUTINE Complex_CaleyMatrix(NLOOPLINE,PCL,M2L,matrix,det)
 
33
    IMPLICIT NONE
 
34
    INTEGER,INTENT(IN)::NLOOPLINE
 
35
    REAL(KIND(1d0)),DIMENSION(1:NLOOPLINE,0:3),INTENT(IN)::PCL
 
36
    COMPLEX(KIND(1d0)),DIMENSION(1:NLOOPLINE),INTENT(IN)::M2L
 
37
    REAL(KIND(1d0)),DIMENSION(0:3)::rij
 
38
    COMPLEX(KIND(1d0)),DIMENSION(NLOOPLINE,NLOOPLINE),INTENT(OUT)::matrix
 
39
    COMPLEX(KIND(1d0)),INTENT(OUT)::det
 
40
    INTEGER::i,j
 
41
    DO i=1,NLOOPLINE
 
42
       DO j=i,NLOOPLINE
 
43
          rij(0:3)=PCL(i,0:3)-PCL(j,0:3)
 
44
          matrix(i,j)=scalarprod(rij(0:3),rij(0:3))-M2L(i)-M2L(j)
 
45
          IF(i.LT.j)THEN
 
46
             matrix(j,i)=matrix(i,j)
 
47
          ENDIF
 
48
       ENDDO
 
49
    ENDDO
 
50
    det=CMNXNDET(NLOOPLINE,matrix)
 
51
    RETURN
 
52
  END SUBROUTINE Complex_CaleyMatrix
 
53
 
 
54
  SUBROUTINE CaleyMatrix2(NLOOPLINE,PijMatrix,M2L,matrix,det)
 
55
    IMPLICIT NONE
 
56
    INTEGER,INTENT(IN)::NLOOPLINE
 
57
    REAL(KIND(1d0)),DIMENSION(NLOOPLINE,NLOOPLINE),INTENT(IN)::PijMatrix ! PijMatrix(i,j)=Pi.Pj
 
58
    REAL(KIND(1d0)),DIMENSION(1:NLOOPLINE),INTENT(IN)::M2L
 
59
    REAL(KIND(1d0)),DIMENSION(NLOOPLINE,NLOOPLINE),INTENT(OUT)::matrix
 
60
    REAL(KIND(1d0)),INTENT(OUT)::det
 
61
    INTEGER::i,j
 
62
    DO i=1,NLOOPLINE
 
63
       DO j=i,NLOOPLINE
 
64
          matrix(i,j)=PijMatrix(i,i)+PijMatrix(j,j)-2d0*PijMatrix(i,j)-M2L(i)-M2L(j)
 
65
          IF(i.LT.j)THEN
 
66
             matrix(j,i)=matrix(i,j)
 
67
          ENDIF
 
68
       ENDDO
 
69
    ENDDO
 
70
    det=MNXNDET(NLOOPLINE,matrix)
 
71
    RETURN
 
72
  END SUBROUTINE CaleyMatrix2
 
73
 
 
74
  SUBROUTINE GramMatrix(NLOOPLINE,PCL,M2L,matrix,det)
 
75
    IMPLICIT NONE
 
76
    INTEGER,INTENT(IN)::NLOOPLINE
 
77
    REAL(KIND(1d0)),DIMENSION(1:NLOOPLINE,0:3),INTENT(IN)::PCL
 
78
    REAL(KIND(1d0)),DIMENSION(1:NLOOPLINE),INTENT(IN)::M2L
 
79
    REAL(KIND(1d0)),DIMENSION(0:3)::rij
 
80
    REAL(KIND(1d0)),DIMENSION(NLOOPLINE,NLOOPLINE),INTENT(OUT)::matrix
 
81
    REAL(KIND(1d0)),INTENT(OUT)::det
 
82
    REAL(KIND(1d0)),DIMENSION(1:NLOOPLINE,0:3)::PE
 
83
    INTEGER::i,j
 
84
    CALL PCL2PE(NLOOPLINE,PCL,PE)
 
85
    DO i=1,NLOOPLINE
 
86
       DO j=i,NLOOPLINE
 
87
          matrix(i,j)=scalarprod(PE(i,0:3),PE(j,0:3))
 
88
          IF(i.LT.j)THEN
 
89
             matrix(j,i)=matrix(i,j)
 
90
          ENDIF
 
91
       ENDDO
 
92
    ENDDO
 
93
    det=MNXNDET(NLOOPLINE,matrix)
 
94
    RETURN
 
95
  END SUBROUTINE GramMatrix
 
96
  ! hep-ph/0303184
 
97
  SUBROUTINE RSMatrices(NLOOPLINE,PCL,M2L,rmatrix,smatrix,rdet,sdet)
 
98
    IMPLICIT NONE
 
99
    INTEGER,INTENT(IN)::NLOOPLINE
 
100
!    INTEGER,INTENT(IN)::NLOOPLINE
 
101
    REAL(KIND(1d0)),DIMENSION(1:NLOOPLINE,0:3),INTENT(IN)::PCL
 
102
    REAL(KIND(1d0)),DIMENSION(1:NLOOPLINE),INTENT(IN)::M2L
 
103
    REAL(KIND(1d0)),DIMENSION(0:3)::rij
 
104
    REAL(KIND(1d0)),DIMENSION(NLOOPLINE,NLOOPLINE),INTENT(OUT)::rmatrix
 
105
    REAL(KIND(1d0)),DIMENSION(0:NLOOPLINE,0:NLOOPLINE),INTENT(OUT)::smatrix
 
106
    REAL(KIND(1d0)),INTENT(OUT)::rdet,sdet
 
107
    INTEGER::i,j
 
108
    LOGICAL::find
 
109
    TYPE(rsmatrices_node),POINTER::item
 
110
    REAL(KIND(1d0)),DIMENSION(NLOOPLINE,0:3)::PCL1
 
111
    IF(RECYCLING)THEN
 
112
       PCL1(1,0:3)=0d0
 
113
       DO i=2,NLOOPLINE
 
114
          PCL1(i,0:3)=PCL(i,0:3)-PCL(1,0:3)
 
115
       ENDDO
 
116
       ALLOCATE(item)
 
117
       item%NLOOPLINE=NLOOPLINE
 
118
       item%M2L(1:NLOOPLINE)=M2L(1:NLOOPLINE)
 
119
       item%PCL(1:NLOOPLINE,0:3)=PCL1(1:NLOOPLINE,0:3)
 
120
       CALL rsmatrices_bt_search(item,rsmatrices_save,find)
 
121
       IF(find)THEN
 
122
          rdet=item%detR
 
123
          sdet=item%detS
 
124
          RMATRIX(1:NLOOPLINE,1:NLOOPLINE)=&
 
125
               item%RMATRIX(1:NLOOPLINE,1:NLOOPLINE)
 
126
          SMATRIX(0:NLOOPLINE,0:NLOOPLINE)=&
 
127
               item%SMATRIX(0:NLOOPLINE,0:NLOOPLINE)
 
128
          DEALLOCATE(item)
 
129
          RETURN
 
130
       ENDIF
 
131
    ENDIF
 
132
    smatrix(0,0)=0d0
 
133
    DO i=1,NLOOPLINE
 
134
       smatrix(0,i)=1d0
 
135
       smatrix(i,0)=1d0
 
136
    ENDDO
 
137
    DO i=1,NLOOPLINE
 
138
       DO j=i,NLOOPLINE
 
139
          rij(0:3)=PCL(i,0:3)-PCL(j,0:3)
 
140
          rmatrix(i,j)=scalarprod(rij(0:3),rij(0:3))-M2L(i)-M2L(j)
 
141
          smatrix(i,j)=rmatrix(i,j)
 
142
          IF(i.LT.j)THEN
 
143
             rmatrix(j,i)=rmatrix(i,j)
 
144
             smatrix(j,i)=smatrix(i,j)
 
145
          ENDIF
 
146
       ENDDO
 
147
    ENDDO
 
148
    rdet=MNXNDET(NLOOPLINE,rmatrix)
 
149
    sdet=MNXNDET(NLOOPLINE+1,smatrix(0:NLOOPLINE,0:NLOOPLINE))
 
150
    IF(RECYCLING)THEN
 
151
       item%detS=sdet
 
152
       item%detR=rdet
 
153
       item%SMATRIX(0:NLOOPLINE,0:NLOOPLINE)=SMATRIX(0:NLOOPLINE,0:NLOOPLINE)
 
154
       item%RMATRIX(1:NLOOPLINE,1:NLOOPLINE)=RMATRIX(1:NLOOPLINE,1:NLOOPLINE)
 
155
    ENDIF
 
156
    RETURN
 
157
  END SUBROUTINE RSMATRICES
 
158
 
 
159
  SUBROUTINE CRSMatrices(NLOOPLINE,PCL,M2L,rmatrix,smatrix,rdet,sdet)
 
160
    IMPLICIT NONE
 
161
    INTEGER,INTENT(IN)::NLOOPLINE
 
162
    REAL(KIND(1d0)),DIMENSION(1:NLOOPLINE,0:3),INTENT(IN)::PCL
 
163
    COMPLEX(KIND(1d0)),DIMENSION(1:NLOOPLINE),INTENT(IN)::M2L
 
164
    REAL(KIND(1d0)),DIMENSION(0:3)::rij
 
165
    COMPLEX(KIND(1d0)),DIMENSION(NLOOPLINE,NLOOPLINE),INTENT(OUT)::rmatrix
 
166
    COMPLEX(KIND(1d0)),DIMENSION(0:NLOOPLINE,0:NLOOPLINE),INTENT(OUT)::smatrix
 
167
    COMPLEX(KIND(1d0)),INTENT(OUT)::rdet,sdet
 
168
    INTEGER::i,j
 
169
    LOGICAL::find
 
170
    TYPE(crsmatrices_node),POINTER::item
 
171
    REAL(KIND(1d0)),DIMENSION(NLOOPLINE,0:3)::PCL1
 
172
    IF(RECYCLING)THEN
 
173
       PCL1(1,0:3)=0d0
 
174
       DO i=2,NLOOPLINE
 
175
          PCL1(i,0:3)=PCL(i,0:3)-PCL(1,0:3)
 
176
       ENDDO
 
177
       ALLOCATE(item)
 
178
       item%NLOOPLINE=NLOOPLINE
 
179
       item%M2L(1:NLOOPLINE)=M2L(1:NLOOPLINE)
 
180
       item%PCL(1:NLOOPLINE,0:3)=PCL1(1:NLOOPLINE,0:3)
 
181
       CALL crsmatrices_bt_search(item,crsmatrices_save,find)
 
182
       IF(find)THEN
 
183
          rdet=item%detR
 
184
          sdet=item%detS
 
185
          RMATRIX(1:NLOOPLINE,1:NLOOPLINE)=&
 
186
               item%RMATRIX(1:NLOOPLINE,1:NLOOPLINE)
 
187
          SMATRIX(0:NLOOPLINE,0:NLOOPLINE)=&
 
188
               item%SMATRIX(0:NLOOPLINE,0:NLOOPLINE)
 
189
          DEALLOCATE(item)
 
190
          RETURN
 
191
       ENDIF
 
192
    ENDIF
 
193
    smatrix(0,0)=DCMPLX(0d0)
 
194
    DO i=1,NLOOPLINE
 
195
       smatrix(0,i)=DCMPLX(1d0)
 
196
       smatrix(i,0)=DCMPLX(1d0)
 
197
    ENDDO
 
198
    DO i=1,NLOOPLINE
 
199
       DO j=i,NLOOPLINE
 
200
          rij(0:3)=PCL(i,0:3)-PCL(j,0:3)
 
201
          rmatrix(i,j)=scalarprod(rij(0:3),rij(0:3))-M2L(i)-M2L(j)
 
202
          smatrix(i,j)=rmatrix(i,j)
 
203
          IF(i.LT.j)THEN
 
204
             rmatrix(j,i)=rmatrix(i,j)
 
205
             smatrix(j,i)=smatrix(i,j)
 
206
          ENDIF
 
207
       ENDDO
 
208
    ENDDO
 
209
    rdet=CMNXNDET(NLOOPLINE,rmatrix)
 
210
    sdet=CMNXNDET(NLOOPLINE+1,smatrix(0:NLOOPLINE,0:NLOOPLINE))
 
211
    IF(RECYCLING)THEN
 
212
       item%detS=sdet
 
213
       item%detR=rdet
 
214
       item%SMATRIX(0:NLOOPLINE,0:NLOOPLINE)=SMATRIX(0:NLOOPLINE,0:NLOOPLINE)
 
215
       item%RMATRIX(1:NLOOPLINE,1:NLOOPLINE)=RMATRIX(1:NLOOPLINE,1:NLOOPLINE)
 
216
    ENDIF
 
217
    RETURN
 
218
  END SUBROUTINE CRSMATRICES
 
219
 
 
220
  SUBROUTINE RSMatrices2(NLOOPLINE,PijMatrix,M2L,rmatrix,smatrix,rdet,sdet)
 
221
    IMPLICIT NONE
 
222
    INTEGER,INTENT(IN)::NLOOPLINE
 
223
    REAL(KIND(1d0)),DIMENSION(1:NLOOPLINE,1:NLOOPLINE),INTENT(IN)::PijMatrix ! PijMatrix(i,j)=Pi.Pj 
 
224
    REAL(KIND(1d0)),DIMENSION(1:NLOOPLINE),INTENT(IN)::M2L
 
225
    REAL(KIND(1d0)),DIMENSION(NLOOPLINE,NLOOPLINE),INTENT(OUT)::rmatrix
 
226
    REAL(KIND(1d0)),DIMENSION(0:NLOOPLINE,0:NLOOPLINE),INTENT(OUT)::smatrix
 
227
    REAL(KIND(1d0)),INTENT(OUT)::rdet,sdet
 
228
    INTEGER::i,j
 
229
    smatrix(0,0)=0d0
 
230
    DO i=1,NLOOPLINE
 
231
       smatrix(0,i)=1d0
 
232
       smatrix(i,0)=1d0
 
233
    ENDDO
 
234
    DO i=1,NLOOPLINE
 
235
       DO j=i,NLOOPLINE
 
236
          rmatrix(i,j)=PijMatrix(i,i)+PijMatrix(j,j)-2d0*PijMatrix(i,j)-M2L(i)-M2L(j)
 
237
          smatrix(i,j)=rmatrix(i,j)
 
238
          IF(i.LT.j)THEN
 
239
             rmatrix(j,i)=rmatrix(i,j)
 
240
             smatrix(j,i)=smatrix(i,j)
 
241
          ENDIF
 
242
       ENDDO
 
243
    ENDDO
 
244
    rdet=MNXNDET(NLOOPLINE,rmatrix)
 
245
    sdet=MNXNDET(NLOOPLINE+1,smatrix(0:NLOOPLINE,0:NLOOPLINE))
 
246
    RETURN
 
247
  END SUBROUTINE RSMATRICES2
 
248
  ! in hep - ph/0509141
 
249
  SUBROUTINE XYZMATRICES(NLOOPLINE,PCL,M2L,XMATRIX,YMATRIX,ZMATRIX,detY,detZ)
 
250
    IMPLICIT NONE
 
251
    INTEGER,INTENT(IN)::NLOOPLINE
 
252
    REAL(KIND(1d0)),DIMENSION(NLOOPLINE,0:3),INTENT(IN)::PCL
 
253
    REAL(KIND(1d0)),DIMENSION(NLOOPLINE),INTENT(IN)::M2L
 
254
    REAL(KIND(1d0)),DIMENSION(2:NLOOPLINE,2:NLOOPLINE),INTENT(OUT)::ZMATRIX
 
255
    REAL(KIND(1d0)),DIMENSION(NLOOPLINE,NLOOPLINE),INTENT(OUT)::YMATRIX,XMATRIX
 
256
    REAL(KIND(1d0)),INTENT(OUT)::detY,detZ
 
257
    INTEGER::i,j
 
258
    REAL(KIND(1d0)),DIMENSION(0:3)::ri1,rj1,rij
 
259
    LOGICAL::find
 
260
    TYPE(xyzmatrices_node),POINTER::item
 
261
    REAL(KIND(1d0)),DIMENSION(NLOOPLINE,0:3)::PCL1
 
262
    IF(RECYCLING)THEN
 
263
       PCL1(1,0:3)=0d0
 
264
       DO i=2,NLOOPLINE
 
265
          PCL1(i,0:3)=PCL(i,0:3)-PCL(1,0:3)
 
266
       ENDDO
 
267
       ALLOCATE(item)
 
268
       item%NLOOPLINE=NLOOPLINE
 
269
       item%M2L(1:NLOOPLINE)=M2L(1:NLOOPLINE)
 
270
       item%PCL(1:NLOOPLINE,0:3)=PCL1(1:NLOOPLINE,0:3)
 
271
       CALL xyzmatrices_bt_search(item,xyzmatrices_save,find)
 
272
       IF(find)THEN
 
273
          detY=item%detY
 
274
          detZ=item%detZ
 
275
          YMATRIX(1:NLOOPLINE,1:NLOOPLINE)=&
 
276
               item%YMATRIX(1:NLOOPLINE,1:NLOOPLINE)
 
277
          XMATRIX(1:NLOOPLINE,1:NLOOPLINE)=&
 
278
               item%XMATRIX(1:NLOOPLINE,1:NLOOPLINE)
 
279
          ZMATRIX(2:NLOOPLINE,2:NLOOPLINE)=&
 
280
               item%ZMATRIX(2:NLOOPLINE,2:NLOOPLINE)
 
281
          DEALLOCATE(item)
 
282
          RETURN
 
283
       ENDIF
 
284
    ENDIF
 
285
    DO i=1,NLOOPLINE
 
286
       DO j=i,NLOOPLINE
 
287
          rij(0:3)=PCL(i,0:3)-PCL(j,0:3)
 
288
          YMATRIX(i,j)=-scalarprod(rij(0:3),rij(0:3))+M2L(i)+M2L(j)
 
289
          IF(i.GE.2)THEN
 
290
             ri1(0:3)=PCL(i,0:3)-PCL(1,0:3)
 
291
             rj1(0:3)=PCL(j,0:3)-PCL(1,0:3)
 
292
             ZMATRIX(i,j)=2d0*scalarprod(ri1(0:3),rj1(0:3))
 
293
             XMATRIX(i,j)=ZMATRIX(i,j)
 
294
          ELSE
 
295
             IF(j.EQ.1)THEN
 
296
                XMATRIX(i,j)=2*M2L(1)
 
297
             ELSE
 
298
                ri1(0:3)=PCL(j,0:3)-PCL(1,0:3)
 
299
                XMATRIX(i,j)=scalarprod(ri1(0:3),ri1(0:3))-M2L(j)+M2L(1)
 
300
             ENDIF
 
301
          ENDIF
 
302
          IF(i.LT.j)THEN
 
303
             YMATRIX(j,i)=YMATRIX(i,j)
 
304
             XMATRIX(j,i)=XMATRIX(i,j)
 
305
             IF(i.GE.2)THEN
 
306
                ZMATRIX(j,i)=ZMATRIX(i,j)
 
307
             ENDIF
 
308
          ENDIF
 
309
       ENDDO
 
310
    ENDDO
 
311
    detZ=MNXNDET(NLOOPLINE-1,ZMATRIX(2:NLOOPLINE,2:NLOOPLINE))
 
312
    detY=MNXNDET(NLOOPLINE,YMATRIX(1:NLOOPLINE,1:NLOOPLINE))
 
313
    IF(RECYCLING)THEN
 
314
       item%detY=detY
 
315
       item%detZ=detZ
 
316
       item%XMATRIX(1:NLOOPLINE,1:NLOOPLINE)=XMATRIX(1:NLOOPLINE,1:NLOOPLINE)
 
317
       item%YMATRIX(1:NLOOPLINE,1:NLOOPLINE)=YMATRIX(1:NLOOPLINE,1:NLOOPLINE)
 
318
       item%ZMATRIX(2:NLOOPLINE,2:NLOOPLINE)=ZMATRIX(2:NLOOPLINE,2:NLOOPLINE)
 
319
    ENDIF
 
320
    RETURN
 
321
  END SUBROUTINE XYZMATRICES
 
322
 
 
323
  SUBROUTINE CXYZMATRICES(NLOOPLINE,PCL,M2L,XMATRIX,YMATRIX,ZMATRIX,detY,detZ)
 
324
    IMPLICIT NONE
 
325
    INTEGER,INTENT(IN)::NLOOPLINE
 
326
    REAL(KIND(1d0)),DIMENSION(NLOOPLINE,0:3),INTENT(IN)::PCL
 
327
    COMPLEX(KIND(1d0)),DIMENSION(NLOOPLINE),INTENT(IN)::M2L
 
328
    COMPLEX(KIND(1d0)),DIMENSION(2:NLOOPLINE,2:NLOOPLINE),INTENT(OUT)::ZMATRIX
 
329
    COMPLEX(KIND(1d0)),DIMENSION(NLOOPLINE,NLOOPLINE),INTENT(OUT)::YMATRIX,XMATRIX
 
330
    COMPLEX(KIND(1d0)),INTENT(OUT)::detY,detZ
 
331
    INTEGER::i,j
 
332
    REAL(KIND(1d0)),DIMENSION(0:3)::ri1,rj1,rij
 
333
    LOGICAL::find
 
334
    TYPE(cxyzmatrices_node),POINTER::item
 
335
    REAL(KIND(1d0)),DIMENSION(NLOOPLINE,0:3)::PCL1
 
336
    IF(RECYCLING)THEN
 
337
       PCL1(1,0:3)=0d0
 
338
       DO i=2,NLOOPLINE
 
339
          PCL1(i,0:3)=PCL(i,0:3)-PCL(1,0:3)
 
340
       ENDDO
 
341
       ALLOCATE(item)
 
342
       item%NLOOPLINE=NLOOPLINE
 
343
       item%M2L(1:NLOOPLINE)=M2L(1:NLOOPLINE)
 
344
       item%PCL(1:NLOOPLINE,0:3)=PCL1(1:NLOOPLINE,0:3)
 
345
       CALL cxyzmatrices_bt_search(item,cxyzmatrices_save,find)
 
346
       IF(find)THEN
 
347
          detY=item%detY
 
348
          detZ=item%detZ
 
349
          YMATRIX(1:NLOOPLINE,1:NLOOPLINE)=&
 
350
               item%YMATRIX(1:NLOOPLINE,1:NLOOPLINE)
 
351
          XMATRIX(1:NLOOPLINE,1:NLOOPLINE)=&
 
352
               item%XMATRIX(1:NLOOPLINE,1:NLOOPLINE)
 
353
          ZMATRIX(2:NLOOPLINE,2:NLOOPLINE)=&
 
354
               item%ZMATRIX(2:NLOOPLINE,2:NLOOPLINE)
 
355
          DEALLOCATE(item)
 
356
          RETURN
 
357
       ENDIF
 
358
    ENDIF
 
359
    DO i=1,NLOOPLINE
 
360
       DO j=i,NLOOPLINE
 
361
          rij(0:3)=PCL(i,0:3)-PCL(j,0:3)
 
362
          YMATRIX(i,j)=-scalarprod(rij(0:3),rij(0:3))+M2L(i)+M2L(j)
 
363
          IF(i.GE.2)THEN
 
364
             ri1(0:3)=PCL(i,0:3)-PCL(1,0:3)
 
365
             rj1(0:3)=PCL(j,0:3)-PCL(1,0:3)
 
366
             ZMATRIX(i,j)=2d0*scalarprod(ri1(0:3),rj1(0:3))
 
367
             XMATRIX(i,j)=ZMATRIX(i,j)
 
368
          ELSE
 
369
             IF(j.EQ.1)THEN
 
370
                XMATRIX(i,j)=2*M2L(1)
 
371
             ELSE
 
372
                ri1(0:3)=PCL(j,0:3)-PCL(1,0:3)
 
373
                XMATRIX(i,j)=scalarprod(ri1(0:3),ri1(0:3))-M2L(j)+M2L(1)
 
374
             ENDIF
 
375
          ENDIF
 
376
          IF(i.LT.j)THEN
 
377
             YMATRIX(j,i)=YMATRIX(i,j)
 
378
             XMATRIX(j,i)=XMATRIX(i,j)
 
379
             IF(i.GE.2)THEN
 
380
                ZMATRIX(j,i)=ZMATRIX(i,j)
 
381
             ENDIF
 
382
          ENDIF
 
383
       ENDDO
 
384
    ENDDO
 
385
    detZ=CMNXNDET(NLOOPLINE-1,ZMATRIX(2:NLOOPLINE,2:NLOOPLINE))
 
386
    detY=CMNXNDET(NLOOPLINE,YMATRIX(1:NLOOPLINE,1:NLOOPLINE))
 
387
    IF(RECYCLING)THEN
 
388
       item%detY=detY
 
389
       item%detZ=detZ
 
390
       item%XMATRIX(1:NLOOPLINE,1:NLOOPLINE)=XMATRIX(1:NLOOPLINE,1:NLOOPLINE)
 
391
       item%YMATRIX(1:NLOOPLINE,1:NLOOPLINE)=YMATRIX(1:NLOOPLINE,1:NLOOPLINE)
 
392
       item%ZMATRIX(2:NLOOPLINE,2:NLOOPLINE)=ZMATRIX(2:NLOOPLINE,2:NLOOPLINE)
 
393
    ENDIF
 
394
    RETURN
 
395
  END SUBROUTINE CXYZMATRICES
 
396
 
 
397
  SUBROUTINE XYZMATRICES2(NLOOPLINE,PijMatrix,M2L,XMATRIX,YMATRIX,ZMATRIX,detY,detZ)
 
398
    IMPLICIT NONE
 
399
    INTEGER,INTENT(IN)::NLOOPLINE
 
400
    REAL(KIND(1d0)),DIMENSION(NLOOPLINE,NLOOPLINE),INTENT(IN)::PijMatrix ! PijMatrix(i,j)=pi.pj
 
401
    REAL(KIND(1d0)),DIMENSION(NLOOPLINE),INTENT(IN)::M2L
 
402
    REAL(KIND(1d0)),DIMENSION(2:NLOOPLINE,2:NLOOPLINE),INTENT(OUT)::ZMATRIX
 
403
    REAL(KIND(1d0)),DIMENSION(NLOOPLINE,NLOOPLINE),INTENT(OUT)::YMATRIX,XMATRIX
 
404
    REAL(KIND(1d0)),INTENT(OUT)::detY,detZ
 
405
    INTEGER::i,j
 
406
    REAL(KIND(1d0))::ri1,rj1,rij
 
407
    DO i=1,NLOOPLINE
 
408
       DO j=i,NLOOPLINE
 
409
          rij=PijMatrix(i,i)+PijMatrix(j,j)-2d0*PijMatrix(i,j)
 
410
          YMATRIX(i,j)=-rij+M2L(i)+M2L(j)
 
411
          IF(i.GE.2)THEN
 
412
             ZMATRIX(i,j)=2d0*(PijMatrix(i,j)+PijMatrix(1,1)-PijMatrix(i,1)-PijMatrix(j,1))
 
413
             XMATRIX(i,j)=ZMATRIX(i,j)
 
414
          ELSE
 
415
             IF(j.EQ.1)THEN
 
416
                XMATRIX(i,j)=2*M2L(1)
 
417
             ELSE
 
418
                ri1=PijMatrix(i,i)+PijMatrix(1,1)-2d0*PijMatrix(i,1)
 
419
                XMATRIX(i,j)=ri1-M2L(j)+M2L(1)
 
420
             ENDIF
 
421
          ENDIF
 
422
          IF(i.LT.j)THEN
 
423
             YMATRIX(j,i)=YMATRIX(i,j)
 
424
             XMATRIX(j,i)=XMATRIX(i,j)
 
425
             IF(i.GE.2)THEN
 
426
                ZMATRIX(j,i)=ZMATRIX(i,j)
 
427
             ENDIF
 
428
          ENDIF
 
429
       ENDDO
 
430
    ENDDO
 
431
    detZ=MNXNDET(NLOOPLINE-1,ZMATRIX(2:NLOOPLINE,2:NLOOPLINE))
 
432
    detY=MNXNDET(NLOOPLINE,YMATRIX(1:NLOOPLINE,1:NLOOPLINE))
 
433
    RETURN
 
434
  END SUBROUTINE XYZMATRICES2
 
435
 
 
436
END MODULE matrices