randoms.f90

Go to the documentation of this file.

 
! Code converted using TO_F90 by Alan Miller
! Date: 2012-03-16  Time: 11:09:33
 
 
MODULE mprand
    USE mpdef
    IMPLICIT NONE
 
    INTEGER(mpi), PARAMETER :: nb=511
    INTEGER(mpi), PARAMETER :: ia=16807
    INTEGER(mpi), PARAMETER :: im=2147483647
    INTEGER(mpi), PARAMETER :: iq=127773
    INTEGER(mpi), PARAMETER :: ir=2836
    REAL(mps), PARAMETER :: aeps=1.0e-10
    REAL(mps), PARAMETER :: scalin=4.6566125e-10
    INTEGER(mpi) :: mbuff(0:nb+1),ian,ic,iboost
    INTEGER(mpi) :: iseed = 4711
    INTEGER(mpi) :: istart = 0
    INTEGER(mpi) :: iwarm = 10
 
END MODULE mprand
 
 
SUBROUTINE gbrshi(n,a)
    USE mprand
 
    IMPLICIT NONE
    INTEGER(mpi) :: i
    INTEGER(mpi) :: idum
    INTEGER(mpi) :: it
    INTEGER(mpi) :: j
    INTEGER(mpi) :: k
 
    INTEGER(mpi), INTENT(IN)                      :: n
    REAL(mps), INTENT(OUT)                        :: a(*)
 
    IF(istart /= 0) GO TO 20
    WRITE(*,*) ' Automatic GBRSHI initialization using:'
    !     initialize buffer
    idum=iseed+9876543          ! prevent damage, if iseed=0
    WRITE(*,*) '           ISEED=',iseed,'   IWARM=',iwarm
    DO j=0,nb+1                 ! fill buffer
        k=idum/iq                  ! minimal standard generator
        idum=ia*(idum-k*iq)-ir*k   !    with Schrages method
        IF(idum < 0) idum=idum+im !
        mbuff(j)=ishft(idum,1)     ! fill in leading bit
    END DO
    ian=iand(ian,nb)            ! mask angle
    ic=1                        ! set pointer
    iboost=0
    DO j=1,iwarm*nb             ! warm up a few times
        it=mbuff(ian)              ! hit ball angle
        mbuff(ian)=ieor(ior(ishft(it,17),ishft(it,-15)),ic)  ! new spin
        ic=it                      ! replace red spin
        ian=iand(it+iboost,nb)     ! boost and mask angle
        iboost=iboost+1            ! increment boost
    END DO
    istart=1                    ! set done-flag
    !     generate array of r.n.
 20 CONTINUE
    DO i=1,n
        it=mbuff(ian)              ! hit ball angle
        mbuff(ian)=ieor(ior(ishft(it,17),ishft(it,-15)),ic)  ! new spin
        ic=it                      ! replace red spin
        ian=iand(it+iboost,nb)     ! boost and mask angle
        a(i)=real(ishft(it,-1),mps)*scalin+aeps ! avoid zero output
        iboost=iboost+1            ! increment boost
    END DO
    iboost=iand(iboost,nb)
    RETURN
END SUBROUTINE gbrshi
 
SUBROUTINE gbrvin(jseed,jwarm)
    USE mprand
 
    IMPLICIT NONE
    INTEGER(mpi) :: idum
    INTEGER(mpi) :: it
    INTEGER(mpi) :: j
    INTEGER(mpi) :: k
 
 
    INTEGER(mpi), INTENT(IN) :: jseed
    INTEGER(mpi), INTENT(IN) :: jwarm
 
    IF(istart == 0) THEN
        WRITE(*,*) ' Gbrshi initialization by GBRVIN-call using:'
        iseed=jseed              ! copy seed and
        iwarm=jwarm              ! warm-up parameter
        istart=-1                ! start flag
        idum=iseed+9876543          ! prevent damage, if iseed=0
        WRITE(*,*) '           ISEED=',iseed,'   IWARM=',iwarm
        DO j=0,nb+1                 ! fill buffer
            k=idum/iq                  ! minimal standard generator
            idum=ia*(idum-k*iq)-ir*k   !    with Schrages method
            IF(idum < 0) idum=idum+im !
            mbuff(j)=ishft(idum,1)     ! fill in leading bit
        END DO
        ian=iand(ian,nb)            ! mask angle
        ic=1                        ! set pointer
        iboost=0
        DO j=1,iwarm*nb             ! warm up a few times
            it=mbuff(ian)              ! hit ball angle
            mbuff(ian)=ieor(ior(ishft(it,17),ishft(it,-15)),ic)  ! new spin
            ic=it                      ! replace red spin
            ian=iand(it+iboost,nb)     ! boost and mask angle
            iboost=iboost+1            ! increment boost
        END DO
    END IF
END SUBROUTINE gbrvin
 
SUBROUTINE gbrtim
    USE mpdef
 
    IMPLICIT NONE
    INTEGER(mpi) :: jseed
    REAL(mps) :: time
 
    LOGICAL :: done
    DATA    done/.false./
    IF(done) RETURN
    jseed=time()
    WRITE(*,*) ' Gbrshi initialialization using Time()'
    CALL gbrvin(jseed,10)
    done=.true.
END SUBROUTINE gbrtim
 
 
REAL(mps) FUNCTION uran()     ! U(0,1)
    USE mpdef
 
    IMPLICIT NONE
    INTEGER(mpi) :: indx
    INTEGER(mpi) :: ndim
 
    parameter(ndim=100)
    REAL(mps) :: buffer(ndim)
    DATA indx/ndim/
    SAVE indx,buffer
    indx=mod(indx,ndim)+1
    IF(indx == 1) CALL gbrshi(ndim,buffer)
    uran=buffer(indx)
END FUNCTION uran
 
 
REAL(mps) FUNCTION gran()     ! N(0,1)
    USE mpdef
 
    IMPLICIT NONE
    REAL(mps) :: al
    REAL(mps) :: cs
    INTEGER(mpi) :: indx
    INTEGER(mpi) :: kn
    INTEGER(mpi) :: ndim
    REAL(mps) :: radsq
    REAL(mps) :: rn1
    REAL(mps) :: rn2
    REAL(mps) :: sn
 
    parameter(ndim=100)
    REAL(mps) :: buffer(ndim)
    DATA indx/ndim/,kn/1/
    SAVE indx,buffer,kn,cs,al
    !     ...
    IF(kn <= 1) THEN
        !        two U(-1,+1) random numbers
10      indx=mod(indx,ndim)+2
        IF(indx == 2) CALL gbrshi(ndim,buffer)
        rn1=buffer(indx-1)-1.0+buffer(indx-1)
        rn2=buffer(indx  )-1.0+buffer(indx)
        radsq=rn1*rn1+rn2*rn2
        IF(radsq > 1.0) GO TO 10 ! test point inside circle?
        !        sine and cosine for random phi
        sn=rn1/sqrt(radsq)
        cs=rn2/sqrt(radsq)
        !        transform to gaussians
        al=sqrt(-2.0*log(radsq))
        kn =2
        gran=sn*al
    ELSE
        kn =1
        gran=cs*al
    END IF
END FUNCTION gran