← Back to branch summary

~maddevelopers/mg5amcnlo/3.0.1

« back to all changes in this revision

Viewing changes to MG4_DECAY/vegas.f

move ./decay to ./mg5decay; resolve unit tests (n.b. __init__ does not check keys of input dictionaries, followed last revision)

Show diffs side-by-side

added added

removed removed

Lines of Context:
MG4_DECAY/vegas.f
1
 
      SUBROUTINE vegas(region,ndim,fxn,init,ncall,itmx,nprn,tgral,sd,
2
 
     *chi2a)
3
 
      INTEGER init,itmx,ncall,ndim,nprn,NDMX,MXDIM
4
 
      REAL tgral,chi2a,sd,region(2*ndim),fxn,ALPH,TINY
5
 
      PARAMETER (ALPH=1.5,NDMX=50,MXDIM=10,TINY=1.e-30)
6
 
      EXTERNAL fxn
7
 
CU    USES fxn,xran1,rebin
8
 
      INTEGER i,idum,it,j,k,mds,nd,ndo,ng,npg,ia(MXDIM),kg(MXDIM)
9
 
      REAL calls,dv2g,dxg,f,f2,f2b,fb,rc,ti,tsi,wgt,xjac,xn,xnd,xo,
10
 
     *d(NDMX,MXDIM),di(NDMX,MXDIM),dt(MXDIM),dx(MXDIM),r(NDMX),x(MXDIM),
11
 
     *xi(NDMX,MXDIM),xin(NDMX),xran1,rndn
12
 
      EXTERNAL XRAN1
13
 
      DOUBLE PRECISION schi,si,swgt
14
 
      DATA idum/1/ ! DOES NOT AFFECT PREVIOUS SETTINGS
15
 
c
16
 
C     GLOBAL
17
 
C
18
 
      INTEGER NDIM1
19
 
      REAL    region1(2*MXDIM)
20
 
      COMMON /VEGAS_PAR1/NG,NDIM1
21
 
      COMMON /VEGAS_PAR2/region1,xi,xnd,dx
22
 
C
23
 
      SAVE ! make variables static
24
 
c
25
 
      if(init.le.0)then         !enter here on a cold start
26
 
        NDIM1=NDIM
27
 
        do j=1,2*ndim
28
 
           region1(j)=region(j)
29
 
        enddo
30
 
        mds=1
31
 
        ndo=1
32
 
        do 11 j=1,ndim
33
 
          xi(1,j)=1.
34
 
11      continue
35
 
      endif
36
 
c
37
 
      if (init.le.1)then !enter here to inherit the grid but not the answers
38
 
        si=0D0
39
 
        swgt=0D0
40
 
        schi=0D0
41
 
      endif 
42
 
c
43
 
      if (init.le.2)then !enter here to inherit the grid and the answers
44
 
        nd=NDMX
45
 
        ng=1
46
 
c
47
 
        if(mds.ne.0)then !setup for stratification
48
 
          ng=(ncall/2.+0.25)**(1./ndim) !number of calls per dimension
49
 
          mds=1
50
 
          if((2*ng-NDMX).ge.0)then
51
 
            mds=-1
52
 
            npg=ng/NDMX+1
53
 
            nd=ng/npg
54
 
            ng=npg*nd
55
 
          endif
56
 
        endif
57
 
        k=ng**ndim
58
 
        npg=max(ncall/k,2)
59
 
        calls=npg*k
60
 
cfax: returns the actual number of calls
61
 
        ncall=int(calls)
62
 
        dxg=1./ng
63
 
        dv2g=(calls*dxg**ndim)**2/npg/npg/(npg-1.)
64
 
        xnd=nd
65
 
        dxg=dxg*xnd
66
 
        xjac=1./calls
67
 
c        WRITE (*,*) 'NG,NDIM,K ',NG,NDIM,K
68
 
c        WRITE (*,*) 'NPG,CALLS,NCALL ',NPG,CALLS,NCALL
69
 
c        WRITE (*,*) 'ND,DXG     ',ND,DXG
70
 
        do 12 j=1,ndim
71
 
          dx(j)=region(j+ndim)-region(j)
72
 
          xjac=xjac*dx(j)
73
 
12      continue
74
 
c-- do the binning if necessary
75
 
        if(nd.ne.ndo)then
76
 
          do 13 i=1,nd
77
 
            r(i)=1.
78
 
13        continue
79
 
          do 14 j=1,ndim
80
 
            call rebin(ndo/xnd,nd,r,xin,xi(1,j))
81
 
14        continue
82
 
          ndo=nd
83
 
        endif
84
 
        if(nprn.ge.0) write(*,200) ndim,calls,it,itmx,nprn,ALPH,mds,nd
85
 
      endif 
86
 
c
87
 
c     MAIN ITERATION LOOP
88
 
c     if init>=3 it enters here
89
 
c
90
 
      do 28 it=1,itmx
91
 
        ti=0.
92
 
        tsi=0.
93
 
        do 16 j=1,ndim
94
 
          kg(j)=1
95
 
          do 15 i=1,nd
96
 
            d(i,j)=0.
97
 
            di(i,j)=0.
98
 
15        continue
99
 
16      continue
100
 
10      continue
101
 
          fb=0.
102
 
          f2b=0.
103
 
          do 19 k=1,npg
104
 
            wgt=xjac
105
 
            do 17 j=1,ndim
106
 
cfax: avoid random numbers exactly equal to 0. or 1.
107
 
 303           rndn=xran1(idum)
108
 
               if(rndn.eq.1e0.or.rndn.eq.0e0) goto 303
109
 
              xn=(kg(j)-rndn)*dxg+1.
110
 
c              write(*,*) k,kg(j),xn
111
 
              ia(j)=max(min(int(xn),NDMX),1)
112
 
              if(ia(j).gt.1)then
113
 
                xo=xi(ia(j),j)-xi(ia(j)-1,j)
114
 
                rc=xi(ia(j)-1,j)+(xn-ia(j))*xo
115
 
              else
116
 
                xo=xi(ia(j),j)
117
 
                rc=(xn-ia(j))*xo
118
 
              endif
119
 
              x(j)=region(j)+rc*dx(j)
120
 
              wgt=wgt*xo*xnd
121
 
17          continue
122
 
            f=wgt*fxn(x,wgt)
123
 
            f2=f*f
124
 
            fb=fb+f
125
 
            f2b=f2b+f2
126
 
            do 18 j=1,ndim
127
 
              di(ia(j),j)=di(ia(j),j)+f
128
 
              if(mds.ge.0) d(ia(j),j)=d(ia(j),j)+f2
129
 
18          continue
130
 
19        continue
131
 
          f2b=sqrt(f2b*npg)
132
 
          f2b=(f2b-fb)*(f2b+fb)
133
 
          if (f2b.le.0.) f2b=TINY
134
 
          ti=ti+fb
135
 
          tsi=tsi+f2b
136
 
          if(mds.lt.0)then
137
 
            do 21 j=1,ndim
138
 
              d(ia(j),j)=d(ia(j),j)+f2b
139
 
21          continue
140
 
          endif
141
 
c this is just a counter: total number of points thrown is ng**ndim
142
 
        do 22 k=ndim,1,-1
143
 
          kg(k)=mod(kg(k),ng)+1
144
 
          if(kg(k).ne.1) goto 10
145
 
22      continue
146
 
c
147
 
c     all points for this iteration have been thrown.
148
 
c     now it calculates the results for this iteration
149
 
c     and the cumulative result,chi^2
150
 
c
151
 
        tsi=tsi*dv2g  
152
 
        wgt=1./tsi
153
 
        si=si+dble(wgt)*dble(ti)
154
 
        schi=schi+dble(wgt)*dble(ti)**2
155
 
        swgt=swgt+dble(wgt)
156
 
        tgral=si/swgt
157
 
        chi2a=max((schi-si*tgral)/(it-.99d0),0.d0)
158
 
        sd=sqrt(1./swgt)
159
 
        tsi=sqrt(tsi)
160
 
        if(nprn.ge.0)then
161
 
           write(*,201) it,ti,tsi
162
 
          if(nprn.ne.0)then
163
 
            do 23 j=1,ndim
164
 
              write(*,202) j,(xi(i,j),di(i,j),i=1+nprn/2,nd,nprn)
165
 
23          continue
166
 
          endif
167
 
        endif
168
 
c
169
 
c     refines the grid after the iteration
170
 
c
171
 
        do 25 j=1,ndim
172
 
          xo=d(1,j)
173
 
          xn=d(2,j)
174
 
          d(1,j)=(xo+xn)/2.
175
 
          dt(j)=d(1,j)
176
 
          do 24 i=2,nd-1
177
 
            rc=xo+xn
178
 
            xo=xn
179
 
            xn=d(i+1,j)
180
 
            d(i,j)=(rc+xn)/3.
181
 
            dt(j)=dt(j)+d(i,j)
182
 
24        continue
183
 
          d(nd,j)=(xo+xn)/2.
184
 
          dt(j)=dt(j)+d(nd,j)
185
 
25      continue
186
 
c-- slow down the evolution
187
 
        do 27 j=1,ndim
188
 
          rc=0.
189
 
          do 26 i=1,nd
190
 
            if(d(i,j).lt.TINY) d(i,j)=TINY
191
 
            r(i)=((1.-d(i,j)/dt(j))/(log(dt(j))-log(d(i,j))))**ALPH
192
 
            rc=rc+r(i)
193
 
26        continue
194
 
          call rebin(rc/xnd,nd,r,xin,xi(1,j))
195
 
27      continue
196
 
c
197
 
28    continue !main iteration loop
198
 
c
199
 
      return
200
 
c
201
 
c     print out formats
202
 
c
203
 
200   FORMAT(/' input param. for vegas:  ndim=',i3,'  calls=',
204
 
     *f8.0/24x,'  it=',i5,'  itmx=',i5/24x,'  nprn=',i3,'  alph=',
205
 
     *f5.2/24x,'  mds=',i3,'   nd=',i4/)
206
 
c200   FORMAT(/' input parameters for vegas:  ndim=',i3,'  ncall=',
207
 
c     *f8.0/28x,'  it=',i5,'  itmx=',i5/28x,'  nprn=',i3,'  alph=',
208
 
c     *f5.2/28x,'  mds=',i3,'   nd=',i4/(30x,'xl(',i2,')= ',g11.4,' xu(',
209
 
c     *i2,')= ',g11.4))
210
 
201   FORMAT(' iter. no.',I3,': ','integral =',g14.7,' +/- ',g9.2)
211
 
c201   FORMAT(/' iteration no.',I3,': ','integral =',g14.7,'+/- ',g9.2/
212
 
c     *' all iterations:   integral =',g14.7,'+/- ',g9.2,' chi**2/it',
213
 
c     *'n =',g9.2)
214
 
202   FORMAT(/' data for axis ',I2/'    X       delta i       ',
215
 
     *'   x       delta i       ','    x       delta i       ',/(1x,
216
 
     *f7.5,1x,g11.4,5x,f7.5,1x,g11.4,5x,f7.5,1x,g11.4))
217
 
c
218
 
      END
219
 
 
220
 
 
221
 
 
222
 
      SUBROUTINE rebin(rc,nd,r,xin,xi)
223
 
      INTEGER nd
224
 
      PARAMETER (NDMX=50)
225
 
      REAL rc,r(NDMX),xi(NDMX),xin(NDMX)
226
 
      INTEGER i,k
227
 
      REAL dr,xn,xo
228
 
      k=0
229
 
      xn=0.
230
 
      dr=0.
231
 
      do 11 i=1,nd-1
232
 
1       if(rc.gt.dr)then
233
 
          k=k+1
234
 
          dr=dr+r(k)
235
 
          xo=xn
236
 
          xn=xi(k)
237
 
        goto 1
238
 
        endif
239
 
        dr=dr-rc
240
 
        xin(i)=xn-(xn-xo)*dr/r(k)
241
 
11    continue
242
 
      do 12 i=1,nd-1
243
 
        xi(i)=xin(i)
244
 
12    continue
245
 
      xi(nd)=1.
246
 
      return
247
 
      END