使用和不使用 SSE 的不同结果（浮点数组乘法）

Question

我有两个二维数组乘法函数。其中之一是 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;
}

Answer 1

根据IEEE standard Formats, 32-bit float can only guanartee 6-7 digits accuracy. Your error is so marginal that no plausible claim can be made on compiler's mechanism. If you want to achieve better precision, it would be wise to choose either 64-bit double(guarentees 15 digits accuracy) or implement your own BigDecimal class 喜欢java 做。

Answer 2

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 的某些信息？它从两个数组加载连续的值。 它也在 xmm3 的 loop 的每次迭代中丢弃了 mulps 结果，所以我认为它只是有问题 .

由于在内部向量循环中循环 j 不会改变 array1[i][k]，只需在循环外广播加载一次 (_mm256_set1_ps)。

然而，这意味着对每个不同的 j 值执行 arraycheck[i][j] 的读取-修改-写入。即 ac[i][j + 0..3] = fma(a1[i][k], a2[k][j + 0..3], ac[i][j + 0..3])。为避免这种情况，您必须先转置其中一个数组。（但这对于 NxN 矩阵来说是 O(N^2)，这仍然比乘法便宜）。

这种方式不使用 horizontal sums，但如果您需要更好的代码，请参阅 link。

它按照与标量 C 相同的顺序执行操作，因此结果应该完全匹配。

另请注意，您需要使用多个累加器来使 CPU 的执行单元饱和。我建议 8，以饱和 Skylake 的 4c 延迟，每 0.5c 吞吐量一个 addps。 Haswell 有 3c 延迟，每 1c addps，但 Skylake 放弃了单独的 FP 添加单元并在 FMA 单元中进行。（见x86 tag wiki, esp. Agner Fog's guides）

实际上，由于我建议的更改根本不使用单个累加器，因此每个循环迭代都访问独立的内存。您将需要一些循环展开来使 FP 执行单元充满两个加载并存储在循环中（即使您只需要两个指针，因为存储返回到与其中一个加载相同的位置）。但是无论如何，如果您的数据适合 L1 缓存，乱序执行应该可以很好地为执行单元提供来自单独迭代的工作。

如果你真的关心性能，你会制作一个 FMA 版本，并且可能为 Sandybridge 制作一个 AVX-without-FMA 版本。您可以每个时钟执行两个 256b FMA，而不是每个时钟执行一个 128b add 和 mul。（当然，你甚至没有得到它，因为你会在延迟上遇到瓶颈，除非循环足够短，以至于乱序 window 可以看到下一次迭代的独立指令）。

您将需要“循环平铺”，又名“缓存阻塞”，以使其不会对大问题造成影响。这是一个矩阵乘法，对吧？有非常好的库，它们针对缓存大小进行了调整，并且可以通过像这样的简单尝试来解决问题。例如ATLAS 上次我检查时还不错，但那是几年前的事了。

使用内在函数，除非你用 asm 编写整个函数。编译器“理解”它们所做的事情，因此可以进行很好的优化，例如在适当的时候展开循环。

使用和不使用 SSE 的不同结果（浮点数组乘法）

different results with and without SSE ( float arrays multiplication)

c++

arrays

floating-point

sse

FP 加法不是完全关联的，因此不同的运算顺序会产生略微不同的舍入误差。