ModRandom.F90

Go to the documentation of this file.

! *********************************************************************
! * Rocstar Simulation Suite                                          *
! * Copyright@2015, Illinois Rocstar LLC. All rights reserved.        *
! *                                                                   *
! * Illinois Rocstar LLC                                              *
! * Champaign, IL                                                     *
! * www.illinoisrocstar.com                                           *
! * sales@illinoisrocstar.com                                         *
! *                                                                   *
! * License: See LICENSE file in top level of distribution package or *
! * http://opensource.org/licenses/NCSA                               *
! *********************************************************************
! *********************************************************************
! * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,   *
! * EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES   *
! * OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND          *
! * NONINFRINGEMENT.  IN NO EVENT SHALL THE CONTRIBUTORS OR           *
! * COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER       *
! * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE,   *
! * Arising FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE    *
! * USE OR OTHER DEALINGS WITH THE SOFTWARE.                          *
! *********************************************************************
!******************************************************************************
!
! Purpose: Portable random number generator
!
! Description: This module provides the following routines:
!
!   RandomSeed: seeds the random number generator
!
!   RandUniform:    generates an array of uniformly distributed random numbers
!                   on the interval [0,1)
!
!   RandNormal:     generates an array of normally distributed random numbers
!
!   RandLogNormal:  generates an array of log-normally distributed numbers
!
!   RandImposedPDF: generates an array of number distributed according to a dscrete PDF
!
!   Rand1Uniform:    function form of RandUniform -- returns a single number
!
!   Rand1Normal:     function form of RandNormal
!
!   Rand1LogNormal:  function form of RandLogNormal
!
!   Rand1ImposedPDF: function form of RandImposedPDF
!
! Notes:
!
!   This is an implementation of the Mersenne Twister random number
!   generator, MT19937, based on the Mersenne Prime 2^19937 - 1.
!
!   It is (somewhat painfully) coded such that it gives the same results
!   whether integers are 32- or 64-bit, and these results are identical
!   to the standard, C version of the Mersenne Twister.
!
!   It outputs random numbers with 53-bit precision.
!
!******************************************************************************
!
! $Id: ModRandom.F90,v 1.6 2008/12/06 08:44:19 mtcampbe Exp $
!
! Copyright: (c) 2003 by the University of Illinois
!
!******************************************************************************

MODULE modrandom

  USE moddatatypes
  IMPLICIT NONE

! Period parameter
  INTEGER, PARAMETER :: RAND_BUFFER_SIZE = 624
  INTEGER, PARAMETER :: RAND_MTI_INDEX   = RAND_BUFFER_SIZE
  INTEGER, PARAMETER :: RAND_SEED_INDEX  = RAND_BUFFER_SIZE + 1
  INTEGER, PARAMETER :: RAND_CALLS_INDEX = RAND_BUFFER_SIZE + 2
  INTEGER, PARAMETER :: RAND_EXTRA_INTS  = 3
  INTEGER, PARAMETER :: RAND_TOTAL_SIZE  = RAND_BUFFER_SIZE + RAND_EXTRA_INTS

  TYPE t_rand_data
    INTEGER :: mt(0:RAND_TOTAL_SIZE-1)
  END TYPE t_rand_data

CONTAINS

! ******************************************************************************
! Initialization subroutine
! ******************************************************************************
!   sets initial seeds to rdata%mt using the generator Line 25 of Table 1 in
!   [KNUTH 1981, The Art of Computer Programming Vol. 2 (2nd Ed.), pp102]

  SUBROUTINE randomseed(seed,rdata)

! ... parameters
    INTEGER,           INTENT(IN)  :: seed
    TYPE(t_rand_data), INTENT(OUT) :: rdata

! ... loop variables
    INTEGER :: i

! ... local variables
    INTEGER, PARAMETER :: one = 1
    INTEGER, PARAMETER :: mult = 1812433253

    INTEGER :: mask32

    mask32 = ishft(one,32) - one

    rdata%mt(0) = iand(seed,mask32)

    DO i=1,rand_buffer_size-1
! --- See Knuth TAOCP Vol2. 3rd Ed. P.106 for multiplier
! --- In the previous versions, MSBs of the seed affect
! --- only MSBs of the array mt[].
! --- 2002/01/09 modified by Makoto Matsumoto

      rdata%mt(i) = mult * (ieor(rdata%mt(i-1), &
                                 ishft(rdata%mt(i-1),-30))) + i

! --- for >32 bit machines
      rdata%mt(i) = iand(rdata%mt(i),mask32)

    END DO ! i

    rdata%mt(rand_mti_index)   = rand_buffer_size
    rdata%mt(rand_seed_index)  = iand(seed,mask32)
    rdata%mt(rand_calls_index) = 0

  END SUBROUTINE randomseed

! ******************************************************************************
! Uniform random number generator on [0,1) with 53-bit precision
! ******************************************************************************

  SUBROUTINE randuniform(a,rdata)

! ... parameters
    REAL(RFREAL),      INTENT(OUT)   :: a(:)
    TYPE(t_rand_data), INTENT(INOUT) :: rdata

! ... loop variables
    INTEGER :: i,j,kk

! ... local variables

! - Period parameters
    INTEGER, PARAMETER :: one = 1
    INTEGER, PARAMETER :: m   = 397
    INTEGER, PARAMETER :: lmask      =  2147483647 ! least significant r bits
    INTEGER, PARAMETER :: pre_mata   = -1727483681
    INTEGER, PARAMETER :: pre_tmaskb = -1658038656
    INTEGER, PARAMETER :: pre_tmaskc = -272236544

! - Normalization parameters
    REAL(RFREAL), PARAMETER :: two32 = 4294967296.0_rfreal
    REAL(RFREAL), PARAMETER :: nfac1 = 67108864.0_rfreal

! - These should be parameters, but cannot be initialized on all platforms
    INTEGER      :: mask32, mata, umask, tmaskb, tmaskc, mag01(0:1)
    REAL(RFREAL) :: nfac2

    INTEGER      :: y,mti,iamin,iamax
    REAL(RFREAL) :: yr(2)

! - Set values of things that should, ideally, be parameters

    umask    = ishft(one,31)           ! most significant w-r bits
    mask32   = ishft(one,32) - one     ! mask to force 32-bit behavior
    mata     = iand(pre_mata,  mask32) ! constant vector
    tmaskb   = iand(pre_tmaskb,mask32) ! tempering parameter
    tmaskc   = iand(pre_tmaskc,mask32) ! tempering parameter
    mag01(0) = 0
    mag01(1) = mata
    nfac2    = 1.0_rfreal/9007199254740992.0_rfreal

! - Begin

    mti = rdata%mt(rand_mti_index)

    iamin = lbound(a,1)
    iamax = ubound(a,1)

    DO i = iamin,iamax

      DO j = 1,2

        IF (mti >= rand_buffer_size) THEN ! make RAND_BUFFER_SIZE integers

          DO kk=0,rand_buffer_size-m-1
            y = ior(iand(rdata%mt(kk),umask),iand(rdata%mt(kk+1),lmask))
            rdata%mt(kk) = ieor(ieor(rdata%mt(kk+m),ishft(y,-1)), &
                                mag01(iand(y,one)))
          END DO ! kk

          DO kk=rand_buffer_size-m,rand_buffer_size-2
            y = ior(iand(rdata%mt(kk),umask),iand(rdata%mt(kk+1),lmask))
            rdata%mt(kk) = ieor(ieor(rdata%mt(kk+(m-rand_buffer_size)), &
                                     ishft(y,-1)), mag01(iand(y,one)))
          END DO ! kk

          y = ior(iand(rdata%mt(rand_buffer_size-1),umask), &
                  iand(rdata%mt(0),lmask))
          rdata%mt(rand_buffer_size-1) = &
            ieor(ieor(rdata%mt(m-1), ishft(y,-1)), mag01(iand(y,one)))

          mti = 0

        END IF ! mti

        y = rdata%mt(mti)

! ----- Tempering

        y = ieor(y,     ishft(y,-11))
        y = ieor(y,iand(ishft(y,  7),tmaskb))
        y = ieor(y,iand(ishft(y, 15),tmaskc))
        y = ieor(y,     ishft(y,-18))

        y = ishft(y,-4-j)

        IF (y < 0) THEN
          yr(j) = REAL(y,KIND=RFREAL) + two32
        ELSE
          yr(j) = REAL(y,kind=rfreal)
        END IF ! y

        mti = mti + 1

      END DO ! j

! --- Form random number from two 32-bit integers

      a(i) = (yr(1)*nfac1 + yr(2)) * nfac2

    END DO ! i

    rdata%mt(rand_mti_index)   = mti
    rdata%mt(rand_calls_index) = rdata%mt(rand_calls_index) + iamax - iamin + 1

  END SUBROUTINE randuniform

! ******************************************************************************
! Normal distribution
! ******************************************************************************

  SUBROUTINE randnormal(a,med,dev,rdata)
! med = median
! dev = standard deviation

    REAL(RFREAL),      INTENT(OUT)   :: a(:)
    REAL(RFREAL),      INTENT(IN)    :: med,dev
    TYPE(t_rand_data), INTENT(INOUT) :: rdata

    INTEGER :: i,iamin,iamax
    REAL(RFREAL) :: r,th
    REAL(RFREAL) :: twopi

    iamin = lbound(a,1)
    iamax = ubound(a,1)
    twopi = 8.0_rfreal*atan(1.0_rfreal)
    CALL randuniform(a,rdata)

    DO i=iamin,iamax-1,2
      th = twopi*a(i)
      r  = dev*sqrt(-2.0_rfreal*log(1.0_rfreal-a(i+1)))
      a(i)   = r*cos(th) + med
      a(i+1) = r*sin(th) + med
    END DO ! i

    IF (modulo(iamax-iamin+1,2)==1) THEN
      th = twopi*a(iamax)
      r  = rand1uniform(rdata)
      r  = dev*sqrt(-2.0_rfreal*log(1.0_rfreal-r))
      a(iamax) = r*cos(th) + med
    END IF

  END SUBROUTINE randnormal

! ******************************************************************************
! Log normal distribution
! ******************************************************************************

  SUBROUTINE randlognormal(a,logmed,gdev,rdata)
! logmed = log of median
! gdev = standard deviation of corresponding Gaussian

    REAL(RFREAL),      INTENT(OUT)   :: a(:)
    REAL(RFREAL),      INTENT(IN)    :: logmed,gdev
    TYPE(t_rand_data), INTENT(INOUT) :: rdata

    CALL randnormal(a,logmed,gdev,rdata)
    a = exp(a)

  END SUBROUTINE randlognormal

! ******************************************************************************
! Discrete PDF distribution from file
! ******************************************************************************

  SUBROUTINE randimposedpdf(a,pdfvalues,locMax,valmax,rdata)

    REAL(RFREAL),      INTENT(IN)    :: pdfvalues(:,:)
    REAL(RFREAL),      INTENT(OUT)   :: a(:)
    TYPE(t_rand_data), INTENT(INOUT) :: rdata
    INTEGER,           INTENT(IN)    :: locmax
    REAL(RFREAL),      INTENT(IN)    :: valmax

    INTEGER :: i,iamin,iamax,k1,k2

    iamin = lbound(a,1)
    iamax = ubound(a,1)
    CALL randuniform(a,rdata)

    DO i=iamin,iamax
       call findpdf(a(i),pdfvalues(:,2),k1,k2)
       a(i) = (a(i) - pdfvalues(k1,2))*(pdfvalues(k2,3)-pdfvalues(k1,3))&
             /(pdfvalues(k2,2) - pdfvalues(k1,2)) + pdfvalues(k1,3)
    END DO ! i

    RETURN
    
  CONTAINS
    SUBROUTINE findpdf(datum,vec,k1,k2)
      REAL(RFREAL),  INTENT(IN)    :: datum
      REAL(RFREAL),  INTENT(IN)    :: vec(:)
      INTEGER,       INTENT(OUT)   :: k1,k2  
      integer :: isgn

      isgn = sign(1.1_rfreal,datum-valmax)
      k1 = locmax
      k2 = k1+isgn
      do while (int(sign(1.1_rfreal,datum-vec(k2)))*isgn > 0)
         k1 = k2
         k2 = k2+isgn
      enddo
    END SUBROUTINE findpdf
      

  END SUBROUTINE randimposedpdf

! ******************************************************************************
! Function form of RandUniform 
! ******************************************************************************

  REAL(RFREAL) FUNCTION rand1uniform(rdata)

! ... parameters
    TYPE(t_rand_data), INTENT(INOUT) :: rdata

! ... local variables
    REAL(RFREAL) :: a(1)

    CALL randuniform(a,rdata)
    rand1uniform = a(1)

  END FUNCTION rand1uniform

! ******************************************************************************
! Function form of RandNormal
! ******************************************************************************

  REAL(RFREAL) FUNCTION rand1normal(med,dev,rdata)

! ... parameters
    REAL(RFREAL),      INTENT(IN)    :: med,dev
    TYPE(t_rand_data), INTENT(INOUT) :: rdata

! ... local variables
    REAL(RFREAL) :: a(1)

    CALL randnormal(a,med,dev,rdata)
    rand1normal = a(1)

  END FUNCTION rand1normal

! ******************************************************************************
! Function form of RandLogNormal
! ******************************************************************************

  REAL(RFREAL) FUNCTION rand1lognormal(logmed,gdev,rdata)

! ... parameters
    REAL(RFREAL),      INTENT(IN)    :: logmed,gdev
    TYPE(t_rand_data), INTENT(INOUT) :: rdata

! ... local variables
    REAL(RFREAL) :: a(1)

    CALL randlognormal(a,logmed,gdev,rdata)
    rand1lognormal = a(1)

  END FUNCTION rand1lognormal

! ******************************************************************************
! Function form of RandUniform 
! ******************************************************************************

  REAL(RFREAL) FUNCTION rand1imposedpdf(rdata,pdfvalues,locMax,valmax)

! ... parameters
    TYPE(t_rand_data), INTENT(INOUT) :: rdata
    REAL(RFREAL),      INTENT(IN)    :: pdfvalues(:,:)
    INTEGER,           INTENT(IN)    :: locmax
    REAL(RFREAL),      INTENT(IN)    :: valmax

! ... local variables
    REAL(RFREAL) :: a(1)

    CALL randimposedpdf(a,pdfvalues,locmax,valmax,rdata)
    rand1imposedpdf = a(1)

  END FUNCTION rand1imposedpdf
END MODULE modrandom

!******************************************************************************
!
! RCS Revision history:
!
! $Log: ModRandom.F90,v $
! Revision 1.6  2008/12/06 08:44:19  mtcampbe
! Updated license.
!
! Revision 1.5  2008/11/19 22:17:30  mtcampbe
! Added Illinois Open Source License/Copyright
!
! Revision 1.4  2005/04/25 18:39:09  luca1
! Imposed PDF from file option for random particle ejection
!
! Revision 1.3  2003/12/08 16:39:47  jferry
! Removed non-standard KIND intrinsic
!
! Revision 1.2  2003/11/21 22:35:51  fnajjar
! Update Random Number Generator
!
! Revision 1.1  2003/02/17 19:31:11  jferry
! Implemented portable random number generator ModRandom
!
!******************************************************************************