!*************************************************************************
!*   Eigenvalues and Eigenvectors of a general complex square matrix     *
!*               using the QR method for complex matrices                *
!* --------------------------------------------------------------------- *
!* SAMPLE RUN:                                                           *
!*                                                                       *
!* Example #1:                                                           *
!*                                                                       *
!* Input matrix (real symmetric for comparison with Jacobi):             *
!*        1,  2,  3, -7,   12                                            *
!*        2,  4,  7,  3,   -1                                            *
!*        3,  7, 10,  8,    4                                            *
!*       -7,  3,  8, -0.75,-9                                            *
!*       12, -1,  4, -9,   10                                            *
!*                                                                       *
!*  Input matrix (real part):                                            *
!*1   1.00000000   2.00000000   3.00000000  -7.00000000  12.00000000     *
!*2   2.00000000   4.00000000   7.00000000   3.00000000  -1.00000000     *
!*3   3.00000000   7.00000000  10.00000000   8.00000000   4.00000000     *
!*4  -7.00000000   3.00000000   8.00000000  -0.75000000  -9.00000000     *
!*5  12.00000000  -1.00000000   4.00000000  -9.00000000  10.00000000     *
!*                                                                       *
!*  Error code =           0                                             *
!*                                                                       *
!*  Eigebvalues:                                                         *
!*  0.463349527981464       0.000000000000000E+000                       *
!*  -7.77457972947830       0.000000000000000E+000                       *
!*  -10.4865451655156       0.000000000000000E+000                       *
!*   18.2918206016854       0.000000000000000E+000                       *
!*   23.7559547653270       0.000000000000000E+000                       *
!*                                                                       *
!*  Eigenvectors (real part in lines):                                   *
!*1   0.20881245   1.00000000  -0.47802717  -0.27500001  -0.21691482     *
!*2   1.00000000  -0.32610016   0.20961973  -0.14768922  -0.81542193     *
!*3   0.46717183  -0.00794702  -0.50782668   1.00000000   0.26443199     *
!*4   0.09305260   0.59491100   1.00000000   0.45566550   0.05074045     *
!*5   0.70581969  -0.01267740   0.13148643  -0.52749949   1.00000000     *
!*                                                                       *
!* Example #2:                                                           *
!*                                                                       *
!* Input matrix (real part):                                             *
!*        1,  5,  3, -7,   12                                            *
!*        2,  4,  7,  3,   -1                                            *
!*        3,  7, 10,  8,    4                                            *
!*       -7,  3,  8, -0.75,-9                                            *
!*       12, -1,  4, -9,   10                                            *
!* Input matrix (imaginary part):                                        *
!*        1,  1,  1,  1,    1                                            *
!*        1,  1,  1,  1,    1                                            *
!*        1,  1,  1,  1,    1                                            *
!*        1,  1,  1,  1     1                                            *
!*        1,  1,  1,  1,    1                                            *
!*                                                                       *
!*  Input Matrix (real  part):                                           *
!*1   1.00000000   5.00000000   3.00000000  -7.00000000  12.00000000     * 
!*2   2.00000000   4.00000000   7.00000000   3.00000000  -1.00000000     *
!*3   3.00000000   7.00000000  10.00000000   8.00000000   4.00000000     *
!*4  -7.00000000   3.00000000   8.00000000  -0.75000000  -9.00000000     *
!*5  12.00000000  -1.00000000   4.00000000  -9.00000000  10.00000000     *
!*  Input Matrix (imag. part):                                           *
!*1   1.00000000   1.00000000   1.00000000   1.00000000   1.00000000     *
!*2   1.00000000   1.00000000   1.00000000   1.00000000   1.00000000     *
!*3   1.00000000   1.00000000   1.00000000   1.00000000   1.00000000     *
!*4   1.00000000   1.00000000   1.00000000   1.00000000   1.00000000     *
!*5   1.00000000   1.00000000   1.00000000   1.00000000   1.00000000     *
!*                                                                       *
!*  Error code =           0                                             *
!*                                                                       *
!*  Eigenvalues:                                                         *
!*  0.888713348723694       3.470700022098407E-002                       *
!*  -8.28426842306568       2.498162361886443E-002                       *
!*  -10.3636733606943       0.879708207660510                            *
!*   18.8947919026401        3.15558227567471                            *
!*   23.1144365323962       0.905020892824929                            *
!*                                                                       *
!*  Eigenvectors (real  part in lines):                                  *
!*1   0.41150176   1.00000000  -0.53250763  -0.21794537  -0.41700119     *
!*2   1.00000000  -0.32418914   0.26914039  -0.27758700  -0.87269351     *
!*3   0.45123301  -0.00728949  -0.51526552   1.00000000   0.27273187     *
!*4   0.17988445   0.60859892   1.00000000   0.45640602   0.06919679     *
!*5   0.68265703  -0.10666293  -0.00445564  -0.61516974   1.00000000     *
!*  Eigenvectors (imag. part in lines):                                  *
!*1   0.00076996   0.00000000  -0.03474965   0.00906216  -0.00518060     *
!*2   0.00000000   0.00167602   0.05162180  -0.07411949  -0.03764759     *
!*3   0.04349160  -0.07014688  -0.05051509   0.00000000  -0.06636677     *
!*4  -0.19851484   0.07550215   0.00000000   0.27102952  -0.31492765     *
!*5   0.04407051   0.19380611   0.29215881   0.17278034   0.00000000     *
!*                                                                       *           
!* --------------------------------------------------------------------- *
!* Reference:   From Numath Library By Tuan Dang Trong in Fortran 77     *
!*              [BIBLI 18].                                              *
!*                                                                       *
!*                               F90 Release 1.1 By J-P Moreau, Paris    *
!*                                        (www.jpmoreau.fr)              *
!* --------------------------------------------------------------------- *
!* Release 1.1: added subroutine MTPRINT and Example #2 (12/04/2004).    *
!*************************************************************************
PROGRAM TEST_CEIGEN
parameter(NM=5,N=5,LMA=5)
  REAL*8 AR(NM,N), AI(NM,N),ZR(LMA,N),ZI(LMA,N)
  REAL*8 WR(N),WI(N),WORK(3*N)
  LOGICAL AVEC
  INTEGER IERR

! For #1:
  DATA AR /1.d0,2.d0,3.d0,-7.d0,12.d0,    &
           2.d0,4.d0,7.d0,3.d0,-1.d0,     &
		   3.d0,7.d0,10.d0,8.d0,4.d0,     &
		   -7.d0,3.d0,8.d0,-0.75d0,-9.d0, &
		   12.d0,-1.d0,4.d0,-9.d0,10.d0/

  AI = 0.d0

! For #2:
! DATA AR /1.d0,2.d0,3.d0,-7.d0,12.d0,    &
!          5.d0,4.d0,7.d0,3.d0,-1.d0,     &
!		   3.d0,7.d0,10.d0,8.d0,4.d0,     &
!		   -7.d0,3.d0,8.d0,-0.75d0,-9.d0, &
!		   12.d0,-1.d0,4.d0,-9.d0,10.d0/

! AI = 1.d0

  CALL MPRINT(N,NM,N,'  Input Matrix (real part):',AR)
! #2 only:
! CALL MPRINT(N,NM,N,'  Input Matrix (imag. part):',AI)
! pause

  CALL CEIGEN(AR,AI,NM,N,.TRUE.,WR,WI,ZR,ZI,WORK,IERR)

  print *,' '
  print *,' Error code =', IERR
  print *,' '
  print *,' Eigenvalues:'
  do i=1,N
    print *, WR(i), WI(i)
  end do
  print *,' '
  CALL MTPRINT(N,NM,N,'  Eigenvectors (real part in lines):',ZR)
! #2 only:
! CALL MTPRINT(N,NM,N,'  Eigenvectors (imag. part in lines):',ZI)
  print *,' '

END !of main program

      !print a matrix with caption
      SUBROUTINE MPRINT(N,MMAX,M,LABEL,A)
      REAL *8 A(MMAX,1)
	  CHARACTER*(*) LABEL
      PRINT 2,LABEL
      DO 1 I=1,N
    1 PRINT 3,I,(A(I,J),J=1,M)
      RETURN
    2 FORMAT(A)
    3 FORMAT(I1,5(F13.8))
      END                 

      !print transpose matrix with caption
      SUBROUTINE MTPRINT(N,MMAX,M,LABEL,A)
      REAL *8 A(MMAX,1)
	  CHARACTER*(*) LABEL
      PRINT 2,LABEL
      DO 1 I=1,N
    1 PRINT 3,I,(A(J,I),J=1,M)
      RETURN
    2 FORMAT(A)
    3 FORMAT(I1,5(F13.8))
      END         

      SUBROUTINE CEIGEN(AR,AI,NM,N,AVEC,WR,WI,ZR,ZI,WORK,IERR)
! --------------------------------------------------------------------
!     This subroutine calculates the eigenvalues and eigenvectors
!     of a general complex matrix. The used method begins with a
!     reduction by unitary transformations to a Hessenberg form.
!     The QR algorithm is then applied to get the eigenvalues.
!     Finally, the corresponding eigenvectors are found by a simple
!     substitution.
! --------------------------------------------------------------------
!     CALLING MODE:
!     CALL CEIGEN (AR,AI,LMA,N,AVEC,WR,WI,ZR,ZI,WORK,IERR)

!     INPUTS:
!     AR   =   TABLE(LMA,N) OF REAL COEF. OF MATRIX A              R*8
!     AI   =   TABLE(LMA,N) OF IMAGINARY COEF. OF MATRIX A         R*8
!     NM   =   1ST DIMENSION OF MATRICES AR,AI,ZR & ZI
!              AS DECLARED IN CALLING PROGRAM                      I*4
!     N    =   2ND DIMENSION OF MATRICES AR,AI,ZR & ZI             I*4
!     AVEC =   CODE TO CALCULATE OR NOT EIGENVECTORS               L*8
!          =   .FALSE. NO EIGENVECTORS
!          =   .TRUE.  WITH EIGENVECTORS

!     OUTPUTS:
!     WR,WI=   TABLEA(N) OF EIGENVALUES                            R*8
!     ZR,ZI=   TABLE(LMA,N) OF EIGENVECTORS                        R*8
!     IERR =   ERROR CODE                                          I*4
!          =   0 NO ERROR
!          < OU = N  NO CONVERGENCE
!          =   9999  DIMENSION ERROR IN TABLES A OR ZR,ZI

!     WORKING SPACE:
!     WORK =  TABLE(3*N)                                           R*8

!     REFERENCE:
!     J.H.WILKINSON & C.REINSCH
!     LINEAR ALGEBRA   SPRINGER-VERLAG 1971
! --------------------------------------------------------------------
      REAL*8 AR(NM,*),AI(NM,*),ZR(NM,*),ZI(NM,*),WR(*),WI(*),WORK(*)
      INTEGER UPP
      LOGICAL AVEC
      IF(N.LE.NM) THEN
        N1=N+1
        N2=N+N1
        CALL CBALAN(NM,N,AR,AI,LOW,UPP,WORK(1))
        CALL CUNITH(NM,N,LOW,UPP,AR,AI,WORK(N1),WORK(N2))
        CALL COMQR2(NM,N,LOW,UPP,WORK(N1),WORK(N2),AR,AI,WR,WI,  &
                    AVEC,ZR,ZI,IERR)
        IF(IERR.EQ.0) THEN
        CALL CBALBK(NM,N,LOW,UPP,N,WORK(1),ZR,ZI)
        CALL NORMAL(NM,N,WR,WI,AVEC,ZR,ZI,WORK(N1),WORK(N2))
        ENDIF
      ELSE
      IERR= 9999
      ENDIF
      RETURN
      END
!
      SUBROUTINE CUNITH(NM,N,LOW,UPP,AR,AI,DR,DI)

!     UNITARY TRANSFORMATIONS TO REDUCE A COMPLEX MATRIX
!     TO A HESSENBERG FORM
      REAL*8 AR(NM,*),AI(NM,*),DR(*),DI(*),XR,XI,G,F,H,S

      INTEGER UPP
      IF(UPP-1.LT.LOW+1) RETURN
      DO M=LOW+1,UPP-1
         DR(M)=0.D0
         DI(M)=0.D0
         S=0.D0
         H=0.D0
         DO I=M,UPP
            S=S+ABS(AR(I,M-1))+ABS(AI(I,M-1))
         ENDDO
!        WRITE(*,'(1X,A,I4,E15.6)') ' M,S  ',M,S
         IF(S.NE.0.D0) THEN
         DO I=UPP,M,-1
            DR(I)=AR(I,M-1)/S
            DI(I)=AI(I,M-1)/S
            H=H+DR(I)*DR(I)+DI(I)*DI(I)
         ENDDO
         G=SQRT(H)
         F=ABSCD(DR(M),DI(M))
         IF(F.EQ.0.D0) THEN
         DR(M)=G
         AR(M,M-1)=S
         ELSE
         H=H+F*G
         G=G/F
         DR(M)=(1.D0+G)*DR(M)
         DI(M)=(1.D0+G)*DI(M)
         ENDIF
         DO J=M,N
            XR=0.D0
            XI=0.D0
            DO I=UPP,M,-1
               XR=XR+DR(I)*AR(I,J)+DI(I)*AI(I,J)
               XI=XI+DR(I)*AI(I,J)-DI(I)*AR(I,J)
            ENDDO
            XR=XR/H
            XI=XI/H
            DO I=M,UPP
               AR(I,J)=AR(I,J)-XR*DR(I)+XI*DI(I)
               AI(I,J)=AI(I,J)-XR*DI(I)-XI*DR(I)
            ENDDO
         ENDDO
         DO I=1,UPP
            XR=0.D0
            XI=0.D0
            DO J=UPP,M,-1
               XR=XR+DR(J)*AR(I,J)-DI(J)*AI(I,J)
               XI=XI+DR(J)*AI(I,J)+DI(J)*AR(I,J)
            ENDDO
            XR=XR/H
            XI=XI/H
            DO J=M,UPP
               AR(I,J)=AR(I,J)-XR*DR(J)-XI*DI(J)
               AI(I,J)=AI(I,J)+XR*DI(J)-XI*DR(J)
            ENDDO
         ENDDO
         DR(M)=DR(M)*S
         DI(M)=DI(M)*S
         AR(M,M-1)=-G*AR(M,M-1)
         AI(M,M-1)=-G*AI(M,M-1)
         ENDIF
      ENDDO
      RETURN
      END

      SUBROUTINE COMQR2(NM,N,LOW,UPP,GR,GI,HR,HI,WR,WI,AVEC,ZR,ZI,IERR)
! ---------------------------------------------------------------------
!     QR ALGORITH APPLIED TO A COMPLEX HESSENBERG MATRIX
! ---------------------------------------------------------------------
      REAL*8 HR(NM,*),HI(NM,*),ZR(NM,*),ZI(NM,*),WR(*),WI(*),GR(*), &
             GI(*),SR,SI,TR,TI,YR,YI,XR,XI,VR,VI,MACHEP,NORM,EPS,EPS1
      INTEGER I,J,K,L,M,N,IN,NA,NM,ITS,LOW,UPP,IERR
      LOGICAL AVEC

!     SEEK FLOATING POINT RELATIVE PRECISION FOR TYPE REAL*8
      EPS=1.D0
    5 EPS=EPS*0.5D0
      EPS1=1.D0+EPS
      IF(EPS1.GT.1.D0) GO TO 5
      MACHEP=2.D0*EPS

      IERR=0
      IF(AVEC) THEN
      DO I=1,N
         DO J=1,N
            ZR(I,J)=0.D0
            ZI(I,J)=0.D0
         ENDDO
         ZR(I,I)=1.D0
      ENDDO

!     ACCUMULATE UNITARY TRANSFORMATIONS
      DO I=UPP-1,LOW+1,-1
         IF((HR(I,I-1).EQ.0.D0.AND.HI(I,I-1).EQ.0.D0)  &
         .OR.(GR(I).EQ.0.D0.AND.GI(I).EQ.0.D0)) GO TO 7
         NORM=HR(I,I-1)*GR(I)+HI(I,I-1)*GI(I)
         DO K=I+1,UPP
            GR(K)=HR(K,I-1)
            GI(K)=HI(K,I-1)
         ENDDO
         DO J=I,UPP
            SR=0.D0
            SI=0.D0
            DO K=I,UPP
               SR=SR+GR(K)*ZR(K,J)+GI(K)*ZI(K,J)
               SI=SI+GR(K)*ZI(K,J)-GI(K)*ZR(K,J)
            ENDDO
            SR=SR/NORM
            SI=SI/NORM
            DO K=I,UPP
               ZR(K,J)=ZR(K,J)+SR*GR(K)-SI*GI(K)
               ZI(K,J)=ZI(K,J)+SR*GI(K)+SI*GR(K)
            ENDDO
         ENDDO
    7 CONTINUE
      ENDDO
      ENDIF

!     CREATE REAL SUBDIAGONAL
      L=LOW+1
      DO I=L,UPP
         M=MIN(I+1,UPP)
         IF(HI(I,I-1).EQ.0.D0) GO TO 9
         NORM=ABSCD(HR(I,I-1),HI(I,I-1))
         YR=HR(I,I-1)/NORM
         YI=HI(I,I-1)/NORM
         HR(I,I-1)=NORM
         HI(I,I-1)=0.D0
         DO J=I,N
            SI=     YR*HI(I,J)-YI*HR(I,J)
            HR(I,J)=YR*HR(I,J)+YI*HI(I,J)
            HI(I,J)=SI
         ENDDO
         DO J=1,M
            SI=     YR*HI(J,I)+YI*HR(J,I)
            HR(J,I)=YR*HR(J,I)-YI*HI(J,I)
            HI(J,I)=SI
         ENDDO
         IF(AVEC) THEN
         DO J=LOW,UPP
            SI=     YR*ZI(J,I)+YI*ZR(J,I)
            ZR(J,I)=YR*ZR(J,I)-YI*ZI(J,I)
            ZI(J,I)=SI
         ENDDO
         ENDIF
    9 CONTINUE
      ENDDO

!     STORE ROOTS ALREADY ISOLATED BY SUBROUTINE BALANC
      DO I=1,N
         IF(I.LT.LOW.OR.I.GT.UPP) THEN
         WR(I)=HR(I,I)
         WI(I)=HI(I,I)
         ENDIF
      ENDDO

      IN=UPP
      TR=0.D0
      TI=0.D0

!     SEEK EIGENVALUES
   10 IF(IN.LT.LOW) GO TO 50
      ITS=0
      NA=IN-1

!     WE LOOK AT REAL PART OF SUBDIAGONAL ELEMENT
   15 DO L=IN,LOW,-1
         SR=MACHEP*(ABS(HR(L-1,L-1))+ABS(HI(L-1,L-1))  &
                    +ABS(HR(L,L))+ABS(HI(L,L)))
         IF (ABS(HR(L,L-1)).LE.SR.OR.L.EQ.LOW) GO TO 20
      ENDDO

!     SHIFTING
   20 IF(L.EQ.IN) GO TO 45
      IF(ITS.EQ.30) GO TO 55
      IF(ITS.EQ.10.OR.ITS.EQ.20) THEN

!     EXTRA SHIFTING
         SR=ABS(HR(IN,NA))+ABS(HR(NA,IN-2))
         SI=0.D0
      ELSE
         SR=HR(IN,IN)
         SI=HI(IN,IN)
         XR=HR(NA,IN)*HR(IN,NA)
         XI=HI(NA,IN)*HR(IN,NA)

         IF(XR.NE.0.D0.OR.XI.NE.0.D0) THEN
            YR=(HR(NA,NA)-SR)*.5D0
            YI=(HI(NA,NA)-SI)*.5D0
            CALL DCSQRT(YR*YR-YI*YI+XR,2.D0*YR*YI+XI,VR,VI)
            IF(YR*VR+YI*VI.LT.0.D0) THEN
               VR=-VR
               VI=-VI
            ENDIF
            CALL DIVCD(XR,XI,YR+VR,YI+VI,VR,VI)
            SR=SR-VR
            SI=SI-VI
         ENDIF
      ENDIF

      DO I=LOW,IN
         HR(I,I)=HR(I,I)-SR
         HI(I,I)=HI(I,I)-SI
      ENDDO

      TR=TR+SR
      TI=TI+SI
      ITS=ITS+1

!     TRIANGULARY DECOMPOSITION H = Q * R
      DO I=L+1,IN
         SR=HR(I,I-1)
         VR=HR(I-1,I-1)
         VI=HI(I-1,I-1)
         NORM=SQRT(VR*VR+VI*VI+SR*SR)
         XR=VR/NORM
         XI=VI/NORM
         WR(I-1)=XR
         WI(I-1)=XI
         HR(I-1,I-1)=NORM
         HI(I-1,I-1)=0.D0
         HR(I,I-1)=0.D0
         HI(I,I-1)=SR/NORM
         DO J=I,N
            YR=HR(I-1,J)
            YI=HI(I-1,J)
            VR=HR(I,J)
            VI=HI(I,J)
            HR(I-1,J)=XR*YR+XI*YI+HI(I,I-1)*VR
            HI(I-1,J)=XR*YI-XI*YR+HI(I,I-1)*VI
            HR(I,J)  =XR*VR-XI*VI-HI(I,I-1)*YR
            HI(I,J)  =XR*VI+XI*VR-HI(I,I-1)*YI
         ENDDO
      ENDDO
      SI=HI(IN,IN)

      IF(SI.EQ.0.D0) GO TO 35
      NORM=ABSCD(HR(IN,IN),SI)
      SR=HR(IN,IN)/NORM
      SI=SI/NORM
      HI(IN,IN)=0.D0
      HR(IN,IN)=NORM
      IF(IN.NE.N) THEN
      DO J=IN+1,N
         YR=HR(IN,J)
         YI=HI(IN,J)
         HR(IN,J)=SR*YR+SI*YI
         HI(IN,J)=SR*YI-SI*YR
      ENDDO
      ENDIF

!     INVERSE OPERATION
   35 DO J=L+1,IN
         XR=WR(J-1)
         XI=WI(J-1)
         DO I=1,J
            YR=HR(I,J-1)
            YI=0.D0
            VR=HR(I,J)
            VI=HI(I,J)
            IF(I.NE.J) THEN
            YI=HI(I,J-1)
            HI(I,J-1)=XR*YI+XI*YR+HI(J,J-1)*VI
            ENDIF
            HR(I,J-1)=XR*YR-XI*YI+HI(J,J-1)*VR
            HR(I,J)  =XR*VR+XI*VI-HI(J,J-1)*YR
            HI(I,J)  =XR*VI-XI*VR-HI(J,J-1)*YI
         ENDDO
         IF(AVEC) THEN
         DO I=LOW,UPP
            YR=ZR(I,J-1)
            YI=ZI(I,J-1)
            VR=ZR(I,J)
            VI=ZI(I,J)
            ZR(I,J-1)=XR*YR-XI*YI+HI(J,J-1)*VR
            ZI(I,J-1)=XR*YI+XI*YR+HI(J,J-1)*VI
            ZR(I,J)  =XR*VR+XI*VI-HI(J,J-1)*YR
            ZI(I,J)  =XR*VI-XI*VR-HI(J,J-1)*YI
         ENDDO
         ENDIF
      ENDDO

      IF(SI.EQ.0.D0) GO TO 15
      DO I=1,IN
         YR=HR(I,IN)
         YI=HI(I,IN)
         HR(I,IN)=SR*YR-SI*YI
         HI(I,IN)=SR*YI+SI*YR
      ENDDO
      IF(AVEC) THEN
      DO I=LOW,UPP
         YR=ZR(I,IN)
         YI=ZI(I,IN)
         ZR(I,IN)=SR*YR-SI*YI
         ZI(I,IN)=SR*YI+SI*YR
      ENDDO
      ENDIF
      GO TO 15

!     ONE ROOT HAS BEEN FOUND
   45 HR(IN,IN)=HR(IN,IN)+TR
      HI(IN,IN)=HI(IN,IN)+TI
      WR(IN)=HR(IN,IN)
      WI(IN)=HI(IN,IN)
      IN=NA
      GO TO 10

!     ALL ROOTS HAVE BEEN FOUND
!     START CALCULATION OF EIGENVECTORS BY SUBSTITUTION
   50 IF(AVEC) THEN
      NORM=0.D0
      DO I=1,N
         DO J=I,N
         NORM=NORM+ABS(HR(I,J))+ABS(HI(I,J))
         ENDDO
      ENDDO
      IF(NORM.EQ.0.D0) GO TO 60
      DO IN=N,2,-1
         XR=WR(IN)
         XI=WI(IN)
         HR(IN,IN)=1.D0
         HI(IN,IN)=0.D0
         DO I=IN-1,1,-1
            VR=HR(I,IN)
            VI=HI(I,IN)
            IF(I.NE.IN-1) THEN
            DO J=I+1,IN-1
               VR=VR+HR(I,J)*HR(J,IN)-HI(I,J)*HI(J,IN)
               VI=VI+HR(I,J)*HI(J,IN)+HI(I,J)*HR(J,IN)
            ENDDO
            ENDIF
            YR=XR-WR(I)
            YI=XI-WI(I)
            IF(YR.EQ.0.D0.AND.YI.EQ.0.D0) YR=NORM*MACHEP
            CALL DIVCD(VR,VI,YR,YI,TR,TI)
            HR(I,IN)=TR
            HI(I,IN)=TI
         ENDDO
      ENDDO
      DO I=1,N-1
         IF(I.LT.LOW.OR.I.GT.UPP) THEN
         DO J=I+1,N
            ZR(I,J)=HR(I,J)
            ZI(I,J)=HI(I,J)
         ENDDO
         ENDIF
      ENDDO

!     END CALCULATION OH HESSENBERG EIGENVECTORS
      DO J=N,LOW+1,-1
         M=MIN(J-1,UPP)
         DO I=LOW,UPP
            VR=ZR(I,J)
            VI=ZI(I,J)
            DO K=LOW,M
               VR=VR+ZR(I,K)*HR(K,J)-ZI(I,K)*HI(K,J)
               VI=VI+ZR(I,K)*HI(K,J)+ZI(I,K)*HR(K,J)
            ENDDO
            ZR(I,J)=VR
            ZI(I,J)=VI
         ENDDO
      ENDDO
      ENDIF

!     END CALCULATION OF EIGENVECTORS OF INPUT MATRIX
      GO TO 60

!     NO CONVERGENCE AFTER 30 ITERATIONS!
   55 IERR = IN
   60 RETURN
      END

      SUBROUTINE CBALAN(NM,N,AR,AI,LOW,UPP,D)

!     BALANCE A COMPLEX MATRIX (VERSION DERIVED FROM BALANC)
      REAL*8 AR(NM,N),AI(NM,N),D(N),B2,C,F,G,R,S
      LOGICAL NOCONV
      INTEGER UPP,B
      DATA B/16/
      B2=B*B
      L=1
      K=N
   10 DO J=K,1,-1
         J1=J
         K1=K
         L1=L
         DO I=1,K
            IF(I.NE.J) THEN
            IF(AR(J,I).NE.0.D0 .OR. AI(J,I).NE.0.D0) GO TO 11
            ENDIF
         ENDDO
         CALL EXCHG1(K,NM,N,AR,AI,D,J1,K1,L1)
         K=K-1
         GO TO 10
   11 CONTINUE
      ENDDO
   20 DO J=L,K
         J1=J
         K1=K
         L1=L
         DO I=L,K
            IF(I.NE.J) THEN
            IF(AR(I,J).NE.0.D0 .OR. AI(I,J).NE.0.D0) GO TO 22
            ENDIF
         ENDDO

         CALL EXCHG1(L,NM,N,AR,AI,D,J1,K1,L1)
         L=L+1
         GO TO 20
   22 CONTINUE
      ENDDO
      LOW=L
      UPP=K
      DO I=L,K
         D(I)=1.D0
      ENDDO
   30 NOCONV=.FALSE.
      DO I=L,K
         C=0.D0
         R=0.D0
         DO J=L,K
            IF(J.NE.I) THEN
            C=C+ABS(AR(J,I))+ABS(AI(J,I))
            R=R+ABS(AR(I,J))+ABS(AI(I,J))
            ENDIF
         ENDDO
         G=R/B
         F=1.D0
         S=C+R
   40    IF(C.LT.G) THEN
         F=F*B
         C=C*B2
         GO TO 40
         ENDIF
         G=R*B
   50    IF(C.GE.G) THEN
         F=F/B
         C=C/B2
         GO TO 50
         ENDIF
         IF((C+R)/F.LT.0.95D0*S) THEN
         G=1.D0/F
         D(I)=D(I)*F
         NOCONV=.TRUE.
         DO J=L,N
            AR(I,J)=AR(I,J)*G
            AI(I,J)=AI(I,J)*G
         ENDDO
         DO J=1,K
            AI(J,I)=AI(J,I)*F
            AR(J,I)=AR(J,I)*F
         ENDDO
      ENDIF
      ENDDO
      IF(NOCONV) GO TO 30
      RETURN
      END

      SUBROUTINE CBALBK(NM,N,LOW,UPP,M,D,Z,T)
!     This subroutine finds the eigenvectors of a complex matrix,
!     as output of subroutine CBALAN.
      REAL*8 Z(NM,*),T(NM,*),D(*),S
      INTEGER UPP
      DO I=LOW,UPP
      S=D(I)
         DO J=1,M
         Z(I,J)=Z(I,J)*S
         T(I,J)=T(I,J)*S
         ENDDO
      ENDDO
      DO L=1,N
      I=L
      IF(I.GE.LOW.AND.I.LE.UPP) GO TO 10
      IF(I.LT.LOW) I=LOW-L
      K=D(I)
         IF(K.NE.I) THEN
         DO J=1,M
         S=Z(I,J)
         Z(I,J)=Z(K,J)
         Z(K,J)=S
         S=T(I,J)
         T(I,J)=T(K,J)
         T(K,J)=S
         ENDDO
         ENDIF
   10 CONTINUE
      ENDDO
      RETURN
      END

      SUBROUTINE NORMAL(NM,N,WR,WI,AVEC,ZR,ZI,TR,TI)
!     Sort eigenvalues by absolute values in ascending order and
!     normalize.
      REAL*8 ZR(NM,*),ZI(NM,*),WR(*),WI(*),TR(*),TI(*),ZM,Z1,Z2,VR,VI, &
      ABSCD
      LOGICAL AVEC

!     SORT SOLUTIONS IN ASCENDING ORDER
      DO J=2,N
         VR=WR(J)
         VI=WI(J)
         IF(AVEC) THEN
         DO K=1,N
            TR(K)=ZR(K,J)
            TI(K)=ZI(K,J)
         ENDDO
         ENDIF
         DO I=J-1,1,-1
            IF( ABSCD(WR(I),WI(I)).LE. ABSCD(VR,VI) ) GO TO 5
            WR(I+1)=WR(I)
            WI(I+1)=WI(I)
            IF(AVEC) THEN
            DO K=1,N
               ZR(K,I+1)=ZR(K,I)
               ZI(K,I+1)=ZI(K,I)
            ENDDO
            ENDIF
         ENDDO
         I=0
    5    WR(I+1)=VR
         WI(I+1)=VI
         IF(AVEC) THEN
         DO K=1,N
            ZR(K,I+1)=TR(K)
            ZI(K,I+1)=TI(K)
         ENDDO
         ENDIF
      ENDDO

!     NORMALIZE WITH RESPECT TO BIGGEST ELEMENT
      IF(AVEC) THEN
      DO J=N,1,-1
         ZM= 0.D0
         DO I=1,N
            VR=ZR(I,J)
            VI=ZI(I,J)
            Z1=ABSCD (VR,VI)
            IF(Z1.GE.ZM) THEN
            IM = I
            ZM = Z1
            ENDIF
         ENDDO
         Z1 = ZR(IM,J)
         Z2 = ZI(IM,J)
         DO I = 1,N
            CALL DIVCD(ZR(I,J),ZI(I,J),Z1,Z2,VR,VI)
            ZR(I,J) = VR
            ZI(I,J) = VI
         ENDDO
      ENDDO
      ENDIF
      END

! Utility routines
      FUNCTION ABSCD(AR,AI)
!     ABSOLUTE VALUE OF A COMPLEX NUMBER C=AR+I*AI
!     ABSCD,AR,AI REAL*8
!     ABSCD=SQRT(AR**2+AI**2)
      REAL*8 ABSCD,AR,AI,XR,XI,W
      XR=ABS(AR)
      XI=ABS(AI)
      IF(XR.LE.XI) THEN
      W=XR
      XR=XI
      XI=W
      ENDIF
      IF(XI.EQ.0.D0) THEN
      ABSCD=XR
      ELSE
      ABSCD=XR*SQRT(1.D0+(XI/XR)**2)
      ENDIF
      RETURN
      END

      SUBROUTINE DCSQRT(X,Y,A,B)

!     SQUARE ROOT OF A COMPLEX NUMBER  A+I*B = RAC(X+I*Y)

      REAL*8 X,Y,A,B,DCABS
      IF(X.EQ.0.D0.AND.Y.EQ.0.D0) THEN
      A=0.D0
      B=0.D0
      ELSE
      A=SQRT((ABS(X)+DCABS(X,Y))*0.5D0)
      IF(X.GE.0.D0) THEN
      B=Y/(A+A)
      ELSE
      IF(Y.LT.0.D0) THEN
      B=-A
      ELSE
      B=A
      A=Y/(B+B)
      ENDIF
      ENDIF
      ENDIF
      END

!     ABSOLUTE VALUE OF A COMPLEX NUMBER X+I*Y

      FUNCTION DCABS(XX,YY)
      REAL*8 DCABS,XX,YY,X,Y,W
      X=ABS(XX)
      Y=ABS(YY)
      IF(X.EQ.0.D0) THEN
      W=Y
      ELSE
      IF(Y.EQ.0) THEN
      W=X
      ELSE
      IF(X.GT.Y) THEN
      W=X*SQRT(1.D0+(Y/X)**2)
      ELSE
      W=Y*SQRT(1.D0+(X/Y)**2)
      ENDIF
      ENDIF
      ENDIF
      DCABS=W
      RETURN
      END

      SUBROUTINE DIVCD(AR,AI,BR,BI,ZR,ZI)

!     EFFECT DIVISION Z=(ZR+I*ZI)=(AR+I*AI)/(BR+I*BI)
!-----DO NOT USE IF BR=BI=0.
      REAL*8 AR,AI,BR,BI,YR,YI,ZR,ZI,W
      YR=BR
      YI=BI
      IF(ABS(YR).GT.ABS(YI)) THEN
      W=YI/YR
      YR=W*YI+YR
      ZR=(AR+W*AI)/YR
      ZI=(AI-W*AR)/YR
      ELSE
      W=YR/YI
      YI=W*YR+YI
      ZR=(W*AR+AI)/YI
      ZI=(W*AI-AR)/YI
      ENDIF
      RETURN
      END

! Exchange elements of a complex matrix
      SUBROUTINE EXCHG1(M,NM,N,A,B,D,J,K,L)
      REAL*8 A(NM,N),B(NM,N),D(N),F
      D(M)=J
      IF(J.NE.M) THEN
      DO I=1,K
         F=A(I,J)
         A(I,J)=A(I,M)
         A(I,M)=F
         F=B(I,J)
         B(I,J)=B(I,M)
         B(I,M)=F
      ENDDO
      DO I=L,N
         F=A(J,I)
         A(J,I)=A(M,I)
         A(M,I)=F
         F=B(J,I)
         B(J,I)=B(M,I)
         B(M,I)=F
      ENDDO
      ENDIF
      RETURN
      END

!end of file tceigen.f90