← Back to branch summary

~madteam/mg5amcnlo/series2.0

« back to all changes in this revision

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

  • Committer: olivier Mattelaer
  • Date: 2015-03-05 00:14:16 UTC
  • mfrom: (258.1.9 2.3)
  • mto: (258.8.1 2.3)
  • mto: This revision was merged to the branch mainline in revision 259.
  • Revision ID: olivier.mattelaer@uclouvain.be-20150305001416-y9mzeykfzwnl9t0j
partial merge

Show diffs side-by-side

added added

removed removed

Lines of Context:
vendor/IREGI/src/gti_reduce.f90
19
19
    LOGICAL,INTENT(OUT)::TSTABLE
20
20
    INTEGER,DIMENSION(1,1)::sol11
21
21
    REAL(KIND(1d0)),DIMENSION(1)::factor1
22
 
    INTEGER::first=0,i,j,numzerp,n,r,nr,init,ntot
 
22
    INTEGER::first=0,i,j,jk,numzerp,n,r,nr,init,ntot
23
23
    INTEGER,DIMENSION(NLOOPLINE)::zerp
24
 
    REAL(KIND(1d0)),DIMENSION(xiarraymax2)::syfactor
25
 
    INTEGER,DIMENSION(xiarraymax2,-1:NLOOPLINE)::sy
 
24
    ! f[i,n]=(i+n-1)!/(n-1)!/i!
 
25
    ! nmax=5,rmax=6,NLOOPLINE=nmax+1
 
26
    ! see syntensor in ti_reduce.f90
 
27
    ! xiarraymax2=(f[0,nmax])+(f[1,nmax])+(f[2,nmax]+f[0,nmax])
 
28
    ! +(f[3,nmax]+f[1,nmax])+(f[4,nmax]+f[2,nmax]+f[0,nmax])
 
29
    ! +(f[5,nmax]+f[3,nmax]+f[1,nmax])+(f[6,nmax]+f[4,nmax]+f[2,nmax]+f[0,nmax])
 
30
    ! when rmax=5,nmax=5 -> xiarraymax2=314
 
31
    ! when rmax=6,nmax=5 -> xiarraymax2=610
 
32
    INTEGER::xiarraymax2,xiarraymax,xiarraymax3
 
33
    !REAL(KIND(1d0)),DIMENSION(xiarraymax2)::syfactor
 
34
    REAL(KIND(1d0)),DIMENSION(:),ALLOCATABLE::syfactor
 
35
    !INTEGER,DIMENSION(xiarraymax2,-1:NLOOPLINE)::sy
 
36
    INTEGER,DIMENSION(:,:),ALLOCATABLE::sy
26
37
!    REAL(KIND(1d0)),PARAMETER::EPS=1d-10
27
38
    REAL(KIND(1d0))::temp
28
39
    INTEGER,DIMENSION(NLOOPLINE)::indices_init2
29
 
    SAVE first
 
40
    SAVE first,xiarraymax,xiarraymax2,xiarraymax3,syfactor,sy
30
41
    MU_R_IREGI=MU
31
42
    IF(first.EQ.0)THEN
32
 
       WRITE(*,*)"#####################################################################################"
33
 
       WRITE(*,*)"#                                                                                   #"
34
 
       WRITE(*,*)"#                           IREGI-alpha-1.0.0                                       #"
35
 
       WRITE(*,*)"#    package for one-loop tensor Integral REduction with General propagator Indices #"
36
 
       WRITE(*,*)"#    Author: Hua-Sheng Shao (huasheng.shao@cern.ch/erdissshaw@gmail.com)            #"
37
 
       WRITE(*,*)"#    a) Physics School, Peking University, Beijing, China                           #"
38
 
       WRITE(*,*)"#    b) PH Department, TH Unit, CERN, Geneva, Switzerland                           #"
39
 
       WRITE(*,*)"#                                                                                   #"
40
 
       WRITE(*,*)"#####################################################################################"
 
43
       IF(.NOT.print_banner)THEN
 
44
          INCLUDE "banner.inc"
 
45
          print_banner=.TRUE.
 
46
       ENDIF
 
47
       xiarraymax2=0
 
48
       DO i=0,MAXRANK_IREGI
 
49
          j=i/2+1
 
50
          DO jk=1,j
 
51
             xiarraymax2=xiarraymax2+&
 
52
                  xiarray_arg1(2*jk-2+MOD(i,2),MAXNLOOP_IREGI-1)
 
53
          ENDDO
 
54
       ENDDO
 
55
       xiarraymax=xiarray_arg1(MAXRANK_IREGI,MAXNLOOP_IREGI-1)
 
56
       xiarraymax3=xiarray_arg1(MAXRANK_IREGI,4)
 
57
       IF(.NOT.ALLOCATED(syfactor))THEN
 
58
          ALLOCATE(syfactor(xiarraymax2))
 
59
       ENDIF
 
60
       IF(.NOT.ALLOCATED(sy))THEN
 
61
          ALLOCATE(sy(xiarraymax2,-1:MAXNLOOP_IREGI))
 
62
       ENDIF
41
63
       ! initialization xiarray and metric, factorial_pair
42
64
       CALL all_Integers(1,1,1,sol11,factor1)
43
65
       CALL calc_factorial_pair
70
92
       ENDIF
71
93
    ENDDO
72
94
    n=NLOOPLINE-numzerp
73
 
    IF(n.GE.6.OR.MAXRANK.GE.6)THEN
74
 
       WRITE(*,*)"ERROR: out of range of general_ti_reduce (N<8,R<7)"
 
95
    IF(n.GE.MAXNLOOP_IREGI.OR.MAXRANK.GT.MAXRANK_IREGI)THEN
 
96
       WRITE(*,100)"ERROR: out of range of general_ti_reduce (N<=",MAXNLOOP_IREGI,",R<=",MAXRANK_IREGI,")"
75
97
       STOP
76
98
    ENDIF
77
 
    CALL general_sytensor(n,indices_init2(1:n),MAXRANK,ntot,sy(1:xiarraymax2,-1:n),syfactor(1:xiarraymax2))
78
 
    CALL general_symmetry(IMODE,NLOOPLINE,idim_init,indices_init,ntot,sy(1:xiarraymax2,-1:n),&
 
99
    CALL general_sytensor(n,indices_init2(1:n),MAXRANK,xiarraymax,xiarraymax2,ntot,&
 
100
         sy(1:xiarraymax2,-1:n),syfactor(1:xiarraymax2))
 
101
    CALL general_symmetry(IMODE,NLOOPLINE,idim_init,indices_init,ntot,xiarraymax,xiarraymax2,xiarraymax3,&
 
102
         sy(1:xiarraymax2,-1:n),&
79
103
         syfactor(1:xiarraymax2),numzerp,zerp,PDEN,M2L,NCOEFS,TICOEFS)
80
104
    TSTABLE=STABLE_IREGI
81
105
    RETURN
 
106
100 FORMAT(2X,A45,I2,A4,I2,A1)
82
107
  END SUBROUTINE general_ti_reduce
83
108
 
84
109
  SUBROUTINE general_ti_reduce2(IMODE,NLOOPLINE,idim_init,indices_init,MAXRANK,NCOEFS,PDEN,PijMatrix,M2L,MU,TICOEFS,TSTABLE)
95
120
    LOGICAL,INTENT(OUT)::TSTABLE
96
121
    INTEGER,DIMENSION(1,1)::sol11
97
122
    REAL(KIND(1d0)),DIMENSION(1)::factor1
98
 
    INTEGER::first=0,i,j,numzerp,n,r,nr,init,ntot
 
123
    INTEGER::first=0,i,j,jk,numzerp,n,r,nr,init,ntot
99
124
    INTEGER,DIMENSION(NLOOPLINE)::zerp
100
 
    REAL(KIND(1d0)),DIMENSION(xiarraymax2)::syfactor
101
 
    INTEGER,DIMENSION(xiarraymax2,-1:NLOOPLINE)::sy
 
125
    ! f[i,n]=(i+n-1)!/(n-1)!/i!
 
126
    ! nmax=5,rmax=6,NLOOPLINE=nmax+1
 
127
    ! see syntensor in ti_reduce.f90
 
128
    ! xiarraymax2=(f[0,nmax])+(f[1,nmax])+(f[2,nmax]+f[0,nmax])
 
129
    ! +(f[3,nmax]+f[1,nmax])+(f[4,nmax]+f[2,nmax]+f[0,nmax])
 
130
    ! +(f[5,nmax]+f[3,nmax]+f[1,nmax])+(f[6,nmax]+f[4,nmax]+f[2,nmax]+f[0,nmax])
 
131
    ! when rmax=5,nmax=5 -> xiarraymax2=314
 
132
    ! when rmax=6,nmax=5 -> xiarraymax2=610
 
133
    INTEGER::xiarraymax2,xiarraymax,xiarraymax3
 
134
    REAL(KIND(1d0)),DIMENSION(:),ALLOCATABLE::syfactor
 
135
    !INTEGER,DIMENSION(xiarraymax2,-1:NLOOPLINE)::sy
 
136
    INTEGER,DIMENSION(:,:),ALLOCATABLE::sy
102
137
    REAL(KIND(1d0))::temp
103
138
    INTEGER,DIMENSION(NLOOPLINE)::indices_init2
104
 
    SAVE first
 
139
    SAVE first,xiarraymax,xiarraymax2,xiarraymax3,syfactor,sy
105
140
    MU_R_IREGI=MU
106
141
    IF(first.EQ.0)THEN
107
 
       WRITE(*,*)"#####################################################################################"
108
 
       WRITE(*,*)"#                                                                                   #"
109
 
       WRITE(*,*)"#                           IREGI-alpha-1.0.0                                       #"
110
 
       WRITE(*,*)"#    package for one-loop tensor Integral REduction with General propagator Indices #"
111
 
       WRITE(*,*)"#    Author: Hua-Sheng Shao (huasheng.shao@cern.ch/erdissshaw@gmail.com)            #"
112
 
       WRITE(*,*)"#    a) Physics School, Peking University, Beijing, China                           #"
113
 
       WRITE(*,*)"#    b) PH Department, TH Unit, CERN, Geneva, Switzerland                           #"
114
 
       WRITE(*,*)"#                                                                                   #"
115
 
       WRITE(*,*)"#####################################################################################"
 
142
       IF(.NOT.print_banner)THEN
 
143
          INCLUDE "banner.inc"
 
144
          print_banner=.TRUE.
 
145
       ENDIF
 
146
       xiarraymax2=0
 
147
       DO i=0,MAXRANK_IREGI
 
148
          j=i/2+1
 
149
          DO jk=1,j
 
150
             xiarraymax2=xiarraymax2+&
 
151
                  xiarray_arg1(2*jk-2+MOD(i,2),MAXNLOOP_IREGI-1)
 
152
          ENDDO
 
153
       ENDDO
 
154
       xiarraymax=xiarray_arg1(MAXRANK_IREGI,MAXNLOOP_IREGI-1)
 
155
       xiarraymax3=xiarray_arg1(MAXRANK_IREGI,4)
 
156
       IF(.NOT.ALLOCATED(syfactor))THEN
 
157
          ALLOCATE(syfactor(xiarraymax2))
 
158
       ENDIF
 
159
       IF(.NOT.ALLOCATED(sy))THEN
 
160
          ALLOCATE(sy(xiarraymax2,-1:MAXNLOOP_IREGI))
 
161
       ENDIF
116
162
       ! initialization xiarray and metric,factorial_pair 
117
163
       CALL all_Integers(1,1,1,sol11,factor1)
118
164
       CALL calc_factorial_pair
145
191
       ENDIF
146
192
    ENDDO
147
193
    n=NLOOPLINE-numzerp
148
 
    IF(n.GE.6.OR.MAXRANK.GE.6)THEN
149
 
       WRITE(*,*)"ERROR: out of range of tensor_integral_reduce2 (N<8,R<7)"
 
194
    IF(n.GE.MAXNLOOP_IREGI.OR.MAXRANK.GT.MAXRANK_IREGI)THEN
 
195
       WRITE(*,100)"ERROR: out of range of general_ti_reduce2 (N<=",MAXNLOOP_IREGI,",R<=",MAXRANK_IREGI,")"
150
196
       STOP
151
197
    ENDIF
152
 
    CALL general_sytensor(n,indices_init2(1:n),MAXRANK,ntot,sy(1:xiarraymax2,-1:n),syfactor(1:xiarraymax2))
153
 
    CALL general_symmetry2(IMODE,NLOOPLINE,idim_init,indices_init,ntot,sy(1:xiarraymax2,-1:n),&
 
198
    CALL general_sytensor(n,indices_init2(1:n),MAXRANK,xiarraymax,xiarraymax2,ntot,&
 
199
         sy(1:xiarraymax2,-1:n),syfactor(1:xiarraymax2))
 
200
    CALL general_symmetry2(IMODE,NLOOPLINE,idim_init,indices_init,ntot,xiarraymax,xiarraymax2,xiarraymax3,&
 
201
         sy(1:xiarraymax2,-1:n),&
154
202
         syfactor(1:xiarraymax2),numzerp,zerp,PDEN,PijMatrix,M2L,NCOEFS,TICOEFS)
155
203
    TSTABLE=STABLE_IREGI
156
204
    RETURN
 
205
100 FORMAT(2X,A46,I2,A4,I2,A1)
157
206
  END SUBROUTINE general_ti_reduce2
158
207
 
159
 
  SUBROUTINE general_sytensor(n,indices_init2,rmax,ntot,sy,syfactor)
 
208
  SUBROUTINE general_sytensor(n,indices_init2,rmax,xiarraymax,xiarraymax2,ntot,sy,syfactor)
160
209
    IMPLICIT NONE
161
 
    INTEGER,INTENT(IN)::n,rmax
 
210
    INTEGER,INTENT(IN)::n,rmax,xiarraymax2,xiarraymax
162
211
    INTEGER,DIMENSION(n),INTENT(IN)::indices_init2
163
212
    INTEGER,INTENT(OUT)::ntot
164
213
    INTEGER,DIMENSION(xiarraymax2,-1:n),INTENT(OUT)::sy
 
214
    !INTEGER,DIMENSION(*,*),INTENT(OUT)::sy
165
215
    INTEGER::r,i,j,k,i0,i0max,rm2x0,nntot,init
166
216
    INTEGER,DIMENSION(xiarraymax,n)::sol
167
217
    REAL(KIND(1d0)),DIMENSION(xiarraymax2),INTENT(OUT)::syfactor
213
263
    RETURN
214
264
  END SUBROUTINE general_sytensor
215
265
 
216
 
  SUBROUTINE general_symmetry2(IMODE,NLOOPLINE,idim_init,indices_init,ntot,sy,syfactor,numzerp,zerp,PDEN,PijMatrix,M2L,NCOEFS,coefs)
 
266
  SUBROUTINE general_symmetry2(IMODE,NLOOPLINE,idim_init,indices_init,ntot,xiarraymax,xiarraymax2,xiarraymax3,&
 
267
       sy,syfactor,numzerp,zerp,PDEN,PijMatrix,M2L,NCOEFS,coefs)
217
268
    ! IMODE=0, IBP reduction
218
269
    ! IMODE=1, PaVe reduction
219
270
    IMPLICIT NONE
220
 
    INTEGER,INTENT(IN)::ntot,numzerp,NLOOPLINE,IMODE,NCOEFS,idim_init
 
271
    INTEGER,INTENT(IN)::ntot,numzerp,NLOOPLINE,IMODE,NCOEFS,idim_init,xiarraymax2,xiarraymax,xiarraymax3
221
272
    INTEGER,DIMENSION(NLOOPLINE),INTENT(IN)::indices_init
222
273
    INTEGER,DIMENSION(xiarraymax2,-1:NLOOPLINE-numzerp),INTENT(IN)::sy
 
274
    !INTEGER,DIMENSION(*,*),INTENT(IN)::sy  
223
275
    REAL(KIND(1d0)),DIMENSION(xiarraymax2),INTENT(IN)::syfactor
224
276
    INTEGER,DIMENSION(NLOOPLINE),INTENT(IN)::zerp
225
277
    REAL(KIND(1d0)),DIMENSION(NLOOPLINE,0:3),INTENT(IN)::PDEN
235
287
    INTEGER::i,j,k,init,r,nntot,nt,nindtot,nindtot0
236
288
    LOGICAL::lzero
237
289
    INTEGER,DIMENSION(0:NLOOPLINE)::nj
238
 
    INTEGER,DIMENSION(0:NLOOPLINE,1:5)::jlist
 
290
    INTEGER,DIMENSION(0:NLOOPLINE,1:MAXRANK_IREGI)::jlist
239
291
    REAL(KIND(1d0))::res,factemp
240
292
    COMPLEX(KIND(1d0)),DIMENSION(1:4)::scalar
 
293
    !REAL(KIND(1d0)),DIMENSION(xiarraymax)::coco
241
294
    REAL(KIND(1d0)),DIMENSION(xiarraymax3)::coco
242
 
    INTEGER,DIMENSION(5)::lor
 
295
    INTEGER,DIMENSION(MAXRANK_IREGI)::lor
243
296
    coefs(0:NCOEFS-1,1:4)=DCMPLX(0d0)
244
297
    j=1
245
298
    nindtot0=0
277
330
             init=init+sol(j,4)
278
331
          ENDIF
279
332
          nj(0:nt)=0
280
 
          jlist(0:nt,1:5)=0
 
333
          jlist(0:nt,1:MAXRANK_IREGI)=0
281
334
          CALL recursive_symmetry(nt,1,r,PCLL(1:nt,0:3),sy(i,-1:nt),lor(1:r),&
282
 
               nj(0:nt),jlist(0:nt,1:5),res)
 
335
               nj(0:nt),jlist(0:nt,1:MAXRANK_IREGI),res)
283
336
          lzero=lzero.AND.(ABS(res).LT.EPS)
284
337
          IF(.NOT.ML5_CONVENTION)THEN
285
338
             coco(j)=res*factor(j)
322
375
    RETURN
323
376
  END SUBROUTINE general_symmetry2
324
377
 
325
 
  SUBROUTINE general_symmetry(IMODE,NLOOPLINE,idim_init,indices_init,ntot,sy,syfactor,numzerp,zerp,PDEN,M2L,NCOEFS,coefs)
 
378
  SUBROUTINE general_symmetry(IMODE,NLOOPLINE,idim_init,indices_init,ntot,xiarraymax,xiarraymax2,xiarraymax3,&
 
379
       sy,syfactor,numzerp,zerp,PDEN,M2L,NCOEFS,coefs)
326
380
    ! IMODE=0, IBP reduction  
327
381
    ! IMODE=1, PaVe reduction 
328
382
    IMPLICIT NONE
329
 
    INTEGER,INTENT(IN)::ntot,numzerp,NLOOPLINE,IMODE,NCOEFS,idim_init
 
383
    INTEGER,INTENT(IN)::ntot,numzerp,NLOOPLINE,IMODE,NCOEFS,idim_init,xiarraymax2,xiarraymax,xiarraymax3
330
384
    INTEGER,DIMENSION(NLOOPLINE),INTENT(IN)::indices_init
331
385
    INTEGER,DIMENSION(xiarraymax2,-1:NLOOPLINE-numzerp),INTENT(IN)::sy
 
386
    !INTEGER,DIMENSION(*,*),INTENT(IN)::sy 
332
387
    REAL(KIND(1d0)),DIMENSION(xiarraymax2),INTENT(IN)::syfactor
333
388
    INTEGER,DIMENSION(NLOOPLINE),INTENT(IN)::zerp
334
389
    REAL(KIND(1d0)),DIMENSION(NLOOPLINE,0:3),INTENT(IN)::PDEN
343
398
    INTEGER::i,j,k,init,r,nntot,nt,nindtot,nindtot0
344
399
    LOGICAL::lzero
345
400
    INTEGER,DIMENSION(0:NLOOPLINE)::nj
346
 
    INTEGER,DIMENSION(0:NLOOPLINE,1:5)::jlist
 
401
    INTEGER,DIMENSION(0:NLOOPLINE,1:MAXRANK_IREGI)::jlist
347
402
    REAL(KIND(1d0))::res,factemp
348
403
    COMPLEX(KIND(1d0)),DIMENSION(1:4)::scalar
 
404
    !REAL(KIND(1d0)),DIMENSION(xiarraymax)::coco
349
405
    REAL(KIND(1d0)),DIMENSION(xiarraymax3)::coco
350
 
    INTEGER,DIMENSION(5)::lor
 
406
    INTEGER,DIMENSION(MAXRANK_IREGI)::lor
351
407
    coefs(0:NCOEFS-1,1:4)=DCMPLX(0d0)
352
408
    j=1
353
409
    nindtot0=0
388
444
             init=init+sol(j,4)
389
445
          ENDIF
390
446
          nj(0:nt)=0
391
 
          jlist(0:nt,1:5)=0
 
447
          jlist(0:nt,1:MAXRANK_IREGI)=0
392
448
          CALL recursive_symmetry(nt,1,r,PCLL(1:nt,0:3),sy(i,-1:nt),lor(1:r),&
393
 
               nj(0:nt),jlist(0:nt,1:5),res)
 
449
               nj(0:nt),jlist(0:nt,1:MAXRANK_IREGI),res)
394
450
          lzero=lzero.AND.(ABS(res).LT.EPS)
395
451
          IF(.NOT.ML5_CONVENTION)THEN
396
452
             coco(j)=res*factor(j)