Vegas integration problem

dani8586 · Oct 15, 2012

Hi,
I have just started to use Vegas for numerical integration and I've found a problem that I can't get rid of...

I'm trying, kust for test, to integrate a naive function, say x**2.
The main function of my program and the vegas function are here in the post.
My problem is the following: with the main and the vegas listed under everithing is working fine, the output are a correct number of iteration with the correct number of calls.
The seed used from the vegas looks like 0, and I found this by putting a print*,iseed just before the end statement at the last line of the subroutine vegas
But if I try to change the seed used in the vegas algorithm, that is to declare my one seed, so if after ncall=1000 i put another declaration like iseed=617647 for example the vegas looks like is stopping after one iteration, and doesn't keep on integrating. The result for that iteration looks correct (is 0.33) but the total result is zero..I can't figure out why, I've tried everything... Can someone help me?

Thanks in advance

**************
*** MAIN ****
**************

program test
implicit real*8 (a-h,o-z)
external fxn

COMMON/BVEG1/NCALL,ITMX,NPRN,NDEV,XL(100),XU(100),ACC
COMMON/BVEG2/IT,NDO,SI,SWGT,SCHI,XI(50,100)
COMMON/BVEG3/ALPH,NDMX,MDS
COMMON/RNDM/ISEED

ndim=1

ncall=1000
itmx=1
nprn=0
acc=1.0d-12
xl(1)=0.0
xu(1)=1.0

call VEGAS(NDIM,FXN,AVGI,SD,CHI2A)
end

function fxn(x,wgt)
implicit real*8 (a-h,o-z)
double precision x(100)

fxn=x(1)**2
return
end

**************
** VEGAS ***
**************

C START VEGAS SECTION
C
BLOCK DATA
IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C MAKES DEFAULT PARAMETER ASSIGNMENTS FOR VEGAS
COMMON/BVEG1/NCALL,ITMX,NPRN,NDEV,XL(100),XU(100),ACC
COMMON/BVEG2/IT,NDO,SI,SWGT,SCHI,XI(50,100)
COMMON/BVEG3/ALPH,NDMX,MDS
DATA NCALL/1000/,ITMX/10/,NPRN/0/,ACC/1.D-2/,
1 XL/100*0.d0/,XU/100*1.d0/,
3 ALPH/1.5d0/,NDMX/50/,MDS/1/,NDEV/6/,
4 NDO/1/,XI/5000*1.d0/,IT/0/,SI,SWGT,SCHI/3*0.d0/
END

SUBROUTINE VEGAS(NDIM,FXN,AVGI,SD,CHI2A)
C
C SUBROUTINE PERFORMS NDIM-DIMENSIONAL MONTE CARLO INTEG'N
C - BY G.P. LEPAGE SEPT 1976/(REV)AUG 1979
C - ALGORITHM DESCRIBED IN J COMP PHYS 27,192(1978)
C
IMPLICIT DOUBLE PRECISION(A-H,O-Z)
COMMON/BVEG1/NCALL,ITMX,NPRN,NDEV,XL(100),XU(100),ACC
COMMON/BVEG2/IT,NDO,SI,SWGT,SCHI,XI(50,100)
COMMON/BVEG3/ALPH,NDMX,MDS
COMMON/BVEG4/CALLS,TI,TSI
COMMON/RESLOCAL/RESL(20),STANDDEVL(20)
DIMENSION D(50,100),DI(50,100),XIN(50),R(50),DX(100),IA(100),
1 KG(100),DT(100),X(100)
DIMENSION RAND(100)
DATA ONE/1.d0/
C
NDO=1
DO 1 J=1,NDIM
1 XI(1,J)=ONE
C
ENTRY VEGAS1(NDIM,FXN,AVGI,SD,CHI2A)
C INITIALIZES CUMMULATIVE VARIABLES, BUT NOT GRID
IT=0
SI=0.d0
SWGT=SI
SCHI=SI
C
ENTRY VEGAS2(NDIM,FXN,AVGI,SD,CHI2A)
C - NO INITIALIZATION
ND=NDMX
NG=1
IF(MDS.EQ.0) GO TO 2
NG=(NCALL/2.d0)**(1.d0/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
DX(J)=XU(J)-XL(J)
3 XJAC=XJAC*DX(J)
C
C REBIN, PRESERVING BIN DENSITY
IF(ND.EQ.NDO) GO TO 8
RC=NDO/XND
DO 7 J=1,NDIM
K=0
XN=0.d0
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 IF(NPRN.GE.0) WRITE(NDEV,200) NDIM,CALLS,IT,ITMX,ACC,NPRN,
1 ALPH,MDS,ND,(XL(J),XU(J),J=1,NDIM)
C
ENTRY VEGAS3(NDIM,FXN,AVGI,SD,CHI2A)
C - MAIN INTEGRATION LOOP
9 IT=IT+1
TI=0.d0
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.d0
F2B=FB
K=0
12 K=K+1
CALL RANDA(NDIM,RAND)
WGT=XJAC
DO 15 J=1,NDIM
XN=(KG(J)-RAND(J))*DXG+ONE
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)
15 WGT=WGT*XO*XND
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=SQRT(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 COMPUTE FINAL RESULTS FOR THIS ITERATION
TSI=TSI*DV2G
TI2=TI*TI
CCCCCCCCCCCCCCCCC MODIFICA 1: PERMETTE INTEGRALE NULLO CCCCCCCCCCCCCCCC`
C
IF(TSI.EQ.0.) THEN
WGT=0.d0
SI=0.d0
SWGT=0.d0
SCHI=0.d0
AVGI=0.d0
CHI2A=0.d0
SD=0.d0
IF(NPRN.LT.0) THEN
CONTINUE
ELSE
TSI=SQRT(TSI)
WRITE(NDEV,201) IT,TI,TSI,AVGI,SD,CHI2A
IF(NPRN.EQ.0) THEN
CONTINUE
ELSE
DO J=1,NDIM
WRITE(NDEV,202) J,(XI(I,J),DI(I,J),I=1,ND,NPRN)
END DO
END IF
END IF
RETURN
END IF
C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC`
C
C PROTECT AGAINST ACCIDENTAL SINGLE PRECISION NEGATIVE VALUES
XXX=TSI/TI**2
IF(XXX.GT.-1.D-4) F2B=ABS(F2B)

WGT=ONE/TSI
SI=SI+TI*WGT
SWGT=SWGT+WGT
SCHI=SCHI+TI2*WGT
AVGI=SI/SWGT
CHI2A=(SCHI-SI*AVGI)/(IT-.9999d0)
SD=SQRT(ONE/SWGT)
C
IF(NPRN.LT.0) GO TO 21
TSI=SQRT(TSI)
WRITE(NDEV,201) IT,TI,TSI,AVGI,SD,CHI2A
RESL(IT)=TI
STANDDEVL(IT)=TSI
IF(NPRN.EQ.0) GO TO 21
DO 20 J=1,NDIM
20 WRITE(NDEV,202) J,(XI(I,J),DI(I,J),I=1,ND,NPRN)
C
C REFINE GRID
21 DO 23 J=1,NDIM
XO=D(1,J)
XN=D(2,J)
D(1,J)=(XO+XN)/2.d0
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.d0
22 DT(J)=DT(J)+D(I,J)
D(ND,J)=(XN+XO)/2.d0
23 DT(J)=DT(J)+D(ND,J)
C
DO 28 J=1,NDIM
RC=0.d0
DO 24 I=1,ND
R(I)=0.d0
IF(D(I,J).LE.0.d0) GO TO 24
XO=DT(J)/D(I,J)
R(I)=((XO-ONE)/XO/LOG(XO))**ALPH
24 RC=RC+R(I)
RC=RC/XND
K=0
XN=0.d0
DR=XN
I=K
25 K=K+1
DR=DR+R(K)
XO=XN
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)
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*ABS(AVGI).LT.SD) GO TO 9
200 FORMAT(/35H INPUT PARAMETERS FOR VEGAS: NDIM=,I3,8H NCALL=,F8.0
1 /28X,5H IT=,I5,7H ITMX=,I5/28X,6H ACC=,G9.3
2 /28X,7H NPRN=,I3,7H ALPH=,F5.2/28X,6H MDS=,I3,6H ND=,I4
3 /28X,10H (XL,XU)=,(T40,2H( ,G12.6,3H , ,G12.6,2H )))
201 FORMAT(///21H INTEGRATION BY VEGAS//14H ITERATION NO.,I3,
1 14H: INTEGRAL =,G14.8/21X,10HSTD DEV =,G10.4/
2 34H ACCUMULATED RESULTS: INTEGRAL =,G14.8/
3 24X,10HSTD DEV =,G10.4/24X,17HCHI**2 PER IT'N =,G10.4)
202 FORMAT(/15H DATA FOR AXIS ,I2/25H X DELTA I ,
1 24H X DELTA I ,18H X DELTA I
2 /(1H ,F7.6,1X,G11.4,5X,F7.6,1X,G11.4,5X,F7.6,1X,G11.4))
RETURN
END

C***************************************************************
C LOAD VEGAS DATA IF DESIRED
SUBROUTINE LOAD_VEGAS(NDIM,NAME)
IMPLICIT DOUBLE PRECISION(A-H,O-Z)
COMMON/BVEG2/IT,NDO,SI,SWGT,SCHI,XI(50,100)
CHARACTER*30 NAME
C OPEN(UNIT=13,FILE=NAME,STATUS='OLD',SHARED)
READ(13,210) IT,NDO,SI,SWGT,SCHI
DO 190 J=1,NDIM
190 READ(13,*) (XI(I,J),I=1,50)
210 FORMAT(2I8,3Z16)
CLOSE(UNIT=13)
RETURN
END

C*******************************************************************
C STORE VEGAS DATA FOR POSSIBLE LATER USE
SUBROUTINE STORE_VEGAS(NDIM,NAME)
IMPLICIT DOUBLE PRECISION(A-H,O-Z)
COMMON/BVEG2/IT,NDO,SI,SWGT,SCHI,XI(50,100)
CHARACTER*30 NAME
OPEN(UNIT=12,FILE=NAME,STATUS='NEW')
WRITE(12,210) IT,NDO,SI,SWGT,SCHI
DO 190 J=1,NDIM
190 WRITE(12,*) (XI(I,J),I=1,50)
210 FORMAT(2I8,3Z16)
CLOSE(UNIT=12)
RETURN
END

SUBROUTINE RANDA(N,RAND)
C
C SUBROUTINE GENERATES UNIFORMLY DISTRIBUTED RANDOM NO'S X(I),I=1,N
IMPLICIT DOUBLE PRECISION(A-H,O-Z)
COMMON/RNDM/ISEED
DIMENSION RAND(N)
DO 1 I=1,N
1 RAND(I)=RAN(ISEED)
RETURN
END

C--------- END THE VEGAS SECTION -------------------------------------

mikrom · Oct 15, 2012

What's not nice, that you initialize global variables on 2 places:
in your main program

Code:

...
ncall=1000
...

and in BLOCK DATA section

Code:

BLOCK DATA
IMPLICIT DOUBLE PRECISION(A-H,O-Z)
COMMON/BVEG1/NCALL,ITMX,NPRN,NDEV,XL(100),XU(100),ACC
...
DATA [COLOR=red]NCALL/1000/[/color],ITMX/10/,NPRN/0/,ACC/1.D-2/,
...
END

If you would post your code properly formatted in [ignore]

Code:

..

[/ignore] marks, maybe someone compile it and looks at it.

dani8586 · Oct 16, 2012

Hi, thanks! Do you mean formatted in this way?
Ah, I use gfortran to compile the code

Thanks

Code:

**************
*** MAIN ****
**************

program test
implicit real*8 (a-h,o-z)
external fxn


COMMON/BVEG1/NCALL,ITMX,NPRN,NDEV,XL(100),XU(100),ACC
COMMON/BVEG2/IT,NDO,SI,SWGT,SCHI,XI(50,100)
COMMON/BVEG3/ALPH,NDMX,MDS
COMMON/RNDM/ISEED

ndim=1

ncall=1000
itmx=1
nprn=0
acc=1.0d-12
xl(1)=0.0
xu(1)=1.0

call VEGAS(NDIM,FXN,AVGI,SD,CHI2A)
end



function fxn(x,wgt)
implicit real*8 (a-h,o-z)
double precision x(100)

fxn=x(1)**2
return
end

**************
** VEGAS ***
**************

C START VEGAS SECTION
C
BLOCK DATA
IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C MAKES DEFAULT PARAMETER ASSIGNMENTS FOR VEGAS
COMMON/BVEG1/NCALL,ITMX,NPRN,NDEV,XL(100),XU(100),ACC
COMMON/BVEG2/IT,NDO,SI,SWGT,SCHI,XI(50,100)
COMMON/BVEG3/ALPH,NDMX,MDS
DATA NCALL/1000/,ITMX/10/,NPRN/0/,ACC/1.D-2/,
1 XL/100*0.d0/,XU/100*1.d0/,
3 ALPH/1.5d0/,NDMX/50/,MDS/1/,NDEV/6/,
4 NDO/1/,XI/5000*1.d0/,IT/0/,SI,SWGT,SCHI/3*0.d0/
END


SUBROUTINE VEGAS(NDIM,FXN,AVGI,SD,CHI2A)
C
C SUBROUTINE PERFORMS NDIM-DIMENSIONAL MONTE CARLO INTEG'N
C - BY G.P. LEPAGE SEPT 1976/(REV)AUG 1979
C - ALGORITHM DESCRIBED IN J COMP PHYS 27,192(1978)
C
IMPLICIT DOUBLE PRECISION(A-H,O-Z)
COMMON/BVEG1/NCALL,ITMX,NPRN,NDEV,XL(100),XU(100),ACC
COMMON/BVEG2/IT,NDO,SI,SWGT,SCHI,XI(50,100)
COMMON/BVEG3/ALPH,NDMX,MDS
COMMON/BVEG4/CALLS,TI,TSI
COMMON/RESLOCAL/RESL(20),STANDDEVL(20)
DIMENSION D(50,100),DI(50,100),XIN(50),R(50),DX(100),IA(100),
1 KG(100),DT(100),X(100)
DIMENSION RAND(100)
DATA ONE/1.d0/
C	
NDO=1
DO 1 J=1,NDIM
1 XI(1,J)=ONE
C
ENTRY VEGAS1(NDIM,FXN,AVGI,SD,CHI2A)
C INITIALIZES CUMMULATIVE VARIABLES, BUT NOT GRID
IT=0
SI=0.d0
SWGT=SI
SCHI=SI
C
ENTRY VEGAS2(NDIM,FXN,AVGI,SD,CHI2A)
C - NO INITIALIZATION
ND=NDMX
NG=1
IF(MDS.EQ.0) GO TO 2
NG=(NCALL/2.d0)**(1.d0/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
DX(J)=XU(J)-XL(J)
3 XJAC=XJAC*DX(J)
C
C REBIN, PRESERVING BIN DENSITY
IF(ND.EQ.NDO) GO TO 8
RC=NDO/XND
DO 7 J=1,NDIM
K=0
XN=0.d0
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 IF(NPRN.GE.0) WRITE(NDEV,200) NDIM,CALLS,IT,ITMX,ACC,NPRN,
1 ALPH,MDS,ND,(XL(J),XU(J),J=1,NDIM)
C
ENTRY VEGAS3(NDIM,FXN,AVGI,SD,CHI2A)
C - MAIN INTEGRATION LOOP
9 IT=IT+1
TI=0.d0
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.d0
F2B=FB
K=0
12 K=K+1
CALL RANDA(NDIM,RAND)
WGT=XJAC
DO 15 J=1,NDIM
XN=(KG(J)-RAND(J))*DXG+ONE
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)
15 WGT=WGT*XO*XND
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=SQRT(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 COMPUTE FINAL RESULTS FOR THIS ITERATION
TSI=TSI*DV2G
TI2=TI*TI
CCCCCCCCCCCCCCCCC MODIFICA 1: PERMETTE INTEGRALE NULLO CCCCCCCCCCCCCCCC`
C
IF(TSI.EQ.0.) THEN
WGT=0.d0
SI=0.d0
SWGT=0.d0
SCHI=0.d0
AVGI=0.d0
CHI2A=0.d0
SD=0.d0
IF(NPRN.LT.0) THEN
CONTINUE
ELSE
TSI=SQRT(TSI) 
WRITE(NDEV,201) IT,TI,TSI,AVGI,SD,CHI2A 
IF(NPRN.EQ.0) THEN
CONTINUE
ELSE
DO J=1,NDIM
WRITE(NDEV,202) J,(XI(I,J),DI(I,J),I=1,ND,NPRN) 
END DO
END IF
END IF
RETURN
END IF
C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC`
C
C PROTECT AGAINST ACCIDENTAL SINGLE PRECISION NEGATIVE VALUES
XXX=TSI/TI**2
IF(XXX.GT.-1.D-4) F2B=ABS(F2B)

WGT=ONE/TSI
SI=SI+TI*WGT
SWGT=SWGT+WGT
SCHI=SCHI+TI2*WGT
AVGI=SI/SWGT
CHI2A=(SCHI-SI*AVGI)/(IT-.9999d0)
SD=SQRT(ONE/SWGT)
C
IF(NPRN.LT.0) GO TO 21
TSI=SQRT(TSI)
WRITE(NDEV,201) IT,TI,TSI,AVGI,SD,CHI2A
RESL(IT)=TI
STANDDEVL(IT)=TSI
IF(NPRN.EQ.0) GO TO 21
DO 20 J=1,NDIM
20 WRITE(NDEV,202) J,(XI(I,J),DI(I,J),I=1,ND,NPRN)
C
C REFINE GRID
21 DO 23 J=1,NDIM
XO=D(1,J)
XN=D(2,J)
D(1,J)=(XO+XN)/2.d0
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.d0
22 DT(J)=DT(J)+D(I,J)
D(ND,J)=(XN+XO)/2.d0
23 DT(J)=DT(J)+D(ND,J)
C
DO 28 J=1,NDIM
RC=0.d0
DO 24 I=1,ND
R(I)=0.d0
IF(D(I,J).LE.0.d0) GO TO 24
XO=DT(J)/D(I,J)
R(I)=((XO-ONE)/XO/LOG(XO))**ALPH
24 RC=RC+R(I)
RC=RC/XND
K=0
XN=0.d0
DR=XN
I=K
25 K=K+1
DR=DR+R(K)
XO=XN
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)
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*ABS(AVGI).LT.SD) GO TO 9
200 FORMAT(/35H INPUT PARAMETERS FOR VEGAS: NDIM=,I3,8H NCALL=,F8.0
1 /28X,5H IT=,I5,7H ITMX=,I5/28X,6H ACC=,G9.3
2 /28X,7H NPRN=,I3,7H ALPH=,F5.2/28X,6H MDS=,I3,6H ND=,I4
3 /28X,10H (XL,XU)=,(T40,2H( ,G12.6,3H , ,G12.6,2H )))
201 FORMAT(///21H INTEGRATION BY VEGAS//14H ITERATION NO.,I3,
1 14H: INTEGRAL =,G14.8/21X,10HSTD DEV =,G10.4/
2 34H ACCUMULATED RESULTS: INTEGRAL =,G14.8/
3 24X,10HSTD DEV =,G10.4/24X,17HCHI**2 PER IT'N =,G10.4)
202 FORMAT(/15H DATA FOR AXIS ,I2/25H X DELTA I ,
1 24H X DELTA I ,18H X DELTA I
2 /(1H ,F7.6,1X,G11.4,5X,F7.6,1X,G11.4,5X,F7.6,1X,G11.4))
RETURN
END

C***************************************************************
C LOAD VEGAS DATA IF DESIRED
SUBROUTINE LOAD_VEGAS(NDIM,NAME)
IMPLICIT DOUBLE PRECISION(A-H,O-Z)
COMMON/BVEG2/IT,NDO,SI,SWGT,SCHI,XI(50,100)
CHARACTER*30 NAME
C OPEN(UNIT=13,FILE=NAME,STATUS='OLD',SHARED)
READ(13,210) IT,NDO,SI,SWGT,SCHI
DO 190 J=1,NDIM
190 READ(13,*) (XI(I,J),I=1,50)
210 FORMAT(2I8,3Z16)
CLOSE(UNIT=13)
RETURN
END

C*******************************************************************
C STORE VEGAS DATA FOR POSSIBLE LATER USE
SUBROUTINE STORE_VEGAS(NDIM,NAME)
IMPLICIT DOUBLE PRECISION(A-H,O-Z)
COMMON/BVEG2/IT,NDO,SI,SWGT,SCHI,XI(50,100)
CHARACTER*30 NAME
OPEN(UNIT=12,FILE=NAME,STATUS='NEW')
WRITE(12,210) IT,NDO,SI,SWGT,SCHI
DO 190 J=1,NDIM
190 WRITE(12,*) (XI(I,J),I=1,50)
210 FORMAT(2I8,3Z16)
CLOSE(UNIT=12)
RETURN
END

SUBROUTINE RANDA(N,RAND)
C
C SUBROUTINE GENERATES UNIFORMLY DISTRIBUTED RANDOM NO'S X(I),I=1,N
IMPLICIT DOUBLE PRECISION(A-H,O-Z)
COMMON/RNDM/ISEED
DIMENSION RAND(N)
DO 1 I=1,N
1 RAND(I)=RAN(ISEED)
RETURN
END

C--------- END THE VEGAS SECTION -------------------------------------

mikrom · Oct 16, 2012

dani8586 said:
Ah, I use gfortran to compile the code

IMO it's not properly formatted. I tried it, but it doesn't compile. What's it fixed form or free form? Which switches are you using to compile it?

gummibaer · Oct 16, 2012

You can see, it is a fixed format file by the C to identify a comment line - but then the statements should start in column 7, continuation identifyers in 6, labels in 1 to 5.

But this one got mixed up. Dani 8586 should proceed to post this file by first setting the (code) and (/code) tags - with [ ] instead of ( ) - and then paste the original file in between them. Otherwise the html of the page will cut off the leading blanks of each line.

Norbert

BTW, mikrom, how did you make the codetags show as typed in your post ??

Norbert

The optimist believes we live in the best of all possible worlds - the pessimist fears this might be true.

dani8586 · Oct 16, 2012

Hi, now should be ok

Code:

	program test
        implicit real*8 (a-h,o-z)
	external fxn


	COMMON/BVEG1/NCALL,ITMX,NPRN,NDEV,XL(100),XU(100),ACC
        COMMON/BVEG2/IT,NDO,SI,SWGT,SCHI,XI(50,100)
	COMMON/BVEG3/ALPH,NDMX,MDS
        COMMON/RNDM/ISEED

	ndim=1
	ncall=10
	itmx=1
	nprn=-1
	acc=1.0d-12
	xl(1)=0.0
	xu(1)=1.0

c       iseed=0
c	call srand(iseed)
	call VEGAS(NDIM,FXN,AVGI,SD,CHI2A)
	end



	function fxn(x,wgt)
        implicit real*8 (a-h,o-z)
	double precision x(100)

	fxn=x(1)**2
	return
	end

C  START VEGAS SECTION
C
      BLOCK DATA
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C   MAKES DEFAULT PARAMETER ASSIGNMENTS FOR VEGAS
      COMMON/BVEG1/NCALL,ITMX,NPRN,NDEV,XL(100),XU(100),ACC
      COMMON/BVEG2/IT,NDO,SI,SWGT,SCHI,XI(50,100)
      COMMON/BVEG3/ALPH,NDMX,MDS
      DATA NCALL/1000/,ITMX/10/,NPRN/0/,ACC/1.D-2/,
     1     XL/100*0.d0/,XU/100*1.d0/,
     3     ALPH/1.5d0/,NDMX/50/,MDS/1/,NDEV/6/,
     4     NDO/1/,XI/5000*1.d0/,IT/0/,SI,SWGT,SCHI/3*0.d0/
      END


      SUBROUTINE VEGAS(NDIM,FXN,AVGI,SD,CHI2A)
C
C   SUBROUTINE PERFORMS NDIM-DIMENSIONAL MONTE CARLO INTEG'N
C      - BY G.P. LEPAGE    SEPT 1976/(REV)AUG 1979
C      - ALGORITHM DESCRIBED IN J COMP PHYS 27,192(1978)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      COMMON/BVEG1/NCALL,ITMX,NPRN,NDEV,XL(100),XU(100),ACC
      COMMON/BVEG2/IT,NDO,SI,SWGT,SCHI,XI(50,100)
      COMMON/BVEG3/ALPH,NDMX,MDS
      COMMON/BVEG4/CALLS,TI,TSI
      COMMON/RESLOCAL/RESL(20),STANDDEVL(20)
      DIMENSION D(50,100),DI(50,100),XIN(50),R(50),DX(100),IA(100),
     1          KG(100),DT(100),X(100)
      DIMENSION RAND(100)
      DATA ONE/1.d0/
C		
      NDO=1
      DO 1 J=1,NDIM
 1    XI(1,J)=ONE
C
      ENTRY VEGAS1(NDIM,FXN,AVGI,SD,CHI2A)
C          INITIALIZES CUMMULATIVE VARIABLES, BUT NOT GRID
      IT=0
      SI=0.d0
      SWGT=SI
      SCHI=SI
C
      ENTRY VEGAS2(NDIM,FXN,AVGI,SD,CHI2A)
C         - NO INITIALIZATION
      ND=NDMX
      NG=1
      IF(MDS.EQ.0) GO TO 2
      NG=(NCALL/2.d0)**(1.d0/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
      DX(J)=XU(J)-XL(J)
 3    XJAC=XJAC*DX(J)
C
C   REBIN, PRESERVING BIN DENSITY
      IF(ND.EQ.NDO) GO TO 8
      RC=NDO/XND
      DO 7 J=1,NDIM
      K=0
      XN=0.d0
      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     IF(NPRN.GE.0) WRITE(NDEV,200) NDIM,CALLS,IT,ITMX,ACC,NPRN,
     1                    ALPH,MDS,ND,(XL(J),XU(J),J=1,NDIM)
C
      ENTRY VEGAS3(NDIM,FXN,AVGI,SD,CHI2A)
C         - MAIN INTEGRATION LOOP
 9    IT=IT+1
      TI=0.d0
      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.d0
      F2B=FB
      K=0
 12   K=K+1
      CALL RANDA(NDIM,RAND)
      WGT=XJAC
      DO 15 J=1,NDIM
      XN=(KG(J)-RAND(J))*DXG+ONE
      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)
 15   WGT=WGT*XO*XND
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=SQRT(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   COMPUTE FINAL RESULTS FOR THIS ITERATION
      TSI=TSI*DV2G
      TI2=TI*TI
CCCCCCCCCCCCCCCCC MODIFICA 1: PERMETTE INTEGRALE NULLO  CCCCCCCCCCCCCCCC`
C
      IF(TSI.EQ.0.) THEN
        WGT=0.d0
        SI=0.d0
        SWGT=0.d0
        SCHI=0.d0
        AVGI=0.d0
        CHI2A=0.d0
        SD=0.d0
        IF(NPRN.LT.0) THEN
          CONTINUE
        ELSE
          TSI=SQRT(TSI) 
          WRITE(NDEV,201) IT,TI,TSI,AVGI,SD,CHI2A       
          IF(NPRN.EQ.0) THEN
            CONTINUE
          ELSE
            DO J=1,NDIM
              WRITE(NDEV,202) J,(XI(I,J),DI(I,J),I=1,ND,NPRN)       
            END DO
          END IF
        END IF
        RETURN
      END IF
C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC`
C
C PROTECT AGAINST ACCIDENTAL SINGLE PRECISION NEGATIVE VALUES
      XXX=TSI/TI**2
      IF(XXX.GT.-1.D-4) F2B=ABS(F2B)

      WGT=ONE/TSI
      SI=SI+TI*WGT
      SWGT=SWGT+WGT
      SCHI=SCHI+TI2*WGT
      AVGI=SI/SWGT
      CHI2A=(SCHI-SI*AVGI)/(IT-.9999d0)
      SD=SQRT(ONE/SWGT)
C
      IF(NPRN.LT.0) GO TO 21
      TSI=SQRT(TSI)
      WRITE(NDEV,201) IT,TI,TSI,AVGI,SD,CHI2A
      RESL(IT)=TI
      STANDDEVL(IT)=TSI
      IF(NPRN.EQ.0) GO TO 21
      DO 20 J=1,NDIM
 20   WRITE(NDEV,202) J,(XI(I,J),DI(I,J),I=1,ND,NPRN)
C
C   REFINE GRID
 21   DO 23 J=1,NDIM
      XO=D(1,J)
      XN=D(2,J)
      D(1,J)=(XO+XN)/2.d0
      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.d0
 22   DT(J)=DT(J)+D(I,J)
      D(ND,J)=(XN+XO)/2.d0
 23   DT(J)=DT(J)+D(ND,J)
C
      DO 28 J=1,NDIM
      RC=0.d0
      DO 24 I=1,ND
      R(I)=0.d0
      IF(D(I,J).LE.0.d0) GO TO 24
      XO=DT(J)/D(I,J)
      R(I)=((XO-ONE)/XO/LOG(XO))**ALPH
 24   RC=RC+R(I)
      RC=RC/XND
      K=0
      XN=0.d0
      DR=XN
      I=K
 25   K=K+1
      DR=DR+R(K)
      XO=XN
      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)
      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*ABS(AVGI).LT.SD) GO TO 9
 200  FORMAT(/35H INPUT PARAMETERS FOR VEGAS:  NDIM=,I3,8H  NCALL=,F8.0
     1  /28X,5H  IT=,I5,7H  ITMX=,I5/28X,6H  ACC=,G9.3
     2  /28X,7H  NPRN=,I3,7H  ALPH=,F5.2/28X,6H  MDS=,I3,6H   ND=,I4
     3  /28X,10H  (XL,XU)=,(T40,2H( ,G12.6,3H , ,G12.6,2H )))
 201  FORMAT(///21H INTEGRATION BY VEGAS//14H ITERATION NO.,I3,
     1  14H:   INTEGRAL =,G14.8/21X,10HSTD DEV  =,G10.4/
     2  34H ACCUMULATED RESULTS:   INTEGRAL =,G14.8/
     3  24X,10HSTD DEV  =,G10.4/24X,17HCHI**2 PER IT'N =,G10.4)
 202  FORMAT(/15H DATA FOR AXIS ,I2/25H    X       DELTA I       ,
     1  24H   X       DELTA I       ,18H   X       DELTA I
     2  /(1H ,F7.6,1X,G11.4,5X,F7.6,1X,G11.4,5X,F7.6,1X,G11.4))
      RETURN
      END

C***************************************************************
C LOAD VEGAS DATA IF DESIRED
      SUBROUTINE LOAD_VEGAS(NDIM,NAME)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      COMMON/BVEG2/IT,NDO,SI,SWGT,SCHI,XI(50,100)
      CHARACTER*30 NAME
C      OPEN(UNIT=13,FILE=NAME,STATUS='OLD',SHARED)
      READ(13,210) IT,NDO,SI,SWGT,SCHI
      DO 190 J=1,NDIM
 190  READ(13,*)  (XI(I,J),I=1,50)
 210  FORMAT(2I8,3Z16)
      CLOSE(UNIT=13)
      RETURN
      END

C*******************************************************************
C STORE VEGAS DATA FOR POSSIBLE LATER USE
      SUBROUTINE STORE_VEGAS(NDIM,NAME)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      COMMON/BVEG2/IT,NDO,SI,SWGT,SCHI,XI(50,100)
      CHARACTER*30 NAME
      OPEN(UNIT=12,FILE=NAME,STATUS='NEW')
      WRITE(12,210) IT,NDO,SI,SWGT,SCHI
      DO 190 J=1,NDIM
 190  WRITE(12,*)  (XI(I,J),I=1,50)
 210  FORMAT(2I8,3Z16)
      CLOSE(UNIT=12)
      RETURN
      END

      SUBROUTINE RANDA(N,RAND)
C
C   SUBROUTINE GENERATES UNIFORMLY DISTRIBUTED RANDOM NO'S X(I),I=1,N
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      COMMON/RNDM/ISEED
      DIMENSION RAND(N)
      DO 1 I=1,N
1     RAND(I)=RAN(iseed)
      RETURN
      END

C--------- END THE VEGAS SECTION -------------------------------------

mikrom · Oct 17, 2012

The problem seems to be witjh the function RAN() - see

http://gcc.gnu.org/onlinedocs/gfortran/RAN.html

RAN(ISEED) doesn't generate for the ISEED > 0 random numbers, but it generates the same numbers:

Code:

$ test_rand
 * iseed =            0
 Output of RAN():
  7.62939453E-06
  0.13153768    
  0.75560522    
  0.45865011    
  0.53276706    
  0.21895909    
  4.70445156E-02
  0.67886448    
  0.67929626    
  0.93469286    
 Output of SRAND() and RAND():
  0.24257827    
  1.34694576E-02
  0.38313866    
  0.41465259    
  6.77688122E-02
  0.99312687    
  0.48430800    
  0.76533771    
  3.18336487E-02
  3.09352875E-02
 * iseed =       617647
 Output of RAN():
[highlight #FCE94F]  0.83393335    
  0.83393335    
  0.83393335    
  0.83393335    
  0.83393335    
  0.83393335    
  0.83393335    
  0.83393335    
  0.83393335    
  0.83393335    [/highlight]    
 Output of SRAND() and RAND():
  0.83393335    
  0.92023325    
  0.36043763    
  0.87803197    
  8.58857632E-02
  0.48455787    
  0.96607089    
  0.75577950    
  0.38717556    
  0.26309896

The Fortran code i used is here
test_rand.f

Code:

      [COLOR=#a020f0]program[/color] test_rand
        [COLOR=#2e8b57][b]integer[/b][/color] iseed, N 
         
        N [COLOR=#804040][b]=[/b][/color] [COLOR=#ff00ff]10[/color]

        iseed [COLOR=#804040][b]=[/b][/color] [COLOR=#ff00ff]0[/color]
        [COLOR=#804040][b]write[/b][/color]([COLOR=#804040][b]*[/b][/color],[COLOR=#804040][b]*[/b][/color]) [COLOR=#ff00ff]'* iseed = '[/color], iseed
        [COLOR=#008080]call[/color] test_random(iseed, N)

        iseed [COLOR=#804040][b]=[/b][/color] [COLOR=#ff00ff]617647[/color]
        [COLOR=#804040][b]write[/b][/color]([COLOR=#804040][b]*[/b][/color],[COLOR=#804040][b]*[/b][/color]) [COLOR=#ff00ff]'* iseed = '[/color], iseed
        [COLOR=#008080]call[/color] test_random(iseed, N)
      [COLOR=#a020f0]end program[/color]

      [COLOR=#a020f0]subroutine[/color] test_random(iseed, N)
        [COLOR=#2e8b57][b]integer[/b][/color] iseed, N
[COLOR=#0000ff]c       function RAN()[/color]
        [COLOR=#804040][b]write[/b][/color]([COLOR=#804040][b]*[/b][/color],[COLOR=#804040][b]*[/b][/color]) [COLOR=#ff00ff]'Output of RAN():'[/color]
        [COLOR=#804040][b]do[/b][/color] i[COLOR=#804040][b]=[/b][/color][COLOR=#ff00ff]1[/color], N
          [COLOR=#804040][b]print[/b][/color] [COLOR=#804040][b]*[/b][/color], ran(iseed)
        [COLOR=#804040][b]end do[/b][/color]
[COLOR=#0000ff]c       functions SRAND() and RAND()[/color]
        [COLOR=#804040][b]write[/b][/color]([COLOR=#804040][b]*[/b][/color],[COLOR=#804040][b]*[/b][/color]) [COLOR=#ff00ff]'Output of SRAND() and RAND():'[/color]
        [COLOR=#008080]call[/color] srand(iseed)
        [COLOR=#804040][b]do[/b][/color] i[COLOR=#804040][b]=[/b][/color][COLOR=#ff00ff]1[/color], N
          [COLOR=#804040][b]print[/b][/color] [COLOR=#804040][b]*[/b][/color], rand()
        [COLOR=#804040][b]end do[/b][/color]
      [COLOR=#a020f0]end subroutine[/color]

mikrom · Oct 17, 2012

To make some observation, I modified the source a little bit.
In the main program I added some prints

Code:

      write(*,*) '*** DBG ISEED =', ISEED
      write(*,*) '*** DBG NCALL = ', NCALL
      call VEGAS(NDIM,FXN,AVGI,SD,CHI2A)
      write(*,*) '*** DBG Results: ', NDIM,AVGI,SD,CHI2A

and in the subroutine RANDA too

Code:

      SUBROUTINE RANDA(N,RAND)
C
C     SUBROUTINE GENERATES UNIFORMLY DISTRIBUTED RANDOM NO'S X(I),I=1,N
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      COMMON/RNDM/ISEED
      DIMENSION RAND(N)
      write(*,*) '*DBG Enter SUBROUTINE RANDA'
      write(*,*) '*DBG N = ', N
      write(*,*) '*DBG ISEED = ', ISEED
      DO 1 I=1,N
      RAND(I)=RAN(ISEED)
 1    write(*,*) '*DBG RAND(', I ,') = ', RAND(I)
      RETURN
      END

Now when we set in the main program

Code:

NCALL = 5 
ISEED = 0

then we get this

Code:

$ dani8586
 *** DBG ISEED =           0
 *** DBG NCALL =            5
 *DBG Enter SUBROUTINE RANDA
 *DBG N =            1
 *DBG ISEED =            0
 *DBG RAND(           1 ) =   7.62939453125000000E-006
 *DBG Enter SUBROUTINE RANDA
 *DBG N =            1
 *DBG ISEED =            0
 *DBG RAND(           1 ) =   0.13153767585754395     
 *DBG Enter SUBROUTINE RANDA
 *DBG N =            1
 *DBG ISEED =            0
 *DBG RAND(           1 ) =   0.75560522079467773     
 *DBG Enter SUBROUTINE RANDA
 *DBG N =            1
 *DBG ISEED =            0
 *DBG RAND(           1 ) =   0.45865011215209961     
 *** DBG Results:            1  0.35490558706289216       5.39359000355652540E-002  3.33066907387583636E-012

RAND(1) got by every call of rubroutine RANDA() other value and the end result is not zero.

Now when we set in the main program

Code:

ISEED = 617647

then RAND(1) obtains by every call the same value and the end result is zero - see:

Code:

$ dani8586
 *** DBG ISEED =      617647
 *** DBG NCALL =            5
 *DBG Enter SUBROUTINE RANDA
 *DBG N =            1
 *DBG ISEED =       617647
 *DBG RAND(           1 ) =   0.83393335342407227     
 *DBG Enter SUBROUTINE RANDA
 *DBG N =            1
 *DBG ISEED =       617647
 *DBG RAND(           1 ) =   0.83393335342407227     
 *DBG Enter SUBROUTINE RANDA
 *DBG N =            1
 *DBG ISEED =       617647
 *DBG RAND(           1 ) =   0.83393335342407227     
 *DBG Enter SUBROUTINE RANDA
 *DBG N =            1
 *DBG ISEED =       617647
 *DBG RAND(           1 ) =   0.83393335342407227     
 *** DBG Results:            1   0.0000000000000000        0.0000000000000000        0.0000000000000000

So it seems that the computation works only with ISEED = 0

dani8586 · Oct 18, 2012

Thanks!
Actually I found a strange thing:
If I compile with gfortran actually looks that it work only with iseed=0; but if before call vegas i put a call srand(iseed) and in the vegas program i take of the argument of Ran(iseed), so if I write Ran() looks like that I can change seed, that is what I'm looking for... Don't know why however!

mikrom · Oct 19, 2012

Don't know why it's so, but as written here

http://gcc.gnu.org/onlinedocs/gfortran/RAN.html

the function RAN() seems to be in gfortran only for compatibility with older compiler.

The better way seems to be to use SRAND() and RAND(), like in example given here

http://gcc.gnu.org/onlinedocs/gfortran/RAND.html#RAND

or RANDOM_SEED and RANDOM_NUMBER like in example here

http://gcc.gnu.org/onlinedocs/gfortran/RANDOM_005fNUMBER.html#RANDOM_005fNUMBER

dani8586 · Oct 19, 2012

thanks !

Tek-Tips is the largest IT community on the Internet today!

Members share and learn making Tek-Tips Forums the best source of peer-reviewed technical information on the Internet!

Vegas integration problem

dani8586

Programmer

mikrom

Programmer

dani8586

Programmer

mikrom

Programmer

gummibaer

Programmer

dani8586

Programmer

mikrom

Programmer

mikrom

Programmer

dani8586

Programmer

mikrom

Programmer

dani8586

Programmer

Similar threads

Part and Inventory Search

Sponsor