/*******************************************************
* Eigenvalues and eigenvectors of a real square matrix *
* by Rutishauser's method and inverse iteration method *
* ---------------------------------------------------- *
* Reference:                                           *
*                                                      *
*   "ALGEBRE Algorithmes et programmes en Pascal       *
*    de Jean-Louis Jardrin - Dunod BO-PRE 1988."       *
*    [BIBLI 10].                                       *
*                                                      *
*                    C++ Release By J-P Moreau, Paris. *
*                            (www.jpmoreau.fr)         *
* ---------------------------------------------------- *
* SAMPLE RUN:                                          *
*                                                      *
* Input file: elpro.dat                                *
*                                                      *
* 5                                                    *
*   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                                *
*                                                      *
* Output to screen:                                    *
*                                                      *
*        *** EIGENVALUES AND EIGENVECTORS ***          *
*              OF A REAL SQUARE MATRIX                 *
*              BY RUTISHAUSER'S METHOD                 *
*            AND INVERSE ITERATION METHOD              *
*                                                      *
*  Matrix A:                                           *
*                                                      *
*  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    *
*                                                      *
*                                                      *
*  Eigenvalue 1: 23.755955                             *
*                                                      *
*  Eigenvector:                                        *
*                                                      *
*  0.705820 -0.012677  0.131486 -0.527499  1.000000    *
*                                                      *
*  Eigenvalue 2:-10.486545                             *
*                                                      *
*  Eigenvector:                                        *
*                                                      *
*  0.467172 -0.007947 -0.507827  1.000000  0.264432    *
*                                                      *
*  Eigenvalue 3: 0.463350                              *
*                                                      *
*  Eigenvector:                                        *
*                                                      *
*  0.208812  1.000000 -0.478027 -0.275000 -0.216915    *
*                                                      *
*  Eigenvalue 4:-7.774580                              *
*                                                      *
*  Eigenvector:                                        *
*                                                      *
*  1.000000 -0.326100  0.209620 -0.147689 -0.815422    *
*                                                      *
*  Eigenvalue 5: 0.463350                              *
*                                                      *
*  Eigenvector:                                        *
*                                                      *
*  0.208812  1.000000 -0.478027 -0.275000 -0.216915    *
*                                                      * 
*******************************************************/
#include <stdio.h>
#include <math.h>

#define  NMAX  20

typedef double MAT[NMAX][NMAX];
typedef double VEC[NMAX];

  FILE *fp;

  void Data(int *n, MAT A, double *eps,double *d1,double *d2,int *m) {
    int i,j;
    
	fp = fopen("elpro6.dat","r");

    fscanf(fp, "%d", n);
    printf("\n  Matrix A:\n\n");
    for (i=1; i<=*n; i++) {
	  for (j=1; j<=*n; j++) { 
        fscanf(fp,"%lf",  &A[i][j]);
        printf(" %10.6f", A[i][j]);
      }
      printf("\n");
    }
    printf("\n");;
    fclose(fp);
    *eps=1e-10;
    *d1=1e-8;
    *d2=1e-8;
    *m=600;
  }

  /**************************************************************
  * Procedure DECCRM determines the lower triangular matrix and *
  * the upper triangukar matrix of Crout's decomposition of a   *
  * given square real matrix, A.                                *
  * ----------------------------------------------------------- *
  * INPUTS:                                                     *
  *         eps: required precision (double)                    *
  *          n : size of matrix A (integer)                     *
  *          A : input matrix (n x n)                           *
  * OUTPUTS:                                                    *
  *          it: flag, =0 if method does not apply              *
  *                    =1 if method is ok.                      *
  *           U: lower triangular matrix.                       *
  *           V: upper triangular matrix.                       *
  **************************************************************/
  void DECCRM(double eps, int n, MAT A, int *it, MAT U, MAT V) {
    int  i,j,k;
    double s;

    if (fabs(A[1][1]) < eps)  *it=0;
    else {
      for (i=1; i<=n; i++)  U[i][1]=A[i][1];
      V[1][1]=1.0;
      for (j=2; j<=n; j++)  V[1][j]=A[1][j]/U[1][1];
      *it=1; k=2;
      while (*it!=0 && k<=n) {
	    for (i=1; i<=n; i++)
	      if (i < k)  
		    U[i][k]=0.0;
		  else {
	        s=0.0;
	        for (j=1; j<k; j++)  s += U[i][j]*V[j][k];
	        U[i][k]=A[i][k] - s;
		  }
	    if (fabs(U[k][k]) < eps) 
          *it=0;
		else {
	      for (j=1; j<=n; j++)
	        if (j < k) 
			  V[k][j]=0.0;
	        else if (j==k) 
			  V[k][j]=1.0;
            else {
	          s=0.0;
	          for (i=1; i<k; i++)  s += U[k][i]*V[i][j];
	          V[k][j]=A[k][j]/U[k][k];
			}
	        k++;
		}
      } //while
    }
  }

  void MatPrint(char *title, MAT A, int n) {
    int i,j;
    printf("\n %s\n", title);
	for (i=1; i<=n; i++) {
      for (j=1; j<=n; j++)
        printf(" %f", A[i][j]);
	  printf("\n");
    }
  }

  void VecPrint(char *title, VEC X, int n) {
    int i;
    printf("\n %s\n", title);
	for (i=1; i<=n; i++) printf(" %f", X[i]);
	printf("\n");
  }

  /********************************************************
  * Calculate the eigenvalues of a real square matrix by  *
  * Rutishauser's Method.                                 *
  * ----------------------------------------------------- *
  * INPUTS:                                               *
  *        eps: absolute precision (double)               *
  *        dta: relative precision (double)               *
  *         m : maximum number of iterations (integer)    *
  *         n : size of matrix A (integer)                *
  *         A : input real square matrix (n x n)          *
  * OUTPUTS:                                              *
  *         it: flag, =-1 if convergence is not obtained  *
  *                   =1 if convergence is ok.            *
  *         R : contains in output the n eigenvalues of A *
  *                                                       *         
  ********************************************************/
  void VAMR(double eps, double dta, int m, int n, MAT A, int *it, VEC R) {
    int i,j,k,l;
    double phi,s,t0;
    MAT U, V;

    t0=0.0;
    l=1;
    while (l<=m && *it!=1) {
      for (i=1; i<=n; i++)  R[i]=A[i][i];

      DECCRM(eps,n,A,it,U,V);

      if (*it==0) {
	    for (i=1; i<=n; i++)  A[i][i]=A[i][i]+1.0;
	    t0 += 1.0;
      }
      else {
	    for (i=1; i<=n; i++)
		  for (j=1; j<=n; j++) {
	        s=0.0;
	        for (k=1; k<=n; k++)  s += V[i][k]*U[k][j];
	        A[i][j]=s;
		  }
	    phi=0.0;
        for (i=1; i<=n; i++) {
	      s=fabs(A[i][i]-R[i]);
	      if (s > phi)  phi=s;
		}
	    if (phi < dta)
	      for (i=1; i<=n; i++)  R[i]=A[i][i]-t0;
		else {
	      l++;
	      *it=-1;
		}
      }
    } //while
  }

  /***********************************************************
  * Procedure IIM calculates a real eigenvalue and the asso- *
  * ciated eigenvector of a real square matrix the inverse   *
  * iteration method.                                        *
  * -------------------------------------------------------- *
  * INPUTS:                                                  *
  *         eps : absolute precision (double)                *
  *         dta : relative precision (double)                *
  *          m  : maximum number of iterations (integer)     *
  *          n  : size of matrix A                           *
  *          A  : input real square matrix (n x n)           *
  * OUTPUTS:                                                 *
  *          it : flag, =-1 if convergence is not obtained   *
  *                     =1 if convergence is ok.             *
  *        Gamma: starting value for the eigenvalue as input *
  *               approximation of the eigenvalue with preci-*
  *               sion dta in output.                        *
  *          X1 : contains in output the associated eigen-   *
  *               vector.                                    *
  *                                                          *
  ***********************************************************/
  void IIM(double eps, double dta, int m, int n, MAT A, int *it,
	       double *gamma, VEC X1)  {
    int i,j,k,l,l0;
    double  p0,phi,s,t0;
    VEC W, X0;
    int LP[NMAX];

    for (i=1; i<=n; i++)  A[i][i]=A[i][i]-(*gamma);
    for (k=1; k<n; k++) {
      p0=A[k][k]; l0=k;
      for (i=k+1; i<=n; i++)
		if (fabs(A[i][k]) > fabs(p0)) {
	      p0=A[i][k]; l0=i;
		}
        LP[k]=l0;
        if (fabs(p0) < eps) {
	      p0=eps; A[l0][k]=eps;
        }
        if (l0 != k)
		  for (j=k; j<=n; j++) {
	        t0=A[k][j]; A[k][j]=A[l0][j]; A[l0][j]=t0;
		  }
	  for (i=k+1; i<=n; i++) {
	    A[i][k]=A[i][k]/p0;
	    for (j=k+1; j<=n; j++)
	      A[i][j]=A[i][j]-A[i][k]*A[k][j];
      }
	} //end k loop

    if (fabs(A[n][n]) < eps)  A[n][n]=eps;
    for (i=1; i<=n; i++)  X0[i]=1.0/sqrt(1.0*i);
    *it=-1; l=1;
    while (*it==-1 && l<=m) {
      for (i=1; i<=n; i++) W[i]=X0[i];
      for (k=1; k<n; k++) {
        l0=LP[k];
        if (l0!=k) {
          t0=W[k]; W[k]=W[l0]; W[l0]=t0;
        }
        for (i=k+1; i<=n; i++)  W[i]=W[i]-A[i][k]*W[k];
      }
      X1[n]=W[n]/A[n][n];
      for (i=n-1; i>0; i--) {
        s=0.0;
        for (j=i+1; j<=n; j++)  s += A[i][j]*X1[j];
	    X1[i]=(W[i]-s)/A[i][i];
      }
      p0=0.0;
      for (i=1; i<=n; i++)
        if (fabs(X1[i]) > fabs(p0))  p0=X1[i];
      for (i=1; i<=n; i++)
        X1[i]=X1[i]/p0;
      phi=0.0;
      for (i=1; i<=n; i++) {
        s=fabs(X1[i]-X0[i]);
        if (s > phi)  phi=s;
      }

      if (phi < dta) {
        *gamma = *gamma + 1.0/p0;
        *it=1;
      }
      else {
        for (i=1; i<=n; i++)  X0[i]=X1[i];
        l++;
      }
	} // while
  }    	  						    


  /*****************************************************
  * INPUTS:                                            *
  * EPS : precision (Double)                           *
  * D1  : precision d1 (Double)                        *
  * D2  : precision d2 (Double)                        *
  * M   : maximum number of iterations (integer)       *
  * N   : order of matrix A (integer)                  *
  * A   : input matrix to study (of MAT type)          *
  * -------------------------------------------------- *
  * OUTPUTS:                                           *
  * IT  : -1 if no convergence is obtained (integer)   *
  * R   : table of eigenvalues (of VEC type)           *
  * VX  : table of eigenvectors (of MAT type)          *
  *****************************************************/
  void EPMRI(double eps, double d1, double d2, int m, int n, MAT A, int *it,
	         VEC R, MAT VX)  {
    int i,j,k; 
    VEC X;
	MAT A1;

    for (i=1; i<=n; i++)
      for (j=1; j<=n; j++)
        A1[i][j] = A[i][j];

    VAMR(eps,d2,m,n,A1,it,R);

// restore A1 after VAMR
    for (i=1; i<=n; i++)
      for (j=1; j<=n; j++)
        A1[i][j] = A[i][j];

    j=1;
    while (*it==1 && j<=n) {
      IIM(eps,d1,m,n,A1,it,&R[j],X);
  
// restore A1 after IIM
      for (i=1; i<=n; i++) {
        for (k=1; k<=n; k++)  
		  A1[i][k] = A[i][k];

        VX[i][j]=X[i];
      }
      j++;     
    }
  }


void main()  {

    int i,it,j,m, n;
    double d1,d2,eps;
    MAT A, VX;
    VEC R;

    printf("\n");;
    printf("    *** EIGENVALUES AND EIGENVECTORS ***\n");
    printf("          OF A REAL SQUARE MATRIX\n");
    printf("          BY RUTISHAUSER'S METHOD\n");
    printf("        AND INVERSE ITERATION METHOD\n\n");

    // read data from input text file
    Data(&n, A, &eps, &d1, &d2, &m);

    EPMRI(eps,d1,d2,m,n,A,&it,R,VX);

    if (it==-1)
      printf("  No convergence!\n\n");
    else
	  for (j=1; j<=n; j++) {
	    printf("\n  Eigenvalue %d: %10.6f\n", j, R[j]);
	    printf("\n  Eigenvector:\n\n");
	    for (i=1; i<=n; i++)  
		  printf(" %10.6f",VX[i][j]);
        printf("\n");
		printf("\n");
      }
    printf("\n");  
}

// End of file elpro.cpp