13
$\begingroup$

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.

$\endgroup$
  • 5
    $\begingroup$ I'm not anywhere near a compiler right now, but consider adding the restrict keyword to your arrays. Aliasing of pointers is usually the difference between Fortran and C function calls on arrays. Also, Fortran stores memory in column-major order and C in row-major. $\endgroup$ – moyner May 13 '12 at 8:46
  • 1
    $\begingroup$ -1 Body of this question talks about implementations, but the title asks which language is faster? How can a language have an attribute of speed? You should edit the question title so it reflects the body of the question. $\endgroup$ – milancurcic May 13 '12 at 17:44
  • $\begingroup$ @IRO-bot, I fixed it. Let me know if it looks ok to you know. $\endgroup$ – Ondřej Čertík May 13 '12 at 17:58
  • 1
    $\begingroup$ Actually the conclusions on "How much better are Fortran compilers really?" are not quite correct in that thread. I had tried the benchmark on a Cray with GCC, PGI, CRAY and Intel compilers and with 3 compilers Fortran was faster than C (b/w 5-40%). Cray compilers produced the fastest Fortran/C code but the Fortran code was 40% faster. I will post detailed results when I get time. Btw anyone with access to Cray machines can verify the benchmark. It is a good platform because 4-5 compilers are available and relevant flags are engaged automatically by the ftn/cc wrapper. $\endgroup$ – stali May 15 '12 at 17:17
  • $\begingroup$ also checked with pgf95/pgcc (11.10) on an Opteron system: #1 and #2 are the fastest (faster than ifort by ~20%), then #6, #8, #7 (in that order). pgf95 was faster than ifort for all your fortran codes, and icpc was faster than pgcpp for all C -- I should mention that for my stuff, I usually find ifort faster, even on the same AMD system. $\endgroup$ – laxxy May 18 '12 at 3:03
12
$\begingroup$

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.

$\endgroup$
  • $\begingroup$ Thanks so much for your time! I was hoping you would reply. :) So first I updated the Makefile to use your flags. Then I put your C code as spectral_norm8.c and updated README. I updated the timings on my machine (github.com/certik/spectral_norm/wiki/Timings) and as you can see, compiler flags did not make the C version faster on my machine (i.e. gfortran still wins), but your SIMD-vectorized version beats gfortran. $\endgroup$ – Ondřej Čertík May 13 '12 at 17:54
  • $\begingroup$ @OndřejČertík: Just out of curiosity, what version of gcc/gfortran are you using? In the previous threads, different versions gave significantly different results. $\endgroup$ – Pedro May 13 '12 at 17:59
  • $\begingroup$ I use 4.6.1-9ubuntu3. Do you have access to intel compilers? My experience with gfortran is that sometimes it doesn't (yet) produce optimal code. IFort usually does. $\endgroup$ – Ondřej Čertík May 13 '12 at 18:03
  • 1
    $\begingroup$ @OndřejČertík: Now the results make more sense! I had overlooked that matmul2 in the Fortran version is semantically equivalent to matmul3 in my C version. The two versions are really now the same and thus gcc/gfortran should produce the same results for both, e.g. no one front-end/language is better than the other in this case. gcc just has the advantage that we could exploit vectorized instructions should we choose to. $\endgroup$ – Pedro May 15 '12 at 16:38
  • 1
    $\begingroup$ @cjordan1: I chose to use the vector_size attribute in order to make the code platform-independent, i.e. using this syntax, gcc should be able to generate vectorized code for other platforms e.g. using AltiVec on the IBM Power architecture. $\endgroup$ – Pedro May 20 '13 at 18:54
7
$\begingroup$

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.

$\endgroup$
  • $\begingroup$ Agreed. For what it's worth, aside from very, very minor edits (-3 characters, +9 characters), I agreed with the main sentiment of his answer. As far as I know, the C++/C/Fortran compiler debate matters only when one has exhausted every other possible avenue for performance enhancement, which is why, for 99.9% of people, these comparisons don't matter. I don't find the discussion particularly enlightening, but I know of at least one person on the site who can attest to choosing Fortran over C and C++ for performance reasons, which is why I can't say it's totally useless. $\endgroup$ – Geoff Oxberry May 17 '12 at 9:25
  • 4
    $\begingroup$ I agree with your main point, but I still think that this discussion is useful as there are a number of people out there who still somehow believe there is some magic that makes one language "faster" than the other, despite the use of identical compiler backends. I contribute to these debates mainly to try to dispel this myth. As for the methodology, there is no "representative case", and in my opinion, taking something as simple as matrix-vector multiplies is a good thing, as it gives the compilers enough space to show what they can do or not. $\endgroup$ – Pedro May 17 '12 at 9:58
  • $\begingroup$ @GeoffOxberry: Sure, you'll always find people who use one language rather than another for more or less well articulated and reasoned causes. My question, though, would be how fast Fortran would be if one were to use the data structures that appear in, say, an unstructured, adaptive finite element meshes. Aside from the fact that this would be awkward to implement in Fortran (everyone who implements this in C++ uses the STL heavily throughout), would Fortran really be faster for this kind of code that has no tight loops, many indirections, lots of ifs? $\endgroup$ – Wolfgang Bangerth May 17 '12 at 13:27
  • $\begingroup$ @WolfgangBangerth: Like I said in my first comment, I agree with you and with user389 (Jonathan Dursi), so asking me that question is pointless. That said, I would invite anyone who does believe that choice of language (among C++/C/Fortran) is important for performance in their application to answer your question. Sadly, I suspect this sort of debate can be had for compiler versions. $\endgroup$ – Geoff Oxberry May 17 '12 at 19:06
  • $\begingroup$ @GeoffOxberry: Yes, and I obviously didn't mean that you needed to answer that question. $\endgroup$ – Wolfgang Bangerth May 17 '12 at 19:22
5
$\begingroup$

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
$\endgroup$
  • 1
    $\begingroup$ Tongue in cheek, I presume? $\endgroup$ – Robert Harvey May 15 '12 at 21:11
  • $\begingroup$ -1 for not answering the question, but I think you know that already. $\endgroup$ – Pedro May 15 '12 at 21:11
  • $\begingroup$ Interesting, what version of gfortran did you use, and did you test out the C code available in the repository with Pedro's flags? $\endgroup$ – Aron Ahmadia May 15 '12 at 21:19
  • 1
    $\begingroup$ Actually, I think it's clearer now, assuming you weren't being sarcastic. $\endgroup$ – Robert Harvey May 15 '12 at 21:23
  • 1
    $\begingroup$ Since this post, and none of the other questions or posts, are being edited by Aron in such a way to better match his opinions, even though my whole point is that all the posts should be labelled with exactly such "these results are meaningless" caveats, I'm just deleting it. $\endgroup$ – user389 May 15 '12 at 21:27
3
$\begingroup$

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 :)

$\endgroup$
  • $\begingroup$ Thanks for checking it! This confirms my experience, that gfortran is not yet fully mature, because for some reason the matmul operation is slow. So the conclusion for me is to simply use matmul and keep the Fortran code simple. $\endgroup$ – Ondřej Čertík May 15 '12 at 16:24
  • $\begingroup$ On the other hand, I think gfortran has a command line option to automatically convert all matmul() calls into BLAS calls (maybe also dot_product(), not sure). Never tried it though. $\endgroup$ – laxxy May 15 '12 at 16:44

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.