回 帖 发 新 帖 刷新版面

主题:谁能帮忙看下程序啊

我的程序不知道哪里出问题了,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

回复列表 (共1个回复)

沙发

找到解决方案了

我来回复

您尚未登录,请登录后再回复。点此登录或注册