主题:谁能帮忙看下程序啊
我的程序不知道哪里出问题了,n=4,m=4时运行正常;但是如果n=24,m=16就提示说stack overflow。是关于求矩阵特征值的,我调试了下,到eigenvalue=eig(h,w=eigenvector)这条语句就出错啦。请大侠们帮忙看下啊,n=4,m=4时候矩阵较小,可能没超出范围,但是我定义的是动态数组啊,数组大小随着m,n变动,m、n变了矩阵也变了,怎么会提示超出范围呢?
program main
use IMSL
! use linear_operators
implicit none
real, parameter::pi=3.14159265
complex, parameter:: fi=(0.0, 1.0)
real t,f,w,temp
complex*8,allocatable::h(:,:)
complex*8,allocatable::h0(:,:)
complex*8,allocatable::r0(:,:)
complex*8,allocatable::l0(:,:)
complex*8,allocatable::eigenvector(:,:)
real*8,allocatable:: eigenvalue(:) !特征值的矩阵
integer ::i,n,m,j,p
n=4
m=4
allocate(h(4*m*n,4*m*n))
allocate(h0(4*n,4*n),r0(4*n,4*n),l0(4*n,4*n))
allocate(eigenvalue(4*n*m))
allocate(eigenvector(4*n*m,4*m*n))
w=0.0
t=1.0
f=0.25
call hmatr(h, n,m)
open(1,file='hmat.txt')
do i=1,4*n*m
do j=1,4*n*m
write(1,*) h(j,i)
end do
end do
close(1)
eigenvalue=eig(h,w=eigenvector)
do j=1,4*n*m-1
p=j
do i=j+1,4*m*n
if(eigenvalue(i)<eigenvalue(p)) p=i
end do
temp=eigenvalue(j)
eigenvalue(j)=eigenvalue(p)
eigenvalue(p)=temp
end do
open(10,file='eeig.txt')
do i=1,4*m*n
write(10,*) eigenvalue(i)
end do
close(10)
stop
contains
subroutine hmatr(h,n,m)
integer n,m
complex h(4*m*n,4*m*n)
complex h0(4*n,4*n)
complex r0(4*n,4*n)
complex l0(4*n,4*n)
integer:: j,n0,x
complex:: c(4,4),a(4,4),b(4,4),r(4,4),l(4,4)
h0=0.0
a=0.0
b=0.0
h=0.0
r0=0.0
l0=0.0
a(1, 4)=t
b=transpose(a)
c=0.0
r=0.0
l=0.0
n0=1
do while(n0<=2*n-1)
r(2,1)=t*exp(fi*(pi-2.0*pi*f/3.0*(3.0*n0/2.0+0.25)))
r(3,4)=t*exp(fi*(pi-2.0*pi*f/3.0*(3.0*(n0+1)/2.0+0.25)))
l(1,2)=t*exp(-fi*(pi-2.0*pi*f/3.0*(3.0*n0/2.0+0.25)))!conjg(r0(2,1))
l(4,3)=t*exp(-fi*(pi-2.0*pi*f/3.0*(3.0*(n0+1)/2.0+0.25)))!conjg(r0(3,4))
r0(2*n0,2*n0-1)=r(2,1)
r0(2*n0+1,2*n0+2)=r(3,4)
l0(2*n0-1,2*n0)=l(1,2)
l0(2*n0+2,2*n0+1)=l(4,3)
c(1,1)=w
c(2, 1) = t*exp(fi*(-pi+2.0*pi/3.0*f*(1.5*n0+0.25)))
c(1, 2) = t*exp(-fi*(-pi+2.0*pi/3.0*f*(1.5*n0+0.25)))
c(2, 2) = w
c(2, 3) = t
c(3, 2) = t
c(3, 3) = w
c(3, 4) = t*exp(fi*(-pi+2.0*pi/3.0*f*(1.5*(n0+1)+0.25)))
c(4, 3) = t*exp(-fi*(-pi+2.0*pi/3.0*f*(1.5*(n0+1)+0.25)))
c(4, 4) = w
h0(2*n0-1,2*n0-1)=c(1,1)
h0(2*n0,2*n0-1)=c(2,1)
h0(2*n0-1,2*n0)=c(1,2)
h0(2*n0,2*n0)=c(2,2)
h0(2*n0,2*n0+1)=c(2,3)
h0(2*n0+1,2*n0)=c(3,2)
h0(2*n0+1,2*n0+1)=c(3,3)
h0(2*n0+1,2*n0+2)=c(3,4)
h0(2*n0+2,2*n0+1)=c(4,3)
h0(2*n0+2,2*n0+2)=c(4,4)
if(2*n0-1==1) then
h0(2*n0-1:2*n0+2,2*n0+3:2*n0+6)=b
h0(2*n0-1:2*n0+2,4*n-3:4*n)=a
else if(2*n0-1<4*n-3) then
h0(2*n0-1:2*n0+2,2*n0+3:2*n0+6)=b
h0(2*n0-1:2*n0+2,2*n0-5:2*n0-2)=a
else
h0(4*n-3:4*n,1:4)=b
h0(2*n0-1:2*n0+2,2*n0-5:2*n0-2)=a
end if
n0=n0+2
end do
do x=1,4*(m-1)*n+1,4*n
h(x:x+4*n-1,x:x+4*n-1)=h0
if( x==1) then
h(x:x+4*n-1,x+4*n:8*n)=r0
h(x:x+4*n-1,(m-1)*4*n+1:m*4*n)=l0
elseif( x>=4*n+1 .and. x<4*(m-1)*n+1) then
h(x:x+4*n-1,x+4*n:x+8*n-1)=r0
h(x:x+4*n-1,x-4*n:x-1)=l0
else
h(x:x+4*n-1,1:4*n)=r0
h(x:x+4*n-1,x-4*n:x-1)=l0
end if
end do
end subroutine hmatr
end program main
program main
use IMSL
! use linear_operators
implicit none
real, parameter::pi=3.14159265
complex, parameter:: fi=(0.0, 1.0)
real t,f,w,temp
complex*8,allocatable::h(:,:)
complex*8,allocatable::h0(:,:)
complex*8,allocatable::r0(:,:)
complex*8,allocatable::l0(:,:)
complex*8,allocatable::eigenvector(:,:)
real*8,allocatable:: eigenvalue(:) !特征值的矩阵
integer ::i,n,m,j,p
n=4
m=4
allocate(h(4*m*n,4*m*n))
allocate(h0(4*n,4*n),r0(4*n,4*n),l0(4*n,4*n))
allocate(eigenvalue(4*n*m))
allocate(eigenvector(4*n*m,4*m*n))
w=0.0
t=1.0
f=0.25
call hmatr(h, n,m)
open(1,file='hmat.txt')
do i=1,4*n*m
do j=1,4*n*m
write(1,*) h(j,i)
end do
end do
close(1)
eigenvalue=eig(h,w=eigenvector)
do j=1,4*n*m-1
p=j
do i=j+1,4*m*n
if(eigenvalue(i)<eigenvalue(p)) p=i
end do
temp=eigenvalue(j)
eigenvalue(j)=eigenvalue(p)
eigenvalue(p)=temp
end do
open(10,file='eeig.txt')
do i=1,4*m*n
write(10,*) eigenvalue(i)
end do
close(10)
stop
contains
subroutine hmatr(h,n,m)
integer n,m
complex h(4*m*n,4*m*n)
complex h0(4*n,4*n)
complex r0(4*n,4*n)
complex l0(4*n,4*n)
integer:: j,n0,x
complex:: c(4,4),a(4,4),b(4,4),r(4,4),l(4,4)
h0=0.0
a=0.0
b=0.0
h=0.0
r0=0.0
l0=0.0
a(1, 4)=t
b=transpose(a)
c=0.0
r=0.0
l=0.0
n0=1
do while(n0<=2*n-1)
r(2,1)=t*exp(fi*(pi-2.0*pi*f/3.0*(3.0*n0/2.0+0.25)))
r(3,4)=t*exp(fi*(pi-2.0*pi*f/3.0*(3.0*(n0+1)/2.0+0.25)))
l(1,2)=t*exp(-fi*(pi-2.0*pi*f/3.0*(3.0*n0/2.0+0.25)))!conjg(r0(2,1))
l(4,3)=t*exp(-fi*(pi-2.0*pi*f/3.0*(3.0*(n0+1)/2.0+0.25)))!conjg(r0(3,4))
r0(2*n0,2*n0-1)=r(2,1)
r0(2*n0+1,2*n0+2)=r(3,4)
l0(2*n0-1,2*n0)=l(1,2)
l0(2*n0+2,2*n0+1)=l(4,3)
c(1,1)=w
c(2, 1) = t*exp(fi*(-pi+2.0*pi/3.0*f*(1.5*n0+0.25)))
c(1, 2) = t*exp(-fi*(-pi+2.0*pi/3.0*f*(1.5*n0+0.25)))
c(2, 2) = w
c(2, 3) = t
c(3, 2) = t
c(3, 3) = w
c(3, 4) = t*exp(fi*(-pi+2.0*pi/3.0*f*(1.5*(n0+1)+0.25)))
c(4, 3) = t*exp(-fi*(-pi+2.0*pi/3.0*f*(1.5*(n0+1)+0.25)))
c(4, 4) = w
h0(2*n0-1,2*n0-1)=c(1,1)
h0(2*n0,2*n0-1)=c(2,1)
h0(2*n0-1,2*n0)=c(1,2)
h0(2*n0,2*n0)=c(2,2)
h0(2*n0,2*n0+1)=c(2,3)
h0(2*n0+1,2*n0)=c(3,2)
h0(2*n0+1,2*n0+1)=c(3,3)
h0(2*n0+1,2*n0+2)=c(3,4)
h0(2*n0+2,2*n0+1)=c(4,3)
h0(2*n0+2,2*n0+2)=c(4,4)
if(2*n0-1==1) then
h0(2*n0-1:2*n0+2,2*n0+3:2*n0+6)=b
h0(2*n0-1:2*n0+2,4*n-3:4*n)=a
else if(2*n0-1<4*n-3) then
h0(2*n0-1:2*n0+2,2*n0+3:2*n0+6)=b
h0(2*n0-1:2*n0+2,2*n0-5:2*n0-2)=a
else
h0(4*n-3:4*n,1:4)=b
h0(2*n0-1:2*n0+2,2*n0-5:2*n0-2)=a
end if
n0=n0+2
end do
do x=1,4*(m-1)*n+1,4*n
h(x:x+4*n-1,x:x+4*n-1)=h0
if( x==1) then
h(x:x+4*n-1,x+4*n:8*n)=r0
h(x:x+4*n-1,(m-1)*4*n+1:m*4*n)=l0
elseif( x>=4*n+1 .and. x<4*(m-1)*n+1) then
h(x:x+4*n-1,x+4*n:x+8*n-1)=r0
h(x:x+4*n-1,x-4*n:x-1)=l0
else
h(x:x+4*n-1,1:4*n)=r0
h(x:x+4*n-1,x-4*n:x-1)=l0
end if
end do
end subroutine hmatr
end program main