Что не так в этом транспонировании SSE2? - PullRequest
Новая
       77

Что не так в этом транспонировании SSE2?

0 голосов
/ 03 января 2019

Я пытаюсь преобразовать этот код: 

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, пытаясь ускорить весь блок (который часто вызывается).Я использую 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];

Но, к сожалению, после суммы "шагов", результаты становятся очень разными для каждого значения фазы,Пытался отладить, но я не смог найти причину проблемы.

Кроме того, на самом деле это не так "быстро", как старая версия.

Вы можете распознатьпроблема?И как вы будете ускорять код?

Вот весь код, если вы хотите проверить два разных выхода: 

#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);
    }
}

1 Ответ

0 голосов
/ 03 января 2019

(обновление: этот ответ был написан до правок, показывающих, что v_phase используется внутри цикла .)

Подождите, подумал я в вашем предыдущем вопросе вам нужно было значение phase на каждом шаге.Да, внутри цикла был комментарий // some other code (that will use phase).

Но похоже, что вас интересует только конечное значение.Так что вы можете изменить порядок вещей, потому что зажим для каждого шага независим.

Это всего лишь сокращение (например, сумма массива) с некоторой обработкой на лету для генерации входных данных длясокращение.

Вы хотите, чтобы 2 элемента v_phase были 2 независимыми частичными суммами для четных / нечетных элементов.Тогда вы в конце суммируете его по горизонтали.(например, _mm_unpackhi_pd(v_phase, v_phase), чтобы привести верхнюю половину к основанию, или см. Самый быстрый способ сделать горизонтальную векторную сумму с плавающей запятой на x86 ).

Затем опционально используйте скаляр fmod для результатачтобы уменьшить диапазон до [0..2Pi).(Случайное уменьшение диапазона во время суммирования может помочь точности, не позволяя значению стать таким большим, если окажется, что точность становится проблемой.)

Если это не так, и выесли для каждого шага i+=2 нужен вектор { phase[i+0], phase[i+1] }, то ваша проблема связана с префиксной суммой .Но только с двумя элементами на вектор, просто избыточно делать все для элементов с невыровненными нагрузками, вероятно, имеет смысл.

Возможно, будет меньше экономии, чем я думал, так как вам нужно фиксировать каждый шаг отдельно: выполнение pB[i+0] + pB[i+1] перед умножениемможет привести к различным зажимам.

Но вы, очевидно, удалили зажим в нашей упрощенной формуле, поэтому вы можете потенциально добавлять элементы перед применением формулы mul / add.

Или, возможно, это победачтобы сделать умножение / сложение для двух шагов одновременно, затем перемешать, чтобы получить нужные вещи.

...