C est-il plus lent que Fortran sur la fusillade de la norme spectrale (en utilisant gcc, intel et d'autres compilateurs)?

La conclusion ici:

Les compilateurs Fortran sont-ils vraiment meilleurs?

est que gfortran et gcc sont aussi rapides pour le code simple. Je voulais donc essayer quelque chose de plus compliqué. J'ai pris l'exemple de la fusillade de la norme spectrale. Je précalcule d'abord la matrice 2D A (:, :), puis calcule la norme. (Je pense que cette solution n'est pas autorisée lors de la fusillade.) J'ai implémenté la version Fortran et C. Voici le code:

https://github.com/certik/spectral_norm

Les versions les plus rapides de gfortran sont spectral_norm2.f90 et spectral_norm6.f90 (l'une utilise les matmul et dot_product intégrés de Fortran, l'autre implémente ces deux fonctions dans le code - sans différence de vitesse). Le code C / C ++ le plus rapide que j'ai pu écrire est spectral_norm7.cpp. Les horaires à partir de la version git 457d9d9 sur mon ordinateur portable sont les suivants:

$ 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

La version de gfortran est donc un peu plus rapide. Pourquoi donc? Si vous envoyez une demande d'extraction avec une implémentation C plus rapide (ou collez simplement un code), je mettrai à jour le référentiel.

Dans Fortran, je passe un tableau 2D, tandis que dans CI, j'utilise un tableau 1D. N'hésitez pas à utiliser un tableau 2D ou tout autre moyen que vous jugerez utile.

Quant aux compilateurs, comparons gcc vs gfortran, icc vs ifort, et ainsi de suite. (Contrairement à la page de shootout, qui compare ifort vs gcc.)

Mise à jour : en utilisant la version 179dae2, qui améliore matmul3 () dans ma version C, ils sont désormais aussi rapides:

$ 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 version vectorisée de Pedro ci-dessous est plus rapide:

$ time ./spectral_norm8 5500
1.274224153

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

Enfin, comme le rapporte laxxy ci-dessous pour les compilateurs Intel, il ne semble pas y avoir de grande différence et même le code Fortran le plus simple (spectral_norm1) est parmi les plus rapides.

Tout d'abord, merci d'avoir posté cette question / défi! Comme avertissement, je suis un programmeur C natif avec une certaine expérience de Fortran, et je me sens plus à l'aise en C, donc en tant que tel, je me concentrerai uniquement sur l'amélioration de la version C. J'invite tous les hacks de Fortran à essayer aussi!

Juste pour rappeler aux nouveaux arrivants de quoi il s'agit: la prémisse de base dans ce fil était que gcc / fortran et icc / ifort devraient, puisqu'ils ont respectivement les mêmes back-ends, produire un code équivalent pour le même programme (sémantiquement identique), indépendamment d'être en C ou Fortran. La qualité du résultat ne dépend que de la qualité des implémentations respectives.

J'ai joué un peu avec le code et sur mon ordinateur (ThinkPad 201x, Intel Core i5 M560, 2,67 GHz), en utilisant gcc4.6.1 et les drapeaux de compilation suivants:

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

Je suis aussi allé de l' avant et écrit un SIMD vectorisé version langage C du code 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;
}

Les trois versions ont été compilées avec les mêmes drapeaux et la même gccversion. Notez que j'ai enveloppé l'appel de fonction principale dans une boucle de 0 à 9 pour obtenir des timings plus précis.

$ 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

Ainsi, avec de «meilleurs» indicateurs de compilateur, la version C ++ surpasse la version Fortran et les boucles vectorisées codées à la main ne fournissent qu'une amélioration marginale. Un rapide coup d'œil à l'assembleur de la version C ++ montre que les boucles principales ont également été vectorisées, bien que déroulées de manière plus agressive.

J'ai également regardé l'assembleur généré par gfortranet voici la grande surprise: pas de vectorisation. J'attribue le fait qu'il n'est que légèrement plus lent au problème de la limitation de la bande passante, du moins sur mon architecture. Pour chacune des multiplications matricielles, 230 Mo de données sont parcourus, ce qui submerge à peu près tous les niveaux de cache. Si vous utilisez une valeur d'entrée plus petite, par exemple100 , les différences de performances augmentent considérablement.

En guise de remarque, au lieu d'obséder par la vectorisation, l'alignement et les drapeaux du compilateur, l'optimisation la plus évidente serait de calculer les premières itérations en arithmétique simple précision, jusqu'à ce que nous ayons ~ 8 chiffres du résultat. Les instructions en simple précision sont non seulement plus rapides, mais la quantité de mémoire à déplacer est également divisée par deux.

La réponse de user389 a été supprimée mais permettez-moi de dire que je suis fermement dans son camp: je ne vois pas ce que nous apprenons en comparant les micro-benchmarks dans différentes langues. Cela ne me surprend pas autant que C et Fortran obtiennent à peu près les mêmes performances sur cette référence, étant donné sa courte durée. Mais la référence est également ennuyeuse car elle peut facilement être écrite dans les deux langues en quelques dizaines de lignes. D'un point de vue logiciel, ce n'est pas un cas représentatif: nous devons nous soucier des logiciels qui ont 10 000 ou 100 000 lignes de code et comment les compilateurs le font. Bien sûr, à cette échelle, on découvrira rapidement d'autres choses: que la langue A nécessite 10 000 lignes tandis que la langue B en requiert 50 000. Ou l'inverse, selon ce que vous voulez faire. Et soudain ça '

En d'autres termes, peu m'importe que mon application puisse être 50% plus rapide si je la développais dans Fortran 77 si au lieu de cela il ne me faudrait qu'un mois pour qu'elle fonctionne correctement alors que cela me prendrait 3 mois en F77. Le problème avec la question ici est qu'elle se concentre sur un aspect (noyaux individuels) qui n'est pas pertinent dans la pratique à mon avis.

Il s'avère que je peux écrire un code Python (en utilisant numpy pour effectuer les opérations BLAS) plus rapidement que le code Fortran compilé avec le compilateur gfortran de mon système.

$ 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"

et sn6a.f90, un spectral_norm6.f90 très légèrement modifié:

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

Vérifié cela avec des compilateurs Intel. Avec 11.1 (-fast, impliquant -O3), et avec 12.0 (-O2) les plus rapides sont 1,2,6,7 et 8 (c'est-à-dire les codes Fortran et C "les plus simples", et le C vectorisé à la main) - ceux-ci sont indiscernables les uns des autres à ~ 1,5 s. Les tests 3 et 5 (avec tableau comme fonction) sont plus lents; # 4 Je n'ai pas pu compiler.

Tout particulièrement, si vous compilez avec 12.0 et -O3, plutôt que -O2, les 2 premiers codes Fortran ("les plus simples") ralentissent BEAUCOUP (1.5 -> 10.2 sec.) - ce n'est pas la première fois que je vois quelque chose comme cela, mais cela peut être l'exemple le plus dramatique. Si c'est toujours le cas dans la version actuelle, je pense que ce serait une bonne idée de le signaler à Intel, car il y a clairement quelque chose qui va très mal avec leurs optimisations dans ce cas plutôt simple.

Sinon, je suis d'accord avec Jonathan que ce n'est pas un exercice particulièrement instructif :)