Ist C langsamer als Fortran bei der Spektralnorm-Schießerei (mit gcc, Intel und anderen Compilern)?

Das Fazit hier:

Wie viel besser sind Fortran-Compiler wirklich?

ist, dass gfortran und gcc für einfachen Code genauso schnell sind. Also wollte ich etwas komplizierteres ausprobieren. Ich nahm das Spektralnorm-Shootout-Beispiel. Ich berechne zuerst die 2D-Matrix A (:, :) und dann die Norm. (Ich denke, diese Lösung ist im Shootout nicht zulässig.) Ich habe die Fortran- und C-Version implementiert. Hier ist der Code:

https://github.com/certik/spectral_norm

Die schnellsten Gfortran-Versionen sind spectral_norm2.f90 und spectral_norm6.f90 (eine verwendet das in Fortran integrierte matmul und dot_product, die andere implementiert diese beiden Funktionen im Code - ohne Geschwindigkeitsunterschied). Der schnellste C / C ++ - Code, den ich schreiben konnte, ist spectral_norm7.cpp. Timings ab der Git-Version 457d9d9 auf meinem Laptop sind:

$ 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

Gfortrans Version ist also etwas schneller. Warum das? Wenn Sie eine Pull-Anfrage mit einer schnelleren C-Implementierung senden (oder einfach einen Code einfügen), aktualisiere ich das Repository.

In Fortran lasse ich ein 2D-Array herumlaufen, während ich in CI ein 1D-Array verwende. Fühlen Sie sich frei, ein 2D-Array oder eine andere Weise zu verwenden, die Sie für richtig halten.

In Bezug auf Compiler vergleichen wir gcc mit gfortran, icc mit ifort und so weiter. (Im Gegensatz zur Shootout-Seite, die ifort mit gcc vergleicht.)

Update : Mit der Version 179dae2, die matmul3 () in meiner C-Version verbessert, sind sie jetzt so schnell:

$ 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

Die vektorisierte Version von Pedro unten ist schneller:

$ time ./spectral_norm8 5500
1.274224153

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

Wie Laxxy weiter unten für Intel-Compiler berichtet, scheint es keinen großen Unterschied zu geben, und selbst der einfachste Fortran-Code (spectral_norm1) gehört zu den schnellsten.

Zunächst einmal vielen Dank für das Posten dieser Frage / Herausforderung! Als Haftungsausschluss bin ich ein gebürtiger C-Programmierer mit Fortran-Erfahrung und fühle mich in C am wohlsten. Daher werde ich mich nur auf die Verbesserung der C-Version konzentrieren. Ich lade alle Fortran-Hacks ein, mitzumachen!

Nur um Neulinge daran zu erinnern, worum es geht: Die Grundvoraussetzung in diesem Thread war, dass gcc / fortran und icc / ifort, da sie jeweils die gleichen Backends haben, für dasselbe (semantisch identische) Programm einen entsprechenden Code erzeugen sollten, unabhängig davon davon in C oder Fortran. Die Qualität des Ergebnisses hängt nur von der Qualität der jeweiligen Implementierungen ab.

Ich habe ein wenig mit dem Code herumgespielt und auf meinem Computer (ThinkPad 201x, Intel Core i5 M560, 2,67 GHz) mit gcc4.6.1 und den folgenden Compiler-Flags gearbeitet:

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

Ich habe auch eine SIMD-vektorisierte C-Sprachversion des C ++ - Codes geschrieben 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;
}

Alle drei Versionen wurden mit denselben Flags und derselben gccVersion kompiliert . Beachten Sie, dass ich den Hauptfunktionsaufruf von 0..9 in eine Schleife eingebunden habe, um genauere Timings zu erhalten.

$ 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

Mit "besseren" Compiler-Flags übertrifft die C ++ - Version die Fortran-Version, und handcodierte vektorisierte Schleifen bieten nur eine marginale Verbesserung. Ein kurzer Blick auf den Assembler für die C ++ - Version zeigt, dass die Hauptschleifen ebenfalls vektorisiert wurden, wenn auch aggressiver abgewickelt.

Ich habe mir auch den Assembler von angeschaut gfortranund hier ist die große Überraschung: keine Vektorisierung. Ich schreibe die Tatsache zu, dass es nur unwesentlich langsamer ist, wenn die Bandbreite begrenzt ist, zumindest in meiner Architektur. Für jede der Matrixmultiplikationen werden 230 MB Daten durchlaufen, wodurch praktisch alle Cache-Ebenen überlastet werden. Wenn Sie beispielsweise einen kleineren Eingabewert verwenden 100, nehmen die Leistungsunterschiede erheblich zu.

Anstatt von Vektorisierung, Ausrichtung und Compiler-Flags besessen zu sein, besteht die offensichtlichste Optimierung darin, die ersten paar Iterationen in Arithmetik mit einfacher Genauigkeit zu berechnen, bis wir ~ 8 Ziffern des Ergebnisses haben. Die Befehle mit einfacher Genauigkeit sind nicht nur schneller, sondern die Menge an Speicher, die verschoben werden muss, halbiert sich auch.

Die Antwort von user389 wurde gelöscht, aber lassen Sie mich feststellen, dass ich fest in seinem Lager bin: Ich sehe nicht, was wir lernen, indem ich Mikro-Benchmarks in verschiedenen Sprachen vergleiche. Es ist für mich keine große Überraschung, dass C und Fortran auf diesem Benchmark in Anbetracht der kurzen Zeit ungefähr die gleiche Leistung erbringen. Der Benchmark ist aber auch langweilig, da er problemlos in zwei Sprachen in ein paar Dutzend Zeilen geschrieben werden kann. Aus Sicht der Software ist dies kein repräsentativer Fall: Wir sollten uns mit Software befassen, die 10.000 oder 100.000 Codezeilen enthält, und wie Compiler dies tun. Auf dieser Skala wird man natürlich schnell andere Dinge herausfinden: Für Sprache A sind 10.000 Zeilen erforderlich, für Sprache B sind 50.000. Oder umgekehrt, je nachdem, was Sie tun möchten. Und plötzlich ist es

Mit anderen Worten, es spielt für mich keine Rolle, dass meine Anwendung möglicherweise 50% schneller wäre, wenn ich sie in Fortran 77 entwickeln würde. Stattdessen würde ich nur 1 Monat brauchen, um sie ordnungsgemäß auszuführen, während ich 3 Monate brauchen würde in F77. Das Problem bei dieser Frage ist, dass sie sich auf einen Aspekt (einzelne Kernel) konzentriert, der aus meiner Sicht in der Praxis nicht relevant ist.

Es stellt sich heraus, dass ich einen Python-Code (mit numpy, um die BLAS-Operationen auszuführen) schneller schreiben kann als den Fortran-Code, der mit dem Gfortran-Compiler meines Systems kompiliert wurde.

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

und sn6a.f90, ein sehr leicht modifizierter 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

Dies wurde mit Intel-Compilern überprüft. Mit 11.1 (-fast, impliziert -O3) und mit 12.0 (-O2) sind die schnellsten 1,2,6,7 und 8 (dh die "einfachsten" Fortran- und C-Codes und das handvektorisierte C) - diese sind bei ~ 1,5s nicht voneinander zu unterscheiden. Die Tests 3 und 5 (mit Array als Funktion) sind langsamer. # 4 Ich konnte nicht kompilieren.

Bemerkenswerterweise verlangsamen beim Kompilieren mit 12.0 und -O3 anstelle von -O2 die ersten 2 ("einfachsten") Fortran-Codes VIEL (1,5 -> 10,2 Sek.) - dies ist nicht das erste Mal, dass ich so etwas sehe dies, aber dies kann das dramatischste Beispiel sein. Wenn dies in der aktuellen Version immer noch der Fall ist, ist es meiner Meinung nach eine gute Idee, dies Intel zu melden, da in diesem eher einfachen Fall eindeutig etwas mit den Optimierungen nicht stimmt.

Ansonsten stimme ich Jonathan zu, dass dies keine besonders informative Übung ist :)