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

Question

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

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的源坐标个数coordinates/displacements。后者适合 L3，因此最内层的循环总是在源点上。

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

Current code

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.戈麦斯

Answer 1

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

98% of the function runtime can be explained by taking the number of square roots and the operation throughput.

或者那也行。

如果是这种情况，您就会使（未完全流水线化的）分频器吞吐量饱和，并且仅仅暴露更多的 ILP 就没有多少好处了。

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

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

vrsqrtps 已在评论中建议。这值得考虑（如果单精度是一个选项），但当您需要 Newton Raphson 迭代以获得足够的精度时，这并不是一个明显的胜利。没有 Newton Raphson 的裸 x * rsqrtps(x) 可能太不准确（并且需要 cmp/AND 来解决 x==0.0），但是 NR 迭代可能需要太多额外的 FMA 微指令是值得的。

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

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

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

How to maximise instruction level parallelism of sqrt-heavy-loop on skylake architecture?

c++

optimization

x86

intrinsics

avx