MATLAB转Fortran的高效乘法算法结果不一致问题排查

借助@veryreverie的思路，我修改实现了块对角矩阵与向量的乘法（矩阵为特定块对角结构）。但出现异常：该算法在MATLAB中运行无问题，直接乘法与C'*x - bn的结果差异因机器加法非关联性控制在1e-15范围内；但迁移至Fortran后结果差异显著，甚至出现NaN值。该问题仅在m与n不相等时出现。

Fortran 实现代码 #

module RHS
implicit none
contains
subroutine RHS_eval(m,n,C,MM,q,LAP,BCterm,ADV,bn)

real*8, intent (in) :: x(:),C(:,:)
real*8, intent (out) :: bn(:)
integer*8 :: row, col
integer*8, intent (in) :: m,n

integer*8 :: offset, blkSz, blkShft, i
integer*8 :: xstart, ystart

bn = 0

bn(1:(m-1)*(n-1)) = bn(1:(m-1)*(n-1)) + x(1:(m-1)*(n-1));
bn(1:(m-1)*(n-1)) = bn(1:(m-1)*(n-1)) - x(m:(m-1)*n);

blkSz = m;
offset = (m-1)*(n);
blkShft = m-1;

do i = 1,(n-1)
    xstart = offset + blkSz*(i-1)+1;
    ystart = (i-1)*(m-1)+1;
    bn(ystart:ystart+blkShft-1) = bn(ystart:ystart+blkShft-1) - x(xstart:xstart+blkShft-1);
    bn(ystart:ystart+blkShft-1) = bn(ystart:ystart+blkShft-1) + x(xstart+1:xstart+blkShft);
enddo

end subroutine RHS_eval

end module RHS

MATLAB 实现代码 #

n = 15; m = 10; 
x = (10*rand((m-1)*n+(n-1)*m,1));
bn = zeros((n-1)*(m-1),1);

bn(1:(m-1)*(n-1)) = bn(1:(m-1)*(n-1)) + x(1:(m-1)*(n-1));
bn(1:(m-1)*(n-1)) = bn(1:(m-1)*(n-1)) - x(m:(m-1)*n);

blkSz = m;
offset = (m-1)*(n);
blkShft = m-1;

for i = 1:(n-1)
    xstart = offset + blkSz*(i-1)+1;
    ystart = (i-1)*(m-1)+1;
    bn(ystart:ystart+blkShft-1) = bn(ystart:ystart+blkShft-1) - x(xstart:xstart+blkShft-1);
    bn(ystart:ystart+blkShft-1) = bn(ystart:ystart+blkShft-1) + x(xstart+1:xstart+blkShft);
end

迭代结果情况 #

我们对矩阵C取转置后与右侧向量相乘，因已知矩阵结构无需存储矩阵。经过小幅修正后，早期迭代结果符合预期——自定义乘法算法与matmul(transpose(C),x)的差值（即maxval(abs(bn - matmul(transpose(C),(x))))）很小；但约第60步开始，结果偏差急剧增大，甚至出现NaN值。

内容的提问来源于stack exchange，提问作者2Napasa