PROGRAM MAIN 
USE IMSL
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      PARAMETER(PI=3.14159265358979D0)
      PARAMETER(THA=0.0D0*PI/180.0D0,THA1=-90.0D0*PI/180.0D0,
     &THA2=90.0D0*PI/180.0D0)
      PARAMETER(G=2.80D6,H=1.47370D4)
   
      DIMENSION A(8,8),EIGENVALUE(8),EIGENVECTOR(8,8)
      
      A(1,1)=THA
      A(1,2)=1.0D0
      A(1,3)=SIN(2.0D0*THA)
      A(1,4)=COS(2.0D0*THA)
      A(1,5)=-A(1,1)
      A(1,6)=-A(1,2)
      A(1,7)=-A(1,3)
      A(1,8)=-A(1,4)
      A(2,1)=1.0D0
      A(2,3)=2.0D0*COS(2.0D0*THA)
      A(2,4)=-2.0D0*SIN(2.0D0*THA)
A(3,5)=1.0D0
            


      A(3,7)=2.0D0*COS(2.0D0*THA)
      A(3,8)=-2.0D0*SIN(2.0D0*THA)
      A(4,1)=1.0D0/G
      A(4,5)=-1.0D0/H
      A(5,1)=THA1
      A(5,2)=1.0D0
      A(5,3)=SIN(2.0D0*THA1)
      A(5,4)=COS(2.0D0*THA1)
      A(6,1)=1.0D0
      A(6,3)=2.0D0*COS(2.0D0*THA1)
      A(6,4)=-2.0D0*SIN(2.0D0*THA1)
      A(7,5)=THA1
      A(7,6)=1.0D0
      A(7,7)=SIN(2.0D0*THA2)
      A(7,8)=COS(2.0D0*THA2)
      A(8,5)=1.0D0
      A(8,7)=2.0D0*COS(2.0D0*THA2)
      A(8,8)=-2.0D0*SIN(2.0D0*THA2)
      EIGENVALUE=EIG(A,V=EIGENVECTOR)
DO I=1,8
WRITE(*,*) 'EIGENVALUE',EIGENVALUE(I)
WRITE(*,'(8E20.10)') EIGENVALUE(I)*EIGENVECTOR(:,I)
WRITE(*,'(8E20.10)') MATMUL(A,EIGENVECTOR(:,I))
ENDDO

END

我是求一个矩阵的特征值特征向量，但是程序可以运行，编译后出现错误，请问是什么原因，谢谢！