[u/mrichter/AliRoot.git] / HIJING / hijing1_36 / vegas.F

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