hijing/src/vegas.f

0001 C*******************************************************************
0002 C
0003 C
0004 C
0005 C
0006 C*******************************************************************
0007 C   SUBROUTINE PERFORMS N-DIMENSIONAL MONTE CARLO INTEG'N
0008 C      - BY G.P. LEPAGE   SEPT 1976/(REV)APR 1978
0009 C*******************************************************************
0010 C
0011       SUBROUTINE VEGAS(FXN,AVGI,SD,CHI2A)
0012       IMPLICIT REAL*8(A-H,O-Z)
0013       COMMON/BVEG1/XL(10),XU(10),ACC,NDIM,NCALL,ITMX,NPRN
0014       SAVE  /BVEG1/
0015       COMMON/BVEG2/XI(50,10),SI,SI2,SWGT,SCHI,NDO,IT
0016       SAVE  /BVEG2/
0017       COMMON/BVEG3/F,TI,TSI
0018       SAVE  /BVEG3/
0019       EXTERNAL FXN
0020       DIMENSION D(50,10),DI(50,10),XIN(50),R(50),DX(10),DT(10),X(10)
0021      1   ,KG(10),IA(10)
0022       REAL*4 QRAN(10)
0023       DATA NDMX/50/,ALPH/1.5D0/,ONE/1.D0/,MDS/-1/
0024
0025       SAVE D, DI, XIN, R, DX, DT, X, KG, IA, QRAN      ! Uzhi
0026       SAVE NDMX, ALPH, ONE, MDS                        ! Uzhi
0027
0028 C
0029       NDO=1
0030       DO 1 J=1,NDIM
0031 1     XI(1,J)=ONE
0032 C
0033       ENTRY VEGAS1(FXN,AVGI,SD,CHI2A)
0034 C         - INITIALIZES CUMMULATIVE VARIABLES, BUT NOT GRID
0035       IT=0
0036       SI=0.
0037       SI2=SI
0038       SWGT=SI
0039       SCHI=SI
0040 C
0041       ENTRY VEGAS2(FXN,AVGI,SD,CHI2A)
0042 C         - NO INITIALIZATION
0043       ND=NDMX
0044       NG=1
0045       IF(MDS.EQ.0) GO TO 2
0046       NG=(NCALL/2.)**(1./NDIM)
0047       MDS=1
0048       IF((2*NG-NDMX).LT.0) GO TO 2
0049       MDS=-1
0050       NPG=NG/NDMX+1
0051       ND=NG/NPG
0052       NG=NPG*ND
0053 2     K=NG**NDIM
0054       NPG=NCALL/K
0055       IF(NPG.LT.2) NPG=2
0056       CALLS=NPG*K
0057       DXG=ONE/NG
0058       DV2G=(CALLS*DXG**NDIM)**2/NPG/NPG/(NPG-ONE)
0059       XND=ND
0060       NDM=ND-1
0061       DXG=DXG*XND
0062       XJAC=ONE/CALLS
0063       DO 3 J=1,NDIM
0064 c***this is the line 50
0065       DX(J)=XU(J)-XL(J)
0066 3     XJAC=XJAC*DX(J)
0067 C
0068 C   REBIN PRESERVING BIN DENSITY
0069 C
0070       IF(ND.EQ.NDO) GO TO 8
0071       RC=NDO/XND
0072       DO 7 J=1,NDIM
0073       K=0
0074       XN=0.
0075       DR=XN
0076       I=K
0077 4     K=K+1
0078       DR=DR+ONE
0079       XO=XN
0080       XN=XI(K,J)
0081 5     IF(RC.GT.DR) GO TO 4
0082       I=I+1
0083       DR=DR-RC
0084       XIN(I)=XN-(XN-XO)*DR
0085       IF(I.LT.NDM) GO TO 5
0086       DO 6 I=1,NDM
0087 6     XI(I,J)=XIN(I)
0088 7     XI(ND,J)=ONE
0089       NDO=ND
0090 C
0091 8     CONTINUE
0092       IF(NPRN.NE.0) WRITE(16,200) NDIM,CALLS,IT,ITMX,ACC,MDS,ND
0093      1                           ,(XL(J),XU(J),J=1,NDIM)
0094 C
0095       ENTRY VEGAS3(FXN,AVGI,SD,CHI2A)
0096 C         - MAIN INTEGRATION LOOP
0097 9     IT=IT+1
0098       TI=0.
0099       TSI=TI
0100       DO 10 J=1,NDIM
0101       KG(J)=1
0102       DO 10 I=1,ND
0103       D(I,J)=TI
0104 10    DI(I,J)=TI
0105 C
0106 11    FB=0.
0107       F2B=FB
0108       K=0
0109 12    K=K+1
0110       CALL ARAN9(QRAN,NDIM)
0111       WGT=XJAC
0112       DO 15 J=1,NDIM
0113       XN=(KG(J)-QRAN(J))*DXG+ONE
0114 c*****this is the line 100
0115       IA(J)=XN
0116       IF(IA(J).GT.1) GO TO 13
0117       XO=XI(IA(J),J)
0118       RC=(XN-IA(J))*XO
0119       GO TO 14
0120 13    XO=XI(IA(J),J)-XI(IA(J)-1,J)
0121       RC=XI(IA(J)-1,J)+(XN-IA(J))*XO
0122 14    X(J)=XL(J)+RC*DX(J)
0123       WGT=WGT*XO*XND
0124 15    CONTINUE
0125 C
0126       F=WGT
0127       F=F*FXN(X,WGT)
0128       F2=F*F
0129       FB=FB+F
0130       F2B=F2B+F2
0131       DO 16 J=1,NDIM
0132       DI(IA(J),J)=DI(IA(J),J)+F
0133 16    IF(MDS.GE.0) D(IA(J),J)=D(IA(J),J)+F2
0134       IF(K.LT.NPG) GO TO 12
0135 C
0136       F2B=DSQRT(F2B*NPG)
0137       F2B=(F2B-FB)*(F2B+FB)
0138       TI=TI+FB
0139       TSI=TSI+F2B
0140       IF(MDS.GE.0) GO TO 18
0141       DO 17 J=1,NDIM
0142 17    D(IA(J),J)=D(IA(J),J)+F2B
0143 18    K=NDIM
0144 19    KG(K)=MOD(KG(K),NG)+1
0145       IF(KG(K).NE.1) GO TO 11
0146       K=K-1
0147       IF(K.GT.0) GO TO 19
0148 C
0149 C   FINAL RESULTS FOR THIS ITERATION
0150 C
0151       TSI=TSI*DV2G
0152       TI2=TI*TI
0153       WGT=TI2/(TSI+1.0d-37)
0154       SI=SI+TI*WGT
0155       SI2=SI2+TI2
0156       SWGT=SWGT+WGT
0157       SWGT=SWGT+1.0D-37
0158       SI2=SI2+1.0D-37
0159       SCHI=SCHI+TI2*WGT
0160       AVGI=SI/(SWGT)
0161       SD=SWGT*IT/(SI2)
0162       CHI2A=SD*(SCHI/SWGT-AVGI*AVGI)/(IT-.999)
0163       SD=DSQRT(ONE/SD)
0164 C****this is the line 150
0165       IF(NPRN.EQ.0) GO TO 21
0166       TSI=DSQRT(TSI)
0167       WRITE(16,201) IT,TI,TSI,AVGI,SD,CHI2A
0168       IF(NPRN.GE.0) GO TO 21
0169       DO 20 J=1,NDIM
0170 20    WRITE(16,202) J,(XI(I,J),DI(I,J),D(I,J),I=1,ND)
0171 C
0172 C   REFINE GRID
0173 C
0174 21    DO 23 J=1,NDIM
0175       XO=D(1,J)
0176       XN=D(2,J)
0177       D(1,J)=(XO+XN)/2.
0178       DT(J)=D(1,J)
0179       DO 22 I=2,NDM
0180       D(I,J)=XO+XN
0181       XO=XN
0182       XN=D(I+1,J)
0183       D(I,J)=(D(I,J)+XN)/3.
0184 22    DT(J)=DT(J)+D(I,J)
0185       D(ND,J)=(XN+XO)/2.
0186 23    DT(J)=DT(J)+D(ND,J)
0187 C
0188       DO 28 J=1,NDIM
0189       RC=0.
0190       DO 24 I=1,ND
0191       R(I)=0.
0192       IF (DT(J).GE.1.0D18) THEN
0193        WRITE(6,*) '************** A SINGULARITY >1.0D18'
0194 C      WRITE(5,1111)
0195 C1111  FORMAT(1X,'**************IMPORTANT NOTICE***************')
0196 C      WRITE(5,1112)
0197 C1112  FORMAT(1X,'THE INTEGRAND GIVES RISE A SINGULARITY >1.0D18')
0198 C      WRITE(5,1113)
0199 C1113  FORMAT(1X,'PLEASE CHECK THE INTEGRAND AND THE LIMITS')
0200 C      WRITE(5,1114)
0201 C1114  FORMAT(1X,'**************END NOTICE*************')
0202       END IF
0203       IF(D(I,J).LE.1.0D-18) GO TO 24
0204       XO=DT(J)/D(I,J)
0205       R(I)=((XO-ONE)/XO/DLOG(XO))**ALPH
0206 24    RC=RC+R(I)
0207       RC=RC/XND
0208       K=0
0209       XN=0.
0210       DR=XN
0211       I=K
0212 25    K=K+1
0213       DR=DR+R(K)
0214       XO=XN
0215 c****this is the line 200
0216       XN=XI(K,J)
0217 26    IF(RC.GT.DR) GO TO 25
0218       I=I+1
0219       DR=DR-RC
0220       XIN(I)=XN-(XN-XO)*DR/(R(K)+1.0d-30)
0221       IF(I.LT.NDM) GO TO 26
0222       DO 27 I=1,NDM
0223 27    XI(I,J)=XIN(I)
0224 28    XI(ND,J)=ONE
0225 C
0226       IF(IT.LT.ITMX.AND.ACC*DABS(AVGI).LT.SD) GO TO 9
0227 200   FORMAT('0INPUT PARAMETERS FOR VEGAS:  NDIM=',I3,'  NCALL=',F8.0
0228      1    /28X,'  IT=',I5,'  ITMX=',I5/28X,'  ACC=',G9.3
0229      2    /28X,'  MDS=',I3,'   ND=',I4/28X,'  (XL,XU)=',
0230      3    (T40,'( ',G12.6,' , ',G12.6,' )'))
0231 201   FORMAT(///' INTEGRATION BY VEGAS' / '0ITERATION NO.',I3,
0232      1    ':   INTEGRAL =',G14.8/21X,'STD DEV  =',G10.4 /
0233      2    ' ACCUMULATED RESULTS:   INTEGRAL =',G14.8 /
0234      3    24X,'STD DEV  =',G10.4 / 24X,'CHI**2 PER IT''N =',G10.4)
0235 202   FORMAT('0DATA FOR AXIS',I2 / ' ',6X,'X',7X,'  DELT I  ',
0236      1    2X,' CONV''CE  ',11X,'X',7X,'  DELT I  ',2X,' CONV''CE  '
0237      2   ,11X,'X',7X,'  DELT I  ',2X,' CONV''CE  ' /
0238      2    (' ',3G12.4,5X,3G12.4,5X,3G12.4))
0239       RETURN
0240       END