Dismiss Notice
Join Physics Forums Today!
The friendliest, high quality science and math community on the planet! Everyone who loves science is here!

[Fortran] Problem with Vegas integration

  1. Oct 15, 2012 #1
    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 --------------------
     
  2. jcsd
Know someone interested in this topic? Share this thread via Reddit, Google+, Twitter, or Facebook

Can you offer guidance or do you also need help?
Draft saved Draft deleted



Similar Discussions: [Fortran] Problem with Vegas integration
  1. Problems In Fortran (Replies: 2)

  2. Fortran problem (Replies: 1)

  3. FORTRAN READ problems (Replies: 5)

Loading...