performance - una - multiplicacion de matrices en fortran 90

Rendimiento de la multiplicación de la matriz Fortran en diferentes optimizaciones

Question

Estoy leyendo el libro "Desarrollo de software científico con Fortran", y hay un ejercicio en él que me parece muy interesante:

"Cree un módulo Fortran llamado MatrixMultiplyModule. Agregue tres subrutinas, llamadas LoopMatrixMultiply, IntrinsicMatrixMultiply y MixMatrixMultiply. Cada rutina debe tomar dos matrices reales como argumento, realizar una multiplicación de matrices y devolver el resultado a través de un tercer argumento. LoopMatrixMultiply debe escribirse por completo con bucles do, y sin operaciones de matriz o procedimientos intrínsecos; IntrinsicMatrixMultiply debe escribirse utilizando la función intrínseca matmul, y MixMatrixMultiply debe escribirse utilizando algunos bucles do y la función intrínseca dot_product. Escribe un pequeño programa para probar el rendimiento de estas tres formas diferentes de realizar la multiplicación de la matriz para diferentes tamaños de matrices ".

Hice algunas pruebas de multiplicar dos matrices de rango 2 y aquí están los resultados, bajo diferentes indicadores de optimización:

compiler:ifort version 13.0.0 on Mac 

Aquí está mi pregunta:

¿Por qué bajo -O0 tienen el mismo rendimiento pero matmul tiene un gran aumento de rendimiento cuando se usa -O3, mientras que el producto explícito de bucle y punto tiene un menor impulso en el rendimiento? Además, ¿por qué dot_product parece tener el mismo rendimiento en comparación con los loops de do explícitos?

El código que uso es el siguiente:

module MatrixMultiplyModule

contains
    subroutine LoopMatrixMultiply(mtx1,mtx2,mtx3)
        real,intent(in)                 ::  mtx1(:,:),mtx2(:,:)
        real,intent(out),allocatable    ::  mtx3(:,:)
        integer                         ::  m,n
        integer                         ::  i,j
        if(size(mtx1,dim=2) /= size(mtx2,dim=1)) stop "input array size not match"
        m=size(mtx1,dim=1)
        n=size(mtx2,dim=2)
        allocate(mtx3(m,n))
        mtx3=0.

        do i=1,m
            do j=1,n
                do k=1,size(mtx1,dim=2)
                    mtx3(i,j)=mtx3(i,j)+mtx1(i,k)*mtx2(k,j)
                end do      
            end do
        end do
    end subroutine

    subroutine IntrinsicMatrixMultiply(mtx1,mtx2,mtx3)
        real,intent(in)                 ::  mtx1(:,:),mtx2(:,:)
        real,intent(out),allocatable    ::  mtx3(:,:)
        integer                         ::  m,n
        integer                         ::  i,j
        if(size(mtx1,dim=2) /= size(mtx2,dim=1)) stop "input array size not match"
        m=size(mtx1,dim=1)
        n=size(mtx2,dim=2)
        allocate(mtx3(m,n))
        mtx3=matmul(mtx1,mtx2)

    end subroutine

    subroutine MixMatrixMultiply(mtx1,mtx2,mtx3)
        real,intent(in)                 ::  mtx1(:,:),mtx2(:,:)
        real,intent(out),allocatable    ::  mtx3(:,:)
        integer                         ::  m,n
        integer                         ::  i,j
        if(size(mtx1,dim=2) /= size(mtx2,dim=1)) stop "input array size not match"
        m=size(mtx1,dim=1)
        n=size(mtx2,dim=2)
        allocate(mtx3(m,n))

        do i=1,m
            do j=1,n
                    mtx3(i,j)=dot_product(mtx1(i,:),mtx2(:,j))
            end do
        end do

    end subroutine

end module




program main
use MatrixMultiplyModule
implicit none

real,allocatable        ::  a(:,:),b(:,:)
real,allocatable    ::  c1(:,:),c2(:,:),c3(:,:)
integer ::  n
integer ::  count, rate
real    ::  timeAtStart, timeAtEnd
real    ::  time(3,10)
do n=100,1000,100
    allocate(a(n,n),b(n,n))

    call random_number(a)
    call random_number(b)

    call system_clock(count = count, count_rate = rate)
    timeAtStart = count / real(rate)
    call LoopMatrixMultiply(a,b,c1)
    call system_clock(count = count, count_rate = rate)
    timeAtEnd = count / real(rate)
    time(1,n/100)=timeAtEnd-timeAtStart

    call system_clock(count = count, count_rate = rate)
    timeAtStart = count / real(rate)
    call IntrinsicMatrixMultiply(a,b,c2)
    call system_clock(count = count, count_rate = rate)
    timeAtEnd = count / real(rate)
    time(2,n/100)=timeAtEnd-timeAtStart

    call system_clock(count = count, count_rate = rate)
    timeAtStart = count / real(rate)
    call MixMatrixMultiply(a,b,c3)
    call system_clock(count = count, count_rate = rate)
    timeAtEnd = count / real(rate)
    time(3,n/100)=timeAtEnd-timeAtStart


    deallocate(a,b)

end do

open(1,file="time.txt")
do n=1,10
    write(1,*) time(:,n)
end do
close(1)
deallocate(c1,c2,c3)
end program

Answer 1

Hay varias cosas que uno debe tener en cuenta al pasar sobre los elementos del conjunto:

• Asegúrese de que el bucle interno contenga elementos consecutivos en la memoria. En su algoritmo de "bucle" actual, el bucle interno está sobre el índice k. Dado que las matrices se presentan en la memoria como columnas (el primer índice varía más rápidamente al pasar por la memoria), el acceso a un nuevo valor de k podría necesitar cargar una nueva página en la memoria caché. En este caso, puede optimizar su algoritmo reordenando los bucles como sigue:

  do j=1,n
      do k=1,size(mtx1,dim=2)
          do i=1,m
              mtx3(i,j)=mtx3(i,j)+mtx1(i,k)*mtx2(k,j)
          end do      
      end do
  end do

  ahora, el ciclo interno está sobre elementos consecutivos en la memoria (el mtx2(k,j) del mtx2(k,j) probablemente será captado solo una vez antes del ciclo interno por el compilador, si no puede almacenarlo en una variable temporal antes del ciclo)

• Asegúrese de que todos los bucles puedan caber en la memoria caché para evitar demasiadas fallas de la memoria caché. Esto se puede hacer bloqueando el algoritmo. En este caso, una solución podría ser, por ejemplo:

  l=size(mtx1,dim=2)
  ichunk=512 ! I have a 3MB cache size (real*4)
  do jj=1,n,ichunk
      do kk=1,l,ichunk
  
  do j=jj,min(jj+ichunk-1,n)
      do k=kk,min(kk+ichunk-1,l)
          do i=1,m
              mtx3(i,j)=mtx3(i,j)+mtx1(i,k)*mtx2(k,j)
          end do      
      end do
  end do
  
      end do
  end do

  en cuyo caso el rendimiento dependerá del tamaño de ichunk , especialmente para matrices suficientemente grandes (incluso se podría bloquear el bucle interno, esto es solo un ejemplo).

• Asegúrese de que el trabajo necesario para realizar el ciclo sea mucho más pequeño que el trabajo dentro del ciclo. Esto se puede resolver mediante el ''desenrollado del bucle'', es decir, combinando varias declaraciones en una iteración del bucle. Por lo general, el compilador puede hacer esto suministrando el -funroll-loops .

Si uso el código anterior y compilo con los flags -O3 -funroll-loops , obtengo un rendimiento ligeramente mejor que con matmul .

Lo importante para recordar de esos tres es el primer punto sobre el orden de bucles, ya que esto es algo que afectará el rendimiento en otros casos de uso, y el compilador generalmente no puede arreglar eso. El desenrollado del bucle, puede dejarlo en el compilador (pero pruébelo, ya que esto no siempre aumenta el rendimiento). En cuanto al segundo punto, dado que esto depende del hardware, usted no debe (generalmente) tratar de implementar una multiplicación de matrices muy eficiente usted mismo y en su lugar considerar usar una biblioteca como p. Ej. Atlas, que puede optimizar el tamaño de la caché, o biblioteca de proveedores como MKL o ACML.