Linear congruential generator issues

[ivan-pi]

I've noticed that BFPack uses some linear congruential generators, e.g.

https://github.com/jomulder/BFpack/blob/cfec7783577f45618457d620be88195c6d6321d5/src/bct_prior.f90#L305-L310

In gfortran-13 the behavior on integer overflow has changed (see also the bottom of the page Porting to GCC 13 or the discussion on comp.lang.fortran). To solve this issue, the procedures which rely on wraparound on overflow should be isolated and compiled with the flag -fwrapv (gfortran-13).

If you don't want to rely on compilation flags, the options remaining are:

• implement the functions in C which has unsigned integers with proper wraparound semantics (recommended here: https://community.intel.com/t5/Intel-Fortran-Compiler/Integer-overflow/m-p/747387)
• use the procedures in the R library (libR.so): https://cran.r-project.org/doc/manuals/r-release/R-exts.html#Random-numbers

I'm happy to provide a pull request with the interfaces needed.

[cjvanlissa]

@ivan-pi just wanted to say thank you for your help with these issues!

[jomulder]

Hi Ivan, thanks for this. About a year ago we actually used the random_seed subroutine but then we read that this was not recommendable anymore for R package, so we switched to linear congruential generators. But apparently this was also not the best idea (maybe due to new updates). Given the options you mention I would rather stay as much as possible in Fortran (rather than using C wrappers for instance) for my understanding. I read that this module is recommended

module rngfuncs
        use iso_c_binding
        interface
          double precision
     *      function unifRand() bind(C, name = "unif_rand")
          end function unifRand

          subroutine getRNGseed() bind(C, name = "GetRNGstate")
          end subroutine getRNGseed

          subroutine putRNGseed() bind(C, name = "PutRNGstate")
          end subroutine putRNGseed
        end interface
end module

So (for my understanding) in my subroutines I can replace

randraw = runiform(iseed)

with

call getRNGseed()
randraw = unifRand()
call putRNGSeed()

So I would need to put call getRNGseed() and call putRNGSeed() every time before and after the function unifRand is used? And how is the seed controlled via the route? Note that in my previous code I used the integer iseed which is given from R as input to the Fortran subroutine. Thanks again for any advice on this!

[ivan-pi]

The R documentation is horribly inconsistent when it comes to describing Fortran features, not just in term of how to specify precision, but also in terms of styles. (Understandably, as I imagine the core R team doesn't keep close track of Fortran development.)

The module, here using free-form syntax, should be:

! rngfuncs.f90
module rngfuncs
implicit none
public
! See the R header `R_ext/Random.h` for details
interface
   ! double unif_rand(void);
   function unif_rand() bind(c,name="unif_rand")
      use, intrinsic :: iso_c_binding, only: c_double
      real(c_double) :: unif_rand
   end function

   ! void GetRNGstate(void);
   subroutine getrngstate() bind(c,name="GetRNGstate")
   end subroutine

   ! void PutRNGstate(void);
   subroutine putrngstate() bind(c,name="PutRNGstate")
   end subroutine
end interface
end module

(You are free to rename the Fortran interfaces to your liking, but IMO this is silly. Fortran syntax is case-insensitive, so we need to use the explicit name=... property.)

I think you don't need to call them every time you call the function, but only once at the beginning before you plan to use it repeatedly. These procedures will set an implicit seed for you. I'm not an R expert, but I imagine you can control the seed from R via set.seed(iseed).

For example to fill an array with uniform variates you would do:

real(c_double) :: a(100)

call getrngstate

do i = 1, size(a)
   a(i) = unif_rand()
end do

call putrngstate

If you are writing the Fortran procedures with the aim to use them from R only, then I think it makes sense to just embrace what the R environment gives you.

[jomulder]

Hi Ivan, thanks a lot for this (again). I tried this and it works fine when setting the seed in R using set.seed(). Btw the draws are also controlled via the seed when I omit the lines call getrngstate and call putrngstate. Not sure if that's supposed to happen though...

[jomulder]

Hi Ivan, these suggestions have resolved many problems for the linux compilers on rhub, so big thanks for this!!!! On a windows check though, I do get a (new) error related to these new function calls of the form:

C:\rtools44\x86_64-w64-mingw32.static.posix\bin/ld.exe: bct_mixedordinal.o:bct_mixedordinal.f90:(.text+0x482): undefined reference todpotrf_'`

Is this familiar to you?

For instance, see link

[ivan-pi]

Those are linker errors. The linker cannot find the library containing the LAPACK symbols (https://cran.r-project.org/doc/FAQ/R-exts.html#index-LAPACK_005fLIBS). Currently the command at which the build fails is

gcc -shared -s -static-libgcc -o BFpack.dll tmp.def BFpack_init.o bct_continuous_final.o bct_mixedordinal.o bct_mixedordinal_old.o bct_prior.o test_kind.o -LC:/rtools44/x86_64-w64-mingw32.static.posix/lib/x64 -LC:/rtools44/x86_64-w64-mingw32.static.posix/lib -lgfortran -lm -lquadmath -LC:/R/bin/x64 -lR

If I understand the R docs correctly, you should add the following line into Makevars:

PKG_LIBS = @LAPACK_LIBS@

This should add -lRlapack to end of the previous command, allowing the linker to find the LAPACK procedures.

[ivan-pi]

Hi Ivan, thanks a lot for this (again). I tried this and it works fine when setting the seed in R using set.seed(). Btw the draws are also controlled via the seed when I omit the lines call getrngstate and call putrngstate. Not sure if that's supposed to happen though...

Interesting. I'm not familiar enough with R to know exactly why. Judging by the implementation (https://github.com/wch/r-source/blob/b7c33a61c4b3ca973b95fe4aedf345fc8b858273/src/main/RNG.c#L409), if no seed has been set, it will initialize the PRNG.

[jomulder]

Thanks again Ivan! Indeed that was the cause. I did need to alter the line slightly to

PKG_LIBS = $(LAPACK_LIBS) $(BLAS_LIBS)

to avoid the error.