NUMAL Section 3.3.1.2.2

BEGIN SECTION : 3.3.1.2.2 (July, 1974)

AUTHORS   : T.J. DEKKER, W. HOFFMANN.

CONTRIBUTORS: W. HOFFMANN, S.P.N. VAN KAMPEN, J.G. VERWER.

INSTITUTE: MATHEMATICAL CENTRE.

RECEIVED: 731205.

BRIEF DESCRIPTION:
    THIS SECTION CONTAINS FIVE PROCEDURES FOR  CALCULATING  EIGENVALUES
    AND / OR EIGENVECTORS OF REAL MATRICES:
    A) REAEIGVAL CALCULATES THE EIGENVALUES OF A MATRIX, PROVIDED  THAT
    ALL EIGENVALUES ARE REAL,
    B) REAEIG1 CALCULATES THE EIGENVALUES, PROVIDED THAT  THEY  ARE ALL
    REAL, AND THE EIGENVECTORS OF A MATRIX,
    C) REAEIG3 CALCULATES THE EIGENVALUES, PROVIDED THAT  THEY  ARE ALL
    REAL, AND THE EIGENVECTORS OF A MATRIX,
    D) COMEIGVAL CALCULATES THE EIGENVALUES OF A MATRIX,
    E) COMEIG1 CALCULATES THE EIGENVALUES AND EIGENVECTORS OF A MATRIX.

KEYWORDS:
    EIGENVALUES,
    EIGENVECTORS.


SUBSECTION: REAEIGVAL.

CALLING SEQUENCE:
    THE HEADING OF THE PROCEDURE IS:
    "INTEGER" "PROCEDURE" REAEIGVAL(A, N, EM, VAL); "VALUE" N;
    "INTEGER" N; "ARRAY" A, EM, VAL;
    "CODE" 34182;

    THE MEANING OF THE FORMAL PARAMETERS IS:
    A:      <ARRAY IDENTIFIER>;
            "ARRAY" A[1:N,1:N];
            ENTRY: THE MATRIX WHOSE EIGENVALUES ARE TO BE CALCULATED;
            EXIT:  THE ARRAY ELEMENTS ARE ALTERED;
    N:      <ARITHMETIC EXPRESSION>;
            THE ORDER OF THE GIVEN MATRIX;
    EM:     <ARRAY IDENTIFIER>;
            "ARRAY" EM[0:5];
            ENTRY: EM[0], THE MACHINE PRECISION;
                   EM[2], THE  RELATIVE TOLERANCE USED  FOR THE  QR
                          ITERATION;
                   EM[4], THE  MAXIMUM  ALLOWED  NUMBER  OF ITERATIONS;
                   FOR THE  CD CYBER 73-28  SUITABLE VALUES  OF THE
                   DATA TO BE GIVEN IN EM ARE:
                   EM[0] = "-14,
                   EM[2] > EM[0] (E.G. EM[2] = "-13),
                   EM[4] = 10 * N;
            EXIT:  EM[1], THE INFINITY NORM OF THE EQUILIBRATED MATRIX;
                   EM[3], THE MAXIMUM ABSOLUTE VALUE OF THE SUBDIAGONAL
                          ELEMENTS NEGLECTED;
                   EM[5], THE  NUMBER OF QR  ITERATIONS  PERFORMED;
                          IF THE  ITERATION  PROCESS IS  NOT  COMPLETED
                          WITHIN EM[4] ITERATIONS, THE VALUE  EM[4] + 1
                          IS DELIVERED  AND IN THIS CASE  ONLY THE LAST
                          N - K ELEMENTS OF VAL  ARE APPROXIMATE EIGEN-
                          VALUES  OF  THE  GIVEN  MATRIX,  WHERE  K  IS
                          DELIVERED IN REAEIGVAL;
    VAL:    <ARRAY IDENTIFIER>;
            "ARRAY" VAL[1:N];
            EXIT: THE  EIGENVALUES OF THE GIVEN  MATRIX  ARE  DELIVERED
                  IN MONOTONICALLY NONINCREASING ORDER;

    MOREOVER:
    REAEIGVAL DELIVERS 0, PROVIDED THAT THE PROCESS IS COMPLETED WITHIN
    EM[4] ITERATIONS; OTHERWISE  REAEIGVAL  DELIVERS  K, THE  NUMBER OF
    EIGENVALUES NOT CALCULATED.

PROCEDURES USED:
    EQILBR    = CP34173,
    TFMREAHES = CP34170,
    REAVALQRI = CP34180.

REQUIRED CENTRAL MEMORY:
    EXECUTION FIELD LENGTH: 3N.

RUNNING TIME: ROUGHLY PROPORTIONAL TO N CUBED.

LANGUAGE: ALGOL 60.

METHOD AND PERFORMANCE:
    THE GIVEN MATRIX IS EQUILIBRATED  BY  CALLING  EQILBR (SEE  SECTION
    3.2.1.1.1) AND TRANSFORMED TO A SIMILAR UPPER-HESSENBERG   MATRIX
    BY CALLING TFMREAHES  (SEE  SECTION 3.2.1.2.1.2). THE   EIGENVALUES
    ARE THEN CALCULATED BY CALLING REAVALQRI, WHICH USES SINGLE  QR
    ITERATION (SEE SECTION 3.3.1.2.1).
    THE PROCEDURE REAEIGVAL SHOULD  BE USED ONLY IF ALL EIGENVALUES ARE
    REAL.
    FOR FURTHER DETAILS SEE REFERENCES [1], [2] AND [3].


SUBSECTION: REAEIG1.

CALLING SEQUENCE:
    THE HEADING OF THE PROCEDURE IS:
    "INTEGER" "PROCEDURE" REAEIG1(A, N, EM, VAL, VEC); "VALUE" N;
    "INTEGER" N; "ARRAY" A, EM, VAL, VEC;
    "CODE" 34184;

    THE MEANING OF THE FORMAL PARAMETERS IS:
    A:      <ARRAY IDENTIFIER>;
            "ARRAY" A[1:N,1:N];
            ENTRY: THE MATRIX WHOSE EIGENVALUES AND EIGENVECTORS ARE TO
                   BE CALCULATED;
            EXIT:  THE ARRAY ELEMENTS ARE ALTERED;
    N:      <ARITHMETIC EXPRESSION>;
            THE ORDER OF THE GIVEN MATRIX;

    EM:     <ARRAY IDENTIFIER>;
            "ARRAY" EM[0:9];
            ENTRY: EM[0], THE MACHINE PRECISION;
                   EM[2], THE RELATIVE  TOLERANCE  USED  FOR THE QR
                          ITERATION;
                   EM[4], THE   MAXIMUM   ALLOWED   NUMBER   OF  QR
                          ITERATIONS;
                   EM[6], THE TOLERANCE  USED  FOR  THE   EIGENVECTORS;
                          FOR EACH EIGENVECTOR  THE  INVERSE  ITERATION
                          ENDS IF THE EUCLIDEAN  NORM  OF THE RESIDUE
                          VECTOR IS SMALLER THAN EM[1] * EM[6];
                   EM[8], THE   MAXIMUM   ALLOWED   NUMBER  OF  INVERSE
                          ITERATIONS  FOR  THE  CALCULATION  OF    EACH
                          EIGENVECTOR;
                   FOR THE  CD CYBER 73-28  SUITABLE VALUES  OF THE
                   DATA TO BE GIVEN IN EM ARE:
                   EM[0] = "-14,
                   EM[2] > EM[0] (E.G. EM[2] = "-13),
                   EM[4] = 10 * N,
                   EM[6] > EM[2] (E.G. EM[6] = "-10),
                   EM[8] = 5;
            EXIT:  EM[1], THE INFINITY NORM OF THE EQUILIBRATED MATRIX;
                   EM[3], THE MAXIMUM ABSOLUTE VALUE OF THE SUBDIAGONAL
                          ELEMENTS NEGLECTED;
                   EM[5], THE  NUMBER OF QR  ITERATIONS  PERFORMED;
                          IF THE  ITERATION  PROCESS IS  NOT  COMPLETED
                          WITHIN EM[4] ITERATIONS, THE VALUE  EM[4] + 1
                          IS DELIVERED  AND IN THIS CASE  ONLY THE LAST
                          N - K ELEMENTS OF VAL AND COLUMNS OF VEC  ARE
                          APPROXIMATE EIGENVALUES AND  EIGENVECTORS  OF
                          THE GIVEN MATRIX, WHERE  K IS  DELIVERED   IN
                          REAEIG1;
                   EM[7], THE  MAXIMUM EUCLIDIAN NORM OF THE RESIDUES
                          OF THE CALCULATED EIGENVECTORS (OF THE TRANS-
                          FORMED MATRIX);
                   EM[9], THE  LARGEST  NUMBER  OF  INVERSE  ITERATIONS
                          PERFORMED FOR THE CALCULATION OF SOME  EIGEN-
                          VECTOR;  IF,  FOR   SOME   EIGENVECTOR    THE
                          EUCLIDEAN  NORM  OF  THE  RESIDUE   REMAINS
                          LARGER   THAN    EM[1] * EM[6],  THE    VALUE
                          EM[8] + 1  IS  DELIVERED;  NEVERTHELESS   THE
                          EIGENVECTORS  MAY  THEN  VERY WELL BE USEFUL,
                          THIS  SHOULD  BE   JUDGED   FROM   THE  VALUE
                          DELIVERED IN EM[7] OR FROM SOME OTHER TEST;
    VAL:    <ARRAY IDENTIFIER>;
            "ARRAY" VAL[1:N];
            EXIT: THE  EIGENVALUES OF THE GIVEN  MATRIX  ARE  DELIVERED
                  IN MONOTONICALLY DECREASING ORDER;
    VEC:    <ARRAY IDENTIFIER>;
            "ARRAY" VEC[1:N,1:N];
            EXIT: THE CALCULATED EIGENVECTORS, CORRESPONDING   TO   THE
                  EIGENVALUES IN ARRAY VAL[1:N], ARE DELIVERED  IN  THE
                  COLUMNS OF ARRAY VEC;

    MOREOVER:
    REAEIG1   DELIVERS 0, PROVIDED THAT THE PROCESS IS COMPLETED WITHIN
    EM[4] ITERATIONS; OTHERWISE  REAEIG1    DELIVERS  K, THE  NUMBER OF
    EIGENVALUES AND EIGENVECTORS NOT CALCULATED.

PROCEDURES USED:
    EQILBR     = CP34173,
    TFMREAHES  = CP34170,
    BAKREAHES2 = CP34172,
    BAKLBR     = CP34174,
    REAVALQRI  = CP34180,
    REAVECHES  = CP34181,
    REASCL     = CP34183.

REQUIRED CENTRAL MEMORY:
    EXECUTION FIELD LENGTH: N * N + 5N.

RUNNING TIME: ROUGHLY PROPORTIONAL TO N CUBED.

OPTIONS: F.

METHOD AND PERFORMANCE:
    THE GIVEN MATRIX IS EQUILIBRATED  BY  CALLING  EQILBR (SEE  SECTION
    3.2.1.1.1) AND TRANSFORMED TO A SIMILAR UPPER-HESSENBERG   MATRIX
    BY CALLING TFMREAHES  (SEE  SECTION 3.2.1.2.1.2). THE   EIGENVALUES
    ARE THEN CALCULATED BY CALLING REAVALQRI, WHICH USES SINGLE  QR
    ITERATION (SEE SECTION 3.3.1.2.1).
    FURTHERMORE, TO FIND THE EIGENVECTORS WILKINSON'S DEVICE IS FIRST
    APPLIED [2, P.328 AND 628]. SUBSEQUENTLY  THE  EIGENVECTORS  OF THE
    UPPER-HESSENBERG MATRIX  ARE  CALCULATED  BY  CALLING  REAVECHES,
    WHICH   USES   INVERSE   ITERATION  (SEE  SECTION  3.3.1.2.1).  THE
    CALCULATED VECTORS ARE THEN BACK-TRANSFORMED TO  THE  CORRESPONDING
    EIGENVECTORS OF THE GIVEN MATRIX BY CALLING BAKREAHES2  AND  BAKLBR
    (SEE SECTIONS 3.2.1.2.1.2 AND 3.2.1.1.1). FINALLY  THE  APPROXIMATE
    EIGENVECTORS ARE NORMALIZED BY  CALLING  REASCL (SEE SECTION 1.1.9)
    SUCH THAT, IN EACH  EIGENVECTOR, AN  ELEMENT  OF  MAXIMUM  ABSOLUTE
    VALUE EQUALS 1.
    THE PROCEDURE REAEIG1 SHOULD BE USED  ONLY IF ALL  EIGENVALUES  ARE
    REAL.
    FOR FURTHER DETAILS SEE THE GIVEN REFERENCES.


SUBSECTION: REAEIG3.

CALLING SEQUENCE:
    THE HEADING OF THE PROCEDURE IS:
    "INTEGER" "PROCEDURE" REAEIG3(A, N, EM, VAL, VEC); "VALUE" N;
    "INTEGER" N; "ARRAY" A, EM, VAL, VEC;
    "CODE" 34187;

    THE MEANING OF THE FORMAL PARAMETERS IS:
    A:      <ARRAY IDENTIFIER>;
            "ARRAY" A[1:N,1:N];
            ENTRY: THE MATRIX WHOSE EIGENVALUES AND EIGENVECTORS ARE TO
                   BE CALCULATED;
            EXIT:  THE ARRAY ELEMENTS ARE ALTERED;
    N:      <ARITHMETIC EXPRESSION>;
            THE ORDER OF THE GIVEN MATRIX;
    EM:     <ARRAY IDENTIFIER>;
            "ARRAY" EM[0:5];
            ENTRY: EM[0], THE MACHINE PRECISION;
                   EM[2], THE RELATIVE  TOLERANCE  USED  FOR THE QR
                          ITERATION;
                   EM[4], THE   MAXIMUM   ALLOWED   NUMBER   OF  QR
                          ITERATIONS;
                   FOR THE  CD CYBER 73-28  SUITABLE VALUES  OF THE
                   DATA TO BE GIVEN IN EM ARE:
                   EM[0] = "-14,
                   EM[2] > EM[0] (E.G. EM[2] = "-13),
                   EM[4] = 10 * N;
            EXIT:  EM[1], THE INFINITY NORM OF THE EQUILIBRATED MATRIX;
                   EM[3], THE MAXIMUM ABSOLUTE VALUE OF THE SUBDIAGONAL
                          ELEMENTS NEGLECTED;
                   EM[5], THE  NUMBER OF QR  ITERATIONS  PERFORMED;
                          IF THE  ITERATION  PROCESS IS  NOT  COMPLETED
                          WITHIN EM[4] ITERATIONS, THE VALUE  EM[4] + 1
                          IS DELIVERED. IN   THIS  CASE  ONLY THE  LAST
                          N - K  ELEMENTS  OF  VAL   ARE    APPROXIMATE
                          EIGENVALUES OF THE GIVEN MATRIX AND NO USEFUL
                          EIGENVECTORS ARE DELIVERED. THE  VALUE  K  IS
                          DELIVERED IN REAEIG3;
    VAL:    <ARRAY IDENTIFIER>;
            "ARRAY" VAL[1:N];
            EXIT: THE  EIGENVALUES OF THE GIVEN  MATRIX ARE  DELIVERED;

    VEC:    <ARRAY IDENTIFIER>;
            "ARRAY" VEC[1:N,1:N];
            EXIT: THE CALCULATED EIGENVECTORS, CORRESPONDING   TO   THE
                  EIGENVALUES IN ARRAY VAL[1:N], ARE DELIVERED  IN  THE
                  COLUMNS OF ARRAY VEC;

    MOREOVER:
    REAEIG3   DELIVERS 0, PROVIDED THAT THE PROCESS IS COMPLETED WITHIN
    EM[4] ITERATIONS; OTHERWISE  REAEIG3    DELIVERS  K, THE  NUMBER OF
    EIGENVALUES NOT CALCULATED.

PROCEDURES USED:
    EQILBR     = CP34173,
    TFMREAHES  = CP34170,
    BAKREAHES2 = CP34172,
    BAKLBR     = CP34174,
    REAQRI     = CP34186,
    REASCL     = CP34183.

REQUIRED CENTRAL MEMORY:
    EXECUTION FIELD LENGTH: 4N.

RUNNING TIME: ROUGHLY PROPORTIONAL TO N CUBED.

LANGUAGE: ALGOL 60.

METHOD AND PERFORMANCE:
    THE GIVEN MATRIX IS EQUILIBRATED  BY  CALLING  EQILBR (SEE  SECTION
    3.2.1.1.1) AND TRANSFORMED TO A SIMILAR UPPER-HESSENBERG   MATRIX
    BY CALLING TFMREAHES  (SEE  SECTION 3.2.1.2.1.2). THE   EIGENVALUES
    AND  EIGENVECTORS  OF  THE   UPPER-HESSENBERG   MATRIX  ARE  THEN
    CALCULATED BY CALLING REAQRI, WHICH USES SINGLE  QR   ITERATION
    FOR THE EIGENVALUES AND A DIRECT METHOD FOR THE  EIGENVECTORS  (SEE
    SECTION 3.3.1.2.1). FINALLY  THE  EIGENVECTORS   OF  THE   UPPER-
    HESSENBERG MATRIX ARE BACK-TRANSFORMED TO THE CORRESPONDING  EIGEN-
    VECTORS OF THE  GIVEN  MATRIX  BY  CALLING  BAKREAHES2 (SEE SECTION
    3.1.2.1.2.1) AND NORMALIZED BY CALLING  REASCL  (SEE SECTION 1.1.9)
    SUCH THAT, IN EACH EIGENVECTOR, AN  ELEMENT  OF  MAXIMUM   ABSOLUTE
    VALUE EQUALS 1.
    THE PROCEDURE REAEIG3  SHOULD  BE USED  ONLY IF ALL EIGENVALUES ARE
    REAL.
    FOR FURTHER DETAILS SEE THE GIVEN REFERENCES.


SUBSECTION: COMEIGVAL.

CALLING SEQUENCE:
    THE HEADING OF THE PROCEDURE IS:
    "INTEGER" "PROCEDURE" COMEIGVAL(A, N, EM, RE, IM); "VALUE" N;
    "INTEGER" N; "ARRAY" A, EM, RE, IM;
    "CODE" 34192;

    THE MEANING OF THE FORMAL PARAMETERS IS:
    A:      <ARRAY IDENTIFIER>;
            "ARRAY" A[1:N,1:N];
            ENTRY: THE MATRIX WHOSE EIGENVALUES ARE TO BE CALCULATED;
            EXIT:  THE ARRAY ELEMENTS ARE ALTERED;
    N:      <ARITHMETIC EXPRESSION>;
            THE ORDER OF THE GIVEN MATRIX;
    EM:     <ARRAY IDENTIFIER>;
            "ARRAY" EM[0:5];
            ENTRY: EM[0], THE MACHINE PRECISION;
                   EM[2], THE  RELATIVE TOLERANCE USED  FOR THE  QR
                          ITERATION;
                   EM[4], THE  MAXIMUM  ALLOWED  NUMBER  OF ITERATIONS;
                   FOR THE  CD CYBER 73-28  SUITABLE VALUES  OF THE
                   DATA TO BE GIVEN IN EM ARE:
                   EM[0] = "-14,
                   EM[2] > EM[0] (E.G. EM[2] = "-13),
                   EM[4] = 10 * N;
            EXIT:  EM[1], THE INFINITY NORM OF THE EQUILIBRATED MATRIX;
                   EM[3], THE MAXIMUM ABSOLUTE VALUE OF THE SUBDIAGONAL
                          ELEMENTS NEGLECTED;
                   EM[5], THE  NUMBER OF QR  ITERATIONS  PERFORMED;
                          IF THE  ITERATION  PROCESS IS  NOT  COMPLETED
                          WITHIN EM[4] ITERATIONS, THE VALUE  EM[4] + 1
                          IS DELIVERED  AND IN THIS CASE  ONLY THE LAST
                          N - K ELEMENTS  OF RE AND IM  ARE APPROXIMATE
                          EIGENVALUES  OF THE GIVEN MATRIX, WHERE  K IS
                          DELIVERED IN COMEIGVAL;
    RE,IM:  <ARRAY IDENTIFIER>;
            "ARRAY" RE, IM[1:N];
            EXIT:  THE  REAL  AND  IMAGINARY  PARTS  OF THE  CALCULATED
                   EIGENVALUES  OF THE  GIVEN  MATRIX  ARE DELIVERED IN
                   ARRAY  RE, IM[1:N], THE  MEMBERS  OF  EACH   NONREAL
                   COMPLEX CONJUGATE PAIR BEING CONSECUTIVE;

    MOREOVER:
    COMEIGVAL DELIVERS 0, PROVIDED THAT THE PROCESS IS COMPLETED WITHIN
    EM[4] ITERATIONS; OTHERWISE  COMEIGVAL  DELIVERS  K, THE  NUMBER OF
    EIGENVALUES NOT CALCULATED.

PROCEDURES USED:
    EQILBR    = CP34173,
    TFMREAHES = CP34170,
    COMVALQRI = CP34190.

REQUIRED CENTRAL MEMORY:
    EXECUTION FIELD LENGTH: 3N.

RUNNING TIME: ROUGHLY PROPORTIONAL TO N CUBED.

LANGUAGE: ALGOL 60.

METHOD AND PERFORMANCE:
    THE GIVEN MATRIX IS EQUILIBRATED  BY  CALLING  EQILBR (SEE  SECTION
    3.2.1.1.1) AND TRANSFORMED TO A SIMILAR UPPER-HESSENBERG   MATRIX
    BY CALLING TFMREAHES  (SEE  SECTION 3.2.1.2.1.2). THE   EIGENVALUES
    ARE THEN CALCULATED BY CALLING COMVALQRI, WHICH USES DOUBLE  QR
    ITERATION (SEE SECTION 3.3.1.2.1).
    FOR FURTHER DETAILS SEE REFERENCES [1], [2] AND [3].


SUBSECTION: COMEIG1.

CALLING SEQUENCE:
    THE HEADING OF THE PROCEDURE IS:
    "INTEGER" "PROCEDURE" COMEIG1(A, N, EM, RE, IM, VEC); "VALUE" N;
    "INTEGER" N; "ARRAY" A, EM, RE, IM, VEC;
    "CODE" 34194;

    THE MEANING OF THE FORMAL PARAMETERS IS:
    A:      <ARRAY IDENTIFIER>;
            "ARRAY" A[1:N,1:N];
            ENTRY: THE MATRIX WHOSE EIGENVALUES AND EIGENVECTORS ARE TO
                   BE CALCULATED;
            EXIT:  THE ARRAY ELEMENTS ARE ALTERED;
    N:      <ARITHMETIC EXPRESSION>;
            THE ORDER OF THE GIVEN MATRIX;
    EM:     <ARRAY IDENTIFIER>;
            "ARRAY" EM[0:9];
            ENTRY: EM[0], THE MACHINE PRECISION;
                   EM[2], THE RELATIVE  TOLERANCE  USED  FOR THE QR
                          ITERATION;
                   EM[4], THE   MAXIMUM   ALLOWED   NUMBER   OF  QR
                          ITERATIONS;
                   EM[6], THE TOLERANCE  USED  FOR  THE   EIGENVECTORS;
                          FOR EACH EIGENVECTOR  THE  INVERSE  ITERATION
                          ENDS IF THE EUCLIDEAN  NORM  OF THE RESIDUE
                          VECTOR IS SMALLER THAN EM[1] * EM[6];
                   EM[8], THE   MAXIMUM   ALLOWED   NUMBER  OF  INVERSE
                          ITERATIONS  FOR  THE  CALCULATION  OF    EACH
                          EIGENVECTOR;
                   FOR THE  CD CYBER 73-28  SUITABLE VALUES  OF THE
                   DATA TO BE GIVEN IN EM ARE:
                   EM[0] = "-14,
                   EM[2] > EM[0] (E.G. EM[2] = "-13),
                   EM[4] = 10 * N,
                   EM[6] > EM[2] (E.G. EM[6] = "-10),
                   EM[8] = 5;

            EXIT:  EM[1], THE INFINITY NORM OF THE EQUILIBRATED MATRIX;
                   EM[3], THE MAXIMUM ABSOLUTE VALUE OF THE SUBDIAGONAL
                          ELEMENTS NEGLECTED;
                   EM[5], THE  NUMBER OF QR  ITERATIONS  PERFORMED;
                          IF THE  ITERATION  PROCESS IS  NOT  COMPLETED
                          WITHIN EM[4] ITERATIONS, THE VALUE  EM[4] + 1
                          IS DELIVERED  AND IN THIS CASE  ONLY THE LAST
                          N - K ELEMENTS OF RE, IM  AND COLUMNS  OF VEC
                          ARE APPROXIMATE EIGENVALUES AND  EIGENVECTORS
                          OF THE GIVEN MATRIX, WHERE  K IS DELIVERED IN
                          COMEIG1;
                   EM[7], THE  MAXIMUM EUCLIDIAN NORM OF THE RESIDUES
                          OF THE CALCULATED EIGENVECTORS (OF THE TRANS-
                          FORMED MATRIX);
                   EM[9], THE  LARGEST  NUMBER  OF  INVERSE  ITERATIONS
                          PERFORMED FOR THE CALCULATION OF SOME  EIGEN-
                          VECTOR;  IF  THE   EUCLIDIAN  NORM  OF  THE
                          RESIDUE FOR ONE OR MORE EIGENVECTORS  REMAINS
                          LARGER THAN EM[1] * EM[6], THE VALUE  EM[8]+1
                          IS  DELIVERED; NEVERTHELESS  THE EIGENVECTORS
                          MAY THEN  VERY WELL BE USEFUL, THIS SHOULD BE
                          JUDGED  FROM THE VALUE  DELIVERED IN EM[7] OR
                          FROM SOME OTHER TEST;
    RE,IM:  <ARRAY IDENTIFIER>;
            "ARRAY" RE, IM[1:N];
            EXIT:  THE  REAL  AND  IMAGINARY  PARTS  OF THE  CALCULATED
                   EIGENVALUES  OF THE  GIVEN  MATRIX  ARE DELIVERED IN
                   ARRAY  RE, IM[1:N], THE  MEMBERS  OF  EACH   NONREAL
                   COMPLEX CONJUGATE PAIR BEING CONSECUTIVE;
    VEC:    <ARRAY IDENTFIER>;
            "ARRAY" VEC[1:N,1:N];
            EXIT:  THE CALCULATED  EIGENVECTORS  ARE  DELIVERED  IN THE
                   COLUMNS OF ARRAY VEC;
                   AN EIGENVECTOR, CORRESPONDING  TO A REAL  EIGENVALUE
                   GIVEN IN ARRAY RE, IS DELIVERED IN THE CORRESPONDING
                   COLUMN OF ARRAY VEC;
                   THE  REAL  AND  IMAGINARY  PART  OF AN  EIGENVECTOR,
                   CORRESPONDING  TO THE  FIRST  MEMBER  OF  A  NONREAL
                   COMPLEX CONJUGATE PAIR  OF EIGENVALUES  GIVEN IN THE
                   ARRAYS RE, IM, ARE DELIVERED  IN THE TWO CONSECUTIVE
                   COLUMNS OF ARRAY VEC CORRESPONDING TO THIS PAIR (THE
                   EIGENVECTORS  CORRESPONDING TO THE SECOND MEMBERS OF
                   NONREAL  COMPLEX CONJUGATE PAIRS  ARE NOT DELIVERED,
                   SINCE THEY ARE SIMPLY THE COMPLEX CONJUGATE OF THOSE
                   CORRESPONDING TO THE FIRST MEMBER OF SUCH PAIRS);

    MOREOVER:
    COMEIG1   DELIVERS 0, PROVIDED THAT THE PROCESS IS COMPLETED WITHIN
    EM[4]  ITERATIONS; OTHERWISE  COMEIG1  DELIVERS  K, THE  NUMBER  OF
    EIGENVALUES AND EIGENVECTORS NOT CALCULATED.

PROCEDURES USED:
    EQILBR     = CP34173,
    TFMREAHES  = CP34170,
    BAKREAHES2 = CP34172,
    BAKLBR     = CP34174,
    REAVECHES  = CP34181,
    COMVALQRI  = CP34190,
    COMVECHES  = CP34191,
    COMSCL     = CP34193.

REQUIRED CENTRAL MEMORY:
    EXECUTION FIELD LENGTH: N * N + 5N.

RUNNING TIME: ROUGHLY PROPORTIONAL TO N CUBED.

LANGUAGE: ALGOL 60.

METHOD AND PERFORMANCE:
    THE GIVEN MATRIX IS EQUILIBRATED  BY  CALLING  EQILBR (SEE  SECTION
    3.2.1.1.1) AND TRANSFORMED TO A SIMILAR UPPER-HESSENBERG   MATRIX
    BY CALLING TFMREAHES  (SEE  SECTION 3.2.1.2.1.2). THE   EIGENVALUES
    ARE THEN CALCULATED BY CALLING COMVALQRI, WHICH USES DOUBLE  QR
    ITERATION (SEE SECTION 3.3.1.2.1).
    FURTHERMORE, TO FIND THE EIGENVECTORS WILKINSON'S DEVICE IS FIRST
    APPLIED [2, P.328 AND 628]. SUBSEQUENTLY  THE  EIGENVECTORS  OF THE
    UPPER-HESSENBERG MATRIX ARE COMPUTED BY CALLING REAVECHES FOR THE
    REAL EIGENVALUES AND  COMVECHES FOR THE OTHERS (SECTION 3.3.1.2.1.)
    THE COMPUTED VECTORS ARE THEN BACK-TRANSFORMED TO THE CORRESPONDING
    EIGENVECTORS OF THE GIVEN MATRIX BY CALLING BAKREAHES2  AND  BAKLBR
    (SEE SECTIONS 3.2.1.2.1.2 AND 3.2.1.1.1). FINALLY  THE  APPROXIMATE
    EIGENVECTORS ARE NORMALIZED BY  CALLING  COMSCL (SEE SECTION 1.1.9)
    SUCH  THAT, IN EACH  EIGENVECTOR, AN  ELEMENT  OF  MAXIMUM  MODULUS
    EQUALS 1.
    FOR FURTHER DETAILS SEE THE GIVEN REFERENCES.

REFERENCES:
    [1]. T.J. DEKKER AND W. HOFFMANN.
         ALGOL 60 PROCEDURES IN NUMERICAL ALGEBRA, PART 2.
         MC TRACT 23, 1968, MATH. CENTR., AMSTERDAM.
    [2]. J.H. WILKINSON.
         THE ALGEBRAIC EIGENVALUE PROBLEM.
         CLARENDON PRESS, OXFORD, 1965.
    [3]. J.G. FRANCIS.
         THE QR TRANSFORMATION, PARTS 1 AND 2.
         COMP. J. 4 (1961), 265 - 271 AND 332 - 345.
    [4]. J.M. VARAH.
         EIGENVECTORS OF A REAL MATRIX BY INVERSE ITERATION.
         STANFORD UNIVERSITY, TECH. REP. NO. CS 34, 1966.

EXAMPLE OF USE:

    IN THIS SECTION  WE  ONLY  GIVE AN EXAMPLE OF USE OF THE PROCEDURES
    REAEIG3  AND COMEIGVAL, BECAUSE  A CALL OF  THE OTHER PROCEDURES IS
    ALMOST SIMILAR.

    THE EIGENVALUES AND CORRESPONDING EIGENVECTORS OF A MATRIX, STORED
    IN ARRAY A, WITH A[I,J]:= "IF" I=1 "THEN" 1 "ELSE" 1 / (I + J - 1),
    MAY BE OBTAINED BY THE PROCEDURE REAEIG3  IN THE FOLLOWING PROGRAM:

    "BEGIN" "INTEGER" I, J, M;
        "ARRAY" A, VEC[1:4,1:4], EM[0:5], VAL[1:4];

        "FOR" I:= 1, 2, 3, 4 "DO" "FOR" J:= 1, 2, 3, 4 "DO"
        A[I,J]:= "IF" I = 1 "THEN" 1 "ELSE" 1 / ( I + J - 1);
        EM[0]:= "-14; EM[2]:= "-13; EM[4]:= 40;
        M:= REAEIG3(A, 4, EM, VAL, VEC);
        OUTPUT(61, "("D, /")", M);
        "FOR" I:= M + 1 "STEP" 1 "UNTIL" 4 "DO"
        OUTPUT(61, "("/, 2(+.13D"+2D, 2B), /, 3(21B, +.13D"+2D, /)")",
        VAL[I], VEC[1,I], VEC[2,I], VEC[3,I], VEC[4,I]);
        OUTPUT(61, "("/, 2(.2D"+2D, /), ZD")", EM[1], EM[3], EM[5])
    "END"

    THE PROGRAM DELIVERS(THE RESULTS ARE CORRECT UP TO TWELVE DIGITS):

    THE NUMBER OF NOT CALCULATED EIGENVALUES: 0

    THE EIGENVALUES AND CORRESPONDING EIGENVECTORS:

    +.1886632138548"+01  +.1000000000000"+01
                         +.3942239850770"+00
                         +.2773202862566"+00
                         +.2150878672143"+00

    -.1980145931103"+00  +.1000000000000"+01
                         -.7388484093937"+00
                         -.3116238593839"+00
                         -.1475423243327"+00

    -.1228293686543"-01  -.4634736456357"+00
                         +.1000000000000"+01
                         -.1542548002737"+00
                         -.3765787365625"+00

    -.1441323817331"-03  +.1095712655340"+00
                         -.6208405341138"+00
                         +.1000000000000"+01
                         -.4887465241876"+00

    EM[1] = .40"+01
    EM[3] = .15"-14
    EM[5] = 5 .

    THE COMPLEX EIGENVALUES  OF A MATRIX  STORED IN ARRAY A  WITH N = 3
    AND THE  ROWS  (8, -1, -5), (-4, 4, -2)  AND  (18, -5, -7), MAY  BE
    OBTAINED BY THE PROCEDURE COMEIGVAL IN THE FOLLOWING PROGRAM:

    "BEGIN" "INTEGER" I, M;
        "ARRAY" A[1:3,1:3], EM[0:5], RE, IM[1:3];

        EM[0]:= "-14; EM[2]:= "-13; EM[4]:= 30;
        A[1,1]:= 8; A[1,2]:= -1; A[1,3]:= -5;
        A[2,1]:= -4; A[2,2]:= 4; A[2,3]:= -2;
        A[3,1]:= 18; A[3,2]:= -5; A[3,3]:= -7;
        M:= COMEIGVAL(A, 3, EM, RE, IM);
        OUTPUT(61, "("D, /")", M);
        "FOR" I:= M + 1 "STEP" 1 "UNTIL" 3 "DO"
        OUTPUT(61, "("2(+.13D"+2D, 2B), /")", RE[I], IM[I]);
        OUTPUT(61, "("/, 2(.2D"+2D, /), ZD")", EM[1], EM[3], EM[5])
    "END"

    THE PROGRAM DELIVERS(THE RESULTS ARE CORRECT UP TO TWELVE DIGITS):

    THE NUMBER OF NOT CALCULATED EIGENVALUES: 0

    THE EIGENVALUES: +.2000000000000"+01  +.4000000000000"+01
                     +.2000000000000"+01  -.4000000000000"+01
                     +.9999999999998"+00  +.0000000000000"+00

    THE ARRAY EM: EM[1] = .30"+02
                  EM[3] = .78"-17
                  EM[5] = 6 .

SOURCE TEXTS:
"CODE" 34182;
"CODE" 34184;
"CODE" 34187;
"CODE" 34192;
"CODE" 34194;