어떻게 numpy가 내 Fortran 루틴보다 훨씬 빠를 수 있습니까?

어떻게 numpy가 내 Fortran 루틴보다 훨씬 빠를 수 있습니까?

---

시뮬레이션 (Fortran으로 작성)에서 온도 분포를 나타내는 512 ^ 3 배열을 얻습니다. 어레이는 약 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)

당연히 더 빨라질 거라고 생각했지만 정말 놀랐습니다. 동일한 조건에서 1 초도 걸리지 않습니다. 평균은 내 Fortran 루틴이 찾은 것 (128 비트 부동 소수점으로 실행되었으므로 어떻게 든 더 많이 신뢰 함)과는 다르지만 7 번째 유효 숫자 정도에서만 나타납니다.

numpy가 어떻게 그렇게 빠를 수 있습니까? 이 값을 찾으려면 배열의 모든 항목을 살펴 봐야합니다. 그렇죠? Fortran 루틴에서 너무 오래 걸리기 위해 매우 어리석은 일을하고 있습니까?

편집하다:

댓글의 질문에 답하려면 :

• 예, 또한 32 비트 및 64 비트 부동으로 Fortran 루틴을 실행했지만 성능에 영향을 미치지 않았습니다.
• 나는 iso_fortran_env128 비트 수레를 제공하는 것을 사용했습니다 .
• 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이것을 어떻게 피합니까?

지금까지 도움을 주셔서 감사합니다.

---

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))

only marginally increases the calculation time. I tried performing the operation element-wise and in slices, but that did only increase the required time at the default optimization level.

At -O3, the element-wise addition performs ~3 % better than the array operation. The difference between double and single precision operations is less than 2% on my machine - on average (the individual runs deviate by far more).

---

Here is a very fast implementation using 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

This uses the single precision matrix 1-norm SLANGE on matrix columns. The run-time is even faster than the approach using single precision array functions - and does not show the precision issue.

---

The numpy is faster because you wrote much more efficient code in python (and much of the numpy backend is written in optimized Fortran and C) and terribly inefficient code in Fortran.

Look at your python code. You load the entire array at once and then call functions that can operate on an array.

Look at your fortran code. You read one value at a time and do some branching logic with it.

The majority of your discrepancy is the fragmented IO you have written in Fortran.

You can write the Fortran just about the same way as you wrote the python and you'll find it runs much faster that way.

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

참고URL : https://stackoverflow.com/questions/33723771/how-can-numpy-be-so-much-faster-than-my-fortran-routine