Optimizing a Fortran subroutine

Question

I've written a minimal implementation for the fast xoroshiro128plus pseudo-random number generator in Fortran to replace the intrinsic random_number. This implementation is quite fast (4X faster than random_number) and the quality is good enough for my purposes, I don't use it in cryptography applications.

My question is how can I optimize this subroutine to get the last drop of performance from my compiler, even 10% improvement is appreciated. This subroutine is to be used in tight loops inside long simulations. I'm interested more in generating a single random number at a time and not big vectors or nD arrays at once.

Here is a test program to give you some context about how my subroutine is used:

program test_xoroshiro128plus
   implicit none
   integer, parameter :: n = 10000
   real*8  :: A(n,n)
   integer :: i, j, t0, t1, count_rate, count_max

   call system_clock(t0, count_rate, count_max)
   do j = 1,n
      do i = 1,n
         call drand128(A(i,j))
      end do
   end do
   ! call drand128(A)  ! works also with 2D 
   call system_clock(t1)

   print *, "Time :", real(t1-t0)/count_rate
   print *, "Mean :", sum(A)/size(A), char(10), A(1:2,1:3)

 contains

   impure elemental subroutine drand128(r)
      real*8, intent(out) :: r
      integer*8 :: s0 = 113, s1 = 19937
      s1 = xor(s0,s1)
      s0 = xor(xor(ior(ishft(s0,55), ishft(s0,-9)),s1), ishft(s1,14))
      s1 = ior(ishft(s1,36), ishft(s1,-28))
      r = ishft(s0+s1, -1) / 9223372036854775808.d0
   end 

end program

Answer 1

Only now I realized you are asking about this particular PRNG. I am using it in Fortran myself https://bitbucket.org/LadaF/elmm/src/eb5b54b9a8eb6af158a38038f72d07865fe23ee3/src/rng_par_zig.f90?at=master&fileviewer=file-view-default

My code in the link is slower than yours, because it calls several subroutines and aims to be more universal. Bet let's try to condense the code I use into a single subroutine.

So let's just compare the performance of your code and the optimized version by @SeverinPappadeux and my optimized code with Gfortran 4.8.5

> gfortran -cpp -O3 -mtune=native xoroshiro.f90 

 Time drand128 sub:   1.80900002    
 Time drand128 fun:   1.80900002    
 Time rng_uni:   1.32900000 

the code is here, remember to let the CPU spin-up, the first iteration of the k loop is just garbage!!!

program test_xoroshiro128plus
   use iso_fortran_env       
   implicit none
   integer, parameter :: n = 30000
   real*8  :: A(n,n)
   real*4  :: B(n,n)
   integer :: i, j, k, t0, t1, count_rate, count_max       

   integer(int64) :: s1 = int(Z'1DADBEEFBAADD0D0', int64), s2 = int(Z'5BADD0D0DEADBEEF', int64)

!let the CPU spin-up                                           
do k = 1, 3                                           
   call system_clock(t0, count_rate, count_max)
   do j = 1,n
      do i = 1,n
         call drand128(A(i,j))
      end do
   end do
   ! call drand128(A)  ! works also with 2D 
   call system_clock(t1)

   print *, "Time drand128 sub:", real(t1-t0)/count_rate

   call system_clock(t0, count_rate, count_max)
   do j = 1,n
      do i = 1,n
         A(i,j) = drand128_fun()
      end do
   end do
   ! call drand128(A)  ! works also with 2D 
   call system_clock(t1)

   print *, "Time drand128 fun:", real(t1-t0)/count_rate


   call system_clock(t0, count_rate, count_max)
   do j = 1,n
      do i = 1,n
         call rng_uni(A(i,j))
      end do
   end do
   call system_clock(t1)

   print *, "Time rng_uni:", real(t1-t0)/count_rate
end do

   print *, "Mean :", sum(A)/size(A), char(10), A(1:2,1:3)

 contains

   impure elemental subroutine drand128(r)
      real*8, intent(out) :: r
      integer*8 :: s0 = 113, s1 = 19937
      s1 = xor(s0,s1)
      s0 = xor(xor(ior(ishft(s0,55), ishft(s0,-9)),s1), ishft(s1,14))
      s1 = ior(ishft(s1,36), ishft(s1,-28))
      r = ishft(s0+s1, -1) / 9223372036854775808.d0
   end 

   impure elemental real*8 function drand128_fun()
     real*8, parameter :: c = 1.0d0/9223372036854775808.d0
     integer*8 :: s0 = 113, s1 = 19937
     s1 = xor(s0,s1)
     s0 = xor(xor(ior(ishft(s0,55), ishft(s0,-9)),s1), ishft(s1,14))
     s1 = ior(ishft(s1,36), ishft(s1,-28))
     drand128_fun = ishft(s0+s1, -1) * c
  end

  impure elemental subroutine rng_uni(fn_val)
    real(real64), intent(inout) ::  fn_val
    integer(int64) :: ival

    ival = s1 + s2

    s2 = ieor(s2, s1)
    s1 = ieor( ieor(rotl(s1, 24), s2), shiftl(s2, 16))
    s2 = rotl(s2, 37)    

    ival  = ior(int(Z'3FF0000000000000',int64), shiftr(ival, 12))
    fn_val = transfer(ival, 1.0_real64) - 1;    
  end subroutine

  function rotl(x, k)
    integer(int64) :: rotl
    integer(int64) :: x
    integer :: k

    rotl = ior( shiftl(x, k), shiftr(x, 64-k))
  end function    

end program

The main difference should come from the faster and better way to convert from integers to reals http://experilous.com/1/blog/post/perfect-fast-random-floating-point-numbers#half-open-range

If you are bored, you could try to inline rotl() manually, but I trust the compiler here.

Answer 2

Ok, here is my attempt. First, I made it to function - in x64 or similar ABI function returning float value do in in register - much faster than parameter transfer. Second, replaced final division by multiplication, though Intel compiler might do it for you.

Timing, Intel i7 6820, WSL, Ubuntu 18.04:

before -   0.850000024
after  -   0.601000011

GNU Fortran 7.3.0, command line

gfortran -std=gnu -O3 -ffast-math -mavx2 /mnt/c/Users/kkk/Documents/CPP/a.for

Code

  program test_xoroshiro128plus
  implicit none
  integer, parameter :: n = 10000
  real*8  :: A(n,n)
  integer :: i, j, t0, t1, count_rate, count_max

  call system_clock(t0, count_rate, count_max)
  do j = 1,n
     do i = 1,n
        A(i,j) = drand128()
     end do
  end do
  A = drand128()  ! works also with 2D
  call system_clock(t1)

  print *, "Time :", real(t1-t0)/count_rate
  print *, "Mean :", sum(A)/size(A), char(10), A(1:2,1:3)

  contains

  impure elemental real*8 function drand128()
     real*8, parameter :: c = 1.0d0/9223372036854775808.d0
     integer*8 :: s0 = 113, s1 = 19937
     s1 = xor(s0,s1)
     s0 = xor(xor(ior(ishft(s0,55), ishft(s0,-9)),s1), ishft(s1,14))
     s1 = ior(ishft(s1,36), ishft(s1,-28))
     drand128 = ishft(s0+s1, -1) * c
  end

  end program