[R-pkg-devel] For reproducibility issue

[Avraham Adler]

Hi.

If it helps, I call the R RNG from Fortran in my Delaporte package
[1], also using iso_c_bindings. Specifically, I have the following C
code [2]:

void F77_SUB(unifrnd) (int *n, double *x){
GetRNGstate();
for (int i = 0; i < *n; ++i){
*(x + i) = unif_rand();
}
PutRNGstate();
}
and call it in Fortran like so [3]:

subroutine rdelap_f(n, a, na, b, nb, l, nl, vars) bind(C, name="rdelap_f_")

external unifrnd

integer(kind = c_int), intent(in), value :: n, na, nb, nl
real(kind = c_double), intent(out), dimension(n) :: vars
real(kind = c_double), intent(in) :: a(na), b(nb), l(nl)
real(kind = c_double), dimension(n) :: p
integer(kind = c_int) :: lg, lt

call unifrnd(n, p)
lt = 1_c_int
lg = 0_c_int
call qdelap_f(p, n, a, na, b, nb, l, nl, lt, lg, vars)

end subroutine rdelap_f

The package passes CRAN tests (at least as of now) on anything between
GCC 4 and GCC9 [4].

Hope that helps,

Avi

[1] https://bitbucket.org/aadler/delaporte/src/master/
[2] https://bitbucket.org/aadler/delaporte/src/master/src/utils_and_wrappers.c
[3] https://bitbucket.org/aadler/delaporte/src/master/src/delaporte.f95
[4] https://cran.r-project.org/web/checks/check_results_Delaporte.html


On Fri, Jan 17, 2020 at 2:39 PM Ivan Krylov <krylov.r00t using gmail.com> wrote:
>
> On Fri, 17 Jan 2020 13:55:39 +0000
> وليد خلف معوض المطيرى <wkmtierie using qu.edu.sa> wrote:
>
> > So, does anyone have an idea of how to solve this issue.
>
> "Writing R Extensions", 1.6. Writing portable packages:
>
> >> Compiled code should not call the system random number generators
> >> such as rand, drand48 and random, but rather use the interfaces to
> >> R’s RNGs described in Random numbers. In particular, if more than
> >> one package initializes the system RNG (e.g. via srand), they will
> >> interfere with each other.
>
> >> Nor should the C++11 random number library be used, nor any other
> >> third-party random number generators such as those in GSL.
>
> It somewhat less convenient to call the R random number generator from
> Fortran than it would be from C or C++, but still possible. There is a
> F77-style example of such use [1], but since you are already using
> iso_c_binding, you should be able to declare the C API [2] right in the
> Fortran source:
>
> subroutine GetRNGState() bind(c)
> end subroutine
>
> subroutine PutRNGstate() bind(c)
> end subroutine
>
> As a bonus, you get to use the R distribution functions [3], without
> the need to implement them yourself from uniformly distributed samples:
>
> function rnorm(mu, sigma) bind(c) result(ret)
>  use intrinsic, iso_c_binding, only: c_double
>  real(c_double), value :: mu, sigma
>  real(c_double) ret
> end function
>
> function rgamma(shape, scale) bind(c) result(ret)
>  use intrinsic, iso_c_binding, only: c_double
>  real(c_double), value :: shape, scale
>  real(c_double) ret
> end function
>
> (The prototypes above are unchecked; I haven't written any Fortran 2003
> in more than a year.)
>
> --
> Best regards,
> Ivan
>
> [1]:
> https://cran.r-project.org/doc/manuals/R-exts.html#index-Random-numbers-in-Fortran
>
> [2]: https://cran.r-project.org/doc/manuals/R-exts.html#Random-numbers
>
> [3]:
> https://cran.r-project.org/doc/manuals/R-exts.html#Distribution-functions
>
> ______________________________________________
> R-package-devel using r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-package-devel