Разные результаты с и без SSE (умножение чисел с плавающей запятой)

У меня есть две функции умножения двумерных массивов. Один из них с SSE. Еще одна функция без какой-либо оптимизации. Обе функции работают хорошо. Но результаты немного отличаются. Например, 20,333334 и 20,333332.

Можете ли вы объяснить, почему результаты отличаются? И что я могу сделать с функциями, чтобы получить тот же результат?

функция с SSE

float** sse_multiplication(float** array1, float** array2, float** arraycheck)
{
    int i, j, k;
    float *ms1, *ms2, result;
    float *end_loop;

    for( i = 0; i < rows1; i++)
    {
        for( j = 0; j < columns2; j++)
        {
            result = 0;
            ms1 = array1[i];
            ms2 = array2[j];
            end_loop = &array1[i][columns1];

            __asm{
                     mov rax, ms1
                     mov rbx, ms2
                     mov rdx, end_loop
                     xorps xmm2, xmm2

                loop:
                     movups xmm0, [rax]
                     movups xmm1, [rbx]
                     movups xmm3, [rax+16]
                     movups xmm4, [rbx+16]

                     mulps xmm0, xmm1
                     mulps xmm3, xmm4

                     addps xmm2, xmm0

                     add rax, 32
                     add rbx, 32

                     cmp rdx, rax
                     jne loop

                     haddps xmm2, xmm2
                     haddps xmm2, xmm2

                     movups result, xmm2
               }

             arraycheck[i][j] = result;
        }
    }
    return arraycheck;
}

функция без какой-либо оптимизации

float** multiplication(float** array1, float** array2, float** arraycheck)
{
    for (int i = 0; i < rows1; i++)
        for (int j = 0; j < columns2; j++)
            for (int k = 0; k < rows1; k++)
                arraycheck[i][j] += array1[i][k] * array2[k][j];

    return arraycheck;
}
Источник

2 ответа

Решение

Добавление FP не совсем ассоциативно, поэтому другой порядок операций приводит к немного отличающимся ошибкам округления.

Ваш С суммирует элементы по порядку. (Если вы не используете -ffast-math чтобы позволить компилятору сделать то же самое предположение, что операции FP достаточно близки к ассоциативным).

Ваша асма суммирует каждый 4-й элемент с 4 различными смещениями, а затем суммирует их по горизонтали. Сумма в каждом элементе вектора округляется по-разному в каждой точке.

Ваша векторизованная версия, похоже, не соответствует версии C. Индексация выглядит иначе. AFAICT, единственный вменяемый способ векторизации arraycheck[i][j] += array1[i][k] * array2[k][j]; кончено j, Зацикливание k потребует больших нагрузок от array2 и зацикливание i потребует больших нагрузок от array1,

Я что-то упускаю из-за вашего асма? Он загружает смежные значения из обоих массивов. Это также выбрасывает mulps результат в xmm3 каждая итерация loop , так что я думаю, что это просто глючит.

С циклом j во внутреннем векторном цикле не меняется array1[i][k] Просто передайте-загрузите его один раз за пределы цикла (_mm256_set1_ps).

Тем не менее, это означает, что делать чтение-изменение-запись arraycheck[i][j] для каждого другого j значение. т.е. ac[i][j + 0..3] = fma(a1[i][k], a2[k][j + 0..3], ac[i][j + 0..3]), Чтобы избежать этого, вам нужно сначала перенести один из массивов. (Но это O(N^2) для матрицы NxN, которая все же дешевле, чем умножение).

Этот способ не использует горизонтальные суммы, но посмотрите эту ссылку, если вы хотите лучший код для этого.

Он выполняет операции в том же порядке, что и скаляр C, поэтому результаты должны точно совпадать.

Также обратите внимание, что вам нужно использовать более одного аккумулятора для насыщения исполнительных блоков ЦП. Я бы предложил 8, чтобы насытить задержку Skylake 4c, по одному на пропускную способность 0.5c addps, Haswell имеет задержку 3c, по одному на 1c addps, но Skylake отбросил отдельный модуль добавления FP и делает это в модуле FMA. (См. Вики-тег x86, особенно руководства Агнера Фога)

На самом деле, так как мое предлагаемое изменение вообще не использует один аккумулятор, каждая итерация цикла обращается к независимой памяти. Вам нужно будет немного развернуть цикл, чтобы насытить исполнительные блоки FP двумя загрузками и сохранить в цикле (даже если вам нужны только два указателя, поскольку хранилище вернулось в то же место, что и одна из загрузок). Но в любом случае, если ваши данные помещаются в кэш-память L1, выполнение в неупорядоченном порядке должно обеспечивать достаточное снабжение исполнительных блоков работой из отдельных итераций.

Если вы действительно заботитесь о производительности, вы сделаете версию FMA и, возможно, версию AVX-без-FMA для Sandybridge. Вы могли бы делать два 256-битных FMA на такт вместо одного 128-битного и мульти на такт. (И, конечно, вы этого даже не получаете, потому что у вас узкое место в задержке, если цикл не достаточно короткий для того, чтобы окно вышло из строя, чтобы увидеть независимые инструкции из следующей итерации).

Вам понадобится "циклическое разбиение по циклам", или "блокировка кеша", чтобы это не помогало при больших размерах проблем. Это матрица умножения, верно? Для этого есть очень хорошие библиотеки, которые настроены на размеры кеша и позволят избежать такой простой попытки, как эта. например, ATLAS был хорош в прошлый раз, когда я проверял, но это было несколько лет назад.

Используйте встроенные функции, если только вы не напишите всю функцию в asm. Компиляторы "понимают", что они делают, поэтому могут делать приятные оптимизации, например, развертывание цикла при необходимости.

Источник

В соответствии со стандартными форматами IEEE32-разрядные числа с плавающей запятой могут гарантировать точность только 6-7 цифр. Ваша ошибка настолько незначительна, что на механизм компилятора невозможно сделать правдоподобное утверждение. Если вы хотите добиться большей точности, было бы разумно выбрать либо 64-битный двойной (гарантия точности 15 цифр), либо реализовать свой собственный класс BigDecimal, такой как java do.

Источник
Другие вопросы по тегам