mersenne_twister.F90

!>\file mersenne_twister.f
!! This file contains the module that calculates random numbers using the 
!! Mersenne twister

!> \defgroup mersenne_ge Mersenne Twister Module
!! Module: mersenne_twister   Modern random number generator
!!\author 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
        use kinddef, only: kind_dbl_prec
        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                         !< Generator state 
          private
          integer:: mti=n+1
          integer:: mt(0:n-1)
          integer:: iset
          real(kind_dbl_prec):: 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
!> This subroutine 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
            size=nrest
          elseif(present(put)) then  ! restore from seed array
            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
!> This subroutine 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
!> This subroutine 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
!> This function generates random numbers in functional mode.
        function random_number_f() result(harvest)
          implicit none
          real(kind_dbl_prec):: harvest
          real(kind_dbl_prec) :: 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
!> This subroutine generates random numbers in interactive mode.
        subroutine random_number_i(harvest,inseed)
          implicit none
          real(kind_dbl_prec),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
!> This subroutine generates random numbers in saved mode; overloads Fortran 90 standard.
        subroutine random_number_s(harvest)
          implicit none
          real(kind_dbl_prec),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
!> This subroutine generates random numbers in thread-safe mode.
        subroutine random_number_t(harvest,stat)
          implicit none
          real(kind_dbl_prec),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,kind_dbl_prec)+2.0_kind_dbl_prec**32)/(2.0_kind_dbl_prec**32-1.0_kind_dbl_prec)
            else
              harvest(j)=real(y,kind_dbl_prec)/(2.0_kind_dbl_prec**32-1.0_kind_dbl_prec)
            endif
            stat%mti=stat%mti+1
          enddo
        end subroutine
!  Subprogram random_gauss_f
!> This subrouitne generates Gaussian random numbers in functional mode.
        function random_gauss_f() result(harvest)
          implicit none
          real(kind_dbl_prec):: harvest
          real(kind_dbl_prec) :: 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
!> This subrouitne generates Gaussian random numbers in interactive mode.
        subroutine random_gauss_i(harvest,inseed)
          implicit none
          real(kind_dbl_prec),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
!> This subroutine generates Gaussian random numbers in saved mode.
        subroutine random_gauss_s(harvest)
          implicit none
          real(kind_dbl_prec),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
!> This subroutine generates Gaussian random numbers in thread-safe mode.
        subroutine random_gauss_t(harvest,stat)
          implicit none
          real(kind_dbl_prec),intent(out):: harvest(:)
          type(random_stat),intent(inout):: stat
          integer mx,my,mz,j
          real(kind_dbl_prec) :: 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
!> This subroutine contains numerical Recipes algorithm to generate Gaussian random numbers.
          subroutine rgauss(r1,r2,r,g1,g2)
            real(kind_dbl_prec),intent(in):: r1,r2
            real(kind_dbl_prec),intent(out):: r,g1,g2
            real(kind_dbl_prec) :: v1,v2,fac
            v1=2.*r1-1._kind_dbl_prec
            v2=2.*r2-1._kind_dbl_prec
            r=v1**2+v2**2
            if(r.lt.1.) then
              fac=sqrt(-2._kind_dbl_prec*log(r)/r)
              g1=v1*fac
              g2=v2*fac
            endif
          end subroutine
        end subroutine
!  Subprogram random_index_f
!> This subroutine 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
!> This subroutine 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
!> This subroutine 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
!> This subroutine  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(kind_dbl_prec) :: 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