Оптимизация подпрограммы на Фортране

Question

Я написал минимальную реализацию для быстрого xoroshiro128plus генератора псевдослучайных чисел в Фортране для замены встроенного random_number.Эта реализация довольно быстрая (в 4 раза быстрее, чем random_number) и качество достаточно хорошее для моих целей, я не использую его в криптографических приложениях.

Мой вопрос заключается в том, как я могу оптимизировать эту подпрограмму, чтобы получить последнее падение производительности от моего компилятора, даже улучшение на 10% приветствуется.Эта подпрограмма должна использоваться в тесных циклах внутри длинных симуляций.Меня больше интересует генерация одного случайного числа за раз, а не больших векторов или nD-массивов одновременно.

Вот тестовая программа, которая даст вам некоторый контекст о том, как используется моя подпрограмма:

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

Только теперь я понял, что вы спрашиваете об этом конкретном PRNG.Я сам использую его в Фортране https://bitbucket.org/LadaF/elmm/src/eb5b54b9a8eb6af158a38038f72d07865fe23ee3/src/rng_par_zig.f90?at=master&fileviewer=file-view-default

Мой код в ссылке медленнее вашего, потому что он вызывает несколько подпрограмм и стремится быть более универсальным.Давайте попробуем сконцентрировать код, который я использую, в одну подпрограмму.

Итак, давайте просто сравним производительность вашего кода и оптимизированной версии @SeverinPappadeux и моего оптимизированного кода с 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 

код здесь, не забудьте позволить процессору раскрутиться, первая итерация цикла k просто мусор !!!

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

Основное отличие должно быть быстрееи лучший способ конвертировать целые числа в действительные http://experilous.com/1/blog/post/perfect-fast-random-floating-point-numbers#half-open-range

Если вам скучно, вы можете попробовать встроить rotl() вручную, но я доверяю компилятору здесь.

Answer 2

Хорошо, вот моя попытка.Во-первых, я сделал это для функции - в x64 или аналогичной функции ABI, возвращающей значение типа float, в регистре - намного быстрее, чем передача параметра.Во-вторых, заменили окончательное деление на умножение, хотя компилятор Intel может сделать это за вас.

Время, Intel i7 6820, WSL, Ubuntu 18.04:

before -   0.850000024
after  -   0.601000011

GNU Fortran 7.3.0, командалиния

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

код

  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