| 1 |
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 |
