/***********************************************************************
*   Find Eigenvalues and Eigenvectors of a real symmetric matrix by    *
*   using Householder reduction and  QL method.                        *
* -------------------------------------------------------------------- *
* SAMPLE RUNS:                                                         *
*                                                                      *
* Example #1:    (symmetric matrix)                                    *
*                                                                      *
*          4 -2 -1  0                                                  *
*     A = -2  4  0 -1                                                  *
*         -1  0  4 -2                                                  *
*          0 -1 -2  4                                                  *
*                                                                      *
*  Output file ttred2.lst contains:                                    *
*                                                                      *
* EIGENVALUES AND EIGENVECTORS OF A SYMMETRIC REAL MATRIX              *
* BY HOUSEHOLDER REDUCTION AND QL METHOD.                              *
*                                                                      *
* N = 4                                                                *
*                                                                      *
* Input symmetric matrix:                                              *
*  4.000000 -2.000000 -1.000000  0.000000                              *
* -2.000000  4.000000  0.000000 -1.000000                              *
* -1.000000  0.000000  4.000000 -2.000000                              *
*  0.000000 -1.000000 -2.000000  4.000000                              *
*                                                                      *
* Error code: 0                                                        *
*                                                                      *
* Eigenvalues:                                                         *
*  1     1.00000000                                                    *
*  2     3.00000000                                                    *
*  3     5.00000000                                                    *
*  4     7.00000000                                                    *
*                                                                      *
* Eigenvectors (in lines):                                             *
*  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                              *
*                                                                      *
*                                                                      *
* Example #2:                                                          *
*                                                                      *
*          1  2  3 -7    12                                            *
*          2  4  7  3    -1                                            *
*     A =  3  7 10  8     4                                            *
*         -7  3  8 -0.75 -9                                            *
*         12 -1  4 -9    10                                            *
*                                                                      *
*                                                                      *
* EIGENVALUES AND EIGENVECTORS OF A SYMMETRIC MATRIX                   *
* BY HOUSEHOLDER REDUCTION AND QL METHOD.                              *
*                                                                      *
* N = 5                                                                *
*                                                                      *
* Input symmetric matrix:                                              *
*  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:                                                         *
*  1   -10.48654517                                                    *
*  2    -7.77457973                                                    *
*  3     0.46334953                                                    *
*  4    18.29182060                                                    *
*  5    23.75595477                                                    *
*                                                                      *
* Eigenvectors (in lines):                                             *
*  0.467172 -0.007947 -0.507827  1.000000  0.264432                    *
*  1.000000 -0.326100  0.209620 -0.147689 -0.815422                    *
*  0.208812  1.000000 -0.478027 -0.275000 -0.216915                    *
*  0.093053  0.594911  1.000000  0.455666  0.050740                    *
*  0.705820 -0.012677  0.131486 -0.527499  1.000000                    *
*                                                                      *
* -------------------------------------------------------------------- *
* Reference:   From Numath Library By Tuan Dang Trong in Fortran 77    *
*              [BIBLI 18].                                             *
*                                                                      *
*                               C++ Release 1.1 By J-P Moreau, Paris   *
* -------------------------------------------------------------------- *
* Release 1.1:  added results of example #2.                           *
***********************************************************************/
#include <basis.h>
#include <vmblock.h> 

      REAL *D, *E;
      REAL **Z;
      REAL zmax;
      int  I,IERR,J,N;

      FILE *fp;
	  void *vmblock = NULL;

      REAL Sign(REAL a, REAL b) {
        if (b <ZERO) return (- ABS(a));
        else return (ABS(a));
      }

	void TQL2(int N, REAL *D, REAL *E, REAL **Z, int *IER)  {
/*--------------------------------------------------------------------------
!     QL METHOD TO DETERMINE THE EIGENVALUES AND EIGENVECTORS OF:
!
!       1)  A SYMMETRIC TRIDIAGONAL MATRIX.
!       2)  A FULL SYMMETRIC MATRIX AFTER A PREVIOUS CALL TO TRED2.
!
!     CALLING MODE:
!               TQL2(N,D,E,Z,&IER)
!     INPUTS:
!     NM  (I4)  1ST DIMENSION OF MATRICES A AND Z IN CALLING PROGRAM
!     N   (I4)  SIZE OF Z
!     D  (R*8)  MAIN DIAGONAL (N) OF THE TRIDIAGONAL MATRIX
!     E  (R*8)  SUB-DIAGONAL (N) OF THE TRIDIAGONAL MATRIX
!     Z  (R*8)  TABLE (NM,N) STORING THE UNITY MATRIX IF THE TRIDIAGONAL
!               MATRIX IS DEFINED BY D AND E, CASE #1.
!               FOR CASE #2, IT CONTAINS THE ELEMENTS OF THE TRANSFORMATION
!               MATRIX AFTER A CALL TO TRED2.
!     OUTPUTS:
!     D  (R*8)  EIGENVALUES
!     Z  (R*8)  EIGENVECTORS
!     IER (I4)  ERROR CODE = 0,  CONVERGENCE OK.
!                          = L,  NO CONVERGENCE FOR THE Lth EIGENVALUE
!
!     REFERENCE:
!     J.H.WILKINSON,-C.REINSCH,R.S.MARTIN
!     HANDBOOK FOR AUTOMATIC COMPUTATION, VOL.2, LINEAR ALGEBRA
!     SPRINGER-VERLAG 1971.
!-------------------------------------------------------------------------*/
//  Labels: e10,e12,e18,e20,e26,e30,e34,e36
      REAL B,C,F,G,H,P,R,S,EPS,EPS1;
      int  I,J,K,L,M,JM;

      EPS=ZERO; JM=30;
      *IER = 0;
      if (N == 1)  return;

//    MACHINE EPSILON

      if (EPS != ZERO) goto e12;
      EPS = 1.0;
e10:  EPS = EPS/TWO;
      EPS1 = ONE + EPS;
      if (EPS1 > ONE) goto e10;

e12:  for (I=2; I<=N; I++)  E[I-1] = E[I];
      E[N] = ZERO;
      F = ZERO;
      B = ZERO;

      for (L=1; L<=N; L++) {
        J = 0;
        H = EPS*(ABS(D[L])+ABS(E[L]));
        if (B < H)  B = H;

//      SEEK SMALLEST ELEMENT OF SUBDIAGONAL

        for (M=L; M<=N; M++)
          if (ABS(E[M]) <= B) goto e18;

e18:    if (M == L)  goto e26;

//      START ITERATION

e20:    if (J == JM) goto e36;
        J = J+1;

//      SHIFT

        G = D[L];
        P = (D[L+1]-G)/(TWO*E[L]);
        R = SQRT(P*P+ONE);
        D[L] = E[L]/(P + Sign(R,P));
        H = G-D[L];
        for (I=L+1; I<=N; I++)  D[I] = D[I] - H;
        F = F + H;

//      QL TRANSFORMATION

        P = D[M];
        C = ONE;
        S = ZERO;
        for (I=M-1; I>=L; I--) {
          G = C*E[I];
          H = C*P;
          if (ABS(P) >= ABS(E[I])) {
            C = E[I]/P;
            R = SQRT(C*C+ONE);
            E[I+1] = S*P*R;
            S = C/R;
            C = ONE/R;
          }
          else {
            C = P/E[I];
            R = SQRT(C*C+ONE);
            E[I+1] = S*E[I]*R;
            S = ONE/R;
            C = C*S;
          }
          P = C*D[I]-S*G;
          D[I+1] = H+S*(C*G+S*D[I]);

//      ELEMENTS OF EIGENVECTORS

          for (K=1; K<=N; K++) {
            H = Z[K][I+1];
            Z[K][I+1] = S*Z[K][I]+C*H;
            Z[K][I] = Z[K][I]*C-S*H;
          }
        }
        E[L] = S*P;
        D[L] = C*P;
        if (ABS(E[L]) > B) goto e20;

//      CONVERGENCE

e26:    D[L] = D[L] + F;
      }

//      SORT EIGENVALUES AND EIGENVECTORS
//      IN ASVENDING ORDER

      for (L=2; L<=N; L++) {
        I = L-1;
        K = I;
        P = D[I];
        for (J=L; J<=N; J++) {
          if (D[J] >= P) goto e30;
          K = J;
          P = D[J];
e30:   ;}
        if (K == I)  goto e34;
        D[K] = D[I];
        D[I] = P;
        for (J=1; J<=N; J++) {
          P = Z[J][I];
          Z[J][I] = Z[J][K];
          Z[J][K] = P;
        }
e34: ;}
      return;

//    NO CONVERGENCE

e36:  *IER = L;
   } //TQL2

void TRED2(double **a, int n, double *d, double *e)  {
/*----------------------------------------------------------------- 
   HOUSEHOLDER REDUCTION OF MATRIX a TO A TRIDIAGONAL FORM.
   Algorithm: Martin et al., Num. Math. 11, 181-195, 1968.
   Ref: Smith et al., Matrix Eigensystem Routines -- EISPACK Guide
        Springer-Verlag, 1976, pp. 489-494.
        W H Press et al., Numerical Recipes in C, Cambridge U P,
        1988, pp. 373-374.
-------------------------------------------------------------------
  INPUTS:
    a  (R*8)  TABLE (N,N) STORING THE REAL ELEMENTS OF INPUT
              SYMMETRIC MATRIX
    n   (I4)  SIZE OF a
  OUTPUTS:
    d  (R*8)  MAIN DIAGONAL(N) OF THE TRIDIAGONAL MATRIX
    e  (R*8)  SUB-DIAGONAL (N) OF THE TRIDIAGONAL MATRIX
    a  (R*8)  TABLE (N,N) STORING THE ELEMENTS OF THE
              TRANSFORMATION MATRIX.
'----------------------------------------------------------------*/
int    l, k, j, i;
double scale, hh, h, g, f;

for (i = n; i >= 2; i--)  {
    l = i - 1;
    h = scale = 0.0;
    if (l > 1)  {
       for (k = 1; k <= l; k++)
           scale += fabs(a[i][k]);
       if (scale == 0.0)
          e[i] = a[i][l];
       else  {
          for (k = 1; k <= l; k++) {
            a[i][k] /= scale;
            h += a[i][k] * a[i][k];
          }
          f = a[i][l];
          g = (f>0 ? -sqrt(h) : sqrt(h));
          e[i] = scale * g;
          h -= f * g;
          a[i][l] = f - g;
          f = 0.0;
          for (j = 1; j <= l; j++) {
            a[j][i] = a[i][j]/h;
            g = 0.0;
            for (k = 1; k <= j; k++)
              g += a[j][k] * a[i][k];
            for (k = j+1; k <= l; k++)
              g += a[k][j] * a[i][k];
            e[j] = g / h;
            f += e[j] * a[i][j];
          }
          hh = f / (h + h);
          for (j = 1; j <= l; j++) {
            f = a[i][j];
            e[j] = g = e[j] - hh * f;
            for (k = 1; k <= j; k++)
              a[j][k] -= (f * e[k] + g * a[i][k]);
          }
       }
    }
    else
        e[i] = a[i][l];
    d[i] = h;
}
d[1] = 0.0;
e[1] = 0.0;

for (i = 1; i <= n; i++)
    {
    l = i - 1;
    if (d[i])
       {
       for (j = 1; j <= l; j++)
           {
           g = 0.0;
           for (k = 1; k <= l; k++)
               g += a[i][k] * a[k][j];
           for (k = 1; k <= l; k++)
               a[k][j] -= g * a[k][i];
           }
       }
       d[i] = a[i][i];
       a[i][i] = 1.0;
       for (j = 1; j <= l; j++)
           a[j][i] = a[i][j] = 0.0;
    }
} //tred2()


void main()  {

    N = 4;

    // allocate memory for matrix A and vectors B, X, D (index 0 not used)
    vmblock = vminit();
    Z   = (REAL **) vmalloc(vmblock, MATRIX, N+1, N+1);
    D   = (REAL *)  vmalloc(vmblock, VEKTOR, N+1, 0);
    E   = (REAL *)  vmalloc(vmblock, VEKTOR, N+1, 0);

    if (! vmcomplete(vmblock))
      LogError ("No Memory", 0, __FILE__, __LINE__);

    // define input symmetric matrix
    Z[1][1]= 4.0; Z[1][2]=-2.0; Z[1][3]=-1.0; Z[1][4]= 0.0;
    Z[2][1]=-2.0; Z[2][2]= 4.0; Z[2][3]= 0.0; Z[2][4]=-1.0;
    Z[3][1]=-1.0; Z[3][2]= 0.0; Z[3][3]= 4.0; Z[3][4]=-2.0;
    Z[4][1]= 0.0; Z[4][2]=-1.0; Z[4][3]=-2.0; Z[4][4]= 4.0;

	fp = fopen("ttred2.lst","w");

    fprintf(fp,"\n EIGENVALUES AND EIGENVECTORS OF A SYMMETRIC REAL MATRIX\n");
    fprintf(fp," BY HOUSEHOLDER REDUCTION AND QL METHOD.\n\n");
    fprintf(fp," N = %d\n\n", N);
    fprintf(fp," Input symmetric matrix:\n");
    for (I=1; I<=N; I++) {
	  for (J=1; J<=N; J++)
	    fprintf(fp,"%10.6f", Z[I][J]);
      fprintf(fp,"\n");
    }

    TRED2(Z,N,D,E);
    TQL2(N,D,E,Z,&IERR);

    fprintf(fp,"\n Error code: %d\n\n", IERR);
    // print eigenvalues
    fprintf(fp," Eigenvalues:\n");
    for (I=1; I<=N; I++)
      fprintf(fp,"%3d%15.8f\n", I, D[I]);
    // print normalized eigenvectors (with respect to unity)
    fprintf(fp,"\n Eigenvectors (in lines):\n");
    for (J=1; J<=N; J++) {
      zmax=Z[1][J];
      for (I=2; I<=N; I++)
	    if (ABS(Z[I][J]) > ABS(zmax))   zmax=Z[I][J];
      for (I=1; I<=N; I++)
		  fprintf(fp,"%10.6f", Z[I][J]/zmax);
      fprintf(fp,"\n");
    }
  
    fclose(fp);
    printf("\n Results in file ttred2.lst\n\n");
      
    free(vmblock);
 
}

// end of file ttred2.cpp