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 |