I'm writing some code in C++ for scientific computing purposes. The linear algebra library I chose to use is Armadillo. My problem is that Armadillo uses either the simple rand() function or a Mersenne Twister rng, while I need better quality numbers, for example ranlux. GSL provides an implementation of ranulx, and Armadillo allows for an external rng by defining the preprocessor macro ARMA_RNG_ALT and writing an appropriate header file that we will call arma_rng_alt.hpp.
Now, I'm no C++ wizard but I managed to do it by more or less copying and pasting from the existing headers arma_rng_cxx98.hpp and arma_rng_cxx11.hpp. My solution might not be the prettiest but I think it works alright. This is an extract of my header:
#if defined(ARMA_RNG_ALT)
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>
#define RNG_TYPE gsl_rng_ranlxd1
class arma_rng_alt
{
public:
typedef unsigned long int seed_type;
inline void set_seed(const seed_type val);
arma_inline int randi_val();
arma_inline double randu_val();
arma_inline double randn_val();
private:
gsl_rng* engine;
};
inline
void
arma_rng_alt::set_seed(const arma_rng_alt::seed_type val)
{
engine = gsl_rng_alloc(RNG_TYPE);
gsl_rng_set(engine, val);
}
arma_inline
int
arma_rng_alt::randi_val()
{
return gsl_rng_get(engine);
}
arma_inline
double
arma_rng_alt::randu_val()
{
return gsl_rng_uniform(engine);
}
arma_inline
double
arma_rng_alt::randn_val()
{
return gsl_ran_gaussian_ziggurat(engine, 1.);
}
#endif
Now, this would be ok if this external rng was treated by armadillo as a class without a private variable (like in the arma_rng_cxx98.hpp header, that uses rand() and doesn't need to allocate an engine). Instead I need to allocate the engine and call it every time, therefore I need to instantiate the class otherwise I get a "cannot call member function without object" error. So I actually needed to modify a header from the Armadillo library itself, namely arma_rng.hpp. To make things clearer, this is how a typical chunk of the original arma_rng.hpp to generate random integers looks like:
template<typename eT>
struct arma_rng::randi
{
arma_inline
operator eT ()
{
#if defined(ARMA_RNG_ALT)
{
return eT( arma_rng_alt::randi_val() );
}
#elif defined(ARMA_USE_EXTERN_CXX11_RNG)
{
return eT( arma_rng_cxx11_instance.randi_val() );
}
#else
{
return eT( arma_rng_cxx98::randi_val() );
}
#endif
}
As you can see, if cxx98 or alt is used, the function randi_val is simply called, while if cxx11 is used, randi_val is called on an actual instance of the calss arma_rng_cxx11 (a thread local variable declared earlier in the header). What I need is something similar to cxx11, so this is how I modified the library header arma_rng.hpp for this specific example:
template<typename eT>
struct arma_rng::randi
{
arma_inline
operator eT ()
{
#if defined(ARMA_RNG_ALT)
{
return eT( arma_rng_alt_instance.randi_val() );
}
#elif defined(ARMA_USE_EXTERN_CXX11_RNG)
{
return eT( arma_rng_cxx11_instance.randi_val() );
}
#else
{
return eT( arma_rng_cxx98::randi_val() );
}
#endif
}
after declaring my own thread local instance of the class arma_rng_alt I've created.
Works like a charm, the catch is that now I need to run my code on a HPC cluster and, while I can create my own arma_rng_alt.hpp, I can't modify the library header arma_rng.hpp itself. So I need a way to make the original header work (so without instantiating the class). Do you have any idea on how to do that? My C++ expertise is frankly stretched to the limit here.
Thanks
Aucun commentaire:
Enregistrer un commentaire