ran_state.F

Go to the documentation of this file.

        MODULE ran_state
!
!git $Id$
!================================================== Hernan G. Arango ===
!  Copyright (c) 2002-2025 The ROMS Group                              !
!    Licensed under a MIT/X style license                              !
!    See License_ROMS.md                                               !
!=======================================================================
!                                                                      !
!  This module supports the random number generation routines. It      !
!  provides each generator with five vectors integers, for use as      !
!  internal random state space.  The  first three integers (iran,      !
!  jran, kran) are maintained as nonnegative  values,  while  the      !
!  last two (mran, nran) have 32-bit nonzero values.  This module      !
!  also includes support for initializing  or  reinitializing the      !
!  state to a desired sequence number, hashing the initial values      !
!  to randon values, and allocating and deallocateing of internal      !
!  workspace.                                                          !
!                                                                      !
!  Adapted from Numerical Recepies.                                    !
!                                                                      !
!  Press, W.H., S.A. Teukolsky, W.T. Vetterling, and B.P. Flannery,    !
!     1996:  Numerical Recipes in Fortran 90,  The Art of Parallel     !
!     Scientific Computing, 2nd Edition, Cambridge Univ. Press.        !
!                                                                      !
!=======================================================================
!
      USE mod_kinds
 
      implicit none
 
      PUBLIC
 
      integer(i8b), parameter :: hg = huge(1_i8b)
      integer(i8b), parameter :: hgm = -hg
      integer(i8b), parameter :: hgng = hgm - 1
 
      integer(i8b), save :: lenran = 0
      integer(i8b), save :: seq = 0
      integer(i8b), save :: iran0
      integer(i8b), save :: jran0
      integer(i8b), save :: kran0
      integer(i8b), save :: nran0
      integer(i8b), save :: mran0
      integer(i8b), save :: rans
 
      integer(i8b), pointer, save :: iran(:)
      integer(i8b), pointer, save :: jran(:)
      integer(i8b), pointer, save :: kran(:)
      integer(i8b), pointer, save :: nran(:)
      integer(i8b), pointer, save :: mran(:)
      integer(i8b), pointer, save :: ranv(:)
 
      integer(i8b), pointer, save :: ranseeds(:,:)
 
      real(r8), save :: amm
!
      INTERFACE ran_hash
        MODULE PROCEDURE ran_hash_s, ran_hash_v
      END INTERFACE
!
      CONTAINS
 
      SUBROUTINE ran_init (length)
!
!=======================================================================
!                                                                      !
!  This routine initializes or reinitializes the random generator      !
!  state space to vectors of size LENGTH.  The saved variable SEQ      !
!  is hashed (via a call to RAN_HASH) to create  unique  starting      !
!  seeds, different for each vector component.                         !
!                                                                      !
!=======================================================================
!
      USE nrutil, ONLY : arth, nrerror, reallocate
!
!  Imported variable declarations.
!
      integer(i8b), intent(in) :: length
!
!  Local variable declarations.
!
      integer(i8b) :: hgt, j, new, sz
!
!-----------------------------------------------------------------------
!  Initialize randon number generator vectors.
!-----------------------------------------------------------------------
!
      IF (length.lt.lenran) RETURN
      hgt=hg
!
!  Check that kind value I8B is in fact a 32-bit integer with the usual
!  properties that we expect it to have (under negation and wrap-around
!  addition). If all these test are satisfied, then the routines that
!  use this module are portable, even though they go beyond F90 integer
!  model.
!
      IF (hg.ne.2147483647)                                             &
     &  CALL nrerror ('RAN_INIT: arith assump 1 fails')
      IF (hgng.ge.0)                                                    &
     &  CALL nrerror ('RAN_INIT: arith assump 2 fails')
      IF ((hgt+1).ne.hgng)                                              &
     &  CALL nrerror ('RAN_INIT: arith assump 3 fails')
      IF (not(hg).ge.0)                                                 &
     &  CALL nrerror ('RAN_INIT: arith assump 4 fails')
      IF (not(hgng).lt.0)                                               &
     &  CALL nrerror ('RAN_INIT: arith assump 5 fails')
      IF ((hg+hgng).ge.0)                                               &
     &  CALL nrerror ('RAN_INIT: arith assump 6 fails')
      IF (not(-1_i8b).lt.0)                                             &
     &  CALL nrerror ('RAN_INIT: arith assump 7 fails')
      IF (not(0_i8b).ge.0)                                              &
     &  CALL nrerror ('RAN_INIT: arith assump 8 fails')
      IF (not(1_i8b).ge.0)                                              &
     &  CALL nrerror ('RAN_INIT: arith assump 9 fails')
!
!  Reallocate or allocate state space.
!
      IF (lenran.gt.0) THEN
        ranseeds => reallocate(ranseeds, length, 5_i8b)
        ranv => reallocate(ranv, length-1_i8b)
        new=lenran+1
      ELSE
        ALLOCATE (ranseeds(length,5))
        ALLOCATE (ranv(length-1))
        new=1
        amm=nearest(1.0_r8,-1.0_r8)/hgng
        IF ((amm*hgng.ge.1.0_r8).or.(amm*hgng.le.0.0_r8))               &
     &    CALL nrerror ('RAN_INIT: arth assump 10 fails')
      END IF
!
!  Set starting values, unique by SEQ and vector component.
!
      ranseeds(new:,1)=seq
      sz=SIZE(ranseeds(new:,1))
      ranseeds(new:,2:5)=spread(arth(new,1_i8b,sz),2,4)
!
!  Hash them.
!
      DO j=1,4
        CALL ran_hash (ranseeds(new:,j), ranseeds(new:,j+1))
      END DO
!
!  Enforce nonnegativity.
!
      WHERE (ranseeds(new:,1:3).lt.0)                                   &
     &  ranseeds(new:,1:3)=not(ranseeds(new:,1:3))
!
!  Enforce nonzero.
!
      WHERE (ranseeds(new:,4:5).eq.0)                                   &
     &   ranseeds(new:,4:5)=1
!
!  Set scalar seeds.
!
      IF (new.eq.1) THEN
        iran0=ranseeds(1,1)
        jran0=ranseeds(1,2)
        kran0=ranseeds(1,3)
        mran0=ranseeds(1,4)
        nran0=ranseeds(1,5)
        rans=nran0
      END IF
!
!  Point to vector seeds.
!
      IF (length.gt.1) THEN
        iran => ranseeds(2:,1)
        jran => ranseeds(2:,2)
        kran => ranseeds(2:,3)
        mran => ranseeds(2:,4)
        nran => ranseeds(2:,5)
        ranv = nran
      END IF
      lenran=length
 
      END SUBROUTINE ran_init
 
      SUBROUTINE ran_deallocate
!
!=======================================================================
!                                                                      !
!  User interface to release the workspace used by random number       !
!  routines.                                                           !
!                                                                      !
!=======================================================================
!
      IF (lenran.gt.0) THEN
        DEALLOCATE (ranseeds, ranv)
        NULLIFY (ranseeds, ranv, iran, jran, kran, mran, nran)
        lenran=0
      END IF
      END SUBROUTINE ran_deallocate
 
      SUBROUTINE ran_seed (sequence, size, put, get)
!
!=======================================================================
!                                                                      !
!  User interface for seeding the random number routines.  Syntax is   !
!  exactly like Fortran 90 RANDOM_SEED, with one additional argument   !
!  keyword: SEQUENCE, set to any integer value,  causes an immediate   !
!  new initialization, seeded by that integer.                         !
!                                                                      !
!=======================================================================
!
!  Imported variable declarations.
!
      integer, optional, intent(in) :: sequence
      integer, optional, intent(out) :: size
 
      integer, optional, intent(in) :: put(:)
      integer, optional, intent(out) :: get(:)
!
!-----------------------------------------------------------------------
!  Set random number seeds.
!-----------------------------------------------------------------------
!
      IF (PRESENT(size)) THEN
        size=5*lenran
      ELSE IF (PRESENT(put)) THEN
        IF (lenran.eq.0) RETURN
        ranseeds=reshape(put,shape(ranseeds))
        WHERE (ranseeds(:,1:3).lt.0)                                    &
     &    ranseeds(:,1:3)=not(ranseeds(:,1:3))
        WHERE (ranseeds(:,4:5).eq.0)                                    &
     &    ranseeds(:,4:5)=1
          iran0=ranseeds(1,1)
          jran0=ranseeds(1,2)
          kran0=ranseeds(1,3)
          mran0=ranseeds(1,4)
          nran0=ranseeds(1,5)
      ELSE IF (present(get)) THEN
        IF (lenran.eq.0) RETURN
        ranseeds(1,1:5)=(/ iran0,jran0,kran0,mran0,nran0 /)
        get=reshape(ranseeds,shape(get))
      ELSE IF (PRESENT(sequence)) THEN
        CALL ran_deallocate
        seq=sequence
      END IF
 
      RETURN
      END SUBROUTINE ran_seed
 
      SUBROUTINE ran_hash_s (il, ir)
!
!=======================================================================
!                                                                      !
!  DES-like hashing of 32-bit integer, using shifts, xor, and adds to  !
!  make the interval nonlinear function. Scalar version.               !
!                                                                      !
!=======================================================================
!
!  Imported variable declarations.
!
      integer(i8b), intent(inout) :: il, ir
!
!  Local variable declarations.
!
      integer(i8b) :: is, j
!
!-----------------------------------------------------------------------
!  Bit mixing. The various constants should not be changed.
!-----------------------------------------------------------------------
!
      DO j=1,4
        is=ir
        ir=ieor(ir,ishft(ir,5))+1422217823
        ir=ieor(ir,ishft(ir,-16))+1842055030
        ir=ieor(ir,ishft(ir,9))+80567781
        ir=ieor(il,ir)
        il=is
      END DO
 
      RETURN
      END SUBROUTINE ran_hash_s
 
      SUBROUTINE ran_hash_v (il, ir)
!
!=======================================================================
!                                                                      !
!  DES-like hashing of 32-bit integer, using shifts, xor, and adds to  !
!  make the interval nonlinear function. Vector version.               !
!                                                                      !
!=======================================================================
!
!  Imported variable declarations.
!
      integer(i8b), intent(inout) :: il(:)
      integer(i8b), intent(inout) :: ir(:)
!
!  Local variable declarations.
!
      integer(i8b) :: j
      integer(i8b), dimension(SIZE(il)) :: is
!
!-----------------------------------------------------------------------
!  Bit mixing. The various constants should not be changed.
!-----------------------------------------------------------------------
!
      DO j=1,4
        is=ir
        ir=ieor(ir,ishft(ir,5))+1422217823
        ir=ieor(ir,ishft(ir,-16))+1842055030
        ir=ieor(ir,ishft(ir,9))+80567781
        ir=ieor(il,ir)
        il=is
      END DO
 
      RETURN
      END SUBROUTINE ran_hash_v
 
      END MODULE ran_state