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

82

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.

user35915
fuente
10
¿Probó la rutina de Fortran sin flotadores de 128 bits? No conozco ningún hardware que realmente los admita, por lo que tendrían que hacerse en software.
user2357112 apoya a Monica el
4
¿Qué pasa si prueba la versión de Fortran usando una matriz (y en particular usando una lectura en lugar de mil millones)?
francescalus
9
¿También consideró usar operadores de matriz en Fortran? A continuación, puede probar minval(), maxval()y sum()? Además, está mezclando IO con las operaciones en Fortran, pero no en Python; esa no es una comparación justa ;-)
Alexander Vogt
4
Al comparar algo que involucre un archivo grande, asegúrese de que se almacene en caché de la misma manera para todas las ejecuciones.
Tom Zych
1
También tenga en cuenta que la precisión es un problema bastante importante en Fortran y tiene un costo. Incluso después de solucionar todos esos problemas obvios con su código de Fortran, es muy posible que se requiera una precisión adicional y cause una pérdida de velocidad significativa.
Luaan

Respuestas:

110

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.

Alexander Vogt
fuente
4
¿Por qué mezclar la entrada con el cálculo lo ralentiza tanto? Ambos tienen que leer el archivo completo, ese será el cuello de botella. Y si el sistema operativo lee con anticipación, el código de Fortran no debería tener que esperar mucho para la E / S.
Barmar
3
@Barmar Aún tendrá la sobrecarga de llamada de función y la lógica para verificar si los datos están en caché cada vez.
Overv el
55

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
casey
fuente
¿La media calculada de esta manera obtiene la misma precisión que numpyla .meanllamada de '? Tengo algunas dudas sobre eso.
Bakuriu
1
@Bakuriu No, no es así. Vea la respuesta de Alexander Vogt y mis ediciones sobre la pregunta.
user35915