/*******************************************************************************
*     Find Eigenvalues and Eigenvectors of a symmetric tridiagonal matrix      *
*     using QL method                                                          *
* ---------------------------------------------------------------------------- *
* SAMPLE RUNS                                                                  *
*                                                                              *
* Example #1:    (tridiagonal matrix)                                          *
*          1  2  0  0     0                                                    *
*          2  4  7  0     0                                                    *
*     A =  0  7 10  8     0                                                    *
*          0  0  8 -0.75 -9                                                    *
*          0  0  0 -9    10                                                    *
*                                                                              *
*  Output file ttql2.lst contains:                                             *
*                                                                              *
*  EIGENVALUES AND EIGENVECTORS OF A SYMMETRIC                                 *
*  TRIDIAGONAL MATRIX BY QL METHOD.                                            *
*                                                                              *
*  N =           5                                                             *    
*                                                                              *
*  Input tridiagonal matrix:                                                   *
*   1.000000  2.000000  0.000000  0.000000  0.000000                           *
*   2.000000  4.000000  7.000000  0.000000  0.000000                           *
*   0.000000  7.000000 10.000000  8.000000  0.000000                           *
*   0.000000  0.000000  8.000000 -0.750000 -9.000000                           *
*   0.000000  0.000000  0.000000 -9.000000 10.000000                           *
*                                                                              *
*  Error code: 0                                                               *
*                                                                              *
*  Eigenvalues:                                                                *
*   1    -9.15659229                                                           *
*   2    -0.78071442                                                           *
*   3     2.53046878                                                           *
*   4    12.53989066                                                           *
*   5    19.11694726                                                           *
*                                                                              *
*  Eigenvectors (in lines]:                                                    *
*  -0.056408  0.286458 -0.522285  1.000000  0.469812                           *
*   1.000000 -0.890357  0.322363  0.344649  0.287721                           *
*   1.000000  0.765234 -0.446362 -0.252815 -0.304616                           *
*   0.097161  0.560616  0.656182 -0.282210  1.000000                           *
*   0.051876  0.469920  1.000000  0.728439 -0.719095                           *
*                                                                              *
* ---------------------------------------------------------------------------- *
* Reference:   From Numath Library By Tuan Dang Trong in Fortran 77            *
*              [BIBLI 18].                                                     *
*                                                                              *
*                                       C++ Release 1.0 By J-P Moreau, Paris   *
*******************************************************************************/
#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:
!               CALL TQL2(NM,N,D,E,Z,IER)
!     INPUTSS:
!     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 main()  {

    N = 5;

    // 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 main diagonal
    D[1]=ONE; D[2]=4.0; D[3]=10.0; D[4]=-0.75; D[5]=10.0;
    // define subdiagonal
    E[1]=ZERO; E[2]=TWO; E[3]=7.0; E[4]=8.0; E[5]=-9.0;

    // Initialize matrix Z to unity matrix
    for (I=1; I<=N; I++)
      for (J=1; J<=N; J++)
        if (J==I)  Z[I][J] = ONE;
        else  Z[I][J] = ZERO;

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

    fprintf(fp,"\n EIGENVALUES AND EIGENVECTORS OF A SYMMETRIC\n");
    fprintf(fp," TRIDIAGONAL MATRIX BY QL METHOD.\n\n");
    fprintf(fp," N = %d\n\n", N);
    fprintf(fp," Input tridiagonal matrix:\n");
    for (I=1; I<=N; I++) {
	  for (J=1; J<=N; J++) {
	    if (J < I-1) fprintf(fp,"%10.6f", ZERO);
	    if (J == I-1) fprintf(fp,"%10.6f", E[I]);
		if (J == I+1) fprintf(fp,"%10.6f", E[I+1]);
	    if (J == I) fprintf(fp,"%10.6f", D[I]);
	    if (J>I+1 && J<=N) fprintf(fp,"%10.6f", ZERO);
      }
      fprintf(fp,"\n");
    }

    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 ttql2.lst\n\n");
      
    free(vmblock);
 
}

// end of file ttql2.cpp