| e74335a4 |
1 | * $Id$ |
| 2 | C******************************************************************* |
| 3 | C |
| 4 | C |
| 5 | C |
| 6 | 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******************************************************************* |
| 11 | C |
| 12 | SUBROUTINE VEGAS(FXN,AVGI,SD,CHI2A) |
| 13 | IMPLICIT REAL*8(A-H,O-Z) |
| 14 | #include "bveg1_hijing.inc" |
| 15 | #include "bveg2_hijing.inc" |
| 16 | #include "bveg3_hijing.inc" |
| 17 | EXTERNAL FXN |
| 18 | DIMENSION D(50,10),DI(50,10),XIN(50),R(50),DX(10),DT(10),X(10) |
| 19 | 1 ,KG(10),IA(10) |
| 20 | REAL*4 QRAN(10) |
| 21 | DATA NDMX/50/,ALPH/1.5D0/,ONE/1.D0/,MDS/-1/ |
| 22 | SAVE |
| 23 | C |
| 24 | NDO=1 |
| 25 | DO 1 J=1,NDIM |
| 26 | 1 XI(1,J)=ONE |
| 27 | C |
| 28 | ENTRY VEGAS1(FXN,AVGI,SD,CHI2A) |
| 29 | C - INITIALIZES CUMMULATIVE VARIABLES, BUT NOT GRID |
| 30 | IT=0 |
| 31 | SI=0. |
| 32 | SI2=SI |
| 33 | SWGT=SI |
| 34 | SCHI=SI |
| 35 | C |
| 36 | ENTRY VEGAS2(FXN,AVGI,SD,CHI2A) |
| 37 | C - NO INITIALIZATION |
| 38 | ND=NDMX |
| 39 | NG=1 |
| 40 | IF(MDS.EQ.0) GO TO 2 |
| 41 | NG=(NCALL/2.)**(1./NDIM) |
| 42 | MDS=1 |
| 43 | IF((2*NG-NDMX).LT.0) GO TO 2 |
| 44 | MDS=-1 |
| 45 | NPG=NG/NDMX+1 |
| 46 | ND=NG/NPG |
| 47 | NG=NPG*ND |
| 48 | 2 K=NG**NDIM |
| 49 | NPG=NCALL/K |
| 50 | IF(NPG.LT.2) NPG=2 |
| 51 | CALLS=NPG*K |
| 52 | DXG=ONE/NG |
| 53 | DV2G=(CALLS*DXG**NDIM)**2/NPG/NPG/(NPG-ONE) |
| 54 | XND=ND |
| 55 | NDM=ND-1 |
| 56 | DXG=DXG*XND |
| 57 | XJAC=ONE/CALLS |
| 58 | DO 3 J=1,NDIM |
| 59 | c***this is the line 50 |
| 60 | DX(J)=XU(J)-XL(J) |
| 61 | 3 XJAC=XJAC*DX(J) |
| 62 | C |
| 63 | C REBIN PRESERVING BIN DENSITY |
| 64 | C |
| 65 | IF(ND.EQ.NDO) GO TO 8 |
| 66 | RC=NDO/XND |
| 67 | DO 7 J=1,NDIM |
| 68 | K=0 |
| 69 | XN=0. |
| 70 | DR=XN |
| 71 | I=K |
| 72 | 4 K=K+1 |
| 73 | DR=DR+ONE |
| 74 | XO=XN |
| 75 | XN=XI(K,J) |
| 76 | 5 IF(RC.GT.DR) GO TO 4 |
| 77 | I=I+1 |
| 78 | DR=DR-RC |
| 79 | XIN(I)=XN-(XN-XO)*DR |
| 80 | IF(I.LT.NDM) GO TO 5 |
| 81 | DO 6 I=1,NDM |
| 82 | 6 XI(I,J)=XIN(I) |
| 83 | 7 XI(ND,J)=ONE |
| 84 | NDO=ND |
| 85 | C |
| 86 | 8 CONTINUE |
| 87 | IF(NPRN.NE.0) WRITE(16,200) NDIM,CALLS,IT,ITMX,ACC,MDS,ND |
| 88 | 1 ,(XL(J),XU(J),J=1,NDIM) |
| 89 | C |
| 90 | ENTRY VEGAS3(FXN,AVGI,SD,CHI2A) |
| 91 | C - MAIN INTEGRATION LOOP |
| 92 | 9 IT=IT+1 |
| 93 | TI=0. |
| 94 | TSI=TI |
| 95 | DO 10 J=1,NDIM |
| 96 | KG(J)=1 |
| 97 | DO 10 I=1,ND |
| 98 | D(I,J)=TI |
| 99 | 10 DI(I,J)=TI |
| 100 | C |
| 101 | 11 FB=0. |
| 102 | F2B=FB |
| 103 | K=0 |
| 104 | 12 K=K+1 |
| 105 | CALL ARAN9(QRAN,NDIM) |
| 106 | WGT=XJAC |
| 107 | DO 15 J=1,NDIM |
| 108 | XN=(KG(J)-QRAN(J))*DXG+ONE |
| 109 | c*****this is the line 100 |
| 110 | IA(J)=XN |
| 111 | IF(IA(J).GT.1) GO TO 13 |
| 112 | XO=XI(IA(J),J) |
| 113 | RC=(XN-IA(J))*XO |
| 114 | GO TO 14 |
| 115 | 13 XO=XI(IA(J),J)-XI(IA(J)-1,J) |
| 116 | RC=XI(IA(J)-1,J)+(XN-IA(J))*XO |
| 117 | 14 X(J)=XL(J)+RC*DX(J) |
| 118 | WGT=WGT*XO*XND |
| 119 | 15 CONTINUE |
| 120 | C |
| 121 | F=WGT |
| 122 | F=F*FXN(X,WGT) |
| 123 | F2=F*F |
| 124 | FB=FB+F |
| 125 | F2B=F2B+F2 |
| 126 | DO 16 J=1,NDIM |
| 127 | DI(IA(J),J)=DI(IA(J),J)+F |
| 128 | 16 IF(MDS.GE.0) D(IA(J),J)=D(IA(J),J)+F2 |
| 129 | IF(K.LT.NPG) GO TO 12 |
| 130 | C |
| 131 | F2B=DSQRT(F2B*NPG) |
| 132 | F2B=(F2B-FB)*(F2B+FB) |
| 133 | TI=TI+FB |
| 134 | TSI=TSI+F2B |
| 135 | IF(MDS.GE.0) GO TO 18 |
| 136 | DO 17 J=1,NDIM |
| 137 | 17 D(IA(J),J)=D(IA(J),J)+F2B |
| 138 | 18 K=NDIM |
| 139 | 19 KG(K)=MOD(KG(K),NG)+1 |
| 140 | IF(KG(K).NE.1) GO TO 11 |
| 141 | K=K-1 |
| 142 | IF(K.GT.0) GO TO 19 |
| 143 | C |
| 144 | C FINAL RESULTS FOR THIS ITERATION |
| 145 | C |
| 146 | TSI=TSI*DV2G |
| 147 | TI2=TI*TI |
| 148 | WGT=TI2/(TSI+1.0d-37) |
| 149 | SI=SI+TI*WGT |
| 150 | SI2=SI2+TI2 |
| 151 | SWGT=SWGT+WGT |
| 152 | SWGT=SWGT+1.0D-37 |
| 153 | SI2=SI2+1.0D-37 |
| 154 | SCHI=SCHI+TI2*WGT |
| 155 | AVGI=SI/(SWGT) |
| 156 | SD=SWGT*IT/(SI2) |
| 157 | CHI2A=SD*(SCHI/SWGT-AVGI*AVGI)/(IT-.999) |
| 158 | SD=DSQRT(ONE/SD) |
| 159 | C****this is the line 150 |
| 160 | IF(NPRN.EQ.0) GO TO 21 |
| 161 | TSI=DSQRT(TSI) |
| 162 | WRITE(16,201) IT,TI,TSI,AVGI,SD,CHI2A |
| 163 | IF(NPRN.GE.0) GO TO 21 |
| 164 | DO 20 J=1,NDIM |
| 165 | 20 WRITE(16,202) J,(XI(I,J),DI(I,J),D(I,J),I=1,ND) |
| 166 | C |
| 167 | C REFINE GRID |
| 168 | C |
| 169 | 21 DO 23 J=1,NDIM |
| 170 | XO=D(1,J) |
| 171 | XN=D(2,J) |
| 172 | D(1,J)=(XO+XN)/2. |
| 173 | DT(J)=D(1,J) |
| 174 | DO 22 I=2,NDM |
| 175 | D(I,J)=XO+XN |
| 176 | XO=XN |
| 177 | XN=D(I+1,J) |
| 178 | D(I,J)=(D(I,J)+XN)/3. |
| 179 | 22 DT(J)=DT(J)+D(I,J) |
| 180 | D(ND,J)=(XN+XO)/2. |
| 181 | 23 DT(J)=DT(J)+D(ND,J) |
| 182 | C |
| 183 | DO 28 J=1,NDIM |
| 184 | RC=0. |
| 185 | DO 24 I=1,ND |
| 186 | R(I)=0. |
| 187 | IF (DT(J).GE.1.0D18) THEN |
| 188 | WRITE(6,*) '************** A SINGULARITY >1.0D18' |
| 189 | C WRITE(5,1111) |
| 190 | C1111 FORMAT(1X,'**************IMPORTANT NOTICE***************') |
| 191 | C WRITE(5,1112) |
| 192 | C1112 FORMAT(1X,'THE INTEGRAND GIVES RISE A SINGULARITY >1.0D18') |
| 193 | C WRITE(5,1113) |
| 194 | C1113 FORMAT(1X,'PLEASE CHECK THE INTEGRAND AND THE LIMITS') |
| 195 | C WRITE(5,1114) |
| 196 | C1114 FORMAT(1X,'**************END NOTICE*************') |
| 197 | END IF |
| 198 | IF(D(I,J).LE.1.0D-18) GO TO 24 |
| 199 | XO=DT(J)/D(I,J) |
| 200 | R(I)=((XO-ONE)/XO/DLOG(XO))**ALPH |
| 201 | 24 RC=RC+R(I) |
| 202 | RC=RC/XND |
| 203 | K=0 |
| 204 | XN=0. |
| 205 | DR=XN |
| 206 | I=K |
| 207 | 25 K=K+1 |
| 208 | DR=DR+R(K) |
| 209 | XO=XN |
| 210 | c****this is the line 200 |
| 211 | XN=XI(K,J) |
| 212 | 26 IF(RC.GT.DR) GO TO 25 |
| 213 | I=I+1 |
| 214 | DR=DR-RC |
| 215 | XIN(I)=XN-(XN-XO)*DR/(R(K)+1.0d-30) |
| 216 | IF(I.LT.NDM) GO TO 26 |
| 217 | DO 27 I=1,NDM |
| 218 | 27 XI(I,J)=XIN(I) |
| 219 | 28 XI(ND,J)=ONE |
| 220 | C |
| 221 | IF(IT.LT.ITMX.AND.ACC*DABS(AVGI).LT.SD) GO TO 9 |
| 222 | 200 FORMAT('0INPUT PARAMETERS FOR VEGAS: NDIM=',I3,' NCALL=',F8.0 |
| 223 | 1 /28X,' IT=',I5,' ITMX=',I5/28X,' ACC=',G9.3 |
| 224 | 2 /28X,' MDS=',I3,' ND=',I4/28X,' (XL,XU)=', |
| 225 | 3 (T40,'( ',G12.6,' , ',G12.6,' )')) |
| 226 | 201 FORMAT(///' INTEGRATION BY VEGAS' / '0ITERATION NO.',I3, |
| 227 | 1 ': INTEGRAL =',G14.8/21X,'STD DEV =',G10.4 / |
| 228 | 2 ' ACCUMULATED RESULTS: INTEGRAL =',G14.8 / |
| 229 | 3 24X,'STD DEV =',G10.4 / 24X,'CHI**2 PER IT''N =',G10.4) |
| 230 | 202 FORMAT('0DATA FOR AXIS',I2 / ' ',6X,'X',7X,' DELT I ', |
| 231 | 1 2X,' CONV''CE ',11X,'X',7X,' DELT I ',2X,' CONV''CE ' |
| 232 | 2 ,11X,'X',7X,' DELT I ',2X,' CONV''CE ' / |
| 233 | 2 (' ',3G12.4,5X,3G12.4,5X,3G12.4)) |
| 234 | RETURN |
| 235 | END |
git://git.uio.no / u / mrichter / AliRoot.git / blame