If I understand what you're doing, your Fortran seed knows nothing
about the R state. Is it possible to switch to using the R random
number generators? Or can you at least use the probability integral
transform to generate uniform [0, 1] via R and then convert those to
observations in Fortran (that's what I did)?
Avi
On Fri, Jan 17, 2020 at 7:10 PM وليد خلف معوض المطيرى
<wkmtie...@qu.edu.sa> wrote:
>
> Hi all,
>
>
>
> I think what I’ve done is something different. Inside the Fortran subroutine,
> I have a subroutine for setting the seed value of the RNG of GNU Fortran
> which I think it is not related to the R RNG like the one below:
>
>
>
> subroutine initrandomseedsinr(temp)
>
> implicit none
>
> integer :: n
>
> integer, intent(in):: temp
>
> integer, dimension(:), allocatable :: seed
>
>
>
> call random_seed(size = n)
>
> allocate(seed(n))
>
> seed = temp
>
> call random_seed(PUT = seed)
>
> deallocate(seed)
>
>
>
> end subroutine initrandomseedsinr
>
> , where temp is an argument of the Fortran subroutine as well as in the
> wrapper R function. This will be related to the RNG method used in the GNU
> Fortran that build on GCC not to the R. I am not sure if I am right on this,
> but tried with using RNGkind(sample.kind = "Rounding") and it doesn’t help.
> The difference in the results were not major. The output at the end of
> running the functions kept having very similar results, but still have the
> issue of reproducing exact results which I need it for relating work that is
> based on the package.
>
>
>
> Many thanks,
>
>
>
> Waleed
>
>
>
>
>
>
>
>
>
> Sent from Mail for Windows 10
>
>
>
> From: Avraham Adler
> Sent: Friday, January 17, 2020 06:30 م
> To: Ivan Krylov
> Cc: وليد خلف معوض المطيرى; R Package Development
> Subject: Re: [R-pkg-devel] For reproducibility issue
>
>
>
> 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.r...@gmail.com> wrote:
> >
> > On Fri, 17 Jan 2020 13:55:39 +0000
> > وليد خلف معوض المطيرى <wkmtie...@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@r-project.org mailing list
> > https://stat.ethz.ch/mailman/listinfo/r-package-devel
>
>
______________________________________________
R-package-devel@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-package-devel
