Sorry. Please see the code below.
Just added some _i8 so gfortran is happier.
gfortran -O3 allinone2.f90 -o test
then
./test
should work.
It is RNG from my PhD advisor, probably box muller. It has been used in many nuclear quantum Monte Carlo code. If you have better ones please advise.
Thank you so much!
module random2
implicit none
integer, private, parameter :: i4=selected_int_kind(9)
integer, private, parameter :: i8=selected_int_kind(15)
integer, private, parameter :: r8=selected_real_kind(15,9)
contains
subroutine ran1(rn,irn)
! multiplicative congruential with additive constant
integer(kind=i8), parameter :: mask24 = ishft(1_8,24)-1
integer(kind=i8), parameter :: mask48 = ishft(1_8,48_8)-1_8
real(kind=r8), parameter :: twom48=2.0d0**(-48)
integer(kind=i8), parameter :: mult1 = 44485709377909_8
integer(kind=i8), parameter :: m11 = iand(mult1,mask24)
integer(kind=i8), parameter :: m12 = iand(ishft(mult1,-24),mask24)
integer(kind=i8), parameter :: iadd1 = 96309754297_8
integer(kind=i8) :: irn
real(kind=r8) :: rn
integer(kind=i8) :: is1,is2
is2 = iand(ishft(irn,-24),mask24)
is1 = iand(irn,mask24)
irn = iand(ishft(iand(is1*m12+is2*m11,mask24),24)+is1*m11+iadd1,mask48)
rn = ior(irn,1_8)*twom48
return
end subroutine ran1
subroutine ran2(rn,irn)
! multiplicative congruential with additive constant
integer(kind=i8), parameter :: mask24 = ishft(1_8,24)-1
integer(kind=i8), parameter :: mask48 = ishft(1_8,48)-1
real(kind=r8), parameter :: twom48=2.0d0**(-48)
integer(kind=i8), parameter :: mult2 = 34522712143931_8
integer(kind=i8), parameter :: m21 = iand(mult2,mask24)
integer(kind=i8), parameter :: m22 = iand(ishft(mult2,-24),mask24)
integer(kind=i8), parameter :: iadd2 = 55789347517_8
integer(kind=i8) :: irn
real(kind=r8) :: rn
integer(kind=i8) :: is1,is2
is2 = iand(ishft(irn,-24),mask24)
is1 = iand(irn,mask24)
irn = iand(ishft(iand(is1*m22+is2*m21,mask24),24)+is1*m21+iadd2,mask48)
rn = ior(irn,1_8)*twom48
return
end subroutine ran2
end module random2
module random
implicit none
integer, private, parameter :: i4=selected_int_kind(9)
integer, private, parameter :: i8=selected_int_kind(15)
integer, private, parameter :: r8=selected_real_kind(15,9)
integer(kind=i8), private, parameter :: mask48 = ishft(1_8,48)-1
integer(kind=i8), private, save :: irn = 1_8,irnsave
contains
function randn(n)
! return an array of random variates (0,1)
use random2
integer(kind=i8) :: n
real(kind=r8), dimension(n) :: randn
integer(kind=i8) :: i
do i=1,n
call ran1(randn(i),irn)
enddo
return
end function randn
function gaussian(n)
! return an array of gaussian random variates with variance 1
integer(kind=i8) :: n
real(kind=r8), dimension(n) :: gaussian
real(kind=r8), dimension(2*((n+1)/2)) :: rn
real(kind=r8) :: rn1,rn2,arg,x1,x2,x3
real(kind=r8), parameter :: pi=4.0d0*atan(1.0_r8)
integer(kind=i8) :: i
rn=randn(size(rn,kind=i8))
do i=1,n/2
rn1=rn(2*i-1)
rn2=rn(2*i)
arg=2.0_r8*pi*rn1
x1=sin(arg)
x2=cos(arg)
x3=sqrt(-2.0_r8*log(rn2))
gaussian(2*i-1)=x1*x3
gaussian(2*i)=x2*x3
enddo
if (mod(n,2).ne.0) then
rn1=rn(n)
rn2=rn(n+1)
arg=2.0_r8*pi*rn1
x2=cos(arg)
x3=sqrt(-2.0_r8*log(rn2))
gaussian(n)=x2*x3
endif
return
end function gaussian
subroutine setrn(irnin)
! set the seed
integer(kind=i8) :: irnin
integer(kind=i8) :: n
integer(kind=i8), allocatable :: iseed(:)
irn=iand(irnin,mask48)
call savern(irn)
!call RANDOM_SEED(SIZE=n)
allocate(iseed(n))
iseed=int(irnin)
!call RANDOM_SEED(PUT=iseed) ! set the internal ran function's seed.
return
end subroutine setrn
subroutine savern(irnin)
! save the seed
integer(kind=i8) :: irnin
irnsave=irnin
return
end subroutine savern
subroutine showirn(irnout)
! show the seed
integer(kind=i8) :: irnout
irnout=irnsave
return
end subroutine showirn
function randomn(n) ! the default ramdom_number subroutine.
! return an array of random variates (0,1)
integer(kind=i8) :: n
real(kind=r8), dimension(n) :: randomn
call RANDOM_NUMBER(randomn)
return
end function randomn
subroutine random1(x) ! generate random number x using the default ramdom_number subroutine.
real(kind=r8) :: x
call RANDOM_NUMBER(x)
return
end subroutine random1
end module random
module constants
implicit none
integer, public, parameter :: i4=selected_int_kind(9)
integer, public, parameter :: i8=selected_int_kind(15)
integer, public, parameter :: r8=selected_real_kind(15,9)
real(kind=r8), public, parameter :: zero=0.0_r8,one=1.0_r8,two=2.0_r8,three=3.0_r8,four=4.0_r8 &
,five=5.0_r8,six=6.0_r8,seven=7.0_r8,eight=8.0_r8,nine=9.0_r8 &
,ten=10.0_r8,tenth=.1_r8,half=.5_r8,third=1.0_r8/3.0_r8,sixth=1.0_r8/6.0_r8 &
,pi=4.0_r8*atan(1.0_r8) &
,normpdf_factor=1.0_r8/sqrt(8.0_r8*atan(1.0_r8)) ! 1/sqrt(2 pi)
complex(kind=r8), public, parameter :: czero=(0.0_r8,0.0_r8),ci=(0.0_r8,1.0_r8),cone=(1.0_r8,0.0_r8)
end module constants
program main
use random
use random2
use constants
implicit none
integer(kind=i8) :: seed
real (kind=r8) :: x1(1),x100(100)
! use the RNG illustration
seed = 123456789
call setrn(int(seed,8))
! generate one uniform randum number
x1=randn(1_8)
write (6,*) 'the one random number is', x1(1)
! generate 100 uniform random number array
x100=randn(100_i8)
write (6,*) 'the 100 randum numbers are ', x100
! generate one gaussian randum number
x1=gaussian(1_i8)
write (6,*) ' the one Gaussian random number is', x1(1)
! generate 100 gaussian random number array
x100=gaussian(100_i8)
write (6,*) 'the 100 Gaussian randum numbers are ', x100
stop ('Program end normally.')
end program main