Ответ 1
Добавление FP не является вполне ассоциативным, поэтому другой порядок операций вызывает несколько разные ошибки округления.
Ваши C суммируют элементы по порядку. (Если вы не используете -ffast-math, чтобы позволить компилятору сделать то же самое предположение, что вы делали это, операции FP достаточно близки к ассоциативному).
Ваш asm суммирует каждый 4-й элемент с 4 различными смещениями, затем горизонтально суммирует их. Сумма в каждом векторном элементе округляется по-разному в каждой точке.
Ваша векторная версия, похоже, не соответствует версии C. Индексирование выглядит по-другому. AFAICT, единственный разумный способ векторизации arraycheck[i][j] += array1[i][k] * array2[k][j]; - над j. Для циклического перехода через k потребуются перемещенные нагрузки от array2, а для циклического перехода через i потребуются перемещенные нагрузки от array1.
Я что-то пропустил в вашем asm? Он загружает смежные значения из обоих массивов. Он также отбрасывает результат mulps в xmm3 каждую итерацию loop, поэтому я думаю, что это просто глючит.
Так как циклирование по j во внутреннем векторном цикле не меняет array1[i][k], просто передайте его один раз за пределы цикла (_mm256_set1_ps).
Однако это означает выполнение read-modify-write 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, поэтому результаты должны точно соответствовать.
Также обратите внимание, что вам необходимо использовать более одного аккумулятора для насыщения исполнительных блоков CPU. Я бы предложил 8, чтобы насытить Skylake 4c латентность, по одной на 0,5c пропускную способность addps. Haswell имеет 3c латентность, по одному на 1c addps, но Skylake сбросил отдельный блок добавления FP и делает это в блоке FMA. (Смотрите x86 tag wiki, esp. Руководство Agner Fog)
На самом деле, поскольку в предложенном мной изменении вообще не используется один аккумулятор, каждая итерация цикла обращается к независимой памяти. Вам нужно немного развернуть цикл, чтобы насытить исполняемые модули FP двумя нагрузками и сохранить в цикле (хотя вам нужны только два указателя, так как хранилище возвращается в то же место, что и одна из нагрузок). Но в любом случае, если ваши данные вписываются в кеш-память L1, выполнение вне очереди должно обеспечивать, чтобы исполняемые модули были хорошо снабжены работой из отдельных итераций.
Если вы действительно заботитесь о производительности, вы сделаете версию FMA и, возможно, версию AVX-без FMA для Sandybridge. Вы можете делать два 256-бит FMA за такт вместо одного 128b add и mul за такт. (И, конечно же, вы даже этого не понимаете, потому что вы отстаете от латентности, если цикл не достаточно короткий, чтобы окно вне порядка отображало независимые инструкции со следующей итерации).
Вам понадобится "петлевая черепица", иначе "блокировка кеша", чтобы это не сосать для больших размеров проблем. Это матрица умножить, правильно? Для этого есть очень хорошие библиотеки, которые настроены на размеры кеша и будут избивать штаны простой попыткой. например ATLAS был хорошим в прошлый раз, когда я проверил, но это было несколько лет назад.
Использовать intrinsics, если вы не напишете всю функцию в asm. Компиляторы "понимают", что они делают, поэтому могут делать приятные оптимизации, например, при необходимости разворачивать цикл.