/*************************************************************
* Calculate the smallest eigenvalue lambda of a real square  *
* matrix and the associated eigenvector. The method consists *
* to inverse the matrix by a classical Gauss-Jordan method,  *
* then to calculate the greatest eigenvalue gamma of the in- *
* verse matrix by the power method, and then we have:        *
* lambda=1/gamma (assuming gamma <> 0).                      * 
* ---------------------------------------------------------- *
* Ref.: "Algèbre - Algorithmes et programmes en Pascal By    *
*        Jean-Louis Jardrin, Dunod Editor, Paris, 1988"      *
*        [BIBLI 10].                                         *
* ---------------------------------------------------------- *
* SAMPLE RUN:                                                *
*                                                            *
* ----------------------------------------------             *
*  Calculate the smallest eigenvalue of a real               *
*  square matrix and the associated eigenvector              *
*  by inversion and 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: 32                          *
*                                                            *
*    DET= -20.000000                                         *
*                                                            *
*    Eigenvalue: 1.000000                                    *
*                                                            *
*    Eigenvector:                                            *
*     0.000000                                               *
*     0.000000                                               *
*     0.750000                                               *
*     1.000000                                               *            
*                                                            *
*                 English C++ version by J-P Moreau, Paris.  *
*                             (www.jpmoreau.fr)              *
**************************************************************
  Exact values are: lambda = 1
                    eigenvector = (0,0,0.75,1)
  NOTE: the method fails if the result is an opposite value
        of another eigenvalue or if the matrix A is singular
        (not inversible) or if lambda is not double.
------------------------------------------------------------*/
#include <stdio.h>
#include <math.h>

#define   NMAX  26
#define   MACH_EPS  1.2e-16

typedef   double MAT[NMAX][NMAX];

int    i,it,m,n;
double dta,eps,lambda;
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()                                                                                                     


/*****************************************************************************
*  SOLVING A LINEAR MATRIX SYSTEM AX = B with Gauss-Jordan method using full *
*  pivoting at each step. During the process, original AA and BB matrices are*
* destroyed to spare storage location.                                       *
* -------------------------------------------------------------------------- *
* INPUTS:                                                                    *
*            N    : size of matrix AA                                        *
*            M    : number of simultaneous second members                    *
*            AA   : MATRIX N*N of type MAT                                   *
*            BB   : MATRIX N*M of type MAT (not used here)                   *
* -------------------------------------------------------------------------- *
* OUTPUTS:   AA   : INVERSE OF AA of type MAT                                *
*            DET  : DETERMINANT OF AA (Double)                               *
*            BB   : SOLUTION MATRIX N*M of type MAT (not used here)          *
* -------------------------------------------------------------------------- *
* NOTA:      If M=0 inversion of AA matrix only (BB is not used here).       *
*                                                                            *
*                               C++ english version by J-P Moreau, Paris.    *
*****************************************************************************/
void MATINV(int N,int M,MAT AA,MAT BB,double *DET)   {                                           
      double PC[NMAX],PL[NMAX],CS[NMAX];                                                  
      double PV,PAV,temp,TT;   
      int I,IK,J,JK,K;
//Put AA starting at index zero
	  for (I=0; I<N; I++)
        for (J=0; J<N; J++)
			A[I][J] = A[I+1][J+1]; 
// INITIALIZATIONS :                                                             
      *DET= 1.0; 
      for (I=0; I<N; I++)  {                                                                
        PC[I]= 0.0;                                                                
        PL[I]= 0.0;                                                                
        CS[I]= 0.0;                                                                
      }              
                                                                       
      for (K=0; K<N; K++) {                                                                  
// SEARCH FOR GREATEST PIVOT :                                               
        PV=AA[K][K];                                                              
        IK=K;                                                                    
        JK=K;                                                                    
        PAV= fabs(PV);                                                            
        for (I=K; I<N; I++)                                                                
          for (J=K; J<N; J++)  {     
            temp = fabs(AA[I][J]);                                                        
            if (temp > PAV) {                                      
              PV=AA[I][J];                                                        
              PAV= fabs(PV);                                                      
              IK=I;                                                              
              JK=J;                                                              
            }                                                               
          }                                                                 
                                                                               
// SEARCH TERMINATED, PIVOT PV IS AT LOCATION I=IK,J=JK.                    
                                                                               
// STORE LOCATION OF PIVOT :                                        
        PC[K]=JK;                                                                
        PL[K]=IK;                                                                
                                                                               
// UPDATE DETERMINANT DET                                              
// IF DET=0, ERROR MESSAGE AND RETURN                                           
// EPSMACH IS MACHINE DEPENDENT                                         
// WE TAKE HERE 1.2e-16 (CF. #define) :                                        
                                                                               
        if (IK!=K) *DET=-*DET;                                                   
        if (JK!=K) *DET=-*DET;                                                   
        *DET=*DET*PV;  
        temp= fabs(*DET);                                                            
        if (temp < MACH_EPS) {                                          
// Error MESSAGE and return                                                     
          printf("\n  THE DETERMINANT EQUALS ZERO, NO SOLUTION !!!\n");                                                              
          return;                                                                  
        }                                                                   
                                                                               
// LOCATE PIVOT IN K,K :                                              
        if(IK!=K)                                                         
          for (I=0; I<N; I++) {                                                              
// EXCHANGE LINES IK AND K OF AA :                                            
            TT=AA[IK][I];                                                         
            AA[IK][I]=AA[K][I];                                                    
            AA[K][I]=TT;                                                          
          }                                                                 
                                                                           
        if(M!=0)                                                           
          for (I=0; I<M; I++) {                                                               
            TT=BB[IK][I];                                                           
            BB[IK][I]=BB[K][I];                                                      
            BB[K][I]=TT;                                                            
          }                                                                   
                                                                           
// PIVOT IS LOCATED AT CORRECT LINE.                                                
        if(JK!=K)                                                         
          for (I=0; I<N; I++) {                                                              
// EXCHANGE COLUMNS JK AND K OF AA :                                          
            TT=AA[I][JK];                                                         
            AA[I][JK]=AA[I][K];                                                    
            AA[I][K]=TT;                                                          
          }                                                                 
                                                                           
// PIVOT IS LOCATED AT CORRECT COLUMN.                                              
// PIVOT IS NOW IN K,K.                                               
                                                                               
// STORE COLUMN K OF AA IN VECTOR CS                             
// THEN SET COLUMN TO ZERO  :                                             
        for (I=0; I<N; I++) {                                                                
          CS[I]=AA[I][K];                                                         
          AA[I][K]= 0.0;                                                          
        }                                                                   
                                                                               
        CS[K]= 0.0;                                                                
        AA[K][K]= 1.0;                                                              
// MODIFY LINE K OF AA :  
        temp= fabs(PV);                                          
        if(temp < MACH_EPS) {                                            
          printf("\n  PIVOT TROP PETIT - ARRET\n");                               
          return;                                                                  
        }                                                                   
        for (I=0; I<N; I++)                                                                
          AA[K][I]=AA[K][I]/PV;                                                    
                                                                           
        if (M!=0)                                                         
          for (I=0; I<M; I++)                                                             
            BB[K][I]=BB[K][I]/PV;                                                  
                                                                           
// MODIFY OTHER LINES OF AA :                                        
        for (J=0; J<N; J++) {                                                                
          if (J==K) J++;                                                  
          for (I=0; I<N; I++)                                                              
// MODIFY LINE J OF AA :                                            
            AA[J][I]=AA[J][I]-CS[J]*AA[K][I];
          if (M!=0)                                                       
            for (I=0; I<M; I++)                                                          
              BB[J][I]=BB[J][I]-CS[J]*BB[K][I];                                     
        }                                                                   
// LINE K IS FINISHED                                                
      } //  end of loop K                                                              
                                                                               
// MATRIX AA IS INVERSED - UNSCRAMBLE MATRIX                         
                                                                               
// EXCHANGE LINES                                                            
      for (I=N-1; I>=0; I--) {                                                               
        IK=(int) PC[I];                                                                
        if (IK==I) goto fin1;                                                   
// EXCHANGE LINES I and PC(I) of AA :                                         
        for (J=0; J<N; J++) {                                                                
          TT=AA[I][J];                                                            
          AA[I][J]=AA[IK][J];                                                      
          AA[IK][J]=TT;                                                           
        }
        if (M!=0)                                                         
          for (J=0; J<M; J++) {                                                             
            TT=BB[I][J];                                                          
            BB[I][J]=BB[IK][J];                                                    
            BB[IK][J]=TT;                                                         
          }                                                                 
        fin1: ;                                                                 
// No more EXCHANGE is NECESSARY                                                      
// Go to next line                                                  
      } // end of i loop                                                                     
                                                                               
// EXCHANGE COLUMNS                                                          
      
      for (J=N-1; J>=0; J--) {                                                                         
        JK=(int) PL[J];                                                                
        if (JK==J) goto fin2;                                                   
// EXCHANGE COLUMNS J and PL(J) of AA :                                       
        for (I=0; I<N; I++) {                                                                
          TT=AA[I][J];                                                            
          AA[I][J]=AA[I][JK];                                                      
          AA[I][JK]=TT;                                                           
        } 
        fin2: ;                                                                  
// No more exchange is necessary,                                           
// go to next column
      }                                                                     
// UNSCRAMBLE TERMINATED
	  
// Put AA starting with index 1
	  for (I=N; I>0; I--)
        for (J=N; J>0; J--)
          A[I][J] = A[I-1][J-1];
      return;  
}               

//input data from screen
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);
}


/***********************************************************
* calculate smallest eigenvalue and associated eigenvector *
* by the Gauss-Jordan and the power methods                *
* -------------------------------------------------------- *
* 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.                        *
*        lambda : smallest eigenvalue (in absolute value)  *
*                 of input matrix A(n,n)                   *
*         X     : associated eigenvector                   *
***********************************************************/
void PWIMGT(double eps,double dta, int m, int n, MAT A, int *it,
			 double *lambda, double *X)  {
  double det,gamma;
  MAT AM1;  // dummy matrix for MATINV
  MATINV(n,0,A,AM1,&det);  // inverse of A is now in A
  printf("\n  DET=%f\n", det);
  if (fabs(det)>eps)  {
    PWM(eps,dta,m,n,A,it,&gamma,X);
    if (*it==1)  *lambda=1.0/gamma;
  }
}                                                                                                     


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

  Read_data();

  PWIMGT(eps,dta,m,n,A,&it,&lambda,X);

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

// end of file tpwimgt.cpp