Go to the documentation of this file.
38 subroutine rand_check_gsl(code, value)
43 integer(kind=c_int) :: value
48 call die_msg(code,
"GSL RNG initialization failed")
50 call die_msg(code,
"GSL RNG has not been successfully initialized")
52 call die_msg(code,
"GSL RNG was already initialized")
54 call die_msg(code,
"Unknown GSL RNG interface return value: " &
58 end subroutine rand_check_gsl
69 integer,
intent(in) :: seed
71 integer,
intent(in) :: offset
73 integer :: i, n, clock
74 integer,
allocatable :: seed_vec(:)
76 integer(kind=c_int) :: c_clock
80 #ifndef DOXYGEN_SKIP_DOC
84 integer(kind=c_int),
value :: seed
92 call system_clock(count = clock)
100 clock = clock + 67 * offset
102 c_clock = int(clock, kind=c_int)
105 call random_seed(
size = n)
106 allocate(seed_vec(n))
108 seed_vec = clock + 37 * (/ (i - 1, i = 1, n) /)
109 call random_seed(put = seed_vec)
122 #ifndef DOXYGEN_SKIP_DOC
141 real(kind=c_double),
target :: rnd
142 type(c_ptr) :: rnd_ptr
148 #ifndef DOXYGEN_SKIP_DOC
150 integer(kind=c_int) function pmc_rand_gsl(harvest)
bind(c)
152 type(c_ptr),
value :: harvest
163 call random_number(rnd)
175 integer,
intent(in) :: n
178 integer(kind=c_int) :: n_c
179 integer(kind=c_int),
target :: harvest
180 type(c_ptr) :: harvest_ptr
184 #ifndef DOXYGEN_SKIP_DOC
188 integer(kind=c_int),
value :: n
189 type(c_ptr),
value :: harvest
195 call assert(669532625, n >= 1)
197 n_c = int(n, kind=c_int)
198 harvest_ptr = c_loc(harvest)
217 real(kind=
dp),
intent(in) :: val
221 if (
pmc_random() < real(ceiling(val), kind=
dp) - val)
then
255 real(kind=
dp),
intent(in) :: mean
258 real(kind=c_double) :: mean_c
259 integer(kind=c_int),
target :: harvest
260 type(c_ptr) :: harvest_ptr
262 real(kind=
dp) :: l, p
267 #ifndef DOXYGEN_SKIP_DOC
272 real(kind=c_double),
value :: mean
273 type(c_ptr),
value :: harvest
279 call assert(368397056, mean >= 0d0)
281 mean_c = real(mean, kind=c_double)
282 harvest_ptr = c_loc(harvest)
283 call rand_check_gsl(353483140, &
287 if (mean <= 10d0)
then
322 integer,
intent(in) :: n
324 real(kind=
dp),
intent(in) :: p
327 integer(kind=c_int) :: n_c
328 real(kind=c_double) :: p_c
329 integer(kind=c_int),
target :: harvest
330 type(c_ptr) :: harvest_ptr
332 real(kind=
dp) :: np, nomp, q, g_real
337 #ifndef DOXYGEN_SKIP_DOC
342 integer(kind=c_int),
value :: n
343 real(kind=c_double),
value :: p
344 type(c_ptr),
value :: harvest
350 call assert(130699849, n >= 0)
351 call assert(754379796, p >= 0d0)
352 call assert(678506029, p <= 1d0)
354 n_c = int(n, kind=c_int)
355 p_c = real(p, kind=c_double)
356 harvest_ptr = c_loc(harvest)
357 call rand_check_gsl(208869397, &
361 np = real(n, kind=
dp) * p
362 nomp = real(n, kind=
dp) * (1d0 - p)
363 if ((np >= 10d0) .and. (nomp >= 10d0))
then
367 elseif (np < 1d-200)
then
369 elseif (nomp < 1d-200)
then
386 if (g_real > real(n - sum, kind=
dp))
exit
410 real(kind=
dp),
intent(in) :: mean
412 real(kind=
dp),
intent(in) :: stddev
415 real(kind=c_double) :: mean_c, stddev_c
416 real(kind=c_double),
target :: harvest
417 type(c_ptr) :: harvest_ptr
419 real(kind=
dp) :: u1, u2, r, theta, z0, z1
423 #ifndef DOXYGEN_SKIP_DOC
428 real(kind=c_double),
value :: mean
429 real(kind=c_double),
value :: stddev
430 type(c_ptr),
value :: harvest
436 call assert(898978929, stddev >= 0d0)
438 mean_c = real(mean, kind=c_double)
439 stddev_c = real(stddev, kind=c_double)
440 harvest_ptr = c_loc(harvest)
441 call rand_check_gsl(102078576, &
449 r = sqrt(-2d0 * log(u1))
450 theta = 2d0 *
const%pi * u2
469 real(kind=
dp),
intent(in) :: mean(:)
471 real(kind=
dp),
intent(in) :: stddev(
size(mean))
473 real(kind=
dp),
intent(out) :: val(
size(mean))
492 real(kind=
dp),
intent(in) :: pdf(:)
494 real(kind=
dp) :: pdf_max
499 pdf_max = maxval(pdf)
500 if (minval(pdf) < 0d0)
then
501 call die_msg(121838078,
'pdf contains negative values')
503 if (pdf_max <= 0d0)
then
504 call die_msg(119208863,
'pdf is not positive')
507 do while (.not. found)
526 integer,
intent(in) :: pdf(:)
528 integer :: pdf_max, k
532 pdf_max = maxval(pdf)
533 if (minval(pdf) < 0)
then
534 call die_msg(598024763,
'pdf contains negative values')
536 if (pdf_max <= 0)
then
537 call die_msg(109961454,
'pdf is not positive')
540 do while (.not. found)
542 if (
pmc_random() < real(pdf(k), kind=
dp) / real(pdf_max, kind=
dp))
then
560 real(kind=
dp),
intent(in) :: vec_cts(:)
562 integer,
intent(in) :: n_samp
564 integer,
intent(out) :: vec_disc(size(vec_cts))
568 call assert(617770167, n_samp >= 0)
572 vec_disc(k) = vec_disc(k) + 1
602 character(len=PMC_UUID_LEN),
intent(out) :: uuid
630 uuid(20:20) = achar(iachar(
'8') + i - 1)
632 uuid(20:20) = achar(iachar(
'A') + i - 3)
integer function prob_round(val)
Round val to floor(val) or ceiling(val) with probability proportional to the relative distance from v...
integer, parameter pmc_rand_gsl_init_fail
Result code indicating initialization failure.
Wrapper functions for MPI.
int pmc_srand_gsl(int seed)
Initialize the random number generator with the given seed.
subroutine sample_vec_cts_to_disc(vec_cts, n_samp, vec_disc)
Convert a real-valued vector into an integer-valued vector by sampling.
integer function rand_poisson(mean)
Generate a Poisson-distributed random number with the given mean.
real(kind=dp) function rand_normal(mean, stddev)
Generates a normally distributed random number with the given mean and standard deviation.
subroutine die_msg(code, error_msg)
Error immediately.
integer function pmc_mpi_rank()
Returns the rank of the current process.
integer, parameter dp
Kind of a double precision real number.
integer, parameter pmc_rand_gsl_not_init
Result code indicating the generator was not initialized when it should have been.
subroutine assert(code, condition_ok)
Errors unless condition_ok is true.
int pmc_rand_poisson_gsl(double mean, int *harvest)
Generate a Poisson-distributed random integer.
integer, parameter pmc_rand_gsl_already_init
Result code indicating the generator was already initialized when an initialization was attempted.
integer function pmc_rand_int(n)
Returns a random integer between 1 and n.
subroutine pmc_srand(seed, offset)
Initializes the random number generator to the state defined by the given seed plus offset....
subroutine pmc_mpi_bcast_integer(val)
Broadcast the given value from process 0 to all other processes.
character(len=pmc_util_convert_string_len) function integer_to_string(val)
Convert an integer to a string format.
type(const_t), save const
Fixed variable for accessing the constant's values.
integer, parameter pmc_rand_gsl_success
Result code indicating successful completion.
integer function sample_cts_pdf(pdf)
Sample the given continuous probability density function.
int pmc_rand_finalize_gsl()
Cleanup and deallocate the random number generator.
int pmc_rand_int_gsl(int n, int *harvest)
Generate a uniform random integer in .
Random number generators.
subroutine rand_normal_array_1d(mean, stddev, val)
Generates a vector of normally distributed random numbers with the given means and standard deviation...
subroutine uuid4_str(uuid)
Generate a version 4 UUID as a string.
integer, parameter pmc_uuid_len
Length of a UUID string.
Common utility subroutines.
int pmc_rand_gsl(double *harvest)
Generate a uniform random number in .
subroutine pmc_rand_finalize()
Cleanup the random number generator.
integer function rand_binomial(n, p)
Generate a Binomial-distributed random number with the given parameters.
character function rand_hex_char()
Generate a random hexadecimal character.
integer function sample_disc_pdf(pdf)
Sample the given discrete probability density function.
real(kind=dp) function pmc_random()
Returns a random number between 0 and 1.
int pmc_rand_binomial_gsl(int n, double p, int *harvest)
Generate a Binomial-distributed random integer.
int pmc_rand_normal_gsl(double mean, double stddev, double *harvest)
Generate a normally-distributed random number.