¿Cómo puede numpy ser mucho más rápido que mi rutina de Fortran?

Obtengo una matriz de 512 ^ 3 que representa una distribución de temperatura de una simulación (escrita en Fortran). La matriz se almacena en un archivo binario de aproximadamente 1 / 2G de tamaño. Necesito saber el mínimo, el máximo y la media de esta matriz y, como pronto necesitaré comprender el código de Fortran de todos modos, decidí intentarlo y se me ocurrió la siguiente rutina muy fácil.

  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)

Esto lleva unos 25 segundos por archivo en la máquina que utilizo. Eso me pareció bastante largo, así que seguí adelante e hice lo siguiente en 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)

Ahora, esperaba que esto fuera más rápido, por supuesto, pero estaba realmente impresionado. Tarda menos de un segundo en condiciones idénticas. La media se desvía de la que encuentra mi rutina de Fortran (que también ejecuté con flotantes de 128 bits, así que de alguna manera confío más en ella) pero solo en el séptimo dígito significativo más o menos.

¿Cómo puede ser tan rápido Numpy? Quiero decir, tienes que mirar cada entrada de una matriz para encontrar estos valores, ¿verdad? ¿Estoy haciendo algo muy estúpido en mi rutina de Fortran para que tarde tanto tiempo?

EDITAR:

Para responder a las preguntas en los comentarios:

• Sí, también ejecuté la rutina Fortran con flotadores de 32 y 64 bits, pero no tuvo ningún impacto en el rendimiento.
• Usé iso_fortran_envque proporciona flotantes de 128 bits.
• Sin embargo, al usar flotadores de 32 bits, mi promedio está bastante desviado, por lo que la precisión es realmente un problema.
• Ejecuté ambas rutinas en archivos diferentes en un orden diferente, por lo que el almacenamiento en caché debería haber sido justo en la comparación, supongo.
• De hecho, intenté abrir MP, pero para leer el archivo en diferentes posiciones al mismo tiempo. Habiendo leído sus comentarios y respuestas, esto suena realmente estúpido ahora y también hizo que la rutina tomara mucho más tiempo. Podría intentarlo en las operaciones de la matriz, pero tal vez ni siquiera sea necesario.
• Los archivos tienen un tamaño de 1 / 2G, eso fue un error tipográfico, gracias.
• Probaré la implementación de la matriz ahora.

EDITAR 2:

Implementé lo que @Alexander Vogt y @casey sugirieron en sus respuestas, y es tan rápido como, numpypero ahora tengo un problema de precisión como @Luaan señaló que podría tener. Usando una matriz flotante de 32 bits, la media calculada por sumes un 20% de descuento. Haciendo

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

Resuelve el problema pero aumenta el tiempo de cálculo (no mucho, pero sí de forma notable). ¿Existe una mejor manera de solucionar este problema? No pude encontrar una manera de leer sencillos del archivo directamente en dobles. ¿Y cómo numpyevita esto?

Gracias por toda la ayuda hasta ahora.

Su implementación de Fortran adolece de dos deficiencias importantes:

• Mezcla IO y cálculos (y lee del archivo entrada por entrada).
• No utiliza operaciones de vector / matriz.

Esta implementación realiza la misma operación que la suya y es más rápida por un factor de 20 en mi máquina:

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

La idea es leer todo el archivo en una matriz tmpde una sola vez. Entonces, puedo usar las funciones MAXVAL, MINVALy SUMen la matriz directamente.

Para el problema de la precisión: simplemente use valores de doble precisión y realice la conversión sobre la marcha como

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

solo aumenta marginalmente el tiempo de cálculo. Intenté realizar la operación por elementos y en porciones, pero eso solo aumentó el tiempo requerido en el nivel de optimización predeterminado.

En -O3, la adición de elementos funciona un 3% mejor que la operación de matriz. La diferencia entre las operaciones de precisión simple y doble es menos del 2% en mi máquina, en promedio (las corridas individuales se desvían mucho más).

Aquí hay una implementación muy rápida usando 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

Esto utiliza la norma de matriz 1 de precisión simple SLANGEen columnas de matriz. El tiempo de ejecución es incluso más rápido que el enfoque que utiliza funciones de matriz de precisión única, y no muestra el problema de la precisión.

El numpy es más rápido porque escribió un código mucho más eficiente en Python (y gran parte del backend numpy está escrito en Fortran y C optimizados) y un código terriblemente ineficiente en Fortran.

Mira tu código de Python. Carga la matriz completa a la vez y luego llama a funciones que pueden operar en una matriz.

Mira tu código de fortran. Lees un valor a la vez y haces una lógica de ramificación con él.

La mayor parte de su discrepancia es el IO fragmentado que ha escrito en Fortran.

Puede escribir el Fortran casi de la misma manera que escribió el pitón y encontrará que corre mucho más rápido de esa manera.

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