双精度矩阵乘法中MATMUL结果与手动循环计算结果不一致问题

核心结论 #

你的代码不存在逻辑错误，MATMUL内置函数也并非仅支持单精度精度，你观测到的微小差值是双精度浮点数运算的正常现象，不属于错误范畴。

原因说明 #

• 浮点数加法不满足严格的数学结合律，不同的累加顺序会导致最终计算结果出现微小的数值差异。你手写的三重循环采用的是m从1到ny的逐次串行累加逻辑，而GNU Fortran的MATMUL为了最大化计算效率，内部会采用分块运算、SIMD向量指令并行累加、缓存友好的运算顺序调整等优化策略，和你手写循环的累加顺序不一致，因此会产生可观测的微小差值。
• 你输出的最大差值仅为5.8207660913467407E-011，远低于双精度浮点数的相对精度上限（约1e-15）乘以你的矩阵元素量级，属于完全可接受的数值误差范围，不会影响后续计算的有效性。
• 整数运算不存在精度损失，所有整数加减乘操作都是精确的，因此你用整数类型测试时两种方法的结果完全一致，差值为0，符合预期。

可选优化建议 #

如果需要进一步缩小两种方法的结果差异，可以在编译时添加浮点精度严格限制参数，比如给gfortran添加-ffloat-store或-fstrict-float编译选项，该选项会限制编译器对浮点运算顺序的自动优化，能让两种方法的计算结果进一步对齐，但代价是MATMUL的计算效率会有所下降。

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