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.