quickrnd.F90

Go to the documentation of this file.

!! Copyright (C) 2002-2024 M. Marques, A. Castro, A. Rubio, G. Bertsch, X. Andrade, A. Buccheri, H. Menke
!!
!! This program is free software; you can redistribute it and/or modify
!! it under the terms of the GNU General Public License as published by
!! the Free Software Foundation; either version 2, or (at your option)
!! any later version.
!!
!! This program is distributed in the hope that it will be useful,
!! but WITHOUT ANY WARRANTY; without even the implied warranty of
!! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
!! GNU General Public License for more details.
!!
!! You should have received a copy of the GNU General Public License
!! along with this program; if not, write to the Free Software
!! Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
!! 02110-1301, USA.
!!
 
#include "global.h"
 
module quickrnd_oct_m
  use, intrinsic :: iso_fortran_env
  use debug_oct_m
  use global_oct_m
  use messages_oct_m
  use profiling_oct_m
 
  implicit none
 
  private
  public ::                     &
    quickrnd,                   &
    shiftseed,                  &
    fisher_yates_shuffle
 
  ! xorshift64 needs a good seed, in particular it cannot be seeded with 0.
  ! These seeds are generated with the splitmix64 generator and have a
  ! reasonable amount of entropy.
  integer(int64), parameter, public ::                &
    splitmix64_123 = int(z'B4DC9BD462DE412B', int64), & ! splitmix64(123)
    splitmix64_321 = int(z'E95F1C4EFCF85DEE', int64)    ! splitmix64(321)
 
  interface quickrnd
    module procedure dquickrnd_single, dquickrnd_array, zquickrnd_array
  end interface quickrnd
 
contains
 
  function xorshift64(x) result(res)
    integer(int64), intent(inout) :: x
    integer(int64) :: res
 
    assert(x /= 0_int64)
 
    x = ieor(x, shiftl(x, 13))
    x = ieor(x, shiftr(x, 7))
    x = ieor(x, shiftl(x, 17))
 
    res = x
  end function xorshift64
 
  function to_double(x) result(d)
    integer(int64), intent(in) :: x
    real(real64) :: d
    real(real64), parameter :: two_pow_minus_53 = real(z'3CA0000000000000', real64)
 
    d = real(shiftr(x, 11), real64) * two_pow_minus_53
  end function to_double
 
  ! ---------------------------------------------------------
 
  subroutine dquickrnd_single(iseed, rnd)
    integer(int64), intent(inout) :: iseed
    real(real64), intent(out) :: rnd
 
    ! no PUSH_SUB, called too often
 
    rnd = to_double(xorshift64(iseed))
 
  end subroutine dquickrnd_single
 
  ! ---------------------------------------------------------
 
  subroutine dquickrnd_array(iseed, nn, rnd)
    integer(int64), intent(inout) :: iseed
    integer, intent(in) :: nn
    real(real64), intent(inout) :: rnd(:)
 
    integer :: ii
 
    push_sub(quickrnd_array)
 
    do ii = 1, nn
      rnd(ii) = to_double(xorshift64(iseed))
    end do
 
    pop_sub(quickrnd_array)
 
  end subroutine dquickrnd_array
 
  ! ---------------------------------------------------------
 
  subroutine zquickrnd_array(iseed, nn, rnd)
    integer(int64), intent(inout) :: iseed
    integer, intent(in) :: nn
    complex(real64), intent(inout) :: rnd(:)
 
    real(real64), parameter :: inv_sqrt2 = 1.0_real64 / sqrt(m_two)
    real(real64) :: re, im
    integer :: ii
 
    push_sub(quickrnd_array)
 
    do ii = 1, nn
      re = to_double(xorshift64(iseed))
      im = to_double(xorshift64(iseed))
      rnd(ii) = cmplx(re, im, real64) * inv_sqrt2
    end do
 
    pop_sub(quickrnd_array)
 
  end subroutine zquickrnd_array
 
  ! ---------------------------------------------------------
 
  subroutine shiftseed(iseed, n)
    integer(int64), intent(inout) :: iseed
    integer(int64), intent(in) :: n
 
    integer(int64) :: discard
    integer(int64) :: ii
 
    push_sub(shiftseed)
 
    do ii = 1, n
      discard = xorshift64(iseed)
    end do
 
    pop_sub(shiftseed)
 
  end subroutine shiftseed
 
 
  function random_integer_in_range(x_min, x_max, seed) result(rand)
    integer(int64), intent(in   ) :: x_min, x_max
    integer(int64), intent(inout) :: seed
    integer(int64)                :: rand
 
    integer(int64) :: range
 
    assert(seed /= 0_int64)
    seed = xorshift64(seed)
    if (seed < 0_int64) seed = -seed  ! Enforce a positive result
    range = x_max - x_min + 1_int64
    rand = x_min + mod(seed, range)
 
  end function random_integer_in_range
 
 
  subroutine fisher_yates_shuffle(m, n, seed, values)
    integer(int32), intent(in   ) :: m
    integer(int64), intent(in   ) :: n
    integer(int64), intent(inout) :: seed
    integer(int64), intent(  out) :: values(:)
 
    integer(int64), allocatable :: indices(:)
    integer(int64) :: i, j, temp
 
    push_sub(fisher_yates_shuffle)
 
    safe_allocate(indices(1:n))
    do j = 1, n
      indices(j) = j
    end do
 
    do i = 1, m
      j = random_integer_in_range(i, n, seed)
      ! Perform the swap
      temp = indices(i)
      indices(i) = indices(j)
      indices(j) = temp
      values(i) = indices(i)
    end do
 
    values(m) = indices(m)
    safe_deallocate_a(indices)
 
    pop_sub(fisher_yates_shuffle)
 
  end subroutine fisher_yates_shuffle
 
end module quickrnd_oct_m
 
!! Local Variables:
!! mode: f90
!! coding: utf-8
!! End: