/**************************************************************
*  Calculate the greatest eigenvalue of a real square matrix  *
*     and the associated eigenvector by the power method      *
* ----------------------------------------------------------- *
* Ref.: "Algèbre - Algorithmes et programmes en Pascal By     *
*        Jean-Louis Jardrin, Dunod Editor, Paris, 1988"       *
*        [BIBLI 10].                                          *
* ----------------------------------------------------------- *
* SAMPLE RUN:                                                 *
*                                                             *
* ----------------------------------------------              *
*  Calculate the greatest eigenvalue of a real                *
*  square matrix and the associated eigenvector               *
*  by the power method.                                       *
* ----------------------------------------------              *
*                                                             *
*  Size of matrix (maximum 25): 4                             *
*                                                             *
*  Line 1                                                     *
*    Element 1: 1                                             *
*    Element 2: 3                                             *
*    Element 3: 0                                             *
*    Element 4: 0                                             *
*                                                             *
*  Line 2                                                     *
*    Element 1: 4                                             *
*    Element 2: 2                                             *
*    Element 3: 0                                             *
*    Element 4: 0                                             *
*                                                             *
*  Line 3                                                     *
*    Element 1: 1                                             *
*    Element 2: -1                                            *
*    Element 3: 5                                             *
*    Element 4: -3                                            *
*                                                             *
*  Line 4                                                     *
*    Element 1: 2                                             *
*    Element 2: 0                                             *
*    Element 3: 4                                             *
*    Element 4: -2                                            *
*                                                             *
*  Precision: 1e-10                                           *
*  Epsilon  : 1e-10                                           *
*  Maximum number of iterations: 27                           *
*                                                             *
*                                                             *
*    Eigenvalue: 5.000000                                     *
*                                                             *
*    Eigenvector:                                             *
*    0.750000                                                 *
*    1.000000                                                 *
*    -0.520833                                                *
*    -0.083333                                                *            
*                                                             *
*                  English C++ version by J-P Moreau, Paris.  *
*                              (www.jpmoreau.fr)              *
***************************************************************
   Exac values are: gamma = 5
                    eigenvector = (1/48)(36,48,-25,-4)
-------------------------------------------------------------*/                                  
#include <stdio.h>
#include <math.h>

#define   NMAX  26

typedef   double MAT[NMAX][NMAX];

int    i,it,m,n;
double dta,eps,gamma;
MAT    A;
double X[NMAX];


/***********************************************************
* calculate greatest eigenvalue and associated eigenvector *
* by the power method                                      *
* -------------------------------------------------------- *
* INPUTS:                                                  *
*         eps   : smallest number in double precision      *
*         dta   : required precision                       *
*         m     : maximum number of iterations             *
*         n     : size of real square matrix A(n,n)        *
*         A     : real square matrix A(n,n)                *
* OUTPUTS:                                                 *
*         it    : error indicator: -1=no convergence,      *
*                 0=method cannot be applied,              *
*                 1=convergence ok.                        *
*         gamma : greatest eigenvalue (in absolute value)  *
*                 of input matrix A(n,n)                   *
*         X1    : associated eigenvector                   *
***********************************************************/
void PWM(double eps,double dta,int m,int n, MAT A, int *it,
		 double *gamma, double *X1)  {
  int i,j,l;
  double phi,s;
  double X0[NMAX];
  for (i=1; i<=n; i++) X0[i]=1.0/sqrt(i);
  *it=-1; l=1;
  while (*it==-1 && l<=m) {
    *gamma=0.0;
    for (i=1; i<=n; i++) {
      X1[i]=0.0;
      for (j=1; j<=n; j++) X1[i] += A[i][j]*X0[j];
      if (fabs(X1[i])>fabs(*gamma))  *gamma=X1[i];
    }
    if (fabs(*gamma) < eps) *it=0;
    else {
      for (i=1; i<=n; i++)  X1[i] /= *gamma;
      phi=0.0;
      for (i=1; i<=n; i++)  {
        s=fabs(X1[i]-X0[i]);
        if (s>phi) phi=s;
      }
      if (phi<dta) *it=1;
      else {
        for (i=1; i<=n; i++) X0[i]=X1[i];
        l++;
      }
    }
  }
} //PWM()                                                                                                     


void Read_data() {
  int i,j;
  printf("\n Size of matrix (maximum %d): ",NMAX-1); scanf("%d", &n);
  for (i=1; i<=n; i++) {
    printf("\n Line %d\n", i);
    for (j=1; j<=n; j++) {
      printf("  Element %d: ", j); scanf("%lf", &A[i][j]);
    }
  }
  printf("\n Precision: "); scanf("%lf", &dta);
  printf(" Epsilon  : "); scanf("%lf", &eps);
  printf(" Maximum number of iterations: "); scanf("%d", &m);
}


void main()  {
  printf(" ----------------------------------------------\n");
  printf("  Calculate the greatest eigenvalue of a real  \n");
  printf("  square matrix and the associated eigenvector \n");
  printf("  by the power method.                         \n");
  printf(" ----------------------------------------------\n");

  Read_data();

  PWM(eps,dta,m,n,A,&it,&gamma,X);

  switch(it+1) {
    case 0: printf("  No convergence !\n"); break;
    case 1: printf("  Method does not apply.\n"); break;
    case 2: printf("\n  Eigenvalue: %f\n\n", gamma);
            printf("  Eigenvector:\n");
            for (i=1; i<=n; i++)  printf("  %f\n", X[i]);
  }
  printf("\n\n");
}

// end of file tpwm.cpp