causes of seg fault in fortran

Question

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.

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.

Answer 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.

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.

Answer 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.