mirror/music-panphasia
1
0
Fork 0
mirror of https://github.com/glatterf42/music-panphasia.git synced 2024-09-19 16:13:46 +02:00
Code Releases Activity
music-panphasia/panphasia/generic_lecuyer.f90
glatterf42 1e32e12349 First commit of original files.
2022-04-29 14:37:23 +02:00

683 lines
17 KiB
Fortran
Raw Permalink Blame History

!=====================================================================================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
!
!}}}
View git blame Copy permalink