为了介绍自己的x86内部函数（以及更次要的缓存友好性），我明确地将我用于基于径向基函数的网格变形的一些代码矢量化。
在发现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的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. Gomes
arith.divider_active与核心时钟周期的性能计数器？如果是这样，你正在饱和（不完全流水线化的）除法器吞吐量，除非你可以避免一些平方根操作，否则就没有太多可获得的了。 - Peter Cordesx * rsqrtps(x)而不使用牛顿迭代法可能太不准确（并且在x==0时会出错），但NR迭代需要太多额外的FMA uops，不值得。Skylake具有惊人的单精度vsqrtps吞吐量：每6个周期一个（相对于每9到12个周期的vsqrtpd）。上次我调整了一个包括多项式逼近的平方根函数用于Skylake，快速逼近倒数并不值得。但是，将其转换为float可以使您在这里获得超过2倍的加速。 - Peter Cordesui、vi和wi设置为零，并在循环结束时添加前一个值（不应该对性能造成任何影响）。 - chtz