← Back to branch summary

~mg5core1/mg5amcnlo/2.6.4

« back to all changes in this revision

Viewing changes to 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:
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