/************************************************************************
*   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.000000    2.000000    3.000000   -7.000000   12.000000           *
*    2.000000    4.000000    7.000000    3.000000   -1.000000           *
*    3.000000    7.000000   10.000000    8.000000    4.000000           *
*   -7.000000    3.000000    8.000000   -0.750000   -9.000000           *
*   12.000000   -1.000000    4.000000   -9.000000   10.000000           *
*                                                                       *
*  Error code = 0                                                       *
*                                                                       *
*  Eigenvalues:                                                         *
*    0.463350    0.000000                                               *
*   -7.774580    0.000000                                               *
*  -10.486545    0.000000                                               *
*   18.291821    0.000000                                               *
*   23.755955    0.000000                                               *
*                                                                       *
*  Eigenvectors (real part in lines):                                   *
*    0.208812    1.000000   -0.478027   -0.275000   -0.216915           *
*    1.000000   -0.326100    0.209620   -0.147689   -0.815422           *
*    0.467172   -0.007947   -0.507827    1.000000    0.264432           *
*    0.093053    0.594911    1.000000    0.455666    0.050740           *
*    0.705820   -0.012677    0.131486   -0.527499    1.000000           *
*                                                                       *
* 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 (real 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.000000    5.000000    3.000000   -7.000000   12.000000           *
*    2.000000    4.000000    7.000000    3.000000   -1.000000           *
*    3.000000    7.000000   10.000000    8.000000    4.000000           *
*   -7.000000    3.000000    8.000000   -0.750000   -9.000000           *
*   12.000000   -1.000000    4.000000   -9.000000   10.000000           *
*  Input Matrix (imag. part):                                           *
*    1.000000    1.000000    1.000000    1.000000    1.000000           *
*    1.000000    1.000000    1.000000    1.000000    1.000000           *
*    1.000000    1.000000    1.000000    1.000000    1.000000           *
*    1.000000    1.000000    1.000000    1.000000    1.000000           *
*    1.000000    1.000000    1.000000    1.000000    1.000000           *
*                                                                       *
*  Error code = 0                                                       *
*                                                                       *
*  Eigenvalues:                                                         *
*    0.888713    0.034707                                               *
*   -8.284268    0.024981                                               *
*  -10.363673    0.879708                                               *
*   18.894792    3.155582                                               *
*   23.114437    0.905021                                               *
*                                                                       *
*  Eigenvectors (real part in lines):                                   *
*    0.411502    1.000000   -0.532508   -0.217945   -0.417001           *
*    1.000000   -0.324189    0.269140   -0.277587   -0.872694           *
*    0.451233   -0.007289   -0.515266    1.000000    0.272732           *
*    0.179884    0.608599    1.000000    0.456406    0.069197           *
*    0.682657   -0.106663   -0.004456   -0.615170    1.000000           *
*                                                                       *
*  Eigenvectors (imag. part in lines):                                  *
*    0.000770    0.000000   -0.034750    0.009062   -0.005181           *
*    0.000000    0.001676    0.051622   -0.074119   -0.037648           *
*    0.043492   -0.070147   -0.050515    0.000000   -0.066367           *
*   -0.198515    0.075502    0.000000    0.271030   -0.314928           *
*    0.044071    0.193806    0.292159    0.172780    0.000000           *
*                                                                       *
* --------------------------------------------------------------------- *
* Reference:   From Numath Library By Tuan Dang Trong in Fortran 77     *
*              [BIBLI 18].                                              *
*                                                                       *
*                               C++ Release 1.1 By J-P Moreau, Paris    *
*                                         (www.jpmoreau.fr)             *
* --------------------------------------------------------------------- *
* Release 1.1: Added function MTPRINT and Example #2 (12/04/2004).      *
*************************************************************************
To link with basis_r.cpp and vmblock.cpp.
----------------------------------------- */
#include <basis.h>
#include <vmblock.h>

      //global variables
      int IERR;        //error code (must be zero)
	  int NM;      //1st dimension of input complex matrix
      int N;           //2nd dimension of input complex matrix
      REAL **AR, **AI; //pointers to real, imag. parts of complex matrix A
      REAL **ZR, **ZI; //complex matrix to store complex eigenvectors
      REAL *WR, *WI;   //real, imag. parts of complex eigenvalues
      REAL *WORK;      //working space, size = N
      boolean AVEC;    //TRUE = eigenvectors asked for 

      int I,J;         //loop index
	  void *vmblock;   //pointer to a working space for dynamic allocations
      char c; 
	  
      // print a matrix with caption
      void MPRINT(int N, int MMAX, int M, char *Name, REAL **A) {
        int I,J;
        printf("\n%s\n",Name);
        for (I=1; I<=N; I++) {
          for (J=1; J<=N; J++) printf("%12.6f", A[I][J]);
          printf("\n");
        }
      }

      // print transpose matrix with caption
      void MTPRINT(int N, int MMAX, int M, char *Name, REAL **A) {
        int I,J;
        printf("\n%s\n",Name);
        for (I=1; I<=N; I++) {
          for (J=1; J<=N; J++) printf("%12.6f", A[J][I]);
          printf("\n");
        }
      }

      //for debug only
      void VPRINT(int N, char *Name, REAL *V) {
        int I;
        printf("\n%s\n",Name);
        for (I=1; I<=N; I++) {
          printf("%12.6f", V[I]);
          if ((I % 5) == 0) printf("\n");
        }
      }

      //Headers of functions defined below
      void CBALAN(int NM,int N,REAL **AR,REAL **AI,int *LOW,int *UPP,REAL *D);
      void CUNITH(int NM,int N,int LOW,int UPP,REAL **AR,REAL **AI,REAL *DR, REAL *DI);
      void COMQR2(int NM,int N,int LOW,int UPP,REAL *GR,REAL *GI, REAL **HR, REAL **HI,
                  REAL *WR, REAL *WI, boolean AVEC, REAL **ZR, REAL **ZI, int *IERR);
      void CBALBK(int NM,int N,int LOW,int UPP,int M,REAL *D,REAL **Z,REAL **T);
      void NORMAL(int NM,int N,REAL *WR,REAL *WI,boolean AVEC,REAL **ZR,REAL **ZI,REAL *TR,REAL *TI);
      REAL ABSCD (REAL AR, REAL AI);
      void DCSQRT(REAL X,REAL Y,REAL *A,REAL *B);
      void DIVCD (REAL AR,REAL AI,REAL BR,REAL BI,REAL *ZR,REAL *ZI);
      void EXCHG1(int M,int NM,int N,REAL **A,REAL **B, REAL *D,int J,int K,int L);
      REAL DCABS (REAL XX,REAL YY);

      //Utility function
      int IMIN(int a,int b) {
        if (a<b) return a;
        else return b;
      }

    // main function
    void CEIGEN(REAL **AR,REAL **AI, int NM,int N, boolean AVEC,
	 	        REAL *WR,REAL *WI, REAL **ZR,REAL **ZI, REAL *WORK, int *IERR) {
/*--------------------------------------------------------------------
!     This Procedure 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           Boolean
!           =   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(N)                                             R*8
!
!     REFERENCE:
!     J.H.WILKINSON & C.REINSCH
!     LINEAR ALGEBRA   SPRINGER-VERLAG 1971
! --------------------------------------------------------------------*/
      int LOW,N1,N2,UPP;
      REAL *TMP1,*TMP2;       //temporary vectors, size = N
	  void *vmblock1 = NULL;

      if (N <= NM) {

		vmblock1 = vminit();
        TMP1 = (REAL *) vmalloc(vmblock1, VEKTOR, N, 0);
		TMP2 = (REAL *) vmalloc(vmblock1, VEKTOR, N, 0);
 
        N1=N+1;
        N2=N+N1;

        CBALAN(NM,N,AR,AI,&LOW,&UPP,WORK);
        CUNITH(NM,N,LOW,UPP,AR,AI,TMP1,TMP2);
        COMQR2(NM,N,LOW,UPP,TMP1,TMP2,AR,AI,WR,WI,AVEC,ZR,ZI,IERR);

        if (*IERR == 0) {
          CBALBK(NM,N,LOW,UPP,N,WORK,ZR,ZI);
          NORMAL(NM,N,WR,WI,AVEC,ZR,ZI,TMP1,TMP2);
        }
      }
      else
        *IERR= 9999;
      vmfree(vmblock1);
    } //CEIGEN

    void CUNITH(int NM,int N,int LOW,int UPP,REAL **AR,REAL **AI,REAL *DR, REAL *DI) {
//  UNITARY TRANSFORMATIONS TO REDUCE A COMPLEX SQUARE MATRIX
//  TO A HESSENBERG FORM.
      REAL XR,XI,G,F,H,S;
      int J,M;
      if (UPP-1 < LOW+1) return;
      for (M=LOW+1; M<UPP; M++) {
        DR[M]=0.0;
        DI[M]=0.0;
        S=0.0;
        H=0.0;
        for (I=M; I<=UPP; I++)
          S += ABS(AR[I][M-1])+ABS(AI[I][M-1]);
        if (S != 0.0)  {
 		  for (I=UPP; I>=M; I--) {
            DR[I]=AR[I][M-1]/S;
            DI[I]=AI[I][M-1]/S;
            H += DR[I]*DR[I]+DI[I]*DI[I];
          }
          G=SQRT(H);
          F=ABSCD(DR[M],DI[M]);
          if (F == 0.0) {
            DR[M]=G;
            AR[M][M-1]=S;
          }
          else {
            H=H+F*G;
            G=G/F;
            DR[M]=(1.0+G)*DR[M];
            DI[M]=(1.0+G)*DI[M];
          }
          for (J=M; J<=N; J++) {
            XR=0.0;
            XI=0.0;
            for (I=UPP; I>=M; I--) {
              XR += (DR[I]*AR[I][J]+DI[I]*AI[I][J]);
              XI += (DR[I]*AI[I][J]-DI[I]*AR[I][J]);
            }
            XR /= H;
            XI /= H;
            for (I=M; I<=UPP; I++) {
              AR[I][J] += (-XR*DR[I]+XI*DI[I]);
              AI[I][J] -= (XR*DI[I]+XI*DR[I]);
            }
          }
          for (I=1; I<=UPP; I++) {
            XR=0.0;
            XI=0.0;
            for (J=UPP; J>=M; J--) {
              XR += (DR[J]*AR[I][J]-DI[J]*AI[I][J]);
              XI += (DR[J]*AI[I][J]+DI[J]*AR[I][J]);
            }
            XR=XR/H;
            XI=XI/H;
            for (J=M; J<=UPP; J++) {
              AR[I][J] -= (XR*DR[J]+XI*DI[J]);
              AI[I][J] += (XR*DI[J]-XI*DR[J]);
            }
          }
          DR[M] *= S;
          DI[M] *= S;
          AR[M][M-1]=-G*AR[M][M-1];
          AI[M][M-1]=-G*AI[M][M-1];
		} // if S<>0
      } //M loop
    } //CUNITH()

    void COMQR2(int NM,int N,int LOW,int UPP,REAL *GR,REAL *GI, REAL **HR, REAL **HI,
		        REAL *WR, REAL *WI, boolean AVEC, REAL **ZR, REAL **ZI, int *IERR) {
/*---------------------------------------------------------------------
      QR ALGORITH APPLIED TO A COMPLEX HESSENBERG MATRIX
      (translated from Fortran 77).
  --------------------------------------------------------------------- */
//    Labels: e5,e7,e10,e15,e20,e35,e45,e50,e55
      REAL SR,SI,TR,TI,YR,YI,XR,XI,VR,VI,MACHEP,NORM,EPS,EPS1;
      int I,J,K,L,M,IN0,NA,ITS;
//    SEEK FLOATING POINT RELATIVE PRECISION FOR REAL
      EPS=1.0;
      do {
        EPS *= 0.5;
        EPS1=1.0+EPS;
      } while (EPS1 > 1.0);
      MACHEP=2.0*EPS;

      *IERR=0;
      if (AVEC)
		for (I=1; I<=N; I++) {
		  for (J=1; J<=N; J++) {
            ZR[I][J]=0.0;
            ZI[I][J]=0.0;
          }
          ZR[I][I]=1.0;
        }

//    ACCUMULATE UNITARY TRANSFORMATIONS
	  for (I=UPP-1; I>LOW; I--) {
         if (((HR[I][I-1] == 0.0) && (HI[I][I-1] == 0.0)) ||
            ((GR[I] == 0.0) && (GI[I] == 0.0)))  goto e5;
         NORM=HR[I][I-1]*GR[I]+HI[I][I-1]*GI[I];
         for (K=I+1; K<=UPP; K++) {
           GR[K]=HR[K][I-1];
           GI[K]=HI[K][I-1];
         }
         for (J=I; J<=UPP; J++) {
           SR=0.0;
           SI=0.0;
           for (K=I; K<=UPP; K++) {
             SR += GR[K]*ZR[K][J]+GI[K]*ZI[K][J];
             SI += GR[K]*ZI[K][J]-GI[K]*ZR[K][J];
           }
           SR /= NORM;
           SI /= NORM;
           for (K=I; K<=UPP; K++) {
             ZR[K][J] += SR*GR[K]-SI*GI[K];
             ZI[K][J] += SR*GI[K]+SI*GR[K];
           }
         }
e5:  ;}

//    CREATE REAL SUBDIAGONAL
      L=LOW+1;
      for (I=L; I<=UPP; I++) {
         M=IMIN(I+1,UPP);
         if (HI[I][I-1] == 0.0) goto e7;
         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.0;
         for (J=I; J<=N; J++) {
           SI= YR*HI[I][J]-YI*HR[I][J];
           HR[I][J]=YR*HR[I][J]+YI*HI[I][J];
           HI[I][J]=SI;
         }
         for (J=1; J<=M; J++) {
           SI= YR*HI[J][I]+YI*HR[J][I];
           HR[J][I]=YR*HR[J][I]-YI*HI[J][I];
           HI[J][I]=SI;
         }
         if (AVEC)
		   for (J=LOW; J<=UPP; J++) {
             SI= YR*ZI[J][I]+YI*ZR[J][I];
             ZR[J][I]=YR*ZR[J][I]-YI*ZI[J][I];
             ZI[J][I]=SI;
           }
 e7: ;}

//    STORE ROOTS ALREADY ISOLATED BY FUNCTION BALANC
      for (I=1; I<=N; I++)
	    if (I < LOW || I > UPP) {
          WR[I]=HR[I][I];
          WI[I]=HI[I][I];
        }

      IN0=UPP;
      TR=0.0;
      TI=0.0;

//    SEEK EIGENVALUES
e10:  if (IN0 < LOW) goto e50;
      ITS=0;
      NA=IN0-1;

//    WE LOOK AT REAL PART OF SUBDIAGONAL ELEMENT
e15:  for (L=IN0; L>=LOW; L--) {
         //here index 0 is temporarily used!
         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]) <= SR || L == LOW) goto e20;
      }

//    SHIFTING
e20:  if (L == IN0) goto e45;
      if (ITS == 30) goto e55;
      if (ITS == 10 || ITS == 20) {
//    EXTRA SHIFTING
        SR=ABS(HR[IN0][NA])+ABS(HR[NA][IN0-2]);
        SI=0.0;
      }
      else {
         SR=HR[IN0][IN0];
         SI=HI[IN0][IN0];
         XR=HR[NA][IN0]*HR[IN0][NA];
         XI=HI[NA][IN0]*HR[IN0][NA];

         if (XR != 0.0 || XI != 0.0) {
           YR=(HR[NA][NA]-SR)*0.5;
           YI=(HI[NA][NA]-SI)*0.5;
           DCSQRT(YR*YR-YI*YI+XR,2.0*YR*YI+XI,&VR,&VI);
           if (YR*VR+YI*VI < 0.0) {
             VR=-VR;
             VI=-VI;
           }
           DIVCD(XR,XI,YR+VR,YI+VI,&VR,&VI);
           SR=SR-VR;
           SI=SI-VI;
         }
      }

      for (I=LOW; I<=IN0; I++) {
        HR[I][I] -= SR;
        HI[I][I] -= SI;
      }

      TR += SR;
      TI += SI;
      ITS++;

//    TRIANGULARY DECOMPOSITION H = Q * R
      for (I=L+1; I<=IN0; I++) {
         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.0;
         HR[I][I-1]=0.0;
         HI[I][I-1]=SR/NORM;
         for (J=I; J<=N; J++) {
           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;
         }
      }

      SI=HI[IN0][IN0];

      if (SI == 0.0) goto e35;
      NORM=ABSCD(HR[IN0][IN0],SI);
      SR=HR[IN0][IN0]/NORM;
      SI /= NORM;
      HI[IN0][IN0]=0.0;
      HR[IN0][IN0]=NORM;
      if (IN0 != N)
		for (J=IN0+1; J<=N; J++) {
          YR=HR[IN0][J];
          YI=HI[IN0][J];
          HR[IN0][J]=SR*YR+SI*YI;
          HI[IN0][J]=SR*YI-SI*YR;
        }

//    INVERSE OPERATION
e35:  for (J=L+1; J<=IN0; J++) {
        XR=WR[J-1];
        XI=WI[J-1];
        for (I=1; I<=J; I++) {
          YR=HR[I][J-1];
          YI=0.0;
          VR=HR[I][J];
          VI=HI[I][J];
          if (I != J) {
            YI=HI[I][J-1];
            HI[I][J-1]=XR*YI+XI*YR+HI[J][J-1]*VI;
          }
          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;
        }
        if (AVEC)
		  for (I=LOW; I<=UPP; I++) {
            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;
          }      
	  }

      if (SI == 0.0) goto e15;
      for (I=1; I<=IN0; I++) {
         YR=HR[I][IN0];
         YI=HI[I][IN0];
         HR[I][IN0]=SR*YR-SI*YI;
         HI[I][IN0]=SR*YI+SI*YR;
      }
      if (AVEC)
	    for (I=LOW; I<=UPP; I++) {
          YR=ZR[I][IN0];
          YI=ZI[I][IN0];
          ZR[I][IN0]=SR*YR-SI*YI;
          ZI[I][IN0]=SR*YI+SI*YR;
        }
      goto e15;

//    ONE ROOT HAS BEEN FOUND
e45:  HR[IN0][IN0]=HR[IN0][IN0]+TR;
      HI[IN0][IN0]=HI[IN0][IN0]+TI;
      WR[IN0]=HR[IN0][IN0];
      WI[IN0]=HI[IN0][IN0];
      IN0=NA;
      goto e10;

//    ALL ROOTS HAVE BEEN FOUND
//    START CALCULATION OF EIGENVECTORS BY SUBSTITUTION
e50:  if (AVEC) {
      NORM=0.0;
      for (I=1; I<=N; I++)
        for (J=I; J<=N; J++)
          NORM += ABS(HR[I][J])+ABS(HI[I][J]);
		if (NORM == 0.0)  {
          printf(" 1 IERR=%d\n", *IERR);
		  return;
        }
      for (IN0=N; IN0>1; IN0--) {
         XR=WR[IN0];
         XI=WI[IN0];
         HR[IN0][IN0]=1.0;
         HI[IN0][IN0]=0.0;
         for (I=IN0-1; I>0; I--) {
            VR=HR[I][IN0];
            VI=HI[I][IN0];
            if (I != IN0-1)
			  for (J=I+1; J<IN0; J++) {
                VR=VR+HR[I][J]*HR[J][IN0]-HI[I][J]*HI[J][IN0];
                VI=VI+HR[I][J]*HI[J][IN0]+HI[I][J]*HR[J][IN0];
              }
            YR=XR-WR[I];
            YI=XI-WI[I];
            if (YR == 0.0 && YI == 0.0)  YR=NORM*MACHEP;
            DIVCD(VR,VI,YR,YI,&TR,&TI);
            HR[I][IN0]=TR;
            HI[I][IN0]=TI;
         }
      }
      for (I=1; I<N; I++) {
         if (I < LOW || I > UPP)
		   for (J=I+1; J<=N; J++) {
             ZR[I][J]=HR[I][J];
             ZI[I][J]=HI[I][J];
           }
      }

//    END CALCULATION OH HESSENBERG EIGENVECTORS
	   for (J=N; J>LOW; J--) {
         M=IMIN(J-1,UPP);
         for (I=LOW; I<=UPP; I++) {
           VR=ZR[I][J];
           VI=ZI[I][J];
           for (K=LOW; K<=M; K++) {
             VR += ZR[I][K]*HR[K][J]-ZI[I][K]*HI[K][J];
             VI += ZR[I][K]*HI[K][J]+ZI[I][K]*HR[K][J];
           }
           ZR[I][J]=VR;
           ZI[I][J]=VI;
         }
       }
      } //if AVEC

//    END CALCULATION OF EIGENVECTORS OF INPUT MATRIX
      return;

//    NO CONVERGENCE AFTER 30 ITERATIONS!
e55:  *IERR = IN0;
    } //COMQR2()


    void CBALAN(int NM,int N,REAL **AR,REAL **AI,int *LOW,int *UPP,REAL *D) {
//  BALANCE A COMPLEX MATRIX (VERSION DERIVED FROM BALANC)
//  No effect for a symmetrical matrix).
//  Labels: e10,e11,e20,e22,e30,e40,e50
      REAL B2,C,F,G,R,S;
      boolean NOCONV;
      int B,J,J1,K,K1,L,L1;
      B=16;
      B2=B*B;
      L=1;
      K=N;
e10:  for (J=K; J>0; J--) {
         J1=J;
         K1=K;
         L1=L;
         for (I=1; I<=K; I++) {
           if (I != J)
             if (AR[J][I] != 0.0 || AI[J][I] != 0.0) goto e11;
         }
         EXCHG1(K,NM,N,AR,AI,D,J1,K1,L1);
         K--;
         goto e10;
e11: ;}
e20:  for (J=L; J<=K; J++) {
         J1=J;
         K1=K;
         L1=L;
         for (I=L; I<=K; I++) {
            if (I != J)
              if (AR[I][J] != 0.0 || AI[I][J] != 0.0) goto e22;
         }
         EXCHG1(L,NM,N,AR,AI,D,J1,K1,L1);
         L++;
         goto e20;
e22: ;}
      *LOW=L;
      *UPP=K;
      for (I=L; I<=K; I++)  D[I]=1.0;
e30:  NOCONV=FALSE;
      for (I=L; I<=K; I++) {
        C=0.0;
        R=0.0;
        for (J=L; J<=K; J++)
		  if (J != I) {
            C += ABS(AR[J][I])+ABS(AI[J][I]);
            R += ABS(AR[I][J])+ABS(AI[I][J]);
          }
        G=R/B;
        F=1.0;
        S=C+R;
e40:    if (C < G) {
          F=F*B;
          C=C*B2;
          goto e40;
        }
        G=R*B;
e50:    if (C >= G) {
          F=F/B;
          C=C/B2;
          goto e50;
        }
        if ((C+R)/F < 0.95*S) {
          G=1.0/F;
          D[I] *= F;
          NOCONV=TRUE;
          for (J=L; J<=N; J++) {
            AR[I][J] *= G;
            AI[I][J] *= G;
          }
          for (J=1; J<=K; J++) {
            AI[J][I] *= F;
            AR[J][I] *= F;
          }
        }
      }
      if (NOCONV) goto e30;
	} //CBALAN()

    void CBALBK(int NM,int N,int LOW,int UPP,int M,REAL *D,REAL **Z,REAL **T) {
//  This Procedure finds the eigenvectors of a complex matrix,
//  as output of function CBALAN.
//  Label: e10
      REAL S;
      int  I,J,K,L;
      for (I=LOW; I<=UPP; I++) {
        S=D[I];
        for (J=1; J<=M; J++) {
          Z[I][J] *= S;
          T[I][J] *= S;
        }
      }
      for (L=1; L<=N; L++) {
        I=L;
        if (I >= LOW && I <= UPP) goto e10;
        if (I < LOW)  I=LOW-L;
        K=int(D[I]);
        if (K != I)
		  for (J=1; J<=M; J++) {
            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;
          }
e10: ;}
    }

    void NORMAL(int NM,int N,REAL *WR,REAL *WI,boolean AVEC,REAL **ZR,REAL **ZI,
		REAL *TR,REAL *TI) {
//  Sort eigenvalues by absolute values in ascending order and
//  normalize.
//  Label: e5
      REAL ZM,Z1,Z2,VR,VI;
      int IM,J,K;
//    SORT SOLUTIONS IN ASCENDING ORDER
      for (J=2; J<=N; J++) {
         VR=WR[J];
         VI=WI[J];
         if (AVEC)
		   for (K=1; K<=N; K++) {
             TR[K]=ZR[K][J];
             TI[K]=ZI[K][J];
           }

		 for (I=J-1; I>0; I--) {
           if (ABSCD(WR[I],WI[I]) <= ABSCD(VR,VI)) goto e5;
           WR[I+1]=WR[I];
           WI[I+1]=WI[I];
           if (AVEC)
			 for (K=1; K<=N; K++) {
               ZR[K][I+1]=ZR[K][I];
               ZI[K][I+1]=ZI[K][I];
             }
         }
         I=0;
e5:      WR[I+1]=VR;
         WI[I+1]=VI;
         if (AVEC)
		   for (K=1; K<=N; K++) {
             ZR[K][I+1]=TR[K];
             ZI[K][I+1]=TI[K];
           }
      }

//    NORMALIZE WITH RESPECT TO BIGGEST ELEMENT
      if (AVEC)
	    for (J=N; J>0; J--) {
          ZM= 0.0;
          for (I=1; I<=N; I++) {
            VR=ZR[I][J];
            VI=ZI[I][J];
            Z1=ABSCD(VR,VI);
            if (Z1 >= ZM) {
              IM = I;
              ZM = Z1;
            }
          }
          Z1 = ZR[IM][J];
          Z2 = ZI[IM][J];
          for (I = 1; I<=N; I++) {
            DIVCD(ZR[I][J],ZI[I][J],Z1,Z2,&VR,&VI);
            ZR[I][J] = VR;
            ZI[I][J] = VI;
          }
        }
    } //NORMAL()

//  Utility routines for complex numbers
    REAL ABSCD(REAL AR,REAL AI) {
//  ABSOLUTE VALUE OF A COMPLEX NUMBER C=AR+I*AI
//    ABSCD,AR,AI REAL*8
//    ABSCD=SQRT(AR**2+AI**2)
      REAL XR,XI,W;
      XR=ABS(AR);
      XI=ABS(AI);
      if (XR <= XI) {
        W=XR;
        XR=XI;
        XI=W;
      }
      if (XI == 0.0)
        return XR;
      else
        return (XR*SQRT(1.0+(XI/XR)*(XI/XR)));
    }

    void DCSQRT(REAL X,REAL Y, REAL *A, REAL *B) {
//  SQUARE ROOT OF A COMPLEX NUMBER  A+I*B = RAC(X+I*Y)
	  if (X == 0.0 && Y == 0.0) {
        *A=0.0;
        *B=0.0;
      }
      else {
        *A=SQRT((ABS(X)+DCABS(X,Y))*0.5);
        if (X >= 0.0)
          *B=Y/(*A+*A);
        else
          if (Y < 0.0)
            *B=-*A;
          else {
            *B=*A;
            *A=Y/(*B+*B);
          }
      }
    }

//  ABSOLUTE VALUE OF A COMPLEX NUMBER X+I*Y
    REAL DCABS(REAL XX,REAL YY) {
      REAL X,Y,W;
      X=ABS(XX);
      Y=ABS(YY);
      if (X == 0.0)
        W=Y;
      else
        if (Y == 0.0)
          W=X;
        else
          if (X > Y)
            W=X*SQRT(1.0+(Y/X)*(Y/X));
          else
            W=Y*SQRT(1.0+(X/Y)*(Y/X));
      return W;
    }

    void DIVCD (REAL AR,REAL AI,REAL BR,REAL BI, REAL *ZR,REAL *ZI) {
//  EFFECT DIVISION Z=(ZR+I*ZI)=(AR+I*AI)/(BR+I*BI)
//  ---- DO NOT USE IF BR=BI=0 ----
      REAL YR,YI,W;
      YR=BR;
      YI=BI;
      if (ABS(YR) > ABS(YI)) {
        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;
      }
    }

// Exchange elements of a complex matrix
    void EXCHG1(int M,int NM,int N, REAL **A,REAL **B, REAL *D,
		        int J,int K,int L)  {
      REAL F;
      D[M]=J;
      if (J != M) {
	    for (I=1; I<=K; I++) {
          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;
        }
        for (I=L; I<=N; I++) {
          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;
        }
      }
    }


void main()  {

  NM = 5;
  N  = 5;   // size of input complex matrix

// dynamic allocations
  vmblock = vminit();
  AR = (REAL **) vmalloc(vmblock, MATRIX,  NM+1, N+1);
  AI = (REAL **) vmalloc(vmblock, MATRIX,  NM+1, N+1);
  ZR = (REAL **) vmalloc(vmblock, MATRIX,  NM+1, N+1);
  ZI = (REAL **) vmalloc(vmblock, MATRIX,  NM+1, N+1);
  WR = (REAL *)  vmalloc(vmblock, VEKTOR,  N+1,  0);
  WI = (REAL *)  vmalloc(vmblock, VEKTOR,  N+1,  0);
  WORK = (REAL *) vmalloc(vmblock, VEKTOR, N+1,  0);

  for (I=0; I<=NM; I++)
    for (J=0; J<=NM; J++) {
      AR[I][J] = 0.0;
      AI[I][J] = 0.0;
    }

  for (I=0; I<=N; I++)  WORK[I] = 0.0;

  //data section #1
  AR[1][1] = 1.0; AR[1][2] = 2.0; AR[1][3] =  3.0; AR[1][4] = -7.0; AR[1][5] = 12.0;
  AR[2][1] = 2.0; AR[2][2] = 4.0; AR[2][3] =  7.0; AR[2][4] =  3.0; AR[2][5] = -1.0;
  AR[3][1] = 3.0; AR[3][2] = 7.0; AR[3][3] = 10.0; AR[3][4] =  8.0; AR[3][5] =  4.0;
  AR[4][1] =-7.0; AR[4][2] = 3.0; AR[4][3] =  8.0; AR[4][4] =-0.75; AR[4][5] = -9.0;
  AR[5][1] =12.0; AR[5][2] =-1.0; AR[5][3] =  4.0; AR[5][4] = -9.0; AR[5][5] = 10.0;

  /* data section #2
  AR[1][1] = 1.0; AR[1][2] = 5.0; AR[1][3] =  3.0; AR[1][4] = -7.0; AR[1][5] = 12.0;
  AR[2][1] = 2.0; AR[2][2] = 4.0; AR[2][3] =  7.0; AR[2][4] =  3.0; AR[2][5] = -1.0;
  AR[3][1] = 3.0; AR[3][2] = 7.0; AR[3][3] = 10.0; AR[3][4] =  8.0; AR[3][5] =  4.0;
  AR[4][1] =-7.0; AR[4][2] = 3.0; AR[4][3] =  8.0; AR[4][4] =-0.75; AR[4][5] = -9.0;
  AR[5][1] =12.0; AR[5][2] =-1.0; AR[5][3] =  4.0; AR[5][4] = -9.0; AR[5][5] = 10.0;

  for (I=1; I<=N; I++)
    for (J=1; J<=N; J++)
      AI[I][J] = 1.0;      */

  MPRINT(N,NM,N,"  Input Matrix (real part):", AR);
  /* #2 only:
  MPRINT(N,NM,N,"  Input Matrix (imag. part):", AI);
  printf(" Any key + Enter to continue... "); scanf("%s", c);  */

  //call main subroutine
  CEIGEN(AR,AI,NM,N,TRUE,WR,WI,ZR,ZI,WORK,&IERR);

  printf("\n  Error code = %d\n\n", IERR);

  if (IERR == 0) {  //no error
    printf("  Eigenvalues:\n");
    for (I=1; I<=N; I++)
      printf("%12.6f%12.6f\n", WR[I], WI[I]);

      MTPRINT(N,NM,N,"  Eigenvectors (real part in lines):", ZR);
	  /* #2 only:
	  MTPRINT(N,NM,N,"  Eigenvectors (imag. part in lines):", ZI);  */
      printf("\n");
  }

  //free memory
  vmfree(vmblock);

} //end of main program

// end of file tceigen.cpp