Фортран. Генератор случайных вещественных чисел

[Sicker]

Как сделать генератор случайных вещественных чисел в фортране?
Функция random_number(r) выдает одно и тоже число $3,92\cdot 10^{-7}$

[Pphantom]

Sicker в сообщении #1068591 писал(а):

Как сделать генератор случайных вещественных чисел в фортране?
Функция random_number(r) выдает одно и тоже число $3,92\cdot 10^{-7}$

Минимальный рецепт - перед первым вызовом random_number выполнить

Используется синтаксис Fortran

 call random_seed

Для части компиляторов этого хватит. Если же результат все равно не меняется при разных запусках программы, можно использовать классический рецепт:

Используется синтаксис Fortran

subroutine init_random_seed()
  integer :: i, n, clock
  integer, dimension(:), allocatable :: seed

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

  call system_clock(count=clock)

  seed = clock + 37 * (/ (i - 1, i = 1, n) /)
  call random_seed(put = seed)

  deallocate(seed)
end subroutine init_random_seed

Тогда перед первым использованием датчика случайных чисел надо вызвать эту процедуру.

[Sicker]

Спасибо, уже сделал как в первом совете

[Gickle]

Доброго времени суток. У меня возникла похожая проблема. Мне нужно для моего массива (длины N) "самодельного" типа закинуть в поле coords случайные вектора длины 3.

Используется синтаксис Fortran

        do while (i<N+1)
        call random_number(coord)
        lattice(i)%coords = coord*L
        i=i+1
    end do

На деле же происходит какая-то чертовщина:
1) Практически (что понимается под практически - см. далее) все элементы получаются одинаковыми. При каждом новом запуске этот одинаковый почти для всех элементов вектор получается при этом различным.
2) Если вывести в конце, например, следующие элементы

Используется синтаксис Fortran

    print*, lattice(N)%coords
    print*, lattice(N-23)%coords
    print*, lattice(28)%coords

, то выяснится, что первые два равны друг другу, а третий отличается. Числа 23 и 28 не играют никакой роли и взяты для примера, N заведомо больше выбранных для примера чисел (~85000 в рассматриваемом случае). Причём, что я заметил, у первой пары и последнего первая координата (то есть lattice(...)%coords(1)) отличаются незначительно и всегда почти что на одну величину.
Да и в целом какая-то крайне хреновая случайность получается, что уж. Пытался пользоваться советами выше - не помогло. На всякий случай даже покоординатно вектора задавал - ничего, разумеется, не поменялось.

В чём моя ошибка, не подскажете?

[Pphantom]

Gickle в сообщении #1069961 писал(а):

В чём моя ошибка, не подскажете?

Я попробовал на всякий случай воспроизвести проблему, дописав недостающие участки кода, и мне это не удалось ни на одном из пяти имеющихся под руками компиляторов. Так что, полагаю, проблема где-то в другом месте.

Судя по проявлениям, у Вас там ошибки при использовании памяти (грубо говоря, Вы почему-то извлекаете данные не оттуда, откуда хотите), но хотелось бы посмотреть на код целиком. Если он небольшой - давайте сюда полностью, если большой - выложите куда-нибудь и давайте ссылку.

[Gickle]

Pphantom
Вот:

Используется синтаксис Fortran

    program frame
    implicit none
    integer N, p 
    integer :: i = 1
    real, parameter :: Pi = 3.141593
    real(8) L, phi, R_avg
    real(8), dimension(3)::coord
    type pore                                               
        real(8), dimension(3) :: coords
        real R
    end type
    type(pore), dimension(:), allocatable::lattice
    
    print*, 'L =' 
    read*, L    
    print*, 'phi ='
    read*, phi
    print*, 'R_avg ='
    read*, R_avg
    N = nint((L**3*log(1/(1-phi)))/((4*Pi*(R_avg)**3)/3))
    allocate(lattice(N))
    do while (i<N+1)
        call init_random_seed()
        call random_number(coord)
        lattice(i)%coords = coord*L
        i=i+1
    end do
    pause
    end program frame

Массив lattice я не очищаю, потому что по планам он дальше будет использован.

[Pphantom]

Ну что ж, видны как минимум две проблемы:
1) Зачем Вы вызов init_random_seed внутрь цикла вставили? Эту штуку нужно использовать один раз в начале программы.
2) Если N у Вас получается в районе 85000, то в стандартный для Вашего компилятора integer (который может оказаться и integer(2)) он может просто не влезть. Результаты попыток его туда втиснуть могут быть весьма разнообразны.

Итого, уберите инициализатор датчика случайных чисел в начало программы и смените целые типы на integer(4). Заодно проверьте результат вычисления N - что он соответствует ожиданиям.

Кстати, я бы и программу переписал, нынешний вариант жутковат, несмотря на небольшие размеры.

[Gickle]

Pphantom
Спасибо большое, помогло. Теперь вроде как всё работает нормально - как минимум среднее значение какой-либо координаты во всех прогнанных циклах оказалось очень близко к мат. ожиданию.

Цитата:

Зачем Вы вызов init_random_seed внутрь цикла вставили? Эту штуку нужно использовать один раз в начале программы.

Думал проверить, вдруг если каждый раз затравку менять, то всё исправится.

Цитата:

Заодно проверьте результат вычисления N - что он соответствует ожиданиям.

Да тут всё в порядке. К слову, как я понимаю, integer(k) отвечает тому, что предельное значение числа этого типа весить не больше k байт? То есть быть не больше $10^{8k-1}$? А то у меня в теории N должны быть сильно побольше.

Кстати, я бы и программу переписал, нынешний вариант жутковат, несмотря на небольшие размеры.

А что бы вы порекомендовали? Просто я с Фортраном только с пятницы знаком, если честно, а до этого весь опыт в программировании был ограничен баловством на матлабе (численные методы) по большей части. Так что и сам подозреваю, что всё жутко коряво.

[Pphantom]

Gickle в сообщении #1070024 писал(а):

Да тут всё в порядке. К слову, как я понимаю, integer(k) отвечает тому, что предельное значение числа этого типа весить не больше k байт?

В частности, и это тоже, хотя исходный смысл - просто количество байт, которое выделяется для хранения этого числа.

Gickle в сообщении #1070024 писал(а):

То есть быть не больше $10^{8k-1}$?

Только $2^{8\,k-1}$.

Gickle в сообщении #1070024 писал(а):

А то у меня в теории N должны быть сильно побольше.

Практически все современные компиляторы поддерживают integer(8). Надеюсь, этого хватит.

Gickle в сообщении #1070024 писал(а):

А что бы вы порекомендовали?

Первый банальный совет - воспользоваться циклом со счетчиком вместо его эрзаца, организованного вручную:

Используется синтаксис Fortran

do i=1,N
 содержимое цикла
end do

Второй, менее банальный (хотя и ненамного) совет - не использовать составные типы без явной и сильной необходимости. В данном случае данные проще хранить не в массиве записей, а в одном двумерном массиве, например, таком: real(8),dimension(0:3,1:N) :: lattice (считая, что первый индекс 0 соответствует полю R). Это, в частности, позволило бы вместо всего пресловутого цикла написать просто:

Используется синтаксис Fortran

 call random_number(lattice(1:3,:))
lattice(1:3,:)=L*lattice(1:3,:)

Такая конструкция, кстати, еще и работать быстрее будет.

[Gickle]

Pphantom

Цитата:

Только $2^{8\,k-1}$.

Ой, да, опечатался.

Цитата:

Практически все современные компиляторы поддерживают integer(8). Надеюсь, этого хватит.

Ну, я тут интересовался скорее из соображений "чтобы знать", нежели действительно беспокоился за то, что у меня памяти не хватит.

Цитата:

Первый банальный совет - воспользоваться циклом со счетчиком вместо его эрзаца, организованного вручную:

Это, вообще говоря, "пережиток прошлого", когда я проверял не совпадут ли координаты двух точек. Но решил, что раз это событие практически невероятное, то смысла в этом особого нет, а вот сменить while на do при этом забыл.

Цитата:

Второй, менее банальный (хотя и ненамного) совет - не использовать составные типы без явной и сильной необходимости. В данном случае данные проще хранить не в массиве записей, а в одном двумерном массиве, например, таком: real(8),dimension(0:3,1:N) :: lattice (считая, что первый индекс 0 соответствует полю R). Это, в частности, позволило бы вместо всего пресловутого цикла написать просто:

А насколько это в целом влияет на производительность, например? Просто это обусловлено в целом соображениями удобства, так как реально дальше полей будет больше и они несут в себе определённый физический смысл. Сильно ли выгоднее заменить массив составного типа на двумерный вида real(8),dimension(1:k,1:N) :: lattice? Конечно, удобство мне важно (я и так-то сильно путаюсь в фортране ), но производительность важнее, пожалуй.

[Pphantom]

Gickle в сообщении #1070042 писал(а):

А насколько это в целом влияет на производительность, например? Просто это обусловлено в целом соображениями удобства, так как реально дальше полей будет больше и они несут в себе определённый физический смысл. Сильно ли выгоднее заменить массив составного типа на двумерный вида real(8),dimension(1:k,1:N) :: lattice? Конечно, удобство мне важно (я и так-то сильно путаюсь в фортране ), но производительность важнее, пожалуй.

Использование массивов означает возможность использования операций над массивами. Собственно, выше как раз был пример этого: вместо организации цикла по блокам из трех чисел можно выполнить все необходимое сразу.

С точки зрения производительности это тоже выгоднее -компилятору удобнее обрабатывать групповые операции, поскольку они заведомо параллельны.

Если полей будет больше, то имеет смысл подумать о том, чтобы хранить именно отдельные массивы полей, а не массивы записей. В Фортране идеология сложных типов работает плохо (хотя и поддерживается), куда удобнее разнести по разным массивам разнотипные данные, относящиеся к однотипным физическим объектам, в этом случае выполнение однообразных операций с этими объектами резко упростится.