fortran语言中seg错误的原因。

I want to find out all eigenvalues and relevant right eigenvectors using lapack. However, the code could just work in small matrix size. While the size of my matrix is over some specific number like 300, it collapsed and all I got is segmentation fault (core dumped). I checked my code several times and still could not find out where the problem is. Following is my code.

我想找出所有的特征值和相关的右特征向量。但是，代码可以在小的矩阵大小中工作。虽然我的矩阵的大小超过了一些特定的数字，比如300，但它崩溃了，我得到的只是分割错误(核心转储)。我检查了我的代码好几次，仍然无法找到问题所在。下面是我的代码。

program Eigenvalue
! finding the eigenvalues of a complex matrix using LAPACK
! http://www.netlib.org/lapack/complex16/zgeev.f
!find out the inverse complex matrix in fortran


Implicit none
!matrix dimension
integer, parameter::N=300
integer,parameter::M=3*(N-1)
! declarations, notice double precision
! A and B are matrices with given value while inver is inverse matrix of A
! S is to check whether inverse matrix is right or not        

real*8,allocatable,dimension(:,:)::A
real*8,allocatable,dimension(:,:)::Inver
real*8,allocatable,dimension(:,:)::B
real*8,allocatable,dimension(:,:)::S
!Mul is the matrix to be calculated in the subroutine and its values are
!product of inverse matrix and B
!eigval represents eigenvalues while Vec are eigenvectors
complex*16,allocatable,dimension(:,:):: Mul,DUMMY,Vec
complex*16,allocatable,dimension(:)::eigval,WORK1
double precision::WORK2(2*M)

integer i,j,ok,temp,error

!specific parameters in this case to define A and B 

real*8 ky,k,lb,ln,kz,ni,to
real*8 x,dx
real*8::x0=-30.d0,x1=30.d0
real*8 co1,co2,co3,co4

allocate(A(M,M),B(M,M),Inver(M,M),S(M,M),stat=error)
if (error.ne.0)then
  print *,"error:not enough memory"
  stop
end if

allocate(Mul(M,M),DUMMY(M,M),Vec(M,M),eigval(M),WORK1(2*M),stat=error)
if (error.ne.0)then
  print *,"error:not enough memory"
  stop
end if
!values of parameters 

ni=5
to=1
ln=10.d0
lb=50.d0
ky=0.4d0
kz=0.002d0
k=(1.d0+ni)/to
dx=(x1-x0)/N
x=x1-dx
co1=(1+ky*ky)*dx*dx-2
co2=ky*(-2*k+(1-k*ky*ky)*dx*dx)
co3=dx*dx*(kz-ky*(x/lb)**2)
co4=ky*k
call mkl_set_num_threads( 2  )

!definition of the initial matrix A

do i=1,M
  do j=1,N-1
        if(j.eq.i)then
        A(i,j)=-2+(1+ky*ky)*dx*dx
        else if (j.eq.(i+1)) then
        A(i,j)=-1
        else if (j.eq.(i-1).and.j.ne.(N-1)) then
        A(i,j)=-1
        else if(j.eq.(i-1).and.j.eq.(N-1))then
         A(i,j)=0
        else
        A(i,j)=(0,0)
        end if
   end do
   do j=N,M
       if(j.eq.i)then
        A(i,j)=1
       else
        A(i,j)=0
       end if
   end do
 end do
call inverse(A,Inver,M)
!definition of matrix B and the multiplication Mul

do i=1,M
      do j=1,N-1
        if(j.eq.i)then
        B(i,j)=ky*((1-k*ky*ky)*dx*dx-2*k)
        else if (j.eq.(i+1)) then
        B(i,j)=ky*k
        else if (j.eq.(i-1).and.j.ne.(N-1)) then
        B(i,j)=ky*k
        else if (j.eq.(i-1).and.j.eq.(N-1))then
        B(i,j)=0
        else if(j.eq.(i-N+1)) then
        B(i,j)=kz-((x0+j*dx)**2)*ky/(lb*lb)
        else if (j.eq.(i-2*N+2)) then
        B(i,j)=k*ky
        else
        B(i,j)=0
        end if
      end do
      do j=N,2*(N-1)
       if(j.eq.(i+N-1))then
        B(i,j)=(kz-((x0+i*dx)**2)*ky/(lb*lb))*dx*dx
       else
        B(i,j)=0
       end if
      end do
      do j=2*N-1,M
        if(j.eq.(i+N-1)) then
         B(i,j)=kz-((x0+(i-N+1)*dx)**2)*ky/(lb*lb)
        else
         B(i,j)=0
      end if
      end do
end do

Mul=matmul(Inver,B)
call ZGEEV('N', 'V', M,Mul, M, eigval, DUMMY, M, Vec, M, WORK1, 2*M, WORK2, ok)
deallocate(A,B,Inver,S,stat=error)
if (error.ne.0)then
  print *,"error:fail to release"
  stop
end if
deallocate(Mul,DUMMY,Vec,eigval,WORK1,stat=error)
if (error.ne.0)then
  print *,"error:failed to release"
  stop
end if
end
subroutine inverse(a,c,n)

!============================================================ 
! Inverse matrix
! Method: Based on Doolittle LU factorization for Ax=b
! Alex G. December 2009
!this is for real numbers
!-----------------------------------------------------------
! input ...
! a(n,n) - array of coefficients for matrix A
! n      - dimension
! output ...
! c(n,n) - inverse matrix of A
! comments ...
! the original matrix a(n,n) will be destroyed 
! during the calculation
!===========================================================

implicit none
integer n
double precision a(n,n), c(n,n)
double precision,allocatable,dimension(:,:):: L, U
double precision,allocatable,dimension(:):: b, d, x
double precision coeff
integer i, j, k,error

allocate(L(n,n),U(n,n),b(n),d(n),x(n),stat=error)
 if (error.ne.0)then
  print *,"error:not enough memory"
  stop
 end if

! step 0: initialization for matrices L and U and b
! Fortran 90/95 allows such operations on matrices

L=0.0
U=0.0
b=0.0

! step 1: forward elimination

do k=1, n-1
 do i=k+1,n
  coeff=a(i,k)/a(k,k)
  L(i,k) = coeff
  do j=k+1,n
     a(i,j) = a(i,j)-coeff*a(k,j)
  end do
 end do
end do

! Step 2: prepare L and U matrices 
! L matrix is a matrix of the elimination coefficient
! + the diagonal elements are 1.0

do i=1,n
  L(i,i) = 1.0
end do
! U matrix is the upper triangular part of A

do j=1,n
  do i=1,j
   U(i,j) = a(i,j)
  end do
end do

! Step 3: compute columns of the inverse matrix C

do k=1,n
  b(k)=1.0
  d(1) = b(1)
! Step 3a: Solve Ld=b using the forward substitution

  do i=2,n
    d(i)=b(i)
    do j=1,i-1
      d(i) = d(i) - L(i,j)*d(j)
    end do
  end do
! Step 3b: Solve Ux=d using the back substitution

  x(n)=d(n)/U(n,n)
  do i = n-1,1,-1
    x(i) = d(i)
    do j=n,i+1,-1
     x(i)=x(i)-U(i,j)*x(j)
    end do
  x(i) = x(i)/u(i,i)
  end do
! Step 3c: fill the solutions x(n) into column k of C

 do i=1,n
  c(i,k) = x(i)
 end do
  b(k)=0.0
end do
  deallocate(L,U,b,d,x,stat=error)
if (error.ne.0)then
  print *,"error:failed to deallocate"
  stop
end if
end subroutine inverse

Hope someone can give me a tip or advice about why the segmentation fault occurs when N is over 300.

希望有人能给我一个建议，告诉我为什么当N超过300的时候会发生分割错误。

2 个解决方案

#1

In your call to ZGEEV, could it be that the dimension of WORK1 (argument after WORK1) is incorrect, shouldn't it be 2*M instead of 3*M? Also, WORK2 should be a double precision array.

在您对ZGEEV的调用中，可能是WORK1的维度(工作1之后的参数)是不正确的，应该是2*M而不是3*M吗?另外，WORK2应该是一个双精度数组。

With the updated information, the problem seems to be in the inverse subroutine. It's a combination of using automatic arrays, and the fact that the Intel compiler puts all automatic arrays on the stack. With large numbers for n, L and U together take up almost 16 MB on the stack. You can either not use automatic arrays, or compile with the option -heap-arrays. ALternatively, on linux, you can use ulimit -s unlimited to allow for unlimited stack size.

有了更新的信息，问题似乎在逆子例程中。这是一个使用自动数组的组合，以及Intel编译器将所有自动数组放在堆栈上的事实。大量的n, L和U一起占用了将近16mb的堆栈。您可以不使用自动数组，也可以使用选项-heap-array编译。或者，在linux上，您可以使用ulimit -s unlimited来允许无限的堆栈大小。

#2

Typical causes of segmentation fault in Fortran are accessing an array beyond its declared size or inconsistencies between actual arguments in the caller and the dummy arguments of the called procedure. You might find a case of the first problem by turning on runtime error checking, especially subscript checking. The second can be found by the compiler when all of your routines have explicit interfaces. The easiest way is to have your procedures in modules and use the modules. This enables the compiler to check consistency of arguments.

Fortran中分割错误的典型原因是访问一个超出其声明大小的数组或调用程序中的实际参数与被调用过程的虚拟参数之间的不一致。您可以通过打开运行时错误检查，特别是下标检查来发现第一个问题。当所有例程都有显式接口时，编译器可以找到第二个。最简单的方法是在模块中设置程序并使用模块。这使编译器能够检查参数的一致性。