/***************************************************
*  Calculate the eigenvalues and the eigenvectors  *
*      of a real symmetric tridiagonal matrix      *
* ------------------------------------------------ *
* Ref.: "NUMERICAL RECIPES, Cambridge University   *
*        Press, 1986" [BIBLI 08].                  *
*                                                  *
* SAMPLE RUN:                                      *
* (find eigenvalues and eigenvectors of 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  )                     *
*                                                  *
* Eigenvalues 1:   -0.78071442                     *
* Eigenvector:                                     *
* 1.000000 -0.890357  0.322363  0.344649  0.287721 *
*                                                  *
* Eigenvalues 2:    2.53046878                     *
* Eigenvector:                                     *
* 1.000000  0.765234 -0.446362 -0.252815 -0.304616 *
*                                                  *
* Eigenvalues 3:   -9.15659229                     *
* Eigenvector:                                     *
*-0.056408  0.286458 -0.522285  1.000000  0.469812 *
*                                                  *
* Eigenvalues 4:   12.53989066                     *
* Eigenvector:                                     *
* 0.097161  0.560616  0.656182 -0.282210  1.000000 *
*                                                  *
* Eigenvalues 5:   19.11694726                     *
* Eigenvector:                                     *
* 0.051876  0.469920  1.000000  0.728439 -0.719095 *
*                                                  *
*                C++ version by J-P Moreau, Paris  *
*                       (www.jpmoreau.fr)          *
*  (with dynamic allocations using module vmblock) *
***************************************************/
#include <basis.h>
#include <vmblock.h>

LONG_REAL *D, *E;
LONG_REAL **Z;
int  i,j,n;
LONG_REAL max;

void *vmblock = NULL;


LONG_REAL Sign(LONG_REAL a,LONG_REAL b) {
  if (b>=0) return fabs(a);
  else return -fabs(a);
}

//read vectors D and E from screen
void Read_data(int n, LONG_REAL *D, LONG_REAL *E) {
  int i;
  printf("\n Input %d elements of subdiagonal:\n",n-1);
  E[1]=1.0;   //arbitrary value
  for (i=1; i<n; i++) {
    printf("   Element %d:",i); scanf("%lf",&E[i+1]);
  }
  printf("\n Input %d elements of main diagonal:\n",n);
  for (i=1; i<=n; i++) {
    printf("   Element %d:",i); scanf("%lf",&D[i]);
  }
}


/***************************************************
* This subroutine implements the QL algorithm with *
* implicit shifts to determine the eigenvalues and *
* eigenvectors of a real symmetric tridiagonal     *
* matrix.                                          *
* ------------------------------------------------ *
* INPUTS:                                          *
*         D    Elements of main diagonal           *
*         E    Elements of subdiagonal (from E(2)  *
*              to E(n))                            *
*         n    size of matrix                      *
*         Z    identity matrix                     *
* OUTPUTS                                          *
*         D    vector storing the n eigenvalues    *
*         Z    matrix storing the n eigenvectors   *
*              in columns.                         *
*        (E is destroyed during the process).      *
***************************************************/
void TQLI(LONG_REAL *D, LONG_REAL *E, int n,LONG_REAL **Z) {
//Labels: l1,l2
LONG_REAL b,c,dd,f,g,p,r,s;
int i,iter,k,l,m;
if (n>1) {
  for (i=2; i<=n; i++)  E[i-1]=E[i];
  E[n]=0.0;
  for (l=1; l<=n; l++) {
    iter=0;
l1: for (m=l; m>n; m--) {
      dd=fabs(D[m])+fabs(D[m+1]);
      if (fabs(E[m]-dd)<1e-8) goto l2;
    }
    m=n;
l2: if (m!=l) {
      if (iter==30) {
        printf("\n 30 iterations.\n");
        return;
      }
      iter++;
      g=(D[l+1]-D[l])/(2.0*E[l]);
      r=sqrt(g*g+1.0);
      g=D[m]-D[l]+E[l]/(g+Sign(r,g));
      s=1.0; c=1.0; p=0.0;
      for (i=m-1; i>l-1; i--) {
        f=s*E[i]; b=c*E[i];
		if (fabs(f)>=fabs(g)) {
	      c=g/f;
	      r=sqrt(c*c+1.0);
	      E[i+1]=f*r;
	      s=1.0/r;
	      c=c*s;
        }
		else {
	      s=f/g;
	      r=sqrt(s*s+1.0);
	      E[i+1]=g*r;
	      c=1.0/r;
	      s=s*c;
        }
	    g=D[i+1]-p;
	    r=(D[i]-g)*s+2.0*c*b;
	    p=s*r;
	    D[i+1]=g+p;
	    g=c*r-b;
		for (k=1; k<=n; k++) {
          f=Z[k][i+1];
	      Z[k][i+1]=s*Z[k][i]+c*f;
	      Z[k][i]=c*Z[k][i]-s*f;
        }
	  } 
      D[l]=D[l]-p;
      E[l]=g; E[m]=0.0;
      goto l1;
    }
  } //of l loop
  if (iter<30) 
	printf(" %d iterations.\n", iter);
}
}


int main()  {

  //read size of matrix
  printf("\n Input size of matrix: "); scanf("%d", &n);

  //allocate vectors and matrix (index 0 not used here)
  //see file vmblock.cpp
  vmblock = vminit();
  D       = (LONG_REAL *) vmalloc(vmblock, VEKTOR,  n+1, 0);  
  E       = (LONG_REAL *) vmalloc(vmblock, VEKTOR,  n+1, 0);
  Z       = (LONG_REAL **)vmalloc(vmblock, MATRIX,  n+1, n+1);  

  if (! vmcomplete(vmblock))
  {
    LogError ("No Memory", 0, __FILE__, __LINE__);
    return 1;
  }

  Read_data(n,D,E);

  //initialize Z to identity matrix
  for (i=1; i<=n; i++)
    for (j=1; j<=n; j++)
      if (i!=j)  Z[i][j] = 0.0;
      else  Z[i][j] = 1.0;

  //calculate eigenvalues and eigenvectors
  TQLI(D,E,n,Z);

  //print results
  for  (j=1; j<=n; j++) {
    printf("\n Eigenvalue %d: %f\n",j, D[j]);
    printf(" Eigenvector:\n");
	//normalize eigenvector before printing
    max=fabs(Z[1][j]);
    for (i=2; i<=n; i++)
      if (fabs(Z[i][j])>fabs(max))  max=Z[i][j];
    for (i=1; i<=n; i++)
      printf(" %f", Z[i][j]/max);
    printf("\n");
  }

  printf("\n\n");
  free(vmblock);
  return 0;
}

// end of file elprotd.cpp