Fortran 90程序中openmp的分段错误

Question

我正在编写一个小程序来练习Fortran 90和openmp。以下示例代码在我使用gfortran -fopenmp heat.f90 -o heat进行编译时返回分段错误，但是在编译后的[[without -fopenmp选项中运行得很好。有人可以帮助我了解为什么会这样吗？

我在shell中设置export OMP_NUM_THREADS=4（在Windows 10上为Ubuntu上的Bash）

program heat

    use omp_lib
    implicit none

    integer, parameter :: nx = 1000
    integer, parameter :: ny = 800
    real, parameter :: H = 1.
    real, parameter :: alpha = 1.e-4
    real, parameter :: dt = 0.001
    real, parameter :: physical_time = 50.

    real aspect, L, dx, dy
    real, dimension(nx,ny) :: T

    aspect = real(nx)/real(ny)
    L = H*aspect
    dx = L/real(nx)
    dy = H/real(ny)

    T = initialize_T()
    call evolve_field()

    contains

        function initialize_T() result(T)
            implicit none
            real, parameter :: sigma = 0.2
            integer i, j
            real x, y
            real, dimension(nx,ny) :: T
            do i=1, nx
                do j=1, ny
                    x = real(i)/real(nx)*L
                    y = real(j)/real(ny)*H
                    T(i,j) = 10. + &
                             2.* ( &
                                   1./(2.*3.14*sigma**2.) * &
                                   exp(-1. * ( (x-L/2.)**2. + (y-H/2.)**2. ) / (2.*sigma**2.)) &
                                 )
                enddo
            enddo
        end function initialize_T

        subroutine heat_eqn()
            implicit none
            real, dimension(nx,ny) :: Tn
            real d2Tdx2, d2Tdy2
            integer i, j

            Tn(:,:) = T(:,:)

            !$omp parallel shared(T) private(i,j,d2Tdx2,d2Tdy2)
            !$omp do
            do i=2, nx-1
                do j=2, ny-1            
                    d2Tdx2 = ( Tn(i+1,j) - 2.*Tn(i,j) + Tn(i-1,j) ) / dx**2.
                    d2Tdy2 = ( Tn(i,j+1) - 2.*Tn(i,j) + Tn(i,j-1) ) / dy**2.
                    T(i,j) =  Tn(i,j) + dt*(alpha*(d2Tdx2 + d2Tdy2))
                end do
            end do
            !$omp end do                     
            !$omp end parallel

            T(:, 1)  = T(:, 2)
            T(:, ny) = T(:, ny-1)
            T(1, :)  = T(2, :)
            T(nx, :) = T(nx-1, :)

        end subroutine heat_eqn

        subroutine evolve_field()
            implicit none
            integer ts, total_ts, frames_total, output_period_ts, pic_counter
            real progress_pct
            character(len=16) :: pic_counter_str

            total_ts = ceiling(physical_time/dt)
            frames_total = 30*20
            output_period_ts = int(total_ts/frames_total)

            pic_counter = 0

            do ts = 0, total_ts

                if (mod(ts,output_period_ts) .eq. 0) then

                    progress_pct = 100.*real(ts)/real(total_ts)
                    print '(I8, F7.2, A)', ts, progress_pct, "%"

                    !! -- for plotting

                    !open(3, file="T.dat", access="stream")
                    !write(3) T(:,:)
                    !close(3)

                    !write (pic_counter_str,'(I5.5)') pic_counter
                    !call system('gnuplot -c plot3d.gnu '//trim(pic_counter_str))
                    !pic_counter = pic_counter + 1

                end if

                ! -----

                call heat_eqn()

                ! -----

            end do

        end subroutine evolve_field

end program heat

Answer 1

如注释中的[[evets所建议，]该程序由于堆栈内存不足而崩溃。也就是说，您创建了太多太大的自动（堆栈分配）阵列。作为指示，您可以在Valgrind中运行该程序

> gfortran -fopenmp heat.f90
> valgrind ./a.out

在打印很多“无效写入”错误之前，它将输出类似的内容

==31723== Warning: client switching stacks?  SP change: 0x1ffefffbc0 --> 0x1ffecf2740
==31723==          to suppress, use: --max-stackframe=3200128 or greater

当堆栈内存不足时，通常会出现有关“客户端切换堆栈”的警告。

在真实程序中，您将使用allocatable数组，这些数组不是人为限制的。但是，在这种简单情况下，我成功使用了标志-fno-automatic，该标志避免了使用自动（堆栈分配的）数组，而是将它们放置在另一种内存中：
> gfortran -fno-automatic -fopenmp heat.f90
> ./a.out
       0   0.00%
     500   0.17%
    1000   0.33%
    1500   0.50%
...

您可以在手册页中阅读有关该选项的详细信息。

注意：我在其他Linux发行版（不是Windows）中测试了该程序。