gpvmc/RanGen.cpp at master · bdallapi/gpvmc · GitHub

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
#include "RanGen.h"
#include <iostream>
#include <cstdlib>

#ifdef USE_RNG_MKL
    #ifdef __cplusplus
        extern "C" {
    #endif//__cplusplus
    #include <mkl_vsl.h>
    #ifdef __cplusplus
        }
    #endif//__cplusplus
#elif defined USE_RNG_GSL
    #include <gsl/gsl_rng.h>
#elif defined USE_RNG_ACML
        namespace acml {
            #include <acml.h>
    }
#elif defined USE_RNG_STD
    #include <cstdlib>
#endif//USE_RNG_MKL

struct state_t{
#ifdef USE_RNG_MKL
        VSLStreamStatePtr state;
#elif defined USE_RNG_ACML
        int *state;
#elif defined USE_RNG_GSL
        gsl_rng* state;
#endif
};


RanGen::RanGen()
{
    m_rng=new state_t;
#ifdef USE_RNG_MKL
    vslNewStream(&m_rng->state,VSL_BRNG_MT19937,1);
#elif defined USE_RNG_ACML
    m_rng->state=new int[650];
#elif defined USE_RNG_GSL
    m_rng->state=gsl_rng_alloc(gsl_rng_ranlxd2);
    gsl_rng_set(m_rng->state,1);
#endif
}

RanGen::~RanGen()
{
#ifdef USE_RNG_MKL
    vslDeleteStream(&m_rng->state);
#elif defined USE_RNG_ACML
    delete [] m_rng->state;
#elif defined USE_RNG_MKL
    gsl_rng_free(m_rng->state);
#endif
    delete m_rng;
}

void RanGen::srand(size_t seed)
{
#ifdef USE_RNG_MKL
    vslDeleteStream(&m_ran.m_rng->state);
    vslNewStream(&m_ran.m_rng->state,VSL_BRNG_MT19937,seed);
#elif defined USE_RNG_ACML
    int sseed(seed),info(0),lseed(1), lstate(650);
    acml::drandinitialize(3,0,&sseed,&lseed,m_ran.m_rng->state,&lstate,&info);
#elif defined USE_RNG_GSL
    gsl_rng_set(m_ran.m_rng->state,seed);
#elif defined USE_RNG_STD
    std::srand(seed);
#else
    std::cerr<<"No random generator set at build time. Abort."<<std::endl;
    std::abort();
#endif
}

double RanGen::uniform()
{
    double out(0);
#ifdef USE_RNG_MKL
    vdRngUniform(VSL_METHOD_DUNIFORM_STD_ACCURATE,m_ran.m_rng->state,1,&out,0,1);
#elif defined USE_RNG_ACML
    int info;
    acml::dranduniform(1,0.0,1.0,m_ran.m_rng->state,&out,&info);
#elif defined USE_RNG_GSL
    out=gsl_rng_uniform_pos(m_ran.m_rng->state);
#elif defined USE_RNG_STD
    out=((double)std::rand())/((double)RAND_MAX);
#else
    std::cerr<<"No random generator set at build time. Abort."<<std::endl;
    std::abort();
#endif
#ifdef DEBUG
    if(out==0){
        std::cerr<<"Random generator returned 0. Looks suspicious."<<std::endl;
        std::abort();
    }
#endif

    return out;
}

RanGen RanGen::m_ran;