'************************************************************** '* 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.: "Algebre - 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 * '* * '* Eigenvalue: 0.999999999638076 * '* * '* Eigenvector: * '* -3.09880974276851D-11 * '* 3.098809742764874D-11 * '* 0.7499999999250131 * '* 1 * '* * '* Number of iterations: 32 * '* * '* ---------------------------------------------------------- * '* Subroutines used: * '* * '* 500 ReadData: Read matrix and some parameters * '* 1000 PWM : Calculate greatest eigenvalue and associated* '* eigenvector by the power method * '* 2000 MATINV : Solve a linear system AX=B by Gauss-Jordan * '* (see also program sysmat.bas). * '* 3000 PWIMGT : Calculate smallest eigenvalue and associated* '* eigenvector by the Gauss-Jordan and the * '* power methods. * '* * '* English Basic 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 real. '-------------------------------------------------------------- 'Program Test_PWIMGT DEFDBL A-H, O-Z DEFINT I-N NMAX = 25 DIM A(NMAX, NMAX), X(NMAX) CLS PRINT " ----------------------------------------------" PRINT " Calculate the greatest eigenvalue of a real " PRINT " square matrix and the associated eigenvector " PRINT " by the power method. " PRINT " ----------------------------------------------" GOSUB 500 'call Read_data(n,A,dta,eps,m) GOSUB 3000 'call PWIMGT(eps,dta,m,n,A,it,gamma,X) IF it + 1 = 0 THEN PRINT PRINT " No convergence !" ELSEIF it + 1 = 1 THEN PRINT PRINT " Method does not apply." ELSE PRINT PRINT " Eigenvalue: "; xlambda PRINT PRINT " Eigenvector:" FOR I = 1 TO N PRINT " "; X(I) NEXT I END IF PRINT PRINT " Number of iterations: "; l PRINT END 'of main program 500 'Subroutine Read_data(n,A,dta,eps,m) PRINT " Size of matrix (maximum "; NMAX; "): "; : INPUT "", N FOR I = 1 TO N PRINT PRINT " Line "; I FOR J = 1 TO N PRINT " Element "; J; : INPUT "", A(I, J) NEXT J NEXT I PRINT INPUT " Precision: ", dta INPUT " Epsilon : ", eps INPUT " Maximum number of iterations: ", MAXIT RETURN '************************************************************ '* 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 * '************************************************************ 1000 'Subroutine PWM(eps,dta,m,n,A,it,gamma,X) DIM X0(NMAX) FOR I = 1 TO N X0(I) = 1# / SQR(I) NEXT I it = -1: l = 1 1010 gamma = 0# FOR I = 1 TO N X(I) = 0# FOR J = 1 TO N X(I) = X(I) + A(I, J) * X0(J) NEXT J IF ABS(X(I)) > ABS(gamma) THEN gamma = X(I) NEXT I IF ABS(gamma) < eps THEN it = 0 ELSE FOR I = 1 TO N X(I) = X(I) / gamma NEXT I phi = 0# FOR I = 1 TO N s = ABS(X(I) - X0(I)) IF s > phi THEN phi = s NEXT I IF phi < dta THEN it = 1 ELSE FOR I = 1 TO N X0(I) = X(I) NEXT I l = l + 1 END IF END IF IF it = -1 AND l <= MAXIT THEN GOTO 1010 RETURN '******************************************* '* SOLVING A LINEAR MATRIX SYSTEM AX = B * '* with Gauss-Jordan method using full * '* pivoting at each step. During the pro- * '* cess, original A and B matrices are * '* destroyed to spare storage location. * '* --------------------------------------- * '* INPUTS: A MATRIX N*N * '* B MATRIX N*M * '* --------------------------------------- * '* OUTPUS: A INVERSE OF A N*N * '* DET DETERMINANT OF A * '* B SOLUTION MATRIX N*M * '* --------------------------------------- * '* NOTA - If M=0 inversion of A matrix * '* only (B not used here). * '******************************************* 2000 'Subroutine MATINV 'Labels: 100, 200, 300 'Initializations det = 1# FOR I = 1 TO N PC(I) = 0# PL(I) = 0# CS(I) = 0# NEXT I 'Main loop FOR K = 1 TO N 'Searching greatest pivot PV = A(K, K) IK = K: JK = K PAV = ABS(PV) FOR I = K TO N FOR J = K TO N temp = ABS(A(I, J)) IF temp > PAV THEN PV = A(I, J) PAV = ABS(PV) IK = I: JK = J END IF NEXT J NEXT I 'Search terminated, the pivot is in location I=IK, J=JK. 'Memorizing pivot location: PC(K) = JK PL(K) = IK 'DETERMINANT DET is actualised 'If DET=0, ERROR MESSAGE and STOP 'Machine dependant EPSMACH equals here 2E-16 IF IK <> K THEN det = -det END IF IF JK <> K THEN det = -det END IF det = det * PV temp = ABS(det) IF temp < EPSMACH THEN PRINT #2, PRINT #2, " The determinant equals ZERO !!!" END 'stop program END IF 'POSITIONNING PIVOT IN K,K:} IF IK <> K THEN FOR I = 1 TO N 'EXCHANGE LINES IK and K of matrix A:} TT = A(IK, I) A(IK, I) = A(K, I) A(K, I) = TT NEXT I END IF IF M <> 0 THEN FOR I = 1 TO M TT = B(IK, I) B(IK, I) = B(K, I) B(K, I) = TT NEXT I END IF 'PIVOT is at correct line IF JK <> K THEN FOR I = 1 TO N 'EXCHANGE COLUMNS JK and K of matrix A: TT = A(I, JK) A(I, JK) = A(I, K) A(I, K) = TT NEXT I END IF 'The PIVOT is at correct column. 'and is located in K,K. 'Column K of matrix A is stored in CS vector 'then column K is set to zero. FOR I = 1 TO N CS(I) = A(I, K) A(I, K) = 0# NEXT I CS(K) = 0# A(K, K) = 1# 'Line K of matrix A is modified: temp = ABS(PV) IF temp < EPSMACH THEN PRINT #2, PRINT #2, " PIVOT TOO SMALL - STOP" END 'stop program END IF FOR I = 1 TO N A(K, I) = A(K, I) / PV NEXT I IF M <> 0 THEN FOR I = 1 TO M B(K, I) = B(K, I) / PV NEXT I END IF 'Other lines of matrix A are modified: FOR J = 1 TO N IF J = K THEN GOTO 100 END IF FOR I = 1 TO N 'Line J of matrix A is modified: A(J, I) = A(J, I) - CS(J) * A(K, I) NEXT I IF M <> 0 THEN FOR I = 1 TO M B(J, I) = B(J, I) - CS(J) * B(K, I) NEXT I END IF 100 NEXT J 'Line K is ready NEXT K 'end K loop 'MATRIX A INVERSION IS DONE - UNSCRAMBLE OF MATRIX A 'EXCHANGE LINES FOR I = N TO 1 STEP -1 IK = INT(PC(I)) IF IK = I THEN GOTO 200 END IF 'EXCHANGE LINES I and PC(I) of matrix A: FOR J = 1 TO N TT = A(I, J) A(I, J) = A(IK, J) A(IK, J) = TT NEXT J IF M <> 0 THEN FOR J = 1 TO M TT = B(I, J) B(I, J) = B(IK, J) B(IK, J) = TT NEXT J END IF 'NO MORE EXCHANGE is NECESSARY 'Go to next line. 200 NEXT I 'End of loop i=N to 1 step -1 'EXCHANGE COLUMNS FOR J = N TO 1 STEP -1 JK = INT(PL(J)) IF JK = J THEN GOTO 300 END IF 'EXCHANGE COLUMNS J ET PL(J) of matrix A: FOR I = 1 TO N TT = A(I, J) A(I, J) = A(I, JK) A(I, JK) = TT NEXT I 'NO MORE EXCHANGE is NECESSARY 'Go to next column. 300 NEXT J ' UNSCRAMBLE TERMINATED RETURN 'End of subroutine 1000 '************************************************************ '* 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 * '************************************************************ 3000 'Subroutine PWIMGT(eps,dta,m,n,A,it,lambda,X) ' double xlambda,det,gamma M = 0 GOSUB 2000 'call MATINV(N,M,A,X,det) inverse of A is now in A PRINT PRINT " DET="; det IF ABS(det) > eps THEN GOSUB 1000 'call PWM(eps,dta,m,n,A,it,gamma,X) IF it = 1 THEN xlambda = 1# / gamma END IF RETURN 'end of file tpwimgt.bas