如何使用点积获得峰值 CPU 性能？
Posted
技术标签:
【中文标题】如何使用点积获得峰值 CPU 性能？【英文标题】：How Do I Attain Peak CPU Performance With Dot Product? 【发布时间】：2014-07-28 22:05:33 【问题描述】：问题
我一直在研究 HPC，特别是使用矩阵乘法作为我的项目（请参阅我在个人资料中的其他帖子）。我在这些方面取得了不错的成绩，但还不够好。我退后一步，看看我在点积计算方面能做多好。
点积与矩阵乘法
点积更简单，可以让我在不处理打包和其他相关问题的情况下测试 HPC 概念。缓存阻塞仍然是一个问题，这是我的第二个问题。
算法
将n两个double数组A和B中的对应元素相乘并将它们相加。组装中的double 点积只是movapd、mulpd、addpd 的系列。展开并以一种巧妙的方式排列，可以让movapd/mulpd/addpd 组在不同的xmm 寄存器上运行，因此是独立的，优化了流水线。当然，事实证明这并不重要，因为我的 CPU 有乱序执行。另请注意，重新排列需要剥离最后一次迭代。
其他假设
我不是为一般的点积编写代码。该代码适用于特定尺寸，我不处理边缘案例。这只是为了测试 HPC 概念，看看我能达到什么类型的 CPU 使用率。
结果
使用gcc -std=c99 -O2 -m32 -mincoming-stack-boundary=2 -msse3 -mfpmath=sse,387 -masm=intel 编译。我在与平常不同的计算机上。这台计算机有一个i5 540m，它可以在两步英特尔睿频加速后获得2.8 GHz * 4 FLOPS/cycle/core = 11.2 GFLOPS/s per core（两个内核现在都打开了，所以它只有两步……如果我关闭一个内核，则可以提高四步） . 32 位 LINPACK 设置为使用一个线程运行时，速度约为 9.5 GFLOPS/s。
N Total Gflops/s Residual
256 5.580521 1.421085e-014
384 5.734344 -2.842171e-014
512 5.791168 0.000000e+000
640 5.821629 0.000000e+000
768 5.814255 2.842171e-014
896 5.807132 0.000000e+000
1024 5.817208 -1.421085e-013
1152 5.805388 0.000000e+000
1280 5.830746 -5.684342e-014
1408 5.881937 -5.684342e-014
1536 5.872159 -1.705303e-013
1664 5.881536 5.684342e-014
1792 5.906261 -2.842171e-013
1920 5.477966 2.273737e-013
2048 5.620931 0.000000e+000
2176 3.998713 1.136868e-013
2304 3.370095 -3.410605e-013
2432 3.371386 -3.410605e-013
问题 1
我怎样才能做得比这更好？我什至没有接近巅峰表现。我已经将汇编代码优化到了天堂。进一步展开可能会稍微提升一点，但较少展开似乎会降低性能。
问题 2
n > 2048 时，您可以看到性能下降。这是因为我的L1缓存是32KB，当n = 2048和A和B是double时，就占满了整个缓存。任何更大的，它们都是从内存中流式传输的。
我尝试了缓存阻塞（源代码中未显示），但也许我做错了。谁能提供一些代码或解释如何阻止缓存的点积？
源代码
#include <stdio.h>
#include <time.h>
#include <stdlib.h>
#include <string.h>
#include <x86intrin.h>
#include <math.h>
#include <omp.h>
#include <stdint.h>
#include <windows.h>
// computes 8 dot products
#define KERNEL(address) \
"movapd xmm4, XMMWORD PTR [eax+"#address"] \n\t" \
"mulpd xmm7, XMMWORD PTR [edx+48+"#address"] \n\t" \
"addpd xmm2, xmm6 \n\t" \
"movapd xmm5, XMMWORD PTR [eax+16+"#address"] \n\t" \
"mulpd xmm4, XMMWORD PTR [edx+"#address"] \n\t" \
"addpd xmm3, xmm7 \n\t" \
"movapd xmm6, XMMWORD PTR [eax+96+"#address"] \n\t" \
"mulpd xmm5, XMMWORD PTR [edx+16+"#address"] \n\t" \
"addpd xmm0, xmm4 \n\t" \
"movapd xmm7, XMMWORD PTR [eax+112+"#address"] \n\t" \
"mulpd xmm6, XMMWORD PTR [edx+96+"#address"] \n\t" \
"addpd xmm1, xmm5 \n\t"
#define PEELED(address) \
"movapd xmm4, XMMWORD PTR [eax+"#address"] \n\t" \
"mulpd xmm7, [edx+48+"#address"] \n\t" \
"addpd xmm2, xmm6 \n\t" \
"movapd xmm5, XMMWORD PTR [eax+16+"#address"] \n\t" \
"mulpd xmm4, XMMWORD PTR [edx+"#address"] \n\t" \
"addpd xmm3, xmm7 \n\t" \
"mulpd xmm5, XMMWORD PTR [edx+16+"#address"] \n\t" \
"addpd xmm0, xmm4 \n\t" \
"addpd xmm1, xmm5 \n\t"
inline double
__attribute__ ((gnu_inline))
__attribute__ ((aligned(64))) ddot_ref(
int n,
const double* restrict A,
const double* restrict B)
double sum0 = 0.0;
double sum1 = 0.0;
double sum2 = 0.0;
double sum3 = 0.0;
double sum;
for(int i = 0; i < n; i+=4)
sum0 += *(A + i ) * *(B + i );
sum1 += *(A + i+1) * *(B + i+1);
sum2 += *(A + i+2) * *(B + i+2);
sum3 += *(A + i+3) * *(B + i+3);
sum = sum0+sum1+sum2+sum3;
return(sum);
inline double
__attribute__ ((gnu_inline))
__attribute__ ((aligned(64))) ddot_asm
( int n,
const double* restrict A,
const double* restrict B)
double sum;
__asm__ __volatile__
(
"mov eax, %[A] \n\t"
"mov edx, %[B] \n\t"
"mov ecx, %[n] \n\t"
"pxor xmm0, xmm0 \n\t"
"pxor xmm1, xmm1 \n\t"
"pxor xmm2, xmm2 \n\t"
"pxor xmm3, xmm3 \n\t"
"movapd xmm6, XMMWORD PTR [eax+32] \n\t"
"movapd xmm7, XMMWORD PTR [eax+48] \n\t"
"mulpd xmm6, XMMWORD PTR [edx+32] \n\t"
"sar ecx, 7 \n\t"
"sub ecx, 1 \n\t" // peel
"L%=: \n\t"
KERNEL(64 * 0)
KERNEL(64 * 1)
KERNEL(64 * 2)
KERNEL(64 * 3)
KERNEL(64 * 4)
KERNEL(64 * 5)
KERNEL(64 * 6)
KERNEL(64 * 7)
KERNEL(64 * 8)
KERNEL(64 * 9)
KERNEL(64 * 10)
KERNEL(64 * 11)
KERNEL(64 * 12)
KERNEL(64 * 13)
KERNEL(64 * 14)
KERNEL(64 * 15)
"lea eax, [eax+1024] \n\t"
"lea edx, [edx+1024] \n\t"
" \n\t"
"dec ecx \n\t"
"jnz L%= \n\t" // end loop
" \n\t"
KERNEL(64 * 0)
KERNEL(64 * 1)
KERNEL(64 * 2)
KERNEL(64 * 3)
KERNEL(64 * 4)
KERNEL(64 * 5)
KERNEL(64 * 6)
KERNEL(64 * 7)
KERNEL(64 * 8)
KERNEL(64 * 9)
KERNEL(64 * 10)
KERNEL(64 * 11)
KERNEL(64 * 12)
KERNEL(64 * 13)
KERNEL(64 * 14)
PEELED(64 * 15)
" \n\t"
"addpd xmm0, xmm1 \n\t" // summing result
"addpd xmm2, xmm3 \n\t"
"addpd xmm0, xmm2 \n\t" // cascading add
"movapd xmm1, xmm0 \n\t" // copy xmm0
"shufpd xmm1, xmm0, 0x03 \n\t" // shuffle
"addsd xmm0, xmm1 \n\t" // add low qword
"movsd %[sum], xmm0 \n\t" // mov low qw to sum
: // outputs
[sum] "=m" (sum)
: // inputs
[A] "m" (A),
[B] "m" (B),
[n] "m" (n)
: //register clobber
"memory",
"eax","ecx","edx","edi",
"xmm0","xmm1","xmm2","xmm3","xmm4","xmm5","xmm6","xmm7"
);
return(sum);
int main()
// timers
LARGE_INTEGER frequency, time1, time2;
double time3;
QueryPerformanceFrequency(&frequency);
// clock_t time1, time2;
double gflops;
int nmax = 4096;
int trials = 1e7;
double sum, residual;
FILE *f = fopen("soddot.txt","w+");
printf("%16s %16s %16s\n","N","Total Gflops/s","Residual");
fprintf(f,"%16s %16s %16s\n","N","Total Gflops/s","Residual");
for(int n = 256; n <= nmax; n += 128 )
double* A = NULL;
double* B = NULL;
A = _mm_malloc(n*sizeof(*A), 64); if (!A) printf("A failed\n"); return(1);
B = _mm_malloc(n*sizeof(*B), 64); if (!B) printf("B failed\n"); return(1);
srand(time(NULL));
// create arrays
for(int i = 0; i < n; ++i)
*(A + i) = (double) rand()/RAND_MAX;
*(B + i) = (double) rand()/RAND_MAX;
// warmup
sum = ddot_asm(n,A,B);
QueryPerformanceCounter(&time1);
// time1 = clock();
for (int count = 0; count < trials; count++)
// sum = ddot_ref(n,A,B);
sum = ddot_asm(n,A,B);
QueryPerformanceCounter(&time2);
time3 = (double)(time2.QuadPart - time1.QuadPart) / frequency.QuadPart;
// time3 = (double) (clock() - time1)/CLOCKS_PER_SEC;
gflops = (double) (2.0*n*trials)/time3/1.0e9;
residual = ddot_ref(n,A,B) - sum;
printf("%16d %16f %16e\n",n,gflops,residual);
fprintf(f,"%16d %16f %16e\n",n,gflops,residual);
_mm_free(A);
_mm_free(B);
fclose(f);
return(0); // successful completion
编辑：组装说明
点积只是两个数字乘积的重复总和：sum += a[i]*b[i]。 sum 必须在第一次迭代之前初始化为 0。向量化，您一次可以做 2 个求和，必须在最后求和：[sum0 sum1] = [a[i] a[i+1]]*[b[i] b[i+1]]，sum = sum0 + sum1。在 (Intel) 汇编中，这是 3 个步骤（在初始化之后）：
pxor xmm0, xmm0 // accumulator [sum0 sum1] = [0 0]
movapd xmm1, XMMWORD PTR [eax] // load [a[i] a[i+1]] into xmm1
mulpd xmm1, XMMWORD PTR [edx] // xmm1 = xmm1 * [b[i] b[i+1]]
addpd xmm0, xmm1 // xmm0 = xmm0 + xmm1
此时你没有什么特别的，编译器可以想出这个。通常，通过展开代码足够多的时间以使用所有可用的xmm 寄存器（32 位模式下的 8 个寄存器），您通常可以获得更好的性能。因此，如果您将其展开 4 次，则可以使用所有 8 个寄存器 xmm0 到 xmm7。您将有 4 个累加器和 4 个寄存器用于存储 movapd 和 addpd 的结果。同样，编译器可以提出这一点。真正思考的部分是想出一种方法来管道代码，即使 MOV/MUL/ADD 组中的每条指令在不同的寄存器上操作，以便所有 3 条指令同时执行（通常在大多数 CPU）。这就是你击败编译器的方式。因此，您必须对 4x 展开的代码进行模式化才能做到这一点，这可能需要提前加载向量并剥离第一次或最后一次迭代。这就是KERNEL(address)。为方便起见，我制作了 4x 展开的流水线代码的宏。这样，我只需更改 address 就可以轻松地将其展开为 4 的倍数。每个 KERNEL 计算 8 个点积。
【问题讨论】：
您可能希望使用compiler intrinsics 而不是内联汇编代码。它看起来更好。 @tangrs，无论标志如何，它们都不会优化人类手动操作的方式。而且它们速度较慢。我已经试过了。 嗯，这很有趣。我一直认为内在函数与下面的程序集是 1-1 映射的。 @tangrs 我也这么认为。它们通常会生成正确的 MOVPD/MULPD/ADDPD 分组，但我似乎从来没有像他们那样进行重新排序以使每个 MOV/MUL/ADD 在不同的寄存器上工作。具有讽刺意味的是，编译器内在函数为我的矩阵乘法项目生成了一个快速内核，它比我从 Goto 复制的其他一些内核运行得更快。尽管如此，在内部函数案例中仍有改进的空间。 为什么是-O2 而不是-O3？为什么不-march=native？
【参考方案1】：
要回答您的整体问题，您无法使用点积实现最佳性能。
问题是您的 CPU 可以在每个时钟周期执行一次 128 位加载，而要执行点积，您需要在每个时钟周期进行两次 128 位加载。
但是对于大 n 来说，情况比这更糟。你的第二个问题的答案是，点积是内存绑定的，而不是计算绑定的，因此它不能并行处理具有快速内核的大 n。这在why-vectorizing-the-loop-does-not-have-performance-improvement 有更好的解释。这是快速内核并行化的一个大问题。我花了一段时间才弄清楚这一点，但学习这一点非常重要。
实际上很少有基本算法可以充分受益于快速内核的并行化。就 BLAS 算法而言，只有诸如矩阵乘法之类的 Level-3 算法 (O(n^3)) 才能真正受益于并行化。慢核心的情况更好，例如使用 GPU 和 Xeon Phi，因为内存速度和核心速度之间的差异要小得多。
如果您想找到一种算法可以接近小 n 的峰值翻牌，请尝试例如标量 * 向量或标量 * 向量的总和。第一种情况应在每个时钟周期执行一次加载、一次乘法和一次存储，而第二种情况应在每个时钟周期执行一次乘法、一次加法和一次加载。
我在 Knoppix 7.3 32 位的 Core 2 Duo P9600@2.67GHz 上测试了以下代码。 我得到了标量积峰值的 75% 和标量积之和的峰值的 75%。 标量积的 flops/cycle 是 2 和标量积之和现在是 4。
用g++ -msse2 -O3 -fopenmp foo.cpp -ffast-math编译
#include <stdio.h>
#include <stdlib.h>
#include <omp.h>
#include <x86intrin.h>
void scalar_product(double * __restrict a, int n)
a = (double*)__builtin_assume_aligned (a, 64);
double k = 3.14159;
for(int i=0; i<n; i++)
a[i] = k*a[i];
void scalar_product_SSE(double * __restrict a, int n)
a = (double*)__builtin_assume_aligned (a, 64);
__m128d k = _mm_set1_pd(3.14159);
for(int i=0; i<n; i+=8)
__m128d t1 = _mm_load_pd(&a[i+0]);
_mm_store_pd(&a[i],_mm_mul_pd(k,t1));
__m128d t2 = _mm_load_pd(&a[i+2]);
_mm_store_pd(&a[i+2],_mm_mul_pd(k,t2));
__m128d t3 = _mm_load_pd(&a[i+4]);
_mm_store_pd(&a[i+4],_mm_mul_pd(k,t3));
__m128d t4 = _mm_load_pd(&a[i+6]);
_mm_store_pd(&a[i+6],_mm_mul_pd(k,t4));
double scalar_sum(double * __restrict a, int n)
a = (double*)__builtin_assume_aligned (a, 64);
double sum = 0.0;
double k = 3.14159;
for(int i=0; i<n; i++)
sum += k*a[i];
return sum;
double scalar_sum_SSE(double * __restrict a, int n)
a = (double*)__builtin_assume_aligned (a, 64);
__m128d sum1 = _mm_setzero_pd();
__m128d sum2 = _mm_setzero_pd();
__m128d sum3 = _mm_setzero_pd();
__m128d sum4 = _mm_setzero_pd();
__m128d k = _mm_set1_pd(3.14159);
for(int i=0; i<n; i+=8)
__m128d t1 = _mm_load_pd(&a[i+0]);
sum1 = _mm_add_pd(_mm_mul_pd(k,t1),sum1);
__m128d t2 = _mm_load_pd(&a[i+2]);
sum2 = _mm_add_pd(_mm_mul_pd(k,t2),sum2);
__m128d t3 = _mm_load_pd(&a[i+4]);
sum3 = _mm_add_pd(_mm_mul_pd(k,t3),sum3);
__m128d t4 = _mm_load_pd(&a[i+6]);
sum4 = _mm_add_pd(_mm_mul_pd(k,t4),sum4);
double tmp[8];
_mm_storeu_pd(&tmp[0],sum1);
_mm_storeu_pd(&tmp[2],sum2);
_mm_storeu_pd(&tmp[4],sum3);
_mm_storeu_pd(&tmp[6],sum4);
double sum = 0;
for(int i=0; i<8; i++) sum+=tmp[i];
return sum;
int main()
//_MM_SET_FLUSH_ZERO_MODE(_MM_FLUSH_ZERO_ON);
//_mm_setcsr(_mm_getcsr() | 0x8040);
double dtime, peak, flops, sum;
int repeat = 1<<18;
const int n = 2048;
double *a = (double*)_mm_malloc(sizeof(double)*n,64);
double *b = (double*)_mm_malloc(sizeof(double)*n,64);
for(int i=0; i<n; i++) a[i] = 1.0*rand()/RAND_MAX;
dtime = omp_get_wtime();
for(int r=0; r<repeat; r++)
scalar_product_SSE(a,n);
dtime = omp_get_wtime() - dtime;
peak = 2*2.67;
flops = 1.0*n/dtime*1E-9*repeat;
printf("time %f, %f, %f\n", dtime,flops, flops/peak);
//for(int i=0; i<n; i++) a[i] = 1.0*rand()/RAND_MAX;
//sum = 0.0;
dtime = omp_get_wtime();
for(int r=0; r<repeat; r++)
scalar_sum_SSE(a,n);
dtime = omp_get_wtime() - dtime;
peak = 2*2*2.67;
flops = 2.0*n/dtime*1E-9*repeat;
printf("time %f, %f, %f\n", dtime,flops, flops/peak);
//printf("sum %f\n", sum);
【讨论】：
我用-msse3为我的电脑编译了这个并更新了peak，我只得到了大约25%的CPU使用率。我还尝试将 sum 函数展开 4 次而不是 3 次，但由于某种原因它运行得非常慢（几乎到了我认为它只是挂起的地步）。至于最初的问题，我在其中一个 cmets 中写道，根据 A Fog，movapd 如何至少需要 2 个周期，因此 movapd/mulpd/addpd 的每个分组（自它们在不同的寄存器上运行）将需要 2 个周期进行 4 个触发器，所以我能做的最好的就是 50% 峰值（对我来说是 5.6 GFLOPS/s）
对不起，我的 unroll3 功能不好。我使用内在函数重新完成了这个，并在与您的非常相似的 Core 2 Duo 系统上对其进行了测试。我现在得到了大约 50% 的峰值。我不知道为什么......我用新代码更新了我的答案。展开三个应该就足够了，但我用四个来说服（所以 n 只需要是 2 的倍数）。
OK thx，我还想知道展开 3 的正确性，因为数组长度是 2 的幂。它不是严格正确的，但对于性能测量来说它很好（我们会错过最后 3 个元素）。我很快就会试一试。
@matmul，该函数还有其他问题。我通常自己对代码进行矢量化。我应该从一开始就这样做。
我将在今天或明天发布它并附上解释，我试图消除打包的开销（这就是阻止我接近 95% 的原因）。但是，如果您想要预览，我基本上在循环 .L16 上解决了 OpenBLAS 内核中的概念 gemm_kernel_2x4_PENRYN.S 并对其进行了一些改进。诀窍似乎是（1）数据重用/消除负载和（2）一旦你发出第一个mul，剩下的muls在循环中的每个周期发出一次。这样你只受adds的限制。以上是关于如何使用点积获得峰值 CPU 性能？的主要内容，如果未能解决你的问题，请参考以下文章