VBF: vegas.f

File vegas.f, 7.4 KB (added by trac, 7 years ago)
Line 
1c     second version of vegas with double precision
2c...............
3
4        SUBROUTINE vegas(region,ndim,fxn,init,ncall,itmx,nprn,
5     >    tgral,sd,chi2a)
6        INTEGER init,itmx,ncall,ndim,nprn,NDMX,MXDIM
7        REAL*8 tgral,chi2a,sd,region(2*ndim),fxn,ALPH,TINY
8        PARAMETER (ALPH=1.5,NDMX=50,MXDIM=10,TINY=1.d-30)
9        EXTERNAL fxn
10C       USES fxn,ran2,rebin
11C
12        INTEGER i,idum,it,j,k,mds,nd,ndo,ng,npg,ia(MXDIM),kg(MXDIM)
13      REAL*8 calls,dv2g,dxg,f,f2,f2b,fb,rc,ti,tsi,wgt,xjac,xn,xnd,xo,
14     >       d(NDMX,MXDIM),di(NDMX,MXDIM),dt(MXDIM),dx(MXDIM),
15     >       r(NDMX),x(MXDIM),xi(NDMX,MXDIM),xin(NDMX),ran2
16        DOUBLE PRECISION schi,si,swgt
17        COMMON /ranno/ idum
18        SAVE
19        if(init.le.0)then
20           mds=1
21           ndo=1
22           do 11 j=1,ndim
23              xi(1,j)=1
24 11        enddo
25        endif
26        if (init.le.1) then
27           si=0
28           swgt=0
29           schi=0
30        endif
31        if (init.le.2)then
32           nd=NDMX
33           ng=1
34           if(mds.ne.0)then
35               ng=(ncall/2.d0+0.25d0)**(1.d0/ndim)
36               mds=1
37               if((2*ng-NDMX).ge.0)then
38                  mds=-1
39                  npg=ng/NDMX+1
40                  nd=ng/npg
41                  ng=npg*nd
42               endif
43           endif
44           k=ng**ndim
45           npg=max(ncall/k,2)
46           calls=npg*k
47           dxg=1.d0/ng
48           dv2g=(calls*dxg**ndim)**2/npg/npg/(npg-1.d0)
49           xnd=nd
50           dxg=dxg*xnd
51           xjac=1.d0/calls
52           do 12 j=1, ndim
53              dx(j)=region(j+ndim)-region(j)
54              xjac=xjac*dx(j)
55 12        enddo
56           if(nd.ne.ndo)then
57              do 13 i=1,nd
58                r(i)=1.d0
59 13           enddo
60              do 14 j=1,ndim
61                call rebin(ndo/xnd,nd,r,xin,xi(1,j))
62 14           enddo
63              ndo=nd
64           endif
65           if(nprn.ge.0) write(*,200) ndim,calls,it,itmx,nprn,
66     *         ALPH,mds,nd,(j,region(j),j,region(j+ndim),j=1,ndim)
67        endif
68        do 28 it=1,itmx
69             ti=0.d0
70             tsi=0.d0
71             do 16 j=1,ndim
72                 kg(j)=1
73                 do 15 i=1,nd
74                     d(i,j)=0.d0
75                     di(i,j)=0.d0
76 15              enddo
77 16          enddo
78 10          continue
79                 fb=0.d0
80                 f2b=0.d0
81                 do 19 k=1,npg
82                    wgt=xjac
83                    do 17 j=1,ndim
84                       xn=(kg(j)-ran2(idum))*dxg+1.d0
85                       ia(j)=max(min(int(xn),NDMX),1)
86                       if(ia(j).gt.1)then
87                          xo=xi(ia(j),j)-xi(ia(j)-1,j)
88                          rc=xi(ia(j)-1,j)+(xn-ia(j))*xo
89                       else
90                          xo=xi(ia(j),j)
91                          rc=(xn-ia(j))*xo
92                       endif
93                       x(j)=region(j)+rc*dx(j)
94                       wgt=wgt*xo*xnd
95 17                 enddo
96                    f=wgt*fxn(x,wgt)
97                    f2=f*f
98                    fb=fb+f
99                    f2b=f2b+f2
100                    do 18 j=1,ndim
101                        di(ia(j),j)=di(ia(j),j)+f
102                        if(mds.ge.0) d(ia(j),j)=d(ia(j),j)+f2
103 18                 enddo
104 19              enddo
105                 f2b=sqrt(f2b*npg)
106                 f2b=(f2b-fb)*(f2b+fb)
107                 if (f2b.le.0) f2b=TINY
108                 ti=ti+fb
109                 tsi=tsi+f2b
110                 if(mds.lt.0)then
111                    do 21 j=1,ndim
112                        d(ia(j),j)=d(ia(j),j)+f2b
113 21                 enddo
114                 endif
115             do 22 k=ndim,1,-1
116                 kg(k)=mod(kg(k),ng)+1
117                 if(kg(k).ne.1) goto 10
118 22          enddo
119             tsi=tsi*dv2g
120             wgt=1.d0/tsi
121             si=si+dble(wgt)*dble(ti)
122             schi=schi+dble(wgt)*dble(ti)**2
123             swgt=swgt+dble(wgt)
124             tgral=si/swgt
125             chi2a=max((schi-si*tgral)/(it-.99d0),0.d0)
126             sd=sqrt(1./swgt)
127             tsi=sqrt(tsi)
128             if(nprn.ge.0)then
129                write(*,201) it,ti,tsi,tgral,sd,chi2a
130                if(nprn.ne.0)then
131                   do 23 j=1,ndim
132                      write(*,202) j,(xi(i,j),di(i,j),
133     *                   i=1+nprn/2,nd,nprn)
134 23                enddo
135                endif
136             endif
137             do 25 j=1,ndim
138                xo=d(1,j)
139                xn=d(2,j)
140                d(1,j)=(xo+xn)/2.d0
141                dt(j)=d(1,j)
142                do 24 i=2,nd-1
143                    rc=xo+xn
144                    xo=xn
145                    xn=d(i+1,j)
146                    d(i,j)=(rc+xn)/3.d0
147                    dt(j)=dt(j)+d(i,j)
148 24             enddo
149                d(nd,j)=(xo+xn)/2.d0
150                dt(j)=dt(j)+d(nd,j)
151 25          enddo
152             do 27 j=1,ndim
153                rc=0.d0
154                do 26 i=1,nd
155                    if(d(i,j).lt.TINY) d(i,j)=TINY
156          r(i)=((1.d0-d(i,j)/dt(j))/(log(dt(j))-log(d(i,j))))**ALPH
157                    rc=rc+r(i)
158 26             enddo
159                call rebin(rc/xnd,nd,r,xin,xi(1,j))
160 27          enddo
161 28     enddo
162        return
163 200  FORMAT(/' input parameters for vegas: ndim=',i3,' ncall=',f8.0
164     *   /28x,' it=',i5,' itmx=',i5
165     *   /28x,' nprn=',i3,' alph=',f5.2/28x,' mds=',i3,' nd=',i4
166     *   /(30x,'xl(',i2,')= ',g11.4,' xu(',i2,')= ',g11.4))
167 201  FORMAT(/' iteration no.',I3,': ','integral =',g14.7,'+/- ',g9.2
168     *   /' all iterations: integral =',g14.7,'+/- ',g9.2,
169     *   ' chi**2/it''n =',g9.2)
170 202  FORMAT(/' data for axis ',I2/' X delta i',
171     *  ' x delta i ',' xdelta i ',
172     *  /(1x,f7.5,1x,g11.4,5x,f7.5,1x,g11.4,5x,f7.5,1x,g11.4))
173        END
174
175c...............................................................
176         SUBROUTINE rebin(rc,nd,r,xin,xi)
177         INTEGER nd
178         REAL*8 rc,r(*),xi(*),xin(*)
179         INTEGER i,k
180         REAL*8 dr,xn,xo
181         k=0
182         xn=0.d0
183         dr=0.d0
184         do 11 i=1,nd-1
185 1           if(rc.gt.dr)then
186                k=k+1
187                dr=dr+r(k)
188                xo=xn
189                xn=xi(k)
190             goto 1
191             endif
192             dr=dr-rc
193             xin(i)=xn-(xn-xo)*dr/r(k)
194 11      enddo
195         do 12 i=1,nd-1
196             xi(i)=xin(i)
197 12      enddo
198         xi(nd)=1.d0
199         return
200         END
201c...............................................................
202       FUNCTION ran2(idum)
203       INTEGER idum,IM1,IM2,IMM1,IA1,IA2,IQ1,IQ2,IR1,IR2,NTAB,NDIV
204       REAL*8 ran2,AM,EPS,RNMX
205       PARAMETER (IM1=2147483563,IM2=2147483399,AM=1./IM1,IMM1=IM1-1,
206     *   IA1=40014,IA2=40692,IQ1=53668,IQ2=52774,IR1=12211,
207     *   IR2=3791,NTAB=32,NDIV=1+IMM1/NTAB,EPS=1.2d-7,RNMX=1.d0-EPS)
208       INTEGER idum2,j,k,iv(NTAB),iy
209       SAVE iv,iy,idum2
210       DATA idum2/123456789/, iv/NTAB*0/, iy/0/
211       if (idum.le.0) then
212           idum=max(-idum,1)
213           idum2=idum
214           do 11 j=NTAB+8,1,-1
215               k=idum/IQ1
216               idum=IA1*(idum-k*IQ1)-k*IR1
217               if (idum.lt.0) idum=idum+IM1
218               if (j.le.NTAB) iv(j)=idum
219 11        enddo
220           iy=iv(1)
221       endif
222       k=idum/IQ1
223       idum=IA1*(idum-k*IQ1)-k*IR1
224       if (idum.lt.0) idum=idum+IM1
225       k=idum2/IQ2
226       idum2=IA2*(idum2-k*IQ2)-k*IR2
227       if (idum2.lt.0) idum2=idum2+IM2
228       j=1+iy/NDIV
229       iy=iv(j)-idum2
230       iv(j)=idum
231       if(iy.lt.1) iy=iy+IMM1
232       ran2=min(AM*iy,RNMX)
233       return
234       END
235c....................................................................
236