Генератор случайных чисел в PGI Fortran не такой уж и случайный

Следующий код просто генерирует простую тройку случайных чисел:

program testrand

integer, parameter :: nz = 160, nf = 160, nlt = 90

real :: tmpidx(3)
integer :: idxarr(3), idx1, idx2, idx3, seed_size, ticks
integer, allocatable :: seed(:)

call random_seed(size=seed_size)
allocate(seed(seed_size))

call system_clock(count=ticks)
seed = ticks+37*(/(i-1, i=1,seed_size)/)
call random_seed(put=seed)

deallocate(seed)

call random_number(tmpidx)
idxarr = tmpidx * (/nz, nf, nlt/)
idx1 = max(1,idxarr(1))
idx2 = max(1,idxarr(2))
idx3 = max(1,idxarr(3))
print *,idx1, idx2, idx3

end program

Я компилирую это с помощью gfortran, запускаю несколько раз и получаю:

> gfortran testrand.f90
> ./a.out
          74          98          86
> ./a.out
         113           3          10
> ./a.out
          44         104          27

Выглядит довольно случайно. Теперь компилирую с PGI Fortran и запускаю несколько раз:

> pgf90 testrand.f90
> ./a.out
            1            1            1
> ./a.out
            1            1            1
> ./a.out
            1            1            1

Конечно, нельзя быть полностью уверенным, но я подозреваю, что это не случайно. :) Кто-нибудь знает, что здесь происходит? Кто-нибудь знает правильный способ получения случайных чисел с помощью PGI Fortran?


person bob.sacamento    schedule 14.01.2016    source источник


Ответы (1)


Почему-то PGI не реализует system_clock, как в компиляторах GNU. Я не знаю почему, я нашел это недавно, делая такие же вещи, как и вы.

Чтобы понять, о чем я говорю, просто напечатайте ticks после вызова system_clock. Скорее всего, вы 0 все время получаете с PGI и разные числа с компиляторами GNU. Чтобы решить вашу проблему, вы можете адаптировать приведенный ниже код. Это слегка измененная версия кода, которую вы можете получить на веб-сайте GNU fortran

program testrand
use iso_fortran_env, only: int64

    integer, parameter :: nz = 160, nf = 160, nlt = 90
    real :: tmpidx(3)
    integer :: idxarr(3), idx1, idx2, idx3, seed_size, ticks
    integer, allocatable :: seed(:)

    call random_seed(size=seed_size)
    allocate(seed(seed_size))

    ! call system_clock(count=ticks)
    ! seed = ticks+37*(/(i-1, i=1,seed_size)/)
    ! call random_seed(put=seed)
    ! 
    ! deallocate(seed)

    call init_random_seed()

    call random_number(tmpidx)
    idxarr = tmpidx * (/nz, nf, nlt/)
    idx1 = max(1,idxarr(1))
    idx2 = max(1,idxarr(2))
    idx3 = max(1,idxarr(3))
    print *,idx1, idx2, idx3

contains
    !
    subroutine init_random_seed()
        implicit none
        integer, allocatable :: seed(:)
        integer :: i, n, istat, dt(8), pid
        integer(int64) :: t

        integer, parameter :: un=703

        call random_seed(size = n)
        allocate(seed(n))
        ! First try if the OS provides a random number generator
        open(unit=un, file="/dev/urandom", access="stream", &
            form="unformatted", action="read", status="old", iostat=istat)
        if (istat == 0) then
            read(un) seed
            close(un)
        else
            ! The PID is
            ! useful in case one launches multiple instances of the same
            ! program in parallel.
            call system_clock(t)
            if (t == 0) then
                call date_and_time(values=dt)
                t = (dt(1) - 1970) * 365_int64 * 24 * 60 * 60 * 1000 &
                + dt(2) * 31_int64 * 24 * 60 * 60 * 1000 &
                + dt(3) * 24_int64 * 60 * 60 * 1000 &
                + dt(5) * 60 * 60 * 1000 &
                + dt(6) * 60 * 1000 + dt(7) * 1000 &
                + dt(8)
            end if
            pid = getpid()
            t = ieor( t, int(pid, kind(t)) )
            do i = 1, n
                seed(i) = lcg(t)
            end do
        end if
        call random_seed(put=seed)
        !print*, "optimal seed = ", seed
    end subroutine init_random_seed
    !
    function lcg(s)
        integer :: lcg
        integer(int64), intent(in out) :: s
        if (s == 0) then
            s = 104729
        else
            s = mod(s, 4294967296_int64)
        end if
        s = mod(s * 279470273_int64, 4294967291_int64)
        lcg = int(mod(s, int(huge(0), 8)), kind(0))
    end function lcg
    !
    !this option is especially used for pgf90 to provide a getpid() function
    !> @brief Returns the process ID of the current process
    !! @todo write the actual code, for now returns a fixed value
    !<
    function getpid()result(pid)
        integer pid
        pid = 53 !just a prime number, no special meaning
    end function getpid
end program
person innoSPG    schedule 14.01.2016
comment
Ага! Я видел это, надеялся избежать этого! Тем не менее date_and_time (), похоже, работает намного лучше. Спасибо! - person bob.sacamento; 15.01.2016