MATMUL 结果与双精度显式计算不相等？

Question

对不起，一个看似愚蠢的问题。在用内在函数替换矩阵上的 for 循环操作时，我正在测试计算效率。当我检查两种方法的矩阵乘积结果时，我很困惑这两种输出不一样。这是我使用的简化代码

program matmultest
    integer,parameter::nx=64,ny=32,nz=16
    real*8::mat1(nx,ny),mat2(ny,nz)
    real*8::result1(nx,nz),result2(nx,nz),diff(nx,nz)
    real*8::localsum
    integer::i,j,m
    do i=1,ny
        do j=1,nx
            mat1(j,i)=dble(j)/7d0+2.65d0*dble(i)
        enddo
    enddo
    do i=1,nz
        do j=1,ny
            mat2(j,i)=5d0*dble(j)-dble(i)*0.45d0
        enddo
    enddo
    do j=1,nz
        do i=1,nx
            localsum=0d0
            do m=1,ny
                localsum=localsum+mat1(i,m)*mat2(m,j)
            enddo
            result1(i,j)=localsum
        enddo
    enddo
    result2=matmul(mat1,mat2)
    diff=result2-result1
    print*,sum(abs(diff)),maxval(diff)
end program matmultest

结果给出

   1.6705598682165146E-008   5.8207660913467407E-011

real8当我integer稍后测试时，差异不为零，但为零。我想知道是因为我的代码在某处出错还是MATMUL()单精度的数值精度？

我使用的编译器是 GNU Fortran (Ubuntu 9.3.0-17ubuntu1~20.04) 9.3.0

谢谢！

Answer 1

francescalus 解释说，操作的重新排序会导致这些差异。让我们试着找出它实际上是如何发生的。

关于矩阵乘积的几句话

考虑矩阵 A(n,p)、B(p,q)、C(n,q) 和 C = A*B。

天真的方法（您使用的一种变体）涉及以下嵌套循环：

c = 0
do i = 1, n
    do j = 1, p
        do k = 1, q
            c(i, j) = c(i, j) + a(i, k) * b(k, j)
        end do
    end do
end do

这些循环可以按 6 个顺序中的任何一个执行，具体取决于您在每个级别选择的变量。在上面的例子中，循环被命名为“ijk”，其他变体“ikj”、“jik”等都是正确的。

由于内存缓存，存在速度差异：当内循环运行在连续的内存元素上时，循环会更快。那是 jki 或 kji 案例。

实际上，由于 Fortran 矩阵按列主序存储，如果最内层循环在 i 上运行，则在指令 c(i, j) = c(i, j) + a(i, k) * c(k, j) 中，值 c(k, j) 是常数，运算等价于 v(i) = v(i) + x * u(i)，其中向量 v 和 u 的元素是连续的。

但是，关于操作顺序，应该没有区别：您可以自己检查 C 的所有元素是否以相同的顺序计算。至少在“更高级别”：编译器可能会以不同的方式优化事物，这才是真正有趣的地方。

MATMUL呢？我相信它通常是一个简单的矩阵产品，基于上面的嵌套循环，比如 jki 循环。

还有其他方法可以乘以矩阵，包括使用Strassen 算法来提高算法复杂度或使用阻塞（即子矩阵的计算乘积）来改进缓存使用。其他可能改变结果的方法是OpenMP（即多线程），或使用FMA 指令。但在这里我们不打算深入研究这些方法。它实际上仅与嵌套循环有关。如果你有兴趣，网上有很多资源，看看这个。

关于优化的几句话

先说三点：

• 在没有SIMD 指令的处理器上，您会得到与 MATMUL 相同的结果（即最后会打印零）。
• 如果您已经实现了上述循环，您也会得到相同的结果。您的代码中存在微小但显着的差异。
• 如果您将循环实现为子例程，您也会得到相同的结果。在这里，我怀疑编译器优化器正在进行一些重新排序，因为我无法使用子例程重现您的“累加器”代码，至少使用 Intel Fortran 是这样。

这是您的实现：

do i = 1, n
    do j = 1, p
        s = 0
        do k = 1, q
            s = s + a(i, k) * b(k, j)
        end do
        c(i, j) = s
    end do
end do

当然也是对的。在这里，您使用的是累加器，并且在最内层循环的末尾，累加器的值写入矩阵 C 中。

优化通常主要与最内层循环相关。为了我们的目的，如果我们去掉所有其他细节，最内层循环中的两个“基本”指令是相关的：

• v(i) = v(i) + x*u(i)（jki 循环）
• s = s + x(k)*y(k)（累加器循环，其中 y 在内存中是连续的，但不是 x）

第一个通常称为“daxpy”（来自 BLAS 例程的名称），对于“AX Plus Y”，“D”表示双精度。第二个只是一个累加器。

在旧的顺序处理器上，没有太多可优化的工作。在具有 SIMD 的现代处理器上，寄存器可以保存多个值，并且可以同时并行地对所有值进行计算。例如，在 x86 上，一个 XMM 寄存器（来自 SSE 指令集）可以保存两个双精度浮点数。一个 YMM 寄存器（来自 AVX2）可以容纳四个数字，一个 ZMM 寄存器（AVX512，在 Xeon 上找到）可以容纳八个数字。

例如，在 YMM 上，最内层的循环将被“展开”以一次处理四个向量元素（如果使用多个寄存器甚至更多）。

以下是基本循环块大致执行的操作：

daxpy 案例：

• 将 u 中的 4 个数字读入寄存器 YMM1
• 将 v 中的 4 个数字读入寄存器 YMM2
• x 是常数并保存在另一个寄存器中
• x与YMM1平行相乘，与YMM2平行相加，结果放入YMM2
• 将结果写回到 v 的对应元素

如果元素在内存中是连续的，则读/写部分会更快，但如果不是，则仍然值得并行执行此操作。

请注意，在这里，我们没有更改高级 Fortran 循环添加的执行顺序。

蓄电池箱

为了使并行性有用，将有一个技巧：在 YMM 寄存器中并行累加四个值，然后将四个累加值相加。

因此，基本循环块是这样做的：

• 累加器保存在YMM3（四个数字）
• 将 X 中的 4 个数字读入寄存器 YMM1
• 将 Y 中的 4 个数字读入寄存器 YMM2
• YMM1 与 YMM2 并行相乘，与 YMM3 并行相加

在最内层循环的末尾，将累加器的四个分量相加，并将其写回作为矩阵元素。

就像我们计算过：

• s1 = x(1)*y(1) + x(5)*y(5) + ... + x(29)*y(29)
• s2 = x(2)*y(2) + x(6)*y(6) + ... + x(30)*y(30)
• s3 = x(3)*y(3) + x(7)*y(7) + ... + x(31)*y(31)
• s4 = x(4)*y(4) + x(8)*y(8) + ... + x(32)*y(32)

然后写入的矩阵元素为 c(i,j) = s1+s2+s3+s4。

此处添加顺序已更改！然后，由于顺序不同，结果很可能不同。