Ok, I think I got the idea.
I don't use c++11, however I think that the ansi c code attached works in the
same way.
I allocated memory for a number = omp_get_max_threads() of generators.
Thank you very much,
Al.
On Sunday, 4 May 2014, 14:14, Klaus Huthmacher <[email protected]>
wrote:
Dear Altro,
I would simply use a wrapper around a C++ std::vector who holds as many
C++11 random number generators (engines) as you use threads in your code.
By calling a random number, the engine in the ith index of the vector is
used, where i is the thread id.
Please see the appended minimal example mb.cpp with the wrapper
rng-omp.{h,cpp}.
Best wishes and keep us informed,
-- Klaus.
> Dear All
>
> I wish to use effectively use the gsl random number generator jointly
> with OpenMP on a HPC cluster.
>
> The simulation I have to do are quite simple: just let evolve an
> "ensemble" of system with the same dynamics. Parallelize such a code
> with OpenMP "seems" straightforward.
> As the GSL does not support parallel processing, I suppose that we must
> use a different random number generator for each element of the ensemble;
> Indeed I don't know if this is a good procedure.
>
> The ensemble could be very big and it would be simply impossible to
> allocate memory for all generators associated with the Systems of the
> ensemble.
> Then my question is:
> Is this still a correct way to use the GSL random generator on OpenMP?
>
> An idea of how the allocation could be done is given in the code below.
> Does it make any sense? If not, which is the correct way?
>
> Best wishes,
>
> Al.
>
>
>
> #include <stdlib.h>
> #include <stdio.h>
> #include <gsl/gsl_rng.h>
> #include <gsl/gsl_randist.h>
> #include <omp.h>
>
> int main (){
> int i;
>
> int DIM_ENSEMBLE = 10000000000;
>
> double *d;
> gsl_rng **r ;
>
> d = malloc(DIM_ENSEMBLE*sizeof(double));
> r = malloc(DIM_ENSEMBLE*sizeof(gsl_rng*))
> ; ;
> for (i=0;i<DIM_ENSEMBLE;i++){
> r[i] = gsl_rng_alloc (gsl_rng_mt19937);
> // r[i] = gsl_rng_alloc (gsl_rng_taus);
> gsl_rng_set(r[i],i);
> }
>
>
> size_t n = gsl_rng_size (r[0]); //size in byte (the size of
> tausworthe is 24 ; the size of mersenne twister is 5000)
>
> #pragma omp parallel for default(none) shared(r,d) shared(DIM_ENSEMBLE)
> for (i =0; i< DIM_ENSEMBLE; i++) {
> d[i] = gsl_rng_uniform(r[i]);
> }
>
> #pragma omp parallel for default(none) shared(d) shared(DIM_ENSEMBLE)
> for (i=0;i<DIM_ENSEMBLE;i++){
> printf("%d %f\n",i,d[i]);
> }
>
> printf("The size in Mb of the vector r is %lu
> Mb\n",n*DIM_ENSEMBLE/1000000);
>
> free(d);
>
> for (i=0;i<DIM_ENSEMBLE;i++){
> gsl_rng_free (r[i]);
> }
> free(r);
>
> return 0;
> }
>
>/*
* gcc gsl-rng-parallel.c -lm -lgsl -lgslcblas -fopenmp
*/
#include <stdlib.h>
#include <stdio.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>
#include <omp.h>
void usage(int argc, char *argv[]){
printf("\n This program demonstrates how (not) to use the GSL Random Number Generator in OpenMP environment\n"
"\n\n"
" This is Free Software - You can use and distribute it under \n"
" the terms of the GNU General Public License, version 3 or later\n\n"
" (c) Al. Tro. ([email protected])\n\n");
printf("Usage: %s [DIM. ENSAMBLE]\n" , argv[0]);
}
int main (int argc, char *argv[]){
if (argc != 2 ){
usage(argc,argv);
exit(1);
}
int i;
int DIM_ENSEMBLE = atoi(argv[1]);
printf("n of threads is %d\n", 4);
int n_threads = 4;
n_threads = omp_get_max_threads();
printf("n of threads is %d\n", n_threads);
double *d;
gsl_rng **r ;
d = malloc(DIM_ENSEMBLE*sizeof(double));
r = malloc(n_threads*sizeof(gsl_rng*));
for (i=0;i<n_threads;i++){
r[i] = gsl_rng_alloc (gsl_rng_mt19937);
// r[i] = gsl_rng_alloc (gsl_rng_taus);
gsl_rng_set(r[i],i); // seed each rng
}
size_t n = gsl_rng_size (r[0]); //size in byte (the size of tausworthe is 24 ; the size of mersenne twister is 5000)
// r = gsl_rng_alloc (gsl_rng_gfsr4);
// n = gsl_rng_size (r); //the size is bery big... 131080
// printf ("%lu\n",n);
// double ciao = gsl_rng_uniform(r);
// n = gsl_rng_size(r);
// printf ("%lu\n",n);
#pragma omp parallel for default(none) shared(r,d) shared(DIM_ENSEMBLE)
for (i =0; i< DIM_ENSEMBLE; i++) {
d[i] = gsl_rng_uniform(r[ omp_get_thread_num() ]);
}
#pragma omp parallel for default(none) shared(d) shared(DIM_ENSEMBLE)
for (i=0;i<DIM_ENSEMBLE;i++){
printf("%d %f\n",i,d[i]);
}
printf("The size of the vector r is %lu Mb\n",n*DIM_ENSEMBLE/1000000);
free(d);
for (i=0; i<n_threads; i++){
gsl_rng_free (r[i]);
}
free(r);
return 0;
}
