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 |