为什么同一程序的 c 和 fortran 版本会产生不同的结果?

Why the c and fortran versions of this same program produce different results?

提问人:Antonio Serrano 提问时间:11/12/2020 更新时间:11/12/2020 访问量:241

问:

我用的是gcc(Ubuntu 9.3.0-17ubuntu1~20.04)9.3.0

c代码为:

// Compile with:
//  gcc -o little_c little.c
#include <stdio.h>  // printf

void main(void) {
    int n = 800;
    float a[n][n], b[n][n], c[n][n];
    int i, j, k;

    for (i = 0; i < n; i++) {
        for (j = 0; j < n; j++) {
            a[i][j] = (float) (i+j);
            b[i][j] = (float) (i-j);
        }
    }

    for (i = 0; i < n; i++) {
        for (j = 0; j < n; j++) {
            float t = (float) 0.0;
            for (k = 0; k < n; k++)
                t += a[i][k] * a[i][k] + b[k][j] * b[k][j];
                //t += a[i][k] + b[k][j]; // If I comment the above line and uncomment this, the c and fortran reults are the same
            c[i][j] = t;
        }
    }
    printf("%f", c[n-1][n-1]); // prints the very last element

}

Fortran 代码:

! Compile with:
!  gfortran -o little_fort little.f90 

program little

implicit none

integer, parameter  :: n = 800
real                :: a(n,n), b(n,n), c(n,n)
real                :: t
integer             :: i, j, k  ! Counters

do i = 1, n
    do j = 1, n
        a(i,j) = real(i-1+j-1)      ! Minus one, for it to be like the c version
        b(i,j) = real(i-1-(j-1))    ! because in c, the index goes from 0 to n-1
    end do
end do

do i = 1, n
    do j = 1, n
        t = 0.0
        do k = 1, n
            t = t + a(i,k) * a(i,k) + b(k,j) * b(k,j)
            !t = t + a(i,k) + b(k,j) ! If I comment the above line and uncomment this, the c and fortran reults are the same
        end do
        c(i,j) = t
    end do
end do
    
write (*,"(F20.4)") c(n,n)  ! This is the same as c[n-1][n-1] in c


end program little

c 程序打印:1362136192.000000

Fortran 程序打印:1362137216.0000

如果我不将每个元素本身相乘,正如我在代码的注释中所述,我会得到两个版本的程序相同的值:

c prigram: 639200.000000

Fortran 计划:639200.0000

为什么当我使用乘法时,c 和 Fortran 代码会产生不同的结果?它是否必须与实数的不同实现有关?

C Fortran 精确 Gfortran

评论

1赞 Weather Vane 11/12/2020
结果在第 7 位有效数字处有所不同,这与 .我建议你使用而不是到处使用,除非有很好的理由不这样做。请参阅为什么浮点数不准确?floatdoublefloat
1赞 Max 11/12/2020
我不确定,但是否有可能一个使用 SSE,一个使用 x87 数学,它们具有不同的精度?你能拆开内部循环,看看它们是否使用不同的FP子系统吗?此外,float 只有 32 位精度,我不确定 real 默认做什么。
1赞 John Alexiou 11/12/2020
只是把这个放在这里:每个计算机科学家都应该知道的关于浮点运算的知识
0赞 tim18 11/13/2020
根据 Max 所说的,如果您应该设置编译选项以允许 sse 优化,由于您使用了 stride 1 内部循环,Fortran 更有可能加速并获得更准确的结果。为了使 Fortraersrersn 计算更像 C 而部分使用括号的另一种方法是使用 SUM 内部函数,因为不再有编译器局限于 Fortran 77。

答:

6赞 rtoijala 11/12/2020 #1

这种差异是由于计算顺序与浮点类型的有限精度相结合所致。

如果将 Fortran 版本更改为

t = t + (a(i,k) * a(i,k) + b(k,j) * b(k,j))

即在术语 和 周围添加括号,您将获得两种语言的相同结果。由于使用了赋值运算符,C 版本已经使用了此计算顺序。ab+=

如注释中所述,这是在可用精度限制下的预期行为。

评论

1赞 Eric Postpischil 11/12/2020
或者,在我的 C 实现中,更改 C 版本中的行以获得与 FORTRAN 版本相同的结果(Apple Clang 11,带有 )。我怀疑可能有多种评估方法可以得到相同的结果。t += b[k][j] * b[k][j] + a[i][k] * a[i][k];-O3
1赞 John Alexiou 11/12/2020
请注意,两个不同量级的加法和最重要的减法会引入明显的舍入误差。
0赞 Jim Rogers 11/12/2020 #2

当我编写该程序的 Ada 版本时,我发现我必须将十进制精度降低到 6 位小数才能获得 Fortran 答案。

Ada 版本是:

与Ada.Text_IO;使用Ada.Text_Io;

procedure Main is
   type Real is digits 6;
   package Real_IO is new Ada.Text_IO.Float_IO(Real);
   use Real_IO;
   
   subtype Index is integer range 1..800;
   type Matrix is array(Index, Index) of Real;
   
   A : Matrix;
   B : Matrix;
   C : Matrix;
   T : Real := 0.0;
begin
   for I in Index loop
      for J in Index loop
         A(I,J) := Real(I - 1 + J - 1);
         B(I,J) := Real(I - 1 - (J - 1));
      end loop;
   end loop;
   
   for I in Index loop
      for J in Index loop
         T := 0.0;
         for K in Index loop
            T := T + A(I,K) * A(I,K) + B(K,J) *B(K,J);
         end loop;
         C(I,J) := T;
      end loop;
   end loop;
   
   Put(Item => C(Index'Last, Index'Last), Exp => 0, Aft => 4);
   New_Line;
end Main;

线定义类型 Real 定义浮点类型的精度:

type Real is digits 6;

使用六位数精度产生的值为

1362137216.0000

使用更高精度的浮点类型导致该值

1362135200.0000

评论

0赞 Vladimir F Героям слава 11/12/2020
那么,这个问题的答案是什么?
0赞 Jim Rogers 11/12/2020
答案是,差异可能由多种因素引起,包括由于不同的计算顺序而导致的不同舍入误差,以及由于精度降低而导致的舍入误差,浮点数。
1赞 Vladimir F Героям слава 11/12/2020
但是问题中的 Fortran 和 C 都使用单精度。证明它在不同语言中的双精度会有所不同有什么意义?很明显,在任何语言的双精度下,它都会有所不同。
0赞 Jim Rogers 11/12/2020
我使用了 GNAT Ada 编译器,它是 GNU 工具链的一部分。它应该产生与GNU C++编译器相同的结果。两种语言共享相同的编译器后端。