EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
vegas.f
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file vegas.f
1 C*******************************************************************
2 C
3 C
4 C
5 C
6 C*******************************************************************
7 C SUBROUTINE PERFORMS N-DIMENSIONAL MONTE CARLO INTEG'N
8 C - BY G.P. LEPAGE SEPT 1976/(REV)APR 1978
9 C*******************************************************************
10 C
11  SUBROUTINE vegas(FXN,AVGI,SD,CHI2A)
12  IMPLICIT REAL*8(a-h,o-z)
13  common/bveg1/xl(10),xu(10),acc,ndim,ncall,itmx,nprn
14  SAVE /bveg1/
15  common/bveg2/xi(50,10),si,si2,swgt,schi,ndo,it
16  SAVE /bveg2/
17  common/bveg3/f,ti,tsi
18  SAVE /bveg3/
19  EXTERNAL fxn
20  dimension d(50,10),di(50,10),xin(50),r(50),dx(10),dt(10),x(10)
21  1 ,kg(10),ia(10)
22  REAL*4 qran(10)
23  DATA ndmx/50/,alph/1.5d0/,one/1.d0/,mds/-1/
24 
25  SAVE d, di, xin, r, dx, dt, x, kg, ia, qran ! Uzhi
26  SAVE ndmx, alph, one, mds ! Uzhi
27 
28 C
29  ndo=1
30  DO 1 j=1,ndim
31 1 xi(1,j)=one
32 C
33  entry vegas1(fxn,avgi,sd,chi2a)
34 C - INITIALIZES CUMMULATIVE VARIABLES, BUT NOT GRID
35  it=0
36  si=0.
37  si2=si
38  swgt=si
39  schi=si
40 C
41  entry vegas2(fxn,avgi,sd,chi2a)
42 C - NO INITIALIZATION
43  nd=ndmx
44  ng=1
45  IF(mds.EQ.0) go to 2
46  ng=(ncall/2.)**(1./ndim)
47  mds=1
48  IF((2*ng-ndmx).LT.0) go to 2
49  mds=-1
50  npg=ng/ndmx+1
51  nd=ng/npg
52  ng=npg*nd
53 2 k=ng**ndim
54  npg=ncall/k
55  IF(npg.LT.2) npg=2
56  calls=npg*k
57  dxg=one/ng
58  dv2g=(calls*dxg**ndim)**2/npg/npg/(npg-one)
59  xnd=nd
60  ndm=nd-1
61  dxg=dxg*xnd
62  xjac=one/calls
63  DO 3 j=1,ndim
64 c***this is the line 50
65  dx(j)=xu(j)-xl(j)
66 3 xjac=xjac*dx(j)
67 C
68 C REBIN PRESERVING BIN DENSITY
69 C
70  IF(nd.EQ.ndo) go to 8
71  rc=ndo/xnd
72  DO 7 j=1,ndim
73  k=0
74  xn=0.
75  dr=xn
76  i=k
77 4 k=k+1
78  dr=dr+one
79  xo=xn
80  xn=xi(k,j)
81 5 IF(rc.GT.dr) go to 4
82  i=i+1
83  dr=dr-rc
84  xin(i)=xn-(xn-xo)*dr
85  IF(i.LT.ndm) go to 5
86  DO 6 i=1,ndm
87 6 xi(i,j)=xin(i)
88 7 xi(nd,j)=one
89  ndo=nd
90 C
91 8 CONTINUE
92  IF(nprn.NE.0) WRITE(16,200) ndim,calls,it,itmx,acc,mds,nd
93  1 ,(xl(j),xu(j),j=1,ndim)
94 C
95  entry vegas3(fxn,avgi,sd,chi2a)
96 C - MAIN INTEGRATION LOOP
97 9 it=it+1
98  ti=0.
99  tsi=ti
100  DO 10 j=1,ndim
101  kg(j)=1
102  DO 10 i=1,nd
103  d(i,j)=ti
104 10 di(i,j)=ti
105 C
106 11 fb=0.
107  f2b=fb
108  k=0
109 12 k=k+1
110  CALL aran9(qran,ndim)
111  wgt=xjac
112  DO 15 j=1,ndim
113  xn=(kg(j)-qran(j))*dxg+one
114 c*****this is the line 100
115  ia(j)=xn
116  IF(ia(j).GT.1) go to 13
117  xo=xi(ia(j),j)
118  rc=(xn-ia(j))*xo
119  go to 14
120 13 xo=xi(ia(j),j)-xi(ia(j)-1,j)
121  rc=xi(ia(j)-1,j)+(xn-ia(j))*xo
122 14 x(j)=xl(j)+rc*dx(j)
123  wgt=wgt*xo*xnd
124 15 CONTINUE
125 C
126  f=wgt
127  f=f*fxn(x,wgt)
128  f2=f*f
129  fb=fb+f
130  f2b=f2b+f2
131  DO 16 j=1,ndim
132  di(ia(j),j)=di(ia(j),j)+f
133 16 IF(mds.GE.0) d(ia(j),j)=d(ia(j),j)+f2
134  IF(k.LT.npg) go to 12
135 C
136  f2b=dsqrt(f2b*npg)
137  f2b=(f2b-fb)*(f2b+fb)
138  ti=ti+fb
139  tsi=tsi+f2b
140  IF(mds.GE.0) go to 18
141  DO 17 j=1,ndim
142 17 d(ia(j),j)=d(ia(j),j)+f2b
143 18 k=ndim
144 19 kg(k)=mod(kg(k),ng)+1
145  IF(kg(k).NE.1) go to 11
146  k=k-1
147  IF(k.GT.0) go to 19
148 C
149 C FINAL RESULTS FOR THIS ITERATION
150 C
151  tsi=tsi*dv2g
152  ti2=ti*ti
153  wgt=ti2/(tsi+1.0d-37)
154  si=si+ti*wgt
155  si2=si2+ti2
156  swgt=swgt+wgt
157  swgt=swgt+1.0d-37
158  si2=si2+1.0d-37
159  schi=schi+ti2*wgt
160  avgi=si/(swgt)
161  sd=swgt*it/(si2)
162  chi2a=sd*(schi/swgt-avgi*avgi)/(it-.999)
163  sd=dsqrt(one/sd)
164 C****this is the line 150
165  IF(nprn.EQ.0) go to 21
166  tsi=dsqrt(tsi)
167  WRITE(16,201) it,ti,tsi,avgi,sd,chi2a
168  IF(nprn.GE.0) go to 21
169  DO 20 j=1,ndim
170 20 WRITE(16,202) j,(xi(i,j),di(i,j),d(i,j),i=1,nd)
171 C
172 C REFINE GRID
173 C
174 21 DO 23 j=1,ndim
175  xo=d(1,j)
176  xn=d(2,j)
177  d(1,j)=(xo+xn)/2.
178  dt(j)=d(1,j)
179  DO 22 i=2,ndm
180  d(i,j)=xo+xn
181  xo=xn
182  xn=d(i+1,j)
183  d(i,j)=(d(i,j)+xn)/3.
184 22 dt(j)=dt(j)+d(i,j)
185  d(nd,j)=(xn+xo)/2.
186 23 dt(j)=dt(j)+d(nd,j)
187 C
188  DO 28 j=1,ndim
189  rc=0.
190  DO 24 i=1,nd
191  r(i)=0.
192  IF (dt(j).GE.1.0d18) THEN
193  WRITE(6,*) '************** A SINGULARITY >1.0D18'
194 C WRITE(5,1111)
195 C1111 FORMAT(1X,'**************IMPORTANT NOTICE***************')
196 C WRITE(5,1112)
197 C1112 FORMAT(1X,'THE INTEGRAND GIVES RISE A SINGULARITY >1.0D18')
198 C WRITE(5,1113)
199 C1113 FORMAT(1X,'PLEASE CHECK THE INTEGRAND AND THE LIMITS')
200 C WRITE(5,1114)
201 C1114 FORMAT(1X,'**************END NOTICE*************')
202  END IF
203  IF(d(i,j).LE.1.0d-18) go to 24
204  xo=dt(j)/d(i,j)
205  r(i)=((xo-one)/xo/dlog(xo))**alph
206 24 rc=rc+r(i)
207  rc=rc/xnd
208  k=0
209  xn=0.
210  dr=xn
211  i=k
212 25 k=k+1
213  dr=dr+r(k)
214  xo=xn
215 c****this is the line 200
216  xn=xi(k,j)
217 26 IF(rc.GT.dr) go to 25
218  i=i+1
219  dr=dr-rc
220  xin(i)=xn-(xn-xo)*dr/(r(k)+1.0d-30)
221  IF(i.LT.ndm) go to 26
222  DO 27 i=1,ndm
223 27 xi(i,j)=xin(i)
224 28 xi(nd,j)=one
225 C
226  IF(it.LT.itmx.AND.acc*dabs(avgi).LT.sd) go to 9
227 200 FORMAT('0INPUT PARAMETERS FOR VEGAS: NDIM=',i3,' NCALL=',f8.0
228  1 /28x,' IT=',i5,' ITMX=',i5/28x,' ACC=',g9.3
229  2 /28x,' MDS=',i3,' ND=',i4/28x,' (XL,XU)=',
230  3 (t40,'( ',g12.6,' , ',g12.6,' )'))
231 201 FORMAT(///' INTEGRATION BY VEGAS' / '0ITERATION NO.',i3,
232  1 ': INTEGRAL =',g14.8/21x,'STD DEV =',g10.4 /
233  2 ' ACCUMULATED RESULTS: INTEGRAL =',g14.8 /
234  3 24x,'STD DEV =',g10.4 / 24x,'CHI**2 PER IT''N =',g10.4)
235 202 FORMAT('0DATA FOR AXIS',i2 / ' ',6x,'X',7x,' DELT I ',
236  1 2x,' CONV''CE ',11x,'X',7x,' DELT I ',2x,' CONV''CE '
237  2 ,11x,'X',7x,' DELT I ',2x,' CONV''CE ' /
238  2 (' ',3g12.4,5x,3g12.4,5x,3g12.4))
239  RETURN
240  END