NUMAL Section 3.3.1.1.1
BEGIN SECTION : 3.3.1.1.1 (July, 1974)
AUTHORS: T.J.DEKKER AND W.HOFFMANN.
CONTRIBUTORS: W.HOFFMANN, J.G.VERWER.
INSTITUTE: MATHEMATICAL CENTRE.
RECEIVED: 730716.
BRIEF DESCRIPTION:
THIS SECTION CONTAINS FOUR PROCEDURES FOR CALCULATING EIGENVALUES
OR EIGENVECTORS OF A SYMMETRIC TRIDIAGONAL MATRIX.
VALSYMTRI CALCULATES ALL, OR SOME CONSECUTIVE, EIGENVALUES OF A
SYMMETRIC TRIDIAGONAL MATRIX BY MEANS OF LINEAR INTERPOLATION USING
A STURM SEQUENCE;
VECSYMTRI CALCULATES THE CORRESPONDING EIGENVECTORS BY MEANS OF
INVERSE ITERATION.
QRIVALSYMTRI CALCULATES ALL EIGENVALUES OF A SYMMETRIC TRIDIAGONAL
MATRIX BY MEANS OF QR ITERATION;
QRISYMTRI CALCULATES THE EIGENVECTORS AS WELL.
WHEN ALL EIGENVALUES HAVE TO BE CALCULATED, QRIVALSYMTRI IS
PREFERABLE WITH RESPECT TO THE RUNNING TIME; WHEN THE EIGENVECTORS
ALSO HAVE TO BE CALCULATED, INVERSE ITERATION IS PREFERABLE.
KEYWORDS:
EIGENVALUES,
EIGENVECTORS,
TRIDIAGONAL MATRIX,
STURM-SEQUENCE,
INVERSE ITERATION,
QR ITERATION.
SUBSECTION: VALSYMTRI.
CALLING SEQUENCE:
THE HEADING OF THE PROCEDURE IS:
"PROCEDURE" VALSYMTRI(D, BB, N, N1, N2, VAL, EM);
"VALUE" N, N1, N2; "INTEGER" N, N1, N2;
"ARRAY" D, BB, VAL, EM;
"CODE" 34151;
THE MEANING OF THE FORMAL PARAMETERS IS:
N: <ARITHMETIC EXPRESSION>;
THE ORDER OF THE GIVEN MATRIX;
D: <ARRAY IDENTIFIER>;
"ARRAY" D[1:N];
ENTRY: THE MAIN DIAGONAL OF THE SYMMETRIC TRIDIAGONAL
MATRIX;
BB: <ARRAY IDENTIFIER>;
"ARRAY" BB[1:N-1];
ENTRY: THE SQUARES OF THE CODIAGONAL ELEMENTS OF THE
SYMMETRIC TRIDIAGONAL MATRIX;
N1,N2: <ARITHMETIC EXPRESSION>;
THE SERIAL NUMBER OF THE FIRST AND LAST EIGENVALUE TO BE
CALCULATED, RESPECTIVELY;
VAL: <ARRAY IDENTIFIER>;
"ARRAY" VAL[N1:N2];
EXIT: THE N2-N1+1 CALCULATED CONSECUTIVE EIGENVALUES IN
NONINCREASING ORDER;
EM: <ARRAY IDENTIFIER>;
"ARRAY" EM[0:3];
ENTRY: EM[0], THE MACHINE PRECISION,
EM[1], AN UPPERBOUND FOR THE MODULI OF THE
EIGENVALUES OF THE GIVEN MATRIX,
EM[2], A RELATIVE TOLERANCE FOR THE EIGENVALUES;
EXIT: EM[3], THE TOTAL NUMBER OF ITERATIONS USED FOR
CALCULATING THE EIGENVALUES.
PROCEDURES USED:
ZEROIN = CP34150.
RUNNING TIME:
DEPENDS STRONGLY ON THE DISTANCE OF SUCCESSIVE EIGENVALUES.
LANGUAGE: ALGOL 60.
METHOD AND PERFORMANCE:
LET T DENOTE THE GIVEN SYMMETRIC TRIDIAGONAL MATRIX OF ORDER N AND
I THE IDENTITY MATRIX. THE EIGENVALUES OF T ARE THE ZEROES OF THE
N-TH DEGREE POLYNOMIAL P(N,X) = DET(T - X*I). INSTEAD OF SEARCHING
FOR THE ZEROES OF P(N,X) WE LOOK FOR THE ZEROES OF THE FUNCTION
F(N,X) = P(N,X) / P(N-1,X). MAINTAINING A LOWER BOUND FOR
ABS(P(N-1,X)) WE DO AVOID OVERFLOW OF THE REAL NUMBER CAPACITY IN
THE COMPUTATION OF F(N,X). THIS FUNCTION CAN BE CALCULATED AS
FOLLOWS:
F(1,X) = D[1] - X,
F(K,X) = D[K] - X - BB[K-1] /
("IF" ABS(F(K-1,X)) > MACHTOL "THEN" F(K-1,X)
"ELSE" "IF" F(K-1,X) <= 0 "THEN" -MACHTOL
"ELSE" MACHTOL), K = 2, . . . ,N,
WHERE MACHTOL EQUALS EM[0] * EM[1].
USING THE STURM SEQUENCE PROPERTY OF (F(K,X)), K=1,2,...,N, WE CAN
LOCATE THE DESIRED EIGENVALUES BY MEANS OF THE PROCEDURE ZEROIN
(SECTION 5.1.1.1). FOR FURTHER DETAILS SEE REF[1], REF[2].
SUBSECTION: VECSYMTRI.
CALLING SEQUENCE:
THE HEADING OF THE PROCEDURE IS:
"PROCEDURE" VECSYMTRI(D, B, N, N1, N2, VAL, VEC, EM);
"VALUE" N, N1, N2; "INTEGER" N, N1, N2;
"ARRAY" D, B, VAL, VEC, EM;
"CODE" 34152;
THE MEANING OF THE FORMAL PARAMETERS IS:
N: <ARITHMETIC EXPRESSION>;
THE ORDER OF THE GIVEN MATRIX;
D: <ARRAY IDENTIFIER>;
"ARRAY" D[1:N],
ENTRY: THE MAIN DIAGONAL OF THE SYMMETRIC TRIDIAGONAL
MATRIX;
B: <ARRAY IDENTIFIER>;
"ARRAY" B[1:N];
ENTRY: THE CODIAGONAL OF THE SYMMETRIC TRIDIAGONAL MATRIX
FOLLOWED BY AN ADDITIONAL ELEMENT 0;
N1, N2: <ARITHMETIC EXPRESSION>;
LOWER AND UPPER BOUND OF THE ARRAY VAL (SEE ALSO METHOD AND
PERFORMANCE);
VAL: <ARRAY IDENTIFIER>;
"ARRAY" VAL[N1:N2];
ENTRY: A ROW OF NONINCREASING EIGENVALUES AS DELIVERED BY
VALSYMTRI;
VEC: <ARRAY IDENTIFIER>;
"ARRAY" VEC[1:N,N1:N2];
EXIT: THE EIGENVECTORS CORRESPONDING WITH THE GIVEN
EIGENVALUES (SEE ALSO METHOD AND PERFORMANCE);
EM: <ARRAY IDENTIFIER>;
"ARRAY" EM[0:9];
ENTRY: EM[0], THE MACHINE PRECISION,
EM[1], A NORM OF THE GIVEN MATRIX,
EM[4], THE ORTHOGONALISATION PARAMETER (SEE ALSO
METHOD AND PERFORMANCE),
EM[6], THE RELATIVE TOLERANCE FOR THE EIGENVECTORS,
EM[8], THE MAXIMUM NUMBER OF ITERATIONS ALLOWED
FOR THE CALCULATION OF EACH EIGENVECTOR;
EXIT: EM[5], THE NUMBER OF EIGENVECTORS INVOLVED IN THE
LAST GRAM-SCHMIDT ORTHOGONALISATION (SEE
METHOD AND PERFORMANCE),
EM[7], THE MAXIMUM EUCLIDEAN NORM OF THE RESIDUES,
EM[9], THE LARGEST NUMBER OF ITERATIONS PERFORMED
FOR THE CALCULATION OF SOME EIGENVECTOR (SEE
METHOD AND PERFORMANCE).
PROCEDURES USED:
VECVEC = CP34010,
TAMVEC = CP34012,
ELMVECCOL= CP34021.
REQUIRED CENTRAL MEMORY:
EXECUTION FIELD LENGTH: FIVE AUXILIARY ONE-DIMENSIONAL REAL ARRAYS
AND ONE BOOLEAN ARRAY, ALL OF LENGTH N, ARE USED.
RUNNING TIME: THE PROCESS IS OF ORDER N FOR EACH EIGENVECTOR.
LANGUAGE: ALGOL 60.
METHOD AND PERFORMANCE:
AN EIGENVECTOR OF A SYMMETRIC TRIDIAGONAL MATRIX T, CORRESPONDING
TO AN EIGENVALUE LAMBDA, IS CALCULATED BY MEANS OF INVERSE
ITERATION; I.E. STARTING FROM SOME INITIAL VECTOR X, THE LINEAR
SYSTEM (T - LAMBDA * I)Y = X IS SOLVED ITERATIVELY, THE SOLUTION Y,
DIVIDED BY ITS EUCLIDEAN NORM REPLACING X EACH TIME.
IF THE DISTANCE BETWEEN SOME APPROXIMATE EIGENVALUES IS SMALLER
THAN MACHTOL (=EM[0] * EM[1]), THEN THEY ARE SLIGHTLY MODIFIED SUCH
THAT THE DISTANCE BETWEEN THEM EQUALS MACHTOL. IF THE DISTANCE
BETWEEN SOME EIGENVALUES IS SMALLER THAN THE ORTHOGONALISATION
PARAMETER (=EM[4]) TIMES EM[1], THEN IN EACH ITERATION STEP GRAM-
SCHMIDT ORTHOGONALISATION IS CARRIED OUT, SO THAT THE EIGENVECTORS
OBTAINED ARE ORTHOGONAL WITHIN WORKING PRECISION. THE ITERATION
ENDS AS SOON AS EITHER THE EUCLIDEAN NORM OF THE RESIDUE IS SMALLER
THAN EM[1] * EM[6], OR THE MAXIMUM ALLOWED NUMBER OF ITERATIONS
(=EM[8]) HAS BEEN PERFORMED. IN THE LATTER CASE EM[9]:= EM[8] + 1.
IF N1 > 1, THEN VECSYMTRI SHOULD BE PRECEDED BY ONE OR MORE CALLS
OF VECSYMTRI PRODUCING A NUMBER OF EIGENVECTORS CORRESPONDING TO
THE PRECEDING EIGENVALUES. MOREOVER ONE MUST GIVE EM[5], AS
PRODUCED BY THE LAST CALL OF VECSYMTRI; THE K-TH TO N2-TH
EIGENVALUES, WHERE K = N1 - EM[5], MUST BE GIVEN IN ARRAY VAL[K:N2]
IN MONOTONICALLY NONINCREASING ORDER (THE K-TH TO (N1-1)-TH
EIGENVALUES BEING NEEDED FOR THE MODIFYING MENTIONED ABOVE), AND
THE CORRESPONDING EIGENVECTORS UP TO THE (N1-1)-TH (WHICH ARE
NEEDED FOR THE GRAM-SCHMIDT ORTHOGONALISATION) IN THE CORRESPONDING
COLUMNS OF ARRAY VEC[1:N,K:N2].
THE TOLERANCES SHOULD SATISFY: EM[0] ( < EM[2]) < EM[6] AND
EM[4] >= EM[0] / EM[6]. FOR FURTHER DETAILS SEE REF[1].
SUBSECTION: QRIVALSYMTRI.
CALLING SEQUENCE:
THE HEADING OF THE PROCEDURE IS:
"INTEGER" "PROCEDURE" QRIVALSYMTRI(D, BB, N, EM);
"VALUE" N; "INTEGER" N; "ARRAY" D, BB, EM;
"CODE" 34160;
THE MEANING OF THE FORMAL PARAMETERS IS:
N: <ARITHMETIC EXPRESSION>;
THE ORDER OF THE GIVEN MATRIX;
D: <ARRAY IDENTIFIER>;
"ARRAY" D[1:N];
ENTRY: THE MAIN DIAGONAL OF THE SYMMETRIC TRIDIAGONAL
MATRIX;
EXIT: THE EIGENVALUES OF THE MATRIX IN SOME ARBITRARY
ORDER;
BB: <ARRAY IDENTIFIER>;
"ARRAY" BB[1:N];
ENTRY: THE SQUARES OF THE CODIAGONAL ELEMENTS OF THE
SYMMETRIC TRIDIAGONAL MATRIX FOLLOWED BY AN
ADDITONAL ELEMENT O;
EXIT: THE SQUARES OF THE CODIAGONAL ELEMENTS OF THE
SYMMETRIC TRIDIAGONAL MATRIX RESULTING FROM THE QR
ITERATION;
EM: <ARRAY IDENTIFIER>;
"ARRAY" EM[0:5];
ENTRY: EM[0], THE MACHINE PRECISION;
EM[1], A NORM OF THE GIVEN MATRIX;
EM[2], A RELATIVE TOLERANCE FOR THE EIGENVALUES;
EM[4], THE MAXIMUM ALLOWED NUMBER OF ITERATIONS;
EXIT: EM[3], THE MAXIMUM ABSOLUTE VALUE OF THE CODIAGONAL
ELEMENTS NEGLECTED;
EM[5], THE NUMBER OF ITERATIONS PERFORMED.
MOREOVER:
QRIVALSYMTRI:= THE NUMBER OF EIGENVALUES NOT CALCULATED.
PROCEDURES USED: NONE.
RUNNING TIME: THE PROCESS IS OF ORDER N SQUARED.
LANGUAGE: ALGOL 60.
METHOD AND PERFORMANCE:
IN QRIVALSYMTRI THE EIGENVALUES OF A SYMMETRIC TRIDIAGONAL MATRIX
ARE CALCULATED BY MEANS OF QR-ITERATION. FOR THIS PROCEDURE WE USED
ESSENTIALLY THE SQUARE-ROOT-FREE VERSION OF THE QR ALGORITHM DUE TO
REINSCH[3].
IN ADDITION TO THE RELATIVE ERROR, WHICH IS SUPPOSED TO BE BOUNDED
BY EM[1] * EM[2] (I.E. MATRIX NORM TIMES RELATIVE TOLERANCE), THE
CALCULATED EIGENVALUES HAVE AN ABSOLUTE ERROR WHICH IS BOUNDED BY
BY EM[0] * EM[1] (I.E. MACHINE PRECISION TIMES MATRIX NORM).
IN PARTICULAR, WHEN SOME EIGENVALUES ARE VERY SMALL COMPARED TO THE
MATRIX NORM, THE ACCURACY OF THE CALCULATED EIGENVALUES CAN BE
INCREASED BY GIVING EM[0] A (POSITIVE) VALUE WHICH IS LESS THAN
THE MACHINE PRECISION.
A PARTICULAR CHOICE OF EM[0] IS HARMLESS FOR THE PROCEDURE
PROVIDED THAT FOR EACH I THE CALCULATION OF BB[I] / EM[0] ** 2
CAUSES NO OVERFLOW AND THE CALCULATION OF (EM[0] * EM[1]) ** 2
CAUSES NO UNDERFLOW.
ONE SHOULD NOTICE THAT THE NUMBER OF QR ITERATIONS INCREASES BY A
SMALLER CHOICE OF EM[0].
FOR FURTHER DETAILS SEE [2], [3].
SUBSECTION: QRISYMTRI.
CALLING SEQUENCE:
THE HEADING OF THE PROCEDURE IS:
"INTEGER" "PROCEDURE" QRISYMTRI(A, N, D, B, BB, EM);
"VALUE" N; "INTEGER" N; "ARRAY" A, D, B, BB, EM;
"CODE" 34161;
THE MEANING OF THE FORMAL PARAMETERS IS:
N: <ARITHMETIC EXPRESSION>;
THE ORDER OF THE GIVEN MATRIX;
D: <ARRAY IDENTIFIER>;
"ARRAY" D[1:N];
ENTRY: THE MAIN DIAGONAL OF THE SYMMETRIC TRIDIAGONAL
MATRIX;
EXIT: THE EIGENVALUES OF THE MATRIX IN SOME ARBITRARY
ORDER;
B: <ARRAY IDENTIFIER>;
"ARRAY" B[1:N];
ENTRY: THE CODIAGONAL OF THE SYMMETRIC TRIDIAGONAL MATRIX
FOLLOWED BY AN ADDITIONAL ELEMENT 0;
EXIT: THE CODIAGONAL OF THE SYMMETRIC TRIDIAGONAL MATRIX
RESULTING FROM THE QR ITERATION, FOLLOWED BY AN
ADDITIONAL ELEMENT 0;
BB: <ARRAY IDENTIFIER>;
"ARRAY" BB[1:N];
ENTRY: THE SQUARED CODIAGONAL ELEMENTS OF THE SYMMETRIC
TRIDIAGONAL MATRIX, FOLLOWED BY AN ADDITIONAL
ELEMENT O;
EXIT: THE SQUARED CODIAGONAL ELEMENTS OF THE SYMMETRIC
TRIDIAGONAL MATRIX RESULTING FROM THE QR ITERATION;
A: <ARRAY IDENTIFIER>;
"ARRAY" A[1:N,1:N];
ENTRY: SOME MATRIX S, SAY, (POSSIBLY THE IDENTITY MATRIX);
EXIT: THE EIGENVECTORS OF THE ORIGINAL SYMMETRIC
TRIDIAGONAL MATRIX, PREMULTIPLIED BY S (SEE METHOD
AND PERFORMANCE);
EM: <ARRAY IDENTIFIER>;
"ARRAY" EM[0:5];
ENTRY: EM[0], THE MACHINE PRECISION;
EM[1], A NORM OF THE GIVEN MATRIX;
EM[2], A RELATIVE TOLERANCE FOR THE QR ITERATION;
EM[4], THE MAXIMUM ALLOWED NUMBER OF ITERATIONS;
EXIT: EM[3], THE MAXIMUM ABSOLUTE VALUE OF THE CODIAGONAL
ELEMENTS NEGLECTED;
EM[5], THE NUMBER OF ITERATIONS PERFORMED.
MOREOVER:
QRISYMTRI:= THE NUMBER OF EIGENVALUES AND -VECTORS NOT
CALCULATED.
PROCEDURES USED:
ROTCOL = CP34040.
RUNNING TIME: THE PROCESS IS OF ORDER N CUBED.
LANGUAGE: ALGOL 60.
METHOD AND PERFORMANCE:
IN QRISYMTRI THE EIGENVALUES AND EIGENVECTORS OF A SYMMETRIC
TRIDIAGONAL MATRIX ARE COMPUTED SIMULTANEOUSLY.
IN MOST APPLICATIONS QRISYMTRI IS USED IN THE COMPUTATION OF
EIGENVALUES AND -VECTORS OF A GENERAL SYMMETRIC MATRIX (SEE QRISYM
SECTION 3.3.1.1.2); IN THAT CASE ARRAY A IS INITIALLY GIVEN THE
VALUE OF THE TRANSFORMING MATRIX (TFMPREVEC SECTION 3.2.1.2.1.1).
FOR THE COMPUTATION OF EIGENVALUES AND EIGENVECTORS OF A SYMMETRIC
TRIDIAGONAL MATRIX, ARRAY A HAS TO BE INITIALIZED TO THE IDENTITY
MATRIX. THE AVERAGE NUMBER OF ITERATIONS IS ABOUT 3N. WHEN THE
PROCESS IS COMPLETED WITHIN EM[4] ITERATIONS, THEN QRISYMTRI:= 0;
OTHERWISE QRISYMTRI:= THE NUMBER, K, OF EIGENVALUES NOT CALCULATED,
EM[5]:= EM[4] + 1 AND ONLY THE LAST N - K ELEMENTS OF D AND THE
LAST N - K COLUMNS OF A ARE APPROXIMATE EIGENVALUES AND -VECTORS
RESPECTIVELY, OF THE ORIGINAL MATRIX.
FOR FURTHER DETAILS SEE REF[1], REF[2].
REFERENCES:
[1] DEKKER, T.J. AND HOFFMANN, W.
ALGOL 60 PROCEDURES IN NUMERICAL ALGEBRA, PART 2,
MATHEMATICAL CENTRE TRACTS 23,
MATHEMATISCH CENTRUM, AMSTERDAM, 1968;
[2] WILKINSON, J.H.
THE ALGEBRAIC EIGENVALUE PROBLEM,
CLARENDON PRESS, OXFORD 1965.
[3] REINSCH, CHR.H.
A STABLE, RATIONAL QR ALGORITHM FOR THE COMPUTATION OF THE
EIGENVALUES OF AN HERMITIAN, TRIDIAGONAL MATRIX.
MATH. OF COMP. VOL 25(1971) PP. 591-597.
EXAMPLE OF USE:
THE FIRST AND SECOND EIGENVALUE IN MONOTONICALLY NON-INCREASING
ORDER AND THE CORRESPONDING EIGENVECTORS OF T, WITH N = 4 AND
T[I,J] = "IF" I = J "THEN" 2 "ELSE" "IF" ABS(I - J) = 1 "THEN" - 1
"ELSE" 0, MAY BE OBTAINED BY THE FOLLOWING PROGRAM:
"BEGIN"
"INTEGER" J;
"ARRAY" B, D[1:4], BB[1:3], VAL[1:2], EM[0:9], VEC[1:4,1:2];
EM[0]:= "-14; EM[1]:= 4; EM[2]:= "-12;
EM[4]:= "-3; EM[6]:= "-10; EM[8]:= 5;
"FOR" J:= 1, 2, 3, 4 "DO" D[J]:= 2; B[4]:= 0;
"FOR" J:= 1, 2, 3 "DO"
"BEGIN" BB[J]:= 1; B[J]:= -1 "END";
VALSYMTRI(D, BB, 4, 1, 2, VAL, EM);
VECSYMTRI(D, B, 4, 1, 2, VAL, VEC, EM);
OUTPUT(61, "("2(+.13D"+2D, 2B), 2/")", VAL[1], VAL[2]);
"FOR" J:= 1, 2, 3, 4 "DO"
OUTPUT(61, "("2(+.13D"+2D, 2B), /")", VEC[J,1], VEC[J,2]);
OUTPUT(61, "("/, .2D"+2D, /, 3(2ZD, /)")",
EM[7], EM[3], EM[5], EM[9])
"END"
THE PROGRAM DELIVERS:
THE EIGENVALUES: +.3618033988751"+01 +.2618033988750"+01
THE EIGENVECTORS: +.3717480344602"+00 +.6015009550075"+00
+.6015009550075"+00 +.3717480344602"+00
+.6015009550075"+00 -.3717480344602"+00
+.3717480344602"+00 -.6015009550075"+00
EM[7] = .15"-11
EM[3] = 24
EM[5] = 1
EM[9] = 1 .
SOURCE TEXT(S):
"CODE" 34151;
"CODE" 34152;
"CODE" 34160;
"CODE" 34161;