| 1 |
chuckv |
4 |
module sprng_mod |
| 2 |
|
|
use definitions, ONLY : DP |
| 3 |
|
|
use parameter, ONLY : iseed |
| 4 |
|
|
#ifdef MPI |
| 5 |
|
|
use mpi_module |
| 6 |
|
|
#endif |
| 7 |
|
|
implicit none |
| 8 |
|
|
PRIVATE |
| 9 |
|
|
|
| 10 |
|
|
#ifdef MPI |
| 11 |
|
|
#define USE_MPI |
| 12 |
|
|
#endif |
| 13 |
|
|
#define CHECK_POINTERS |
| 14 |
|
|
|
| 15 |
|
|
|
| 16 |
|
|
#include "headers/sprng_f.h" |
| 17 |
|
|
|
| 18 |
|
|
|
| 19 |
|
|
public :: random_gauss |
| 20 |
|
|
public :: init_gauss_stream |
| 21 |
|
|
public :: gran |
| 22 |
|
|
public :: init_random_stream |
| 23 |
|
|
public :: get_random_sprng |
| 24 |
|
|
|
| 25 |
|
|
integer, parameter :: gtype = 3 ! sprng random generator type |
| 26 |
|
|
! Combined Multiple Recursive generator |
| 27 |
|
|
logical :: is_init_gauss = .false. |
| 28 |
|
|
logical :: is_init_ran = .false. |
| 29 |
|
|
#ifdef mpi |
| 30 |
|
|
integer, parameter :: gauss_stream_number = node + nprocs |
| 31 |
|
|
integer, parameter :: ran_stream_number = node |
| 32 |
|
|
integer, parameter :: n_streams = 2*nprocs |
| 33 |
|
|
#else |
| 34 |
|
|
integer, parameter :: gauss_stream_number = 1 |
| 35 |
|
|
integer, parameter :: ran_stream_number = 0 |
| 36 |
|
|
integer, parameter :: n_streams = 2 |
| 37 |
|
|
#endif |
| 38 |
|
|
|
| 39 |
|
|
|
| 40 |
|
|
SPRNG_POINTER gauss_stream |
| 41 |
|
|
SPRNG_POINTER ran_stream |
| 42 |
|
|
contains |
| 43 |
|
|
|
| 44 |
|
|
subroutine init_random_stream () |
| 45 |
|
|
integer :: junk |
| 46 |
|
|
if (is_init_ran) return |
| 47 |
|
|
ran_stream = init_sprng(gtype,ran_stream_number, & |
| 48 |
|
|
n_streams,iseed,SPRNG_DEFAULT) |
| 49 |
|
|
is_init_ran = .true. |
| 50 |
|
|
! junk = print_sprng(ran_stream) |
| 51 |
|
|
end subroutine init_random_stream |
| 52 |
|
|
|
| 53 |
|
|
subroutine init_gauss_stream () |
| 54 |
|
|
integer :: junk |
| 55 |
|
|
! . if stream is already initialized then return |
| 56 |
|
|
if (is_init_gauss) return |
| 57 |
|
|
|
| 58 |
|
|
gauss_stream = init_sprng(gtype,gauss_stream_number, & |
| 59 |
|
|
n_streams,iseed,SPRNG_DEFAULT) |
| 60 |
|
|
is_init_gauss = .true. |
| 61 |
|
|
|
| 62 |
|
|
! junk = print_sprng(gauss_stream) |
| 63 |
|
|
end subroutine init_gauss_stream |
| 64 |
|
|
|
| 65 |
|
|
function get_random_sprng () result(ran_n) |
| 66 |
|
|
real( kind = DP ) :: ran_n |
| 67 |
|
|
if(.not. is_init_ran) then |
| 68 |
|
|
call init_random_stream() |
| 69 |
|
|
endif |
| 70 |
|
|
|
| 71 |
|
|
ran_n = sprng(ran_stream) |
| 72 |
|
|
return |
| 73 |
|
|
end function get_random_sprng |
| 74 |
|
|
|
| 75 |
|
|
function random_gauss () result(X) |
| 76 |
|
|
use constants, ONLY : twopi |
| 77 |
|
|
real( kind = dp ) :: X,Y |
| 78 |
|
|
real( kind = dp ) :: R1,R2 |
| 79 |
|
|
real( kind = dp ) :: rn1,rn2 |
| 80 |
|
|
! . if the gauss stream has not been initialized then init |
| 81 |
|
|
if(.not. is_init_gauss) then |
| 82 |
|
|
call init_gauss_stream() |
| 83 |
|
|
endif |
| 84 |
|
|
|
| 85 |
|
|
|
| 86 |
|
|
rn1 = sprng(gauss_stream) |
| 87 |
|
|
rn2 = sprng(gauss_stream) |
| 88 |
|
|
|
| 89 |
|
|
|
| 90 |
|
|
R1 = sqrt(-2.0_DP * log(rn1)) |
| 91 |
|
|
R2 = twopi * rn2 |
| 92 |
|
|
|
| 93 |
|
|
X = R1*cos(R2) |
| 94 |
|
|
|
| 95 |
|
|
end function random_gauss |
| 96 |
|
|
|
| 97 |
|
|
!*********************************************************************** |
| 98 |
|
|
! qran is a (quick) normal, linear congruential random number |
| 99 |
|
|
! generator with "well-chosen" constants taken from Numerical |
| 100 |
|
|
! Recipes. Call it with a positive integer seed, it'll change that |
| 101 |
|
|
! seed and return a number in [0,1]...sjs 3/18/92 |
| 102 |
|
|
!*********************************************************************** |
| 103 |
|
|
! Note: Numerical Recipes doesn't SAVE the variables, which it |
| 104 |
|
|
! should. |
| 105 |
|
|
!*********************************************************************** |
| 106 |
|
|
function qran(iseed) result(random) |
| 107 |
|
|
! |
| 108 |
|
|
integer, parameter :: imod = 134456 |
| 109 |
|
|
integer, parameter :: imult = 3887 |
| 110 |
|
|
integer, parameter :: iadd = 29573 |
| 111 |
|
|
integer :: iseed |
| 112 |
|
|
|
| 113 |
|
|
real( kind = DP ) :: random |
| 114 |
|
|
! |
| 115 |
|
|
save |
| 116 |
|
|
! |
| 117 |
|
|
iseed = mod(imult * iseed + iadd, imod) |
| 118 |
|
|
random = dble(iseed) / dble(imod) |
| 119 |
|
|
return |
| 120 |
|
|
end function qran |
| 121 |
|
|
!*********************************************************************** |
| 122 |
|
|
! gran is a (good) more reliable random number generator based on 3 |
| 123 |
|
|
! linear congruential cycles, with shuffling. It was taken straight |
| 124 |
|
|
! from Numerical Recipes....sjs 3/18/92 |
| 125 |
|
|
!*********************************************************************** |
| 126 |
|
|
! Note: Numerical Recipes doesn't SAVE the variables, which it |
| 127 |
|
|
! should. |
| 128 |
|
|
!*********************************************************************** |
| 129 |
|
|
function gran(iseed) result(random) |
| 130 |
|
|
|
| 131 |
|
|
logical init |
| 132 |
|
|
integer isize, imod1, imult1, iadd1 |
| 133 |
|
|
integer imod2, imult2, iadd2 |
| 134 |
|
|
integer imod3, imult3, iadd3 |
| 135 |
|
|
integer iseed, iseed1, iseed2, iseed3, islot |
| 136 |
|
|
real( kind = DP ) :: div1, div2 |
| 137 |
|
|
real( kind = DP ) :: random |
| 138 |
|
|
! |
| 139 |
|
|
parameter (isize = 97) |
| 140 |
|
|
parameter (imod1 = 259200, imult1 = 7141, iadd1 = 54773, & |
| 141 |
|
|
div1 = 1.E0_DP / imod1) |
| 142 |
|
|
parameter (imod2 = 134456, imult2 = 8121, iadd2 = 28411, & |
| 143 |
|
|
div2 = 1.E0_DP / imod2) |
| 144 |
|
|
parameter (imod3 = 243000, imult3 = 4561, iadd3 = 51349) |
| 145 |
|
|
! |
| 146 |
|
|
save |
| 147 |
|
|
! |
| 148 |
|
|
real( kind = DP ) :: store(isize) |
| 149 |
|
|
! |
| 150 |
|
|
data init/.true./ |
| 151 |
|
|
! |
| 152 |
|
|
if (iseed .lt. 0 .or. init) then |
| 153 |
|
|
init = .false. |
| 154 |
|
|
iseed1 = mod(iadd1 - iseed, imod1) |
| 155 |
|
|
iseed1 = mod(imult1 * iseed1 + iadd1, imod1) |
| 156 |
|
|
iseed2 = mod(iseed1, imod2) |
| 157 |
|
|
iseed1 = mod(imult1 * iseed1 + iadd1, imod1) |
| 158 |
|
|
iseed3 = mod(iseed1, imod3) |
| 159 |
|
|
do islot = 1, isize |
| 160 |
|
|
iseed1 = mod(imult1 * iseed1 + iadd1, imod1) |
| 161 |
|
|
iseed2 = mod(imult2 * iseed2 + iadd2, imod2) |
| 162 |
|
|
store(islot) = (dble(iseed1) + dble(iseed2) * div2) * div1 |
| 163 |
|
|
end do |
| 164 |
|
|
! iseed = 1 |
| 165 |
|
|
end if |
| 166 |
|
|
iseed1 = mod(imult1 * iseed1 + iadd1, imod1) |
| 167 |
|
|
iseed2 = mod(imult2 * iseed2 + iadd2, imod2) |
| 168 |
|
|
iseed3 = mod(imult3 * iseed3 + iadd3, imod3) |
| 169 |
|
|
islot = 1 + (isize * iseed3) / imod3 |
| 170 |
|
|
random = store(islot) |
| 171 |
|
|
store(islot) = (dble(iseed1) + dble(iseed2) * div2) * div1 |
| 172 |
|
|
return |
| 173 |
|
|
end function gran |
| 174 |
|
|
|
| 175 |
|
|
|
| 176 |
|
|
|
| 177 |
|
|
end module sprng_mod |
