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

arrays performance

Question

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 que tiene aproximadamente 1 / 2G de tamaño. Necesito saber el mínimo, el máximo y la media de esta matriz y, como de todos modos pronto tendré que entender el código de Fortran, decidí probarlo 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 toma alrededor de 25 segundos por archivo en la máquina que uso. 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 realmente me impresionó. Se tarda menos de un segundo en condiciones idénticas. La media se desvía de la que encuentra mi rutina Fortran (que también ejecuté con flotantes de 128 bits, por lo que de alguna manera confío más) pero solo en el séptimo dígito significativo más o menos.

¿Cómo puede ser numpy tan rápido? Quiero decir que 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 me lleve mucho más tiempo?

EDITAR:

Para responder 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.
• iso_fortran_env que proporciona flotantes de 128 bits.
• Sin embargo, con flotadores de 32 bits, mi media está bastante apagada, por lo que la precisión es realmente un problema.
• Ejecuté ambas rutinas en diferentes archivos en diferente orden, por lo que supongo que el almacenamiento en caché debería haber sido justo en la comparación.
• En realidad intenté abrir MP, pero leer el archivo en diferentes posiciones al mismo tiempo. Después de leer 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 con las operaciones de matriz, pero tal vez eso ni siquiera sea necesario.
• Los archivos son en realidad 1 / 2G de tamaño, 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 numpy pero 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 sum es 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 computación (no mucho, pero notablemente). ¿Hay una mejor manera de solucionar este problema? No pude encontrar una manera de leer los singles del archivo directamente a los dobles. ¿Y cómo numpy evita esto?

Gracias por toda la ayuda hasta ahora.

Answer 1

El numpy es más rápido porque escribiste 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 toda la matriz de una vez y luego llama a las funciones que pueden operar en una matriz.

Mira tu código fortran. Lees un valor a la vez y haces algo de 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 verá que funciona 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

Answer 2

Su implementación de Fortran sufre dos deficiencias importantes:

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

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 tmp de una vez. Entonces, puedo usar las funciones MAXVAL , MINVAL y SUM directamente en la matriz.

Para el problema de precisión: simplemente usando valores de doble precisión y haciendo 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 segmentos, pero eso solo aumentó el tiempo requerido en el nivel de optimización predeterminado.

En -O3 , la adición de elementos sabios funciona ~ 3% mejor que la operación de matriz. La diferencia entre las operaciones de precisión doble y simple es menos del 2% en mi máquina, en promedio (las ejecuciones 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 usa la matriz de precisión simple 1-norma SLANGE en 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 precisión.