!$$$ Module Documentation Block
!
! Module: mersenne_twister Modern random number generator
! Prgmmr: Iredell Org: W/NX23 date: 2005-06-14
!
! Abstract: This module calculates random numbers using the Mersenne twister.
! (It has been adapted to a Fortran 90 module from open source software.
! The comments from the original software are given below in the remarks.)
! The Mersenne twister (aka MT19937) is a state-of-the-art random number
! generator based on Mersenne primes and originally developed in 1997 by
! Matsumoto and Nishimura. It has a period before repeating of 2^19937-1,
! which certainly should be good enough for geophysical purposes. :-)
! Considering the algorithm's robustness, it runs fairly speedily.
! (Some timing statistics are given below in the remarks.)
! This adaptation uses the standard Fortran 90 random number interface,
! which can generate an arbitrary number of random numbers at one time.
! The random numbers generated are uniformly distributed between 0 and 1.
! The module also can generate random numbers from a Gaussian distribution
! with mean 0 and standard deviation 1, using a Numerical Recipes algorithm.
! The module also can generate uniformly random integer indices.
! There are also thread-safe versions of the generators in this adaptation,
! necessitating the passing of generator states which must be kept private.
!
! Program History Log:
! 2005-06-14 Mark Iredell
!
! Usage:
! The module can be compiled with 4-byte reals or with 8-byte reals, but
! 4-byte integers are required. The module should be endian-independent.
! The Fortran 90 interfaces random_seed and random_number are overloaded
! and can be used as in the standard by adding the appropriate use statement
! use mersenne_twister
! In the below use cases, harvest is a real array of arbitrary size,
! and iharvest is an integer array of arbitrary size.
! To generate uniformly distributed random numbers between 0 and 1,
! call random_number(harvest)
! To generate Gaussian distributed random numbers with 0 mean and 1 sigma,
! call random_gauss(harvest)
! To generate uniformly distributed random integer indices between 0 and n,
! call random_index(n,iharvest)
! In standard "saved" mode, the random number generator can be used without
! setting a seed. But to set a seed, only 1 non-zero integer is required, e.g.
! call random_setseed(4357) ! set default seed
! The full generator state can be set via the standard interface random_seed,
! but it is recommended to use this method only to restore saved states, e.g.
! call random_seed(size=lsave) ! get size of generator state seed array
! allocate isave(lsave) ! allocate seed array
! call random_seed(get=isave) ! fill seed array (then maybe save to disk)
! call random_seed(put=isave) ! restore state (after read from disk maybe)
! Locally kept generator states can also be saved in a seed array, e.g.
! type(random_stat):: stat
! call random_seed(get=isave,stat=stat) ! fill seed array
! call random_seed(put=isave,stat=stat) ! restore state
! To generate random numbers in a threaded region, the "thread-safe" mode
! must be used where generator states of type random_state are passed, e.g.
! type(random_stat):: stat(8)
! do i=1,8 ! threadable loop
! call random_setseed(7171*i,stat(i)) ! thread-safe call
! enddo
! do i=1,8 ! threadable loop
! call random_number(harvest,stat(i)) ! thread-safe call
! enddo
! do i=1,8 ! threadable loop
! call random_gauss(harvest,stat(i)) ! thread-safe call
! enddo
! do i=1,8 ! threadable loop
! call random_index(n,iharvest,stat(i))! thread-safe call
! enddo
! There is also a relatively inefficient "interactive" mode available, where
! setting seeds and generating random numbers are done in the same call.
! There is also a functional mode available, returning one value at a time.
!
! Public Defined Types:
! random_stat Generator state (private contents)
!
! Public Subprograms:
! random_seed determine size or put or get state
! size optional integer output size of seed array
! put optional integer(:) input seed array
! get optional integer(:) output seed array
! stat optional type(random_stat) (thread-safe mode)
! random_setseed set seed (thread-safe mode)
! inseed integer seed input
! stat type(random_stat) output
! random_setseed set seed (saved mode)
! inseed integer seed input
! random_number get mersenne twister random numbers (thread-safe mode)
! harvest real(:) numbers output
! stat type(random_stat) input
! random_number get mersenne twister random numbers (saved mode)
! harvest real(:) numbers output
! random_number get mersenne twister random numbers (interactive mode)
! harvest real(:) numbers output
! inseed integer seed input
! random_number_f get mersenne twister random number (functional mode)
! harvest real number output
! random_gauss get gaussian random numbers (thread-safe mode)
! harvest real(:) numbers output
! stat type(random_stat) input
! random_gauss get gaussian random numbers (saved mode)
! harvest real(:) numbers output
! random_gauss get gaussian random numbers (interactive mode)
! harvest real(:) numbers output
! inseed integer seed input
! random_gauss_f get gaussian random number (functional mode)
! harvest real number output
! random_index get random indices (thread-safe mode)
! imax integer maximum index input
! iharvest integer(:) numbers output
! stat type(random_stat) input
! random_index get random indices (saved mode)
! imax integer maximum index input
! iharvest integer(:) numbers output
! random_index get random indices (interactive mode)
! imax integer maximum index input
! iharvest integer(:) numbers output
! inseed integer seed input
! random_index_f get random index (functional mode)
! imax integer maximum index input
! iharvest integer number output
!
! Remarks:
! (1) Here are the comments in the original open source code:
! A C-program for MT19937: Real number version
! genrand() generates one pseudorandom real number (double)
! which is uniformly distributed on [0,1]-interval, for each
! call. sgenrand(seed) set initial values to the working area
! of 624 words. Before genrand(), sgenrand(seed) must be
! called once. (seed is any 32-bit integer except for 0).
! Integer generator is obtained by modifying two lines.
! Coded by Takuji Nishimura, considering the suggestions by
! Topher Cooper and Marc Rieffel in July-Aug. 1997.
! This library is free software; you can redistribute it and/or
! modify it under the terms of the GNU Library General Public
! License as published by the Free Software Foundation; either
! version 2 of the License, or (at your option) any later
! version.
! This library is distributed in the hope that it will be useful,
! but WITHOUT ANY WARRANTY; without even the implied warranty of
! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
! See the GNU Library General Public License for more details.
! You should have received a copy of the GNU Library General
! Public License along with this library; if not, write to the
! Free Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA
! 02111-1307 USA
! Copyright (C) 1997 Makoto Matsumoto and Takuji Nishimura.
! When you use this, send an email to: matumoto@math.keio.ac.jp
! with an appropriate reference to your work.
! Fortran translation by Hiroshi Takano. Jan. 13, 1999.
!
! (2) On a single IBM Power4 processor on the NCEP operational cluster (2005)
! each Mersenne twister random number takes less than 30 ns, about 3 times
! slower than the default random number generator, and each random number
! from a Gaussian distribution takes less than 150 ns.
!
! Attributes:
! Language: Fortran 90
!
!$$$
module mersenne_twister
private
! Public declarations
public random_stat
public random_seed
public random_setseed
public random_number
public random_number_f
public random_gauss
public random_gauss_f
public random_index
public random_index_f
! Parameters
integer,parameter:: n=624
integer,parameter:: m=397
integer,parameter:: mata=-1727483681 ! constant vector a
integer,parameter:: umask=-2147483648 ! most significant w-r bits
integer,parameter:: lmask =2147483647 ! least significant r bits
integer,parameter:: tmaskb=-1658038656 ! tempering parameter
integer,parameter:: tmaskc=-272236544 ! tempering parameter
integer,parameter:: mag01(0:1)=(/0,mata/)
integer,parameter:: iseed=4357
integer,parameter:: nrest=n+6
! Defined types
type random_stat
private
integer:: mti=n+1
integer:: mt(0:n-1)
integer:: iset
real:: gset
end type
! Saved data
type(random_stat),save:: sstat
! Overloaded interfaces
interface random_setseed
module procedure random_setseed_s
module procedure random_setseed_t
end interface
interface random_number
module procedure random_number_i
module procedure random_number_s
module procedure random_number_t
end interface
interface random_gauss
module procedure random_gauss_i
module procedure random_gauss_s
module procedure random_gauss_t
end interface
interface random_index
module procedure random_index_i
module procedure random_index_s
module procedure random_index_t
end interface
! All the subprograms
contains
! Subprogram random_seed
! Sets and gets state; overloads Fortran 90 standard.
subroutine random_seed(size,put,get,stat)
implicit none
integer,intent(out),optional:: size
integer,intent(in),optional:: put(nrest)
integer,intent(out),optional:: get(nrest)
type(random_stat),intent(inout),optional:: stat
if(present(size)) then ! return size of seed array
! if(present(put).or.present(get))&
! call errmsg('RANDOM_SEED: more than one option set - some ignored')
size=nrest
elseif(present(put)) then ! restore from seed array
! if(present(get))&
! call errmsg('RANDOM_SEED: more than one option set - some ignored')
if(present(stat)) then
stat%mti=put(1)
stat%mt=put(2:n+1)
stat%iset=put(n+2)
stat%gset=transfer(put(n+3:nrest),stat%gset)
if(stat%mti.lt.0.or.stat%mti.gt.n.or.any(stat%mt.eq.0).or.
& stat%iset.lt.0.or.stat%iset.gt.1) then
call random_setseed_t(iseed,stat)
! call errmsg('RANDOM_SEED: invalid seeds put - default seeds used')
endif
else
sstat%mti=put(1)
sstat%mt=put(2:n+1)
sstat%iset=put(n+2)
sstat%gset=transfer(put(n+3:nrest),sstat%gset)
if(sstat%mti.lt.0.or.sstat%mti.gt.n.or.any(sstat%mt.eq.0)
& .or.sstat%iset.lt.0.or.sstat%iset.gt.1) then
call random_setseed_t(iseed,sstat)
! call errmsg('RANDOM_SEED: invalid seeds put - default seeds used')
endif
endif
elseif(present(get)) then ! save to seed array
if(present(stat)) then
if(stat%mti.eq.n+1) call random_setseed_t(iseed,stat)
get(1)=stat%mti
get(2:n+1)=stat%mt
get(n+2)=stat%iset
get(n+3:nrest)=transfer(stat%gset,get,nrest-(n+3)+1)
else
if(sstat%mti.eq.n+1) call random_setseed_t(iseed,sstat)
get(1)=sstat%mti
get(2:n+1)=sstat%mt
get(n+2)=sstat%iset
get(n+3:nrest)=transfer(sstat%gset,get,nrest-(n+3)+1)
endif
else ! reset default seed
if(present(stat)) then
call random_setseed_t(iseed,stat)
else
call random_setseed_t(iseed,sstat)
endif
endif
end subroutine
! Subprogram random_setseed_s
! Sets seed in saved mode.
subroutine random_setseed_s(inseed)
implicit none
integer,intent(in):: inseed
call random_setseed_t(inseed,sstat)
end subroutine
! Subprogram random_setseed_t
! Sets seed in thread-safe mode.
subroutine random_setseed_t(inseed,stat)
implicit none
integer,intent(in):: inseed
type(random_stat),intent(out):: stat
integer ii,mti
ii=inseed
if(ii.eq.0) ii=iseed
stat%mti=n
stat%mt(0)=iand(ii,-1)
do mti=1,n-1
stat%mt(mti)=iand(69069*stat%mt(mti-1),-1)
enddo
stat%iset=0
stat%gset=0.
end subroutine
! Subprogram random_number_f
! Generates random numbers in functional mode.
function random_number_f() result(harvest)
implicit none
real:: harvest
real h(1)
if(sstat%mti.eq.n+1) call random_setseed_t(iseed,sstat)
call random_number_t(h,sstat)
harvest=h(1)
end function
! Subprogram random_number_i
! Generates random numbers in interactive mode.
subroutine random_number_i(harvest,inseed)
implicit none
real,intent(out):: harvest(:)
integer,intent(in):: inseed
type(random_stat) stat
call random_setseed_t(inseed,stat)
call random_number_t(harvest,stat)
end subroutine
! Subprogram random_number_s
! Generates random numbers in saved mode; overloads Fortran 90 standard.
subroutine random_number_s(harvest)
implicit none
real,intent(out):: harvest(:)
if(sstat%mti.eq.n+1) call random_setseed_t(iseed,sstat)
call random_number_t(harvest,sstat)
end subroutine
! Subprogram random_number_t
! Generates random numbers in thread-safe mode.
subroutine random_number_t(harvest,stat)
implicit none
real,intent(out):: harvest(:)
type(random_stat),intent(inout):: stat
integer j,kk,y
integer tshftu,tshfts,tshftt,tshftl
tshftu(y)=ishft(y,-11)
tshfts(y)=ishft(y,7)
tshftt(y)=ishft(y,15)
tshftl(y)=ishft(y,-18)
do j=1,size(harvest)
if(stat%mti.ge.n) then
do kk=0,n-m-1
y=ior(iand(stat%mt(kk),umask),iand(stat%mt(kk+1),lmask))
stat%mt(kk)=ieor(ieor(stat%mt(kk+m),ishft(y,-1)),
& mag01(iand(y,1)))
enddo
do kk=n-m,n-2
y=ior(iand(stat%mt(kk),umask),iand(stat%mt(kk+1),lmask))
stat%mt(kk)=ieor(ieor(stat%mt(kk+(m-n)),ishft(y,-1)),
& mag01(iand(y,1)))
enddo
y=ior(iand(stat%mt(n-1),umask),iand(stat%mt(0),lmask))
stat%mt(n-1)=ieor(ieor(stat%mt(m-1),ishft(y,-1)),
& mag01(iand(y,1)))
stat%mti=0
endif
y=stat%mt(stat%mti)
y=ieor(y,tshftu(y))
y=ieor(y,iand(tshfts(y),tmaskb))
y=ieor(y,iand(tshftt(y),tmaskc))
y=ieor(y,tshftl(y))
if(y.lt.0) then
harvest(j)=(real(y)+2.0**32)/(2.0**32-1.0)
else
harvest(j)=real(y)/(2.0**32-1.0)
endif
stat%mti=stat%mti+1
enddo
end subroutine
! Subprogram random_gauss_f
! Generates Gaussian random numbers in functional mode.
function random_gauss_f() result(harvest)
implicit none
real:: harvest
real h(1)
if(sstat%mti.eq.n+1) call random_setseed_t(iseed,sstat)
call random_gauss_t(h,sstat)
harvest=h(1)
end function
! Subprogram random_gauss_i
! Generates Gaussian random numbers in interactive mode.
subroutine random_gauss_i(harvest,inseed)
implicit none
real,intent(out):: harvest(:)
integer,intent(in):: inseed
type(random_stat) stat
call random_setseed_t(inseed,stat)
call random_gauss_t(harvest,stat)
end subroutine
! Subprogram random_gauss_s
! Generates Gaussian random numbers in saved mode.
subroutine random_gauss_s(harvest)
implicit none
real,intent(out):: harvest(:)
if(sstat%mti.eq.n+1) call random_setseed_t(iseed,sstat)
call random_gauss_t(harvest,sstat)
end subroutine
! Subprogram random_gauss_t
! Generates Gaussian random numbers in thread-safe mode.
subroutine random_gauss_t(harvest,stat)
implicit none
real,intent(out):: harvest(:)
type(random_stat),intent(inout):: stat
integer mx,my,mz,j
real r2(2),r,g1,g2
mz=size(harvest)
if(mz.le.0) return
mx=0
if(stat%iset.eq.1) then
mx=1
harvest(1)=stat%gset
stat%iset=0
endif
my=(mz-mx)/2*2+mx
do
call random_number_t(harvest(mx+1:my),stat)
do j=mx,my-2,2
call rgauss(harvest(j+1),harvest(j+2),r,g1,g2)
if(r.lt.1.) then
harvest(mx+1)=g1
harvest(mx+2)=g2
mx=mx+2
endif
enddo
if(mx.eq.my) exit
enddo
if(my.lt.mz) then
do
call random_number_t(r2,stat)
call rgauss(r2(1),r2(2),r,g1,g2)
if(r.lt.1.) exit
enddo
harvest(mz)=g1
stat%gset=g2
stat%iset=1
endif
contains
! Numerical Recipes algorithm to generate Gaussian random numbers.
subroutine rgauss(r1,r2,r,g1,g2)
real,intent(in):: r1,r2
real,intent(out):: r,g1,g2
real v1,v2,fac
v1=2.*r1-1.
v2=2.*r2-1.
r=v1**2+v2**2
if(r.lt.1.) then
fac=sqrt(-2.*log(r)/r)
g1=v1*fac
g2=v2*fac
endif
end subroutine
end subroutine
! Subprogram random_index_f
! Generates random indices in functional mode.
function random_index_f(imax) result(iharvest)
implicit none
integer,intent(in):: imax
integer:: iharvest
integer ih(1)
if(sstat%mti.eq.n+1) call random_setseed_t(iseed,sstat)
call random_index_t(imax,ih,sstat)
iharvest=ih(1)
end function
! Subprogram random_index_i
! Generates random indices in interactive mode.
subroutine random_index_i(imax,iharvest,inseed)
implicit none
integer,intent(in):: imax
integer,intent(out):: iharvest(:)
integer,intent(in):: inseed
type(random_stat) stat
call random_setseed_t(inseed,stat)
call random_index_t(imax,iharvest,stat)
end subroutine
! Subprogram random_index_s
! Generates random indices in saved mode.
subroutine random_index_s(imax,iharvest)
implicit none
integer,intent(in):: imax
integer,intent(out):: iharvest(:)
if(sstat%mti.eq.n+1) call random_setseed_t(iseed,sstat)
call random_index_t(imax,iharvest,sstat)
end subroutine
! Subprogram random_index_t
! Generates random indices in thread-safe mode.
subroutine random_index_t(imax,iharvest,stat)
implicit none
integer,intent(in):: imax
integer,intent(out):: iharvest(:)
type(random_stat),intent(inout):: stat
integer,parameter:: mh=n
integer i1,i2,mz
real h(mh)
mz=size(iharvest)
do i1=1,mz,mh
i2=min((i1-1)+mh,mz)
call random_number_t(h(:i2-(i1-1)),stat)
iharvest(i1:i2)=max(ceiling(h(:i2-(i1-1))*imax),1)
enddo
end subroutine
end module