Fortran 中"end program"语句中的分段错误

Question

我有一些Fortran代码，它由一个子例程和对它的调用组成。它应该使用定义的窗口大小来计算矩阵中元素的平均值。例如，使用winsize＝2在(10，10)数组上调用子例程将返回(5，5)数组。

代码如图所示：

SUBROUTINE avgwin(ts, sizelat,sizelon,winsize,size2,size3,ts_new)
implicit none
double precision, dimension(10,sizelat,sizelon) :: ts
double precision, dimension(winsize,winsize) :: store
double precision, dimension(10,size2,size3) :: ts_new
double precision :: par,ave
integer :: sizelat, sizelon,i,j,k,winsize,size2,size3
integer :: A, B,p,m,numb
A=0
B=0
par = 11 !Hypothetical value to be excluded
do i=1,10 !Looping through time
do j=1,sizelat !Looping through latitude
if ((j+winsize) > sizelat) then !Checks if it will exceed bounds
exit !If it'll exceed, discard remaining cells
end if
do k=1,sizelon !Looping through longitude
if ((k+winsize)>sizelon) then
exit
end if
store = ts(i,j:j+winsize,k:k+winsize) !Gets the values for that window
where (store == par) store = -99 !Replaces masked with -99
ave = 0
numb = 0 !Variable to count 
do p=1,winsize
do m=1,winsize
if (store(p,m)==-99) then !Evaluates if it's masked, i.e., =-99
ave = ave
else
ave = ave + store(p,m) !Sum of existent values
numb = numb +1 !Updates counting variable
end if
end do
end do
ave = ave/numb !Calculates the mean
ts_new(i,A,B) = ave 
B=B+1
end do
B=0
A=A+1
end do
A=0
B=0
end do
END SUBROUTINE
program testefor
implicit none
double precision, dimension(10,10,10) :: teste
double precision, dimension(10,5,5) :: oi
integer :: i,j,k
do i=1,10
do j=1,10
do k=1,10
teste(i,j,k)=i
end do
end do
end do
CALL avgwin(teste,10,10,2,5,5,oi)
print*, oi(1,5,5)
end program testefor

然而，当我运行它时，我会遇到分段故障。我试着用GDB调试它，令我惊讶的是，它在退出程序时返回了正确的结果，但出现了segfault。我从gdb得到的可以在下面找到：

Breakpoint 1, testefor () at testefor.f90:56
56      do i=1,10
(gdb) cont
Continuing.
1.0000000000000000     
Program received signal SIGSEGV, Segmentation fault.
0x0000000000400f73 in testefor () at testefor.f90:67
67  end program testefor

因此，程序返回正确的(1,5,5)元素=1.0，但在其他地方出错。

有人能帮我找出问题吗？感谢

Answer 1

有一件事可以帮助您在调试模式下编译程序，并激活许多调试开关。

在我的案例中，我使用gfortran进行编译，并使用：

1. -g包含调试符号
2. -fbacktrace以获得更好的堆栈跟踪
3. -Wall启用所有编译器警告
4. -fcheck=all来启用运行时检查(使程序速度变慢，但在调试过程中非常有用)

这些运行时检查立即发现错误：

At line 21 of file teste.F
Fortran runtime error: Array bound mismatch for dimension 1 of array 'store' (2/3)

store的大小是(winsize, winsize)，但您将一个数组复制到大小为(windsize+1, winsize+1)的数组中。如果在Fortran中对数组进行切片，则它包括起始索引和结束索引：(1:10)是从1到10。如果您有(1:1+10)，它从1变为11，这意味着它的大小是11。

如果没有运行时数组边界检查(在本例中由-fcheck=all激活)，这真的很难调试，并且可能导致各种意外行为。

如果您使用的编译器与gfortran不同，您需要了解如何在编译器上打开此类测试，因为开关不是标准化的。

Answer 2

这并不是你问题的答案(我已经回答过了)，而是一些对你有帮助的提示。

1. 如果你在多维数组上循环，你应该总是在第一个索引上有最内部的循环，以此类推，所以

   ! Inefficient way to do it:
   do i = 1, 10
   do j = 1, 10
   do k = 1, 10
   a(i, j, k) = i*j+k
   end do
   end do
   end do
   !Efficient way to do it:
   do k = 1, 10                   !  <-+
   do j = 1, 10               !    | swapped
   do i = 1, 10           !  <-+
   a(i, j, k) = i*j+k
   end do
   end do
   end do
   

   Fortran存储多维数组的方式是第一个索引变化最快。因此，高效的方法读取连续的元素，而低效的方法需要在内存中跳来跳去。

2. 代替这个复杂的循环，您可以同样轻松地使用SUM和COUNT:

   • COUNT计算数组中.TRUE.的实例，因此numb的值可以像numb=COUNT(store/=-99)一样容易地计算出来
   • SUM有一个可选参数MASK，您可以使用它来省略某些值。要汇总所有非-99的值，可以使用：s = SUM(store, MASK=(store /= -99))

   所以你可以用替换大约15行代码

   ave = sum(store, MASK=(store/=par)) / count(store/=par)
   

3. store和par都是double precision，比较浮点变量的绝对相等性是一个棘手的