Как numpy может быть намного быстрее, чем моя программа Fortran?

Я получаю массив 512 ^ 3, представляющий распределение температуры от моделирования (написанное на языке Fortran). Массив хранится в двоичном файле размером около 1/2G. Мне нужно знать минимальный, максимальный и средний из этого массива, и поскольку мне скоро понадобится понять код Fortran, я решил дать ему повод и придумал следующую очень легкую процедуру.

  integer gridsize,unit,j
  real mini,maxi
  double precision mean

  gridsize=512
  unit=40
  open(unit=unit,file='T.out',status='old',access='stream',&
       form='unformatted',action='read')
  read(unit=unit) tmp
  mini=tmp
  maxi=tmp
  mean=tmp
  do j=2,gridsize**3
      read(unit=unit) tmp
      if(tmp>maxi)then
          maxi=tmp
      elseif(tmp<mini)then
          mini=tmp
      end if
      mean=mean+tmp
  end do
  mean=mean/gridsize**3
  close(unit=unit)

Это занимает около 25 секунд на файл на машине, которую я использую. Это показалось мне довольно длинным, поэтому я пошел дальше и сделал следующее на Python:

    import numpy

    mmap=numpy.memmap('T.out',dtype='float32',mode='r',offset=4,\
                                  shape=(512,512,512),order='F')
    mini=numpy.amin(mmap)
    maxi=numpy.amax(mmap)
    mean=numpy.mean(mmap)

Теперь я ожидал, что это будет быстрее, но я действительно сдулся. В то же время требуется меньше секунды. Среднее отклоняется от одного моего подпрограммы "Fortran" (который я также запускал с 128-битными поплавками, поэтому я как-то доверяю ему больше), но только на 7-й значащей цифре или около того.

Как numpy может быть настолько быстрым? Я имею в виду, что вы должны смотреть на каждую запись массива, чтобы найти эти значения, не так ли? Я делаю что-то очень глупо в моей программе Fortran, чтобы она занимала намного больше времени?

EDIT:

Чтобы ответить на вопросы в комментариях:

  • Да, я также запускал подпрограмму Fortran с 32-битными и 64-битными поплавками, но не влиял на производительность.
  • Я использовал iso_fortran_env, который обеспечивает 128-битные поплавки.
  • Используя 32-битные поплавки, мое среднее значение довольно велико, поэтому точность действительно является проблемой.
  • Я запускал обе процедуры в разных файлах в другом порядке, поэтому кеширование должно быть справедливым в сравнении, я думаю?
  • Я на самом деле пробовал открытый MP, но читал из файла на разных позициях одновременно. Прочитав ваши комментарии и ответы, это звучит действительно глупо сейчас, и это заставило процедуру заняться намного дольше. Я мог бы дать ему попробовать массивные операции, но, возможно, это даже не понадобится.
  • Файлы на самом деле размером 1/2G, это была опечатка, спасибо.
  • Теперь я попробую реализацию массива.

ИЗМЕНИТЬ 2:

Я реализовал то, что @Alexander Vogt и @casey предложили в своих ответах, и это так же быстро, как numpy, но теперь у меня есть проблема с точностью, которую, как заметил @Luaan, я могу получить. Используя 32-битный массив с плавающей точкой, среднее значение, вычисленное на sum, составляет 20%. Выполнение

...
real,allocatable :: tmp (:,:,:)
double precision,allocatable :: tmp2(:,:,:)
...
tmp2=tmp
mean=sum(tmp2)/size(tmp)
...

Решает проблему, но увеличивает вычислительное время (не очень, но заметно). Есть ли лучший способ обойти эту проблему? Я не мог найти способ прочитать синглы из файла непосредственно в парном разряде. И как numpy избежать этого?

Спасибо за всю помощь до сих пор.

Ответы

Ответ 1

В вашей реализации Fortran есть два основных недостатка:

  • Вы смешиваете IO и вычисления (и читаете из записи файла по записи).
  • Вы не используете векторные/матричные операции.

Эта реализация выполняет ту же операцию, что и ваша, и работает быстрее на 20 раз на моей машине:

program test
  integer gridsize,unit
  real mini,maxi,mean
  real, allocatable :: tmp (:,:,:)

  gridsize=512
  unit=40

  allocate( tmp(gridsize, gridsize, gridsize))

  open(unit=unit,file='T.out',status='old',access='stream',&
       form='unformatted',action='read')
  read(unit=unit) tmp

  close(unit=unit)

  mini = minval(tmp)
  maxi = maxval(tmp)
  mean = sum(tmp)/gridsize**3
  print *, mini, maxi, mean

end program

Идея состоит в том, чтобы прочитать весь файл в один массив tmp за один раз. Затем я могу использовать функции MAXVAL, MINVAL и SUM непосредственно в массиве.

Для проблемы точности: просто используйте значения двойной точности и выполняйте преобразование "на лету" как

mean = sum(real(tmp, kind=kind(1.d0)))/real(gridsize**3, kind=kind(1.d0))

лишь незначительно увеличивает время вычисления. Я попытался выполнить операцию с элементами и в срезах, но это только увеличило требуемое время на уровне оптимизации по умолчанию.

В -O3 элементное добавление выполняет на 3% лучше, чем операция массива. Разница между операциями двойной и одиночной точности составляет менее 2% на моей машине - в среднем (индивидуальные прогоны отклоняются намного больше).

Вот очень быстрая реализация с использованием LAPACK:

program test
  integer gridsize,unit, i, j
  real mini,maxi
  integer  :: t1, t2, rate
  real, allocatable :: tmp (:,:,:)
  real, allocatable :: work(:)
!  double precision :: mean
  real :: mean
  real :: slange

  call system_clock(count_rate=rate)
  call system_clock(t1)
  gridsize=512
  unit=40

  allocate( tmp(gridsize, gridsize, gridsize), work(gridsize))

  open(unit=unit,file='T.out',status='old',access='stream',&
       form='unformatted',action='read')
  read(unit=unit) tmp

  close(unit=unit)

  mini = minval(tmp)
  maxi = maxval(tmp)

!  mean = sum(tmp)/gridsize**3
!  mean = sum(real(tmp, kind=kind(1.d0)))/real(gridsize**3, kind=kind(1.d0))
  mean = 0.d0
  do j=1,gridsize
    do i=1,gridsize
      mean = mean + slange('1', gridsize, 1, tmp(:,i,j),gridsize, work)
    enddo !i
  enddo !j
  mean = mean / gridsize**3

  print *, mini, maxi, mean
  call system_clock(t2)
  print *,real(t2-t1)/real(rate)

end program

Это использует матрицу однократной точности 1-norm SLANGE для матричных столбцов. Время выполнения еще быстрее, чем подход, использующий функции одиночной точности, и не показывает проблему точности.

Ответ 2

Целое число быстрее, потому что вы пишете гораздо более эффективный код в python (и большая часть базы данных numpy написана в оптимизированном Fortran и C) и ужасно неэффективный код в Fortran.

Посмотрите на свой код на Python. Вы загружаете весь массив сразу, а затем вызываете функции, которые могут работать с массивом.

Посмотрите на свой код fortran. Вы читаете одно значение за раз и делаете с ним некоторую ветвящуюся логику.

Большая часть вашего несоответствия - фрагментированный IO, написанный вами в Fortran.

Вы можете написать Fortran примерно так же, как написано на питоне, и вы обнаружите, что он работает намного быстрее таким образом.

program test
  implicit none
  integer :: gridsize, unit
  real :: mini, maxi, mean
  real, allocatable :: array(:,:,:)

  gridsize=512
  allocate(array(gridsize,gridsize,gridsize))
  unit=40
  open(unit=unit, file='T.out', status='old', access='stream',&
       form='unformatted', action='read')
  read(unit) array    
  maxi = maxval(array)
  mini = minval(array)
  mean = sum(array)/size(array)
  close(unit)
end program test