{------------------------------------------------------------------------* * * * Turbo Pascal Numerical Methods Toolbox * * Copyright (c) 1986, 87 by Borland International, Inc. * * * * This unit provides the Jacobi PROCEDURE for finding eigenvalues and * * eigenvectors of a real symmetric matrix. Version with pointers for * * use with Turbo Pascal for Windows by J-P Moreau (Paris). * * (www.jpmoreau.fr) * * * *------------------------------------------------------------------------} UNIT Ujacobi; INTERFACE TYPE Float = Double; CONST TNNearlyZero = 1.2E-16; TNArraySize = 50; { Maximum size of matrix } TYPE PTNV = ^TNvector; PTNM = ^TNmatrix; PTNI = ^TNIntVector; TNvector = ARRAY[1..TNArraySize] OF Float; TNmatrix = ARRAY[1..TNArraySize] OF TNvector; TNIntVector = ARRAY[1..TNArraySize] OF integer; VAR Mat : PTNM; { given square symmetric matrix } Eigenvalues : PTNV; { eigenvalues solution ARRAY } Eigenvectors : PTNM; { eigenvectors solution ARRAY Nota 1 - These pointers must be initialized and disposed of after use by the calling PROGRAM. Nota 2 - Also the optional sorting of the eigenvalues must be provided by the calling PROGRAM. The normalization of the eigenvectors to unity is here provided. } PROCEDURE Jacobi(Dimen : integer; MaxIter : integer; Tolerance : Float; VAR Iter : integer; VAR Error : byte); {----------------------------------------------------------------------------} {- -} {- Input: Dimen, Mat, MaxIter, tolerance -} {- Output: Eigenvalues, Eigenvectors, Iter, Error -} {- -} {- Purpose: The eigensystem of a symmetric matrix can be -} {- computed much more simply than the -} {- eigensystem of an arbitrary matrix. The -} {- cyclic Jacobi method is an iterative -} {- technique for approximating the complete -} {- eigensystem of a symmetric matrix to a given -} {- tolerance. The method consists of multiplying -} {- the matrix, Mat, by a series of rotation -} {- matrices, rot[i]. The rotation matrices are -} {- chosen so that the elements of the upper -} {- triangular part of Mat are systematically -} {- annihilated. That is, rot[1] is chosen so -} {- that Mat[1, 1] is identically zero; rot[2] is -} {- chosen so that Mat[1, 2] is identically zero; -} {- etc. Since each operation will probably -} {- change the value of elements annihilated IN -} {- previous operations, the method is iterative. -} {- Eventually, the matrix will be diagonal. The -} {- eigenvalues will be the elements along the -} {- diagonal of the matrix and the eigenvectors -} {- will be the rows of the matrix created by the -} {- product of all the rotation matrices rot[i]. -} {- -} {- The user inputs the matrix, tolerance and maximum -} {- number of iterations. The PROCEDURE returns the -} {- eigenvalues and eigenvectors (or error code) of the -} {- matrix. -} {- -} {- User-Defined Types: TNvector = ARRAY[1..TNArraySize] of real; -} {- TNmatrix = ARRAY[1..TNArraySize] of TNvector; -} {- -} {- Global Variables: Dimen : integer Dimension of square matrix -} {- Mat : TNmatrix Square matrix -} {- MaxIter : integer Maximum number of -} {- Iterations -} {- tolerance : real tolerance IN answer -} {- Eigenvalues : TNvector Eigenvalues of Mat -} {- Eigenvectors : TNmatrix Eigenvectors of Mat -} {- Iter : integer Number of iterations -} {- Error : byte Error code -} {- -} {- Errors: 0: No error -} {- 1: Dimen < 1 -} {- 2: tolerance < TNNearlyZero -} {- 3: MaxIter < 1 -} {- 4: Mat NOT symmetric -} {- 5: Iter > MaxIter -} {- -} {----------------------------------------------------------------------------} IMPLEMENTATION PROCEDURE Jacobi{Dimen : integer; MaxIter : integer; tolerance : Float; VAR Iter : integer; VAR Error : byte}; VAR Row, Column, Diag : integer; SinTheta, CosTheta : Float; SumSquareDiag : Float; Done : boolean; PROCEDURE TestData(Dimen : integer; MaxIter : integer; tolerance : Float; VAR Error : byte); {---------------------------------------------------} {- Input: Dimen, Mat, MaxIter, tolerance -} {- Output: Error -} {- -} {- This PROCEDURE tests the input data for errors. -} {---------------------------------------------------} VAR Row, Column : integer; BEGIN Error := 0; IF Dimen < 1 THEN Error := 1; IF tolerance <= TNNearlyZero THEN Error := 2; IF MaxIter < 1 THEN Error := 3; IF Error = 0 THEN FOR Row := 1 to Dimen - 1 DO FOR Column := Row + 1 to Dimen DO IF ABS(Mat^[Row, Column] - Mat^[Column, Row]) > TNNearlyZero THEN Error := 4 { Matrix NOT symmetric } END; { PROCEDURE TestData } PROCEDURE Initialize(Dimen : integer; VAR Iter : integer); {--------------------------------------------} {- Input: Dimen -} {- Output: Iter, Eigenvectors -} {- -} {- This PROCEDURE initializes Iter to zero -} {- and Eigenvectors to the identity matrix. -} {--------------------------------------------} VAR i,j : integer; BEGIN Iter := 0; FOR i:=1 to Dimen DO FOR j:= 1 to Dimen DO Eigenvectors^[i,j]:=0.0; FOR i := 1 to Dimen DO Eigenvectors^[i, i] := 1.0; END; { PROCEDURE Initialize } PROCEDURE CalculateRotation(RowRow : Float; RowCol : Float; ColCol : Float; VAR SinTheta : Float; VAR CosTheta : Float); {-----------------------------------------------------------} {- Input: RowRow, RowCol, ColCol -} {- Output: SinTheta, CosTheta -} {- -} {- This PROCEDURE calculates the sine and cosine of the -} {- angle Theta through which to rotate the matrix Mat. -} {- Given the tangent of 2-Theta, the tangent of Theta can -} {- be calculated WITH the quadratic formula. The cosine -} {- and sine are easily calculable from the tangent. The -} {- rotation must be such that the Row, Column element is -} {- zero. RowRow is the Row,Row element; RowCol is the -} {- Row,Column element; ColCol is the Column,Column element -} {- of Mat. -} {-----------------------------------------------------------} VAR TangentTwoTheta, TangentTheta, Dummy : Float; BEGIN IF ABS(RowRow - ColCol) > TNNearlyZero THEN BEGIN TangentTwoTheta := (RowRow - ColCol) / (2 * RowCol); Dummy := Sqrt(Sqr(TangentTwoTheta) + 1); IF TangentTwoTheta < 0 THEN { Choose the root nearer to zero } TangentTheta := -TangentTwoTheta - Dummy ELSE TangentTheta := -TangentTwoTheta + Dummy; CosTheta := 1 / Sqrt(1 + Sqr(TangentTheta)); SinTheta := CosTheta * TangentTheta; END ELSE BEGIN CosTheta := Sqrt(1/2); IF RowCol < 0 THEN SinTheta := -Sqrt(1/2) ELSE SinTheta := Sqrt(1/2); END; END; { PROCEDURE CalculateRotation } PROCEDURE RotateMatrix(Dimen : integer; SinTheta : Float; CosTheta : Float; Row : integer; Col : integer); {--------------------------------------------------------------} {- Input: Dimen, SinTheta, CosTheta, Row, Col -} {- Output: Mat -} {- -} {- This PROCEDURE rotates the matrix Mat through an angle -} {- Theta. The rotation matrix is the identity matrix execept -} {- FOR the Row,Row; Row,Col; Col,Col; and Col,Row elements. -} {- The rotation will make the Row,Col element of Mat -} {- to be zero. -} {--------------------------------------------------------------} VAR CosSqr, SinSqr, SinCos : Float; MatRowRow, MatColCol, MatRowCol, MatRowIndex, MatColIndex : Float; Index : integer; BEGIN CosSqr := Sqr(CosTheta); SinSqr := Sqr(SinTheta); SinCos := SinTheta * CosTheta; MatRowRow := Mat^[Row, Row] * CosSqr + 2 * Mat^[Row, Col] * SinCos + Mat^[Col, Col] * SinSqr; MatColCol := Mat^[Row, Row] * SinSqr - 2 * Mat^[Row, Col] * SinCos + Mat^[Col, Col] * CosSqr; MatRowCol := (Mat^[Col, Col] - Mat^[Row, Row]) * SinCos + Mat^[Row, Col] * (CosSqr - SinSqr); FOR Index := 1 to Dimen DO IF NOT(Index IN [Row, Col]) THEN BEGIN MatRowIndex := Mat^[Row, Index] * CosTheta + Mat^[Col, Index] * SinTheta; MatColIndex := -Mat^[Row, Index] * SinTheta + Mat^[Col, Index] * CosTheta; Mat^[Row, Index] := MatRowIndex; Mat^[Index, Row] := MatRowIndex; Mat^[Col, Index] := MatColIndex; Mat^[Index, Col] := MatColIndex; END; Mat^[Row, Row] := MatRowRow; Mat^[Col, Col] := MatColCol; Mat^[Row, Col] := MatRowCol; Mat^[Col, Row] := MatRowCol; END; { PROCEDURE RotateMat^rix } PROCEDURE RotateEigenvectors(Dimen : integer; SinTheta : Float; CosTheta : Float; Row : integer; Col : integer); {--------------------------------------------------------------} {- Input: Dimen, SinTheta, CosTheta, Row, Col -} {- Output: Eigenvectors -} {- -} {- This PROCEDURE rotates the Eigenvectors Matrix through an -} {- angle Theta. The rotation Matrix is the identity Matrix -} {- except FOR the Row,Row; Row,Col; Col,Col; and Col,Row -} {- elements. The Eigenvectors Matrix will be the product of -} {- all the rotation Matrices which operate on Mat. -} {--------------------------------------------------------------} VAR EigenvectorsRowIndex, EigenvectorsColIndex : Float; Index : integer; BEGIN { Transform eigenvector Mat^rix } FOR Index := 1 to Dimen DO BEGIN EigenvectorsRowIndex := CosTheta * Eigenvectors^[Row, Index] + SinTheta * Eigenvectors^[Col, Index]; EigenvectorsColIndex := -SinTheta * Eigenvectors^[Row, Index] + CosTheta * Eigenvectors^[Col, Index]; Eigenvectors^[Row, Index] := EigenvectorsRowIndex; Eigenvectors^[Col, Index] := EigenvectorsColIndex; END; END; { PROCEDURE RotateEigenvectors } PROCEDURE NormalizeEigenvectors(Dimen : integer); {---------------------------------------------------} {- Input: Dimen, Eigenvectors -} {- Output: Eigenvectors -} {- -} {- This PROCEDURE normalizes the eigenvectors so -} {- that the largest element IN each vector is one. -} {---------------------------------------------------} VAR Row : integer; Largest : Float; PROCEDURE FindLargest(Dimen : integer; Eigenvector : TNvector; VAR Largest : Float); {---------------------------------------} {- Input: Dimen, Eigenvectors -} {- Output: Largest -} {- -} {- This PROCEDURE returns the value of -} {- the largest element of the vector. -} {---------------------------------------} VAR Term : integer; BEGIN Largest := Eigenvector[1]; FOR Term := 2 to Dimen DO IF ABS(Eigenvector[Term]) > ABS(Largest) THEN Largest := Eigenvector[Term]; END; { PROCEDURE FindLargest } PROCEDURE DivVecConst(Dimen : integer; VAR ChangingRow : TNvector; Divisor : Float); {--------------------------------------------------------} {- Input: Dimen, ChangingRow -} {- Output: Divisor -} {- -} {- elementary row operation - dividing by a constant -} {--------------------------------------------------------} VAR Term : integer; BEGIN FOR Term := 1 to Dimen DO ChangingRow[Term] := ChangingRow[Term] / Divisor; END; { PROCEDURE DivVecConst } BEGIN { PROCEDURE NormalizeEigenvectors } FOR Row := 1 to Dimen DO BEGIN FindLargest(Dimen, Eigenvectors^[Row], Largest); DivVecConst(Dimen, Eigenvectors^[Row], Largest); END; END; { PROCEDURE NormalizeEigenvectors } BEGIN { PROCEDURE Jacobi } TestData(Dimen, MaxIter, tolerance, Error); IF Error = 0 THEN BEGIN Initialize(Dimen, Iter); REPEAT Iter := Succ(Iter); SumSquareDiag := 0; FOR Diag := 1 to Dimen DO SumSquareDiag := SumSquareDiag + Sqr(Mat^[Diag, Diag]); Done := true; FOR Row := 1 to Dimen - 1 DO FOR Column := Row + 1 to Dimen DO IF ABS(Mat^[Row, Column]) > tolerance * SumSquareDiag THEN BEGIN Done := false; CalculateRotation(Mat^[Row, Row], Mat^[Row, Column], Mat^[Column, Column], SinTheta, CosTheta); RotateMatrix(Dimen, SinTheta, CosTheta, Row, Column); RotateEigenvectors(Dimen, SinTheta, CosTheta, Row, Column); END; UNTIL Done OR (Iter > MaxIter); FOR Diag := 1 to Dimen DO Eigenvalues^[Diag] := Mat^[Diag, Diag]; NormalizeEigenvectors(Dimen); IF Iter > MaxIter THEN Error := 5 END; END; { PROCEDURE Jacobi } END. { Ujacobi } {end of file ujacobi.pas}