程序收到信号 SIGSEGV：分段错误 - 大数组中的无效内存引用
Program received signal SIGSEGV: Segmentation fault - invalid memory reference in arrays with big sizes
我在 运行 此代码时遇到错误。当我 运行 带有小 L 的代码（如 L=16 或 L=32）时，我没有收到任何错误，但在 L = 128 或 L=96 中，经过 7000-8000 步后，我收到以下错误：
Program received signal SIGSEGV: Segmentation fault - invalid memory reference.
Backtrace for this error:
#0 0x7FBA5CAC3E08
#1 0x7FBA5CAC2F90
#2 0x7FBA5C1E84AF
#3 0x402769 in MAIN__ at newhys.f90:?
Segmentation fault (core dumped)
这是完整代码：
SUBROUTINE init_random_seed()
implicit none
INTEGER :: i, n, clock
INTEGER, DIMENSION(:), ALLOCATABLE :: seed
CALL RANDOM_SEED(size = n)
ALLOCATE(seed(n))
CALL SYSTEM_CLOCK(COUNT=clock)
seed = clock + 37 * (/ (i - 1, i = 1, n) /)
CALL RANDOM_SEED(PUT = seed)
DEALLOCATE(seed)
END SUBROUTINE
!end module
Program Activ_mater
USE OMP_LIB
Implicit none
Integer,parameter :: time=1000000000, L=128, N = L**2*2
Integer,parameter:: n_thread = 8
Real(8),parameter :: pi = 3.14159265359
Real(8),parameter :: v0 = 0.50, alpha = 1.0/36.0
real(8)START,END_1 ,eta
type block_p
Integer :: partical_N
Integer :: particle_ad(10*L)
end type
Type(block_p) ,pointer,dimension(:,:) :: C
Real(8),allocatable :: x(:), y(:) , phi,angle_new(:),angle(:)
Real(8) :: sum_a, sum_b,x_in, y_in, x_out, y_out, avrage_t, r,ra,eta1(5)
Integer :: i,j,t,n_p,I_b,J_b,b_l,neighbor_i(9),neighbor_j(9),A,n_p_b,ne = 1,stateta=0,ot=0,op=175
character(len=10)::name1
call omp_set_num_threads(n_thread)
call init_random_seed()
eta1=(/2.100,2.116,2.120,2.126,2.128/) ! The value of ETA
allocate(x(2*n), y(2*n) , phi,angle_new(2*N),angle(2*N))
allocate(C(2*L,2*L))
C(:,:)%partical_N=0
do i =1,N
call random_number(ra)
x(i)=ra*L
I_b = int(x(i))+1
call random_number(ra)
y(i)=ra*L
J_b = int(y(i))+1
call random_number(ra)
angle(i)=ra*2.0*pi
C(I_b,J_b)%partical_N = C(I_b,J_b)%partical_N + 1 !Number of particle in block C(I_b,J_b)
C(I_b,J_b)%particle_ad( C(I_b,J_b)%partical_N ) = i ! The particle number in block C(I_b,J_b)
end do
! loop for eta
eta= 0.0
write(name1,'(f5.3)')eta
open(unit=10, file='Hysteresis,'//trim(name1)//'.dat')
!=====================explanation of system====================================
print*,'==========================================================================='
print*, 'eta = ', etA ,' ',' alpha = ',alpha
print*,'L=',L ,' ', 'Particle Number=', N,' ','Density=', N/L**2
print*,'==========================================================================='
!==============================================================================
START = omp_get_wtime()
do t =1,time
if (ot == 300000 )then
stateta = 0
ot = 0
op = op + 1
end if
if (stateta == 0 )then
eta = eta + ((1.0/3.0) * (10E-6))
endif
if (int(eta * 100) == op) then
stateta = 1
end if
angle_new(:)=0
!$OMP PARALLEL DEFAULT(PRIVATE) SHARED(x,y,angle,angle_new,c) firstprivate(eta)
!$OMP DO schedule(runtime)
do i =1, N
sum_a=0; sum_b=0;n_p=0
I_b = int(x(i))+1; J_b = int(y(i))+1 ! The block of particle i
! Now I should find nine neighbor of particle i-----------------------------------------------
neighbor_i=(/I_b+1, I_b, I_b-1, I_b, I_b+1, I_b-1, I_b-1, I_b+1, I_b/)
neighbor_j=(/J_b, J_b+1, J_b, J_b-1, J_b+1, J_b+1, J_b-1, J_b-1, J_b/)
do b_l = 1, 9
I_b = neighbor_i(b_l) ; J_b=neighbor_j(b_l)
if (I_b >L )I_b=1
if (I_b <1 )I_b=L
if (J_b >L )J_b=1
if (J_b <1 )J_b=L
!neighbor_i(b_l)=I_b; neighbor_j(b_l)=J_b
A = C( I_b, J_b )%partical_N ! number of particle in block C( neighbor_i(b_l), neighbor_j(b_l) )
!=============================================================================================
do n_p_b =1, A
j = C( I_b, J_b)%particle_ad(n_p_b) !particle j in the block C
if (i /= j )then
X_in = abs(max(x(i),x(j)) - min(x(i),x(j)));
Y_in = abs(max(y(i),y(j)) - min(y(i),y(j)));
X_out =L-X_in
Y_out =L-Y_in
r = sqrt(min(X_in,X_out)**2 + min(Y_in,Y_out)**2)
else
r=0.0
end if
if ( r < 1 )then
if ( j <= i )then
sum_A = sum_A + sin(angle(j));
sum_B = sum_B + cos(angle(j));
else
sum_A = sum_A + alpha*sin(angle(j));
sum_B = sum_B + alpha*cos(angle(j));
endif
n_p = n_p + 1;
endif
enddo
enddo
sum_A = sum_A/n_p; sum_B = sum_B/n_p
!if (int(sum_A*1e10) ==0 .and. int(sum_B*1e10) ==0 )print*,'zerrooo'
avrage_t=atan2(sum_A,sum_B);
if (avrage_t<0.0) then
avrage_t=avrage_t+2.0*pi;
endif
call random_number(ra)
angle_new(i)=avrage_t+eta*(ra-0.50)
if( angle_new(i)>=2*pi) angle_new(i)= angle_new(i)-2*pi
if( angle_new(i)<0) angle_new(i)= 2*pi+angle_new(i)
end do
!$OMP END DO
!$OMP END PARALLEL
angle = angle_new
C(:,:)%partical_N=0
! phi=0.0
do i=1, N
x(i) = x(i) + v0*sin(angle(i));
if (x(i)<1) x(i)=L+x(i)
if (x(i)>L) x(i)=x(i)-L
I_b = int(x(i))+1
y(i) = y(i)+ v0*cos(angle(i));
if (y(i)<1) y(i)=L+y(i)
if (y(i)>L) y(i)=y(i)-L
J_b = int(y(i))+1
C(I_b,J_b)%partical_N = C(I_b,J_b)%partical_N + 1
C(I_b,J_b)%particle_ad( C(I_b,J_b)%partical_N ) = i
end do
if (stateta == 1 )then
phi= sqrt((sum(sin(angle))**2+sum(cos(angle))**2))/N;
ot = ot + 1
end if
write(10,*)phi
if (mod(t,10)==0)then
! ave4_phi=sum(phi**4)/t;
! ave2_phi =sum(phi**2)/t;
!print* ,ave4_phi,ave2_phi
print*,'Time=',t,' ==== Eta : ',eta,"Ot : " , ot
end if
end do
END_1 = omp_get_wtime()
print*,'Run Time = ',end_1-start
End Program
P.S.(1) ：我使用 omp lib 运行 我的程序并行
P.S.(2) : 我使用 gfortran 编译代码
P.S.(3) ：使用 -g -fcheck=all 编译的代码并给我这个错误：
At line 155 of file z.f90 Fortran runtime error: Index '1281' of dimension 1 of array 'c%particle_ad' above upper bound of 1280
谢谢大家
您的 particle_ad 阵列只有 space 个用于 10*L 个粒子，但您似乎试图在其中存储最多 N=2*L**2 个粒子（取决于随机数下降）。平均而言，每个人都会有足够的 space，但是如果太多粒子落入一个块中，您的代码将失败，并且（如果我没记错的话）发生这种情况的机会会随着 L增加。
您可以通过将 Integer :: particle_ad(10*L) 替换为 Integer :: particle_ad(N) 来解决此问题，但这会浪费大量的 space。
一个更好的解决方案是动态 re-size particle_ad 数组，使它们在每次变满时变大。为了方便起见，您可以将此行为包装在 block_p class.
的方法中
例如，
type block_p
Integer :: partical_N
Integer, allocatable :: particle_ad(:)
contains
subroutine add_particle(this, index)
class(block_p), intent(inout) :: this
integer, intent(in) :: index
integer, allocatable :: temp
! Update partical_N.
this%partical_N = this%partical_N + 1
! Resize particle_ad if it is full.
if (size(this%particle_ad)<this%partical_N) then
temp = this%particle_ad
deallocate(this%particle_ad)
allocate(this%particle_ad(2*this%partical_N))
this%particle_ad(:size(temp)) = temp
endif
! Add the new index to particle_ad.
this%particle_ad(this%partical_N) = index
end subroutine
end type
然后您可以替换这些行
C(I_b,J_b)%partical_N = C(I_b,J_b)%partical_N + 1
C(I_b,J_b)%particle_ad( C(I_b,J_b)%partical_N ) = i
和
call C(I_b,J_b)%add_particle(i)
请注意，由于每个 particle_ad 数组现在是 allocatable，您需要先初始化每个数组，然后才能调用 add_particle.
我在 运行 此代码时遇到错误。当我 运行 带有小 L 的代码（如 L=16 或 L=32）时，我没有收到任何错误，但在 L = 128 或 L=96 中，经过 7000-8000 步后，我收到以下错误：
Program received signal SIGSEGV: Segmentation fault - invalid memory reference.
Backtrace for this error:
#0 0x7FBA5CAC3E08
#1 0x7FBA5CAC2F90
#2 0x7FBA5C1E84AF
#3 0x402769 in MAIN__ at newhys.f90:?
Segmentation fault (core dumped)
这是完整代码：
SUBROUTINE init_random_seed()
implicit none
INTEGER :: i, n, clock
INTEGER, DIMENSION(:), ALLOCATABLE :: seed
CALL RANDOM_SEED(size = n)
ALLOCATE(seed(n))
CALL SYSTEM_CLOCK(COUNT=clock)
seed = clock + 37 * (/ (i - 1, i = 1, n) /)
CALL RANDOM_SEED(PUT = seed)
DEALLOCATE(seed)
END SUBROUTINE
!end module
Program Activ_mater
USE OMP_LIB
Implicit none
Integer,parameter :: time=1000000000, L=128, N = L**2*2
Integer,parameter:: n_thread = 8
Real(8),parameter :: pi = 3.14159265359
Real(8),parameter :: v0 = 0.50, alpha = 1.0/36.0
real(8)START,END_1 ,eta
type block_p
Integer :: partical_N
Integer :: particle_ad(10*L)
end type
Type(block_p) ,pointer,dimension(:,:) :: C
Real(8),allocatable :: x(:), y(:) , phi,angle_new(:),angle(:)
Real(8) :: sum_a, sum_b,x_in, y_in, x_out, y_out, avrage_t, r,ra,eta1(5)
Integer :: i,j,t,n_p,I_b,J_b,b_l,neighbor_i(9),neighbor_j(9),A,n_p_b,ne = 1,stateta=0,ot=0,op=175
character(len=10)::name1
call omp_set_num_threads(n_thread)
call init_random_seed()
eta1=(/2.100,2.116,2.120,2.126,2.128/) ! The value of ETA
allocate(x(2*n), y(2*n) , phi,angle_new(2*N),angle(2*N))
allocate(C(2*L,2*L))
C(:,:)%partical_N=0
do i =1,N
call random_number(ra)
x(i)=ra*L
I_b = int(x(i))+1
call random_number(ra)
y(i)=ra*L
J_b = int(y(i))+1
call random_number(ra)
angle(i)=ra*2.0*pi
C(I_b,J_b)%partical_N = C(I_b,J_b)%partical_N + 1 !Number of particle in block C(I_b,J_b)
C(I_b,J_b)%particle_ad( C(I_b,J_b)%partical_N ) = i ! The particle number in block C(I_b,J_b)
end do
! loop for eta
eta= 0.0
write(name1,'(f5.3)')eta
open(unit=10, file='Hysteresis,'//trim(name1)//'.dat')
!=====================explanation of system====================================
print*,'==========================================================================='
print*, 'eta = ', etA ,' ',' alpha = ',alpha
print*,'L=',L ,' ', 'Particle Number=', N,' ','Density=', N/L**2
print*,'==========================================================================='
!==============================================================================
START = omp_get_wtime()
do t =1,time
if (ot == 300000 )then
stateta = 0
ot = 0
op = op + 1
end if
if (stateta == 0 )then
eta = eta + ((1.0/3.0) * (10E-6))
endif
if (int(eta * 100) == op) then
stateta = 1
end if
angle_new(:)=0
!$OMP PARALLEL DEFAULT(PRIVATE) SHARED(x,y,angle,angle_new,c) firstprivate(eta)
!$OMP DO schedule(runtime)
do i =1, N
sum_a=0; sum_b=0;n_p=0
I_b = int(x(i))+1; J_b = int(y(i))+1 ! The block of particle i
! Now I should find nine neighbor of particle i-----------------------------------------------
neighbor_i=(/I_b+1, I_b, I_b-1, I_b, I_b+1, I_b-1, I_b-1, I_b+1, I_b/)
neighbor_j=(/J_b, J_b+1, J_b, J_b-1, J_b+1, J_b+1, J_b-1, J_b-1, J_b/)
do b_l = 1, 9
I_b = neighbor_i(b_l) ; J_b=neighbor_j(b_l)
if (I_b >L )I_b=1
if (I_b <1 )I_b=L
if (J_b >L )J_b=1
if (J_b <1 )J_b=L
!neighbor_i(b_l)=I_b; neighbor_j(b_l)=J_b
A = C( I_b, J_b )%partical_N ! number of particle in block C( neighbor_i(b_l), neighbor_j(b_l) )
!=============================================================================================
do n_p_b =1, A
j = C( I_b, J_b)%particle_ad(n_p_b) !particle j in the block C
if (i /= j )then
X_in = abs(max(x(i),x(j)) - min(x(i),x(j)));
Y_in = abs(max(y(i),y(j)) - min(y(i),y(j)));
X_out =L-X_in
Y_out =L-Y_in
r = sqrt(min(X_in,X_out)**2 + min(Y_in,Y_out)**2)
else
r=0.0
end if
if ( r < 1 )then
if ( j <= i )then
sum_A = sum_A + sin(angle(j));
sum_B = sum_B + cos(angle(j));
else
sum_A = sum_A + alpha*sin(angle(j));
sum_B = sum_B + alpha*cos(angle(j));
endif
n_p = n_p + 1;
endif
enddo
enddo
sum_A = sum_A/n_p; sum_B = sum_B/n_p
!if (int(sum_A*1e10) ==0 .and. int(sum_B*1e10) ==0 )print*,'zerrooo'
avrage_t=atan2(sum_A,sum_B);
if (avrage_t<0.0) then
avrage_t=avrage_t+2.0*pi;
endif
call random_number(ra)
angle_new(i)=avrage_t+eta*(ra-0.50)
if( angle_new(i)>=2*pi) angle_new(i)= angle_new(i)-2*pi
if( angle_new(i)<0) angle_new(i)= 2*pi+angle_new(i)
end do
!$OMP END DO
!$OMP END PARALLEL
angle = angle_new
C(:,:)%partical_N=0
! phi=0.0
do i=1, N
x(i) = x(i) + v0*sin(angle(i));
if (x(i)<1) x(i)=L+x(i)
if (x(i)>L) x(i)=x(i)-L
I_b = int(x(i))+1
y(i) = y(i)+ v0*cos(angle(i));
if (y(i)<1) y(i)=L+y(i)
if (y(i)>L) y(i)=y(i)-L
J_b = int(y(i))+1
C(I_b,J_b)%partical_N = C(I_b,J_b)%partical_N + 1
C(I_b,J_b)%particle_ad( C(I_b,J_b)%partical_N ) = i
end do
if (stateta == 1 )then
phi= sqrt((sum(sin(angle))**2+sum(cos(angle))**2))/N;
ot = ot + 1
end if
write(10,*)phi
if (mod(t,10)==0)then
! ave4_phi=sum(phi**4)/t;
! ave2_phi =sum(phi**2)/t;
!print* ,ave4_phi,ave2_phi
print*,'Time=',t,' ==== Eta : ',eta,"Ot : " , ot
end if
end do
END_1 = omp_get_wtime()
print*,'Run Time = ',end_1-start
End Program
P.S.(1) ：我使用 omp lib 运行 我的程序并行
P.S.(2) : 我使用 gfortran 编译代码
P.S.(3) ：使用 -g -fcheck=all 编译的代码并给我这个错误：
At line 155 of file z.f90 Fortran runtime error: Index '1281' of dimension 1 of array 'c%particle_ad' above upper bound of 1280
谢谢大家
您的 particle_ad 阵列只有 space 个用于 10*L 个粒子，但您似乎试图在其中存储最多 N=2*L**2 个粒子（取决于随机数下降）。平均而言，每个人都会有足够的 space，但是如果太多粒子落入一个块中，您的代码将失败，并且（如果我没记错的话）发生这种情况的机会会随着 L增加。
您可以通过将 Integer :: particle_ad(10*L) 替换为 Integer :: particle_ad(N) 来解决此问题，但这会浪费大量的 space。
一个更好的解决方案是动态 re-size particle_ad 数组，使它们在每次变满时变大。为了方便起见，您可以将此行为包装在 block_p class.
例如，
type block_p
Integer :: partical_N
Integer, allocatable :: particle_ad(:)
contains
subroutine add_particle(this, index)
class(block_p), intent(inout) :: this
integer, intent(in) :: index
integer, allocatable :: temp
! Update partical_N.
this%partical_N = this%partical_N + 1
! Resize particle_ad if it is full.
if (size(this%particle_ad)<this%partical_N) then
temp = this%particle_ad
deallocate(this%particle_ad)
allocate(this%particle_ad(2*this%partical_N))
this%particle_ad(:size(temp)) = temp
endif
! Add the new index to particle_ad.
this%particle_ad(this%partical_N) = index
end subroutine
end type
然后您可以替换这些行
C(I_b,J_b)%partical_N = C(I_b,J_b)%partical_N + 1
C(I_b,J_b)%particle_ad( C(I_b,J_b)%partical_N ) = i
和
call C(I_b,J_b)%add_particle(i)
请注意，由于每个 particle_ad 数组现在是 allocatable，您需要先初始化每个数组，然后才能调用 add_particle.