Is C slower than Fortran on the spectral norm shootout (using gcc, intel and other compilers)?

Question

The conclusion here:

How much better are Fortran compilers really?

is that gfortran and gcc are as fast for simple code. So I wanted try something more complicated. I took the spectral norm shootout example. I first precalculate the 2D matrix A(:, :), and then calculate the norm. (This solution is not allowed on the shootout I think.) I have implemented Fortran and C version. Here is the code:

https://github.com/certik/spectral_norm

The fastest gfortran versions are spectral_norm2.f90 and spectral_norm6.f90 (one uses the Fortran's built-in matmul and dot_product, the other implements these two functions in the code -- with no difference in speed). The fastest C/C++ code that I was able to write is spectral_norm7.cpp. Timings as of the git version 457d9d9 on my laptop are:

$ 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

So gfortran's version is a little faster. Why is that? If you send a pull request with a faster C implementation (or just paste a code), I'll update the repository.

In Fortran I pass a 2D array around, while in C I use a 1D array. Feel free to use a 2D array or any other way that you see fit.

As to compilers, let's compare gcc vs gfortran, icc vs ifort, and so on. (Unlike the shootout page, which compares ifort vs gcc.)

Update: using the version 179dae2, which improves matmul3() in my C version, they are now as fast:

$ 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

Pedro's vectorized version below is faster:

$ time ./spectral_norm8 5500
1.274224153

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

Finally, as laxxy reports below for Intel compilers, there doesn't seem to be a big difference there and even the simplest Fortran code (spectral_norm1) is among the fastest.

Answer 1

First of all, thanks for posting this question/challenge! As a disclaimer, I'm a native C programmer with some Fortran experience, and feel most at home in C, so as such, I will focus only on improving the C version. I invite all Fortran hacks to have their go too!

Just to remind newcomers about what this is about: The basic premise in this thread was that gcc/fortran and icc/ifort should, since they have the same back-ends respectively, produce equivalent code for the same (semantically identical) program, irrespective of it being in C or Fortran. The quality of the result depends only on the quality of the respective implementations.

I played around with the code a bit, and on my computer (ThinkPad 201x, Intel Core i5 M560, 2.67 GHz), using gcc 4.6.1 and the following compiler flags:

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

I also went ahead and wrote a SIMD-vectorized C-language version of the C++ code, 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;
}

All three versions were compiled with the same flags and the same gcc version. Note that I wrapped the main function call in a loop from 0..9 to get more accurate timings.

$ 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

So with "better" compiler flags, the C++ version out-performs the Fortran version and hand-coded vectorized loops only provide a marginal improvement. A quick look at the assembler for the C++ version shows that the main loops have also been vectorized, albeit unrolled more aggressively.

I also had a look at the assembler generated by gfortran and here's the big surprise: no vectorization. I attribute the fact that it's only marginally slower to the problem being bandwidth limited, at least on my architecture. For each of the matrix multiplications, 230MB of data are traversed, which pretty-much swamps all levels of cache. If you use a smaller input value, e.g. 100, the performance differences grow considerably.

As a side-note, instead of obsessing about vectorization, alignment and compiler flags, the most obvious optimization would be to compute the first few iterations in single-precision arithmetic, until we have ~8 digits of the result. The single-precision instructions are not only faster, but the amount of memory that has to be moved around is also halved.

Answer 2

user389's answer has been deleted but let me state that I'm firmly in his camp: I fail to see what we learn by comparing micro-benchmarks in different languages. It doesn't come as much of a surprise to me that C and Fortran get pretty much the same performance on this benchmark given how short it is. But the benchmark is also boring since it can easily be written in both languages in a couple dozen lines. From a software viewpoint, that is not a representative case: we should care about software that has 10,000 or 100,000 lines of code and how compilers do on that. Of course, at that scale, one will quickly find out other things: that language A requires 10,000 lines while language B requires 50,000. Or the other way around, depending on what you want to do. And suddenly it's not so much about the difference in milliseconds upon execution anymore, but a difference of months of developer time.

In other words, it doesn't matter much to me that maybe my application could be 50% faster if I developed it in Fortran 77 if instead it will only take me 1 month to get it to run correctly while it would take me 3 months in F77. The problem with the question here is that it focuses on an aspect (individual kernels) that isn't relevant in practice in my view.

Answer 3

It turns out that I can write a Python code (using numpy to do the BLAS operations) faster than the Fortran code compiled with my system's gfortran compiler.

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

and sn6a.f90, a very lightly modified 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

Answer 4

Checked this with Intel compilers. With 11.1 (-fast, implying -O3), and with 12.0 (-O2) the fastest ones are 1,2,6,7,and 8 (i.e. the "simplest" Fortran and C codes, and the hand-vectorized C) -- these are indistinguishable from each other at ~1.5s. Tests 3 and 5 (with array as a function) are slower; #4 I couldn't compile.

Quite notably, if compiling with 12.0 and -O3, rather than -O2, the first 2 ("simplest") Fortran codes slow down A LOT (1.5 -> 10.2 sec.) -- this is not the first time I see something like this, but this may be the most dramatic example. If this still is the case in the current release, I think it would be a good idea to report it to Intel, as there is clearly something going very wrong with their optimizations in this rather simple case.

Otherwise I agree with Jonathan that this is not a particularly informative exercise :)