PROGRAM MAIN
    USE MSIMSL
    IMPLICIT REAL*8(A-H,O-Z)
    PARAMETER(MNEM=8,NA=10,NB=NA,MNPNT=MNEM*(NA+2)-(MNEM-1),
     &MNODS=2*MNPNT)
    PARAMETER(NODV=2,NOD=(NA+2),NODE=NODV*NOD,NOD2=2*NOD)
    PARAMETER(PI=3.141592654)
    DIMENSION NDEM(MNEM,NOD),NEVALC(50),NEVALR(50)
    DIMENSION A(MNODS,MNODS),B(MNODS,MNODS),C(MNODS,MNODS)
    DIMENSION SE1(NOD2,NOD2),SE2(NOD2,NOD2),SE3(NOD2,NOD2)
    DIMENSION AINV(MNODS,MNODS),SM(2*MNODS,2*MNODS)
    DIMENSION AINVB(MNODS,MNODS),AINVC(MNODS,MNODS)
    DIMENSION YONG(MNEM),RNU(MNEM),THTA(MNPNT),ETHTA(NOD)
    COMPLEX(8) EVAL(2*MNODS),EVEC(2*MNODS,2*MNODS)
    COMMON/GAUSSP/GAUSP(17),GAUSW(17)
    COMMON/NG/NGAUSP
    COMMON/VAR/IEM
    COMMON/PROBM/NTYPE
    OPEN(9,FILE='C:\mshing\eignv\r1180IM11.dat')
    OPEN(90,FILE='C:\mshing\eignv\r1180Re24.dat')
    NGAUSP=17
    NTYPE=1
      TH1=180.0D0*PI/180.0D0
    TH2=180.0D0*PI/180.0D0
    NEM=MNEM
    NPOINT=MNPNT
    NODS=MNODS
    NEMH=NEM/2
    DO 5 IE=1,NEMH
    YONG(IE)=100.0D0
    RNU(IE)=0.1D0
    5 CONTINUE
    DO 6 IE=NEMH+1,NEM
    YONG(NEMH+IE)=100D0
    RNU(NEMH+IE)=0.10D0
    6 CONTINUE
      DTH=(TH1+TH2)/(MNPNT-1)
    DO 8 I=1,NPOINT
    THTA(I)=(-TH1)+(I-1)*DTH
    8 CONTINUE
      DO 15 INOD=1,NOD
    NDEM(1,INOD)=INOD
   15 CONTINUE
      DO 18 NE=2,NEM
    DO 19 INOD=1,NOD
    NDEM(NE,INOD)=NDEM(NE-1,INOD)+NOD-1
   19 CONTINUE
   18 CONTINUE
      DO 900 IE=1,NEM
    DO 900 INOD=1,NOD
    WRITE(*,*) IE,NDEM(IE,INOD)
  900 CONTINUE
    DO 20 IRS=1,NODS
    DO 20 ICS=1,NODS
    A(IRS,ICS)=0.0D0
    B(IRS,ICS)=0.0D0
    C(IRS,ICS)=0.0D0
   20 CONTINUE
C
      CALL GAUSS(-1.0D0,1.0D0,GAUSP,GAUSW,NGAUSP)
C     
C     SUM OVER ALL ELEMENTS IN SYSTEM
C
      DO 30 IEM=1,NEM
    DO 32 ID=1,NOD
    NNDEM=NDEM(IEM,ID)
    ETHTA(ID)=THTA(NNDEM)
   32 CONTINUE
    CALL ELESTF(SE1,SE2,SE3,YONG,RNU,ETHTA,NEM,NOD)                 
    DO 35 INOD=1,NOD
    IRS=(NDEM(IEM,INOD)-1)*NODV
    DO 40 INDV=1,NODV
    IRS=IRS+1
    IRE=(INOD-1)*NODV+INDV
    DO 45 JNOD=1,NOD
    ICS=(NDEM(IEM,JNOD)-1)*NODV
    DO 50 JNDV=1,NODV
    ICS=ICS+1
    ICE=(JNOD-1)*NODV+JNDV 
    A(IRS,ICS)=A(IRS,ICS)+SE1(IRE,ICE)
    B(IRS,ICS)=B(IRS,ICS)+SE2(IRE,ICE)
    C(IRS,ICS)=C(IRS,ICS)+SE3(IRE,ICE)
   50 CONTINUE
   45    CONTINUE
   40 CONTINUE
   35 CONTINUE
   30 CONTINUE
C    
C     CALCULATE THE INVERSION OF MATRIX A
C     
    !CALL DLSGRR(NODS,NODS,A,NODS,1.0D-14,IRANK,AINV,NODS)
    CALL DLINRG(NODS,A,NODS,AINV,NODS)
C     
C     CALL THE SUBROUTINE MULTIPLIES THE A-MATRIX BY THE B-MATRIX
C  
      CALL ABMUT(AINV,NODS,NODS,B,NODS,AINVB)
    CALL ABMUT(AINV,NODS,NODS,C,NODS,AINVC)
      DO 60 I=1,2*NODS
    DO 65 J=1,2*NODS
    SM(I,J)=0.0D0
   65 CONTINUE
   60 CONTINUE
    DO 70 I=1,NODS
    DO 70 J=1,NODS
    SM(I,J)=-AINVB(I,J)
   70 CONTINUE
    DO 80 I=1,NODS
    DO 80 J=1,NODS
    SM(I,NODS+J)=-AINVC(I,J)
   80 CONTINUE
      DO 90 I=1,NODS
    SM(NODS+I,I)=1.0D0
   90 CONTINUE
      CALL DEVCRG(2*NODS,SM,2*NODS,EVAL,EVEC,2*NODS)
c    WRITE(1,300)
c 300 FORMAT(28X,'EIGENVALUES:'/)
c     WRITE(1,400)
c 400 FORMAT(15X,'REAL',37X,'IMAGE')
C
C     WRITE ONLY THE EIGENVALUES BETWWEEN -1<Re(Ramad)<0
C
      JCONTC=0
    JCONTR=0
    DO 100 IEVAL=1,2*NODS
    IF(REAL(EVAL(IEVAL)).GT.-1.0D0.AND.
    &REAL(EVAL(IEVAL)).LT.5.1D0.AND.AIMAG(EVAL(IEVAL)).GT.0.0D0) THEN
      JCONTC=JCONTC+1
    NEVALC(JCONTC)=IEVAL
    !WRITE(9,*) 'JCONTC=',JCONTC
c    WRITE(9,200) IEVAL,REAL(EVAL(IEVAL)),IEVAL,AIMAG(EVAL(IEVAL))
      WRITE(9,200) REAL(EVAL(IEVAL)),AIMAG(EVAL(IEVAL))
    END IF
  100 CONTINUE
      DO 105 IEVAL=1,2*NODS
     IF(REAL(EVAL(IEVAL)).GT.-1.0D0.AND.
    &REAL(EVAL(IEVAL)).LT.5.1D0.AND.DABS(AIMAG(EVAL(IEVAL)))
     &.LT.1.0D-14) THEN
      JCONTR=JCONTR+1
    NEVALR(JCONTR)=IEVAL
c    !WRITE(90,*) 'JCONTR=',JCONTR
c    WRITE(90,200) IEVAL,REAL(EVAL(IEVAL)),IEVAL,AIMAG(EVAL(IEVAL))
      WRITE(90,200) REAL(EVAL(IEVAL)),AIMAG(EVAL(IEVAL))
    END IF
  105 CONTINUE
c  200 FORMAT(2X,4HXin(,I3,2H)=,E15.6,2X,5HEtan(,I3,2H)=,E15.6)
  200 FORMAT(2X,E15.6,5X,5E15.6)
      !WRITE(1,500)
C  500 FORMAT(//,28X,'EIGENVECTORS:'/)
C
C     WRITE ONLY EIGENVECTORS CORRESPONDING THE  EIGENVALUES
C
      IF(JCONTC.NE.0) THEN
    DO 600 J=1,JCONTC 
    JEVEC=NEVALC(J)
    DO 600 IEVAL=1,NODS,2
c      WRITE(9,700) IEVAL,JEVEC,REAL(EVEC(NODS+IEVAL,JEVEC)),
c     &IEVAL,JEVEC,AIMAG(EVEC(NODS+IEVAL,JEVEC))
c      WRITE(9,800) IEVAL+1,JEVEC,REAL(EVEC(NODS+IEVAL+1,JEVEC)),
c     &IEVAL+1,JEVEC,AIMAG(EVEC(NODS+IEVAL+1,JEVEC))
C      WRITE(9,700) REAL(EVEC(NODS+IEVAL,JEVEC)),
C     &AIMAG(EVEC(NODS+IEVAL,JEVEC))
C    WRITE(9,800) REAL(EVEC(NODS+IEVAL+1,JEVEC)),
C     &AIMAG(EVEC(NODS+IEVAL+1,JEVEC))
  600 CONTINUE
     END IF
    IF(JCONTR.NE.0) THEN
      DO 605 J=1,JCONTR 
     JEVEC=NEVALR(J)
    DO 605 IEVAL=1,NODS,2
c      WRITE(90,700) IEVAL,JEVEC,REAL(EVEC(NODS+IEVAL,JEVEC)),
c     &IEVAL,JEVEC,AIMAG(EVEC(NODS+IEVAL,JEVEC))
c      WRITE(90,800) IEVAL+1,JEVEC,REAL(EVEC(NODS+IEVAL+1,JEVEC)),
c     &IEVAL+1,JEVEC,AIMAG(EVEC(NODS+IEVAL+1,JEVEC))
C    WRITE(90,700) REAL(EVEC(NODS+IEVAL,JEVEC)),
C     &AIMAG(EVEC(NODS+IEVAL,JEVEC))
C      WRITE(90,800) REAL(EVEC(NODS+IEVAL+1,JEVEC)),
C     &AIMAG(EVEC(NODS+IEVAL+1,JEVEC))
  605 CONTINUE
      END IF
c  700 FORMAT(2X,3Hfr(,I3,1H,I3,2H)=,E15.6,2X,3Hgr(,I3,1H,I3,2H)=,E15.6)
c  800 FORMAT(2X,3Hft(,I3,1H,I3,2H)=,E15.6,2X,3Hgt(,I3,1H,I3,2H)=,
c     &E15.6/)
  700 FORMAT(2X,E15.6,5X,E15.6)
  800 FORMAT(2X,E15.6,5X,E15.6)
      STOP
    END