¿C es más lento que Fortran en el tiroteo de la norma espectral (usando gcc, intel y otros compiladores)?

La conclusión aquí:

¿Cuánto mejor son realmente los compiladores de Fortran?

es que gfortran y gcc son tan rápidos para un código simple. Entonces quería probar algo más complicado. Tomé el ejemplo de tiroteo de la norma espectral. Primero precalculo la matriz 2D A (:, :), y luego calculo la norma. (Creo que esta solución no está permitida en el tiroteo). He implementado la versión Fortran y C. Aquí está el código:

https://github.com/certik/spectral_norm

Las versiones más rápidas de gfortran son spectral_norm2.f90 y spectral_norm6.f90 (una usa el matmul incorporado y dot_product de Fortran, la otra implementa estas dos funciones en el código, sin diferencia de velocidad). El código C / C ++ más rápido que pude escribir es spectral_norm7.cpp. Los tiempos a partir de la versión git 457d9d9 en mi computadora portátil son:

$ time ./spectral_norm6 5500
1.274224153

real    0m2.675s
user    0m2.520s
sys 0m0.132s


$ time ./spectral_norm7 5500
1.274224153

real    0m2.871s
user    0m2.724s
sys 0m0.124s

Entonces la versión de gfortran es un poco más rápida. ¿Porqué es eso? Si envía una solicitud de extracción con una implementación de C más rápida (o simplemente pega un código), actualizaré el repositorio.

En Fortran paso una matriz 2D, mientras que en CI uso una matriz 1D. Siéntase libre de usar una matriz 2D o cualquier otra forma que le parezca adecuada.

En cuanto a los compiladores, comparemos gcc vs gfortran, icc vs ifort, etc. (A diferencia de la página de tiroteo, que compara ifort vs gcc.)

Actualización : usando la versión 179dae2, que mejora matmul3 () en mi versión C, ahora son tan rápidos:

$ time ./spectral_norm6 5500
1.274224153

real    0m2.669s
user    0m2.500s
sys 0m0.144s

$ time ./spectral_norm7 5500
1.274224153

real    0m2.665s
user    0m2.472s
sys 0m0.168s

La versión vectorizada de Pedro a continuación es más rápida:

$ time ./spectral_norm8 5500
1.274224153

real    0m2.523s
user    0m2.336s
sys 0m0.156s

Finalmente, como laxxy informa a continuación para los compiladores Intel, no parece haber una gran diferencia allí e incluso el código Fortran más simple (spectral_norm1) se encuentra entre los más rápidos.

En primer lugar, ¡gracias por publicar esta pregunta / desafío! Como descargo de responsabilidad, soy un programador nativo de C con algo de experiencia en Fortran, y me siento como en casa en C, por lo que me centraré solo en mejorar la versión C. ¡Invito a todos los piratas de Fortran a probar también!

Solo para recordarles a los recién llegados de qué se trata: la premisa básica en este hilo era que gcc / fortran e icc / ifort deberían, dado que tienen los mismos back-end respectivamente, producir código equivalente para el mismo programa (semánticamente idéntico), independientemente de que sea en C o Fortran. La calidad del resultado depende solo de la calidad de las implementaciones respectivas.

Jugué un poco con el código y en mi computadora (ThinkPad 201x, Intel Core i5 M560, 2.67 GHz), usando gcc4.6.1 y los siguientes indicadores de compilación:

GCCFLAGS= -O3 -g -Wall -msse2 -march=native -funroll-loops -ffast-math -fomit-frame-pointer -fstrict-aliasing

También seguí adelante y escribí una versión SIMD-vectorizada en lenguaje C del código C ++ spectral_norm_vec.c:

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>

/* Define the generic vector type macro. */  
#define vector(elcount, type)  __attribute__((vector_size((elcount)*sizeof(type)))) type

double Ac(int i, int j)
{
    return 1.0 / ((i+j) * (i+j+1)/2 + i+1);
}

double dot_product2(int n, double u[], double v[])
{
    double w;
    int i;
    union {
        vector(2,double) v;
        double d[2];
        } *vu = u, *vv = v, acc[2];

    /* Init some stuff. */
    acc[0].d[0] = 0.0; acc[0].d[1] = 0.0;
    acc[1].d[0] = 0.0; acc[1].d[1] = 0.0;

    /* Take in chunks of two by two doubles. */
    for ( i = 0 ; i < (n/2 & ~1) ; i += 2 ) {
        acc[0].v += vu[i].v * vv[i].v;
        acc[1].v += vu[i+1].v * vv[i+1].v;
        }
    w = acc[0].d[0] + acc[0].d[1] + acc[1].d[0] + acc[1].d[1];

    /* Catch leftovers (if any) */
    for ( i = n & ~3 ; i < n ; i++ )
        w += u[i] * v[i];

    return w;

}

void matmul2(int n, double v[], double A[], double u[])
{
    int i, j;
    union {
        vector(2,double) v;
        double d[2];
        } *vu = u, *vA, vi;

    bzero( u , sizeof(double) * n );

    for (i = 0; i < n; i++) {
        vi.d[0] = v[i];
        vi.d[1] = v[i];
        vA = &A[i*n];
        for ( j = 0 ; j < (n/2 & ~1) ; j += 2 ) {
            vu[j].v += vA[j].v * vi.v;
            vu[j+1].v += vA[j+1].v * vi.v;
            }
        for ( j = n & ~3 ; j < n ; j++ )
            u[j] += A[i*n+j] * v[i];
        }

}


void matmul3(int n, double A[], double v[], double u[])
{
    int i;

    for (i = 0; i < n; i++)
        u[i] = dot_product2( n , &A[i*n] , v );

}

void AvA(int n, double A[], double v[], double u[])
{
    double tmp[n] __attribute__ ((aligned (16)));
    matmul3(n, A, v, tmp);
    matmul2(n, tmp, A, u);
}


double spectral_game(int n)
{
    double *A;
    double u[n] __attribute__ ((aligned (16)));
    double v[n] __attribute__ ((aligned (16)));
    int i, j;

    /* Aligned allocation. */
    /* A = (double *)malloc(n*n*sizeof(double)); */
    if ( posix_memalign( (void **)&A , 4*sizeof(double) , sizeof(double) * n * n ) != 0 ) {
        printf( "spectral_game:%i: call to posix_memalign failed.\n" , __LINE__ );
        abort();
        }


    for (i = 0; i < n; i++) {
        for (j = 0; j < n; j++) {
            A[i*n+j] = Ac(i, j);
        }
    }


    for (i = 0; i < n; i++) {
        u[i] = 1.0;
    }
    for (i = 0; i < 10; i++) {
        AvA(n, A, u, v);
        AvA(n, A, v, u);
    }
    free(A);
    return sqrt(dot_product2(n, u, v) / dot_product2(n, v, v));
}

int main(int argc, char *argv[]) {
    int i, N = ((argc >= 2) ? atoi(argv[1]) : 2000);
    for ( i = 0 ; i < 10 ; i++ )
        printf("%.9f\n", spectral_game(N));
    return 0;
}

Las tres versiones fueron compiladas con las mismas banderas y la misma gccversión. Tenga en cuenta que envolví la llamada a la función principal en un bucle de 0..9 para obtener tiempos más precisos.

$ time ./spectral_norm6 5500
1.274224153
...
real    0m22.682s
user    0m21.113s
sys 0m1.500s

$ time ./spectral_norm7 5500
1.274224153
...
real    0m21.596s
user    0m20.373s
sys 0m1.132s

$ time ./spectral_norm_vec 5500
1.274224153
...
real    0m21.336s
user    0m19.821s
sys 0m1.444s

Entonces, con indicadores de compilador "mejores", la versión de C ++ supera a la versión de Fortran y los bucles vectorizados codificados a mano solo proporcionan una mejora marginal. Un vistazo rápido al ensamblador para la versión C ++ muestra que los bucles principales también se han vectorizado, aunque se han desenrollado de forma más agresiva.

También eché un vistazo al ensamblador generado por gfortrany aquí está la gran sorpresa: sin vectorización. Atribuyo el hecho de que es solo un poco más lento al problema de que el ancho de banda es limitado, al menos en mi arquitectura. Para cada una de las multiplicaciones de la matriz, se atraviesan 230 MB de datos, lo que prácticamente inunda todos los niveles de caché. Si utiliza un valor de entrada más pequeño, por ejemplo 100, las diferencias de rendimiento aumentan considerablemente.

Como nota al margen, en lugar de obsesionarse con las banderas de vectorización, alineación y compilación, la optimización más obvia sería calcular las primeras iteraciones en aritmética de precisión simple, hasta que tengamos ~ 8 dígitos del resultado. Las instrucciones de precisión simple no solo son más rápidas, sino que la cantidad de memoria que debe moverse también se reduce a la mitad.

La respuesta del usuario 389 ha sido eliminada, pero permítanme decir que estoy firmemente en su campo: no veo lo que aprendemos al comparar micro-puntos de referencia en diferentes idiomas. No me sorprende tanto que C y Fortran obtengan el mismo rendimiento en este punto de referencia dado lo corto que es. Pero el punto de referencia también es aburrido, ya que se puede escribir fácilmente en ambos idiomas en un par de docenas de líneas. Desde el punto de vista del software, ese no es un caso representativo: deberíamos preocuparnos por el software que tiene 10,000 o 100,000 líneas de código y cómo los compiladores lo hacen. Por supuesto, a esa escala, uno descubrirá rápidamente otras cosas: que el lenguaje A requiere 10,000 líneas mientras que el lenguaje B requiere 50,000. O al revés, dependiendo de lo que quieras hacer. Y de repente '

En otras palabras, no me importa mucho que tal vez mi aplicación podría ser un 50% más rápida si la desarrollé en Fortran 77 si en cambio solo me tomaría 1 mes para que funcione correctamente, mientras que me llevaría 3 meses en F77. El problema con la pregunta aquí es que se centra en un aspecto (núcleos individuales) que, en mi opinión, no es relevante en la práctica.

Resulta que puedo escribir un código Python (usando numpy para hacer las operaciones BLAS) más rápido que el código Fortran compilado con el compilador gfortran de mi sistema.

$ gfortran -o sn6a sn6a.f90 -O3 -march=native
    
    $ ./sn6a 5500
1.274224153
1.274224153
1.274224153
   1.9640001      sec per iteration

$ python ./foo1.py
1.27422415279
1.27422415279
1.27422415279
1.20618661245 sec per iteration

foo1.py:

import numpy
import scipy.linalg
import timeit

def specNormDot(A,n):
    u = numpy.ones(n)
    v = numpy.zeros(n)

    for i in xrange(10):
        v  = numpy.dot(numpy.dot(A,u),A)
        u  = numpy.dot(numpy.dot(A,v),A)

    print numpy.sqrt(numpy.vdot(u,v)/numpy.vdot(v,v))

    return

n = 5500

ii, jj = numpy.meshgrid(numpy.arange(1,n+1), numpy.arange(1,n+1))
A  = (1./((ii+jj-2.)*(ii+jj-1.)/2. + ii))

t = timeit.Timer("specNormDot(A,n)", "from __main__ import specNormDot,A,n")
ntries = 3

print t.timeit(ntries)/ntries, "sec per iteration"

y sn6a.f90, un spectral_norm6.f90 muy ligeramente modificado:

program spectral_norm6
! This uses spectral_norm3 as a starting point, but does not use the
! Fortrans
! builtin matmul and dotproduct (to make sure it does not call some
! optimized
! BLAS behind the scene).
implicit none

integer, parameter :: dp = kind(0d0)
real(dp), allocatable :: A(:, :), u(:), v(:)
integer :: i, j, n
character(len=6) :: argv
integer :: calc, iter
integer, parameter :: niters=3

call get_command_argument(1, argv)
read(argv, *) n

allocate(u(n), v(n), A(n, n))
do j = 1, n
    do i = 1, n
        A(i, j) = Ac(i, j)
    end do
end do

call tick(calc)

do iter=1,niters
    u = 1
    do i = 1, 10
        v = AvA(A, u)
        u = AvA(A, v)
    end do

    write(*, "(f0.9)") sqrt(dot_product2(u, v) / dot_product2(v, v))
enddo

print *, tock(calc)/niters, ' sec per iteration'

contains

pure real(dp) function Ac(i, j) result(r)
integer, intent(in) :: i, j
r = 1._dp / ((i+j-2) * (i+j-1)/2 + i)
end function

pure function matmul2(v, A) result(u)
! Calculates u = matmul(v, A), but much faster (in gfortran)
real(dp), intent(in) :: v(:), A(:, :)
real(dp) :: u(size(v))
integer :: i
do i = 1, size(v)
    u(i) = dot_product2(A(:, i), v)
end do
end function

pure real(dp) function dot_product2(u, v) result(w)
! Calculates w = dot_product(u, v)
real(dp), intent(in) :: u(:), v(:)
integer :: i
w = 0
do i = 1, size(u)
    w = w + u(i)*v(i)
end do
end function

pure function matmul3(A, v) result(u)
! Calculates u = matmul(v, A), but much faster (in gfortran)
real(dp), intent(in) :: v(:), A(:, :)
real(dp) :: u(size(v))
integer :: i, j
u = 0
do j = 1, size(v)
    do i = 1, size(v)
        u(i) = u(i) + A(i, j)*v(j)
    end do
end do
end function

pure function AvA(A, v) result(u)
! Calculates u = matmul2(matmul3(A, v), A)
! In gfortran, this function is sligthly faster than calling
! matmul2(matmul3(A, v), A) directly.
real(dp), intent(in) :: v(:), A(:, :)
real(dp) :: u(size(v))
u = matmul2(matmul3(A, v), A)
end function

subroutine tick(t)
    integer, intent(OUT) :: t

    call system_clock(t)
end subroutine tick

! returns time in seconds from now to time described by t 
real function tock(t)
    integer, intent(in) :: t
    integer :: now, clock_rate

    call system_clock(now,clock_rate)

    tock = real(now - t)/real(clock_rate)
end function tock
end program

Lo comprobé con compiladores Intel. Con 11.1 (-fast, lo que implica -O3), y con 12.0 (-O2), los más rápidos son 1,2,6,7 y 8 (es decir, los códigos "más simples" de Fortran y C, y el C vectorizado a mano) - estos son indistinguibles entre sí a ~ 1.5s. Las pruebas 3 y 5 (con la matriz como función) son más lentas; # 4 No pude compilar.

En particular, si compila con 12.0 y -O3, en lugar de -O2, los primeros 2 ("más simples") códigos Fortran se ralentizan MUCHO (1.5 -> 10.2 segundos) - esta no es la primera vez que veo algo así esto, pero este puede ser el ejemplo más dramático. Si este sigue siendo el caso en la versión actual, creo que sería una buena idea informarlo a Intel, ya que claramente hay algo que va muy mal con sus optimizaciones en este caso bastante simple.

De lo contrario, estoy de acuerdo con Jonathan en que este no es un ejercicio particularmente informativo :)