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

Question 1

시뮬레이션 (Fortran으로 작성)에서 온도 분포를 나타내는 512 ^ 3 배열을 얻습니다. 어레이는 약 1 / 2G 크기의 이진 파일에 저장됩니다. 이 배열의 최소, 최대 및 평균을 알아야합니다. 어쨌든 포트란 코드를 곧 이해할 필요가있을 것이므로 시도하기로 결정하고 다음과 같은 매우 쉬운 루틴을 생각해 냈습니다.

  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이것을 어떻게 피합니까?

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

Question 2

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

이것은 SLANGE행렬 열 에서 단 정밀도 행렬 1- 노름 을 사용합니다 . 런타임은 단 정밀도 배열 함수를 사용하는 접근 방식보다 훨씬 빠르며 정밀도 문제를 나타내지 않습니다.

Question 3

numpy는 파이썬에서 훨씬 더 효율적인 코드를 작성하고 (그리고 numpy 백엔드의 대부분은 최적화 된 Fortran 및 C로 작성 됨) Fortran에서 매우 비효율적 인 코드를 작성했기 때문에 더 빠릅니다.

파이썬 코드를보십시오. 전체 배열을 한 번에로드 한 다음 배열에서 작동 할 수있는 함수를 호출합니다.

포트란 코드를보십시오. 한 번에 하나의 값을 읽고이 값으로 분기 논리를 수행합니다.

대부분의 불일치는 Fortran에서 작성한 조각난 IO입니다.

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

답변

답변

답변