mirror of
https://github.com/cosmo-sims/MUSIC.git
synced 2024-09-20 18:13:45 +02:00
684 lines
17 KiB
Fortran
684 lines
17 KiB
Fortran
|
!=====================================================================================c
|
||
|
!
|
||
|
! The code below was written by: Stephen Booth
|
||
|
! Edinburgh Parallel Computing Centre
|
||
|
! The University of Edinburgh
|
||
|
! JCMB
|
||
|
! Mayfield Road
|
||
|
! Edinburgh EH9 3JZ
|
||
|
! United Kingdom
|
||
|
!
|
||
|
! This file is part of the software made public in
|
||
|
! Jenkins and Booth 2013 - arXiv:1306.XXXX
|
||
|
!
|
||
|
! The software computes the Panphasia Gaussian white noise field
|
||
|
! realisation described in detail in Jenkins 2013 - arXiv:1306.XXXX
|
||
|
!
|
||
|
!
|
||
|
!
|
||
|
! This software is free, subject to a agreeing licence conditions:
|
||
|
!
|
||
|
!
|
||
|
! (i) you will publish the phase descriptors and reference Jenkins (13)
|
||
|
! for any new simulations that use Panphasia phases. You will pass on this
|
||
|
! condition to others for any software or data you make available publically
|
||
|
! or privately that makes use of Panphasia.
|
||
|
!
|
||
|
! (ii) that you will ensure any publications using results derived from Panphasia
|
||
|
! will be submitted as a final version to arXiv prior to or coincident with
|
||
|
! publication in a journal.
|
||
|
!
|
||
|
!
|
||
|
! (iii) that you report any bugs in this software as soon as confirmed to
|
||
|
! A.R.Jenkins@durham.ac.uk
|
||
|
!
|
||
|
! (iv) that you understand that this software comes with no warranty and that is
|
||
|
! your responsibility to ensure that it is suitable for the purpose that
|
||
|
! you intend.
|
||
|
!
|
||
|
!=====================================================================================c
|
||
|
!{{{Rand_base (define kind types)
|
||
|
MODULE Rand_base
|
||
|
! This module just declares the base types
|
||
|
! we may have to edit this to match to the target machine
|
||
|
! we really need a power of 2 selected int kind in fortran-95 we could
|
||
|
! do this with a PURE function I think.
|
||
|
|
||
|
!
|
||
|
! 10 decimal digits will hold 2^31
|
||
|
!
|
||
|
|
||
|
INTEGER, PARAMETER :: Sint = SELECTED_INT_KIND(9)
|
||
|
! INTEGER, PARAMETER :: Sint = SELECTED_INT_KIND(10)
|
||
|
! INTEGER, PARAMETER :: Sint = 4
|
||
|
|
||
|
!
|
||
|
! 18-19 decimal digits will hold 2^63
|
||
|
! but all 19 digit numbers require 2^65 :-(
|
||
|
!
|
||
|
|
||
|
INTEGER, PARAMETER :: Dint = SELECTED_INT_KIND(17)
|
||
|
! INTEGER, PARAMETER :: Dint = SELECTED_INT_KIND(18)
|
||
|
! INTEGER, PARAMETER :: Dint = 8
|
||
|
|
||
|
! type for index counters must hold Nstore
|
||
|
INTEGER, PARAMETER :: Ctype = SELECTED_INT_KIND(3)
|
||
|
END MODULE Rand_base
|
||
|
!}}}
|
||
|
|
||
|
!{{{Rand_int (random integers mod 2^31-1)
|
||
|
|
||
|
MODULE Rand_int
|
||
|
USE Rand_base
|
||
|
IMPLICIT NONE
|
||
|
! The general approach of this module is two have
|
||
|
! two types Sint and Dint
|
||
|
!
|
||
|
! Sint should have at least 31 bits
|
||
|
! dint shouldhave at least 63
|
||
|
|
||
|
!{{{constants
|
||
|
|
||
|
INTEGER(KIND=Ctype), PARAMETER :: Nstate=5_Ctype
|
||
|
INTEGER(KIND=Ctype), PRIVATE, PARAMETER :: Nbatch=128_Ctype
|
||
|
INTEGER(KIND=Ctype), PRIVATE, PARAMETER :: Nstore=Nstate+Nbatch
|
||
|
|
||
|
INTEGER(KIND=Sint), PRIVATE, PARAMETER :: M = 2147483647_Sint
|
||
|
INTEGER(KIND=Dint), PRIVATE, PARAMETER :: Mask = 2147483647_Dint
|
||
|
INTEGER(KIND=Dint), PRIVATE, PARAMETER :: A1 = 107374182_Dint
|
||
|
INTEGER(KIND=Dint), PRIVATE, PARAMETER :: A5 = 104480_Dint
|
||
|
LOGICAL, PARAMETER :: Can_step_int=.TRUE.
|
||
|
LOGICAL, PARAMETER :: Can_reverse_int=.TRUE.
|
||
|
|
||
|
!}}}
|
||
|
|
||
|
!{{{Types
|
||
|
!
|
||
|
! This type holds the state of the generator
|
||
|
!
|
||
|
!{{{TYPE RAND_state
|
||
|
|
||
|
TYPE RAND_state
|
||
|
PRIVATE
|
||
|
INTEGER(KIND=Sint) :: state(Nstore)
|
||
|
! do we need to re-fill state table this is reset when we initialise state.
|
||
|
LOGICAL :: need_fill
|
||
|
! position of the next state variable to output
|
||
|
INTEGER(KIND=Ctype) :: pos
|
||
|
END TYPE RAND_state
|
||
|
|
||
|
!}}}
|
||
|
|
||
|
!
|
||
|
! This type defines the offset type used for stepping.
|
||
|
!
|
||
|
!{{{TYPE RAND_offset
|
||
|
|
||
|
TYPE RAND_offset
|
||
|
PRIVATE
|
||
|
INTEGER(KIND=Sint) :: poly(Nstate)
|
||
|
END TYPE RAND_offset
|
||
|
|
||
|
!}}}
|
||
|
|
||
|
!}}}
|
||
|
|
||
|
!{{{interface and overloads
|
||
|
!
|
||
|
! Allow automatic conversion between integers and offsets
|
||
|
!
|
||
|
INTERFACE ASSIGNMENT(=)
|
||
|
MODULE PROCEDURE Rand_set_offset
|
||
|
MODULE PROCEDURE Rand_load
|
||
|
MODULE PROCEDURE Rand_save
|
||
|
MODULE PROCEDURE Rand_seed
|
||
|
END INTERFACE
|
||
|
INTERFACE OPERATOR(+)
|
||
|
MODULE PROCEDURE Rand_add_offset
|
||
|
END INTERFACE
|
||
|
INTERFACE OPERATOR(*)
|
||
|
MODULE PROCEDURE Rand_mul_offset
|
||
|
END INTERFACE
|
||
|
|
||
|
!
|
||
|
! overload + as the boost/stepping operator
|
||
|
!
|
||
|
INTERFACE OPERATOR(+)
|
||
|
MODULE PROCEDURE Rand_step
|
||
|
MODULE PROCEDURE Rand_boost
|
||
|
END INTERFACE
|
||
|
!}}}
|
||
|
|
||
|
|
||
|
!{{{PUBLIC/PRIVATE
|
||
|
PRIVATE reduce,mod_saxpy,mod_sdot,p_saxpy,p_sdot,poly_mult
|
||
|
PRIVATE poly_square, poly_power
|
||
|
PRIVATE fill_state, repack_state
|
||
|
|
||
|
PUBLIC Rand_sint, Rand_sint_vec
|
||
|
|
||
|
PUBLIC Rand_save, Rand_load
|
||
|
PUBLIC Rand_set_offset, Rand_add_offset, Rand_mul_offset
|
||
|
PUBLIC Rand_step, Rand_boost, Rand_seed
|
||
|
!}}}
|
||
|
|
||
|
CONTAINS
|
||
|
!{{{Internals
|
||
|
!{{{RECURSIVE FUNCTION reduce(A)
|
||
|
RECURSIVE FUNCTION reduce(A)
|
||
|
!
|
||
|
! Take A Dint and reduce to Sint MOD M
|
||
|
!
|
||
|
INTEGER(KIND=Dint), INTENT(IN) :: A
|
||
|
INTEGER(KIND=Sint) reduce
|
||
|
INTEGER(KIND=Dint) tmp
|
||
|
|
||
|
tmp = A
|
||
|
DO WHILE( ISHFT(tmp, -31) .GT. 0 )
|
||
|
tmp = IAND(tmp,Mask) + ISHFT(tmp, -31)
|
||
|
END DO
|
||
|
IF( tmp .GE. M ) THEN
|
||
|
reduce = tmp - M
|
||
|
ELSE
|
||
|
reduce = tmp
|
||
|
END IF
|
||
|
END FUNCTION reduce
|
||
|
!}}}
|
||
|
!{{{RECURSIVE SUBROUTINE fill_state(x)
|
||
|
RECURSIVE SUBROUTINE fill_state(x)
|
||
|
TYPE(RAND_state), INTENT(INOUT) :: x
|
||
|
INTEGER(KIND=Ctype) i
|
||
|
INTRINSIC IAND, ISHFT
|
||
|
INTEGER(KIND=Dint) tmp
|
||
|
DO i=Nstate+1,Nstore
|
||
|
tmp = (x%state(i-5) * A5) + (x%state(i-1)*A1)
|
||
|
!
|
||
|
! now reduce down to mod M efficiently
|
||
|
! really hope the compiler in-lines this
|
||
|
!
|
||
|
! x%state(i) = reduce(tmp)
|
||
|
DO WHILE( ISHFT(tmp, -31) .GT. 0 )
|
||
|
tmp = IAND(tmp,Mask) + ISHFT(tmp, -31)
|
||
|
END DO
|
||
|
IF( tmp .GE. M ) THEN
|
||
|
x%state(i) = tmp - M
|
||
|
ELSE
|
||
|
x%state(i) = tmp
|
||
|
END IF
|
||
|
|
||
|
END DO
|
||
|
x%need_fill = .FALSE.
|
||
|
END SUBROUTINE fill_state
|
||
|
!}}}
|
||
|
!{{{RECURSIVE SUBROUTINE repack_state(x)
|
||
|
RECURSIVE SUBROUTINE repack_state(x)
|
||
|
TYPE(RAND_state), INTENT(INOUT) :: x
|
||
|
INTEGER(KIND=Ctype) i
|
||
|
DO i=1,Nstate
|
||
|
x%state(i) = x%state(i+x%pos-(Nstate+1))
|
||
|
END DO
|
||
|
x%pos = Nstate + 1
|
||
|
x%need_fill = .TRUE.
|
||
|
END SUBROUTINE repack_state
|
||
|
!}}}
|
||
|
!{{{RECURSIVE SUBROUTINE mod_saxpy(y,a,x)
|
||
|
RECURSIVE SUBROUTINE mod_saxpy(y,a,x)
|
||
|
INTEGER(KIND=Ctype) i
|
||
|
INTEGER(KIND=Sint) y(Nstate)
|
||
|
INTEGER(KIND=Sint) a
|
||
|
INTEGER(KIND=Sint) x(Nstate)
|
||
|
INTEGER(KIND=Dint) tx,ty,ta
|
||
|
|
||
|
IF( a .EQ. 0_Sint ) RETURN
|
||
|
|
||
|
! We use KIND=Dint temporaries here to ensure
|
||
|
! that we don't overflow in the expression
|
||
|
|
||
|
ta = a
|
||
|
DO i=1,Nstate
|
||
|
ty=y(i)
|
||
|
tx=x(i)
|
||
|
y(i) = reduce(ty + ta * tx)
|
||
|
END DO
|
||
|
|
||
|
END SUBROUTINE
|
||
|
!}}}
|
||
|
!{{{RECURSIVE SUBROUTINE mod_sdot(res,x,y)
|
||
|
RECURSIVE SUBROUTINE mod_sdot(res,x,y)
|
||
|
INTEGER(KIND=Sint), INTENT(OUT) :: res
|
||
|
INTEGER(KIND=Sint), INTENT(IN) :: x(Nstate) , y(Nstate)
|
||
|
INTEGER(KIND=Dint) dx, dy, dtmp
|
||
|
INTEGER(KIND=Sint) tmp
|
||
|
INTEGER(KIND=Ctype) i
|
||
|
|
||
|
tmp = 0
|
||
|
DO i=1,Nstate
|
||
|
dx = x(i)
|
||
|
dy = y(i)
|
||
|
dtmp = tmp
|
||
|
tmp = reduce(dtmp + dx * dy)
|
||
|
END DO
|
||
|
res = tmp
|
||
|
END SUBROUTINE
|
||
|
!}}}
|
||
|
!{{{RECURSIVE SUBROUTINE p_saxpy(y,a)
|
||
|
RECURSIVE SUBROUTINE p_saxpy(y,a)
|
||
|
! Calculates mod_saxpy(y,a,P)
|
||
|
INTEGER(KIND=Sint), INTENT(INOUT) :: y(Nstate)
|
||
|
INTEGER(KIND=Sint), INTENT(IN) :: a
|
||
|
INTEGER(KIND=Dint) tmp, dy, da
|
||
|
dy = y(1)
|
||
|
da = a
|
||
|
tmp = dy + da*A5
|
||
|
y(1) = reduce(tmp)
|
||
|
dy = y(5)
|
||
|
da = a
|
||
|
tmp = dy + da*A1
|
||
|
y(5) = reduce(tmp)
|
||
|
|
||
|
END SUBROUTINE
|
||
|
!}}}
|
||
|
!{{{RECURSIVE SUBROUTINE p_sdot(res,n,x)
|
||
|
RECURSIVE SUBROUTINE p_sdot(res,x)
|
||
|
INTEGER(KIND=Sint), INTENT(OUT) :: res
|
||
|
INTEGER(KIND=Sint), INTENT(IN) :: x(Nstate)
|
||
|
INTEGER(KIND=Dint) dx1, dx5, dtmp
|
||
|
dx1 = x(1)
|
||
|
dx5 = x(5)
|
||
|
|
||
|
dtmp = A1*dx5 + A5*dx1
|
||
|
res = reduce(dtmp)
|
||
|
END SUBROUTINE
|
||
|
!}}}
|
||
|
!{{{RECURSIVE SUBROUTINE poly_mult(a,b)
|
||
|
RECURSIVE SUBROUTINE poly_mult(a,b)
|
||
|
INTEGER(KIND=Sint), INTENT(INOUT) :: a(Nstate)
|
||
|
INTEGER(KIND=Sint), INTENT(IN) :: b(Nstate)
|
||
|
INTEGER(KIND=Sint) tmp((2*Nstate) - 1)
|
||
|
INTEGER(KIND=Ctype) i
|
||
|
|
||
|
tmp = 0_Sint
|
||
|
|
||
|
DO i=1,Nstate
|
||
|
CALL mod_saxpy(tmp(i:Nstate+i-1),a(i), b)
|
||
|
END DO
|
||
|
DO i=(2*Nstate)-1, Nstate+1, -1
|
||
|
CALL P_SAXPY(tmp(i-Nstate:i-1),tmp(i))
|
||
|
END DO
|
||
|
a = tmp(1:Nstate)
|
||
|
END SUBROUTINE
|
||
|
!}}}
|
||
|
!{{{RECURSIVE SUBROUTINE poly_square(a)
|
||
|
RECURSIVE SUBROUTINE poly_square(a)
|
||
|
INTEGER(KIND=Sint), INTENT(INOUT) :: a(Nstate)
|
||
|
INTEGER(KIND=Sint) tmp((2*Nstate) - 1)
|
||
|
INTEGER(KIND=Ctype) i
|
||
|
|
||
|
tmp = 0_Sint
|
||
|
|
||
|
DO i=1,Nstate
|
||
|
CALL mod_saxpy(tmp(i:Nstate+i-1),a(i), a)
|
||
|
END DO
|
||
|
DO i=(2*Nstate)-1, Nstate+1, -1
|
||
|
CALL P_SAXPY(tmp(i-Nstate:i-1),tmp(i))
|
||
|
END DO
|
||
|
a = tmp(1:Nstate)
|
||
|
END SUBROUTINE
|
||
|
!}}}
|
||
|
!{{{RECURSIVE SUBROUTINE poly_power(poly,n)
|
||
|
RECURSIVE SUBROUTINE poly_power(poly,n)
|
||
|
INTEGER(KIND=Sint), INTENT(INOUT) :: poly(Nstate)
|
||
|
INTEGER, INTENT(IN) :: n
|
||
|
INTEGER nn
|
||
|
INTEGER(KIND=Sint) x(Nstate), out(Nstate)
|
||
|
|
||
|
IF( n .EQ. 0 )THEN
|
||
|
poly = 0_Sint
|
||
|
poly(1) = 1_Sint
|
||
|
RETURN
|
||
|
ELSE IF( n .LT. 0 )THEN
|
||
|
poly = 0_Sint
|
||
|
RETURN
|
||
|
END IF
|
||
|
|
||
|
out = 0_sint
|
||
|
out(1) = 1_Sint
|
||
|
x = poly
|
||
|
nn = n
|
||
|
DO WHILE( nn .GT. 0 )
|
||
|
IF( MOD(nn,2) .EQ. 1 )THEN
|
||
|
call poly_mult(out,x)
|
||
|
END IF
|
||
|
nn = nn/2
|
||
|
IF( nn .GT. 0 )THEN
|
||
|
call poly_square(x)
|
||
|
END IF
|
||
|
END DO
|
||
|
poly = out
|
||
|
|
||
|
END SUBROUTINE poly_power
|
||
|
!}}}
|
||
|
!}}}
|
||
|
|
||
|
!{{{RECURSIVE SUBROUTINE Rand_seed( state, n )
|
||
|
RECURSIVE SUBROUTINE Rand_seed( state, n )
|
||
|
TYPE(Rand_state), INTENT(OUT) :: state
|
||
|
INTEGER, INTENT(IN) :: n
|
||
|
! initialise the genrator using a single integer
|
||
|
! fist initialise to an arbitrary state then boost by a multiple
|
||
|
! of a long distance
|
||
|
!
|
||
|
! state is moved forward by P^n steps
|
||
|
! we want this to be ok for seperating parallel sequences on MPP machines
|
||
|
! P is taken as a prime number as this should prevent strong correlations
|
||
|
! when the generators are operated in tight lockstep.
|
||
|
! equivalent points on different processors will also be related by a
|
||
|
! primative polynomial
|
||
|
! P is 2^48-59
|
||
|
TYPE(Rand_state) tmp
|
||
|
TYPE(Rand_offset), PARAMETER :: P = &
|
||
|
Rand_offset( (/ 1509238949_Sint ,2146167999_Sint ,1539340803_Sint , &
|
||
|
1041407428_Sint ,666274987_Sint /) )
|
||
|
|
||
|
CALL Rand_load( tmp, (/ 5, 4, 3, 2, 1 /) )
|
||
|
state = Rand_boost( tmp, Rand_mul_offset(P, n ))
|
||
|
|
||
|
END SUBROUTINE Rand_seed
|
||
|
!}}}
|
||
|
!{{{RECURSIVE SUBROUTINE Rand_load( state, input )
|
||
|
RECURSIVE SUBROUTINE Rand_load( state, input )
|
||
|
TYPE(RAND_state), INTENT(OUT) :: state
|
||
|
INTEGER, INTENT(IN) :: input(Nstate)
|
||
|
|
||
|
INTEGER(KIND=Ctype) i
|
||
|
|
||
|
state%state = 0_Sint
|
||
|
DO i=1,Nstate
|
||
|
state%state(i) = MOD(INT(input(i),KIND=Sint),M)
|
||
|
END DO
|
||
|
state%need_fill = .TRUE.
|
||
|
state%pos = Nstate + 1
|
||
|
END SUBROUTINE Rand_load
|
||
|
!}}}
|
||
|
!{{{RECURSIVE SUBROUTINE Rand_save( save_vec,state )
|
||
|
RECURSIVE SUBROUTINE Rand_save( save_vec, x )
|
||
|
INTEGER, INTENT(OUT) :: save_vec(Nstate)
|
||
|
TYPE(RAND_state), INTENT(IN) :: x
|
||
|
|
||
|
INTEGER(KIND=Ctype) i
|
||
|
DO i=1,Nstate
|
||
|
save_vec(i) = x%state(x%pos-(Nstate+1) + i)
|
||
|
END DO
|
||
|
END SUBROUTINE Rand_save
|
||
|
!}}}
|
||
|
|
||
|
!{{{RECURSIVE SUBROUTINE Rand_set_offset( offset, n )
|
||
|
RECURSIVE SUBROUTINE Rand_set_offset( offset, n )
|
||
|
TYPE(Rand_offset), INTENT(OUT) :: offset
|
||
|
INTEGER, INTENT(IN) :: n
|
||
|
|
||
|
offset%poly = 0_Sint
|
||
|
IF ( n .GE. 0 ) THEN
|
||
|
offset%poly(2) = 1_Sint
|
||
|
call poly_power(offset%poly,n)
|
||
|
ELSE
|
||
|
!
|
||
|
! This is X^-1
|
||
|
!
|
||
|
offset%poly(4) = 858869107_Sint
|
||
|
offset%poly(5) = 1840344978_Sint
|
||
|
call poly_power(offset%poly,-n)
|
||
|
END IF
|
||
|
END SUBROUTINE Rand_set_offset
|
||
|
!}}}
|
||
|
!{{{TYPE(Rand_offset) RECURSIVE FUNCTION Rand_add_offset( a, b )
|
||
|
TYPE(Rand_offset) RECURSIVE FUNCTION Rand_add_offset( a, b )
|
||
|
TYPE(Rand_offset), INTENT(IN) :: a, b
|
||
|
|
||
|
Rand_add_offset = a
|
||
|
CALL poly_mult(Rand_add_offset%poly,b%poly)
|
||
|
RETURN
|
||
|
END FUNCTION Rand_add_offset
|
||
|
!}}}
|
||
|
!{{{TYPE(Rand_offset) RECURSIVE FUNCTION Rand_mul_offset( a, n )
|
||
|
TYPE(Rand_offset) RECURSIVE FUNCTION Rand_mul_offset( a, n )
|
||
|
TYPE(Rand_offset), INTENT(IN) :: a
|
||
|
INTEGER, INTENT(IN) :: n
|
||
|
Rand_mul_offset = a
|
||
|
CALL poly_power(Rand_mul_offset%poly,n)
|
||
|
RETURN
|
||
|
END FUNCTION Rand_mul_offset
|
||
|
!}}}
|
||
|
!{{{RECURSIVE FUNCTION Rand_boost(x, offset)
|
||
|
RECURSIVE FUNCTION Rand_boost(x, offset)
|
||
|
TYPE(Rand_state) Rand_boost
|
||
|
TYPE(Rand_state), INTENT(IN) :: x
|
||
|
TYPE(Rand_offset), INTENT(IN) :: offset
|
||
|
INTEGER(KIND=Sint) tmp(2*Nstate-1), res(Nstate)
|
||
|
INTEGER(KIND=Ctype) i
|
||
|
|
||
|
DO i=1,Nstate
|
||
|
tmp(i) = x%state(x%pos-(Nstate+1) + i)
|
||
|
END DO
|
||
|
tmp(Nstate+1:) = 0_Sint
|
||
|
|
||
|
DO i=1,Nstate-1
|
||
|
call P_SDOT(tmp(i+Nstate),tmp(i:Nstate+i-1))
|
||
|
END DO
|
||
|
|
||
|
DO i=1,Nstate
|
||
|
call mod_sdot(res(i),offset%poly,tmp(i:Nstate+i-1))
|
||
|
END DO
|
||
|
Rand_boost%state = 0_Sint
|
||
|
DO i=1,Nstate
|
||
|
Rand_boost%state(i) = res(i)
|
||
|
END DO
|
||
|
Rand_boost%need_fill = .TRUE.
|
||
|
Rand_boost%pos = Nstate + 1
|
||
|
|
||
|
END FUNCTION Rand_boost
|
||
|
!}}}
|
||
|
!{{{RECURSIVE FUNCTION Rand_step(x, n)
|
||
|
RECURSIVE FUNCTION Rand_step(x, n)
|
||
|
TYPE(Rand_state) Rand_step
|
||
|
TYPE(RAND_state), INTENT(IN) :: x
|
||
|
INTEGER, INTENT(IN) :: n
|
||
|
TYPE(Rand_offset) tmp
|
||
|
|
||
|
CALL Rand_set_offset(tmp,n)
|
||
|
Rand_step=Rand_boost(x,tmp)
|
||
|
|
||
|
END FUNCTION
|
||
|
!}}}
|
||
|
|
||
|
!{{{RECURSIVE FUNCTION Rand_sint(x)
|
||
|
RECURSIVE FUNCTION Rand_sint(x)
|
||
|
TYPE(RAND_state), INTENT(INOUT) :: x
|
||
|
INTEGER(KIND=Sint) Rand_sint
|
||
|
IF( x%pos .GT. Nstore )THEN
|
||
|
CALL repack_state(x)
|
||
|
END IF
|
||
|
IF( x%need_fill ) CALL fill_state(x)
|
||
|
Rand_sint = x%state(x%pos)
|
||
|
x%pos = x%pos + 1
|
||
|
RETURN
|
||
|
END FUNCTION Rand_sint
|
||
|
!}}}
|
||
|
!{{{RECURSIVE SUBROUTINE Rand_sint_vec(iv,x)
|
||
|
RECURSIVE SUBROUTINE Rand_sint_vec(iv,x)
|
||
|
INTEGER(KIND=Sint), INTENT(OUT) :: iv(:)
|
||
|
TYPE(RAND_state), INTENT(INOUT) :: x
|
||
|
INTEGER left,start, chunk, i
|
||
|
|
||
|
start=1
|
||
|
left=SIZE(iv)
|
||
|
DO WHILE( left .GT. 0 )
|
||
|
IF( x%pos .GT. Nstore )THEN
|
||
|
CALL repack_state(x)
|
||
|
END IF
|
||
|
IF( x%need_fill ) CALL fill_state(x)
|
||
|
|
||
|
chunk = MIN(left,Nstore-x%pos+1)
|
||
|
DO i=0,chunk-1
|
||
|
iv(start+i) = x%state(x%pos+i)
|
||
|
END DO
|
||
|
start = start + chunk
|
||
|
x%pos = x%pos + chunk
|
||
|
left = left - chunk
|
||
|
END DO
|
||
|
|
||
|
RETURN
|
||
|
END SUBROUTINE Rand_sint_vec
|
||
|
!}}}
|
||
|
|
||
|
|
||
|
END MODULE Rand_int
|
||
|
|
||
|
!}}}
|
||
|
|
||
|
!{{{Rand (use Rand_int to make random reals)
|
||
|
|
||
|
MODULE Rand
|
||
|
USE Rand_int
|
||
|
IMPLICIT NONE
|
||
|
|
||
|
!{{{Parameters
|
||
|
|
||
|
INTEGER, PARAMETER :: RAND_kind1 = SELECTED_REAL_KIND(10)
|
||
|
INTEGER, PARAMETER :: RAND_kind2 = SELECTED_REAL_KIND(6)
|
||
|
|
||
|
INTEGER, PARAMETER, PRIVATE :: Max_block=100
|
||
|
INTEGER(KIND=Sint), PRIVATE, PARAMETER :: M = 2147483647
|
||
|
REAL(KIND=RAND_kind1), PRIVATE, PARAMETER :: INVMP1_1 = ( 1.0_RAND_kind1 / 2147483647.0_RAND_kind1 )
|
||
|
REAL(KIND=RAND_kind2), PRIVATE, PARAMETER :: INVMP1_2 = ( 1.0_RAND_kind2 / 2147483647.0_RAND_kind2 )
|
||
|
|
||
|
LOGICAL, PARAMETER :: Can_step = Can_step_int
|
||
|
LOGICAL, PARAMETER :: Can_reverse = Can_reverse_int
|
||
|
|
||
|
!}}}
|
||
|
PUBLIC Rand_real
|
||
|
|
||
|
|
||
|
INTERFACE Rand_real
|
||
|
MODULE PROCEDURE Rand_real1
|
||
|
MODULE PROCEDURE Rand_real2
|
||
|
MODULE PROCEDURE Rand_real_vec1
|
||
|
MODULE PROCEDURE Rand_real_vec2
|
||
|
END INTERFACE
|
||
|
|
||
|
|
||
|
CONTAINS
|
||
|
|
||
|
!{{{RECURSIVE SUBROUTINE Rand_real1(y,x)
|
||
|
RECURSIVE SUBROUTINE Rand_real1(y,x)
|
||
|
REAL(KIND=RAND_kind1), INTENT(OUT) :: y
|
||
|
TYPE(RAND_state), INTENT(INOUT) :: x
|
||
|
INTEGER(KIND=Sint) Z
|
||
|
|
||
|
Z = Rand_sint(x)
|
||
|
IF (Z .EQ. 0) Z = M
|
||
|
|
||
|
y = ((Z-0.5d0)*INVMP1_1)
|
||
|
RETURN
|
||
|
END SUBROUTINE Rand_real1
|
||
|
!}}}
|
||
|
!{{{RECURSIVE SUBROUTINE Rand_real2(y,x)
|
||
|
RECURSIVE SUBROUTINE Rand_real2(y,x)
|
||
|
REAL(KIND=RAND_kind2), INTENT(OUT) :: y
|
||
|
TYPE(RAND_state), INTENT(INOUT) :: x
|
||
|
INTEGER(KIND=Sint) Z
|
||
|
|
||
|
Z = Rand_sint(x)
|
||
|
IF (Z .EQ. 0) Z = M
|
||
|
|
||
|
y = ((Z-0.5d0)*INVMP1_1) ! generate in double and truncate.
|
||
|
RETURN
|
||
|
END SUBROUTINE Rand_real2
|
||
|
!}}}
|
||
|
|
||
|
!{{{RECURSIVE SUBROUTINE Rand_real_vec1(rv,x)
|
||
|
RECURSIVE SUBROUTINE Rand_real_vec1(rv,x)
|
||
|
TYPE(RAND_state), INTENT(INOUT) :: x
|
||
|
REAL(KIND=RAND_kind1) rv(:)
|
||
|
INTEGER left,start, chunk, i
|
||
|
INTEGER(KIND=Sint) Z
|
||
|
INTEGER(KIND=Sint) temp(MIN(SIZE(rv),Max_block))
|
||
|
|
||
|
start=0
|
||
|
left=SIZE(rv)
|
||
|
DO WHILE( left .GT. 0 )
|
||
|
chunk = MIN(left,Max_block)
|
||
|
CALL Rand_sint_vec(temp(1:chunk),x)
|
||
|
DO i=1,chunk
|
||
|
Z = temp(i)
|
||
|
IF (Z .EQ. 0) Z = M
|
||
|
rv(start+i) = (Z-0.5d0)*INVMP1_1
|
||
|
END DO
|
||
|
start = start + chunk
|
||
|
left = left - chunk
|
||
|
END DO
|
||
|
|
||
|
RETURN
|
||
|
END SUBROUTINE Rand_real_vec1
|
||
|
!}}}
|
||
|
!{{{RECURSIVE SUBROUTINE Rand_real_vec2(rv,x)
|
||
|
RECURSIVE SUBROUTINE Rand_real_vec2(rv,x)
|
||
|
TYPE(RAND_state), INTENT(INOUT) :: x
|
||
|
REAL(KIND=RAND_kind2) rv(:)
|
||
|
INTEGER left,start, chunk, i
|
||
|
INTEGER(KIND=Sint) Z
|
||
|
INTEGER(KIND=Sint) temp(MIN(SIZE(rv),Max_block))
|
||
|
|
||
|
start=0
|
||
|
left=SIZE(rv)
|
||
|
DO WHILE( left .GT. 0 )
|
||
|
chunk = MIN(left,Max_block)
|
||
|
CALL Rand_sint_vec(temp(1:chunk),x)
|
||
|
DO i=1,chunk
|
||
|
Z = temp(i)
|
||
|
IF (Z .EQ. 0) Z = M
|
||
|
rv(start+i) = (Z-0.5d0)*INVMP1_2
|
||
|
END DO
|
||
|
start = start + chunk
|
||
|
left = left - chunk
|
||
|
END DO
|
||
|
|
||
|
RETURN
|
||
|
END SUBROUTINE Rand_real_vec2
|
||
|
!}}}
|
||
|
END MODULE Rand
|
||
|
|
||
|
!}}}
|
||
|
|
||
|
!{{{test program
|
||
|
! PROGRAM test_random
|
||
|
! use Rand
|
||
|
! TYPE(RAND_state) x
|
||
|
! REAL y
|
||
|
! CALL Rand_load(x,(/5,4,3,2,1/))
|
||
|
! DO I=0,10
|
||
|
! CALL Rand_real(y,x)
|
||
|
! WRITE(*,10) I,y
|
||
|
! END DO
|
||
|
!
|
||
|
!10 FORMAT(I10,E25.16)
|
||
|
!
|
||
|
! END
|
||
|
|
||
|
! 0 0.5024326127022505E-01
|
||
|
! 1 0.8260946767404675E-01
|
||
|
! 2 0.2123264316469431E-01
|
||
|
! 3 0.6926658791489899E+00
|
||
|
! 4 0.2076155943796039E+00
|
||
|
! 5 0.4327449947595596E-01
|
||
|
! 6 0.2204052871093154E-01
|
||
|
! 7 0.1288446951657534E+00
|
||
|
! 8 0.4859915426932275E+00
|
||
|
! 9 0.5721384193748236E-01
|
||
|
! 10 0.7996825082227588E+00
|
||
|
!
|
||
|
|
||
|
|
||
|
!}}}
|
||
|
|