/*
 * 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. (address@hidden)\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;

}