rngs.f90

Go to the documentation of this file.

!#############################################################################
!#                                                                           #
!# fosite - 3D hydrodynamical simulation program                             #
!# module: rngs.f90                                                          #
!#                                                                           #
!# Copyright (C) 2015 Manuel Jung <mjung@astrophysik.uni-kiel.de>            #
!#                                                                           #
!# This program is free software; you can redistribute it and/or modify      #
!# it under the terms of the GNU General Public License as published by      #
!# the Free Software Foundation; either version 2 of the License, or (at     #
!# your option) any later version.                                           #
!#                                                                           #
!# This program 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, GOOD TITLE or        #
!# NON INFRINGEMENT.  See the GNU General Public License for more            #
!# details.                                                                  #
!#                                                                           #
!# You should have received a copy of the GNU General Public License         #
!# along with this program; if not, write to the Free Software               #
!# Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.                 #
!#                                                                           #
!#############################################################################

!----------------------------------------------------------------------------!
!----------------------------------------------------------------------------!
MODULE rngs
  IMPLICIT NONE
  !--------------------------------------------------------------------------!
  PRIVATE
#if defined(NECSXAURORA)
  INTEGER, PARAMETER :: i8 = 8
#else
  INTEGER, PARAMETER :: i8 = selected_int_kind(18)
#endif
  !--------------------------------------------------------------------------!
  PUBLIC :: &
       i8, &
       kiss64, &
       superkiss64, &
       dkiss64
  !--------------------------------------------------------------------------!

CONTAINS

  FUNCTION kiss64(seed) RESULT(res)
    IMPLICIT NONE
    !------------------------------------------------------------------------!
    INTEGER(KIND=I8)           :: res
    INTEGER(KIND=I8), OPTIONAL :: seed
    !------------------------------------------------------------------------!
    INTEGER(KIND=I8), save     :: x, y, z, c
    INTEGER(KIND=I8)           :: t
    data x, y, z, c &
      / 1234567890987654321_i8, &
        362436362436362436_i8, &
        1066149217761810_i8, &
        123456123456123456_i8 /
    !------------------------------------------------------------------------!
    IF(PRESENT(seed)) THEN
      z = seed
      res = 0
    ELSE
      t = ishft(x, 58) + c
      IF (ishft(x,-63_i8) .EQ. ishft(t,-63_i8)) THEN
        c = ishft(x, -6) + ishft(x, -63_i8)
      ELSE
        c = ishft(x, -6) + 1 - ishft(x+t, -63_i8)
      END IF
      x = t + x
      y = ieor(y, ishft(y,13_i8))
      y = ieor(y, ishft(y,-17_i8))
      y = ieor(y, ishft(y,43_i8))

      z = 6906969069_i8 * z + 1234567
      res = x + y + z
    END IF
  END FUNCTION kiss64

  FUNCTION superkiss64() RESULT(x)
    IMPLICIT NONE
    !------------------------------------------------------------------------!
    INTEGER(KIND=I8)       :: x
    !------------------------------------------------------------------------!
    INTEGER(KIND=I8), save :: i,q(20632),carry,xcng,xs,indx
    data i,carry,xcng,xs,indx &
      / 0_i8, &
        36243678541_i8, &
        12367890123456_i8,&
        521288629546311_i8,&
        20633_i8 /
    !------------------------------------------------------------------------!
! Nec SX internal! compiler error:
! "f90 fatal: Internal error in optimization phase."
#if !defined(NECSXAURORA)
    !Initalize if i=0
    IF(i.EQ.0) THEN
      DO i=1,20632
        !fill Q with Congruential+Xorshift
        xcng = xcng*6906969069_i8+123
        xs = ieor(xs,ishft(xs,13))
        xs = ieor(xs,ishft(xs,-17))
        xs = ieor(xs,ishft(xs,43))
        q(i) = xcng + xs
      END DO
    END IF

    IF(indx <= 20632) THEN
      x=q(indx)
      indx=indx+1
    ELSE
      x=refill()
    END IF

    xcng=xcng*6906969069_i8+123
    xs=ieor(xs,ishft(xs,13))
    xs=ieor(xs,ishft(xs,-17))
    xs=ieor(xs,ishft(xs,43))
    x=x+xcng+xs
  CONTAINS

    FUNCTION refill() RESULT(s)
      IMPLICIT NONE
      !------------------------------------------------------------------------!
      INTEGER(KIND=I8) :: i,s,z,h
      !------------------------------------------------------------------------!
      DO i=1,20632
        h=iand(carry,1_i8)
        z = ishft(ishft(q(i),41),-1) &
            + ishft(ishft(q(i),39),-1) &
            + ishft(carry,-1)
        carry=ishft(q(i),-23)+ishft(q(i),-25)+ishft(z,-63)
        q(i)=not(ishft(z,1)+h)
      END DO
      indx=2
      s=q(1)
      RETURN
    END FUNCTION refill
#else
    x = 0_i8
#endif
  END FUNCTION superkiss64

  FUNCTION i8tod(x) RESULT(res)
    IMPLICIT NONE
    !------------------------------------------------------------------------!
    INTEGER(KIND=I8) :: x
    REAL             :: res
#if defined(NECSXAURORA)
    INTEGER(KIND=8)  :: long
#endif
    !------------------------------------------------------------------------!
#if defined(NECSXAURORA)
    ! No idea what is happeing on the SX.., but HUGE(x)!=HUGE(long)
    ! The 2048? This is 2.*1024. No idea, why the 1024 is needed. This shouldn't
    ! be required.
    res = (REAL(x) / huge(long)) / 2048. + 0.5
#else
    res = 0.5e+0 * (REAL(x) / REAL(huge(x))) + 0.5e+0
#endif
  END FUNCTION i8tod


  FUNCTION dkiss64() RESULT(res)
    IMPLICIT NONE
    !------------------------------------------------------------------------!
    REAL             :: res
    !------------------------------------------------------------------------!
    res = i8tod(kiss64())
  END FUNCTION dkiss64


END MODULE rngs