C é mais lento que Fortran no combate às normas espectrais (usando gcc, intel e outros compiladores)?

A conclusão aqui:

Quão melhores são os compiladores Fortran realmente?

é que o gfortran e o gcc são tão rápidos quanto o código simples. Então, eu queria tentar algo mais complicado. Peguei o exemplo do tiroteio da norma espectral. Primeiro eu pré-calculo a matriz 2D A (:, :) e depois calculo a norma. (Esta solução não é permitida na disputa, eu acho.) Eu implementei o Fortran e a versão C. Aqui está o código:

https://github.com/certik/spectral_norm

As versões mais rápidas do gfortran são spectral_norm2.f90 e spectral_norm6.f90 (uma usa o matmul e o dot_product integrados do Fortran, a outra implementa essas duas funções no código - sem diferença de velocidade). O código C / C ++ mais rápido que eu pude escrever é spectral_norm7.cpp. Os tempos da versão git 457d9d9 no meu laptop são:

$ 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

Portanto, a versão do gfortran é um pouco mais rápida. Por que é que? Se você enviar uma solicitação pull com uma implementação C mais rápida (ou apenas colar um código), atualizarei o repositório.

No Fortran, passo um array 2D, enquanto no CI use um array 1D. Sinta-se à vontade para usar uma matriz 2D ou qualquer outra maneira que achar melhor.

Quanto aos compiladores, vamos comparar gcc vs gfortran, icc vs ifort e assim por diante. (Diferentemente da página de tiroteios, que compara ifort x gcc.)

Atualização : usando a versão 179dae2, que melhora o matmul3 () na minha versão C, eles agora são mais rápidos:

$ 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

A versão vetorizada de Pedro abaixo é mais rápida:

$ time ./spectral_norm8 5500
1.274224153

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

Finalmente, como os relatórios laxxy abaixo para os compiladores Intel, não parece haver uma grande diferença, e mesmo o código Fortran mais simples (spectral_norm1) está entre os mais rápidos.

Antes de tudo, obrigado por postar esta pergunta / desafio! Como isenção de responsabilidade, sou programador nativo em C com alguma experiência em Fortran e me sinto mais à vontade em C; portanto, vou me concentrar apenas em melhorar a versão em C. Convido todos os hacks do Fortran a fazerem o mesmo!

Apenas para lembrar aos recém-chegados sobre o que é isso: A premissa básica neste segmento era que gcc / fortran e icc / ifort deveriam, uma vez que possuem os mesmos back-ends, respectivamente, produzir código equivalente para o mesmo programa (semanticamente idêntico), independentemente sendo em C ou Fortran. A qualidade do resultado depende apenas da qualidade das respectivas implementações.

Eu brinquei um pouco com o código e no meu computador (ThinkPad 201x, Intel Core i5 M560, 2,67 GHz), usando gcc4.6.1 e os seguintes sinalizadores do compilador:

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

Também fui adiante e escrevi uma versão em linguagem C do código C ++ vetorizada pelo SIMD 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;
}

Todas as três versões foram compiladas com os mesmos sinalizadores e os mesmos gcc versão. Observe que envolvi a chamada de função principal em um loop de 0 a 9 para obter tempos mais precisos.

$ 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

Portanto, com sinalizadores "melhores" do compilador, a versão C ++ supera a versão Fortran e os loops vetorizados codificados à mão fornecem apenas uma melhoria marginal. Uma rápida olhada no assembler para a versão C ++ mostra que os loops principais também foram vetorizados, embora desenrolados de forma mais agressiva.

Também dei uma olhada no assembler gerado pelo gfortran e aqui está a grande surpresa: sem vetorização. Atribuo o fato de que é apenas um pouco mais lento ao problema de a largura de banda ser limitada, pelo menos na minha arquitetura. Para cada uma das multiplicações de matriz, são percorridos 230 MB de dados, o que praticamente inverte todos os níveis de cache. Se você usar um valor de entrada menor, por exemplo 100, as diferenças de desempenho aumentam consideravelmente.

Como observação, em vez de ficar obcecado com sinalizações de vetorização, alinhamento e compilador, a otimização mais óbvia seria calcular as primeiras iterações na aritmética de precisão única, até obtermos ~ 8 dígitos do resultado. As instruções de precisão única não são apenas mais rápidas, mas a quantidade de memória que precisa ser movida também é reduzida pela metade.

A resposta do user389 foi excluída, mas deixe-me afirmar que estou firmemente no campo dele: não consigo ver o que aprendemos comparando micro-benchmarks em diferentes idiomas. Não me surpreende que C e Fortran obtenham praticamente o mesmo desempenho neste benchmark, dada a sua baixa. Mas o benchmark também é chato, pois pode ser facilmente escrito nos dois idiomas em algumas dezenas de linhas. Do ponto de vista do software, esse não é um caso representativo: devemos nos preocupar com software que tenha 10.000 ou 100.000 linhas de código e como os compiladores fazem isso. Certamente, nessa escala, descobriremos rapidamente outras coisas: o idioma A requer 10.000 linhas, enquanto o idioma B requer 50.000. Ou o contrário, dependendo do que você deseja fazer. E de repente é '

Em outras palavras, não importa muito para mim que talvez meu aplicativo pudesse ser 50% mais rápido se eu o desenvolvesse no Fortran 77. Em vez disso, levará apenas um mês para que ele funcione corretamente, enquanto isso levaria três meses. em F77. O problema com a pergunta aqui é que ela se concentra em um aspecto (núcleos individuais) que não é relevante na prática, na minha opinião.

Acontece que eu posso escrever um código Python (usando numpy para fazer as operações BLAS) mais rápido que o código Fortran compilado com o compilador gfortran do meu sistema.

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

e sn6a.f90, um spectral_norm6.f90 muito levemente modificado:

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

Verifiquei isso com os compiladores Intel. Com 11.1 (-fast, implicando -O3) e com 12.0 (-O2), os mais rápidos são 1,2,6,7 e 8 (ou seja, os códigos "mais simples" de Fortran e C e o C vetorizado à mão) - estes são indistinguíveis um do outro a ~ 1.5s. Os testes 3 e 5 (com matriz em função) são mais lentos; # 4 Não pude compilar.

Notavelmente, se compilar com 12.0 e -O3, em vez de -O2, os 2 primeiros códigos (mais simples) do Fortran desaceleram MUITO (1,5 -> 10,2 segundos) - essa não é a primeira vez que vejo algo como isso, mas esse pode ser o exemplo mais dramático. Se esse ainda for o caso no release atual, acho que seria uma boa ideia denunciá-lo à Intel, pois há claramente algo de errado com suas otimizações nesse caso bastante simples.

Caso contrário, concordo com Jonathan que este não é um exercício particularmente informativo :)