/************************************************************* * 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 #include #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 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=0; I--) { IK=(int) PC[I]; if (IK==I) goto fin1; // EXCHANGE LINES I and PC(I) of AA : for (J=0; J=0; J--) { JK=(int) PL[J]; if (JK==J) goto fin2; // EXCHANGE COLUMNS J and PL(J) of AA : for (I=0; I0; 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