如果在 Xeon Phi 上编译时不知道循环计数，则性能下降
Posted
技术标签:
【中文标题】如果在 Xeon Phi 上编译时不知道循环计数，则性能下降【英文标题】：Performance degradation if loop count is not known at compile time on Xeon Phi 【发布时间】：2014-11-04 16:39:10 【问题描述】：我正在创建一个简单的矩阵乘法程序，在英特尔至强融核架构上运行。
在多次尝试自动矢量化后，为了获得更好的性能，我不得不使用 Intel Intrinsics。
到目前为止，矩阵大小是由源代码中的#define 给出的，但是当我尝试在运行时给出它时，性能会大幅下降。
源码如下：
#include <stdio.h>
#include <stdlib.h>
#include <time.h>
#include <math.h>
#include <stddef.h>
#include <chrono>
#include <ctime>
#include <mmintrin.h>
#include <xmmintrin.h> // SSE
#include <pmmintrin.h> // SSE2
#include <emmintrin.h> // SSE3
#include <immintrin.h>
#include <zmmintrin.h>
#define ALIGNMENT 64
#ifndef SIZE
#define SIZE 960
#endif
#define vZero(c) (c) = _mm512_setzero_pd();
#define start_time() \
auto start = std::chrono::high_resolution_clock::now();
/** Shows the elapsed time. See start_time for usage*/
#define elapsed_time(STRING) \
auto elapsed = std::chrono::high_resolution_clock::now() - start; \
long long microseconds = std::chrono::duration_cast<std::chrono::microseconds>(elapsed).count(); \
printf(#STRING":%lld\n", microseconds);
void recTranspose(double *__restrict__ a, double *__restrict__ aT, const int n, const int k, const int lda, const int ldat)
if (n*k <= 128)
for(int i = 0; i < n; i++)
for(int j = 0; j < k; j++)
aT[j*ldat+i] = a[i*lda+j];
//printf("Reached _|_");
return;
if(k > n)
recTranspose(a, aT, n, (k+1)/2, lda, ldat);
recTranspose(&a[(k+1)/2], &aT[(k+1)/2*ldat], n, k-((k+1)/2), lda, ldat);
else
recTranspose(a, aT, (n+1)/2, k, lda, ldat);
recTranspose(&a[(n+1)/2*lda], &aT[(n+1)/2], n- (n+1)/2, k, lda, ldat);
/** Calculates 8 cols and 30 rows of c.*/
inline void eightbythirty(double *__restrict__ a, double *__restrict__ b, double * __restrict__ c, const int size)
__m512d c0, c1, c2, c3, c4, c5, c6, c7, c8, c9;
__m512d c10, c11, c12, c13, c14, c15, c16, c17, c18, c19;
__m512d c20, c21, c22, c23, c24, c25, c26, c27, c28, c29;
vZero(c0); vZero(c1); vZero(c2); vZero(c3); vZero(c4); vZero(c5);
vZero(c6); vZero(c7); vZero(c8); vZero(c9); vZero(c10); vZero(c11);
vZero(c12); vZero(c13); vZero(c14); vZero(c15); vZero(c16); vZero(c17);
vZero(c18); vZero(c19); vZero(c20); vZero(c21); vZero(c22); vZero(c23);
vZero(c24); vZero(c25); vZero(c26); vZero(c27); vZero(c28); vZero(c29);
__assume_aligned(a, ALIGNMENT);
__assume_aligned(b, ALIGNMENT);
__assume_aligned(c, ALIGNMENT);
__assume(size%16==0);
for(int i = 0; i < size; i++)
const __m512d bv = _mm512_load_pd(b+i*size);
c0 = _mm512_fmadd_pd(_mm512_extload_pd(a+i*size+0, _MM_UPCONV_PD_NONE, _MM_BROADCAST_1X8, 0), bv, c0);
c1 = _mm512_fmadd_pd(_mm512_extload_pd(a+i*size+1, _MM_UPCONV_PD_NONE, _MM_BROADCAST_1X8, 0), bv, c1);
c2 = _mm512_fmadd_pd(_mm512_extload_pd(a+i*size+2, _MM_UPCONV_PD_NONE, _MM_BROADCAST_1X8, 0), bv, c2);
c3 = _mm512_fmadd_pd(_mm512_extload_pd(a+i*size+3, _MM_UPCONV_PD_NONE, _MM_BROADCAST_1X8, 0), bv, c3);
c4 = _mm512_fmadd_pd(_mm512_extload_pd(a+i*size+4, _MM_UPCONV_PD_NONE, _MM_BROADCAST_1X8, 0), bv, c4);
c5 = _mm512_fmadd_pd(_mm512_extload_pd(a+i*size+5, _MM_UPCONV_PD_NONE, _MM_BROADCAST_1X8, 0), bv, c5);
c6 = _mm512_fmadd_pd(_mm512_extload_pd(a+i*size+6, _MM_UPCONV_PD_NONE, _MM_BROADCAST_1X8, 0), bv, c6);
c7 = _mm512_fmadd_pd(_mm512_extload_pd(a+i*size+7, _MM_UPCONV_PD_NONE, _MM_BROADCAST_1X8, 0), bv, c7);
c8 = _mm512_fmadd_pd(_mm512_extload_pd(a+i*size+8, _MM_UPCONV_PD_NONE, _MM_BROADCAST_1X8, 0), bv, c8);
c9 = _mm512_fmadd_pd(_mm512_extload_pd(a+i*size+9, _MM_UPCONV_PD_NONE, _MM_BROADCAST_1X8, 0), bv, c9);
c10 = _mm512_fmadd_pd(_mm512_extload_pd(a+i*size+10, _MM_UPCONV_PD_NONE, _MM_BROADCAST_1X8, 0),bv, c10);
c11 = _mm512_fmadd_pd(_mm512_extload_pd(a+i*size+11, _MM_UPCONV_PD_NONE, _MM_BROADCAST_1X8, 0),bv, c11);
c12 = _mm512_fmadd_pd(_mm512_extload_pd(a+i*size+12, _MM_UPCONV_PD_NONE, _MM_BROADCAST_1X8, 0), bv, c12);
c13 = _mm512_fmadd_pd(_mm512_extload_pd(a+i*size+13, _MM_UPCONV_PD_NONE, _MM_BROADCAST_1X8, 0), bv, c13);
c14 = _mm512_fmadd_pd(_mm512_extload_pd(a+i*size+14, _MM_UPCONV_PD_NONE, _MM_BROADCAST_1X8, 0), bv, c14);
c15 = _mm512_fmadd_pd(_mm512_extload_pd(a+i*size+15, _MM_UPCONV_PD_NONE, _MM_BROADCAST_1X8, 0), bv, c15);
c16 = _mm512_fmadd_pd(_mm512_extload_pd(a+i*size+16, _MM_UPCONV_PD_NONE, _MM_BROADCAST_1X8, 0), bv, c16);
c17 = _mm512_fmadd_pd(_mm512_extload_pd(a+i*size+17, _MM_UPCONV_PD_NONE, _MM_BROADCAST_1X8, 0), bv, c17);
c18 = _mm512_fmadd_pd(_mm512_extload_pd(a+i*size+18, _MM_UPCONV_PD_NONE, _MM_BROADCAST_1X8, 0), bv, c18);
c19 = _mm512_fmadd_pd(_mm512_extload_pd(a+i*size+19, _MM_UPCONV_PD_NONE, _MM_BROADCAST_1X8, 0), bv, c19);
c20 = _mm512_fmadd_pd(_mm512_extload_pd(a+i*size+20, _MM_UPCONV_PD_NONE, _MM_BROADCAST_1X8, 0), bv, c20);
c21 = _mm512_fmadd_pd(_mm512_extload_pd(a+i*size+21, _MM_UPCONV_PD_NONE, _MM_BROADCAST_1X8, 0), bv, c21);
c22 = _mm512_fmadd_pd(_mm512_extload_pd(a+i*size+22, _MM_UPCONV_PD_NONE, _MM_BROADCAST_1X8, 0), bv, c22);
c23 = _mm512_fmadd_pd(_mm512_extload_pd(a+i*size+23, _MM_UPCONV_PD_NONE, _MM_BROADCAST_1X8, 0), bv, c23);
c24 = _mm512_fmadd_pd(_mm512_extload_pd(a+i*size+24, _MM_UPCONV_PD_NONE, _MM_BROADCAST_1X8, 0), bv, c24);
c25 = _mm512_fmadd_pd(_mm512_extload_pd(a+i*size+25, _MM_UPCONV_PD_NONE, _MM_BROADCAST_1X8, 0), bv, c25);
c26 = _mm512_fmadd_pd(_mm512_extload_pd(a+i*size+26, _MM_UPCONV_PD_NONE, _MM_BROADCAST_1X8, 0), bv, c26);
c27 = _mm512_fmadd_pd(_mm512_extload_pd(a+i*size+27, _MM_UPCONV_PD_NONE, _MM_BROADCAST_1X8, 0), bv, c27);
c28 = _mm512_fmadd_pd(_mm512_extload_pd(a+i*size+28, _MM_UPCONV_PD_NONE, _MM_BROADCAST_1X8, 0), bv, c28);
c29 = _mm512_fmadd_pd(_mm512_extload_pd(a+i*size+29, _MM_UPCONV_PD_NONE, _MM_BROADCAST_1X8, 0), bv, c29);
_mm512_storenr_pd(c+0*size, c0);
_mm512_storenr_pd(c+1*size, c1);
_mm512_storenr_pd(c+2*size, c2);
_mm512_storenr_pd(c+3*size, c3);
_mm512_storenr_pd(c+4*size, c4);
_mm512_storenr_pd(c+5*size, c5);
_mm512_storenr_pd(c+6*size, c6);
_mm512_storenr_pd(c+7*size, c7);
_mm512_storenr_pd(c+8*size, c8);
_mm512_storenr_pd(c+9*size, c9);
_mm512_storenr_pd(c+10*size, c10);
_mm512_storenr_pd(c+11*size, c11);
_mm512_storenr_pd(c+12*size, c12);
_mm512_storenr_pd(c+13*size, c13);
_mm512_storenr_pd(c+14*size, c14);
_mm512_storenr_pd(c+15*size, c15);
_mm512_storenr_pd(c+16*size, c16);
_mm512_storenr_pd(c+17*size, c17);
_mm512_storenr_pd(c+18*size, c18);
_mm512_storenr_pd(c+19*size, c19);
_mm512_storenr_pd(c+20*size, c20);
_mm512_storenr_pd(c+21*size, c21);
_mm512_storenr_pd(c+22*size, c22);
_mm512_storenr_pd(c+23*size, c23);
_mm512_storenr_pd(c+24*size, c24);
_mm512_storenr_pd(c+25*size, c25);
_mm512_storenr_pd(c+26*size, c26);
_mm512_storenr_pd(c+27*size, c27);
_mm512_storenr_pd(c+28*size, c28);
_mm512_storenr_pd(c+29*size, c29);
int main(int argc, const char ** argv)
#ifdef SIZES
const int size = SIZE;
#else
const int size = atoi(argv[1]);
#endif
void* p = malloc((sizeof(double)*5*size*size) + ALIGNMENT-1);
double *__restrict__ a = (double*)(((size_t)p + ALIGNMENT-1) / ALIGNMENT * ALIGNMENT);
double *__restrict__ aT = (double*) a+size*size;
double *__restrict__ b = aT+size*size;
double *__restrict__ c = b+size*size;
double *__restrict__ d = c+size*size;
srand(time(NULL));
for(int i = 0; i < size; i++)
for(int j = 0; j < size; j++)
a[i*size+j] = (double) (rand()%20);
for(int j2=0; j2<size; j2++)
c[i*size+j2] = 0.0;
for(int i = 0; i < size; i++)
for(int j = 0; j < size; j++)
b[i*size+j] = (double) (rand()%20);
start_time();
recTranspose(a, aT, size, size, size, size);
for(int i = 0; i < size; i+=30)
for(int j = 0; j < size; j+=8)
eightbythirty(&aT[i], &b[j], &c[i*size+j], size);
elapsed_time();
double gflops = 2.0*size*size*size*1.0e-03/(microseconds);
printf("Gflops: %f\n", gflops);
for(int i = 0; i < size; i++)
for(int j = 0; j < size; j++)
double s = 0;
for(int u = 0; u < size; u++)
s += a[i*size+u] * b[u*size+j];
d[i*size+j] = s;
int error = 0;
for(int i = 0; i < size; i++)
for(int j = 0; j < size; j++)
if(abs(c[i*size+j] - d[i*size+j]) > 1)
printf("Error at %d %d , %f instead of %f\n", i, j, c[i*size+j], d[i*size+j]);
error++;
if(error > 16) return 0;
printf("OK\n");
例如，尺寸为 960（目前它仅适用于尺寸为 30*8 的倍数）：
如果我使用给定大小的编译时间进行编译：icc -mmic -O3 -restrict -std=c++11 -DSIZES -DSIZE=960 mmul.cpp -o mmul.o
经过时间：0.460745s 触发器：3.840458
如果我使用运行时给定大小进行编译：icc -mmic -O3 -restrict -std=c++11 mmul.cpp -o mmul.o
经过时间：2.204564s Gflops：0.802640
我认为这可能是 icc 无法识别内存访问模式的预取问题。查看生成的 asm 源代码，“编译时”版本中 vprefetch 指令的数量要高得多。
有趣的事实：在编译时版本中检查乘法的正确结果（代码末尾的两个 for 循环，第 178-197 行）要慢得多！
有什么想法吗？我尝试了#pragma loop_count，但它似乎没用，而且手动内部预取似乎也不是很有效。
提前感谢您的任何回答。
问候， 卢卡
【问题讨论】：
【参考方案1】：计算机科学的基本定理指出，任何问题都可以通过另一层间接来解决。
我们的想法是将您的代码保留为固定大小的循环，并编写代码以调度到正确的固定大小的循环。
首先将eightbythirty 改为这样：
template<int size>
inline void eightbythirty(double *__restrict__ a, double *__restrict__ b, double * __restrict__ c)
里面有相同的实现。您可以将其放入 namespace details，因为它通常不面向用户。
接下来，包装它：
inline void eightbythirty(double *__restrict__ a, double *__restrict__ b, double * __restrict__ c, const int size_divided_by_240)
int size=size_divided_by_240;
switch( size&7 )
case 0: break;
case 01: eightbythirty<01>(a,b,c); break;
case 02: eightbythirty<02>(a,b,c); break;
case 03: eightbythirty<03>(a,b,c); break;
case 04: eightbythirty<04>(a,b,c); break;
case 05: eightbythirty<05>(a,b,c); break;
case 06: eightbythirty<06>(a,b,c); break;
case 07: eightbythirty<07>(a,b,c); break;
a+=(size&7)*8*30;
b+=(size&7)*8*30;
c+=(size&7)*8*30;
switch( (size>>3)&7 )
case 0: break;
case 01: eightbythirty<1*8>(a,b,c); break;
case 02: eightbythirty<2*8>(a,b,c); break;
case 03: eightbythirty<3*8>(a,b,c); break;
case 04: eightbythirty<4*8>(a,b,c); break;
case 05: eightbythirty<5*8>(a,b,c); break;
case 06: eightbythirty<6*8>(a,b,c); break;
case 07: eightbythirty<7*8>(a,b,c); break;
a += (size&(7<<3))*8*30;
b += (size&(7<<3))*8*30;
c += (size&(7<<3))*8*30;
switch( (size>>6)&7 )
case 0: break;
case 01: eightbythirty<1*8*8>(a,b,c); break;
case 02: eightbythirty<2*8*8>(a,b,c); break;
case 03: eightbythirty<3*8*8>(a,b,c); break;
case 04: eightbythirty<4*8*8>(a,b,c); break;
case 05: eightbythirty<5*8*8>(a,b,c); break;
case 06: eightbythirty<6*8*8>(a,b,c); break;
case 07: eightbythirty<7*8*8>(a,b,c); break;
default:
a += (size&(7<<6))*8*30;
b += (size&(7<<6))*8*30;
c += (size&(7<<6))*8*30;
int steps = size/8/8/8;
for( int i = 0; i < steps; ++i )
eightbythirty<512>(a+512*i, b+512*i, c+512*i);
这会将您的输入大小分成 3 位块。然后它调用固定大小的实现。上面代码中出现了 4 个分支，其中大部分是简单的跳转表，对于小于 512*8*30 的值。对于大于该值的值，事情大多以 512*8*30 的块完成。
7*3+1 = 实现了你的原始函数的 22 个实现，每个都有一个常量 size，因此编译器可以完全优化它们。
这通常可以通过元编程完成，但不值得一次性使用。
当我调用<int size> 版本的eightbythirty 时，我可能在上面的代码中遗漏了一些*(8*30)。
【讨论】：
The fundamental theorem of computer science states that any problem can be solved with another layer of indirection。除了间接的问题太多;)
lol :) 无论如何，我将“接受”标记为这个答案，因为使用相同的技术（但不是完全相同的代码）它有点解决了我的问题。唯一剩下的问题是，在此之后我无法运行该程序。不幸的是，我的提交非常接近，所以我现在无法向您展示代码中的所有更改。当我在接下来的几天里找到一些时间时，我会更多地澄清我在做什么（总是假设有人感兴趣：P）【参考方案2】：
要检查它是否是预取优化结果，您可以尝试编译您定义大小的版本：
icc -mmic -O3 -restrict -std=c++11 -no-opt-prefetch -DSIZES -DSIZE=960 mmul.cpp -o mmul.o
这将禁用软件预取插入。这是一个合理的假设，因为知道数据集的大小将允许更有效的预取。
如果您的问题是预取，您应该尝试更积极的 compiler option -opt-prefetch=4。另一个想法是使用 #pragma prefetch。请参阅 here 了解更多信息。
您是否也尝试过#pragma loop_count。检查here。
您可以使用 vTune 分析您的代码并检查您的 L1 缓存命中率。对于矩阵乘法，它应该非常高。还要确保使用缓存阻塞。
我还建议您使用大量迭代来运行您的基准测试。这将涵盖缓存预热。
我这边的另一个建议是使用数学库来完成这项任务。我看到你的代码是单核的，所以这将导致所有内核的性能约为 250 GFLOPS。那将是大约 25-30% 的效率，这不是很好。 使用 MKL 会带来更高的效率。
但是，出于教育目的，这个练习是可以的。
【讨论】：
感谢您的回答。是的，减少预取会导致性能崩溃（4.310070 秒，0.410544 GFLOPS）。但我不知道如何以用户定义的大小运行程序并获得良好的性能。更高数量的迭代似乎没有任何区别。关于在所有内核上运行它，遗憾的是可扩展性并不是最好的，所以我得到了大约 210 GFLOPS（238 个线程），但还有另一个故事......当然，正如你猜到的那样，这是出于教育目的，所以我不能使用数学库。 您是否尝试过 -opt-prefetch=4。在我的帖子中也添加了此信息。【参考方案3】：这个问题也发布在https://software.intel.com/en-us/forums/topic/534856。 software.intel.com 上的回复来自不同的专家组。
【讨论】：
是的，我希望不是问题。我还在那里发布了 asm 源代码，但我没有像这里那样搜索解决方法，而是试图了解英特尔编译器是否有问题。 发帖到两个社区都很棒！更多地关注问题。确保感兴趣的各方了解其他线程也很好。 :-)以上是关于如果在 Xeon Phi 上编译时不知道循环计数，则性能下降的主要内容，如果未能解决你的问题，请参考以下文章