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
15 common/bveg2/xi(50,10),si,si2,swgt,schi,ndo,
it
20 dimension d(50,10),di(50,10),xin(50),
r(50),dx(10),
dt(10),
x(10)
23 DATA ndmx/50/,alph/1.5d0/,one/1.d0/,mds/-1/
25 SAVE d, di, xin,
r, dx,
dt,
x, kg, ia, qran
26 SAVE ndmx, alph, one, mds
33 entry vegas1(fxn,avgi,sd,chi2a)
41 entry vegas2(fxn,avgi,sd,chi2a)
46 ng=(ncall/2.)**(1./ndim)
48 IF((2*ng-ndmx).LT.0)
go to 2
58 dv2g=(calls*dxg**ndim)**2/npg/npg/(npg-one)
81 5
IF(rc.GT.dr)
go to 4
92 IF(nprn.NE.0)
WRITE(16,200) ndim,calls,
it,itmx,acc,mds,nd
93 1 ,(xl(
j),xu(
j),
j=1,ndim)
95 entry vegas3(fxn,avgi,sd,chi2a)
110 CALL
aran9(qran,ndim)
113 xn=(kg(
j)-qran(
j))*dxg+one
116 IF(ia(
j).GT.1)
go to 13
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)
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
137 f2b=(f2b-fb)*(f2b+fb)
140 IF(mds.GE.0)
go to 18
142 17 d(ia(
j),
j)=d(ia(
j),
j)+f2b
144 19 kg(
k)=mod(kg(
k),ng)+1
145 IF(kg(
k).NE.1)
go to 11
153 wgt=ti2/(tsi+1.0d-37)
162 chi2a=sd*(schi/swgt-avgi*avgi)/(
it-.999)
165 IF(nprn.EQ.0)
go to 21
167 WRITE(16,201)
it,ti,tsi,avgi,sd,chi2a
168 IF(nprn.GE.0)
go to 21
170 20
WRITE(16,202)
j,(xi(
i,
j),di(
i,
j),d(
i,
j),
i=1,nd)
183 d(
i,
j)=(d(
i,
j)+xn)/3.
192 IF (
dt(
j).GE.1.0d18)
THEN
193 WRITE(6,*)
'************** A SINGULARITY >1.0D18'
203 IF(d(
i,
j).LE.1.0d-18)
go to 24
205 r(
i)=((xo-one)/xo/dlog(xo))**alph
217 26
IF(rc.GT.dr)
go to 25
220 xin(
i)=xn-(xn-xo)*dr/(
r(
k)+1.0d-30)
221 IF(
i.LT.ndm)
go to 26
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 /28
x,
' IT=',
i5,
' ITMX=',
i5/28
x,
' ACC=',g9.3
229 2 /28
x,
' MDS=',
i3,
' ND=',
i4/28
x,
' (XL,XU)=',
230 3 (t40,
'( ',g12.6,
' , ',g12.6,
' )'))
231 201
FORMAT(///
' INTEGRATION BY VEGAS' /
'0ITERATION NO.',
i3,
232 1
': INTEGRAL =',g14.8/21
x,
'STD DEV =',g10.4 /
233 2
' ACCUMULATED RESULTS: INTEGRAL =',g14.8 /
234 3 24
x,
'STD DEV =',g10.4 / 24
x,
'CHI**2 PER IT''N =',g10.4)
235 202
FORMAT(
'0DATA FOR AXIS',
i2 /
' ',6
x,
'X',7
x,
' DELT I ',
236 1 2
x,
' CONV''CE ',11
x,
'X',7
x,
' DELT I ',2
x,
' CONV''CE '
237 2 ,11
x,
'X',7
x,
' DELT I ',2
x,
' CONV''CE ' /
238 2 (
' ',3g12.4,5
x,3g12.4,5
x,3g12.4))