首先，感谢您发布这个问题/挑战！作为免责声明，我是一名具有一定 Fortran 经验的本地 C 程序员，并且对 C 感到最熟悉，因此，我将只专注于改进 C 版本。我也邀请所有 Fortran 黑客一起去！

只是为了提醒新手这是关于什么的：这个线程的基本前提是 gcc/fortran 和 icc/ifort 应该，因为它们分别具有相同的后端，为相同的（语义相同的）程序生成相同的代码，不管它在 C 或 Fortran 中。结果的质量仅取决于各自实现的质量。

我在我的电脑（ThinkPad 201x，Intel Core i5 M560，2.67 GHz）上玩了一下代码，使用gcc了 4.6.1 和以下编译器标志：

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

我还继续编写了 C++ 代码的 SIMD 向量化 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;
}

所有三个版本都使用相同的标志和相同的gcc版本编译。请注意，我将主函数调用包装在从 0..9 开始的循环中，以获得更准确的计时。

$ 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

因此，使用“更好”的编译器标志，C++ 版本的性能优于 Fortran 版本，而手动编码的矢量化循环仅提供边际改进。快速浏览一下 C++ 版本的汇编器会发现，主循环也已被矢量化，尽管展开得更积极。

我还查看了由 生成的汇编程序gfortran，这是一个很大的惊喜：没有矢量化。我将它只是稍微慢一点的事实归因于带宽有限的问题，至少在我的架构上是这样。对于每个矩阵乘法，要遍历 230MB 的数据，这几乎淹没了所有级别的缓存。如果您使用较小的输入值，例如100，则性能差异会显着增加。

作为旁注，最明显的优化是在单精度算术中计算前几次迭代，而不是痴迷于矢量化、对齐和编译器标志，直到我们得到大约 8 位结果。单精度指令不仅速度更快，而且必须移动的内存量也减半。

user389 的回答已被删除，但让我声明我坚定地支持他的阵营：我看不到通过比较不同语言的微基准测试我们学到了什么。考虑到 C 和 Fortran 在这个基准测试中的性能非常短，我并不感到惊讶。但是基准测试也很无聊，因为它可以很容易地用两种语言在几十行中编写。从软件的角度来看，这不是一个具有代表性的案例：我们应该关心具有 10,000 或 100,000 行代码的软件以及编译器如何处理这些代码。当然，在这种规模下，人们很快就会发现其他事情：语言 A 需要 10,000 行，而语言 B 需要 50,000 行。或者反过来，取决于你想做什么。突然间'

换句话说，如果我在 Fortran 77 中开发我的应用程序可能会快 50%，这对我来说并不重要，如果我只需要 1 个月就可以让它正确运行，而我需要 3 个月在 F77 中。这里的问题的问题在于，它侧重于我认为在实践中不相关的方面（单个内核）。

事实证明，我可以比使用系统的 gfortran 编译器编译的 Fortran 代码更快地编写 Python 代码（使用 numpy 执行 BLAS 操作）。

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

和 sn6a.f90，一个非常轻微修改的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

用英特尔编译器检查了这一点。对于 11.1（-fast，意味着 -O3）和 12.0 (-O2)，最快的是 1、2、6、7 和 8（即“最简单的”Fortran 和 C 代码，以及手动矢量化的 C） - 这些在约 1.5 秒时彼此无法区分。测试 3 和 5（将数组作为函数）较慢；#4 我无法编译。

值得注意的是，如果使用 12.0 和 -O3 而不是 -O2 进行编译，前 2 个（“最简单”）Fortran 代码会减慢很多（1.5 -> 10.2 秒）——这不是我第一次看到类似的东西这个，但这可能是最戏剧性的例子。如果在当前版本中仍然是这种情况，我认为向英特尔报告是个好主意，因为在这个相当简单的情况下，他们的优化显然出了问题。

否则我同意乔纳森的观点，这不是一个特别有用的练习:)