我用于估算英特尔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。
矩阵乘法的FLOPs数量大约为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缓存中，因此缓存似乎不是限制因素（虽然也许我可以更有效地使用它）。
编辑2： Eigen基准测试显示其与MKL一样好（对于SSE）。 它们使用Intel（R）Core（TM）2 Quad CPU Q9400 @ 2.66GHz。 因此，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
编辑3: 这是我的代码，供有兴趣的人参考。我可以获得比这更好一点的结果，但不会好多少。我正在使用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]);
}