C     PROGRAM FOR COMPUTING EIGENVALUES AND (LEFT & RIGHT) EIGENVECTORS
C     OF A SYMMETRIC OR WELL-CONDITIONED ASYMMETRIC REAL (SQUARE) MATRIX
C     FOR ILL-CONDITIONED (UNSTABLE) ASYMMETRIC MATRICES OR THE MATRICES
C     WITH REPEATED (OR VERY CLOSELY SPACED) EIGENVALUES USE SOME
C     OTHER METHODS (VIZ. SCHUR's DECOMPOSITION, JORDAN CANONICAL FORM
C     ETCETERA). THE LEFT EIGENVECTOR IS SOLUTION OF YA=YL AND THE RIGHT
C     EIGENVECTOR IS THE SOLUTION OF AX=LX. FOR SYMMETRIC MATRICES WE
C     HAVE Y=X WHILE L (EIGENVALUES) ARE COMMON.
C     -----------------------------------------------------------------
C     METHOD:POWER METHOD IN CONJUNCTION WITH EXTENDED DEFLATION METHOD
C     BY WILSON,EB, DECIUS,JC AND CROSS,PC (1955) MOLECULAR VIBRATIONS.
C     MCGRAW-HILL, NEW YORK
C     -----------------------------------------------------------------
C     PROGRAM BY SK MISHRA, DEPT. OF ECONOMICS, NEHU, SHILLONG (INDIA)
C     -----------------------------------------------------------------
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      PARAMETER (NH=100) ! LARGEST DIMENSION
      PARAMETER (QMAX0=1.0D50, EPS0=1.0D-08) ! EPS0 IS ACCURACY
      DIMENSION A(NH,NH),PR(NH),PL(NH),QR(NH),QL(NH)
      DIMENSION C(NH,NH),AL(NH),B(NH,NH)
      CHARACTER *30 IFIL, OFIL
C     -----------------------------------------------------------------
      WRITE(*,*)'PROGRAM TO COMPUTE EIGEN-VALUES AND EIGEN-VECTORS OF A'
      WRITE(*,*)'REAL SYMMETRIC MATRIX OR (LEFT & RIGHT) EIGEN-VECTORS'
      WRITE(*,*)'OF A WELL-CONDITIONED ASYMMETRIC REAL (SQUARE) MATRIX'
      WRITE(*,*)'------------------------------------------------------'
      WRITE(*,*)' '
C     INPUT FILE NAME WITH DATA : N, A(N,N), AND OUTPUT FILE
      WRITE(*,*)'FILE NAMES OF INPUT MATRIX AND OUTPUT ?'
      READ(*,*) IFIL,OFIL
      OPEN(7,FILE=IFIL)
      READ(7,*) N
      DO I=1,N
      READ(7,*)(A(I,J),J=1,N)
      ENDDO
      CLOSE(7)
C     COMPUTATION BEGINS
      TRACE=0.D0
      DO I=1,N
      TRACE=TRACE+A(I,I)
      ENDDO
      WRITE(*,*)'NO. OF ITERATIONS TO PERFORM, SAY 100 OR MORE'
      READ(*,*) NIT
      OPEN(7,FILE=OFIL)
C     INITIALIZE B BATRIX
      DO I=1,N
      DO J=1,N
      B(I,J)=0.D0
      ENDDO
      ENDDO

      DO K=1,N
C     INITIALIZE PR AND PL AS UNIT VECTORS
      DO I=1,N
      PR(I)=1.D0
      PL(I)=1.D0
      ENDDO
C     RIGHT MULTIPLY AND LEFT MULTIPLY A BY PR AND PL - GET QR, QL
      QMAXR=QMAX0 ! INITIALIZE BY VERY LARGE FIGURE
      QMAXL=QMAX0 ! INITIALIZE BY VERY LARGE FIGURE
      EPS=EPS0  ! A SMALL FIGURE
   1  DO KK=1,NIT ! NO. OF ITERATIONS TO PERFORM
      DO I=1,N
        QR(I)=0.D0
        QL(I)=0.D0
        DO J=1,N
          QR(I)=QR(I)+A(I,J)*PR(J)
        ENDDO
        DO J=1,N
          QL(I)=QL(I)+PL(J)*A(J,I)
        ENDDO
      ENDDO
C     FIND MAX MAGNITUDE ELEMENTS OF QR AND QL, PRESERVING THEIR SIGNS
      QRMAX=DABS(QR(1))
      SIGR=ANINT(QR(1)/QRMAX)
      QLMAX=DABS(QL(1))
      SIGL=ANINT(QL(1)/QLMAX)
      DO J=2,N
      IF(QRMAX.LT.DABS(QR(J))) THEN
      QRMAX=DABS(QR(J))
      SIGR=ANINT(QR(J)/QRMAX)
      ENDIF
      IF(QLMAX.LT.DABS(QL(J))) THEN
      QLMAX=DABS(QL(J))
      SIGL=ANINT(QL(J)/QLMAX)
      ENDIF
      ENDDO
C     RESUME THE SIGN OF QRMAX AND QLMAX
      QRMAX=QRMAX*SIGR
      QLMAX=QLMAX*SIGL
C     DIVIDE QR AND QL BY THEIR RESPECTIVE QRMAX AND QLMAX AND STORE
C     THEM IN PR AND PL RESPECTIVELY

      DO J=1,N
      PR(J)=QR(J)/QRMAX
      PL(J)=QL(J)/QLMAX
      ENDDO


C     COMPARE QRMAX AND QLMAX (CONSIDERING SIGNS AS WELL)
      IF(QRMAX-QLMAX.GT.EPS.OR.(DABS(DABS(QRMAX)-QMAXR)).GT.EPS.OR.
     & (DABS(DABS(QLMAX)-QMAXL)).GT.EPS) THEN
      QMAXR=QRMAX
      QMAXL=QLMAX
      ENDIF
      ENDDO

C     STANDARDIZE PR AND PL TO HAVE UNITY NORM
      SL=0.D0
      SR=0.D0
      DO I=1,N
      SL=SL+PL(I)**2
      SR=SR+PR(I)**2
      ENDDO

      DO I=1,N
      PR(I)=PR(I)/SR
      PL(I)=PL(I)/SL
      ENDDO
      WRITE(7,*)'EIGENVALUE= ',QRMAX
      WRITE(7,*)'RIGHT EIGENVECTOR '
      WRITE(7,*)(PR(I),I=1,N)
      WRITE(7,*)'LEFT EIGENVECTOR '
      WRITE(7,*)(PL(I),I=1,N)
C     STORE EIGENVALUES IN AL
      AL(K)=QRMAX
C     CONSTRUCT C AS PR X QRMAX PL
      SRL=0.D0
      DO I=1,N
      SRL=SRL+PL(I)*PR(I)
      ENDDO
      DO I=1,N
      DO J=1,N
      C(I,J)=PR(I)*QRMAX*PL(J)/SRL

      ENDDO
      ENDDO
      DO I=1,N
      DO J=1,N
      B(I,J)=B(I,J)+C(I,J) ! RECONSTRUCT THE ORIGINAL A MATRIX
      A(I,J)=A(I,J)-C(I,J)
      ENDDO
      ENDDO
      WRITE(7,*)'-------- RESIDUAL MATRIX AT THIS STAGE --------------'
      DO I=1,N
      WRITE(7,10)(A(I,J),J=1,N)
      ENDDO
      ENDDO
  10  FORMAT(8F10.3)

      SL=0.D0
      DET=1.0D0
      DO I=1,N
      SL=SL+AL(I)
      DET=DET*AL(I)
      ENDDO
      WRITE(7,*)'TRACE OF MATRIX, SUM OF EIGENVALUES AND DETERMINANT='
      WRITE(7,*) TRACE, SL,DET
      WRITE(7,*)'EIGENVALUES ARE :'
      WRITE(7,*)(AL(K),K=1,N)
      CLOSE(7)
      WRITE(*,*)'END'
      END