Re: [Blitz-support] parallel MT19937

[Patrik J. <pa...@uc...>]

Hi Patrick with a c, :-)

Patrick Guio wrote:
>
> Dear Patrik,
>
> I have noticed that the Mersenne Twister random generator was patched
> for parallel use.
> I have been using the Blitz implementation of the Mersenne Twister
> random generator with MPI for sometime before your patch was added.
> My approach is that each process has its own random generator seeded
> with a different value. The seed value is calculated as the sum of a
> predefined seed value "seed" and the rank "rankProc" of the MPI process.
> I have attached a code snippet to illustrate this.
> The result is different random series drawn for each process.
That's how I used to do it, too. But since I'm basing my research 
results on these random numbers (using them for Monte Carlo) I thought 
it would be prudent to heed the warnings of the MT19937 authors about 
parallel use.

I don't know if you've seen the paper on parallel MT 
(http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/DC/dc.html) but the 
idea is that if you use the different seed approach, you are just 
starting at different points in the fixed (but long) MT19937 sequence. 
While the sequence is randomly distributed to a large degree, there are 
no guarantees that there aren't correlations between r(i) and r(i+N). 
(I'm summarizing this mostly for the benefit of those who didn't read my 
original email about the patch back this summer.)
>
> Now I would like to try the patched version but I cannot figure out how
> to use it.
>
> I just replaced the constructor and seeding statements
>
> <   ranlib::Uniform<double> r;
> <   r.seed(seed + rankProc);
> ---
>  >   ranlib::Uniform<double> r(rankProc);
>  >   r.seed(seed);
>
> First the code does not compile and I get an error that there is no
> suitable constructor for
>
> /usr/local/include/random/uniform.h(79): error: no suitable constructor
> exists to convert from "unsigned int" to
> "ranlib::IRNGWrapper<ranlib::defaultIRNG, ranlib::defaultState>"
>        IRNGWrapper<IRNG,stateTag>::IRNGWrapper(i) {};
>                                                ^
> I guess the problem is that the defaultState is sharedState and not
> independentState but I cannot figure out how to proceed.
>
> Could you provide also a small example of how to use the parallel 
> version?
The idea is that you instantiate a different random number generator. I 
do something like:

ranlib::Uniform<double, ranlib::MersenneTwister, 
ranlib::independentState> rng(i);

If you want the IRNG directly, then the above just instantiates  
IRNGWrapper<ranlib::MersenneTwister, ranlib::independentState>(i), so 
that should work. At least I can compile it ok. After you've 
instantiated the generator, you use it like the normal ones.

It doesn't make sense to generate the parallel generators with 
sharedState, because they have different recursions and will clobber 
each others state vectors and who knows what will happen... I don't 
remember if I thought about that and that's why it doesn't compile, but 
change it to independentState and it should work.

Also note that, unlike most of blitz functionality, you need to link the 
blitz library for this to work, because that's where the parameter sets are.

hope that helps,

/Patrik