/************************************************************** * 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 #include #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