[Rd] RNG stuck via Fortran call

[Duncan Murdoch]

On 11/30/2005 10:06 AM, Gilles GUILLOT wrote:
> Having not much success with my previous question I try to reformulate it:
> 
> I'm simulating a Markow chain in Fortran interfaced with R.
> Each loop of my Fortran calls various functions of the R RNG through 
> the wrapper given below.
> 
> In a run of 100 iterations of the Markov kernel, 
>  after 20 iterations, the RNG seems like frozen.
> 
> For example, the first call to the RNG in my loop is:
> 
>  rpostlamb = ggrgam(dble(m),1.d0)
>  write(*,*) 'rpostlamb=',rpostlamb
> 
> 
> after the 20th iteration,  it will return the same value 
> rpostlamb=  1.24634557
> 
> for all the following interations.
> 
> Is there any suggestion of explanation for this strange fact ?

You need to put together code that someone else could compile and run, 
or show us your real code if it's short enough.  It sounds as though 
you're not handling the RNG state properly.

Please don't send this code to me, post it to the list.  (Actually, in 
putting together a short reproducible example, you're very likely to 
discover the bug yourself, and won't need to post anything.)

Duncan Murdoch
> 
> 
> Thanks in advance
> 
> Gilles 
> 
> R Version 2.2.0 compiled under Mandrake 10.1
> 
> 
> #include <R.h>
> #include <Rmath.h>
> 
> /******************/
> /* random numbers */
> /******************/
> 
> void F77_SUB(rndstart)(void) 
> { GetRNGstate(); }
> 
> void F77_SUB(rndend)(void) 
> { PutRNGstate(); }
> 
> double F77_SUB(ggrnorm)(double *mu, double *sigma) 
> { return rnorm(*mu, *sigma); }
> 
> double F77_SUB(ggrexp)(double *scale)
> {return rexp(*scale);}
> 
> double F77_SUB(ggrgam)(double *a, double *scale)
> {return rgamma(*a, *scale);}
> 
> double F77_SUB(ggrunif)(double *a, double *b)
> {return runif(*a, *b);}
> 
> double F77_SUB(ggrbinom)(double *n, double *p)
> {return rbinom(*n, *p);}
> 
> double F77_SUB(ggrpois)(double *lambda)
> {return rpois(*lambda);}
> 
> ______________________________________________
> R-devel at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-devel