2 C*******************************************************************
7 C*******************************************************************
8 C SUBROUTINE PERFORMS N-DIMENSIONAL MONTE CARLO INTEG'N
9 C - BY G.P. LEPAGE SEPT 1976/(REV)APR 1978
10 C*******************************************************************
12 SUBROUTINE VEGAS(FXN,AVGI,SD,CHI2A)
13 IMPLICIT REAL*8(A-H,O-Z)
15 #include "bveg1_hijing.inc"
16 #include "bveg2_hijing.inc"
17 #include "bveg3_hijing.inc"
19 DIMENSION D(50,10),DI(50,10),XIN(50),R(50),DX(10),DT(10),X(10)
22 DATA NDMX/50/,ALPH/1.5D0/,ONE/1.D0/,MDS/-1/
29 ENTRY VEGAS1(FXN,AVGI,SD,CHI2A)
30 C - INITIALIZES CUMMULATIVE VARIABLES, BUT NOT GRID
37 ENTRY VEGAS2(FXN,AVGI,SD,CHI2A)
42 NG=(NCALL/2.)**(1./NDIM)
44 IF((2*NG-NDMX).LT.0) GO TO 2
54 DV2G=(CALLS*DXG**NDIM)**2/NPG/NPG/(NPG-ONE)
60 c***this is the line 50
64 C REBIN PRESERVING BIN DENSITY
77 5 IF(RC.GT.DR) GO TO 4
88 IF(NPRN.NE.0) WRITE(16,200) NDIM,CALLS,IT,ITMX,ACC,MDS,ND
89 1 ,(XL(J),XU(J),J=1,NDIM)
91 ENTRY VEGAS3(FXN,AVGI,SD,CHI2A)
92 C - MAIN INTEGRATION LOOP
106 CALL ARAN9(QRAN,NDIM)
109 XN=(KG(J)-QRAN(J))*DXG+ONE
110 c*****this is the line 100
112 IF(IA(J).GT.1) GO TO 13
116 13 XO=XI(IA(J),J)-XI(IA(J)-1,J)
117 RC=XI(IA(J)-1,J)+(XN-IA(J))*XO
118 14 X(J)=XL(J)+RC*DX(J)
128 DI(IA(J),J)=DI(IA(J),J)+F
129 16 IF(MDS.GE.0) D(IA(J),J)=D(IA(J),J)+F2
130 IF(K.LT.NPG) GO TO 12
133 F2B=(F2B-FB)*(F2B+FB)
136 IF(MDS.GE.0) GO TO 18
138 17 D(IA(J),J)=D(IA(J),J)+F2B
140 19 KG(K)=MOD(KG(K),NG)+1
141 IF(KG(K).NE.1) GO TO 11
145 C FINAL RESULTS FOR THIS ITERATION
149 WGT=TI2/(TSI+1.0d-37)
158 CHI2A=SD*(SCHI/SWGT-AVGI*AVGI)/(IT-.999)
160 C****this is the line 150
161 IF(NPRN.EQ.0) GO TO 21
163 WRITE(16,201) IT,TI,TSI,AVGI,SD,CHI2A
164 IF(NPRN.GE.0) GO TO 21
166 20 WRITE(16,202) J,(XI(I,J),DI(I,J),D(I,J),I=1,ND)
179 D(I,J)=(D(I,J)+XN)/3.
180 22 DT(J)=DT(J)+D(I,J)
182 23 DT(J)=DT(J)+D(ND,J)
188 IF (DT(J).GE.1.0D18) THEN
189 WRITE(6,*) '************** A SINGULARITY >1.0D18'
191 C1111 FORMAT(1X,'**************IMPORTANT NOTICE***************')
193 C1112 FORMAT(1X,'THE INTEGRAND GIVES RISE A SINGULARITY >1.0D18')
195 C1113 FORMAT(1X,'PLEASE CHECK THE INTEGRAND AND THE LIMITS')
197 C1114 FORMAT(1X,'**************END NOTICE*************')
199 IF(D(I,J).LE.1.0D-18) GO TO 24
201 C R(I)=((XO-ONE)/XO/DLOG(XO))**ALPH
202 C The line above doesn't work on Opteron, probably because ALPH is
203 C used in DATA. As a temporary solution we use 1.5D0 directly (PH)
204 R(I)=((XO-ONE)/XO/DLOG(XO))**(1.5D0)
214 c****this is the line 200
216 26 IF(RC.GT.DR) GO TO 25
219 XIN(I)=XN-(XN-XO)*DR/(R(K)+1.0d-30)
220 IF(I.LT.NDM) GO TO 26
225 IF(IT.LT.ITMX.AND.ACC*DABS(AVGI).LT.SD) GO TO 9
226 200 FORMAT('0INPUT PARAMETERS FOR VEGAS: NDIM=',I3,' NCALL=',F8.0
227 1 /28X,' IT=',I5,' ITMX=',I5/28X,' ACC=',G9.3
228 2 /28X,' MDS=',I3,' ND=',I4/28X,' (XL,XU)=',
229 3 (T40,'( ',G12.6,' , ',G12.6,' )'))
230 201 FORMAT(///' INTEGRATION BY VEGAS' / '0ITERATION NO.',I3,
231 1 ': INTEGRAL =',G14.8/21X,'STD DEV =',G10.4 /
232 2 ' ACCUMULATED RESULTS: INTEGRAL =',G14.8 /
233 3 24X,'STD DEV =',G10.4 / 24X,'CHI**2 PER IT''N =',G10.4)
234 202 FORMAT('0DATA FOR AXIS',I2 / ' ',6X,'X',7X,' DELT I ',
235 1 2X,' CONV''CE ',11X,'X',7X,' DELT I ',2X,' CONV''CE '
236 2 ,11X,'X',7X,' DELT I ',2X,' CONV''CE ' /
237 2 (' ',3G12.4,5X,3G12.4,5X,3G12.4))