这个 SSE2 转置有什么问题？

Question

我正在尝试转换此代码：

double *pB = b[voiceIndex];
double *pC = c[voiceIndex];
double phase = mPhase;
double bp0 = mNoteFrequency * mHostPitch;

for (int sampleIndex = 0; sampleIndex < blockSize; sampleIndex++) {
    // some other code (that will use phase, like sin(phase))

    phase += std::clamp(radiansPerSample * (bp0 * pB[sampleIndex] + pC[sampleIndex]), 0.0, PI);
}

mPhase = phase;

在 SSE2 中，试图加速整个块（经常调用）。我正在使用带有 Fast optimizazion 标志的 MSVC，但自动矢量化非常糟糕。因为我也在学习矢量化，所以我觉得这是一个很好的挑战。

所以我对上面的公式进行了简化，比如：

radiansPerSampleBp0 = radiansPerSample * bp0;
phase += std::clamp(radiansPerSampleBp0 * pB[sampleIndex] + radiansPerSample * pC[sampleIndex]), 0.0, PI);

可以将其静音为串行依赖项，例如：

phase[0] += (radiansPerSampleBp0 * pB[0] + radiansPerSample * pC[0])
phase[1] += (radiansPerSampleBp0 * pB[1] + radiansPerSample * pC[1]) + (radiansPerSampleBp0 * pB[0] + radiansPerSample * pC[0])

phase[2] += (radiansPerSampleBp0 * pB[2] + radiansPerSample * pC[2]) + (radiansPerSampleBp0 * pB[1] + radiansPerSample * pC[1])
phase[3] += (radiansPerSampleBp0 * pB[3] + radiansPerSample * pC[3]) + (radiansPerSampleBp0 * pB[2] + radiansPerSample * pC[2])

phase[4] += (radiansPerSampleBp0 * pB[4] + radiansPerSample * pC[4]) + (radiansPerSampleBp0 * pB[3] + radiansPerSample * pC[3])
phase[5] += (radiansPerSampleBp0 * pB[5] + radiansPerSample * pC[5]) + (radiansPerSampleBp0 * pB[4] + radiansPerSample * pC[4])

因此，我做的代码：

double *pB = b[voiceIndex];
double *pC = c[voiceIndex];
double phase = mPhase;
double bp0 = mNoteFrequency * mHostPitch;

__m128d v_boundLower = _mm_set1_pd(0.0);
__m128d v_boundUpper = _mm_set1_pd(PI);
__m128d v_radiansPerSampleBp0 = _mm_set1_pd(mRadiansPerSample * bp0);
__m128d v_radiansPerSample = _mm_set1_pd(mRadiansPerSample);

__m128d v_pB0 = _mm_load_pd(pB);
v_pB0 = _mm_mul_pd(v_pB0, v_radiansPerSampleBp0);
__m128d v_pC0 = _mm_load_pd(pC);
v_pC0 = _mm_mul_pd(v_pC0, v_radiansPerSample);

__m128d v_pB1 = _mm_setr_pd(0.0, pB[0]);
v_pB1 = _mm_mul_pd(v_pB1, v_radiansPerSampleBp0);
__m128d v_pC1 = _mm_setr_pd(0.0, pC[0]);
v_pC1 = _mm_mul_pd(v_pC1, v_radiansPerSample);

__m128d v_phase = _mm_set1_pd(phase);
__m128d v_phaseAcc;

for (int sampleIndex = 0; sampleIndex < blockSize; sampleIndex += 2, pB += 2, pC += 2) {
    // some other code (that will use phase, like sin(phase))

    v_phaseAcc = _mm_add_pd(v_pB0, v_pC0);
    v_phaseAcc = _mm_max_pd(v_phaseAcc, v_boundLower);
    v_phaseAcc = _mm_min_pd(v_phaseAcc, v_boundUpper);
    v_phaseAcc = _mm_add_pd(v_phaseAcc, v_pB1);
    v_phaseAcc = _mm_add_pd(v_phaseAcc, v_pC1);
    v_phase = _mm_add_pd(v_phase, v_phaseAcc);

    v_pB0 = _mm_load_pd(pB + 2);
    v_pB0 = _mm_mul_pd(v_pB0, v_radiansPerSampleBp0);
    v_pC0 = _mm_load_pd(pC + 2);
    v_pC0 = _mm_mul_pd(v_pC0, v_radiansPerSample);

    v_pB1 = _mm_load_pd(pB + 1);
    v_pB1 = _mm_mul_pd(v_pB1, v_radiansPerSampleBp0);
    v_pC1 = _mm_load_pd(pC + 1);
    v_pC1 = _mm_mul_pd(v_pC1, v_radiansPerSample);
}

mPhase = v_phase.m128d_f64[blockSize % 2 == 0 ? 1 : 0];

但是，不幸的是，在求和 "steps" 之后，每个相位值的结果变得非常不同。试过调试，但我真的找不到问题所在。

此外，它不是真的 "fast" 而不是旧版本。

你能识别故障吗？你将如何加速代码？

这是完整的代码，如果你想检查两个不同的输出：

#include <iostream>
#include <algorithm>
#include <immintrin.h>
#include <emmintrin.h>

#define PI 3.14159265358979323846

constexpr int voiceSize = 1;
constexpr int bufferSize = 256;

class Param
{
public:
    alignas(16) double mPhase = 0.0;
    alignas(16) double mPhaseOptimized = 0.0;
    alignas(16) double mNoteFrequency = 10.0;
    alignas(16) double mHostPitch = 1.0;
    alignas(16) double mRadiansPerSample = 1.0;

    alignas(16) double b[voiceSize][bufferSize];
    alignas(16) double c[voiceSize][bufferSize];

    Param() { }

    inline void Process(int voiceIndex, int blockSize) {
        double *pB = b[voiceIndex];
        double *pC = c[voiceIndex];
        double phase = mPhase;
        double bp0 = mNoteFrequency * mHostPitch;

        for (int sampleIndex = 0; sampleIndex < blockSize; sampleIndex++) {
            // some other code (that will use phase, like sin(phase))

            phase += std::clamp(mRadiansPerSample * (bp0 * pB[sampleIndex] + pC[sampleIndex]), 0.0, PI);

            std::cout << sampleIndex << ": " << phase << std::endl;
        }

        mPhase = phase;
    }
    inline void ProcessOptimized(int voiceIndex, int blockSize) {
        double *pB = b[voiceIndex];
        double *pC = c[voiceIndex];
        double phase = mPhaseOptimized;
        double bp0 = mNoteFrequency * mHostPitch;

        __m128d v_boundLower = _mm_set1_pd(0.0);
        __m128d v_boundUpper = _mm_set1_pd(PI);
        __m128d v_radiansPerSampleBp0 = _mm_set1_pd(mRadiansPerSample * bp0);
        __m128d v_radiansPerSample = _mm_set1_pd(mRadiansPerSample);

        __m128d v_pB0 = _mm_load_pd(pB);
        v_pB0 = _mm_mul_pd(v_pB0, v_radiansPerSampleBp0);
        __m128d v_pC0 = _mm_load_pd(pC);
        v_pC0 = _mm_mul_pd(v_pC0, v_radiansPerSample);

        __m128d v_pB1 = _mm_setr_pd(0.0, pB[0]);
        v_pB1 = _mm_mul_pd(v_pB1, v_radiansPerSampleBp0);
        __m128d v_pC1 = _mm_setr_pd(0.0, pC[0]);
        v_pC1 = _mm_mul_pd(v_pC1, v_radiansPerSample);

        __m128d v_phase = _mm_set1_pd(phase);
        __m128d v_phaseAcc;

        for (int sampleIndex = 0; sampleIndex < blockSize; sampleIndex += 2, pB += 2, pC += 2) {
            // some other code (that will use phase, like sin(phase))

            v_phaseAcc = _mm_add_pd(v_pB0, v_pC0);
            v_phaseAcc = _mm_max_pd(v_phaseAcc, v_boundLower);
            v_phaseAcc = _mm_min_pd(v_phaseAcc, v_boundUpper);
            v_phaseAcc = _mm_add_pd(v_phaseAcc, v_pB1);
            v_phaseAcc = _mm_add_pd(v_phaseAcc, v_pC1);
            v_phase = _mm_add_pd(v_phase, v_phaseAcc);

            v_pB0 = _mm_load_pd(pB + 2);
            v_pB0 = _mm_mul_pd(v_pB0, v_radiansPerSampleBp0);
            v_pC0 = _mm_load_pd(pC + 2);
            v_pC0 = _mm_mul_pd(v_pC0, v_radiansPerSample);

            v_pB1 = _mm_load_pd(pB + 1);
            v_pB1 = _mm_mul_pd(v_pB1, v_radiansPerSampleBp0);
            v_pC1 = _mm_load_pd(pC + 1);
            v_pC1 = _mm_mul_pd(v_pC1, v_radiansPerSample);

            std::cout << sampleIndex << ": " << v_phase.m128d_f64[0] << std::endl;
            std::cout << sampleIndex + 1 << ": " << v_phase.m128d_f64[1] << std::endl;
        }

        mPhaseOptimized = v_phase.m128d_f64[blockSize % 2 == 0 ? 1 : 0];
    }
};

class MyPlugin
{
public: 
    Param mParam1;

    MyPlugin() {
        // fill b
        for (int voiceIndex = 0; voiceIndex < voiceSize; voiceIndex++) {
            for (int sampleIndex = 0; sampleIndex < bufferSize; sampleIndex++) {
                double value = (sampleIndex / ((double)bufferSize - 1));

                mParam1.b[voiceIndex][sampleIndex] = value;
            }
        }

        // fill c
        for (int voiceIndex = 0; voiceIndex < voiceSize; voiceIndex++) {
            for (int sampleIndex = 0; sampleIndex < bufferSize; sampleIndex++) {
                double value = 0.0;

                mParam1.c[voiceIndex][sampleIndex] = value;
            }
        }
    }
    ~MyPlugin() { }

    void Process(int blockSize) {
        for (int voiceIndex = 0; voiceIndex < voiceSize; voiceIndex++) {
            mParam1.Process(voiceIndex, blockSize);
        }
    }
    void ProcessOptimized(int blockSize) {
        for (int voiceIndex = 0; voiceIndex < voiceSize; voiceIndex++) {
            mParam1.ProcessOptimized(voiceIndex, blockSize);
        }
    }
};

int main() {
    MyPlugin myPlugin;

    long long numProcessing = 1;
    long long counterProcessing = 0;

    // I'll only process once block, just for analysis
    while (counterProcessing++ < numProcessing) {
        // variable blockSize (i.e. it can vary, being even or odd)
        int blockSize = 256;

        // process data
        myPlugin.Process(blockSize);
        std::cout << "#########" << std::endl;
        myPlugin.ProcessOptimized(blockSize);
    }
}

Answer 1

（更新：这个答案是在显示 v_phase 在循环 中使用的编辑之前写的。）

等一下，我以为你在每一步都需要 phase 的值。是的，循环中有一个 // some other code (that will use phase) 注释。

但看起来您只对最终值感兴趣。所以你可以自由地重新排序，因为每一步的夹紧都是独立的。

这只是一个缩减（如数组的总和），并进行了一些动态处理以生成缩减的输入。

您希望 v_phase 的 2 个元素是偶数/奇数元素的 2 个独立部分和。然后你在最后横向求和。（例如 _mm_unpackhi_pd(v_phase, v_phase) 将高半部分置于底部，或参见 Fastest way to do horizontal float vector sum on x86）。

然后可以选择对结果使用标量 fmod 以将范围缩小到 [0..2Pi) 范围内。（如果事实证明精度成为问题，则在求和期间偶尔缩小范围可以通过阻止值变得太大来帮助提高精度。）

如果情况并非如此，并且您在每个 i+=2 步骤确实需要一个 { phase[i+0], phase[i+1] } 向量，那么您的问题似乎与 prefix sum 有关.但是每个向量只有 2 个元素，对具有未对齐负载的元素冗余地执行所有操作可能是有意义的。

可能比我想象的要少节省，因为您需要分别限制每个步骤：在乘法之前执行 pB[i+0] + pB[i+1] 可能会导致不同的限制。

但是您显然已经删除了我们简化公式中的限制，因此您可以在应用 mul/add 公式之前添加元素。

或者一次完成两个步骤的 multiply/add 东西可能是一个胜利，然后将其洗牌以添加正确的东西。

这个 SSE2 转置有什么问题？

What's wrong in this SSE2 transposition?

c++

vectorization

intrinsics

sse2