矩阵乘法Intel/AMD CPU的最大FLOPS

Question

我估计Intel CPU的最大FLOPs/s的公式是

Max SP FLOPs/s = frequencey * 4 SSE(8AVX) * 2 (MAC) * number of cores (not HW threads)
Max DP FLOPs/s = 0.5 * Max SP FLOPs/s

所谓MAC，我的意思是CPU可以同时进行一次SSE(AVX)乘法和加法运算。在我使用的系统上，负载下的最大频率是2.66 GHz。它只有SSE，核心数(非硬件线程)为4。这意味着:最大SP FLOPs/s = 85.12 GFLOPs/s。

矩阵乘法的FLOP数约为2*n*m*k。对于一个带有n=1000的方阵，它是2*10E9个FLOPs (20亿个FLOPs)。一旦我知道了时间，我就可以估计FLOPs/s。

然而，对于我自己的代码，我能得到的最好的结果是大约40 SP GFLOPs/s，例如使用n=1000。我和Eigen得到了大致相同的结果。这大约是45%的效率。我对最大值的计算是错误的吗？英特尔CPU处理大型密集矩阵乘法的最佳效率是多少？有没有人有一篇论文来描述这一点？

我知道在GPU上的效率可以超过60%。http://www.anandtech.com/show/6774/nvidias-geforce-gtx-titan-part-2-titans-performance-unveiled/3

编辑:我得到了与n=500类似的结果，它可以很容易地放入我系统的12MB L3缓存中，所以缓存似乎不是限制因素(尽管我可能可以更有效地使用它)。

Edit2:特征基准测试表明它与MKL (用于上交所)一样好。它们使用英特尔(R)酷睿(TM)2四核中央处理器Q9400 @2.66 They。因此，2.66 * 2 (DP SSE) *2 MAC *4核心= 42.25 DP GFLOPs/s。您可以在图表上看到它们都小于20。像我这样的大约45%。http://eigen.tuxfamily.org/index.php?title=Benchmark

http://ark.intel.com/products/35365/Intel-Core2-Quad-Processor-Q9400-6M-Cache-2_66-GHz-1333-MHz-FSB

Edit3:这是我为关心的人编写的代码。我可以得到比这个稍微好一点的结果，但也不会好到哪里去。我使用Agner Fog的vectorclass来表示SEE/AVX。并将Vec8f设置为float8，将Vec4d设置为double4

//SGEMM and AVX call MM_tile<float, float8>(nthreads, a, b, c, n, m, k);
template <typename ftype, typename floatn>
void GEMM_tile(const int nthreads, const ftype*A , const ftype* B, ftype* C, const int N, const int M, const int K) {       
    for(int i=0; i<N; i++) {
       for(int j=0; j<K; j++) {
           C[K*i + j] = 0;
       }
    }   
    const int nc = 32;
    const int kc = 32;
    const int mc = 32;
    omp_set_num_threads(nthreads);
    #pragma omp parallel for if(nthreads>1)
    for(int ii=0; ii<N; ii+=nc) {
        for(int jj=0; jj<K; jj+=kc)
            for(int ll=0; ll<M; ll+=mc) {
                const int nb = min(N-ii, nc);
                const int kb = min(K-jj, kc);
                const int mb = min(M-ll, mc);
                MM_block<ftype, floatn>(nb, mb, kb, &A[M*ii+ll], N, &B[K*ll+jj], K, &C[K*ii+jj], K );
            }
     }
}

template <typename ftype, typename floatn>
void MM_block(int n, int m, int k, const ftype *a, const int stridea, 
                                   const ftype *b, const int strideb,
                                   ftype *c, const int stridec ) {
    const int vec_size = sizeof(floatn)/sizeof(ftype);
    for(int i=0; i<n; i+=4) {
        for(int j=0; j<k; j+=vec_size) {
            Dot4x4_vec_block<ftype, floatn>(m, &a[strideb*i], &b[j], &c[stridec*i + j], stridea, strideb, stridec);
    }
}

template <typename ftype, typename floatn>
inline void Dot4x4_vec_block(const int n, const ftype *a, const ftype *b, ftype *c, const int stridea, const int strideb, const int stridec) {
    floatn tmp0, tmp1, tmp2, tmp3;
    load(tmp0, &c[stridec*0]);
    load(tmp1, &c[stridec*1]);
    load(tmp2, &c[stridec*2]);
    load(tmp3, &c[stridec*3]);

    ftype *a0_ptr = (ftype*)&a[stridea*0];
    ftype *a1_ptr = (ftype*)&a[stridea*1];
    ftype *a2_ptr = (ftype*)&a[stridea*2];
    ftype *a3_ptr = (ftype*)&a[stridea*3];
    for(int i=0; i<n; i++) {
        floatn breg = floatn().load(&b[i*strideb + 0]);

        floatn areg0 = *a0_ptr++;
        floatn areg1 = *a1_ptr++;
        floatn areg2 = *a2_ptr++;
        floatn areg3 = *a3_ptr++;

        tmp0 += areg0 * breg;
        tmp1 += areg1 * breg;
        tmp2 += areg2 * breg;
        tmp3 += areg3 * breg;
}
    tmp0.store(&c[stridec*0]);
    tmp1.store(&c[stridec*1]);
    tmp2.store(&c[stridec*2]);
    tmp3.store(&c[stridec*3]);
}

Answer 1

通常，处理吞吐量的限制因素是内存带宽，特别是在您的工作集不能放入CPU缓存的情况下(您的float的1000x1000矩阵将占用大约4MB，而您的CPU可能有2MB的L3缓存)。在这种情况下，算法的结构可能会对它的执行方式产生很大的影响，但您通常会在某些时候遇到困难，因为您正在等待来自内存层次结构中更高级别的值。

此外，您的理论数字假设您有足够的指令，而不存在数据依赖性，以便在每个周期内保持所有执行单元的任务。在实践中，这可能很难做到。我不确定一般矩阵乘法的最佳吞吐量是多少，但请查看this previous question，了解如何最大限度地提高指令吞吐量。