1

我想使用 lapack 找出所有特征值和相关的右特征向量。但是，代码只能在小矩阵大小下工作。虽然我的矩阵的大小超过了 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

希望有人能给我一个提示或建议，说明为什么当 N 超过 300 时会发生分段错误。

4

2 回答 2

1

在您对 ZGEEV 的调用中，WORK1(argument after WORK1) 的维度是否不正确，不应该是2*M吗3*M？此外，WORK2应该是一个双精度数组。

有了更新的信息，问题似乎出在inverse子程序上。它结合了使用自动数组和英特尔编译器将所有自动数组放入堆栈这一事实。大量的n,L和U一起在堆栈上占用了将近 16 MB。您可以不使用自动数组，也可以使用选项编译-heap-arrays。或者，在 linux 上，您可以使用ulimit -s unlimited允许无限的堆栈大小。

于 2013-09-11T08:07:30.800 回答
0

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

于 2013-09-11T08:29:03.100 回答