1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
|
/*************************************************************************
* c_gsl_wrap - C front ends for bits of GNU Scientific Library *
* *
* Copyright (C) 2016, 2017 by Boud Roukema *
* boud cosmo.torun.pl *
* *
* This program is free software; you can redistribute it and/or modify*
* it under the terms of the GNU General Public License as published by*
* the Free Software Foundation; either version 2 of the License, or *
* (at your option) any later version. *
* *
* This program is distributed in the hope that it will be useful, *
* but WITHOUT ANY WARRANTY; without even the implied warranty of *
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
* GNU General Public License for more details. *
* *
* You should have received a copy of the GNU General Public License *
* along with this program; if not, write to the *
* Free Software Foundation, Inc., *
* 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA *
* or see https://www.gnu.org/licenses/licenses.html#GPL . *
*************************************************************************/
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>
#include <gsl/gsl_qrng.h>
#ifdef __GNUC__
#include <stdio.h>
#include <strings.h>
#endif
#define N_SIGMA_GAUSS 1.0
void wrap_gsl_rng_uniform(long iseed,
int i_free, /* 0 = initialise, 1 = free the memory */
double *ran_double);
void wrap_gsl_qrng_3D(int i_free, /* 0 = initialise, 1 = free the memory */
double ran_double_3D[3]);
void gauss_gsl(long Seed,
double *RandNum,
int free_gauss_gsl);
/* random double precision float generator
Especially written for convenient calls
from fortran90, requiring only a single fortran-side
wrapper for initialising, calling, and freeing.
*/
void wrap_gsl_rng_uniform(long iseed,
int i_free, /* 0 = initialise, 1 = free the memory */
double *ran_double)
{
static const gsl_rng_type * T_gsl;
static gsl_rng * r_gsl;
static int have_initialised = 0;
/* MT19937 - Makoto Matsumoto and Takuji Nishimura;
Mersenne prime period of 2^19937 - 1;
second revision 2002 - info gsl for details */
/* T_gsl = gsl_rng_mt19937; */
/* RANLUX algorithm of L\"uscher;
period \sim 10^171;
second generation algorithm;
double-precision, highest luxury level;
slow: about 7-10 times slower than MT19937;
claims of well-tested randomness - info gsl for more
*/
T_gsl = gsl_rng_ranlxd2;
if(!have_initialised){
gsl_rng_default_seed = (unsigned long) iseed;
r_gsl = gsl_rng_alloc(T_gsl);
have_initialised = 1;
};
if(!i_free){
*ran_double = gsl_rng_uniform(r_gsl);
}else if(have_initialised){
gsl_rng_free(r_gsl);
have_initialised = 0;
};
}
void wrap_gsl_qrng_3D(int i_free, /* 0 = initialise, 1 = free the memory */
double ran_double_3D[3])
{
static gsl_qrng * r_gsl;
static int have_initialised = 0;
if(!have_initialised){
r_gsl = gsl_qrng_alloc(gsl_qrng_niederreiter_2, 3);
have_initialised = 1;
};
if(!i_free){
gsl_qrng_get(r_gsl, ran_double_3D);
}else if(have_initialised){
gsl_qrng_free(r_gsl);
have_initialised = 0;
};
}
/* Independent mpi threads should create independent instances of this
sequence; providing different Seed's will give different sequences. */
void gauss_gsl(long Seed,
double *RandNum,
int free_gauss_gsl){
/* The static attribute should work for mpi (independent memory per
thread). If porting to openmp (shared memory) is needed, test this
before trusting it.
*/
static int first_time = 1;
static const gsl_rng_type * T_gsl;
static gsl_rng * r_gsl = NULL;
static const double n_sigma_gauss = N_SIGMA_GAUSS;
if(!first_time && free_gauss_gsl){
gsl_rng_free(r_gsl);
first_time = 1; /* in case this function gets called again after
freeing*/
}else{
if(first_time){
first_time = 0; /* no longer first time */
/* initialise the random number generator */
T_gsl = gsl_rng_ranlxd2;
/* A faster option could be:
T_gsl = gsl_rng_mt19937;
See the GSL documentation for more generators.
*/
gsl_rng_default_seed = (unsigned long) Seed;
r_gsl = gsl_rng_alloc(T_gsl); /* uses gsl_rng_default_seed */
};
/* generate double-precision gaussian random value
https://en.wikipedia.org/wiki/Ziggurat_algorithm
*/
*RandNum = gsl_ran_gaussian_ziggurat(r_gsl,
n_sigma_gauss);
};
}
|