如何在 Skylake 架构上最大化 sqrt-heavy loop 的指令级并行性?

为了向自己介绍x86内部函数(并在较小程度上缓存友好性(，我明确矢量化了一些用于基于RBF(径向基函数(的网格变形的代码。发现 vsqrtpd 是主要瓶颈后，我想知道是否/如何进一步掩盖其延迟。这是标量计算内核：

for(size_t i=0; i<nPt; ++i)
{
double xi = X[i], yi = X[i+nPt], zi = X[i+2*nPt];
for(size_t j=0; j<nCP; ++j)
{
// compute distance from i to j
double d = sqrt(pow(xi-Xcp[   j   ],2)+
pow(yi-Xcp[ j+nCP ],2)+
pow(zi-Xcp[j+2*nCP],2));
// compute the RBF kernel coefficient
double t = max(0.0,1.0-d);
t = pow(t*t,2)*(1.0+4.0*d);
// update coordinates
for(size_t k=0; k<nDim; ++k) X[i+k*nPt] += t*Ucp[j+k*nCP];
}
}

nPt 是目标坐标的数量，它远大于源坐标/位移的数量 nCP。后者适合 L3，因此最里面的环路总是在源点上。

第一个优化步骤是同时处理 4 个目标点。源点数据仍然通过标量加载访问，然后通过广播进行访问。
第二步是通过阻塞环路来靶向 L1，阻塞 i 环路在某种程度上比阻塞 j 环路重要得多，后者只提供了轻微的改进。最内层循环仍然超过 j 以减少负载/存储。
第三是加载 4 个控制点并使用随机/排列来遍历 i-j 的 4 种组合，而不是使用广播。
第四，在观察到省略平方根可使速度提高 1.5 倍(大约是 i7-7700 上大型 LLT 的 FP 性能的 70%(之后，将 4 个寄存器专用于 4 平方根的计算，以(也许？(允许进行其他一些计算......与第三步相比，改善 1%。

当前代码

void deform(size_t nPt, size_t nCP, const double* Xcp, const double* Ucp, double* X)
{
const size_t SIMDLEN = 4;
// tile ("cache block") sizes
const size_t TILEH = 512;
const size_t TILEW = 256;
// fill two registers with the constants we need
__m256d vone  = _mm256_set1_pd(1.0),
vfour = _mm256_set1_pd(4.0);
// explicitly vectorized (multiple i's at a time) and blocked
// outer most loop over sets of #TILEH points
for(size_t i0=0; i0<nPt; i0+=TILEH)
{
// displacement buffer, due to tiling, coordinates cannot be modified in-place
alignas(64) double U[3*TILEH*sizeof(double)];
// zero the tile displacements
for(size_t k=0; k<3*TILEH; k+=SIMDLEN)
_mm256_store_pd(&U[k], _mm256_setzero_pd());
// stop point for inner i loop
size_t iend = min(i0+TILEH,nPt);
// second loop over sets of #TILEW control points
for(size_t j0=0; j0<nCP; j0+=TILEW)
{
// stop point for inner j loop
size_t jend = min(j0+TILEW,nCP);
// inner i loop, over #TILEH points
// vectorized, operate on #SIMDLEN points at a time
for(size_t i=i0; i<iend; i+=SIMDLEN)
{
// coordinates and displacements of points i
__m256d wi,
xi = _mm256_load_pd(&X[   i   ]),
yi = _mm256_load_pd(&X[ i+nPt ]),
zi = _mm256_load_pd(&X[i+2*nPt]),
ui = _mm256_load_pd(&U[    i-i0    ]),
vi = _mm256_load_pd(&U[ i-i0+TILEH ]);
wi = _mm256_load_pd(&U[i-i0+2*TILEH]);
// inner j loop, over #TILEW control points, vectorized loads
for(size_t j=j0; j<jend; j+=SIMDLEN)
{
// coordinates of points j, and an aux var
__m256d t,
xj = _mm256_load_pd(&Xcp[   j   ]),
yj = _mm256_load_pd(&Xcp[ j+nCP ]),
zj = _mm256_load_pd(&Xcp[j+2*nCP]);
// compute the possible 4 distances from i to j...
#define COMPUTE_DIST(D) __m256d                         
D = _mm256_sub_pd(xi,xj);  D = _mm256_mul_pd(D,D);      
t = _mm256_sub_pd(yi,yj);  D = _mm256_fmadd_pd(t,t,D);  
t = _mm256_sub_pd(zi,zj);  D = _mm256_fmadd_pd(t,t,D);  
D = _mm256_sqrt_pd(D)
// ...by going through the different permutations
#define SHUFFLE(FUN,IMM8)   
xj = FUN(xj,xj,IMM8);       
yj = FUN(yj,yj,IMM8);       
zj = FUN(zj,zj,IMM8)
COMPUTE_DIST(d0);
SHUFFLE(_mm256_shuffle_pd,0b0101);
COMPUTE_DIST(d1);
SHUFFLE(_mm256_permute2f128_pd,1);
COMPUTE_DIST(d2);
SHUFFLE(_mm256_shuffle_pd,0b0101);
COMPUTE_DIST(d3);
// coordinate registers now hold the displacements
xj = _mm256_load_pd(&Ucp[   j   ]),
yj = _mm256_load_pd(&Ucp[ j+nCP ]);
zj = _mm256_load_pd(&Ucp[j+2*nCP]);
// coefficients for each set of distances...
#define COMPUTE_COEFF(C)                                
t = _mm256_min_pd(vone,C);  t = _mm256_sub_pd(vone,t);  
t = _mm256_mul_pd(t,t);     t = _mm256_mul_pd(t,t);     
C = _mm256_fmadd_pd(vfour,C,vone);                      
C = _mm256_mul_pd(t,C)
// ...+ update i point displacements
#define UPDATE_DISP(C)          
COMPUTE_COEFF(C);               
ui = _mm256_fmadd_pd(C,xj,ui);  
vi = _mm256_fmadd_pd(C,yj,vi);  
wi = _mm256_fmadd_pd(C,zj,wi)
UPDATE_DISP(d0);
SHUFFLE(_mm256_shuffle_pd,0b0101);
UPDATE_DISP(d1);
SHUFFLE(_mm256_permute2f128_pd,1);
UPDATE_DISP(d2);
SHUFFLE(_mm256_shuffle_pd,0b0101);
UPDATE_DISP(d3);
}
// store updated displacements
_mm256_store_pd(&U[    i-i0    ], ui);
_mm256_store_pd(&U[ i-i0+TILEH ], vi);
_mm256_store_pd(&U[i-i0+2*TILEH], wi);
}
}
// add tile displacements to the coordinates
for(size_t k=0; k<3; ++k)
{
for(size_t i=i0; i<iend; i+=SIMDLEN)
{
__m256d
x = _mm256_load_pd(&X[i+k*nPt]),
u = _mm256_load_pd(&U[i-i0+k*TILEH]);
x = _mm256_add_pd(x,u);
_mm256_stream_pd(&X[i+k*nPt], x);
}
}
}
}

那么我还能对它做些什么呢？或者，我做错了什么？

谢谢 P. 戈麦斯

首先检查性能计数器是否arith.divider_active为 ~= 内核时钟周期。

函数运行时间的 98% 可以通过取平方根数和运算吞吐量来解释。

或者这也行得通。

如果是这种情况，那么您正在饱和(未完全流水线化(分频器吞吐量，并且仅公开更多 ILP 就没有多少收益。

算法更改是您获得任何东西的唯一真正机会，例如避免一些sqrt操作或使用单精度。

单精度可免费为每个矢量提供 2 倍的工作量。但对于 sqrt 密集型工作负载，还有一个额外的好处：每个向量vsqrtps吞吐量通常优于vsqrtpd。 Skylake就是这种情况：每6个周期一个，而vsqrtpd每9到12个周期一个。这可能会将瓶颈从sqrt/divide单元转移到前端或FMA单元。

评论中提出了vrsqrtps。这是值得考虑的(如果可以选择单精度(，但是当您需要牛顿拉夫森迭代来获得足够的精度时，这不是一个明显的胜利。没有Newton Raphson的裸x * rsqrtps(x)可能太不准确了(并且需要一个cmp/AND来解决x==0.0(，但是NR迭代可能需要太多额外的FMA uops而不值得。

(带vrsqrt14ps/pd的AVX512在近似方面具有更高的精度，但通常仍然不足以在没有牛顿的情况下使用。但有趣的是，它确实存在双精度。当然，如果你在Xeon Phi上，sqrt非常慢，你打算使用AVX512ERvrsqrt28pd+ Newton，或者只是单独vrsqrt28ps。

上次我为 Skylake 调整了一个包含多项式近似 sqrt 的函数时，快速近似倒数是不值得的。硬件单精度sqrt是最佳选择，它为我们提供了所需的精度(我们甚至没有考虑需要double(。不过，sqrt 操作之间的工作比你多。

相关内容

最新更新

热门标签：