Czy C jest wolniejszy niż Fortran podczas strzelania do normy widmowej (przy użyciu gcc, intel i innych kompilatorów)?

Wniosek tutaj:

O ile lepsze są naprawdę kompilatory Fortran?

jest to, że gfortran i gcc są tak szybkie, jak prosty kod. Chciałem więc spróbować czegoś bardziej skomplikowanego. Wziąłem przykład strzelania do normy spektralnej. Najpierw wstępnie obliczam macierz 2D A (:, :), a następnie obliczam normę. (Myślę, że to rozwiązanie nie jest dozwolone podczas strzelaniny). Zaimplementowałem wersję Fortran i C. Oto kod:

https://github.com/certik/spectral_norm

Najszybsze wersje gfortran to spectral_norm2.f90 i spectral_norm6.f90 (jedna wykorzystuje wbudowane matmul i dot_product Fortrana, druga implementuje te dwie funkcje w kodzie - bez różnicy prędkości). Najszybszy kod C / C ++, jaki udało mi się napisać, to spectral_norm7.cpp. Czasy od wersji git 457d9d9 na moim laptopie wynoszą:

$ 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

Więc wersja gfortran jest trochę szybsza. Dlaczego? Jeśli wyślesz żądanie ściągnięcia z szybszą implementacją C (lub po prostu wkleisz kod), zaktualizuję repozytorium.

W Fortranie przekazuję tablicę 2D, natomiast w CI używam tablicy 1D. Możesz swobodnie korzystać z tablicy 2D lub w dowolny inny sposób, który uznasz za odpowiedni.

Jeśli chodzi o kompilatory, porównajmy gcc vs gfortran, icc vs ifort i tak dalej. (W przeciwieństwie do strony z rzutami, która porównuje ifort z gcc.)

Aktualizacja : używając wersji 179dae2, która poprawia matmul3 () w mojej wersji C, są teraz tak szybkie:

$ 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

Wersja wektoryzacyjna Pedro poniżej jest szybsza:

$ time ./spectral_norm8 5500
1.274224153

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

Wreszcie, jak podają poniższe raporty Laxxy dla kompilatorów Intela, wydaje się, że nie ma tu dużej różnicy, a nawet najprostszy kod Fortrana (spectral_norm1) jest jednym z najszybszych.

Przede wszystkim dziękuję za opublikowanie tego pytania / wyzwania! Jako zrzeczenie się, jestem rodzimym programistą C z pewnym doświadczeniem w Fortranie i czuję się jak w domu w C, dlatego skupię się tylko na ulepszeniu wersji C. Zapraszam też wszystkich hacków z Fortranu!

Aby przypomnieć początkującym o tym, o co w tym chodzi: Podstawowym założeniem w tym wątku było to, że gcc / fortran i icc / ifort powinny, ponieważ mają one odpowiednio te same zaplecza, produkować równoważny kod dla tego samego (semantycznie identycznego) programu, niezależnie od tego jest w C lub Fortran. Jakość wyniku zależy tylko od jakości odpowiednich wdrożeń.

Grałem trochę z kodem i na moim komputerze (ThinkPad 201x, Intel Core i5 M560, 2,67 GHz), używając wersji gcc4.6.1 i następujących flag kompilatora:

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

Napisałem też wektoryzowaną przez SIMD wersję C języka 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;
}

Wszystkie trzy wersje zostały skompilowane z tymi samymi flagami i tą samą gccwersją. Zauważ, że zawinąłem wywołanie funkcji głównej w pętli od 0..9, aby uzyskać dokładniejsze czasy.

$ 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

Dzięki „lepszym” flagom kompilatora wersja C ++ przewyższa wersję Fortrana, a ręcznie kodowane wektoryzowane pętle zapewniają jedynie nieznaczną poprawę. Szybkie spojrzenie na asembler dla wersji C ++ pokazuje, że główne pętle również zostały wektoryzowane, aczkolwiek rozwinięte bardziej agresywnie.

Spojrzałem również na asembler wygenerowany przez gfortrani oto wielka niespodzianka: brak wektoryzacji. Przypisuję fakt, że jest on tylko nieznacznie wolniejszy, ponieważ problem ma ograniczoną przepustowość, przynajmniej w mojej architekturze. Dla każdego z mnożeń macierzy przemierzanych jest 230 MB danych, co prawie zapełnia wszystkie poziomy pamięci podręcznej. Jeśli użyjesz mniejszej wartości wejściowej, np. 100Różnice wydajności znacznie wzrosną.

Na marginesie, zamiast obsesji na punkcie wektoryzacji, wyrównywania i flag kompilatora, najbardziej oczywistą optymalizacją byłoby obliczenie pierwszych kilku iteracji w arytmetyce pojedynczej precyzji, dopóki nie uzyskamy ~ 8 cyfr wyniku. Instrukcje pojedynczej precyzji są nie tylko szybsze, ale ilość pamięci, którą należy przenieść, jest również zmniejszona o połowę.

odpowiedź użytkownika389 została usunięta, ale pozwól mi stwierdzić, że jestem mocno w jego obozie: nie widzę tego, czego się uczymy, porównując mikro-testy w różnych językach. Nie jest dla mnie zaskoczeniem, że C i Fortran uzyskują prawie taką samą wydajność na tym teście, biorąc pod uwagę jego krótki czas. Benchmark jest również nudny, ponieważ można go łatwo napisać w obu językach w kilkudziesięciu wierszach. Z punktu widzenia oprogramowania nie jest to reprezentatywny przypadek: powinniśmy dbać o oprogramowanie, które ma 10 000 lub 100 000 wierszy kodu i jak kompilują to. Oczywiście w tej skali można szybko dowiedzieć się innych rzeczy: ten język A wymaga 10 000 wierszy, a język B 50 000. Lub na odwrót, w zależności od tego, co chcesz zrobić. I nagle to

Innymi słowy, nie ma dla mnie większego znaczenia, że ​​może moja aplikacja mogłaby być o 50% szybsza, jeśli opracowałbym ją w Fortran 77, jeśli zamiast tego zajmie mi to tylko 1 miesiąc, aby uruchomić ją poprawnie, podczas gdy zajmie mi to 3 miesiące w F77. Problem w tym pytaniu polega na tym, że koncentruje się on na aspekcie (poszczególnych jądrach), który moim zdaniem nie jest istotny w praktyce.

Okazuje się, że mogę napisać kod Pythona (używając Numpy do wykonywania operacji BLAS) szybciej niż kod Fortran skompilowany z kompilatorem gfortran mojego systemu.

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

i sn6a.f90, bardzo lekko zmodyfikowany spectral_norm6.f90:

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

Sprawdzono to za pomocą kompilatorów Intel. Przy 11.1 (-szybkim, co sugeruje -O3), a przy 12.0 (-O2) najszybsze to 1,2,6,7 i 8 (tj. „Najprostsze” kody Fortran i C oraz ręcznie wektoryzowane C) - są nie do odróżnienia od siebie po ~ 1,5 s. Testy 3 i 5 (z tablicą jako funkcją) przebiegają wolniej; # 4 Nie mogłem skompilować.

W szczególności, jeśli kompilujemy z 12.0 i -O3, a nie -O2, pierwsze 2 („najprostsze”) kody Fortran spowalniają DUŻO (1,5 -> 10,2 sek.) - to nie pierwszy raz, kiedy widzę coś takiego to, ale może to być najbardziej dramatyczny przykład. Jeśli nadal tak jest w bieżącej wersji, myślę, że dobrym pomysłem byłoby zgłoszenie tego do Intela, ponieważ najwyraźniej coś jest nie tak z ich optymalizacjami w tej dość prostej sprawie.

W przeciwnym razie zgadzam się z Jonathanem, że nie jest to szczególnie pouczające ćwiczenie :)