NUMAL Section 3.3.1.1.2

BEGIN SECTION : 3.3.1.1.2 (July, 1974)

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

CONTRIBUTORS: W.HOFFMANN, J.G.VERWER.

INSTITUTE:    MATHEMATICAL CENTRE.

RECEIVED:     730924.

BRIEF DESCRIPTION:
    THIS SECTION  CONTAINS SEVEN PROCEDURES.
    A) EIGVALSYM1  AND  EIGVALSYM2  CALCULATE  ALL EIGENVALUES, OR SOME
    CONSECUTIVE  EIGENVALUES  INCLUDING  THE  LARGEST, OF  A  SYMMETRIC
    MATRIX  USING  LINEAR  INTERPOLATION ON A FUNCTION  DERIVED  FROM A
    STURM SEQUENCE,
    B) EIGSYM1  AND  EIGSYM2  CALCULATE THE CORRESPONDING  EIGENVECTORS
    AS WELL, BY MEANS OF INVERSE ITERATION,
    C) QRIVALSYM1  AND  QRIVALSYM2   CALCULATE  ALL  EIGENVALUES  OF  A
    SYMMETRIC MATRIX BY MEANS OF QR ITERATION,
    D) QRISYM CALCULATES ALL EIGENVECTORS AS WELL IN THE SAME ITERATION
    PROCESS.
    EIGVALSYM1, EIGSYM1 AND QRIVALSYM1 USE IONAL  ARRAY FOR
    THE GIVEN SYMMETRIC MATRIX; THE OTHER  PROCEDURES EXPECT THE MATRIX
    TO BE STORED IN "ARRAY".
    QRISYM DELIVERS THE EIGENVECTORS IN THE ARRAY THAT WAS USED FOR THE
    ORIGINAL MATRIX IN CONTRAST WITH EIGSYM1 AND EIGSYM2 WHICH  DELIVER
    THE EIGENVECTORS IN AN EXTRA ARRAY.
    WHEN ALL EIGENVALUES HAVE TO BE CALCULATED, THE PROCEDURES USING QR
    ITERATION ARE  PREFERABLE WITH  RESPECT TO THEIR RUNNING TIME. WHEN
    ALSO  THE EIGENVECTORS  HAVE TO BE CALCULATED THE PROCEDURES  USING
    INVERSE  ITERATION ARE FASTER; HOWEVER, THE ONE USING  QR ITERATION
    USES LESS MEMORY SPACE.

KEYWORDS:
    EIGENVALUES,
    EIGENVECTORS,
    SYMMETRIC MATRIX,
    STURM-SEQUENCE,
    INVERSE ITERATION,
    QR ITERATION.


SUBSECTION: EIGVALSYM2.

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

    THE MEANING OF THE FORMAL PARAMETERS IS:
    N:      <ARITHMETIC EXPRESSION>;
            THE ORDER OF THE GIVEN MATRIX;
    NUMVAL: <ARITHMETIC EXPRESSION>;
            THE SERIAL  NUMBER OF THE LAST EIGENVALUE TO BE CALCULATED;
    A:      <ARRAY IDENTIFIER>;
            "ARRAY" A[1:N,1:N];
            ENTRY:  THE UPPER  TRIANGLE OF THE SYMMETRIC MATRIX MUST BE
                    GIVEN  IN  THE  UPPER  TRIANGULAR  PART  OF A  (THE
                    ELEMENTS A[I,J], I<= J);
            EXIT:   THE  DATA  FOR  HOUSEHOLDER'S  BACK  TRANSFORMATION
                    (WHICH ISN'T USED BY THIS PROCEDURE)  IS  DELIVERED
                    IN THE UPPER TRIANGULAR PART OF A;
            THE ELEMENTS A[I,J] FOR I > J ARE NEITHER USED NOR CHANGED;
    VAL:    <ARRAY IDENTIFIER>;
            "ARRAY" VAL[1:NUMVAL];
            EXIT:   THE  NUMVAL  LARGEST  EIGENVALUES IN  MONOTONICALLY
                    NON-INCREASING ORDER;
    EM:     <ARRAY IDENTIFIER>;
            "ARRAY" EM[0:3];
            ENTRY:  EM[0], THE MACHINE PRECISION,
                    EM[2], THE RELATIVE TOLERANCE FOR THE EIGENVALUES;
            EXIT:   EM[1], THE INFINITY NORM OF THE ORIGINAL MATRIX,
                    EM[3], THE   NUMBER   OF   ITERATIONS    USED   FOR
                           CALCULATING THE NUMVAL EIGENVALUES.

PROCEDURES USED:
    TFMSYMTRI2      =      CP34140,
    VALSYMTRI       =      CP34151.

REQUIRED CENTRAL MEMORY:
    EXECUTION FIELD LENGTH:
            THREE ONE-DIMENSIONAL REAL ARRAYS OF LENGTH N ARE USED.

RUNNING TIME:
    ROUGHLY PROPORTIONAL TO N CUBED.

LANGUAGE:   ALGOL 60.

METHOD AND PERFORMANCE:
    THE BODY OF  EIGVALSYM2  CONSISTS OF TWO PROCEDURE  STATEMENTS; THE
    FIRST IS A CALL OF  TFMSYMTRI2 TO  TRANSFORM THE  SYMMETRIC  MATRIX
    INTO  A  SIMILAR  TRIDIAGONAL  MATRIX  BY  MEANS  OF  HOUSEHOLDER'S
    TRANSFORMATION; THE SECOND IS A CALL OF VALSYMTRI TO  CALCULATE THE
    DESIRED EIGENVALUES. OPERATION DETAILS OF BOTH PROCEDURES ARE GIVEN
    IN THEIR DESCRIPTION.


SUBSECTION: EIGVALSYM1.

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

    THE MEANING OF THE FORMAL PARAMETERS IS:
    N:      <ARITHMETIC EXPRESSION>;
            THE ORDER OF THE GIVEN MATRIX;
    NUMVAL: <ARITHMETIC EXPRESSION>;
            THE SERIAL  NUMBER OF THE LAST EIGENVALUE TO BE CALCULATED;
    A:      <ARRAY IDENTIFIER>;
            "ARRAY" A[1:(N+1)*N//2];
            ENTRY:  THE UPPER  TRIANGLE OF THE SYMMETRIC MATRIX MUST BE
                    GIVEN  IN SUCH A WAY  THAT THE (I,J)-TH  ELEMENT OF
                    THE MATRIX IS A[(J-1)*J//2+I],  1 <= I <= J <= N;
            EXIT:   THE  DATA  FOR  HOUSEHOLDER'S  BACK  TRANSFORMATION
                    (WHICH ISN'T USED BY THIS PROCEDURE).
    VAL:    <ARRAY IDENTIFIER>;
            "ARRAY" VAL[1:NUMVAL];
            EXIT:   THE  NUMVAL  LARGEST  EIGENVALUES IN  MONOTONICALLY
                    NON-INCREASING ORDER;
    EM:     <ARRAY IDENTIFIER>;
            "ARRAY" EM[0:3];
            ENTRY:  EM[0], THE MACHINE PRECISION,
                    EM[2], THE RELATIVE TOLERANCE FOR THE EIGENVALUES;
            EXIT:   EM[1], THE INFINITY NORM OF THE ORIGINAL MATRIX,
                    EM[3], THE   NUMBER   OF   ITERATIONS    USED   FOR
                           CALCULATING THE NUMVAL EIGENVALUES.

PROCEDURES USED:
    TFMSYMTRI1      =      CP34143,
    VALSYMTRI       =      CP34151.

REQUIRED CENTRAL MEMORY:
    EXECUTION FIELD LENGTH:
            THREE ONE-DIMENSIONAL REAL ARRAYS OF LENGTH N ARE USED.

RUNNING TIME:
    ROUGHLY PROPORTIONAL TO N CUBED.

LANGUAGE:   ALGOL 60.

METHOD AND PERFORMANCE:
    THE BODY OF  EIGVALSYM1  CONSISTS OF TWO PROCEDURE  STATEMENTS; THE
    FIRST IS A CALL OF  TFMSYMTRI1 TO  TRANSFORM THE  SYMMETRIC  MATRIX
    INTO  A  SIMILAR  TRIDIAGONAL  MATRIX  BY  MEANS  OF  HOUSEHOLDER'S
    TRANSFORMATION; THE SECOND IS A CALL OF VALSYMTRI TO  CALCULATE THE
    DESIRED EIGENVALUES. OPERATION DETAILS OF BOTH PROCEDURES ARE GIVEN
    IN THEIR DESCRIPTION.


SUBSECTION: EIGSYM2.

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

    THE MEANING OF THE FORMAL PARAMETERS IS:
    N:      <ARITHMETIC EXPRESSION>;
            THE ORDER OF THE GIVEN MATRIX;
    NUMVAL: <ARITHMETIC EXPRESSION>;
            THE SERIAL  NUMBER OF THE LAST EIGENVALUE TO BE CALCULATED;
    A:      <ARRAY IDENTIFIER>;
            "ARRAY" A[1:N,1:N];
            ENTRY:  THE UPPER  TRIANGLE OF THE SYMMETRIC MATRIX MUST BE
                    GIVEN  IN  THE  UPPER  TRIANGULAR  PART  OF A  (THE
                    ELEMENTS A[I,J], I<= J);
            EXIT:   THE  DATA  FOR  HOUSEHOLDER'S  BACK  TRANSFORMATION
                    IS DELIVERED IN THE UPPER TRIANGULAR PART OF A;
            THE ELEMENTS A[I,J] FOR I > J ARE NEITHER USED NOR CHANGED;
    VAL:    <ARRAY IDENTIFIER>;
            "ARRAY" VAL[1:NUMVAL];
            EXIT:   THE  NUMVAL  LARGEST  EIGENVALUES IN  MONOTONICALLY
                    NON-INCREASING ORDER;
    VEC:    <ARRAY IDENTIFIER>;
            "ARRAY" VEC[1:N,1:NUMVAL];
            EXIT:   THE NUMVAL CALCULATED  EIGENVECTORS, STORED COLUMN-
                    WISE, CORRESPONDING TO THE CALCULATED  EIGENVALUES;
    EM:     <ARRAY IDENTIFIER>;
            "ARRAY" EM[0:9];
            ENTRY:  EM[0], THE MACHINE PRECISION,
                    EM[2], THE RELATIVE TOLERANCE FOR THE EIGENVALUES,
                    EM[4], THE ORTHOGONALISATION  PARAMETER (SEE METHOD
                           AND PERFORMANCE),
                    EM[6], THE TOLERANCE FOR THE EIGENVECTORS,
                    EM[8], THE  MAXIMUM  NUMBER OF  INVERSE  ITERATIONS
                           ALLOWED  FOR THE CALCULATION OF EACH  EIGEN-
                           VECTOR;
            EXIT:   EM[1], THE INFINITY NORM OF THE MATRIX,
                    EM[3], THE   NUMBER   OF   ITERATIONS    USED   FOR
                           CALCULATING THE NUMVAL EIGENVALUES,
                    EM[5], THE  NUMBER OF EIGENVECTORS  INVOLVED IN THE
                           LAST GRAM-SCHMIDT ORTHOGONALISATION,
                    EM[7], THE MAXIMUM  EUCLIDEAN  NORM OF THE RESIDUES
                           OF THE CALCULATED EIGENVECTORS,
                    EM[9], THE  LARGEST  NUMBER  OF INVERSE  ITERATIONS
                           PERFORMED FOR THE CALCULATION OF SOME EIGEN-
                           VECTOR.

PROCEDURES USED:
    TFMSYMTRI2      =      CP34140,
    VALSYMTRI       =      CP34151,
    VECSYMTRI       =      CP34152,
    BAKSYMTRI2      =      CP34141.

REQUIRED CENTRAL MEMORY:
    EXECUTION FIELD LENGTH:
            THREE ONE-DIMENSIONAL REAL ARRAYS OF LENGTH N ARE DECLARED;
            MOREOVER,  VECSYMTRI USES FIVE ONE-DIMENSIONAL  REAL ARRAYS
            OF LENGTH N AND ONE BOOLEAN ARRAY OF LENGTH N.

RUNNING TIME:
    ROUGHLY PROPORTIONAL TO N CUBED.

LANGUAGE:   ALGOL 60.

METHOD AND PERFORMANCE:
    THE  BODY  OF  EIGSYM2  CONSISTS  OF  FOUR  PROCEDURE   STATEMENTS;
    THE  FIRST  IS  A  CALL  OF TFMSYMTRI2 TO  TRANSFORM THE  SYMMETRIC
    MATRIX  INTO   A   SIMILAR   TRIDIAGONAL   MATRIX   BY   MEANS   OF
    HOUSEHOLDERS TRANSFORMATION,
    THE  SECOND  IS  A  CALL  OF  VALSYMTRI  TO  CALCULATE THE  DESIRED
    EIGENVALUES,
    THE  THIRD  IS A CALL OF VECSYMTRI  TO CALCULATE  THE CORRESPONDING
    EIGENVECTORS AND
    THE   FOURTH   IS  A   CALL  OF  BAKSYMTRI2  TO  PERFORM  THE  BACK
    TRANSFORMATION.
    THE  PARAMETERS  EM[5], EM[7] AND EM[9] ARE  GIVEN ITS VALUE IN THE
    PROCEDURE  VECSYMTRI.  FOR A POSSIBLY  SUBSEQUENT CALL OF VECSYMTRI
    THE VALUE OF EM[5] IS NEEDED.  WHEN CONSECUTIVE EIGENVALUES ARE TOO
    CLOSE TOGETHER, THE CORRESPONDING  EIGENVECTORS ARE NOT NECESSARILY
    DELIVERED ORTHOGONAL BY INVERSE ITERATION (THE METHOD WHICH IS USED
    IN VECSYMTRI). THEREFORE GRAM-SCHMIDT  ORTHOGONALISATION IS APPLIED
    ON  THE  EIGENVECTORS  WHEN  THE  DISTANCE  BETWEEN TWO CONSECUTIVE
    EIGENVALUES IS SMALLER THAN EM[4].
    FOR FURTHER  DETAILS  ONE IS REFERRED  TO  THE  SPECIFIC  PROCEDURE
    DESCRIPTIONS.


SUBSECTION: EIGSYM1.

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

    THE MEANING OF THE FORMAL PARAMETERS IS:
    N:      <ARITHMETIC EXPRESSION>;
            THE ORDER OF THE GIVEN MATRIX;
    NUMVAL: <ARITHMETIC EXPRESSION>;
            THE SERIAL  NUMBER OF THE LAST EIGENVALUE TO BE CALCULATED;
    A:      <ARRAY IDENTIFIER>;
            "ARRAY" A[1:(N+1)*N//2];
            ENTRY:  THE UPPER  TRIANGLE OF THE SYMMETRIC MATRIX MUST BE
                    GIVEN  IN SUCH A WAY  THAT THE (I,J)-TH  ELEMENT OF
                    THE MATRIX IS A[(J-1)*J//2+I],  1 <= I <= J <= N;
            EXIT:   THE  DATA  FOR  HOUSEHOLDER'S  BACK TRANSFORMATION;
    VAL:    <ARRAY IDENTIFIER>;
            "ARRAY" VAL[1:NUMVAL];
            EXIT:   THE  NUMVAL  LARGEST  EIGENVALUES IN  MONOTONICALLY
                    NON-INCREASING ORDER;
    VEC:    <ARRAY IDENTIFIER>;
            "ARRAY" VEC[1:N,1:NUMVAL];
            EXIT:   THE NUMVAL CALCULATED  EIGENVECTORS, STORED COLUMN-
                    WISE, CORRESPONDING TO THE CALCULATED  EIGENVALUES;
    EM:     <ARRAY IDENTIFIER>;
            "ARRAY" EM[0:9];
            ENTRY:  EM[0], THE MACHINE PRECISION,
                    EM[2], THE RELATIVE TOLERANCE FOR THE EIGENVALUES,
                    EM[4], THE ORTHOGONALISATION  PARAMETER (SEE METHOD
                           AND PERFORMANCE),
                    EM[6], THE TOLERANCE FOR THE EIGENVECTORS,
                    EM[8], THE  MAXIMUM  NUMBER OF  INVERSE  ITERATIONS
                           ALLOWED  FOR THE CALCULATION OF EACH  EIGEN-
                           VECTOR;
            EXIT:   EM[1], THE INFINITY NORM OF THE MATRIX,
                    EM[3], THE   NUMBER   OF   ITERATIONS    USED   FOR
                           CALCULATING THE NUMVAL EIGENVALUES,
                    EM[5], THE  NUMBER OF EIGENVECTORS  INVOLVED IN THE
                           LAST GRAM-SCHMIDT ORTHOGONALISATION,
                    EM[7], THE MAXIMUM  EUCLIDEAN  NORM OF THE RESIDUES
                           OF THE CALCULATED EIGENVECTORS,
                    EM[9], THE  LARGEST  NUMBER  OF INVERSE  ITERATIONS
                           PERFORMED FOR THE CALCULATION OF SOME EIGEN-
                           VECTOR.

PROCEDURES USED:
    TFMSYMTRI1      =      CP34143,
    VALSYMTRI       =      CP34151,
    VECSYMTRI       =      CP34152,
    BAKSYMTRI1      =      CP34144.

REQUIRED CENTRAL MEMORY:
    EXECUTION FIELD LENGTH:
            THREE ONE-DIMENSIONAL REAL ARRAYS OF LENGTH N ARE DECLARED;
            MOREOVER,  VECSYMTRI  AND  BAKSYMTRI1 USE A TOTAL AMOUNT OF
            SIX ONE-DIMENSIONAL REAL ARRAYS OF LENGTH N AND ONE BOOLEAN
            ARRAY OF LENGTH N.

RUNNING TIME:
    ROUGHLY PROPORTIONAL TO N CUBED.

LANGUAGE:   ALGOL 60.

METHOD AND PERFORMANCE:
    THE  BODY  OF  EIGSYM1  CONSISTS  OF  FOUR  PROCEDURE   STATEMENTS;
    THE  FIRST  IS  A  CALL  OF TFMSYMTRI1 TO  TRANSFORM THE  SYMMETRIC
    MATRIX  INTO   A   SIMILAR   TRIDIAGONAL   MATRIX   BY   MEANS   OF
    HOUSEHOLDERS TRANSFORMATION,
    THE  SECOND  IS  A  CALL  OF  VALSYMTRI  TO  CALCULATE THE  DESIRED
    EIGENVALUES,
    THE  THIRD  IS A CALL OF VECSYMTRI  TO CALCULATE  THE CORRESPONDING
    EIGENVECTORS AND
    THE   FOURTH   IS  A   CALL  OF  BAKSYMTRI1  TO  PERFORM  THE  BACK
    TRANSFORMATION.
    FOR DETAILS ONE IS REFERRED  TO EIGSYM2  OR TO THE  DESCRIPTIONS OF
    THE FOUR PROCEDURES USED.


SUBSECTION: QRIVALSYM2.

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

    THE MEANING OF THE FORMAL PARAMETERS IS:
    N:      <ARITHMETIC EXPRESSION>;
            THE ORDER OF THE GIVEN MATRIX;
    A:      <ARRAY IDENTIFIER>;
            "ARRAY" A[1:N,1:N];
            ENTRY:  THE UPPER  TRIANGLE OF THE SYMMETRIC MATRIX MUST BE
                    GIVEN  IN  THE  UPPER  TRIANGULAR  PART  OF A  (THE
                    ELEMENTS A[I,J], I<= J);
            EXIT:   THE  DATA  FOR  HOUSEHOLDER'S  BACK  TRANSFORMATION
                    (WHICH ISN'T USED BY THIS PROCEDURE)  IS  DELIVERED
                    IN THE UPPER TRIANGULAR PART OF A;
            THE ELEMENTS A[I,J] FOR I > J ARE NEITHER USED NOR CHANGED;
    VAL:    <ARRAY IDENTIFIER>;
            "ARRAY" VAL[1:N];
            EXIT:   THE EIGENVALUES OF THE  MATRIX  IN  SOME  ARBITRARY
                    ORDER;
    EM:     <ARRAY IDENTIFIER>;
            "ARRAY" EM[0:5];
            ENTRY:  EM[0], THE MACHINE PRECISION,
                    EM[2], THE RELATIVE  TOLERANCE FOR THE EIGENVALUES,
                    EM[4], THE MAXIMUM ALLOWED NUMBER OF ITERATIONS;
            EXIT:   EM[1], THE INFINITY NORM OF THE MATRIX,
                    EM[3], THE MAXIMUM ABSOLUTE VALUE OF THE CODIAGONAL
                           ELEMENTS NEGLECTED,
                    EM[5], THE NUMBER OF ITERATIONS PERFORMED;
    MOREOVER:
            QRIVALSYM2:=   THE NUMBER OF  EIGENVALUES  NOT  CALCULATED.

PROCEDURES USED:
    TFMSYMTRI2      =      CP34140,
    QRIVALSYMTRI    =      CP34160.

REQUIRED CENTRAL MEMORY:
    EXECUTION FIELD LENGTH:
            IONAL REAL ARRAYS OF LENGTH N ARE USED.

RUNNING TIME:
    ROUGHLY PROPORTIONAL TO N CUBED.

LANGUAGE:   ALGOL 60.

METHOD AND PERFORMANCE:
    THE BODY OF  QRIVALSYM2  CONSISTS OF TWO PROCEDURE  STATEMENTS; THE
    FIRST IS A CALL OF  TFMSYMTRI2 TO  TRANSFORM THE  SYMMETRIC  MATRIX
    INTO  A  SIMILAR  TRIDIAGONAL  MATRIX  BY  MEANS  OF  HOUSEHOLDER'S
    TRANSFORMATION; THE SECOND IS A CALL OF  QRIVALSYMTRI TO  CALCULATE
    THE  EIGENVALUES.  WHEN  THE  PROCESS  IS  COMPLETED  WITHIN  EM[4]
    ITERATIONS  THEN QRIVALSYM2:= 0; OTHERWISE QRIVALSYM2:= THE NUMBER,
    K, OF  EIGENVALUES  NOT  CALCULATED, EM[5]:= EM[4] + 1 AND ONLY THE
    LAST N - K ELEMENTS OF VAL ARE APPROXIMATE EIGENVALUES OF THE GIVEN
    MATRIX. OPERATION  DETAILS OF BOTH  PROCEDURES  USED  ARE  GIVEN IN
    THEIR DESCRIPTION.


SUBSECTION: QRIVALSYM1.

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

    THE MEANING OF THE FORMAL PARAMETERS IS:
    N:      <ARITHMETIC EXPRESSION>;
            THE ORDER OF THE GIVEN MATRIX;
    A:      <ARRAY IDENTIFIER>;
            "ARRAY" A[1:(N+1)*N//2];
            ENTRY:  THE UPPER  TRIANGLE OF THE SYMMETRIC MATRIX MUST BE
                    GIVEN  IN SUCH A WAY  THAT THE (I,J)-TH  ELEMENT OF
                    THE MATRIX IS A[(J-1)*J//2+I],  1 <= I <= J <= N;
            EXIT:   THE  DATA  FOR  HOUSEHOLDER'S  BACK  TRANSFORMATION
                    (WHICH ISN'T USED BY THIS PROCEDURE).
    VAL:    <ARRAY IDENTIFIER>;
            "ARRAY" VAL[1:N];
            EXIT:   THE EIGENVALUES OF THE  MATRIX  IN  SOME  ARBITRARY
                    ORDER;
    EM:     <ARRAY IDENTIFIER>;
            "ARRAY" EM[0:5];
            ENTRY:  EM[0], THE MACHINE PRECISION,
                    EM[2], THE RELATIVE  TOLERANCE FOR THE EIGENVALUES,
                    EM[4], THE MAXIMUM ALLOWED NUMBER OF ITERATIONS;
            EXIT:   EM[1], THE INFINITY NORM OF THE MATRIX,
                    EM[3], THE MAXIMUM ABSOLUTE VALUE OF THE CODIAGONAL
                           ELEMENTS NEGLECTED,
                    EM[5], THE NUMBER OF ITERATIONS PERFORMED;
    MOREOVER:
            QRIVALSYM1:=   THE NUMBER OF  EIGENVALUES  NOT  CALCULATED.

PROCEDURES USED:
    TFMSYMTRI1      =      CP34143,
    QRIVALSYMTRI    =      CP34160.

REQUIRED CENTRAL MEMORY:
    EXECUTION FIELD LENGTH:
            TWO ONE-DIMENSIONAL REAL ARRAYS OF LENGTH N ARE USED.

RUNNING TIME:
    ROUGHLY PROPORTIONAL TO N CUBED.

LANGUAGE:   ALGOL 60.

METHOD AND PERFORMANCE:
    THE BODY OF  QRIVALSYM1  CONSISTS OF TWO PROCEDURE  STATEMENTS; THE
    FIRST IS A CALL OF  TFMSYMTRI1 TO  TRANSFORM THE  SYMMETRIC  MATRIX
    INTO  A  SIMILAR  TRIDIAGONAL  MATRIX  BY  MEANS  OF  HOUSEHOLDER'S
    TRANSFORMATION; THE SECOND IS A CALL OF  QRIVALSYMTRI TO  CALCULATE
    THE  EIGENVALUES.  WHEN  THE  PROCESS  IS  COMPLETED  WITHIN  EM[4]
    ITERATIONS  THEN QRIVALSYM1:= 0; OTHERWISE QRIVALSYM1:= THE NUMBER,
    K, OF  EIGENVALUES  NOT  CALCULATED, EM[5]:= EM[4] + 1 AND ONLY THE
    LAST N - K ELEMENTS OF VAL ARE APPROXIMATE EIGENVALUES OF THE GIVEN
    MATRIX. OPERATION  DETAILS OF BOTH  PROCEDURES  USED  ARE  GIVEN IN
    THEIR DESCRIPTION.


SUBSECTION: QRISYM.

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

    THE MEANING OF THE FORMAL PARAMETERS IS:
    N:      <ARITHMETIC EXPRESSION>;
            THE ORDER OF THE GIVEN MATRIX;
    A:      <ARRAY IDENTIFIER>;
            "ARRAY" A[1:N,1:N];
            ENTRY:  THE UPPER  TRIANGLE OF THE SYMMETRIC MATRIX MUST BE
                    GIVEN  IN  THE  UPPER  TRIANGULAR  PART  OF A  (THE
                    ELEMENTS A[I,J], I<= J);
            EXIT:   THE  EIGENVECTORS OF THE  SYMMETRIC  MATRIX, STORED
                    COLUMNWISE;
    VAL:    <ARRAY IDENTIFIER>;
            "ARRAY" VAL[1:N];
            EXIT:   THE EIGENVALUES OF THE MATRIX  CORRESPONDING TO THE
                    CALCULATED EIGENVECTORS;
    EM:     <ARRAY IDENTIFIER>;
            "ARRAY" EM[0:5];
            ENTRY:  EM[0], THE MACHINE PRECISION,
                    EM[2], THE RELATIVE TOLERANCE FOR THE QR ITERATION,
                    EM[4], THE MAXIMUM ALLOWED NUMBER OF ITERATIONS;
            EXIT:   EM[1], THE INFINITY NORM OF THE MATRIX,
                    EM[3], THE MAXIMUM ABSOLUTE VALUE OF THE CODIAGONAL
                           ELEMENTS NEGLECTED,
                    EM[5], THE NUMBER OF ITERATIONS PERFORMED;
    MOREOVER:
            QRISYM  :=     THE NUMBER OF  EIGENVALUES AND  -VECTORS NOT
                           CALCULATED.

PROCEDURES USED:
    TFMSYMTRI2      =      CP34140,
    TFMPREVEC       =      CP34142,
    QRISYMTRI       =      CP34161.

REQUIRED CENTRAL MEMORY:

    EXECUTION FIELD LENGTH:
            TWO ONE-DIMENSIONAL REAL ARRAYS OF LENGTH N ARE USED.

RUNNING TIME:
    PROPORTIONAL TO N CUBED.

LANGUAGE:   ALGOL 60.

METHOD AND PERFORMANCE:
    THE BODY OF  QRISYM  CONSISTS OF  THREE  PROCEDURE  STATEMENTS; THE
    FIRST IS A CALL OF  TFMSYMTRI2 TO  TRANSFORM THE  SYMMETRIC  MATRIX
    INTO  A  SIMILAR  TRIDIAGONAL  MATRIX  BY  MEANS  OF  HOUSEHOLDER'S
    TRANSFORMATION,  THE  SECOND IS A CALL OF TFMPREVEC  TO PERFORM THE
    DESIRED  BACK TRANSFORMATION ON THE EIGENVECTORS IN ADVANCE AND THE
    THIRD IS A CALL OF  QRISYMTRI  TO  CALCULATE  THE  EIGENVALUES  AND
    THE EIGENVECTORS.  WHEN  THE  PROCESS  IS  COMPLETED  WITHIN  EM[4]
    ITERATIONS  THEN  QRISYM:= 0;  OTHERWISE QRISYM:= THE NUMBER, K, OF
    EIGENVALUES AND -VECTORS NOT CALCULATED, EM[5]:= EM[4] + 1 AND ONLY
    THE LAST N - K  ELEMENTS OF VAL AND THE LAST N - K COLUMNS OF A ARE
    APPROXIMATE EIGENVALUES AND EIGENVECTORS  RESPECTIVELY OF THE GIVEN
    MATRIX. OPERATION  DETAILS OF  THE  PROCEDURES  USED  ARE  GIVEN IN
    THEIR DESCRIPTION.

EXAMPLES OF USE:

    THE TWO LARGEST  EIGENVALUES IN MONOTONICALLY NON INCREASING  ORDER
    AND THE  CORRESPONDING  EIGENVECTORS OF M, WITH N = 4 AND  M[I,J] =
    1 / (I + J - 1), MAY BE OBTAINED BY THE FOLLOWING PROGRAM:

    "BEGIN"
        "INTEGER" I, J;
        "ARRAY" A[1:10], VAL[1:2], EM[0:9], VEC[1:4,1:2];

        EM[0]:= "-14; EM[2]:= "-12; EM[4]:= "-3;
        EM[6]:= 10"-10; EM[8]:= 5;
        "FOR" I:= 1 "STEP" 1 "UNTIL" 4 "DO"
        "FOR" J:= I "STEP" 1 "UNTIL" 4 "DO"
        A[(J * J - J) / 2 + I]:= 1 / (I + J - 1);
        EIGSYM1(A, 4, 2, VAL, VEC, EM);
        OUTPUT(61, "("2(+.13D"+2D, 2B), 2/")", VAL[1], VAL[2]);
        "FOR" I:= 1, 2, 3, 4 "DO"
        OUTPUT(61, "("2(+.13D"+2D, 2B), /")", VEC[I,1], VEC[I,2]);
        OUTPUT(61, "("2(.2D"+2D, /), 3(2ZD, /)")",
        EM[1], EM[7], EM[3], EM[5], EM[9])
    "END"

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

    THE EIGENVALUES:  +.1500214280059"+01  +.1691412202214"+00

    THE EIGENVECTORS: -.7926082911638"+00  +.5820756994972"+00
                      -.4519231209016"+00  -.3705021850671"+00
                      -.3224163985818"+00  -.5095786345018"+00
                      -.2521611696882"+00  -.5140482722222"+00

    EM[1] = .21"+01
    EM[7] = .92"-14
    EM[3] = 32
    EM[5] = 1
    EM[9] = 1 .

    THE  TWO  LARGEST  EIGENVALUES  OF  M, WITH  N = 4   AND   M[I,J] =
    1 / (I + J - 1), MAY  BE  OBTAINED IN  MONOTONICALLY NON INCREASING
    ORDER BY THE FOLLOWING PROGRAM:

    "BEGIN"
        "INTEGER" I, J;
        "ARRAY" A[1:4,1:4], VAL[1:2], EM[0:3];

        EM[0]:= "-14; EM[2]:= "-12;
        "FOR" I:= 1 "STEP" 1 "UNTIL" 4 "DO"
        "FOR" J:= I "STEP" 1 "UNTIL" 4 "DO" A[I,J]:= 1 / (I + J -1);
        EIGVALSYM2(A, 4, 2, VAL, EM);
        OUTPUT(61, "("2(+.13D"+2D, 2B)")", VAL[1], VAL[2]);
        OUTPUT(61, "("2/, .2D"+2D, /, 2ZD")", EM[1], EM[3])
    "END"

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

    THE EIGENVALUES: +.1500214280059"+01  +.1691412202214"+00
    EM[3] = .21"+01
    EM[1] = 32  .

SOURCE TEXT(S) :
"CODE" 34153;
"CODE" 34154;
"CODE" 34155;
"CODE" 34156;
"CODE" 34162;
"CODE" 34163;
"CODE" 34164;