1SECTION:3.2.1.2.1.1          (DECEMBER 1979)                     PAGE 1
 
 
 
 AUTHORS:      T.J.DEKKER, W.HOFFMANN.
 
 
 CONTRIBUTORS: W.HOFFMANN, J.G.VERWER.
 
 
 INSTITUTE:    MATHEMATICAL CENTRE.
 
 
 RECEIVED:     730705.
 
 
 BRIEF DESCRIPTION:
     THIS SECTION CONTAINS FIVE PROCEDURES.
     A) TFMSYMTRI2 AND TFMSYMTRI1 TRANSFORM A REAL SYMMETRIC MATRIX INTO
     A SIMILAR TRIDIAGONAL ONE BY MEANS OF HOUSEHOLDER'S TRANSFORMATION,
     B) BAKSYMTRI2   AND   BAKSYMTRI1  PERFORM  THE  CORRESPONDING  BACK
     TRANSFORMATION  AND  FINALLY,
     C) TFMPREVEC (WHICH IS TO BE USED IN COMBINATION  WITH  TFMSYMTRI2)
     CALCULATES THE TRANSFORMING MATRIX.
     TFMSYMTRI2  AND  BAKSYMTRI2  USE  THE  UPPER  TRIANGLE  OF  A  TWO-
     DIMENSIONAL  ARRAY  FOR THE UPPER  TRIANGLE OF THE GIVEN  SYMMETRIC
     MATRIX  (TFMSYMTRI2)  OR  FOR  THE  DATA  FOR   HOUSEHOLDER'S  BACK
     TRANSFORMATION  (BAKSYMTRI2).  THE OTHER  ELEMENTS ARE NEITHER USED
     NOR CHANGED.
     TFMSYMTRI1 AND BAKSYMTRI1 USE AN ARRAY A[1:(N+1)*N//2] FOR THE GIVEN
     SYMMETRIC  MATRIX (TFMSYMTRI1)  OR  FOR THE DATA FOR  HOUSEHOLDER'S
     TRANSFORMATION (BAKSYMTRI1).
 
 KEYWORDS:
     HOUSEHOLDER'S TRANSFORMATION,
     TRIANGULARIZATION.
 
 
1SECTION:3.2.1.2.1.1          (DECEMBER 1979)                     PAGE 2
 
 
 
 SUBSECTION: TFMSYMTRI2.
 
 CALLING SEQUENCE:
     THE HEADING OF THIS PROCEDURE IS:
     "PROCEDURE" TFMSYMTRI2(A, N, D, B, BB, EM); "VALUE" N;
     "INTEGER" N; "ARRAY" A, D, B, BB, EM; "CODE" 34140;
 
     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 IS
                     DELIVERED IN THE  UPPER  TRIANGULAR  PART OF A. THE
                     ELEMENTS  A[I,J], I>J ARE NEITHER USED NOR CHANGED;
     D:      <ARRAY IDENTIFIER>;
             "ARRAY"D[1:N];
             EXIT:   THE  MAIN  DIAGONAL  OF THE  SYMMETRIC  TRIDIAGONAL
                     MATRIX   T  (SAY),   PRODUCED   BY  HOUSEHOLDER'S
                     TRANSFORMATION;
     B:      <ARRAY IDENTIFIER>;
             "ARRAY" B[1:N];
             EXIT:   THE CODIAGONAL  ELEMENTS OF T ARE DELIVERED IN B[1]
                     THROUGH B[N-1]; B[N] IS SET EQUAL TO ZERO;
     BB:     <ARRAY IDENTIFIER>;
             "ARRAY"BB[1:N];
             EXIT:   THE  SQUARES  OF THE  CODIAGONAL  ELEMENTS OF T ARE
                     DELIVERED  IN  BB[1] THROUGH BB[N-1];  BB[N] IS SET
                     EQUAL TO ZERO;
     EM:     <ARRAY IDENTIFIER>;
             "ARRAY"EM[0:1];
             ENTRY:  EM[0], THE MACHINE PRECISION;
             EXIT:   EM[1], THE  INFINITY  NORM OF THE ORIGINAL  MATRIX.
 
 PROCEDURES USED:
     TAMVEC   =      CP34012,
     MATMAT   =      CP34013,
     TAMMAT   =      CP34014,
     ELMVECCOL=      CP34021,
     ELMCOLVEC=      CP34022,
     ELMCOL   =      CP34023.
 
 
 RUNNING TIME: ROUGHLY PROPORTIONAL TO N CUBED.
 
 
 
1SECTION:3.2.1.2.1.1          (DECEMBER 1979)                     PAGE 3
 
 
 
 METHOD AND PERFORMANCE:
     A  GIVEN  SYMMETRIC  MATRIX M IS TRANSFORMED INTO A TRIDIAGONAL ONE
     BY  MEANS OF N-1 ORTHOGONAL  SIMILARITY  TRANSFORMATIONS; THE  P-TH
     TRANSFORMATION  IS  CHOSEN  IN  SUCH  A  WAY THAT IN THE (N-P+1)-TH
     COLUMN AND ROW OF M THE DESIRED  ZEROES  ARE  INTRODUCED.  HOWEVER,
     IF, IN THIS COLUMN AND ROW, ALL ELEMENTS OUTSIDE THE MAIN  DIAGONAL
     AND THE  ADJACENT  CODIAGONALS  ARE SMALLER IN ABSOLUTE  VALUE THAN
     THE INFINITY NORM OF M TIMES THE MACHINE  PRECISION, THEN THE  P-TH
     TRANSFORMATION IS SKIPPED.
     FOR FURTHER DETAILS SEE REF[1] AND REF[2].
 
 
 SUBSECTION: BAKSYMTRI2.
 
 CALLING SEQUENCE:
     THE HEADING OF THE PROCEDURE IS:
     "PROCEDURE" BAKSYMTRI2(A, N, N1, N2, VEC); "VALUE" N, N1, N2;
     "INTEGER" N, N1, N2; "ARRAY" A, VEC; "CODE" 34141;
 
     THE MEANING OF THE FORMAL PARAMETERS IS:
     A:      <ARRAY IDENTIFIER>;
             "ARRAY"A[1:N,1:N];
             ENTRY:  THE DATA FOR THE BACK  TRANSFORMATION,  AS PRODUCED
                     BY  TFMSYMTRI2,   MUST   BE   GIVEN  IN  THE  UPPER
                     TRIANGULAR PART OF A;
     N:      <ARITHMETIC EXPRESSION>;
             THE ORDER OF THE GIVEN MATRIX;
     N1,N2:  <ARITHMETIC EXPRESSION>;
             LOWER AND UPPER BOUND, RESPECTIVELY, OF THE COLUMN NUMBERS
             OF VEC;
     VEC:    <ARRAY IDENTIFIER>;
             "ARRAY"VEC[1:N,N1:N2];
             ENTRY:  THE VECTORS ON WHICH THE BACK TRANSFORMATION HAS TO
                     BE PERFORMED;
             EXIT:   THE TRANSFORMED VECTORS.
 
 PROCEDURES USED:
     TAMMAT  =       CP34014,
     ELMCOL  =       CP34023.
 
 
 RUNNING TIME: ROUGHLY PROPORTIONAL TO N SQUARED TIMES (N2-N1+1).
 
 
 METHOD AND PERFORMANCE: SEE REF[1].
 
 
1SECTION:3.2.1.2.1.1          (DECEMBER 1979)                     PAGE 4
 
 SUBSECTION: TFMPREVEC.
 
 CALLING SEQUENCE:
     THE HEADING OF THIS PROCEDURE IS:
     "PROCEDURE" TFMPREVEC(A, N); "VALUE" N; "INTEGER" N; "ARRAY" A;
     "CODE" 34142;
 
     THE MEANING OF THE FORMAL PARAMETERS IS:
     A:      <ARRAY IDENTIFIER>;
             "ARRAY"A[1:N,1:N];
             ENTRY:  THE DATA  FOR THE BACK  TRANSFORMATION, AS PRODUCED
                     BY   TFMSYMTRI2,   MUST   BE  GIVEN  IN  THE  UPPER
                     TRIANGULAR PART OF A;
             EXIT:   THE MATRIX  WHICH  TRANSFORMS  THE ORIGINAL  MATRIX
                     INTO A SIMILAR TRIDIAGONAL ONE.
     N:      <ARITHMETIC EXPRESSION>;
             THE ORDER OF THE MATRIX;
 
 PROCEDURES USED:
     TAMMAT  =       CP34014,
     ELMCOL  =       CP34023.
 
 
 RUNNING TIME: ROUGHLY PROPORTIONAL TO N CUBED.
 
 
 METHOD AND PERFORMANCE: SEE REF[1].
 
 
1SECTION:3.2.1.2.1.1          (DECEMBER 1979)                     PAGE 5
 
 
 
 SUBSECTION: TFMSYMTRI1.
 
 
 CALLING SEQUENCE:
     THE HEADING OF THIS PROCEDURE IS:
     "PROCEDURE" TFMSYMTRI1(A, N, D, B, BB, EM); "VALUE" N;
     "INTEGER" N; "ARRAY" A, D, B, BB, EM; "CODE" 34143;
 
     THE MEANING OF THE FORMAL PARAMETERS IS:
     A:      <ARRAY IDENTIFIER>;
             "ARRAY"A[1:(N+1)*N//2];
             ENTRY:  THE  UPPER  TRIANGLE OF THE GIVEN  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 AS
                     USED BY BAKSYMTRI1;
     N:      <ARITHMETIC EXPRESSION>;
             THE ORDER OF THE GIVEN MATRIX;
     D:      <ARRAY IDENTIFIER>;
             "ARRAY"D[1:N];
             EXIT:   THE  MAIN  DIAGONAL  OF THE  SYMMETRIC  TRIDIAGONAL
                     MATRIX   T  (SAY),   PRODUCED   BY  HOUSEHOLDER'S
                     TRANSFORMATION;
     B:      <ARRAY IDENTIFIER>;
             "ARRAY" B[1:N];
             EXIT:   THE CODIAGONAL  ELEMENTS OF T ARE DELIVERED IN B[1]
                     THROUGH B[N-1]; B[N] IS SET EQUAL TO ZERO;
     BB:     <ARRAY IDENTIFIER>;
             "ARRAY"BB[1:N];
             EXIT:   THE  SQUARES  OF THE  CODIAGONAL  ELEMENTS OF T ARE
                     DELIVERED  IN  BB[1] THROUGH BB[N-1];  BB[N] IS SET
                     EQUAL TO ZERO;
     EM:     <ARRAY IDENTIFIER>;
             "ARRAY"EM[0:1];
             ENTRY:  EM[0], THE MACHINE PRECISION;
             EXIT:   EM[1], THE  INFINITY  NORM OF THE  ORIGINAL MATRIX.
 
 PROCEDURES USED:
     VECVEC  =       CP34010,
     SEQVEC  =       CP34016,
     ELMVEC  =       CP34020.
 
 
 RUNNING TIME: ROUGHLY PROPORTIONAL TO N CUBED.
 
 
 METHOD AND PERFORMANCE: SEE TFMSYMTRI2 (THIS SECTION).
 
 
1SECTION:3.2.1.2.1.1          (DECEMBER 1979)                     PAGE 6
 
 
 
 SUBSECTION: BAKSYMTRI1.
 
 CALLING SEQUENCE:
     THE HEADING OF THE PROCEDURE IS:
     "PROCEDURE" BAKSYMTRI1(A, N, N1, N2, VEC); "VALUE" N, N1, N2;
     "INTEGER" N, N1, N2; "ARRAY" A, VEC; "CODE" 34144;
 
     THE MEANING OF THE FORMAL PARAMETERS IS:
     A:      <ARRAY IDENTIFIER>;
             "ARRAY"A[1:(N+1)*N//2];
             ENTRY:  THE DATA FOR THE BACK  TRANSFORMATION, AS  PRODUCED
                     BY  TFMSYMTRI1;
     N:      <ARITHMETIC EXPRESSION>;
             THE ORDER OF THE GIVEN MATRIX;
     N1,N2:  <ARITHMETIC EXPRESSION>;
             LOWER AND UPPER BOUND, RESPECTIVELY, OF THE COLUMN NUMBERS
             OF VEC;
     VEC:    <ARRAY IDENTIFIER>;
             "ARRAY"VEC[1:N,N1:N2];
             ENTRY:  THE VECTORS ON WHICH THE BACK TRANSFORMATION HAS TO
                     BE PERFORMED;
             EXIT:   THE TRANSFORMED VECTORS.
 
 PROCEDURES USED:
     VECVEC  =       CP34010,
     ELMVEC  =       CP34020.
 
 REQUIRED CENTRAL MEMORY:
     AN AUXILIARY  ONE-DIMENSIONAL REAL ARRAY OF LENGTH N IS USED.
 
 RUNNING TIME: ROUGHLY PROPORTIONAL TO N SQUARED TIMES (N2-N1+1).
 
 METHOD AND PERFORMANCE: SEE REF[1].
 
 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.
 
 EXAMPLE OF USE:
 
     THE   FIVE   PROCEDURES  OF  THIS   SECTION  ARE  USED  IN
     SECTION 3.3.1.1.2:
             EIGSYM2 USES TFMSYMTRI2 AND BAKSYMTRI2;
             EIGSYM1 USES TFMSYMTRI1 AND BAKSYMTRI1;
             QRISYM  USES TFMSYMTRI2 AND TFMPREVEC.
1SECTION:3.2.1.2.1.1          (JUNE 1974)                         PAGE 7
 
 
 
 SOURCE TEXT(S):
 
0"CODE" 34140;
     "COMMENT" MCA 2300;
     "PROCEDURE" TFMSYMTRI2(A, N, D, B, BB, EM); "VALUE" N;"INTEGER" N;
      "ARRAY" A, B, BB, D, EM;
     "BEGIN" "INTEGER" I, J, R, R1;
         "REAL" W, X, A1, B0, BB0, D0, MACHTOL, NORM;
 
         "REAL" "PROCEDURE" TAMMAT(L, U, I, J, A, B); "CODE" 34014;
         "REAL" "PROCEDURE" MATMAT(L, U, I, J, A, B); "CODE" 34013;
         "PROCEDURE" ELMVECCOL(L, U, I, A, B, X); "CODE" 34021;
         "REAL" "PROCEDURE" TAMVEC(L, U, I, A, B); "CODE" 34012;
         "PROCEDURE" ELMCOL(L, U, I, J, A, B, X); "CODE" 34023;
         "PROCEDURE" ELMCOLVEC(L, U, I, A, B, X); "CODE" 34022;
 
         NORM:= 0;
         "FOR" J:= 1 "STEP" 1 "UNTIL" N "DO"
         "BEGIN" W:= 0;
             "FOR" I:= 1 "STEP" 1 "UNTIL" J "DO" W:= ABS(A[I,J]) + W;
             "FOR" I:= J + 1 "STEP" 1 "UNTIL" N "DO" W:= ABS(A[J,I]) +
             W; "IF" W > NORM "THEN" NORM:= W
         "END";
         MACHTOL:= EM[0] * NORM; EM[1]:= NORM; R:= N;
         "FOR" R1:= N - 1 "STEP" -1 "UNTIL" 1 "DO"
         "BEGIN" D[R]:= A[R,R]; X:= TAMMAT(1, R - 2, R, R, A, A);
             A1:= A[R1,R]; "IF" SQRT(X) <= MACHTOL "THEN"
             "BEGIN" B0:= B[R1]:= A1; BB[R1]:= B0 * B0;A[R,R]:= 1 "END"
             "ELSE"
             "BEGIN" BB0:= BB[R1]:= A1 * A1 + X;
                 B0:= "IF" A1 > 0 "THEN" -SQRT(BB0) "ELSE" SQRT(BB0);
                 A1:= A[R1,R]:= A1 - B0; W:= A[R,R]:= 1 / (A1 * B0);
                 "FOR" J:= 1 "STEP" 1 "UNTIL" R1 "DO" B[J]:= (TAMMAT(1,
                 J, J, R, A, A) + MATMAT(J + 1, R1, J, R, A, A)) * W;
                 ELMVECCOL(1, R1, R, B, A, TAMVEC(1, R1, R, A, B) *
                 W * .5); "FOR" J:= 1 "STEP" 1 "UNTIL" R1 "DO"
                 "BEGIN" ELMCOL(1, J, J, R, A, A, B[J]);
                     ELMCOLVEC(1, J, J, A, B, A[J,R])
                 "END"; B[R1]:= B0
             "END"; R:= R1
         "END";
         D[1]:= A[1,1]; A[1,1]:= 1; B[N]:= BB[N]:= 0
     "END" TFMSYMTRI2;
         "EOP"
 
1SECTION:3.2.1.2.1.1          (JUNE 1974)                         PAGE 8
 
 
0"CODE" 34141;
     "COMMENT" MCA 2301;
     "PROCEDURE" BAKSYMTRI2(A, N, N1, N2, VEC); "VALUE" N, N1, N2;
     "INTEGER" N, N1, N2; "ARRAY" A, VEC;
     "BEGIN" "INTEGER" I, J, K; "REAL" W;
 
         "REAL" "PROCEDURE" TAMMAT(L, U, I, J, A, B); "CODE" 34014;
         "PROCEDURE" ELMCOL(L, U, I, J, A, B, X); "CODE" 34023;
 
         "FOR" J:= 2 "STEP" 1 "UNTIL" N "DO"
         "BEGIN" W:= A[J,J]; "IF" W < 0 "THEN"
             "FOR" K:= N1 "STEP" 1 "UNTIL" N2 "DO"
             ELMCOL(1, J - 1, K, J, VEC, A,
             TAMMAT(1, J - 1, J, K, A, VEC) * W)
         "END"
     "END" BAKSYMTRI2;
         "EOP"
 
0"CODE" 34142;
     "COMMENT" MCA 2302;
     "PROCEDURE" TFMPREVEC(A, N); "VALUE" N; "INTEGER" N; "ARRAY" A;
     "BEGIN" "INTEGER" I, J, J1, K; "REAL" AB;
 
         "REAL" "PROCEDURE" TAMMAT(L, U, I, J, A, B); "CODE" 34014;
         "PROCEDURE" ELMCOL(L, U, I, J, A, B, X); "CODE" 34023;
 
         J1:= 1;
         "FOR" J:= 2 "STEP" 1 "UNTIL" N "DO"
         "BEGIN" "FOR" I:= 1 "STEP" 1 "UNTIL" J1 - 1 ,
             J "STEP" 1 "UNTIL" N "DO" A[I,J1]:= 0;
             A[J1,J1]:= 1; AB:= A[J,J];
             "IF" AB < 0 "THEN"
             "FOR" K:= 1 "STEP" 1 "UNTIL" J1 "DO"
             ELMCOL(1, J1, K, J, A, A,
             TAMMAT(1, J1, J, K, A, A) * AB); J1:= J
         "END";
         "FOR" I:= N - 1 "STEP" -1 "UNTIL" 1 "DO"
         A[I,N]:= 0; A[N,N]:= 1
     "END" TFMPREVEC;
         "EOP"
 
1SECTION:3.2.1.2.1.1          (DECEMBER 1979)                     PAGE 9
 
 
0"CODE" 34143;
     "COMMENT" MCA 2305;
     "PROCEDURE" TFMSYMTRI1(A, N, D, B, BB, EM); "VALUE" N;"INTEGER" N;
      "ARRAY" A, B, BB, D, EM;
     "BEGIN" "INTEGER" I, J, R, R1, P, Q, TI, TJ;
         "REAL" S, W, X, A1, B0, BB0, D0, NORM, MACHTOL;
 
         "REAL" "PROCEDURE" VECVEC(L, U, SHIFT, A, B); "CODE" 34010;
         "REAL" "PROCEDURE" SEQVEC(L, U, IL, SHIFT, A, B);"CODE" 34016;
         "PROCEDURE" ELMVEC(L, U, SHIFT, A, B, X); "CODE" 34020;
 
         NORM:= 0; TJ:= 0;
         "FOR" J:= 1 "STEP" 1 "UNTIL" N "DO"
         "BEGIN" W:= 0;
             "FOR" I:= 1 "STEP" 1 "UNTIL" J "DO" W:= ABS(A[I + TJ]) +W;
             TJ:= TJ + J; TI:= TJ + J;
             "FOR" I:= J + 1 "STEP" 1 "UNTIL" N "DO"
             "BEGIN" W:= ABS(A[TI]) + W; TI:= TI + I "END";
             "IF" W > NORM "THEN" NORM:= W
         "END";
         MACHTOL:= EM[0] * NORM; EM[1]:= NORM; Q:= (N + 1) * N // 2;
         R:= N; "FOR" R1:= N - 1 "STEP" -1 "UNTIL" 1 "DO"
         "BEGIN" P:= Q - R; D[R]:= A[Q];
             X:= VECVEC(P + 1, Q - 2, 0, A, A);
             A1:= A[Q - 1]; "IF" SQRT(X) <= MACHTOL "THEN"
             "BEGIN" B0:= B[R1]:= A1; BB[R1]:= B0 * B0; A[Q]:= 1 "END"
             "ELSE"
             "BEGIN" BB0:= BB[R1]:= A1 * A1 + X;
                 B0:= "IF" A1 > 0 "THEN" -SQRT(BB0) "ELSE" SQRT(BB0);
                 A1:= A[Q - 1]:= A1 - B0; W:= A[Q]:= 1 / (A1 * B0);
                 TJ:= 0; "FOR" J:= 1 "STEP" 1 "UNTIL" R1 "DO"
                 "BEGIN" TI:= TJ + J; S:= VECVEC(TJ + 1, TI, P - TJ,
                     A, A); TJ:= TI + J;
                     B[J]:= (SEQVEC(J + 1, R1, TJ, P, A, A) + S) * W;
                     TJ:= TI
                 "END";
                 ELMVEC(1, R1, P, B, A, VECVEC(1,R1,P,B,A)* W *.5);
                 TJ:= 0; "FOR" J:= 1 "STEP" 1 "UNTIL" R1 "DO"
                 "BEGIN" TI:= TJ + J; ELMVEC(TJ + 1, TI, P - TJ, A, A,
                     B[J]);ELMVEC(TJ + 1, TI, -TJ, A, B, A[J + P]);
                     TJ:= TI
                 "END"; B[R1]:= B0
             "END";
             Q:= P; R:= R1
         "END";
         D[1]:= A[1]; A[1]:= 1; B[N]:= BB[N]:= 0
     "END" TFMSYMTRI1;
         "EOP"
 
1SECTION:3.2.1.2.1.1          (JUNE 1974)                        PAGE 10
 
 
0"CODE" 34144;
     "COMMENT" MCA 2306;
     "PROCEDURE" BAKSYMTRI1(A, N, N1, N2, VEC); "VALUE" N, N1, N2;
     "INTEGER" N, N1, N2; "ARRAY" A, VEC;
     "BEGIN" "INTEGER" J, J1, K, TI, TJ;
         "REAL" W; "ARRAY" AUXVEC[1:N];
 
         "REAL" "PROCEDURE" VECVEC(L, U, SHIFT, A, B); "CODE" 34010;
         "PROCEDURE" ELMVEC(L, U, SHIFT, A, B, X); "CODE" 34020;
 
         "FOR" K:= N1 "STEP" 1 "UNTIL" N2 "DO"
         "BEGIN" "FOR" J:= 1 "STEP" 1 "UNTIL" N "DO"
             AUXVEC[J]:= VEC[J,K]; TJ:= J1:= 1;
             "FOR" J:= 2 "STEP" 1 "UNTIL" N "DO"
             "BEGIN" TI:= TJ + J; W:= A[TI];
                 "IF" W < 0 "THEN" ELMVEC(1, J1, TJ, AUXVEC,A,VECVEC(1,
                 J1, TJ, AUXVEC, A) * W); J1:= J; TJ:= TI
             "END";
             "FOR" J:= 1 "STEP" 1 "UNTIL" N "DO" VEC[J,K]:= AUXVEC[J]
         "END"
     "END" BAKSYMTRI1;
         "EOP"
1SECTION:3.2.1.2.1.2          (JUNE 1974)                         PAGE 1      
        
        
        
 AUTHORS:      T.J.DEKKER AND W.HOFFMANN.       
        
        
 CONTRIBUTORS: W.HOFFMANN, J.G.VERWER.
        
        
 INSTITUTE:    MATHEMATICAL CENTRE.   
        
        
 RECEIVED:     731112.      
        
        
 BRIEF DESCRIPTION:         
        
     THIS SECTION  CONTAINS THREE PROCEDURES.   
     A) TFMREAHES TRANSFORMS A MATRIX INTO A SIMILAR  UPPER-HESSENBERG        
     MATRIX BY MEANS OF WILKINSON'S TRANSFORMATION,       
     B) BAKREAHES1  PERFORMS  THE  CORRESPONDING  BACK TRANSFORMATION ON      
     A VECTOR AND SHOULD BE CALLED AFTER TFMREAHES,       
     C) BAKREAHES2  PERFORMS  THE  CORRESPONDING  BACK TRANSFORMATION ON      
     THE COLUMNS OF A MATRIX AND SHOULD BE CALLED AFTER TFMREAHES.  
        
 KEYWORDS:        
        
     SIMILARITY TRANSFORMATION,       
     UPPER-HESSENBERG MATRIX.         
        
        
1SECTION:3.2.1.2.1.2          (DECEMBER 1975)                     PAGE 2      
        
        
        
 SUBSECTION: TFMREAHES.     
        
 CALLING SEQUENCE:
        
     THE HEADING OF THE PROCEDURE IS: 
     "PROCEDURE" TFMREAHES(A, N, EM, INT); "VALUE" N;     
     "INTEGER" N; "ARRAY" A, EM; "INTEGER" "ARRAY" INT;   
        
     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 MATRIX TO BE TRANSFORMED;        
             EXIT:   THE UPPER-HESSENBERG MATRIX IS  DELIVERED IN  THE        
                     UPPER TRIANGLE AND THE FIRST  SUBDIAGONAL OF A, THE      
                     (NONTRIVIAL ELEMENTS OF THE) TRANSFORMING   MATRIX,      
                     L, IN  THE  REMAINING  PART OF A, I.E. A[I,J] =
                     L[I,J + 1], FOR I = 3,...,N AND J = 1,...,I - 2;         
     EM:     <ARRAY IDENTIFIER>;      
             "ARRAY"EM[0:1];
             ENTRY:  EM[0], THE MACHINE PRECISION;        
             EXIT:   EM[1], THE INFINITY NORM OF THE ORIGINAL MATRIX;         
     INT:    <ARRAY IDENTIFIER>;      
             "INTEGER""ARRAY" INT[1:N];         
             EXIT:   THE PIVOTAL INDICES DEFINING  THE  STABILIZING  ROW      
                     AND COLUMN INTERCHANGES;   
        
 PROCEDURES USED: 
     MATVEC          =      CP34011,  
     MATMAT          =      CP34013,  
     ICHCOL          =      CP34031,  
     ICHROW          =      CP34032.  
        
 REQUIRED CENTRAL MEMORY:   
     EXECUTION FIELD LENGTH:
             A  ONE-DIMENSIONAL REAL ARRAY OF LENGTH N IS DECLARED. 
        
 RUNNING TIME: ROUGHLY PROPORTIONAL TO N CUBED. 
        
 LANGUAGE:   ALGOL 60.      
        
 METHOD AND PERFORMANCE:    
     WILKINSON'S   TRANSFORMATION   IS   A TRIANGULAR       SIMILARITY        
     TRANSFORMATION   WITH   STABILIZING  ROW  AND  COLUMN  INTERCHANGES      
     TRANSFORMING  A  MATRIX, M, INTO  AN UPPER-HESSENBERG  MATRIX, H.        
     THE TRANSFORMING MATRIX IS THE PRODUCT OF A PERMUTATION  MATRIX, P,      
     AND A UNIT LOWER-TRIANGULAR MATRIX, L. THE NONDIAGONAL ELEMENTS  IN      
     THE FIRST COLUMN OF L ARE 0, AND THE  ROW  AND  COLUMN INTERCHANGES      
     ARE CHOSEN IN SUCH A WAY THAT THE  ABSOLUTE  VALUE OF EACH  ELEMENT      
     OF L IS AT MOST 1.     
     BECAUSE  OF  THE SPECIAL  FORM  OF L, THE  MATRICES H AND L CAN  BE      
     STORED TOGETHER  IN  THE  ARRAY  USED FOR THE MATRIX M (SEE CALLING      
     SEQUENCE). FOR FURTHER DETAILS SEE REFERENCE [1] AND [2].      
1SECTION:3.2.1.2.1.2          (JUNE 1974)                         PAGE 3      
        
        
        
 SUBSECTION: BAKREAHES1.    
        
        
 CALLING SEQUENCE:
     THE HEADING OF THE PROCEDURE IS: 
     "PROCEDURE" BAKREAHES1(A, N, INT, V); "VALUE" N;     
     "INTEGER" N; "ARRAY" A, V; "INTEGER" "ARRAY" INT;    
        
     THE MEANING OF THE FORMAL PARAMETERS IS:   
     N:      <ARITHMETIC EXPRESSION>; 
             THE LENGTH OF THE VECTOR TO BE TRANSFORMED;  
     A:      <ARRAY IDENTIFIER>;      
             "ARRAY"A[1:N,1:N];       
             ENTRY:  THE (NONTRIVIAL  ELEMENTS  OF  THE)    TRANSFORMING      
                     MATRIX, L, AS PRODUCED BY TFMREAHES MUST  BE  GIVEN      
                     IN THE PART  BELOW  THE  FIRST  SUBDIAGONAL  OF  A,      
                     I.E. A[I,J] = L[I,J + 1], FOR I = 3,...,N   AND
                     J = 1,...,I - 2; 
     INT:    <ARRAY IDENTIFIER>;      
             "INTEGER""ARRAY" INT[1:N];         
             ENTRY:  PIVOTAL INDICES DEFINING  THE  STABILIZING  ROW AND      
                     COLUMN INTERCHANGES AS PRODUCED BY TFMREAHES;  
     V:      <ARRAY IDENTIFIER>;      
             "ARRAY"V[1:N]; 
             ENTRY:  THE VECTOR TO BE TRANSFORMED;        
             EXIT:   THE TRANSFORMED VECTOR.    
        
        
 PROCEDURES USED: 
     MATVEC          =      CP34011.  
        
        
 RUNNING TIME: ROUGHLY PROPORTIONAL TO N SQUARED.         
        
        
 LANGUAGE:   ALGOL 60.      
        
        
 METHOD AND PERFORMANCE:    
     THE    BACK  TRANSFORMATION  WHICH  CORRESPONDS  TO   WILKINSON'S        
     TRANSFORMATION AS PERFORMED BY TFMREAHES  TRANSFORMS  A  VECTOR, X,      
     INTO THE VECTOR PLX, WHERE PL IS THE TRANSFORMING MATRIX.      
        
        
1SECTION:3.2.1.2.1.2          (JUNE 1974)                         PAGE 4      
        
        
        
 SUBSECTION: BAKREAHES2.    
        
 CALLING SEQUENCE:
     THE HEADING OF THE PROCEDURE IS: 
     "PROCEDURE" BAKREAHES2(A, N, N1, N2, INT, VEC); "VALUE" N, N1, N2;       
     "INTEGER" N, N1, N2; "ARRAY" A, VEC; "INTEGER" "ARRAY" INT;    
        
     THE MEANING OF THE FORMAL PARAMETERS IS:   
     N:      <ARITHMETIC EXPRESSION>; 
             THE LENGTH OF THE VECTORS TO BE TRANSFORMED; 
     N1, N2: <ARITHMETIC EXPRESSION>; 
             THE  COLUMN  NUMBERS  OF  THE  FIRST  AND LAST VECTOR TO BE      
             TRANSFORMED;   
     A:      <ARRAY IDENTIFIER>;      
             "ARRAY"A[1:N,1:N];       
             ENTRY:  THE (NONTRIVIAL  ELEMENTS  OF  THE)    TRANSFORMING      
                     MATRIX, L, AS PRODUCED BY TFMREAHES MUST  BE  GIVEN      
                     IN THE PART  BELOW  THE  FIRST  SUBDIAGONAL  OF  A,      
                     I.E. A[I,J] = L[I,J + 1], FOR I = 3,...,N   AND
                     J = 1,...,I - 2; 
     INT:    <ARRAY IDENTIFIER>;      
             "INTEGER""ARRAY" INT[1:N];         
             ENTRY:  PIVOTAL INDICES DEFINING  THE  STABILIZING  ROW AND      
                     COLUMN INTERCHANGES AS PRODUCED BY TFMREAHES;  
     VEC:    <ARRAY IDENTIFIER>;      
             "ARRAY"VEC[1:N,N1:N2];   
             ENTRY:  THE   N2 - N1 + 1   VECTORS   OF  LENGTH  N  TO  BE      
                     TRANSFORMED;     
             EXIT:   THE N2 - N1 + 1 VECTORS OF LENGTH N RESULTING  FROM      
                     THE BACK TRANSFORMATION;   
        
 PROCEDURES USED: 
     TAMVEC          =      CP34012,  
     ICHROW          =      CP34032.  
        
 REQUIRED CENTRAL MEMORY:   
     EXECUTION FIELD LENGTH:
             IONAL REAL ARRAY OF LENGTH N IS DECLARED.    
        
 RUNNING TIME: ROUGHLY PROPORTIONAL TO (N2 - N1 + 1) * N * N.       
        
 LANGUAGE:   ALGOL 60.      
        
 METHOD AND PERFORMANCE:    
     SEE SUBSECTION BAKREAHES1.       
        
 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]      J. H. WILKINSON, THE ALGEBRAIC EIGENVALUE PROBLEM,    
              CLARENDON PRESS, OXFORD, 1965.    
1SECTION:3.2.1.2.1.2          (JUNE 1974)                         PAGE 5      
        
        
        
 EXAMPLES OF USE: 
        
     EXAMPLES  OF  USE OF  TFMREAHES, BAKREAHES1  AND  BAKREAHES2 CAN BE      
     FOUND  IN  THE  PROCEDURES  FOR   CALCULATING    EIGENVALUES    AND      
     EIGENVECTORS AS DESCRIBED IN SECTION 3.3.1.2.2.      
        
 SOURCE TEXT(S) : 
        
0"CODE" 34170;    
     "COMMENT" MCA 2400;    
     "PROCEDURE" TFMREAHES(A, N, EM, INT); "VALUE" N; "INTEGER" N;  
      "ARRAY" A, EM; "INTEGER" "ARRAY" INT;     
     "BEGIN" "INTEGER" I, J, J1, K, L;
         "REAL" S, T, MACHTOL, MACHEPS, NORM;   
          "ARRAY" B[0:N - 1];         
        
         "REAL" "PROCEDURE" MATVEC(L, U, I, A, B); "CODE" 34011;    
         "REAL" "PROCEDURE" MATMAT(L, U, I, J, A, B); "CODE" 34013; 
         "PROCEDURE" ICHCOL(L, U, I, J, A); "CODE" 34031; 
         "PROCEDURE" ICHROW(L, U, I, J, A); "CODE" 34032; 
        
         MACHEPS:= EM[0]; NORM:= 0;   
         "FOR" I:= 1 "STEP" 1 "UNTIL" N "DO"    
         "BEGIN" S:= 0;     
             "FOR" J:= 1 "STEP" 1 "UNTIL" N "DO" S:= S + ABS(A[I,J]);         
             "IF" S > NORM "THEN" NORM:= S      
         "END";   
         EM[1]:= NORM; MACHTOL:= NORM * MACHEPS; INT[1]:= 0;        
         "FOR" J:= 2 "STEP" 1 "UNTIL" N "DO"    
         "BEGIN" J1:= J - 1; L:= 0; S:= MACHTOL;
             "FOR" K:= J + 1 "STEP" 1 "UNTIL" N "DO"      
             "BEGIN" T:= ABS(A[K,J1]); "IF" T > S "THEN"  
                 "BEGIN" L:= K; S:= T "END"     
             "END";         
             "IF" L ^= 0 "THEN"       
             "BEGIN" "IF" ABS(A[J,J1]) < S "THEN"         
                 "BEGIN" ICHROW(1, N, J, L, A); 
                     ICHCOL(1, N, J, L, A)      
                 "END"      
                 "ELSE" L:= J; T:= A[J,J1];     
                 "FOR" K:= J + 1 "STEP" 1 "UNTIL" N "DO"  
                 A[K,J1]:= A[K,J1] / T
             "END"
             "ELSE"         
             "FOR" K:= J + 1 "STEP" 1 "UNTIL" N "DO" A[K,J1]:= 0;   
             "FOR" I:= 1 "STEP" 1 "UNTIL" N "DO"
             B[I - 1]:= A[I,J]:= A[I,J] +       
             ("IF" L = 0 "THEN" 0 "ELSE" MATMAT(J + 1, N, I, J1, A, A))-      
             MATVEC(1, "IF" J1 < I - 2 "THEN" J1 "ELSE" I - 2, I, A, B);      
             INT[J]:= L     
         "END"    
     "END" TFMREAHES;       
         "EOP"    
1SECTION:3.2.1.2.1.2          (JUNE 1974)                         PAGE 6      
        
        
        
0"CODE" 34171;    
     "COMMENT" MCA 2401;    
     "PROCEDURE" BAKREAHES1(A, N, INT, V); "VALUE" N; "INTEGER" N;  
      "ARRAY" A, V; "INTEGER" "ARRAY" INT;      
     "BEGIN" "INTEGER" I, L;
         "REAL" W; "ARRAY" X[1:N];    
        
         "REAL" "PROCEDURE" MATVEC(L, U, I, A, B); "CODE" 34011;    
        
         "FOR" I:= 2 "STEP" 1 "UNTIL" N "DO" X[I - 1]:= V[I];       
         "FOR" I:= N "STEP" -1 "UNTIL" 2 "DO"   
         "BEGIN" V[I]:= V[I] + MATVEC(1, I - 2, I, A, X); 
             L:= INT[I]; "IF" L > I "THEN"      
             "BEGIN" W:= V[I]; V[I]:= V[L]; V[L]:= W "END"
         "END"    
     "END" BAKREAHES1;      
         "EOP"    
        
0"CODE" 34172;    
     "COMMENT" MCA 2402;    
     "PROCEDURE" BAKREAHES2(A, N, N1, N2, INT, VEC); "VALUE" N, N1, N2;       
     "INTEGER" N, N1, N2; "ARRAY" A, VEC; "INTEGER" "ARRAY" INT;    
     "BEGIN" "INTEGER" I, L, K; "ARRAY" U[1:N]; 
        
         "REAL" "PROCEDURE" TAMVEC(L, U, I, A, B); "CODE" 34012;    
         "PROCEDURE" ICHROW(L, U, I, J, A); "CODE" 34032; 
        
         "FOR" I:= N "STEP" -1 "UNTIL" 2 "DO"   
         "BEGIN" "FOR" K:= I - 2 "STEP" -1 "UNTIL" 1 "DO" 
             U[K + 1]:= A[I,K];       
             "FOR" K:= N1 "STEP" 1 "UNTIL" N2 "DO"        
             VEC[I,K]:= VEC[I,K] + TAMVEC(2 , I - 1, K, VEC, U);    
             L:= INT[I]; "IF" L > I "THEN" ICHROW(N1, N2, I, L, VEC)
         "END"    
     "END" BAKREAHES2;      
         "EOP"    
1SECTION:3.2.1.2.2.1          (JUNE 1974)                         PAGE 1      
        
        
        
 AUTHOR   : C.G. VAN DER LAAN.        
        
        
 CONTRIBUTORS : H.FIOLET, C.G. VAN DER LAAN.    
        
        
 INSTITUTE : MATHEMATICAL CENTRE.     
        
        
 RECEIVED : 730903.         
        
        
 BRIEF DESCRIPTION :        
        
     THIS SECTION CONTAINS THREE PROCEDURES;    
     A) HSHHRMTRI TRANSFORMS THE HERMITIAN MATRIX M INTO A SIMILAR  
         REAL SYMMETRIC TRIDIAGONAL MATRIX S;   
     B) BAKHRMTRI PERFORMS A BACK TRANSFORMATION CORRESPONDING TO   
     HSHHRMTRI;   
     C) HSHHRMTRIVAL DELIVERS THE MAIN DIAGONAL ELEMENTS AND THE    
     SQUARES OF THE CODIAGONAL ELEMENTS OF A HERMITIAN TRIDIAGONAL  
     MATRIX WHICH IS UNITARY SIMILAR WITH A GIVEN HERMITIAN MATRIX. 
        
        
 KEYWORDS :       
        
     HERMITIAN MATRIX ,     
     TRIDIAGONALIZATION ,   
     COMPLEX HOUSEHOLDER,S TRANSFORMATION .     
        
        
 SUBSECTION : HSHHRMTRI.    
        
        
 CALLING SEQUENCE :         
        
     THE HEADING OF THE PROCEDURE READS :       
     "PROCEDURE" HSHHRMTRI(A, N, D, B, BB, EM, TR, TI); "VALUE" N;  
     "INTEGER" N; "ARRAY" A, D, B, BB, EM, TR, TI;        
        
     THE MEANING OF THE FORMAL PARAMETERS IS :  
     A,TR,TI:   <ARRAY IDENTIFIER>;   
                "ARRAY"A[1:N,1:N];    
                "ARRAY" TR,TI[1:N-1]; 
                ENTRY: THE  REAL  PART  OF  THE  UPPER  TRIANGLE OF  THE      
                       HERMITIAN   MATRIX  MUST  BE  GIVEN IN THE  UPPER      
                       TRIANGULAR PART OF A (THE ELEMENTS A[I,J], I<=J);      
                       THE  IMAGINARY  PART OF THE STRICT LOWER TRIANGLE      
                       OF  THE   HERMITIAN  MATRIX  MUST BE GIVEN IN THE      
                       STRICT LOWER PART OF A (THE ELEMENTS A[I,J],I>J);      
                EXIT: DATA FOR THE BACKTRANSFORMATION;    
1SECTION:3.2.1.2.2.1          (JUNE 1974)                         PAGE 2      
        
        
        
     N:         <ARITHMETIC EXPRESSION>;        
                THE ORDER OF THE MATRIX;        
     D:         <ARRAY IDENTIFIER>;   
                "ARRAY"D[1:N];        
                EXIT : THE  MAIN  DIAGONAL  OF THE  RESULTING  SYMMETRIC      
                       TRIDIAGONAL MATRIX;      
     B:         <ARRAY IDENTIFIER>;   
                "ARRAY"B[1:N-1];      
                EXIT: THE CODIAGONAL ELEMENTS OF THE RESULTING SYMMETRIC      
                      TRIDIAGONAL MATRIX;       
     BB:        <ARRAY IDENTIFIER>;   
                "ARRAY"BB[1:N-1];     
                EXIT : THE  SQUARES  OF  THE  MODULI  OF  THE CODIAGONAL      
                       ELEMENTS  OF THE RESULTING  SYMMETRIC TRIDIAGONAL      
                       MATRIX;        
     EM:        <ARRAY IDENTIFIER>;   
                "ARRAY"EM[0:1];       
                ENTRY: EM[0], THE MACHINE PRECISION;      
                EXIT:  EM[1], AN ESTIMATE FOR A NORM OF THE ORIGINAL
                              MATRIX. 
        
        
 PROCEDURES USED :
        
     MATVEC    = CP34011 ,  
     TAMVEC    = CP34012 ,  
     MATMAT    = CP34013 ,  
     TAMMAT    = CP34014 ,  
     MATTAM    = CP34015 ,  
     ELMVECCOL = CP34021 ,  
     ELMCOLVEC = CP34022 ,  
     ELMCOL    = CP34023 ,  
     ELMROW    = CP34024 ,  
     ELMVECROW = CP34026 ,  
     ELMROWVEC = CP34027 ,  
     ELMROWCOL = CP34028 ,  
     ELMCOLROW = CP34029 ,  
     CARPOL    = CP34344 .  
        
        
 RUNNING TIME : PROPORTIONAL TO N CUBED .       
        
        
 LANGUAGE: ALGOL 60.        
        
        
 METHOD AND PERFORMANCE: SEE HSHHRMTRIVAL (THIS SECTION). 
        
        
1SECTION:3.2.1.2.2.1          (JUNE 1974)                         PAGE 3      
        
        
        
 SUBSECTION : BAKHRMTRI.    
        
        
 CALLING SEQUENCE :         
        
     THE HEADING OF THE PROCEDURE READS :       
     "PROCEDURE" BAKHRMTRI(A, N, N1, N2, VECR, VECI, TR, TI);       
     "VALUE" N, N1, N2; "INTEGER" N, N1, N2;    
        
     THE MEANING OF THE FORMAL PARAMETERS IS :  
     A,TR,TI:   <ARRAY IDENTIFIER>;   
                "ARRAY"A[1:N,1:N];    
                "ARRAY" TR,TI[1:N-1]; 
                ENTRY: THE DATA FOR THE BACKTRANSFORMATION AS PRODUCED        
                       BY HSHHRMTRI;  
     N:         <ARITHMETIC EXPRESSION>;        
                THE ORDER OF THE MATRIX OF WHICH THE EIGENVECTORS ARE         
                CALCULATED; 
     N1,N2:     <ARITHMETIC EXPRESSION>;        
                THE EIGENVECTORS CORRESPONDING TO THE  EIGENVALUES  WITH      
                INDICES N1,...,N2 ARE TO BE TRANSFORMED;  
     VECR,VECI: <ARRAY IDENTIFIER>;   
                "ARRAY" VECR,VECI[1:N,N1:N2];   
                ENTRY:      
                THE BACK TRANSFORMATION IS PERFORMED ON THE REAL    
                EIGENVECTORS GIVEN IN THE COLUMNS OF ARRAY VECR;    
                EXIT:       
                VECR : REAL PART OF THE TRANSFORMED EIGENVECTORS;   
                VECI : IMAGINARY PART OF THE TRANSFORMED EIGENVECTORS.        
        
        
 PROCEDURES USED :
        
     MATMAT    = CP34013 ,  
     TAMMAT    = CP34014 ,  
     ELMCOL    = CP34023 ,  
     ELMCOLROW = CP34029 ,  
     COMMUL    = CP34341 ,  
     COMROWCST = CP34353 .  
        
 RUNNING TIME: PROPORTIONAL TO  (N2-N1+1)*N**2 .
        
        
 LANGUAGE: ALGOL 60.        
        
        
 METHOD AND PERFORMANCE: SEE HSHHRMTRIVAL (THIS SECTION). 
        
        
1SECTION:3.2.1.2.2.1          (JUNE 1974)                         PAGE 4      
        
        
        
 SUBSECTION : HSHHRMTRIVAL. 
        
        
 CALLING SEQUENCE :         
        
     THE HEADING OF THE PROCEDURE READS :       
     "PROCEDURE" HSHHRMTRIVAL(A, N, D, BB, EM); "VALUE" N; "INTEGER" N;       
      "ARRAY" A, D, BB, EM; 
        
     THE MEANING OF THE FORMAL PARAMETERS IS :  
     A:         <ARRAY IDENTIFIER>;   
                "ARRAY"A[1:N,1:N];    
                ENTRY: THE  REAL  PART  OF  THE  UPPER  TRIANGLE OF  THE      
                       HERMITIAN   MATRIX  MUST  BE  GIVEN IN THE  UPPER      
                       TRIANGULAR PART OF A (THE ELEMENTS A[I,J], I<=J);      
                       THE  IMAGINARY  PART OF THE STRICT LOWER TRIANGLE      
                       OF  THE   HERMITIAN  MATRIX  MUST BE GIVEN IN THE      
                       STRICT LOWER PART OF A (THE ELEMENTS A[I,J],I>J);      
                       THE ELEMENTS OF A ARE ALTERED;     
     N:         <ARITHMETIC EXPRESSION>;        
                THE ORDER OF THE GIVEN MATRIX;  
     D:         <ARRAY IDENTIFIER>;   
                "ARRAY"D[1:N];        
                EXIT: THE  MAIN  DIAGONAL  OF THE  RESULTING HERMITIAN        
                      TRIDIAGONAL MATRIX;       
     BB:        <ARRAY IDENTIFIER>;   
                "ARRAY"BB[1:N-1];     
                EXIT: THE   SQUARES   OF  THE  MODULI OF THE  CODIAGONAL      
                      ELEMENTS OF THE RESULTING HERMITIAN  TRIDIAGONAL        
                      MATRIX;         
     EM:        <ARRAY IDENTIFIER>;   
                "ARRAY"EM[0:1];       
                ENTRY: EM[0],THE MACHINE PRECISION;       
                EXIT:  EM[1], AN ESTIMATE FOR A NORM OF THE ORIGINAL
                              MATRIX; 
        
        
 PROCEDURES USED :
        
     MATVEC    = CP34011 ,  
     TAMVEC    = CP34012 ,  
     MATMAT    = CP34013 ,  
     TAMMAT    = CP34014 ,  
     MATTAM    = CP34015 ,  
     ELMVECCOL = CP34021 ,  
     ELMCOLVEC = CP34022 ,  
     ELMCOL    = CP34023 ,  
     ELMROW    = CP34024 ,  
     ELMVECROW = CP34026 ,  
     ELMROWVEC = CP34027 ,  
     ELMROWCOL = CP34028 ,  
     ELMCOLROW = CP34029 .  
        
1SECTION:3.2.1.2.2.1          (JUNE 1974)                         PAGE 5      
        
        
        
        
 RUNNING TIME : PROPORTIONAL TO N CUBED .       
        
        
 LANGUAGE : ALGOL 60.       
        
        
 THE FOLLOWING HOLDS FOR THE THREE PROCEDURES : 
        
        
 METHOD AND PERFORMANCE :   
        
     HSHHRMTRIVAL TRANSFORMS A HERMITIAN MATRIX INTO A SIMILAR      
     HERMITIAN TRIDIAGONAL MATRIX BY MEANS OF HOUSEHOLDER'S         
     TRANSFORMATION.        
     HSHHRMTRI TRANSFORMS A HERMITIAN MATRIX INTO A SIMILAR         
     REAL  TRIDIAGONAL MATRIX BY MEANS OF HOUSEHOLDER'S   
     TRANSFORMATION FOLLOWED BY A COMPLEX DIAGONAL UNITARY
     SIMILARITY TRANSFORMATION IN ORDER TO MAKE THE RESULTING       
     TRIDIAGONAL MATRIX REAL SYMMETRIC;         
     HOUSEHOLDER'S TRANSFORMATION FOR COMPLEX HERMITIAN MATRICES    
     IS A UNITARY SIMILARITY TRANSFORMATION,TRANSFORMING A HERMITIAN
     MATRIX INTO A SIMILAR COMPLEX TRIDIAGONAL ONE (SEE WILKINSON,  
     1965, P. 342-343). LET M BE A GIVEN HERMITIAN MATRIX OF ORDER  
     N, WITH REAL PART MR AND IMAGINARY PART MI, P THE    
     TRANSFORMING MATRIX AND T THE RESULTING HERMITIAN TRIDIAGONAL  
     MATRIX. SINCE P IS UNITARY, WE HAVE T = P"MP, WHERE  
     " STANDS  FOR CONJUGATING AND TRANSPOSING. THE MATRIX  P IS THE
     PRODUCT OF N-2 HOUSEHOLDER MATRICES, THESE BEING UNITARY       
     HERMITIAN MATRICES OF THE FORM I-U"U/T, WHERE T IS A SCALAR    
     (>0), AND U A COMPLEX VECTOR. THE K-TH HOUSEHOLDER MATRIX,     
     K=1,...,N-2, IS CHOSEN IN SUCH A WAY THAT THE LAST K ELEMENTS OF U       
     VANISH, AND THE DESIRED ZEROS ARE INTRODUCED IN THE (N-K+1)-TH 
     COLUMN AND ROW OF THE MATRIX M. HOWEVER, IF THE EUCLIDIAN NORM 
     OF THE FIRST N-K-1 ELEMENTS OF COLUMN N-K+1 OF THE MATRIX M IS 
     SMALLER THAN THE MACHINE PRECISION TIMES THE INFINITY NORM OF THE        
     MATRIX ( NORM(MR) + NORM(MI) ), THEN THE K-TH TRANSFORMATION   
     IS SKIPPED (I.E. THE K-TH HOUSEHOLDER MATRIX IS REPLACED BY I).
1SECTION:3.2.1.2.2.1          (JUNE 1974)                         PAGE 6      
        
        
        
     THE COMPLEX DIAGONAL SIMILARITY TRANSFORMATION D TRANSFORMS THE
     HERMITIAN TRIDIAGONAL MATRIX T INTO A REAL SYMMETRIC TRIDIAGONAL         
     MATRIX S (MUELLER, 1966). THE DIAGONAL OF D IS CHOSEN IN SUCH  
     A WAY THAT THE CODIAGONAL ELEMENTS OF T ARE TRANSFORMED INTO   
     THEIR ABSOLUTE VALUES. 
     BAKHRMTRI PERFORMS THE BACK TRANSFORMATION TO REPLACE THE      
     EIGENVECTORS OF THE TRIDIAGONAL SYMMETRIC MATRIX S BY THE      
     EIGENVECTORS OF THE ORIGINAL HERMITIAN MATRIX M. IF X IS AN    
     EIGENVECTOR OF S THEN PDX IS THE CORRESPONDING EIGENVECTOR OF  
     M. STARTING FROM THE VECTOR V=DX, THE VECTOR PDX IS  
     OBTAINED BY SUCCESSIVELY REPLACING V BY THE K-TH HOUSEHOLDER   
     MATRIX TIMES V, FOR K=N-2,...,1. THE RESULTING VECTOR V THEN EQUALS      
     PDX.         
        
        
 REFERENCES :     
        
     MUELLER, D.J. (1966),  
     HOUSEHOLDER,S METHOD  FOR  COMPLEX MATRICES AND EIGENSYSTEMS   OF        
     HERMITIAN MATRICES,    
     NUMER.MATH., 8, P.72-92;         
        
     WILKINSON, J.H. (1965),
     THE ALGEBRAIC EIGENVALUE PROBLEM,
     CLARENDON PRESS, OXFORD.         
        
        
 EXAMPLE OF USE : 
        
     THE  PROCEDURES HSHHRMTRIVAL  AND  BAKHRMTRI  ARE  USED  IN SECTION      
     3.3.2.1. :   
     EIGVALHRM AND QRIVALHRM USE HSHHRMTRIVAL,  
     EIGHRM AND QRIHRM USE BAKHRMTRI .
     AS A FORMAL  TEST OF THE PROCEDURE HSHHRMTRI, THE FOLLOWING  MATRIX      
     WAS USED (SEE GREGORY AND KARNEY, CHAPTER 6, EXAMPLE 6.6) :    
          3    1    0   +2I 
          1    3   -2I   0  
          0   +2I   1    1  
         -2I   0    1    1  
        
1SECTION:3.2.1.2.2.1          (JUNE 1974)                         PAGE 7      
        
        
        
     "BEGIN"      
     "COMMENT" GREGORY AND KARNEY,CHAPTER 6, EXAMPLE 6.6; 
     "PROCEDURE" HSHHRMTRI(A,N,D,B,BB,EM,TR,TI);"CODE" 34363;       
     "PROCEDURE" INIMAT(LR,UR,LC,UC,A,X);"CODE" 31011;    
     "REAL" "ARRAY" A[1:4,1:4],D,B,BB[1:4],TR,TI[1:3],EM[0:1];      
     "INTEGER" I,J;         
     "PROCEDURE" OUT(S,A,N);
     "VALUE" N;"INTEGER" N;"ARRAY" A;"STRING" S;
     "BEGIN" "INTEGER" I,J; 
         OUTPUT(61,"("10S")",S);      
         "FOR" I:=1 "STEP" 1 "UNTIL" N "DO"     
         OUTPUT(61,"("+D.3DBB")",A[I]);         
         OUTPUT(61,"("/")") 
     "END" OUT;   
        
     INIMAT(1,4,1,4,A,0);   
     A[1,1]:=A[2,2]:=3;     
     A[1,2]:=A[3,3]:=A[3,4]:=A[4,4]:=1;         
     A[3,2]:=2;A[4,1]:=-2;  
     EM[0]:="-14; 
     OUTPUT(61,"(""("INITIAL MATRIX GIVEN IN ARRAY A[1:4,1:4]:")",/")");      
     "FOR" I:=1 "STEP" 1 "UNTIL" 4 "DO"         
     "BEGIN" "FOR" J:=1 "STEP" 1 "UNTIL" 4 "DO" 
         OUTPUT(61,"("-DBBB")",A[I,J]);         
         OUTPUT(61,"("/")") 
     "END";       
     OUTPUT(61,"("/,"("HSHHRMTRI DELIVERS:")",//")");     
     HSHHRMTRI(A,4,D,B,BB,EM,TR,TI);  
     OUT("("D[1:4]: ")",D,4);         
     OUT("("B[1:3]: ")",B,3);         
     OUT("("BB[1:3]: ")",BB,3);       
     OUT("("EM[1]: ")",EM,1);         
     "END"        
        
     OUTPUT :     
     INITIAL MATRIX GIVEN IN ARRAY A[1:4,1:4]:  
      3    1    0    0      
      0    3    0    0      
      0    2    1    1      
     -2    0    0    1      
        
     HSHHRMTRI DELIVERS:    
        
     D[1:4]:   +3.000  +1.400  +2.600  +1.000   
     B[1:3]:   +2.236  +0.800  +2.236 
     BB[1:3]:  +5.000  +0.640  +5.000 
     EM[1]:    +6.000       
1SECTION:3.2.1.2.2.1          (JUNE 1974)                         PAGE 8      
        
        
        
 SOURCE TEXT(S) : 
        
0"CODE" 34363;    
     "PROCEDURE" HSHHRMTRI(A, N, D, B, BB, EM, TR, TI); "VALUE" N;  
     "INTEGER" N; "ARRAY" A, D, B, BB, EM, TR, TI;        
     "BEGIN" "INTEGER" I, J, J1, JM1, R, RM1;   
         "REAL" NRM, W, TOL2, X, AR, AI, MOD, C, S, H, K, T, Q,     
         AJR, ARJ, BJ, BBJ; 
         "REAL" "PROCEDURE" MATVEC(L,U,I,A,B);"CODE" 34011;         
         "REAL" "PROCEDURE" TAMVEC(L,U,I,A,B);"CODE" 34012;         
         "REAL" "PROCEDURE" MATMAT(L,U,I,J,A,B);"CODE" 34013;       
         "REAL" "PROCEDURE" TAMMAT(L,U,I,J,A,B);"CODE" 34014;       
         "REAL" "PROCEDURE" MATTAM(L,U,I,J,A,B);"CODE" 34015;       
         "PROCEDURE" ELMVECCOL(L,U,I,A,B,X);"CODE" 34021; 
         "PROCEDURE" ELMCOLVEC(L,U,I,A,B,X);"CODE" 34022; 
         "PROCEDURE" ELMCOL(L,U,I,J,A,B,X);"CODE" 34023;  
         "PROCEDURE" ELMROW(L,U,I,J,A,B,X);"CODE" 34024;  
         "PROCEDURE" ELMVECROW(L,U,I,A,B,X);"CODE" 34026; 
         "PROCEDURE" ELMROWVEC(L,U,I,A,B,X);"CODE" 34027; 
         "PROCEDURE" ELMROWCOL(L,U,I,J,A,B,X);"CODE" 34028;         
         "PROCEDURE" ELMCOLROW(L,U,I,J,A,B,X);"CODE" 34029;         
         "PROCEDURE" CARPOL(AR,AI,R,C,S);"CODE" 34344;    
         NRM:= 0; 
         "FOR" I:= 1 "STEP" 1 "UNTIL" N "DO"    
         "BEGIN" W:= ABS(A[I,I]);     
             "FOR" J:= I - 1 "STEP" - 1 "UNTIL" 1, I + 1 "STEP" 1   
             "UNTIL" N "DO" W:= W + ABS(A[I,J]) + ABS(A[J,I]);      
             "IF" W > NRM "THEN" NRM:= W        
         "END" I; 
         TOL2:= (EM[0] * NRM) ** 2; EM[1]:= NRM; R:= N;   
         "FOR" RM1:= N - 1 "STEP" - 1 "UNTIL" 1 "DO"      
         "BEGIN" X:= TAMMAT(1, R - 2, R, R, A, A) + MATTAM(1, R -   
             2, R, R, A, A); AR:= A[RM1,R]; AI:= - A[R,RM1];        
             D[R]:= A[R,R]; CARPOL(AR, AI, MOD, C, S);    
             "IF" X < TOL2 "THEN"     
             "BEGIN" A[R,R]:= - 1; B[RM1]:= MOD;
                 BB[RM1]:= MOD * MOD  
             "END"
        
1SECTION:3.2.1.2.2.1          (JUNE 1974)                         PAGE 9      
        
        
        
             "ELSE"         
             "BEGIN" H:= MOD * MOD + X; K:= SQRT(H);      
                 T:= A[R,R]:= H + MOD * K;      
                 "IF" AR = 0 "AND" AI = 0 "THEN" A[RM1,R]:= K "ELSE"
                 "BEGIN" A[RM1,R]:= AR + C * K; 
                     A[R,RM1]:= - AI - S * K; S:= - S     
                 "END";     
                 C:= - C; J:= 1; JM1:= 0;       
                 "FOR" J1:= 2 "STEP" 1 "UNTIL" R "DO"     
                 "BEGIN" B[J]:= (TAMMAT(1, J, J, R, A, A) +         
                     MATMAT(J1, RM1, J, R, A, A) + MATTAM(1,        
                     JM1, J, R, A, A) - MATMAT(J1, RM1, R, J,       
                     A, A)) / T;      
                     BB[J]:= (MATMAT(1, JM1, J, R, A, A) -
                     TAMMAT(J1, RM1, J, R, A, A) - MATMAT(1, J,     
                     R, J, A, A) - MATTAM(J1, RM1, J, R, A, A))     
                     / T; JM1:= J; J:= J1       
                 "END" J1;  
                 Q:= (TAMVEC(1, RM1, R, A, B) - MATVEC(1, RM1,      
                 R, A, BB)) / T / 2;  
                 ELMVECCOL(1, RM1, R, B, A, - Q);         
                 ELMVECROW(1, RM1, R, BB, A, Q); J:= 1;   
                 "FOR" J1:= 2 "STEP" 1 "UNTIL" R "DO"     
                 "BEGIN" AJR:= A[J,R]; ARJ:= A[R,J]; BJ:= B[J];     
                     BBJ:= BB[J];     
                     ELMROWVEC(J, RM1, J, A, B, - AJR);   
                     ELMROWVEC(J, RM1, J, A, BB, ARJ);    
                     ELMROWCOL(J, RM1, J, R, A, A, - BJ); 
                     ELMROW(J, RM1, J, R, A, A, BBJ);     
                     ELMCOLVEC(J1, RM1, J, A, B, - ARJ);  
                     ELMCOLVEC(J1, RM1, J, A, BB, - AJR); 
                     ELMCOL(J1, RM1, J, R, A, A, BBJ);    
                     ELMCOLROW(J1, RM1, J, R, A, A, BJ); J:= J1;    
                 "END" J1;  
                 BB[RM1]:= H; B[RM1]:= K;       
             "END";         
             TR[RM1]:= C; TI[RM1]:= S; R:= RM1; 
         "END" RM1;         
         D[1]:= A[1,1];     
     "END" HSHHRMTRI;       
         "EOP"    
        
1SECTION:3.2.1.2.2.1          (JUNE 1974)                        PAGE 10      
        
        
0"CODE" 34365;    
     "PROCEDURE" BAKHRMTRI(A, N, N1, N2, VECR, VECI, TR, TI);       
     "VALUE" N, N1, N2; "INTEGER" N, N1, N2;    
      "ARRAY" A, VECR, VECI, TR, TI;  
     "BEGIN" "INTEGER" I, J, R, RM1;  
         "REAL" C, S, T, QR, QI;      
         "REAL" "PROCEDURE" MATMAT(L,U,I,J,A,B);"CODE" 34013;       
         "REAL" "PROCEDURE" TAMMAT(L,U,I,J,A,B);"CODE" 34014;       
         "PROCEDURE" ELMCOL(L,U,I,J,A,B,X);"CODE" 34023;  
         "PROCEDURE" ELMCOLROW(L,U,I,J,A,B,X);"CODE" 34029;         
         "PROCEDURE" COMMUL(AR,AI,BR,BI,RR,RI);"CODE" 34341;        
         "PROCEDURE" COMROWCST(L,U,I,AR,AI,XR,XI);"CODE" 34353;     
         "FOR" I:= 1 "STEP" 1 "UNTIL" N "DO"    
         "FOR" J:= N1 "STEP" 1 "UNTIL" N2 "DO" VECI[I,J]:= 0; C:= 1;
         S:= 0;   
         "FOR" J:= N - 1 "STEP" - 1 "UNTIL" 1 "DO"        
         "BEGIN" COMMUL(C, S, TR[J], TI[J], C, S);        
             COMROWCST(N1, N2, J, VECR, VECI, C, S)       
         "END" J; 
         RM1:= 2; 
         "FOR" R:= 3 "STEP" 1 "UNTIL" N "DO"    
         "BEGIN" T:= A[R,R]; "IF" T > 0 "THEN"  
             "FOR" J:= N1 "STEP" 1 "UNTIL" N2 "DO"        
             "BEGIN" QR:= (TAMMAT(1, RM1, R, J, A, VECR) -
                 MATMAT(1, RM1, R, J, A, VECI)) / T;      
                 QI:= (TAMMAT(1, RM1, R, J, A, VECI) +    
                 MATMAT(1, RM1, R, J, A, VECR)) / T;      
                 ELMCOL(1, RM1, J, R, VECR, A, - QR);     
                 ELMCOLROW(1, RM1, J, R, VECR, A, - QI);  
                 ELMCOLROW(1, RM1, J, R, VECI, A, QR);    
                 ELMCOL(1, RM1, J, R, VECI, A, - QI)      
             "END";         
             RM1:= R;       
         "END" R  
     "END" BAKHRMTRI;       
         "EOP"    
        
1SECTION:3.2.1.2.2.1          (JUNE 1974)                        PAGE 11      
        
        
0"CODE" 34364;    
     "PROCEDURE" HSHHRMTRIVAL(A, N, D, BB, EM); "VALUE" N; "INTEGER" N;       
      "ARRAY" A, D, BB, EM; 
     "BEGIN" "INTEGER" I, J, J1, JM1, R, RM1;   
         "REAL" NRM, W, TOL2, X, AR, AI, H, T, Q, AJR, ARJ, DJ,     
         BBJ, MOD2;         
         "REAL" "PROCEDURE" MATVEC(L,U,I,A,B);"CODE" 34011;         
         "REAL" "PROCEDURE" TAMVEC(L,U,I,A,B);"CODE" 34012;         
         "REAL" "PROCEDURE" MATMAT(L,U,I,J,A,B);"CODE" 34013;       
         "REAL" "PROCEDURE" TAMMAT(L,U,I,J,A,B);"CODE" 34014;       
         "REAL" "PROCEDURE" MATTAM(L,U,I,J,A,B);"CODE" 34015;       
         "PROCEDURE" ELMVECCOL(L,U,I,A,B,X);"CODE" 34021; 
         "PROCEDURE" ELMCOLVEC(L,U,I,A,B,X);"CODE" 34022; 
         "PROCEDURE" ELMCOL(L,U,I,J,A,B,X);"CODE" 34023;  
         "PROCEDURE" ELMROW(L,U,I,J,A,B,X);"CODE" 34024;  
         "PROCEDURE" ELMVECROW(L,U,I,A,B,X);"CODE" 34026; 
         "PROCEDURE" ELMROWVEC(L,U,I,A,B,X);"CODE" 34027; 
         "PROCEDURE" ELMROWCOL(L,U,I,J,A,B,X);"CODE" 34028;         
         "PROCEDURE" ELMCOLROW(L,U,I,J,A,B,X);"CODE" 34029;         
         NRM:= 0; 
         "FOR" I:= 1 "STEP" 1 "UNTIL" N "DO"    
         "BEGIN" W:= ABS(A[I,I]);     
             "FOR" J:= I - 1 "STEP" - 1 "UNTIL" 1, I + 1 "STEP" 1   
             "UNTIL" N "DO" W:= W + ABS(A[I,J]) + ABS(A[J,I]);      
             "IF" W > NRM "THEN" NRM:= W        
         "END" I; 
         TOL2:= (EM[0] * NRM) ** 2; EM[1]:= NRM; R:= N;   
         "FOR" RM1:= N - 1 "STEP" - 1 "UNTIL" 1 "DO"      
         "BEGIN" X:= TAMMAT(1, R - 2, R, R, A, A) + MATTAM(1, R -   
             2, R, R, A, A); AR:= A[RM1,R]; AI:= - A[R,RM1];        
             D[R]:= A[R,R]; 
             "IF" X < TOL2 "THEN" BB[RM1]:= AR * AR + AI * AI "ELSE"
             "BEGIN" MOD2:= AR * AR + AI * AI; "IF" MOD2 = 0 "THEN" 
                 "BEGIN" A[RM1,R]:= SQRT(X); T:= X "END"  
                 "ELSE"     
                 "BEGIN" X:= X + MOD2; H:= SQRT(MOD2 * X);
                     T:= X + H; H:= 1 + X / H;  
                     A[R,RM1]:= - AI * H; A[RM1,R]:= AR * H;        
                 "END";     
                                                               "COMMENT"      
1SECTION:3.2.1.2.2.1          (JUNE 1974)                        PAGE 12      
                                                                 ;  
        
        
                 J:= 1; JM1:= 0;      
                 "FOR" J1:= 2 "STEP" 1 "UNTIL" R "DO"     
                 "BEGIN" D[J]:= (TAMMAT(1, J, J, R, A, A) +         
                     MATMAT(J1, RM1, J, R, A, A) + MATTAM(1,        
                     JM1, J, R, A, A) - MATMAT(J1, RM1, R, J,       
                     A, A)) / T;      
                     BB[J]:= (MATMAT(1, JM1, J, R, A, A) -
                     TAMMAT(J1, RM1, J, R, A, A) - MATMAT(1, J,     
                     R, J, A, A) - MATTAM(J1, RM1, J, R, A, A))     
                     / T; JM1:= J; J:= J1       
                 "END" J1;  
                 Q:= (TAMVEC(1, RM1, R, A, D) - MATVEC(1, RM1,      
                 R, A, BB)) / T / 2;  
                 ELMVECCOL(1, RM1, R, D, A, - Q);         
                 ELMVECROW(1, RM1, R, BB, A, Q); J:= 1;   
                 "FOR" J1:= 2 "STEP" 1 "UNTIL" R "DO"     
        
                 "BEGIN" AJR:= A[J,R]; ARJ:= A[R,J]; DJ:= D[J];     
                     BBJ:= BB[J];     
                     ELMROWVEC(J, RM1, J, A, D, - AJR);   
                     ELMROWVEC(J, RM1, J, A, BB, ARJ);    
                     ELMROWCOL(J, RM1, J, R, A, A, - DJ); 
                     ELMROW(J, RM1, J, R, A, A, BBJ);     
                     ELMCOLVEC(J1, RM1, J, A, D, - ARJ);  
                     ELMCOLVEC(J1, RM1, J, A, BB, - AJR); 
                     ELMCOL(J1, RM1, J, R, A, A, BBJ);    
                     ELMCOLROW(J1, RM1, J, R, A, A, DJ); J:= J1;    
                 "END" J1;  
                 BB[RM1]:= X;         
             "END";         
             R:= RM1;       
         "END" RM1;         
         D[1]:= A[1,1];     
     "END" HSHHRMTRIVAL;    
         "EOP"    
1SECTION:3.2.1.2.2.2          (JUNE 1974)                         PAGE 1      
        
        
        
 AUTHOR   : C.G. VAN DER LAAN.        
        
        
 CONTRIBUTORS : H.FIOLET, C.G. VAN DER LAAN.    
        
        
 INSTITUTE : MATHEMATICAL CENTRE.     
        
        
 RECEIVED: 731016.
        
        
 BRIEF DESCRIPTION :        
        
     THIS SECTION CONTAINS THE PROCEDURES HSHCOMHES AND BAKCOMHES.  
     HSHCOMHES TRANSFORMS A COMPLEX MATRIX BY MEANS OF HOUSEHOLDER'S
     TRANSFORMATION FOLLOWED BY A COMPLEX DIAGONAL TRANSFORMATION INTO        
     A SIMILAR UNITARY UPPER-HESSENBERG MATRIX WITH A REAL NONNEGATIVE        
     SUBDIAGONAL. 
     BAKCOMHES PERFORMS THE CORRESPONDING BACK TRANSFORMATION.      
        
        
 KEYWORDS:        
        
     COMPLEX EIGENPROBLEM,  
     REDUCTION HESSENBERG FORM,       
     HOUSEHOLDER'S TRANSFORMATION.    
        
        
 SUBSECTION: HSHCOMHES.     
        
        
 CALLING SEQUENCE:
        
     THE HEADING OF THE PROCEDURE READS:        
     "PROCEDURE" HSHCOMHES(AR, AI, N, EM, B, TR, TI, DEL); "VALUE" N;         
     "INTEGER" N; "ARRAY" AR, AI, EM, B, TR, TI, DEL;     
        
     THE MEANING OF THE FORMAL PARAMETERS IS:   
     AR,AI:     <ARRAY IDENTIFIER>;   
                "ARRAY" AR,AI[1:N,1:N];         
                ENTRY:      
                THE REAL PART AND THE IMAGINARY PART OF THE MATRIX TO BE      
                TRANSFORMED  MUST BE GIVEN IN THE ARRAYS AR AND AI, 
                RESPECTIVELY;         
                EXIT:       
                THE REAL PART AND THE IMAGINARY PART OF THE UPPER   
                TRIANGLE OF THE RESULTING UPPER-HESSENBERG MATRIX ARE         
                DELIVERED IN THE CORRESPONDING PARTS OF THE ARRAYS AR         
                AND AI, RESPECTIVELY; DATA FOR THE HOUSEHOLDER BACK-
                TRANSFORMATION ARE DELIVERED IN THE STRICT LOWER    
                TRIANGLES OF THE ARRAYS AR AND AI;        
1SECTION:3.2.1.2.2.2          (JUNE 1974)                         PAGE 2      
        
        
        
     N:         <ARITHMETIC EXPRESSION>;        
                THE ORDER OF THE GIVEN MATRIX;  
     EM:        <ARRAY IDENTIFIER>;   
                "ARRAY"EM[0:1];       
                ENTRY:      
                EM[0]: THE MACHINE PRECISION;   
                EM[1]: AN ESTIMATE OF THE NORM OF THE COMPLEX MATRIX;         
                  (OR, E.G. THE SUM OF THE INFINITY NORMS OF THE REAL         
                  (PART AND IMAGINARY PART OF THE MATRIX);
     B:         <ARRAY IDENTIFIER>;   
                "ARRAY"B[1:N-1];      
                EXIT:       
                THE REAL NONNEGATIVE SUBDIAGONAL OF THE RESULTING   
                UPPER-HESSENBERG MATRIX;        
     TR,TI:     <ARRAY IDENTIFIER>;   
                "ARRAY" TR,TI[1:N];   
                EXIT:       
                THE REAL PART AND THE IMAGINARY PART OF THE DIAGONAL
                ELEMENTS OF A DIAGONAL SIMILARITY TRANSFORMATION ARE
                DELIVERED IN THE ARRAYS TR AND TI, RESPECTIVELY; BY THIS      
                INFORMATION THE COMPLEX UPPER-HESSENBERG MATRIX IS  
                TRANSFORMED INTO A UPPER-HESSENBERG MATRIX WITH A REAL        
                SUBDIAGONAL;
     DEL:       <ARRAY IDENTIFIER>;   
                "ARRAY"DEL[1:N-2];    
                EXIT:       
                INFORMATION CONCERNING THE SEQUENCE OF HOUSEHOLDER  
                MATRICES.   
        
        
 PROCEDURES USED: 
        
     HSHCOMCOL = CP34355,   
     MATMAT    = CP34013,   
     ELMROWCOL = CP34028,   
     HSHCOMPRD = CP34356,   
     CARPOL    = CP34344,   
     COMMUL    = CP34341,   
     COMCOLCST = CP34352,   
     COMROWCST = CP34353.   
        
        
 RUNNING TIME: ROUGHLY PROPORTIONAL TO N CUBED. 
        
        
 LANGUAGE: ALGOL 60.        
        
        
 METHOD AND PERFORMANCE: SEE BAKCOMHES (THIS SECTION).    
        
        
1SECTION:3.2.1.2.2.2          (JUNE 1974)                         PAGE 3      
        
        
        
 SUBSECTION: BAKCOMHES.     
        
        
 CALLING SEQUENCE:
        
     THE HEADING OF THE PROCEDURE READS:        
     "PROCEDURE" BAKCOMHES(AR, AI, TR, TI, DEL, VR, VI, N, N1, N2); 
     "VALUE" N, N1, N2; "INTEGER" N, N1, N2;    
      "ARRAY" AR, AI, TR, TI, DEL, VR, VI;      
        
     THE MEANING OF THE FORMAL PARAMETERS IS:   
     AR,AI,TR,TI,DEL:       
                <ARRAY IDENTIFIER>;   
                "ARRAY" AR,AI[1:N,1:N];         
                "ARRAY" TR,TI[1:N];   
                "ARRAY"DEL[1:N-2];    
                ENTRY: THE DATA FOR THE BACKTRANSFORMATION AS PRODUCED        
                       BY HSHCOMHES;  
     VR,VI:     <ARRAY IDENTIFIER>;   
                "ARRAY" VR,VI[1:N,N1:N2];       
                ENTRY:      
                THE BACK TRANSFORMATION IS PERFORMED ON THE EIGENVECTORS      
                WITH THE REAL PARTS GIVEN IN ARRAY VR AND THE IMAGINARY       
                PARTS GIVEN IN ARRAY VI;        
                EXIT:       
                THE REAL PARTS AND IMAGINARY PARTS OF THE RESULTING 
                EIGENVECTORS ARE DELIVERED IN THE COLUMNS OF THE ARRAYS       
                VR AND VI, RESPECTIVELY;        
     N:         <ARITHMETIC EXPRESSION>;        
                THE ORDER OF THE MATRIX OF WHICH THE EIGENVECTORS ARE         
                CALCULATED; 
     N1,N2:     <ARITHMETIC EXPRESSION>;        
                THE EIGENVECTORS CORRESPONDING TO THE EIGENVALUES WITH        
                INDICES N1,...,N2 ARE TO BE TRANSFORMED;  
        
        
 PROCEDURES USED: 
        
     COMROWCST = CP34353,   
     HSHCOMPRD = CP34356.   
        
        
 RUNNING TIME: PROPORTIONAL TO (N2-N1) * N**2.  
        
        
 LANGUAGE : ALGOL 60.       
        
1SECTION:3.2.1.2.2.2          (JUNE 1974)                         PAGE 4      
        
        
        
 THE FOLLOWING HOLDS FOR BOTH PROCEDURES:       
        
        
 METHOD AND PERFORMANCE:    
        
     HSHCOMHES:   
     HOUSEHOLDER'S TRANSFORMATION (FOR COMPLEX MATRICES) IS A UNITARY         
     SIMILARITY TRANSFORMATION, WHICH TRANSFORMS A COMPLEX MATRIX INTO A      
     SIMILAR UPPER-HESSENBERG MATRIX (SEE WILKINSON, 1965, P. 347-349).       
     LET M BE A GIVEN COMPLEX MATRIX OF ORDER N, P THE TRANSFORMING 
     MATRIX AND H THE RESULTING UPPER-HESSENBERG MATRIX. SINCE P IS 
     UNITARY, WE THEN HAVE H = P''MP, WHERE '' STANDS FOR CONJUGATING         
     AND TRANSPOSING.  THE MATRIX P IS  THE PRODUCT OF  N-2 HOUSEHOLDER       
     MATRICES,  THESE  BEING  UNITARY  HERMITEAN  MATRICES  OF      
     THE FORM I - UU''/T, WHERE T IS A SCALAR (>0), AND U A COMPLEX 
     VECTOR. THE R-TH HOUSEHOLDER MATRIX, R=1,...,N-2, IS CHOSEN IN SUCH      
     A WAY THAT THE FIRST R ELEMENTS OF U VANISH, AND THE DESIRED ZEROS       
     ARE INTRODUCED IN THE LAST N-R-1 ELEMENTS OF THE R-TH COLUMN OF THE      
     MATRIX M. HOWEVER, IF THE EUCLIDEAN NORM OF THE LAST N-R-1 ELEMENTS      
     OF COLUMN R OF THE MATRIX M IS SMALLER THAN THE MACHINE PRECISION        
     TIMES A NORM OF THE MATRIX THEN THE R-TH TRANSFORMATION IS SKIPPED       
     (I.E. THE R-TH HOUSEHOLDER MATRIX IS REPLACED BY I). THE COMPLEX         
     DIAGONAL SIMILARITY TRANSFORMATION D TRANSFORMS THE UPPER-     
     HESSENBERG MATRIX H INTO AN UPPER-HESSENBERG MATRIX HR, WITH REAL        
     NONNEGATIVE ELEMENTS. THE DIAGONAL OF D IS CHOSEN IN SUCH A WAY
     THAT SUBDIAGONAL ELEMENTS OF H ARE TRANSFORMED INTO THEIR ABSOLUTE       
     VALUES (SEE MUELLER, 1966).      
     BAKCOMHES:   
     THE BACK TRANSFORMATION TRANSFORMS A COMPLEX VECTOR X INTO THE 
     COMPLEX VECTOR PDX. IF X IS AN EIGENVECTOR OF H  THEN PDX IS THE         
     CORRESPONDING EIGENVECTOR OF M. STARTING FROM THE VECTOR V=DX, THE       
     VECTOR PDX IS OBTAINED BY SUCCESSIVELY REPLACING V BY THE R-TH 
     HOUSEHOLDER MATRIX TIMES V, FOR R=N-2,...,1. THE RESULTING VECTOR        
     THEN EQUALS PDX.       
        
        
 REFERENCES:      
        
     MUELLER, D.J. (1966),  
     HOUSEHOLDER'S METHOD  FOR  COMPLEX MATRICES AND EIGENSYSTEMS   OF        
     HERMITIAN MATRICES,    
     NUMER.MATH., 8, P.72-92;         
        
     WILKINSON, J.H. (1965),
     THE ALGEBRAIC EIGENVALUE PROBLEM,
     CLARENDON PRESS, OXFORD;         
        
        
 EXAMPLE OF USE:  
        
     HSHCOMHES IS USED IN THE PROCEDURES EIGVALCOM AND EIGCOM.      
     BAKCOMHES IS USED IN THE PROCEDURE EIGCOM. 
     (SEE SECTION 3.3.2.2.2.).        
1SECTION:3.2.1.2.2.2          (JUNE 1974)                         PAGE 5      
        
        
        
 SOURCE TEXT(S) : 
        
0"CODE" 34366;    
     "PROCEDURE" HSHCOMHES(AR, AI, N, EM, B, TR, TI, DEL); "VALUE" N;         
     "INTEGER" N; "ARRAY" AR, AI, EM, B, TR, TI, DEL;     
     "BEGIN" "INTEGER" R, RM1, I, J, NM1;       
         "REAL" TOL, T, XR, XI;       
         "REAL" "PROCEDURE" MATMAT(L,U,I,J,A,B);"CODE" 34013;       
         "PROCEDURE" ELMROWCOL(L,U,I,J,A,B,X);"CODE" 34028;         
         "PROCEDURE" HSHCOMPRD(I,II,L,U,J,AR,AI,BR,BI,T);"CODE" 34356;        
         "PROCEDURE" COMCOLCST(L,U,J,AR,AI,XR,XI);"CODE" 34352;     
         "PROCEDURE" COMROWCST(L,U,I,AR,AI,XR,XI);"CODE" 34353;     
         "PROCEDURE" CARPOL(AR,AI,R,C,S);"CODE" 34344;    
         "PROCEDURE" COMMUL(AR,AI,BR,BI,RR,RI);"CODE" 34341;        
         "BOOLEAN" "PROCEDURE" HSHCOMCOL(L,U,J,AR,AI,TOL,K,C,S,T);  
         "CODE" 34355;      
         NM1:= N - 1; TOL:= (EM[0] * EM[1]) ** 2; RM1:= 1;
         "FOR" R:= 2 "STEP" 1 "UNTIL" NM1 "DO"  
         "BEGIN" "IF" HSHCOMCOL(R, N, RM1, AR, AI, TOL, B[RM1],     
             TR[R], TI[R], T) "THEN"  
             "BEGIN" "FOR" I:= 1 "STEP" 1 "UNTIL" N "DO"  
                 "BEGIN" XR:= (MATMAT(R, N, I, RM1, AI, AI) -       
                     MATMAT(R, N, I, RM1, AR, AR)) / T;   
                     XI:= ( - MATMAT(R, N, I, RM1, AR, AI) -        
                     MATMAT(R, N, I, RM1, AI, AR)) / T;   
                     ELMROWCOL(R, N, I, RM1, AR, AR, XR); 
                     ELMROWCOL(R, N, I, RM1, AR, AI, XI); 
                     ELMROWCOL(R, N, I, RM1, AI, AR, XI); 
                     ELMROWCOL(R, N, I, RM1, AI, AI, - XR)
                 "END";     
                 HSHCOMPRD(R, N, R, N, RM1, AR, AI, AR, AI, T);     
             "END";         
             DEL[RM1]:= T; RM1:= R    
         "END" FORR;        
         "IF" N > 1 "THEN" CARPOL(AR[N,NM1], AI[N,NM1], B[NM1],     
         TR[N], TI[N]); RM1:= 1; TR[1]:= 1; TI[1]:= 0;    
         "FOR" R:= 2 "STEP" 1 "UNTIL" N "DO"    
         "BEGIN" COMMUL(TR[RM1], TI[RM1], TR[R], TI[R], TR[R],      
             TI[R]); COMCOLCST(1, RM1, R, AR, AI, TR[R], TI[R]);    
             COMROWCST(R + 1, N, R, AR, AI, TR[R], - TI[R]);        
             RM1:= R        
         "END";   
     "END" HSHCOMHES;       
         "EOP"    
        
1SECTION:3.2.1.2.2.2          (JUNE 1974)                         PAGE 6      
        
        
0"CODE" 34367;    
     "PROCEDURE" BAKCOMHES(AR, AI, TR, TI, DEL, VR, VI, N, N1, N2); 
     "VALUE" N, N1, N2; "INTEGER" N, N1, N2;    
      "ARRAY" AR, AI, TR, TI, DEL, VR, VI;      
     "BEGIN" "INTEGER" I, R, RM1;     
         "REAL" H;
         "PROCEDURE" HSHCOMPRD(I,II,L,U,J,AR,AI,BR,BI,T);"CODE" 34356;        
         "PROCEDURE" COMROWCST(L,U,I,AR,AI,XR,XI);"CODE" 34353;     
         "FOR" I:= 2 "STEP" 1 "UNTIL" N "DO" COMROWCST(N1, N2, I, VR,         
         VI, TR[I], TI[I]); R:= N - 1;
         "FOR" RM1:= N - 2 "STEP" - 1 "UNTIL" 1 "DO"      
         "BEGIN" H:= DEL[RM1];        
             "IF" H > 0 "THEN" HSHCOMPRD(R, N, N1, N2, RM1, VR, VI, 
             AR, AI, H); R:= RM1      
         "END"    
     "END" BAKCOMHES;       
         "EOP"    
1SECTION:3.2.2.1.1            (DECEMBER 1975)                     PAGE 1      
        
        
        
 AUTHOR : D.T.WINTER        
        
        
 INSTITUTE : MATHEMATICAL CENTRE      
        
 RECEIVED : 731217
        
 BRIEF DESCRIPTION :        
        
     THIS SECTION CONTAINS THREE PROCEDURES :   
     1.  HSHREABID.         
         THIS PROCEDURE  TRANSFORMS  A GIVEN MATRIX  TO BIDIAGONAL FORM,      
         BY PREMULTIPLYING  AND  POSTMULTIPLYING  THE GIVEN MATRIX  WITH      
         ORTHOGONAL MATRICES.         
     2.  PSTTFMMAT.         
         THIS PROCEDURE CALCULATES THE  POSTMULTIPLYING  MATRIX FROM THE      
         DATA GENERATED BY HSHREABID. 
     3.  PRETFMMAT.         
         THIS PROCEDURE  CALCULATES THE  PREMULTIPLYING  MATRIX FROM THE      
         DATA GENERATED BY HSHREABID. 
        
 KEYWORDS :       
     HOUSEHOLDER'S TRANSFORMATION     
     BIDIAGONALISATION      
        
 SUBSECTION : HSHREABID     
        
 CALLING SEQUENCE :         
     THE HEADING OF THE PROCEDURE IS :
     "PROCEDURE" HSHREABID(A, M, N, D, B, EM);  
     "VALUE" M, N; "INTEGER" M, N; "ARRAY" A, D, B, EM;   
        
     THE MEANING OF THE FORMAL PARAMETERS IS :  
     A:  <ARRAY IDENTIFIER>;
         "ARRAY" A[1:M,1:N];
         ENTRY: THE GIVEN MATRIX;     
         EXIT: DATA CONCERNING THE  PREMULTIPLYING  AND  POSTMULTIPLYING      
             MATRICES;      
     M:  <ARITHMETIC EXPRESSION>;     
         THE NUMBER OF ROWS OF THE GIVEN MATRIX;
     N:  <ARITHMETIC EXPRESSION>;     
         THE NUMBER OF COLUMNS OF THE GIVEN MATRIX,       
         N SHOULD SATISFY N <= M;     
     D:  <ARRAY IDENTIFIER>;
         "ARRAY"D[1:N];     
         EXIT: THE DIAGONAL OF THE BIDIAGONAL MATRIX;     
     B:  <ARRAY IDENTIFIER>;
         "ARRAY"B[1:N];     
         EXIT: THE SUPERDIAGONAL OF THE BIDIAGONAL MATRIX IS DELIVERED        
         IN B[1:N-1];       
     EM: <ARRAY IDENTIFIER>;
         "ARRAY"EM[0:1];    
         ENTRY: EM[0]: THE MACHINE-PRECISION;   
         EXIT: EM[1]: THE INFINITY NORM OF THE GIVEN MATRIX.        
1SECTION:3.2.2.1.1            (JUNE 1974)                         PAGE 2      
        
        
        
 PROCEDURES USED :
     TAMMAT = CP34014       
     MATTAM = CP34015       
     ELMCOL = CP34023       
     ELMROW = CP34024       
        
        
 RUNNING TIME :   
     RUNNING TIME IS ROUGHLY PROPORTIONAL TO (M + N) * N * N        
        
        
 METHOD AND PERFORMANCE :   
     LET  US  ASSUME  A  GIVEN  MATRIX  A[ 1:M , 1:N ],  WITH  M  >=  N.      
     FIRSTLY WE PREMULTIPLY A WITH A HOUSEHOLDER MATRIX,  CHOSEN IN SUCH      
     A WAY THAT THE FIRST COLUMN OF THE RESULTING MATRIX A' IS ZERO WITH      
     THE EXCEPTION  OF THE FIRST  ELEMENT.  SECONDLY WE  POSTMULTIPLY A'      
     WITH A  HOUSEHOLDER  MATRIX SO  THAT THE FIRST ROW OF THE RESULTING      
     MATRIX  IS ZERO  WITH  THE EXCEPTION  OF  THE FIRST  TWO  ELEMENTS.      
     NOW  WE REMOVE  THE FIRST ROW AND COLUMN,  AND REPEAT  THIS PROCESS      
     UNTIL  THE  MATRIX  IS  TOTALLY  TRANSFORMED  TO  BIDIAGONAL  FORM.      
     THIS  PROCEDURE  IS A  REWRITING  OF A PART  OF  THE PROCEDURE  SVD      
     PUBLISHED  BY  G.H.GOLUB  AND  C.REINSCH[1]. HOWEVER IN CONTRAST TO      
     THEIR PROCEDURE, HERE WE SKIP A TRANSFORMATION IF THE COLUMN OR ROW      
     ON WHICH  OUR  ATTENTION  IS FOCUSSED  IS ALREADY  (NEARLY)  IN THE      
     DESIRED FORM,  I.E.  IF THE SUM OF THE SQUARES OF THE ELEMENTS THAT      
     OUGHT TO BE  ZERO  IS SMALLER THAN  A CERTAIN CONSTANT,  IN SVD THE      
     TRANSFORMATION  IS  SKIPPED  ONLY  IF THE  NORM  OF THE FULL ROW OR      
     COLUMN  IS SMALL ENOUGH.  OUR WAY SEEMS TO GIVE BETTER RESULTS,  AS      
     SOME  ILL-DEFINED  TRANSFORMATIONS  ARE  SKIPPED.  MOREOVER,  IF  A      
     TRANSFORMATION  IS  SKIPPED,   WE  DO  NOT  STORE  A  ZERO  IN  THE      
     DIAGONAL  OR  SUPERDIAGONAL,  BUT  WE  STORE  THE VALUE  THAT WOULD      
     HAVE  BEEN  FOUND  IF THE  COLUMN  OR  ROW  WAS IN THE DESIRED FORM      
     ALREADY.     
        
        
 LANGUAGE : ALGOL-60        
        
        
1SECTION:3.2.2.1.1            (DECEMBER 1975)                     PAGE 3      
        
        
        
 SUBSECTION : PSTTFMMAT     
        
        
 CALLING SEQUENCE :         
     THE HEADING OF THE PROCEDURE IS :
     "PROCEDURE" PSTTFMMAT(A, N, V, B);         
     "VALUE" N; "INTEGER" N; "ARRAY" A, V, B;   
        
     THE MEANING OF THE FORMAL PARAMETERS IS:   
     A:  <ARRAY IDENTIFIER>;
         "ARRAY"A[1:N,1:N]; 
         THE DATA  CONCERNING THE  POSTMULTIPLYING MATRIX,  AS GENERATED      
         BY HSHREABID;      
     N:  <ARITHMETIC EXPRESSION>;     
         THE NUMBER OF COLUMNS AND ROWS OF A;   
     V:  <ARRAY IDENTIFIER>;
         "ARRAY"V[1:N,1:N]; 
         EXIT: THE POSTMULTIPLYING MATRIX;      
     B:  <ARRAY IDENTIFIER>;
         "ARRAY"B[1:N];     
         THE SUPERDIAGONAL AS GENERATED BY HSHREABID.     
        
        
 PROCEDURES USED :
     MATMAT = CP34013       
     ELMCOL = CP34023       
        
        
 RUNNING TIME :   
     THE RUNNING TIME IS ABOUT PROPORTIONAL TO N ** 3     
        
        
 LANGUAGE : ALGOL 60        
        
        
1SECTION:3.2.2.1.1            (JUNE 1974)                         PAGE 4      
        
        
        
 SUBSECTION : PRETFMMAT     
        
        
 CALLING SEQUENCE :         
     THE HEADING OF THE PROCEDURE IS :
     "PROCEDURE" PRETFMMAT(A, M, N, D);         
     "VALUE" M, N; "INTEGER" M, N; "ARRAY" A, D;
        
     THE MEANING OF THE FORMAL PARAMETERS IS :  
     A:  <ARRAY IDENTIFIER>;
         "ARRAY"A[1:M,1:N]; 
         ENTRY:  THE  DATA   CONCERNING  THE  PREMULTIPLYING  MATRIX  AS      
             GENERATED BY HSHREABID   
         EXIT : THE PREMULTIPLYING MATRIX.      
     M:  <ARITHMETIC EXPRESSION>;     
         THE NUMBER OF ROWS OF A.     
     N:  <ARITHMETIC EXPRESSION>;     
         THE NUMBER OF COLUMNS OF A, N SHOULD SATISFY N <= M.       
     D:  <ARRAY IDENTIFIER>;
         "ARRAY"D[1:N];     
         THE DIAGONAL AS GENERATED BY HSHREABID.
        
        
 PROCEDURES USED :
     TAMMAT = CP34014       
     ELMCOL = CP34023       
        
        
 RUNNING TIME :   
     THE RUNNING TIME IS ABOUT PROPORTIONAL TO M * N * N  
        
        
 LANGUAGE : ALGOL-60        
        
        
 REFERENCES :     
     [1] WILKINSON, J.H. AND C.REINSCH
         HANDBOOK FOR AUTOMATIC COMPUTATION, VOL. 2       
         LINEAR ALGEBRA     
         HEIDELBERG (1971)  
        
        
 EXAMPLE OF USE : 
        
     FOR AN EXAMPLE OF USE ONE IS REFERRED TO SECTION 3.5.1.2       
        
1SECTION:3.2.2.1.1            (JUNE 1974)                         PAGE 5      
        
        
        
 SOURCE TEXT(S) : 
0"CODE" 34260;    
 "PROCEDURE" HSHREABID(A, M, N, D, B, EM);      
 "VALUE" M, N; "INTEGER" M, N; "ARRAY" A, D, B, EM;       
 "BEGIN" "INTEGER" I, J, I1;
     "REAL" NORM, MACHTOL, W, S, F, G, H;       
        
     "REAL" "PROCEDURE" TAMMAT(L, U, I, J, A, B);         
     "VALUE" L, U, I, J; "INTEGER" L, U, I, J; "ARRAY" A, B;        
     "CODE" 34014;
     "REAL" "PROCEDURE" MATTAM(L, U, I, J, A, B);         
     "VALUE" L, U, I, J; "ARRAY" A, B;
     "CODE" 34015;
        
     "PROCEDURE" ELMCOL(L, U, I, J, A, B, X);   
     "VALUE" L, U, I, J, X; "INTEGER" L, U, I, J; "REAL" X;         
      "ARRAY" A, B;         
     "CODE" 34023;
     "PROCEDURE" ELMROW(L, U, I, J, A, B, X);   
     "VALUE" L, U, I, J, X; "INTEGER" L, U, I, J; "REAL" X;         
      "ARRAY" A, B;         
     "CODE" 34024;
        
     NORM:= 0;    
     "FOR" I:= 1 "STEP" 1 "UNTIL" M "DO"        
     "BEGIN" W:= 0;         
         "FOR" J:= 1 "STEP" 1 "UNTIL" N "DO" W:= ABS(A[I,J]) + W;   
         "IF" W > NORM "THEN" NORM:= W
     "END";       
     MACHTOL:= EM[0] * NORM; EM[1]:= NORM;      
     "FOR" I:= 1 "STEP" 1 "UNTIL" N "DO"        
     "BEGIN" I1:= I + 1; S:= TAMMAT(I1, M, I, I, A, A);   
         "IF" S < MACHTOL "THEN" D[I]:= A[I,I] "ELSE"     
         "BEGIN" F:= A[I,I]; S:= F * F + S;     
             D[I]:= G:= "IF" F < 0 "THEN" SQRT(S) "ELSE" - SQRT(S); 
             H:= F * G - S; A[I,I]:= F - G;     
             "FOR" J:= I1 "STEP" 1 "UNTIL" N "DO"         
             ELMCOL(I, M, J, I, A, A, TAMMAT(I, M, I, J, A, A) / H) 
         "END";   
         "IF" I < N "THEN"  
         "BEGIN" S:= MATTAM(I1 + 1, N, I, I, A, A);       
             "IF" S < MACHTOL "THEN" B[I]:= A[I,I1] "ELSE"
             "BEGIN" F:= A[I,I1]; S:= F * F + S;
                 B[I]:= G:= "IF" F < 0 "THEN" SQRT(S) "ELSE" - SQRT(S);       
                 H:= F * G - S; A[I,I1]:= F - G;
                 "FOR" J:= I1 "STEP" 1 "UNTIL" M "DO"     
                 ELMROW(I1, N, J, I, A, A, MATTAM(I1, N, I, J, A, A) /        
                 H)         
             "END"
         "END"    
     "END"        
 "END" HSHREABID; 
         "EOP"    
1SECTION:3.2.2.1.1            (JUNE 1974)                         PAGE 6      
        
        
        
0"CODE" 34261;    
 "PROCEDURE" PSTTFMMAT(A, N, V, B);   
 "VALUE" N; "INTEGER" N; "ARRAY" A, V, B;       
 "BEGIN" "INTEGER" I, I1, J;
     "REAL" H;    
     "REAL" "PROCEDURE" MATMAT(L, U, I, J, A, B);         
     "VALUE" L, U, I, J; "INTEGER" L, U, I, J; "ARRAY" A, B;        
     "CODE" 34013;
     "PROCEDURE" ELMCOL(L, U, I, J, A, B, X);   
     "VALUE" L, U, I, J, X; "INTEGER" L, U, I, J; "REAL" X;         
      "ARRAY" A, B;         
     "CODE" 34023;
        
     I1:= N; V[N,N]:= 1;    
     "FOR" I:= N - 1 "STEP" - 1 "UNTIL" 1 "DO"  
     "BEGIN" H:= B[I] * A[I,I1]; "IF" H < 0 "THEN"        
         "BEGIN" "FOR" J:= I1 "STEP" 1 "UNTIL" N "DO" V[J,I]:= A[I,J] /       
             H;   
             "FOR" J:= I1 "STEP" 1 "UNTIL" N "DO"         
             ELMCOL(I1, N, J, I, V, V, MATMAT(I1, N, I, J, A, V))   
         "END";   
         "FOR" J:= I1 "STEP" 1 "UNTIL" N "DO" V[I,J]:= V[J,I]:= 0;  
         V[I,I]:= 1; I1:= I 
     "END"        
 "END" PSTTFMMAT; 
         "EOP"    
0"CODE" 34262;    
 "PROCEDURE" PRETFMMAT(A, M, N, D);   
 "VALUE" M, N; "INTEGER" M, N; "ARRAY" A, D;    
 "BEGIN" "INTEGER" I, I1, J;
     "REAL" G, H; 
     "REAL" "PROCEDURE" TAMMAT(L, U, I, J, A, B);         
     "VALUE" L, U, I, J; "INTEGER" L, U, I, J; "ARRAY" A, B;        
     "CODE" 34014;
     "PROCEDURE" ELMCOL(L, U, I, J, A, B, X);   
     "VALUE" L, U, I, J, X; "INTEGER" L, U, I, J; "REAL" X;         
      "ARRAY" A, B;         
     "CODE" 34023;
        
     "FOR" I:= N "STEP" - 1 "UNTIL" 1 "DO"      
     "BEGIN" I1:= I + 1; G:= D[I]; H:= G * A[I,I];        
         "FOR" J:= I1 "STEP" 1 "UNTIL" N "DO" A[I,J]:= 0; 
         "IF" H < 0 "THEN"  
         "BEGIN" "FOR" J:= I1 "STEP" 1 "UNTIL" N "DO"     
             ELMCOL(I, M, J, I, A, A, TAMMAT(I1, M, I, J, A, A) / H);         
             "FOR" J:= I "STEP" 1 "UNTIL" M "DO" A[J,I]:= A[J,I] / G
         "END"    
         "ELSE"   
         "FOR" J:= I "STEP" 1 "UNTIL" M "DO" A[J,I]:= 0;  
         A[I,I]:= A[I,I] + 1
     "END"        
 "END" PRETFMMAT; 
         "EOP"    
1SECTION 3.3.1.1.1            (JULY 1974)                         PAGE 1      
        
        
        
 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.
        
        
1SECTION 3.3.1.1.1            (DECEMBER 1975)                     PAGE 2      
        
        
        
 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;
        
     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:     
1SECTION 3.3.1.1.1            (JULY 1974)                         PAGE 3      
        
        
        
             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;      
        
     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;      
1SECTION 3.3.1.1.1            (JULY 1974)                         PAGE 4      
        
        
        
             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].      
        
        
        
1SECTION 3.3.1.1.1            (JULY 1974)                         PAGE 5      
        
        
        
 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; 
        
     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.      
        
        
1SECTION 3.3.1.1.1            (DECEMBER 1975)                     PAGE 6      
        
        
        
 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;     
        
     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;      
1SECTION 3.3.1.1.1            (JULY 1974)                         PAGE 7      
        
        
        
     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].  
1SECTION 3.3.1.1.1            (JULY 1974)                         PAGE 8      
        
        
        
 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];         
         "PROCEDURE" VALSYMTRI(D, BB, N, N1, N2, VAL, EM); "CODE" 34151;      
         "PROCEDURE" VECSYMTRI(D,B,N,N1,N2, VAL, VEC, EM); "CODE" 34152;      
        
         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  . 
        
        
1SECTION 3.3.1.1.1            (JULY 1974)                         PAGE 9      
        
        
        
 SOURCE TEXT(S):  
0"CODE" 34151;    
     "COMMENT" MCA 2311;    
     "PROCEDURE" VALSYMTRI(D, BB, N, N1, N2, VAL, EM);    
     "VALUE" N, N1, N2;     
     "INTEGER" N, N1, N2; "ARRAY" D, BB, VAL, EM;         
         "BEGIN" "INTEGER" K, COUNT;  
         "REAL" MAX, X, Y, MACHEPS, NORM, RE, MACHTOL, UB, LB, LAMBDA;        
        
         "REAL" "PROCEDURE" STURM;    
         "BEGIN" "INTEGER" P, I; "REAL" F;      
             COUNT:= COUNT + 1;       
             P:= K; F:= D[1] - X;     
             "FOR" I:= 2 "STEP" 1 "UNTIL" N "DO"
             "BEGIN" "IF" F <= 0 "THEN"         
                 "BEGIN" P:= P + 1;   
                     "IF" P > N "THEN" "GOTO" OUT         
                 "END"      
                 "ELSE" "IF" P < I - 1 "THEN"   
                 "BEGIN" LB:= X; "GOTO" OUT "END";        
                 "IF" ABS(F) < MACHTOL "THEN"   
                 F:= "IF" F <= 0 "THEN" - MACHTOL "ELSE" MACHTOL;   
                 F:= D[I] - X - BB[I - 1] / F   
             "END";         
             "IF" P = N "OR" F <= 0 "THEN"      
             "BEGIN" "IF" X < UB "THEN" UB:= X "END" "ELSE" LB:= X; 
        OUT: STURM:= "IF" P = N "THEN" F "ELSE" (N - P) * MAX       
         "END" STURM;       
        
         "BOOLEAN" "PROCEDURE" ZEROIN(X, Y, FX, TOLX); "CODE"34150; 
        
         MACHEPS:= EM[0]; NORM:= EM[1]; RE:= EM[2];       
         MACHTOL:= NORM * MACHEPS; MAX:= NORM / MACHEPS; COUNT:= 0; 
         UB:= 1.1 * NORM; LB:= - UB; LAMBDA:= UB;         
         "FOR" K:= N1 "STEP" 1 "UNTIL" N2 "DO"  
         "BEGIN" X:= LB; Y:= UB; LB:= -1.1  * NORM;       
             ZEROIN(X, Y, STURM, ABS(X) * RE + MACHTOL);  
             VAL[K]:= LAMBDA:= "IF" X > LAMBDA "THEN" LAMBDA "ELSE" X;        
             "IF" UB > X "THEN" UB:= "IF" X > Y "THEN" X "ELSE" Y   
         "END";   
         EM[3]:= COUNT      
     "END" VALSYMTRI;       
         "EOP"    
1SECTION 3.3.1.1.1            (JULY 1974)                        PAGE 10      
        
        
        
0"CODE" 34152;    
     "COMMENT" MCA 2312;    
     "PROCEDURE" VECSYMTRI(D, B, N, N1, N2, VAL, VEC, EM);
     "VALUE" N, N1, N2;     
     "INTEGER" N, N1, N2; "ARRAY" D, B, VAL, VEC, EM;     
     "BEGIN" "INTEGER" I, J, K, COUNT, MAXCOUNT, COUNTLIM, ORTH, IND;         
         "REAL" BI, BI1, U, W, Y, MI1, LAMBDA, OLDLAMBDA, ORTHEPS,  
         VALSPREAD, SPR, RES, MAXRES, OLDRES, NORM, NEWNORM, OLDNORM,         
         MACHTOL, VECTOL;   
         "ARRAY" M, P, Q, R, X[1:N];  
         "BOOLEAN" "ARRAY" INT[1:N];  
        
         "REAL" "PROCEDURE" VECVEC(L, U, SHIFT, A, B); "CODE" 34010;
         "PROCEDURE" ELMVECCOL(L, U, I, A, B, X); "CODE" 34021;     
         "REAL" "PROCEDURE" TAMVEC(L, U, I, A, B); "CODE" 34012;    
        
         NORM:= EM[1]; MACHTOL:= EM[0] * NORM; VALSPREAD:= EM[4] * NORM;      
         VECTOL:= EM[6] * NORM; COUNTLIM:= EM[8]; ORTHEPS:= SQRT(EM[0]);      
         MAXCOUNT:= IND:= 0; MAXRES:= 0;        
         "IF" N1 > 1 "THEN" 
         "BEGIN" ORTH:= EM[5]; OLDLAMBDA:= VAL[N1 - ORTH];
             "FOR" K:= N1 - ORTH + 1 "STEP" 1 "UNTIL" N1 - 1  "DO"  
             "BEGIN" LAMBDA:= VAL[K]; SPR:= OLDLAMBDA - LAMBDA;     
                 "IF" SPR < MACHTOL "THEN" LAMBDA:= OLDLAMBDA - MACHTOL;      
                 OLDLAMBDA:= LAMBDA   
             "END"
         "END" "ELSE" ORTH:= 1;       
         "FOR" K:= N1 "STEP" 1 "UNTIL" N2 "DO"  
         "BEGIN" LAMBDA:= VAL[K]; "IF" K > 1 "THEN"       
             "BEGIN" SPR:= OLDLAMBDA - LAMBDA;  
                 "IF" SPR < VALSPREAD "THEN"    
                 "BEGIN" "IF" SPR < MACHTOL "THEN"        
                     LAMBDA:= OLDLAMBDA - MACHTOL;        
                     ORTH:= ORTH +1   
                 "END" "ELSE" ORTH:= 1
             "END";         
             COUNT:= 0; U:= D[1] - LAMBDA; BI:= W:= B[1]; 
             "IF" ABS(BI) < MACHTOL "THEN" BI:= MACHTOL;  
             "FOR" I:= 1 "STEP" 1 "UNTIL" N - 1  "DO"     
             "BEGIN" BI1:= B[I + 1];  
                 "IF" ABS(BI1) < MACHTOL "THEN" BI1:= MACHTOL;      
                 "IF" ABS(BI) >= ABS(U) "THEN"  
                 "BEGIN" MI1:= M[I + 1]:= U / BI; P[I]:= BI;        
                     Y:= Q[I]:= D[I + 1] - LAMBDA; R[I]:= BI1;      
                     U:= W - MI1 * Y; W:= - MI1 * BI1; INT[I]:= "TRUE"        
                 "END"      
                 "ELSE"     
                 "BEGIN" MI1:= M[I + 1]:= BI / U; P[I]:= U; Q[I]:= W;         
                     R[I]:= 0; U:= D[I + 1] - LAMBDA - MI1 * W;W:= BI1;       
                     INT[I]:= "FALSE" 
                 "END";     
                 X[I]:= 1; BI:= BI1   
             "END" TRANSFORM
1SECTION 3.3.1.1.1            (JULY 1974)                        PAGE 11      
                                                                 ;  
        
        
             P[N]:= "IF" ABS(U) < MACHTOL "THEN" MACHTOL "ELSE" U;  
             Q[N]:= R[N]:= 0; X[N]:= 1; "GOTO" ENTRY;     
         ITERATE: W:= X[1]; 
             "FOR" I:= 2 "STEP" 1 "UNTIL" N "DO"
             "BEGIN" "IF" INT[I - 1] "THEN"     
                 "BEGIN" U:= W; W:= X[I - 1]:= X[I] "END" 
                 "ELSE" U:= X[I]; W:= X[I]:= U - M[I] * W 
             "END" ALTERNATE;         
         ENTRY: U:= W:= 0;  
             "FOR" I:= N "STEP" -1 "UNTIL" 1 "DO"         
             "BEGIN" Y:= U; U:= X[I]:= (X[I] - Q[I] * U - R[I] * W) /         
                 P[I]; W:= Y
             "END" NEXT ITERATION;    
             NEWNORM:= SQRT(VECVEC(1, N, 0, X, X)); "IF" ORTH > 1"THEN"       
             "BEGIN" OLDNORM:= NEWNORM;         
                 "FOR" J:= K - ORTH + 1 "STEP" 1 "UNTIL" K - 1 "DO" 
                 ELMVECCOL(1, N, J, X, VEC, -TAMVEC(1, N, J, VEC, X));        
                 NEWNORM:= SQRT(VECVEC(1, N, 0, X, X));   
                 "IF" NEWNORM < ORTHEPS * OLDNORM "THEN"  
                 "BEGIN" IND:= IND + 1; COUNT:= 1;        
                     "FOR" I:= 1 "STEP" 1 "UNTIL" IND - 1,
                     IND + 1 "STEP" 1 "UNTIL" N "DO" X[I]:= 0;      
                     X[IND]:= 1; "IF" IND = N "THEN" IND:= 0;       
                     "GOTO" ITERATE   
                 "END" NEW START      
             "END" ORTHOGONALISATION; 
             RES:= 1 / NEWNORM; "IF" RES > VECTOL "OR" COUNT = 0 "THEN"       
             "BEGIN" COUNT:= COUNT + 1; "IF" COUNT <= COUNTLIM "THEN"         
                 "BEGIN" "FOR" I:= 1 "STEP" 1 "UNTIL" N "DO"        
                     X[I]:= X[I] * RES; "GOTO" ITERATE    
                 "END"      
             "END";         
             "FOR" I:= 1 "STEP" 1 "UNTIL" N "DO" VEC[I,K]:= X[I] * RES;       
             "IF" COUNT > MAXCOUNT "THEN" MAXCOUNT:= COUNT;         
             "IF" RES > MAXRES "THEN" MAXRES:= RES; OLDLAMBDA:= LAMBDA        
         "END";   
         EM[5]:= ORTH; EM[7]:= MAXRES; EM[9]:= MAXCOUNT   
     "END" VECSYMTRI;       
         "EOP"    
1SECTION 3.3.1.1.1            (JULY 1974)                        PAGE 12      
        
        
        
     "CODE" 34160;
     "INTEGER" "PROCEDURE" QRIVALSYMTRI(D, BB, N, EM); "VALUE" N;   
     "INTEGER" N; "ARRAY" D, BB, EM;  
     "BEGIN" "INTEGER" I, I1, LOW, OLDLOW, N1, COUNT, MAX;
         "REAL" BBTOL, BBMAX, BBI, BBN1, MACHTOL, DN, DELTA, F, NUM,
         SHIFT, G, H, T, P, R, S, C, OLDG;      
         BBTOL:= (EM[2] * EM[1]) ** 2; MACHTOL:= EM[0] * EM[1];     
         MAX:= EM[4]; BBMAX:= 0; COUNT:= 0; OLDLOW:= N;   
         "FOR" N1:= N - 1 "WHILE" N > 0 "DO"    
         "BEGIN"  
             "FOR" I:= N, I - 1 "WHILE" ("IF" I >= 1 "THEN"         
             BB[I] > BBTOL "ELSE" "FALSE") "DO" LOW:= I;  
             "IF" LOW > 1 "THEN" "BEGIN" "IF" BB[LOW-1] > BBMAX "THEN"        
             BBMAX:= BB[LOW-1] "END"; 
             "IF" LOW = N "THEN" N:= N1 "ELSE"  
             "BEGIN" DN:= D[N]; DELTA:= D[N1] - DN;       
                 BBN1:= BB[N1];       
                 "IF" ABS(DELTA) < MACHTOL "THEN" R:= SQRT(BBN1) "ELSE"       
                 "BEGIN"    
                     F:= 2 / DELTA; NUM:= BBN1 * F;       
                     R:= -NUM / (SQRT(NUM * F + 1) + 1)   
                 "END";     
                 "IF" LOW = N1 "THEN" 
                 "BEGIN" D[N]:= DN + R; D[N1]:= D[N1] - R; N:= N - 2
                 "END"      
                 "ELSE"     
                 "BEGIN" COUNT:= COUNT + 1;     
                     "IF" COUNT > MAX "THEN" "GOTO" END;  
                     "IF" LOW < OLDLOW "THEN"   
                     "BEGIN" SHIFT:= 0; OLDLOW:= LOW "END"
                     "ELSE"  SHIFT:= DN + R;    
                     H:= D[LOW] - SHIFT;        
                     "IF" ABS(H) < MACHTOL "THEN" H:= "IF" H <= 0 "THEN"      
                     -MACHTOL "ELSE" MACHTOL;   
                     G:= H; T:= G * H;
                     BBI:= BB[LOW]; P:= T + BBI; I1:= LOW;
                     "FOR" I:= LOW + 1 "STEP" 1 "UNTIL" N "DO"      
                     "BEGIN" S:= BBI / P; C:= T / P;      
                         H:= D[I] - SHIFT - BBI / H;      
                         "IF" ABS(H) < MACHTOL "THEN" H:= "IF" H <= 0         
                         "THEN" -MACHTOL "ELSE" MACHTOL;  
                         OLDG:= G; G:= H * C; T:= G * H;  
                         D[I1]:= OLDG - G + D[I];         
                         BBI:= "IF" I = N "THEN" 0 "ELSE" BB[I];    
                         P:= T + BBI; BB[I1]:= S * P; I1:= I        
                     "END"; 
                     D[N]:= G + SHIFT 
                 "END" QRSTEP         
             "END"
         "END";   
      END: EM[3]:= SQRT(BBMAX); EM[5]:= COUNT; QRIVALSYMTRI:= N     
     "END" QRIVALSYMTRI;    
          "EOP"   
1SECTION 3.3.1.1.1            (JULY 1974)                        PAGE 13      
        
        
        
0"CODE" 34161;    
     "COMMENT" MCA 2321;    
     "INTEGER" "PROCEDURE" QRISYMTRI(A, N, D, B, BB, EM); "VALUE" N;
     "INTEGER" N; "ARRAY" A, D, B, BB, EM;      
     "BEGIN" "INTEGER" I, J, J1, K, M, M1, COUNT, MAX;    
         "REAL" BBMAX, R, S, SIN, T, C, COS, OLDCOS, G, P, W, TOL, TOL2,      
         LAMBDA, DK1, A0, A1;         
        
         "PROCEDURE" ROTCOL(L, U, I, J, A, C, S); "CODE" 34040;     
        
         TOL:= EM[2] * EM[1]; TOL2:= TOL * TOL; COUNT:= 0; BBMAX:= 0;         
         MAX:= EM[4]; M:= N;
      IN: K:= M; M1:= M - 1;
      NEXT: K:= K - 1; "IF" K > 0 "THEN"        
         "BEGIN" "IF" BB[K] >= TOL2 "THEN" "GOTO" NEXT;   
             "IF" BB[K] > BBMAX "THEN" BBMAX:= BB[K]      
         "END";   
         "IF" K = M1 "THEN" M:= M1 "ELSE"       
         "BEGIN"  
             T:= D[M] - D[M1]; R:= BB[M1];      
             "IF" ABS(T) < TOL "THEN" S:= SQRT(R) "ELSE"  
             "BEGIN" W:= 2 / T; S:= W * R / (SQRT(W * W * R + 1) + 1)         
             "END"; "IF" K = M - 2 "THEN"       
             "BEGIN" D[M]:= D[M] + S; D[M1]:= D[M1] - S;  
                 T:= - S / B[M1]; R:= SQRT(T * T + 1); COS:= 1 / R; 
                 SIN:= T / R; ROTCOL(1,N,M1,M,A,COS,SIN); M:= M - 2 
             "END"
             "ELSE"         
             "BEGIN" COUNT:= COUNT + 1;         
                 "IF" COUNT > MAX "THEN" "GOTO" END;      
                 LAMBDA:= D[M] + S; "IF" ABS(T) < TOL "THEN"        
                 "BEGIN" W:= D[M1] - S;         
                     "IF" ABS(W) < ABS(LAMBDA) "THEN" LAMBDA:= W    
                 "END";     
                 K:= K + 1; T:= D[K] - LAMBDA; COS:= 1; W:= B[K];   
                 P:= SQRT(T * T + W * W); J1:= K;         
                 "FOR" J:= K + 1 "STEP" 1 "UNTIL" M "DO"  
                 "BEGIN" OLDCOS:= COS; COS:= T / P; SIN:= W / P;    
                     DK1:= D[J] - LAMBDA; T:= OLDCOS * T; 
                     D[J1]:= (T + DK1) * SIN * SIN + LAMBDA + T;    
                     T:= COS * DK1 - SIN * W * OLDCOS; W:= B[J];    
                     P:= SQRT(T * T + W * W); G:= B[J1]:= SIN * P;  
                     BB[J1]:= G * G; ROTCOL(1, N, J1, J, A, COS, SIN);        
                     J1:= J 
                 "END";     
                 D[M]:= COS * T + LAMBDA; "IF" T < 0 "THEN" B[M1]:= - G       
             "END" QRSTEP   
         "END";   
         "IF" M > 0 "THEN" "GOTO" IN; 
      END: EM[3]:= SQRT(BBMAX); EM[5]:= COUNT; QRISYMTRI:= M        
     "END" QRISYMTRI;       
         "EOP"    
1SECTION 3.3.1.1.2            (JULY 1974)                         PAGE 1      
        
        
        
 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.
        
        
1SECTION 3.3.1.1.2            (JULY 1974)                         PAGE 2      
        
        
        
 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;    
        
     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.  
        
        
1SECTION 3.3.1.1.2            (JULY 1974)                         PAGE 3      
        
        
        
 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;    
        
     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.  
        
        
1SECTION 3.3.1.1.2            (JULY 1974)                         PAGE 4      
        
        
        
 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;         
        
     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.   
        
1SECTION 3.3.1.1.2            (JULY 1974)                         PAGE 5      
        
        
        
 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.
        
        
1SECTION 3.3.1.1.2            (JULY 1974)                         PAGE 6      
        
        
        
 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;         
        
     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.   
        
        
1SECTION 3.3.1.1.2            (JULY 1974)                         PAGE 7      
        
        
        
 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.        
        
        
1SECTION 3.3.1.1.2            (JULY 1974)                         PAGE 8      
        
        
        
 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;
        
     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.      
        
        
1SECTION 3.3.1.1.2            (JULY 1974)                         PAGE 9      
        
        
        
 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;
        
     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.  
        
        
1SECTION 3.3.1.1.2            (JULY 1974)                        PAGE 10      
        
        
        
 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;
        
     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.         
        
1SECTION 3.3.1.1.2            (JULY 1974)                        PAGE 11      
        
        
        
 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];
         "PROCEDURE" EIGSYM1(A, N, NUMVAL, VAL, VEC, EM); "CODE" 34156;       
        
         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"        
        
1SECTION 3.3.1.1.2            (JULY 1974)                        PAGE 12      
        
        
        
     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]; 
         "PROCEDURE" EIGVALSYM2(A, N, NUMVAL, VAL, EM); "CODE" 34153;         
        
         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  .
        
        
1SECTION 3.3.1.1.2            (JULY 1974)                        PAGE 13      
        
        
        
 SOURCE TEXT(S) : 
0"CODE" 34153;    
 "COMMENT" MCA 2313;        
 "PROCEDURE" EIGVALSYM2(A, N, NUMVAL, VAL, EM); "VALUE" N, NUMVAL;  
 "INTEGER" N, NUMVAL; "ARRAY" A, VAL, EM;       
 "BEGIN" "ARRAY" B, BB, D[1:N];       
        
         "PROCEDURE" TFMSYMTRI2(A, N, D, B, BB, EM); "CODE" 34140;  
         "PROCEDURE" VALSYMTRI(D, BB, N, N1, N2, VAL, EM); "CODE" 34151;      
        
         TFMSYMTRI2(A, N, D, B, BB, EM);        
         VALSYMTRI(D, BB, N, 1, NUMVAL, VAL, EM)
     "END" EIGVALSYM2;      
         "EOP"    
0"CODE" 34154;    
 "COMMENT" MCA 2314;        
 "PROCEDURE" EIGSYM2(A, N, NUMVAL, VAL, VEC, EM); "VALUE" N, NUMVAL;
 "INTEGER" N, NUMVAL; "ARRAY" A, VAL, VEC, EM;  
 "BEGIN" "ARRAY" B, BB, D[1:N];       
        
         "PROCEDURE" TFMSYMTRI2(A, N, D, B, BB, EM); "CODE" 34140;  
         "PROCEDURE" VALSYMTRI(D, BB, N, N1, N2, VAL, EM); "CODE" 34151;      
         "PROCEDURE" VECSYMTRI(D, B, N, N1, N2, VAL, VEC, EM);      
         "CODE" 34152;      
         "PROCEDURE" BAKSYMTRI2(A, N, N1, N2, VEC); "CODE" 34141;   
        
         TFMSYMTRI2(A, N, D, B, BB, EM);        
         VALSYMTRI(D, BB, N, 1, NUMVAL, VAL, EM);         
         VECSYMTRI(D, B, N, 1, NUMVAL, VAL, VEC, EM);     
         BAKSYMTRI2(A, N, 1, NUMVAL, VEC)       
     "END" EIGSYM2;         
         "EOP"    
0"CODE" 34155;    
 "COMMENT" MCA 2318;        
 "PROCEDURE" EIGVALSYM1(A, N, NUMVAL, VAL, EM); "VALUE" N, NUMVAL;  
 "INTEGER" N, NUMVAL; "ARRAY" A, VAL, EM;       
 "BEGIN" "ARRAY" B, BB, D[1:N];       
        
         "PROCEDURE" TFMSYMTRI1(A, N, D, B, BB, EM); "CODE" 34143;  
         "PROCEDURE" VALSYMTRI(D, BB, N, N1, N2, VAL, EM); "CODE" 34151;      
        
         TFMSYMTRI1(A, N, D, B, BB, EM);        
         VALSYMTRI(D, BB, N, 1, NUMVAL, VAL, EM)
     "END" EIGVALSYM1;      
         "EOP"    
1SECTION 3.3.1.1.2            (JULY 1974)                        PAGE 14      
        
        
        
0"CODE" 34156;    
 "COMMENT" MCA 2319;        
 "PROCEDURE" EIGSYM1(A, N, NUMVAL, VAL, VEC, EM); "VALUE" N, NUMVAL;
 "INTEGER" N, NUMVAL; "ARRAY" A, VAL, VEC, EM;  
 "BEGIN" "ARRAY" B, BB, D[1:N];       
        
         "PROCEDURE" TFMSYMTRI1(A, N, D, B, BB, EM); "CODE" 34143;  
         "PROCEDURE" VALSYMTRI(D, BB, N, N1, N2, VAL, EM); "CODE" 34151;      
         "PROCEDURE" VECSYMTRI(D, B, N, N1, N2, VAL, VEC, EM);      
         "CODE" 34152;      
         "PROCEDURE" BAKSYMTRI1(A, N, N1, N2, VEC); "CODE" 34144;   
        
         TFMSYMTRI1(A, N, D, B, BB, EM);        
         VALSYMTRI(D, BB, N, 1, NUMVAL, VAL, EM);         
         VECSYMTRI(D, B, N, 1, NUMVAL, VAL, VEC, EM);     
         BAKSYMTRI1(A, N, 1, NUMVAL, VEC)       
     "END" EIGSYM1;         
         "EOP"    
0"CODE" 34162;    
 "COMMENT" MCA 2322;        
 "INTEGER" "PROCEDURE" QRIVALSYM2(A, N, VAL, EM); "VALUE" N;        
 "INTEGER" N; "ARRAY" A, VAL, EM;     
 "BEGIN" "ARRAY" B, BB[1:N];
        
         "PROCEDURE" TFMSYMTRI2(A, N, D, B, BB, EM); "CODE" 34140;  
         "INTEGER" "PROCEDURE" QRIVALSYMTRI(D, BB, N, EM);
         "CODE" 34160;      
        
         TFMSYMTRI2(A, N, VAL, B, BB, EM);      
         QRIVALSYM2:= QRIVALSYMTRI(VAL, BB, N, EM)        
     "END" QRIVALSYM2;      
         "EOP"    
1SECTION 3.3.1.1.2            (JULY 1974)                        PAGE 15      
        
        
        
0"CODE" 34163;    
 "COMMENT" MCA 2323;        
 "INTEGER" "PROCEDURE" QRISYM(A, N, VAL, EM); "VALUE" N;  
 "INTEGER" N; "ARRAY" A, VAL, EM;     
 "BEGIN" "ARRAY" B, BB[1:N];
        
         "PROCEDURE" TFMSYMTRI2(A, N, D, B, BB, EM); "CODE" 34140;  
         "PROCEDURE" TFMPREVEC(A, N); "CODE" 34142;       
         "INTEGER" "PROCEDURE" QRISYMTRI(A, N, D, B, BB, EM);       
         "CODE" 34161;      
        
         TFMSYMTRI2(A, N, VAL, B, BB, EM); TFMPREVEC(A, N);         
         QRISYM:= QRISYMTRI(A, N, VAL, B, BB, EM)         
     "END" QRISYM;
         "EOP"    
0"CODE" 34164;    
 "COMMENT" MCA 2327;        
 "INTEGER" "PROCEDURE" QRIVALSYM1(A, N, VAL, EM); "VALUE" N;        
 "INTEGER" N; "ARRAY" A, VAL, EM;     
 "BEGIN" "ARRAY" B, BB[1 : N];        
        
         "PROCEDURE" TFMSYMTRI1(A, N, D, B, BB, EM); "CODE" 34143;  
         "INTEGER" "PROCEDURE" QRIVALSYMTRI(D, BB, N, EM);
         "CODE" 34160;      
        
         TFMSYMTRI1(A, N, VAL, B, BB, EM);      
         QRIVALSYM1:= QRIVALSYMTRI(VAL, BB, N, EM)        
     "END" QRIVALSYM1;      
         "EOP"    
1SECTION 3.3.1.2.1            (JULY 1974)                         PAGE 1      
        
        
        
 AUTHORS: T.J. DEKKER, W. HOFFMANN.   
        
        
 CONTRIBUTORS: W. HOFFMANN, S.P.N. VAN KAMPEN.  
        
        
 INSTITUTE: MATHEMATICAL CENTRE.      
        
        
 RECEIVED: 731115.
        
        
 BRIEF DESCRIPTION:         
        
     THIS SECTION CONTAINS FIVE PROCEDURES:     
     A) REAVALQRI CALCULATES  THE EIGENVALUES  OF A REAL UPPER-HESSEN-        
     BERG  MATRIX, PROVIDED THAT  ALL EIGENVALUES ARE REAL, BY MEANS  OF      
     SINGLE QR ITERATION;   
     B) REAVECHES CALCULATES  THE EIGENVECTOR  CORRESPONDING TO  A GIVEN      
     REAL EIGENVALUE  OF A  REAL UPPER-HESSENBERG MATRIX  BY MEANS  OF        
     INVERSE ITERATION;     
     C) REAQRI  CALCULATES  THE EIGENVALUES AND  EIGENVECTORS OF  A REAL      
     UPPER-HESSENBERG MATRIX, PROVIDED THAT  ALL EIGENVALUES ARE REAL,        
     BY MEANS OF SINGLE QR ITERATION; 
     D) COMVALQRI CALCULATES THE REAL  AND COMPLEX EIGENVALUES OF A REAL      
     UPPER-HESSENBERG MATRIX BY MEANS OF DOUBLE QR ITERATION;       
     E) COMVECHES CALCULATES THE  EIGENVECTOR  CORRESPONDING  TO A GIVEN      
     COMPLEX EIGENVALUE OF A REAL UPPER-HESSENBERG MATRIX  BY MEANS OF        
     INVERSE ITERATION.     
        
        
 KEYWORDS:        
        
     EIGENVALUE,  
     EIGENVECTOR, 
     UPPER-HESSENBERG MATRIX,         
     QR ITERATION,
     INVERSE ITERATION.     
        
        
1SECTION 3.3.1.2.1            (JULY 1974)                         PAGE 2      
        
        
        
 SUBSECTION: REAVALQRI.     
        
        
 CALLING SEQUENCE:
        
     THE HEADING OF THE PROCEDURE IS: 
     "INTEGER" "PROCEDURE" REAVALQRI(A, N, EM, VAL); "VALUE" N;     
     "INTEGER" N; "ARRAY" A, EM, VAL; 
        
     THE MEANING OF THE FORMAL PARAMETERS IS:   
     A:      <ARRAY IDENTIFIER>;      
             "ARRAY" A[1:N,1:N];      
             ENTRY: THE ELEMENTS OF THE  REAL UPPER-HESSENBERG  MATRIX        
                    MUST BE  GIVEN  IN THE UPPER TRIANGLE  AND THE FIRST      
                    SUBDIAGONAL OF ARRAY A;     
             EXIT:  THE HESSENBERG PART OF ARRAY A IS ALTERED;      
     N:      <ARITHMETIC EXPRESSION>; 
             THE ORDER OF THE GIVEN MATRIX;     
     EM:     <ARRAY IDENTIFIER>;      
             "ARRAY" EM[0:5];         
             ENTRY: EM[0], THE MACHINE PRECISION;         
                    EM[1], A NORM OF THE GIVEN MATRIX;    
                    EM[2], THE  RELATIVE TOLERANCE USED  FOR THE  QR
                           ITERATION; 
                           IF  THE ABSOLUTE  VALUE  OF SOME  SUBDIAGONAL      
                           ELEMENT IS  SMALLER THAN  EM[1] * EM[2], THEN      
                           THIS ELEMENT IS  NEGLECTED AND  THE MATRIX IS      
                           PARTITIONED;         
                    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[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 REAVALQRI;        
     VAL:    <ARRAY IDENTIFIER>;      
             "ARRAY" VAL[1:N];        
             THE EIGENVALUES OF  THE GIVEN MATRIX  ARE DELIVERED IN VAL.      
        
     MOREOVER:    
     REAVALQRI DELIVERS 0, PROVIDED THAT THE PROCESS IS COMPLETED WITHIN      
     EM[4] ITERATIONS; OTHERWISE REAVALQRI DELIVERS THE NUMBER OF EIGEN-      
     VALUES NOT CALCULATED. 
        
        
1SECTION 3.3.1.2.1            (JULY 1974)                         PAGE 3      
        
        
        
 PROCEDURES USED: 
        
     ROTCOL = CP34040,      
     ROTROW = CP34041.      
        
        
 RUNNING TIME: ROUGHLY PROPORTIONAL TO N CUBED. 
        
        
 LANGUAGE: ALGOL 60.        
        
        
 METHOD AND PERFORMANCE:    
        
     THE METHOD  USED IN  THE PROCEDURE  REAVALQRI  IS THE SINGLE QR
     ITERATION OF FRANCIS (SEE REF[1], P. 54, REF[2] P. 515 - 543  AND        
     REF[3]). THE EIGENVALUES  OF A REAL  UPPER-HESSENBERG MATRIX  ARE        
     CALCULATED, PROVIDED THAT  THE MATRIX HAS  REAL  EIGENVALUES  ONLY.      
        
        
 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.     
        
        
 EXAMPLE OF USE:  
        
     THE PROCEDURE REAVALQRI  IS USED  IN REAEIGVAL, SECTION  3.3.1.2.2.      
        
        
1SECTION 3.3.1.2.1            (JULY 1974)                         PAGE 4      
        
        
        
 SUBSECTION: REAVECHES.     
        
        
 CALLING SEQUENCE:
        
     THE HEADING OF THE PROCEDURE IS: 
     "PROCEDURE" REAVECHES(A, N, LAMBDA, EM, V); "VALUE" N, LAMBDA; 
     "INTEGER" N; "REAL" LAMBDA; "ARRAY" A, EM, V;        
        
     THE MEANING OF THE FORMAL PARAMETERS IS:   
     A:      <ARRAY IDENTIFIER>;      
             "ARRAY" A[1:N,1:N];      
             ENTRY: THE ELEMENTS OF THE  REAL UPPER-HESSENBERG  MATRIX        
                    MUST BE  GIVEN  IN THE UPPER TRIANGLE  AND THE FIRST      
                    SUBDIAGONAL OF ARRAY A;     
             EXIT:  THE HESSENBERG PART OF ARRAY A IS ALTERED;      
     N:      <ARITHMETIC EXPRESSION>; 
             THE ORDER OF THE GIVEN MATRIX;     
     LAMBDA: <ARITHMETIC EXPRESSION>; 
             THE GIVEN REAL EIGENVALUE OF THE UPPER-HESSENBERG MATRIX;        
     EM:     <ARRAY IDENTIFIER>;      
             "ARRAY" EM[0:9];         
             ENTRY: EM[0], THE MACHINE PRECISION;         
                    EM[1], A NORM OF THE GIVEN MATRIX;    
                    EM[6], THE TOLERANCE  USED FOR THE  EIGENVECTOR; THE      
                           INVERSE  ITERATION  ENDS IF  THE  EUCLIDIAN        
                           NORM  OF THE  RESIDUE VECTOR  IS SMALLER THAN      
                           EM[1] * EM[6];       
                    EM[8], 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[6] = "-10,     
                    EM[8] = 5;        
             EXIT:  EM[7], THE EUCLIDIAN NORM OF THE RESIDUE VECTOR OF        
                           THE CALCULATED EIGENVECTOR;    
                    EM[9], THE NUMBER OF  INVERSE ITERATIONS  PERFORMED;      
                           IF  EM[7] REMAINS  LARGER THAN  EM[1] * EM[6]      
                           DURING EM[8] ITERATIONS, THE VALUE  EM[8] + 1      
                           IS DELIVERED;        
     V:      <ARRAY IDENTIFIER>;      
             "ARRAY" V[1:N];
             THE CALCULATED EIGENVECTOR IS DELIVERED IN V.
1SECTION 3.3.1.2.1            (JULY 1974)                         PAGE 5      
        
        
        
 PROCEDURES USED: 
        
     VECVEC = CP34010,      
     MATVEC = CP34011.      
        
        
 RUNNING TIME: ROUGHLY PROPORTIONAL TO N SQUARED.         
        
        
 LANGUAGE: ALGOL 60.        
        
        
 METHOD AND PERFORMANCE:    
        
     THE PROCEDURE REAVECHES CALCULATES  AN EIGENVECTOR CORRESPONDING TO      
     A GIVEN  APPROXIMATE REAL  EIGENVALUE OF A REAL  UPPER-HESSENBERG        
     MATRIX, BY MEANS  OF INVERSE  ITERATION (SEE REF[1], P. 55, REF[2],      
     P. 619 - 629 AND REF[3]).        
        
        
 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.M. VARAH.       
          EIGENVECTORS OF A REAL MATRIX BY INVERSE ITERATION.       
          STANFORD UNIVERSITY, TECH. REP. NO. CS 34, 1966.
        
        
 EXAMPLE OF USE:  
        
     THE PROCEDURE  REAVECHES  IS  USED  IN REAEIG1, SECTION  3.3.1.2.2.      
        
        
1SECTION 3.3.1.2.1            (JULY 1974)                         PAGE 6      
        
        
        
 SUBSECTION: REAQRI.        
        
        
 CALLING SEQUENCE:
        
     THE HEADING OF THE PROCEDURE IS: 
     "INTEGER" "PROCEDURE" REAQRI(A, N, EM, VAL, VEC); "VALUE" N;   
     "INTEGER" N; "ARRAY" A, EM, VAL, VEC;      
        
     THE MEANING OF THE FORMAL PARAMETERS IS:   
     A:      <ARRAY IDENTIFIER>;      
             "ARRAY" A[1:N,1:N];      
             ENTRY: THE ELEMENTS OF THE  REAL UPPER-HESSENBERG  MATRIX        
                    MUST BE  GIVEN  IN THE UPPER TRIANGLE  AND THE FIRST      
                    SUBDIAGONAL OF ARRAY A;     
             EXIT:  THE HESSENBERG PART OF ARRAY A IS ALTERED;      
     N:      <ARITHMETIC EXPRESSION>; 
             THE ORDER OF THE GIVEN MATRIX;     
     EM:     <ARRAY IDENTIFIER>;      
             "ARRAY" EM[0:5];         
             ENTRY: EM[0], THE MACHINE PRECISION;         
                    EM[1], A NORM OF THE GIVEN MATRIX;    
                    EM[2], THE  RELATIVE TOLERANCE USED  FOR THE  QR
                           ITERATION; 
                           IF  THE ABSOLUTE  VALUE  OF SOME  SUBDIAGONAL      
                           ELEMENT IS  SMALLER THAN  EM[1] * EM[2], THEN      
                           THIS ELEMENT IS  NEGLECTED AND  THE MATRIX IS      
                           PARTITIONED;         
                    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[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  AND  THE  LAST N - K      
                           COLUMNS OF VEC  ARE APPROXIMATED  EIGENVALUES      
                           AND EIGENVECTORS OF THE GIVEN MATRIX, WHERE K      
                           IS DELIVERED IN REAQRI;        
     VAL:    <ARRAY IDENTIFIER>;      
             "ARRAY" VAL[1:N];        
             THE EIGENVALUES OF  THE GIVEN MATRIX  ARE DELIVERED IN VAL;      
     VEC:    <ARRAY IDENTIFIER>;      
             "ARRAY" VEC[1:N,1:N];    
             THE  CALCULATED  EIGENVECTORS, CORRESPONDING  TO THE EIGEN-      
             VALUES  IN ARRAY VAL[1:N], ARE DELIVERED  IN THE COLUMNS OF      
             ARRAY VEC.     
        
1SECTION 3.3.1.2.1            (JULY 1974)                         PAGE 7      
        
        
        
     MOREOVER:    
     REAQRI DELIVERS 0, PROVIDED  THAT THE  PROCESS IS COMPLETED  WITHIN      
     EM[4] ITERATIONS; OTHERWISE  REAQRI DELIVERS  THE NUMBER OF  EIGEN-      
     VALUES AND EIGENVECTORS NOT CALCULATED.    
        
        
 PROCEDURES USED: 
        
     MATVEC = CP34011,      
     ROTCOL = CP34040,      
     ROTROW = CP34041.      
        
        
 RUNNING TIME: ROUGHLY PROPORTIONAL TO N CUBED. 
        
        
 LANGUAGE: ALGOL 60.        
        
        
 METHOD AND PERFORMANCE:    
        
     THE  PROCEDURE  REAQRI  CALCULATES  THE  EIGENVALUES  OF AN  UPPER-      
     HESSENBERG MATRIX BY MEANS OF SINGLE QR ITERATION (SEE METHOD  
     AND PERFORMANCE OF REAVALQRI, THIS SECTION).         
     THE EIGENVECTORS  ARE CALCULATED  BY A  DIRECT  METHOD (SEE REF[1],      
     P. 55-56), IN CONTRAST WITH REAVECHES WHICH USES INVERSE ITERATION.      
     IF THE HESSENBERG MATRIX  IS NOT TOO ILL-CONDITIONED WITH RESPECT        
     TO ITS  EIGENVALUE  PROBLEM, THEN  THIS METHOD  YIELDS  NUMERICALLY      
     INDEPENDENT EIGENVECTORS AND IS COMPETITIVE WITH  INVERSE ITERATION      
     AS TO ACCURACY AND COMPUTATION TIME.       
     IF THE QR ITERATION PROCESS  IS NOT COMPLETED  WITHIN THE GIVEN
     NUMBER  OF  ITERATIONS, NOT ALL  EIGENVALUES  AND  EIGENVECTORS ARE      
     DELIVERED.   
     THE PROCEDURE REAQRI  SHOULD  BE USED ONLY  IF ALL EIGENVALUES  ARE      
     REAL.        
        
        
 REFERENCES:      
     [1]. T.J. DEKKER AND W. HOFFMANN.
          ALGOL 60 PROCEDURES IN NUMERICAL ALGEBRA, PART 2.         
          MC TRACT 23, 1968, MATH. CENTR., AMSTERDAM.     
        
        
 EXAMPLE OF USE:  
        
     THE PROCEDURE REAQRI IS USED IN REAEIG3, SECTION 3.3.1.2.2.    
        
        
1SECTION 3.3.1.2.1            (JULY 1974)                         PAGE 8      
        
        
        
 SUBSECTION: COMVALQRI.     
        
        
 CALLING SEQUENCE:
        
     THE HEADING OF THE PROCEDURE IS: 
     "INTEGER" "PROCEDURE" COMVALQRI(A, N, EM, RE, IM); "VALUE" N;  
     "INTEGER" N; "ARRAY" A, EM, RE, IM;        
        
     THE MEANING OF THE FORMAL PARAMETERS IS:   
     A:      <ARRAY IDENTIFIER>;      
             "ARRAY" A[1:N,1:N];      
             ENTRY: THE ELEMENTS OF THE  REAL UPPER-HESSENBERG  MATRIX        
                    MUST BE  GIVEN  IN THE UPPER TRIANGLE  AND THE FIRST      
                    SUBDIAGONAL OF ARRAY A;     
             EXIT:  THE HESSENBERG PART OF ARRAY A IS ALTERED;      
     N:      <ARITHMETIC EXPRESSION>; 
             THE ORDER OF THE GIVEN MATRIX;     
     EM:     <ARRAY IDENTIFIER>;      
             "ARRAY" EM[0:5];         
             ENTRY: EM[0], THE MACHINE PRECISION;         
                    EM[1], A NORM OF THE GIVEN MATRIX;    
                    EM[2], THE  RELATIVE TOLERANCE USED  FOR THE  QR
                           ITERATION; 
                           IF  THE ABSOLUTE  VALUE  OF SOME  SUBDIAGONAL      
                           ELEMENT IS  SMALLER THAN  EM[1] * EM[2], THEN      
                           THIS ELEMENT IS  NEGLECTED AND  THE MATRIX IS      
                           PARTITIONED;         
                    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[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 COMVALQRI;        
     RE,IM:  <ARRAY IDENTIFIER>;      
             "ARRAY" RE, IM[1:N];     
             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:    
     COMVALQRI DELIVERS 0, PROVIDED THAT THE PROCESS IS COMPLETED WITHIN      
     EM[4] ITERATIONS; OTHERWISE COMVALQRI DELIVERS THE NUMBER OF EIGEN-      
     VALUES NOT CALCULATED. 
        
        
1SECTION 3.3.1.2.1            (JULY 1974)                         PAGE 9      
        
        
        
 PROCEDURES USED: NONE.     
        
        
 RUNNING TIME: ROUGHLY PROPORTIONAL TO N CUBED. 
        
        
 LANGUAGE: ALGOL 60.        
        
        
        
        
 METHOD AND PERFORMANCE:    
        
     THE METHOD USED IN THE PROCEDURE COMVALQRI FOR CALCULATING THE REAL      
     AND COMPLEX EIGENVALUES  OF A REAL UPPER-HESSENBERG MATRIX IS THE        
     DOUBLE QR  ITERATION  OF  FRANCIS (SEE REF[1], P. 74,  REF[2]  
     P. 528 - 537 AND REF[3]).        
        
        
 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.     
        
        
 EXAMPLE OF USE:  
        
     THE COMPLEX EIGENVALUES  AND -VECTORS OF H, WITH N = 4 AND H[I,J] =      
     "IF" I = 1 "THEN" -1 "ELSE" "IF" I - J = 1 "THEN" 1 "ELSE" 0,   MAY      
     BE OBTAINED BY THE FOLLOWING PROGRAM:      
        
     "BEGIN" "INTEGER" I, J, M;       
         "ARRAY" A[1:4,1:4], RE, IM[1:4], EM[0:9];        
         "INTEGER" "PROCEDURE" COMVALQRI(A, N, EM, RE, IM);         
         "CODE" 34190;      
         "PROCEDURE" COMVECHES(A, N, LAMBDA, MU, EM, U, V);         
         "CODE" 34191;      
        
         EM[0]:= "-14; EM[2]:= "-13; EM[1]:= 4; EM[4]:= 40;         
         EM[6]:= "-10; EM[8]:= 5;     
         "FOR" I:= 1, 2, 3, 4 "DO" "FOR" J:= 1, 2, 3, 4 "DO" A[I,J]:=         
         "IF" I = 1 "THEN" -1 "ELSE" "IF" I - J = 1 "THEN" 1 "ELSE" 0;        
         M:= COMVALQRI(A, 4, EM, RE, IM); OUTPUT(61, "("D, /")", M);
1SECTION 3.3.1.2.1            (JULY 1974)                        PAGE 10      
        
        
        
         "FOR" J:= M + 1 "STEP" 1 "UNTIL" 4 "DO"
         "BEGIN" "INTEGER" K; "ARRAY" U, V[1:4];
             "FOR" I:= 1, 2, 3, 4 "DO" "FOR" K:= 1, 2, 3, 4 "DO"    
             A[I,K]:= "IF" I = 1 "THEN" -1 "ELSE"         
             "IF" I - K = 1 "THEN" 1 "ELSE" 0;  
             COMVECHES(A, 4, RE[J], IM[J], EM, U, V);     
             OUTPUT(61, "("/, 2(+.13D"+2D, 2B), 2/")", RE[J], IM[J]);         
             "FOR" I:= 1, 2, 3, 4 "DO"
             OUTPUT(61, "("21B, 2(+.13D"+2D, 2B), /")", U[I], V[I]) 
         "END";   
         OUTPUT(61, "("/, 2(.2D"+2D, /), 2(ZD, /)")",     
         EM[3], EM[7], EM[5], EM[9])  
     "END"        
        
     THE PROGRAM DELIVERS (THE RESULTS ARE CORRECT UP TO TWELVE DIGITS):      
        
     THE NUMBER OF NOT CALCULATED EIGENVALUES: 0
        
     THE EIGENVALUES AND -VECTORS:    
        
     +.3090169943750"+00  +.9510565162952"+00   
        
                          -.2527643931136"+00  -.4314048696688"+00  
                          -.4883989055049"+00  +.1070817869743"+00  
                          -.4908273055667"-01  +.4975850535950"+00  
                          +.4580641097602"+00  +.2004426884413"+00  
        
     +.3090169943750"+00  -.9510565162952"+00   
        
                          -.2527643931136"+00  +.4314048696688"+00  
                          -.4883989055049"+00  -.1070817869743"+00  
                          -.4908273055667"-01  -.4975850535950"+00  
                          +.4580641097602"+00  -.2004426884413"+00  
        
     -.8090169943749"+00  +.5877852522924"+00   
        
                          +.1095191711534"+00  -.4878581260468"+00  
                          -.3753586823743"+00  +.3303117611685"+00  
                          +.4978239349006"+00  -.4659753040772"-01  
                          -.4301373647081"+00  -.2549153731770"+00  
        
     -.8090169943749"+00  -.5877852522924"+00   
        
                          +.1095191711534"+00  +.4878581260468"+00  
                          -.3753586823743"+00  -.3303117611685"+00  
                          +.4978239349006"+00  +.4659753040772"-01  
                          -.4301373647081"+00  +.2549153731770"+00  
        
    THE ARRAY EM: EM[3] = .67"-22     
                  EM[7] = .17"-13     
                  EM[5] = 9 
                  EM[9] = 1 .         
        
        
1SECTION 3.3.1.2.1            (JULY 1974)                        PAGE 11      
        
        
        
 SUBSECTION: COMVECHES.     
        
        
 CALLING SEQUENCE:
        
     THE HEADING OF THE PROCEDURE IS: 
     "PROCEDURE" COMVECHES(A, N, LAMBDA, MU, EM, U, V);   
     "VALUE" N, LAMBDA, MU; 
     "INTEGER" N; "REAL" LAMBDA, MU; "ARRAY" A, EM, U, V; 
        
     THE MEANING OF THE FORMAL PARAMETERS IS:   
     A:      <ARRAY IDENTIFIER>;      
             "ARRAY" A[1:N,1:N];      
             ENTRY: THE ELEMENTS OF THE  REAL UPPER-HESSENBERG  MATRIX        
                    MUST BE  GIVEN  IN THE UPPER TRIANGLE  AND THE FIRST      
                    SUBDIAGONAL OF ARRAY A;     
             EXIT:  THE HESSENBERG PART OF ARRAY A IS ALTERED;      
     N:      <ARITHMETIC EXPRESSION>; 
             THE ORDER OF THE GIVEN MATRIX;     
     LAMBDA, MU:  
             <ARITHMETIC EXPRESSION>; 
             THE REAL AND IMAGINARY PART OF THE GIVEN EIGENVALUE;   
     EM:     <ARRAY IDENTIFIER>;      
             "ARRAY" EM[0:9];         
             ENTRY: EM[0], THE MACHINE PRECISION;         
                    EM[1], A NORM OF THE GIVEN MATRIX;    
                    EM[6], THE TOLERANCE  USED FOR THE  EIGENVECTOR; THE      
                           INVERSE  ITERATION  ENDS IF  THE  EUCLIDIAN        
                           NORM  OF THE  RESIDUE VECTOR  IS SMALLER THAN      
                           EM[1] * EM[6];       
                    EM[8], 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[6] = "-10,     
                    EM[8] = 5;        
             EXIT:  EM[7], THE EUCLIDIAN NORM OF THE RESIDUE VECTOR OF        
                           THE CALCULATED EIGENVECTOR;    
                    EM[9], THE NUMBER OF  INVERSE ITERATIONS  PERFORMED;      
                           IF  EM[7] REMAINS  LARGER THAN  EM[1] * EM[6]      
                           DURING EM[8] ITERATIONS, THE VALUE  EM[8] + 1      
                           IS DELIVERED;        
     U, V:   <ARRAY IDENTIFIER>;      
             "ARRAY" U, V[1:N];       
             THE REAL AND IMAGINARY  PARTS OF THE CALCULATED EIGENVECTOR      
             ARE DELIVERED IN THE ARRAYS U, V[1:N].       
        
1SECTION 3.3.1.2.1            (JULY 1974)                        PAGE 12      
        
 PROCEDURES USED: 
        
     VECVEC = CP34010,      
     MATVEC = CP34011,      
     TAMVEC = CP34012.      
        
        
 RUNNING TIME: ROUGHLY PROPORTIONAL TO N SQUARED.         
        
        
 LANGUAGE: ALGOL 60.        
        
        
        
        
 METHOD AND PERFORMANCE:    
        
     THE PROCEDURE COMVECHES CALCULATES  AN EIGENVECTOR CORRESPONDING TO      
     A GIVEN APPROXIMATE EIGENVALUE OF A REAL UPPER-HESSENBERG MATRIX,        
     BY   MEANS  OF   INVERSE  ITERATION  (SEE  REF[1],  P. 75,  REF[2],      
     P. 629 - 633 AND REF[3]).        
        
        
 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.M. VARAH.       
          EIGENVECTORS OF A REAL MATRIX BY INVERSE ITERATION.       
          STANFORD UNIVERSITY, TECH. REP. NO. CS 34, 1966.
        
        
 EXAMPLE OF USE:  
        
     SEE EXAMPLE OF USE OF COMVALQRI, THIS SECTION.       
        
        
 SOURCE TEXT(S) : 
0"CODE" 34180;    
     "COMMENT" MCA 2410;    
     "INTEGER" "PROCEDURE" REAVALQRI(A, N, EM, VAL); "VALUE" N;     
     "INTEGER" N; "ARRAY" A, EM, VAL; 
     "BEGIN" "INTEGER" N1, I, I1, J, Q, MAX, COUNT;       
         "REAL" DET, W, SHIFT, KAPPA, NU, MU, R, TOL, DELTA, MACHTOL, S;      
        
         "PROCEDURE" ROTCOL(L, U, I, J, A, C, S); "CODE" 34040;     
         "PROCEDURE" ROTROW(L, U, I, J, A, C, S); "CODE" 34041;     
         "COMMENT"
1SECTION 3.3.1.2.1            (JULY 1974)                        PAGE 13      
                                                                 ;  
        
        
         MACHTOL:= EM[0] * EM[1]; TOL:= EM[1] * EM[2]; MAX:= EM[4]; 
         COUNT:= 0; R:= 0;  
      IN: N1:= N - 1;       
         "FOR" I:= N, I - 1 "WHILE" ("IF" I >= 1 "THEN"   
         ABS(A[I + 1,I]) > TOL "ELSE" "FALSE") "DO" Q:= I;
         "IF" Q > 1 "THEN"  
         "BEGIN" "IF" ABS(A[Q,Q - 1]) > R "THEN"
             R:= ABS(A[Q,Q - 1])      
         "END";   
         "IF" Q = N "THEN"  
         "BEGIN" VAL[N]:= A[N,N]; N:= N1 "END"  
         "ELSE"   
         "BEGIN" DELTA:= A[N,N] - A[N1,N1]; DET:= A[N,N1] * A[N1,N];
             "IF" ABS(DELTA) < MACHTOL "THEN" S:= SQRT(DET) "ELSE"  
             "BEGIN" W:= 2 / DELTA; S:= W * W * DET + 1;  
                 S:= "IF" S <= 0 "THEN" -DELTA * .5 "ELSE"
                 W * DET / (SQRT(S) + 1)        
             "END";         
             "IF" Q = N1 "THEN"       
             "BEGIN" VAL[N]:= A[N,N] + S;       
                 VAL[N1]:= A[N1,N1] - S; N:= N - 2        
             "END"
             "ELSE"         
             "BEGIN" COUNT:= COUNT + 1;         
                 "IF" COUNT > MAX "THEN" "GOTO" OUT;      
                 SHIFT:= A[N,N] + S; "IF" ABS(DELTA) < TOL "THEN"   
                 "BEGIN" W:= A[N1,N1] - S;      
                     "IF" ABS(W) < ABS(SHIFT) "THEN" SHIFT:= W      
                 "END";     
                 A[Q,Q]:= A[Q,Q] - SHIFT;       
                 "FOR" I:= Q "STEP" 1 "UNTIL" N - 1 "DO"  
                 "BEGIN" I1:= I + 1; A[I1,I1]:= A[I1,I1] - SHIFT;   
                     KAPPA:= SQRT(A[I,I] ** 2 + A[I1,I] ** 2);      
                     "IF" I > Q "THEN"
                     "BEGIN" A[I,I - 1]:= KAPPA * NU;     
                         W:= KAPPA * MU         
                     "END"  
                     "ELSE" W:= KAPPA; MU:= A[I,I] / KAPPA;         
                     NU:= A[I1,I] / KAPPA; A[I,I]:= W;    
                     ROTROW(I1, N, I, I1, A, MU, NU);     
                     ROTCOL(Q, I, I, I1, A, MU, NU);      
                     A[I,I]:= A[I,I] + SHIFT    
                 "END";     
                 A[N,N - 1]:= A[N,N] * NU; A[N,N]:= A[N,N] * MU + SHIFT       
             "END"
         "END";   
         "IF" N > 0 "THEN" "GOTO" IN; 
      OUT: EM[3]:= R; EM[5]:= COUNT; REAVALQRI:= N        
     "END" REAVALQRI;       
         "EOP"    
1SECTION 3.3.1.2.1            (JULY 1974)                        PAGE 14      
        
        
        
0"CODE" 34181;    
     "COMMENT" MCA 2411;    
     "PROCEDURE" REAVECHES(A, N, LAMBDA, EM, V); "VALUE" N, LAMBDA; 
     "INTEGER" N; "REAL" LAMBDA; "ARRAY" A, EM, V;        
     "BEGIN" "INTEGER" I, I1, J, COUNT, MAX;    
         "REAL" M, R, NORM, MACHTOL, TOL;       
         "BOOLEAN" "ARRAY" P[1:N];    
        
         "REAL" "PROCEDURE" VECVEC(L, U, SHIFT, A, B); "CODE" 34010;
         "REAL" "PROCEDURE" MATVEC(L, U, I, A, B); "CODE" 34011;    
        
         NORM:= EM[1]; MACHTOL:= EM[0] * NORM; TOL:= EM[6] * NORM;  
         MAX:= EM[8]; A[1,1]:= A[1,1] - LAMBDA; 
      GAUSS: "FOR" I:= 1 "STEP" 1 "UNTIL" N - 1 "DO"      
         "BEGIN" I1:= I + 1; R:= A[I,I]; M:= A[I1,I];     
             "IF" ABS(M) < MACHTOL "THEN" M:= MACHTOL;    
             P[I]:= ABS(M) <= ABS(R); 
             "IF" P[I] "THEN"         
             "BEGIN" A[I1,I]:= M:= M / R;       
                 "FOR" J:= I1 "STEP" 1 "UNTIL" N "DO"     
                 A[I1,J]:= ("IF" J > I1 "THEN" A[I1,J]    
                 "ELSE" A[I1,J] - LAMBDA) - M * A[I,J]    
             "END"
             "ELSE"         
             "BEGIN" A[I,I]:= M; A[I1,I]:= M:= R / M;     
                 "FOR" J:= I1 "STEP" 1 "UNTIL" N "DO"     
                 "BEGIN" R:= ("IF" J > I1 "THEN" A[I1,J] "ELSE"     
                     A[I1,J] - LAMBDA);         
                     A[I1,J]:= A[I,J] - M * R; A[I,J]:= R 
                 "END"      
             "END"
         "END" GAUSS;       
         "IF" ABS(A[N,N]) < MACHTOL "THEN" A[N,N]:= MACHTOL;        
         "FOR" J:= 1 "STEP" 1 "UNTIL" N "DO" V[J]:= 1; COUNT:= 0;   
      FORWARD: COUNT:= COUNT + 1; "IF" COUNT > MAX "THEN" "GOTO" OUT;         
         "FOR" I:= 1 "STEP" 1 "UNTIL" N - 1 "DO"
         "BEGIN" I1:= I + 1;
             "IF" P[I] "THEN" V[I1]:= V[I1] - A[I1,I] * V[I] "ELSE" 
             "BEGIN" R:= V[I1]; V[I1]:= V[I] - A[I1,I] * R;         
                 V[I]:=R    
             "END"
         "END" FORWARD;     
      BACKWARD: "FOR" I:= N "STEP" -1 "UNTIL" 1 "DO"      
         V[I]:= (V[I] - MATVEC(I + 1, N, I, A, V)) / A[I,I];        
         R:= 1 / SQRT(VECVEC(1, N, 0, V, V));   
         "FOR" J:= 1 "STEP" 1 "UNTIL" N "DO" V[J]:= V[J] * R;       
         "IF" R > TOL "THEN" "GOTO" FORWARD;    
      OUT: EM[7]:= R; EM[9]:= COUNT   
     "END" REAVECHES;       
         "EOP"    
1SECTION 3.3.1.2.1            (JULY 1974)                        PAGE 15      
        
        
        
0"CODE" 34186;    
     "COMMENT" MCA 2416;    
     "INTEGER" "PROCEDURE" REAQRI(A, N, EM, VAL, VEC); "VALUE" N;   
     "INTEGER" N; "ARRAY" A, EM, VAL, VEC;      
     "BEGIN" "INTEGER" M1, I, I1, M, J, Q, MAX, COUNT;    
         "REAL" W, SHIFT, KAPPA, NU, MU, R, TOL, S, MACHTOL,        
         ELMAX, T, DELTA, DET;        
         "ARRAY" TF[1:N];   
        
         "REAL" "PROCEDURE" MATVEC(L, U, I, A, B); "CODE" 34011;    
         "PROCEDURE" ROTCOL(L, U, I, J, A, C, S); "CODE" 34040;     
         "PROCEDURE" ROTROW(L, U, I, J, A, C, S); "CODE" 34041;     
        
         MACHTOL:= EM[0] * EM[1]; TOL:= EM[1] * EM[2]; MAX:= EM[4]; 
         COUNT:= 0; ELMAX:= 0; M:= N; 
         "FOR" I:= 1 "STEP" 1 "UNTIL" N "DO"    
         "BEGIN" VEC[I,I]:= 1;        
             "FOR" J:= I + 1 "STEP" 1 "UNTIL" N "DO"      
             VEC[I,J]:= VEC[J,I]:= 0  
         "END";   
      IN: M1:= M - 1;       
         "FOR" I:= M, I - 1 "WHILE" ("IF" I >= 1 "THEN"   
         ABS(A[I + 1,I]) > TOL "ELSE" "FALSE") "DO" Q:= I;
         "IF" Q > 1 "THEN"  
         "BEGIN" "IF" ABS(A[Q,Q - 1]) > ELMAX "THEN"      
             ELMAX:= ABS(A[Q, Q - 1]) 
         "END";   
         "IF" Q = M "THEN"  
         "BEGIN" VAL[M]:= A[M,M]; M:= M1 "END"  
         "ELSE"   
         "BEGIN" DELTA:= A[M,M] - A[M1,M1]; DET:= A[M,M1] * A[M1,M];
             "IF" ABS(DELTA) < MACHTOL "THEN" S:= SQRT(DET) "ELSE"  
             "BEGIN" W:= 2 / DELTA; S:= W * W * DET + 1;  
                 S:= "IF" S <= 0 "THEN" -DELTA * .5 "ELSE"
                 W * DET / (SQRT(S) + 1)        
             "END";         
             "IF" Q = M1 "THEN"       
             "BEGIN" A[M,M]:= VAL[M]:= A[M,M] + S;        
                 A[Q,Q]:= VAL[Q]:= A[Q,Q] - S;  
                 T:= "IF" ABS(S) < MACHTOL "THEN"         
                 (S + DELTA) / A[M,Q] "ELSE" A[Q,M] / S;  
                 R:= SQRT(T * T + 1); NU:= 1 / R;         
                 MU:= -T * NU; A[Q,M]:= A[Q,M] - A[M,Q];  
                 ROTROW(Q + 2, N, Q, M, A, MU, NU);       
                 ROTCOL(1, Q - 1, Q, M, A, MU, NU);       
                 ROTCOL(1, N, Q, M, VEC, MU, NU); M:= M - 2         
             "END"
1SECTION 3.3.1.2.1            (JULY 1974)                        PAGE 16      
        
        
        
             "ELSE"         
             "BEGIN" COUNT:= COUNT + 1;         
                 "IF" COUNT > MAX "THEN" "GOTO" END;      
                 SHIFT:= A[M,M] + S; "IF" ABS(DELTA) < TOL "THEN"   
                 "BEGIN" W:= A[M1,M1] - S;      
                     "IF" ABS(W) < ABS(SHIFT) "THEN" SHIFT:= W      
                 "END";     
                 A[Q,Q]:= A[Q,Q] - SHIFT;       
                 "FOR" I:= Q "STEP" 1 "UNTIL" M1 "DO"     
                 "BEGIN" I1:= I + 1; A[I1,I1]:= A[I1,I1] - SHIFT;   
                     KAPPA:= SQRT(A[I,I] ** 2 + A[I1,I] ** 2);      
                     "IF" I > Q "THEN"
                     "BEGIN" A[I,I - 1]:= KAPPA * NU;     
                         W:= KAPPA * MU         
                     "END"  
                     "ELSE" W:= KAPPA; MU:= A[I,I] / KAPPA;         
                     NU:= A[I1,I] / KAPPA; A[I,I]:= W;    
                     ROTROW(I1, N, I, I1, A, MU, NU);     
                     ROTCOL(1, I, I, I1, A, MU, NU);      
                     A[I,I]:= A[I,I] + SHIFT;   
                     ROTCOL(1, N, I, I1, VEC, MU, NU)     
                 "END";     
                 A[M,M1]:= A[M,M] * NU; A[M,M]:= A[M,M] * MU + SHIFT
             "END"
         "END";   
         "IF" M > 0 "THEN" "GOTO" IN; 
         "FOR" J:= N "STEP" -1 "UNTIL" 2 "DO"   
         "BEGIN" TF[J]:= 1; T:= A[J,J];         
             "FOR" I:= J - 1 "STEP" -1 "UNTIL" 1 "DO"     
             "BEGIN" DELTA:= T - A[I,I];        
                 TF[I]:= MATVEC(I + 1, J, I, A, TF) /     
                 ("IF" ABS(DELTA) < MACHTOL "THEN" MACHTOL "ELSE" DELTA)      
             "END";         
             "FOR" I:= 1 "STEP" 1 "UNTIL" N "DO"
             VEC[I,J]:= MATVEC(1, J, I, VEC, TF)
         "END";   
      END: EM[3]:= ELMAX; EM[5]:= COUNT; REAQRI:= M       
     "END" REAQRI;
         "EOP"    
1SECTION 3.3.1.2.1            (JULY 1974)                        PAGE 17      
        
        
        
0"CODE" 34190;    
     "COMMENT" MCA 2420;    
     "INTEGER" "PROCEDURE" COMVALQRI(A, N, EM, RE, IM); "VALUE" N;  
     "INTEGER" N; "ARRAY" A, EM, RE, IM;        
     "BEGIN" "INTEGER" I, J, P, Q, MAX, COUNT, N1, P1, P2, IMIN1,   
         I1, I2, I3;        
         "REAL" DISC, SIGMA, RHO, G1, G2, G3, PSI1, PSI2, AA, E, K, 
         S, NORM, MACHTOL2, TOL, W;   
         "BOOLEAN" B;       
        
         NORM:= EM[1]; MACHTOL2:= (EM[0] * NORM) ** 2;    
         TOL:= EM[2] * NORM; MAX:= EM[4]; COUNT:= 0; W:= 0;         
      IN: "FOR" I:= N, I - 1 "WHILE"  
         ("IF" I >= 1 "THEN" ABS(A[I + 1,I]) > TOL "ELSE" "FALSE")  
         "DO" Q:= I; "IF" Q > 1 "THEN"
         "BEGIN" "IF" ABS(A[Q,Q - 1]) > W "THEN" W:= ABS(A[Q,Q - 1])
         "END";   
         "IF" Q >= N - 1 "THEN"       
         "BEGIN" N1:= N - 1; "IF" Q = N "THEN"  
             "BEGIN" RE[N]:= A[N,N]; IM[N]:= 0; N:= N1 "END"        
             "ELSE"         
             "BEGIN" SIGMA:= A[N,N] - A[N1,N1]; 
                 RHO:= -A[N,N1] * A[N1,N];      
                 DISC:= SIGMA ** 2 - 4 * RHO; "IF" DISC > 0 "THEN"  
                 "BEGIN" DISC:= SQRT(DISC);     
                     S:= -2 * RHO / (SIGMA + ("IF" SIGMA >= 0       
                     "THEN" DISC "ELSE" -DISC));
                     RE[N]:= A[N,N] + S;        
                     RE[N1]:= A[N1,N1] - S; IM[N]:= IM[N1]:= 0      
                 "END"      
                 "ELSE"     
                 "BEGIN" RE[N]:= RE[N1]:= (A[N1,N1] + A[N,N]) / 2;  
                     IM[N1]:= SQRT( -DISC) / 2; IM[N]:= -IM[N1]     
                 "END";     
                 N:= N - 2  
             "END"
         "END"    
         "ELSE"   
         "BEGIN" COUNT:= COUNT + 1; "IF" COUNT > MAX "THEN"         
             "GOTO" OUT; N1:= N - 1;  
             SIGMA:= A[N,N] + A[N1,N1] + SQRT(ABS(A[N1,N - 2] * A[N,N1])      
             * EM[0]); RHO:= A[N,N] * A[N1,N1] - A[N,N1] * A[N1,N]; 
             "FOR" I:= N - 1, I - 1 "WHILE"     
             ("IF" I - 1 >= Q "THEN" ABS(A[I,I - 1] *     
             A[I1,I] * (ABS(A[I,I] + A[I1,I1] - SIGMA) +  
             ABS(A[I + 2,I1]))) > ABS(A[I,I] * ((A[I,I] - SIGMA) +  
             A[I,I1] * A[I1,I] + RHO)) * TOL    
             "ELSE" "FALSE") "DO" P1:= I1:= I; P:= P1 - 1;
             P2:= P + 2;    
                                                               "COMMENT"      
1SECTION 3.3.1.2.1            (JULY 1974)                        PAGE 18      
      @ @     @   @                                              ;  
        
        
             "FOR" I:= P "STEP" 1 "UNTIL" N - 1 "DO"      
             "BEGIN" IMIN1:= I - 1; I1:= I + 1; I2:= I + 2;         
                 "IF" I = P "THEN"    
                 "BEGIN" G1:= A[P,P] * (A[P,P] - SIGMA) + A[P,P1] * 
                     A[P1,P] + RHO;   
                     G2:= A[P1,P] * (A[P,P] + A[P1,P1] - SIGMA);    
                     "IF" P1 <= N1 "THEN"       
                     "BEGIN" G3:= A[P1,P] * A[P2,P1]; A[P2,P]:= 0 "END"       
                     "ELSE" G3:= 0    
                 "END"      
                 "ELSE"     
                 "BEGIN" G1:= A[I,IMIN1]; G2:= A[I1,IMIN1];         
                     G3:= "IF" I2 <= N "THEN" A[I2,IMIN1] "ELSE" 0  
                 "END";     
                 K:= "IF" G1 >= 0 "THEN"        
                 SQRT(G1 ** 2 + G2 ** 2 + G3 ** 2) "ELSE" 
                 -SQRT(G1 ** 2 + G2 ** 2 + G3 ** 2);      
                 B:= ABS(K) > MACHTOL2;         
                 AA:= "IF" B "THEN" G1 / K + 1 "ELSE" 2;  
                 PSI1:= "IF" B "THEN" G2 / (G1 + K) "ELSE" 0;       
                 PSI2:= "IF" B "THEN" G3 / (G1 + K) "ELSE" 0;       
                 "IF" I ^= Q "THEN" A[I,IMIN1]:= "IF" I = P "THEN"  
                 -A[I,IMIN1] "ELSE" -K;         
                 "FOR" J:= I "STEP" 1 "UNTIL" N "DO"      
                 "BEGIN" E:= AA * (A[I,J] + PSI1 * A[I1,J] +        
                     ("IF" I2 <= N "THEN" PSI2 * A[I2,J] "ELSE" 0));
                     A[I,J]:= A[I,J] - E; A[I1,J]:= A[I1,J] - PSI1 * E;       
                     "IF" I2 <= N "THEN" A[I2,J]:= A[I2,J] - PSI2 * E         
                 "END";     
                 "FOR" J:= Q "STEP" 1 "UNTIL"   
                 ("IF" I2 <= N "THEN" I2 "ELSE" N) "DO"   
                 "BEGIN" E:= AA * (A[J,I] + PSI1 * A[J,I1] +        
                     ("IF" I2 <= N "THEN" PSI2 * A[J,I2] "ELSE" 0));
                     A[J,I]:= A[J,I] - E; A[J,I1]:= A[J,I1] - PSI1 * E;       
                     "IF" I2 <= N "THEN" A[J,I2]:= A[J,I2] - PSI2 * E         
                 "END";     
                 "IF" I2 <= N1 "THEN" 
                 "BEGIN" I3:= I + 3; E:= AA * PSI2 * A[I3,I2];      
                     A[I3,I]:= -E;    
                     A[I3,I1]:= -PSI1 * E;      
                     A[I3,I2]:= A[I3,I2] - PSI2 * E       
                 "END"      
             "END"
         "END";   
         "IF" N > 0 "THEN" "GOTO" IN; 
     OUT: EM[3]:= W; EM[5]:= COUNT; COMVALQRI:= N         
     "END" COMVALQRI;       
         "EOP"    
1SECTION 3.3.1.2.1            (JULY 1974)                        PAGE 19      
        
        
        
0"CODE" 34191;    
     "COMMENT" MCA 2421;    
     "PROCEDURE" COMVECHES(A, N, LAMBDA, MU, EM, U, V);   
     "VALUE" N, LAMBDA, MU; 
     "INTEGER" N; "REAL" LAMBDA, MU; "ARRAY" A, EM, U, V; 
     "BEGIN" "INTEGER" I, I1, J, COUNT, MAX;    
         "REAL" AA, BB, D, M, R, S, W, X, Y, NORM, MACHTOL, TOL;    
         "ARRAY" G, F[1:N]; 
         "BOOLEAN" "ARRAY" P[1:N];    
        
         "REAL" "PROCEDURE" VECVEC(L, U, SHIFT, A, B); "CODE" 34010;
         "REAL" "PROCEDURE" MATVEC(L, U, I, A, B); "CODE" 34011;    
         "REAL" "PROCEDURE" TAMVEC(L, U, I, A, B); "CODE" 34012;    
        
         NORM:= EM[1]; MACHTOL:= EM[0] * NORM; TOL:= EM[6] * NORM;  
         MAX:= EM[8];       
         "FOR" I:= 2 "STEP" 1 "UNTIL" N "DO"    
         "BEGIN" F[I - 1]:= A[I,I - 1]; A[I,1]:= 0 "END"; 
         AA:= A[1,1] - LAMBDA; BB:= -MU;        
         "FOR" I:= 1 "STEP" 1 "UNTIL" N - 1 "DO"
         "BEGIN" I1:= I + 1; M:= F[I];
             "IF" ABS(M) < MACHTOL "THEN" M:= MACHTOL;    
             A[I,I]:= M; D:= AA ** 2 + BB ** 2; P[I]:= ABS(M) < SQRT(D);      
             "IF" P[I] "THEN"         
             "BEGIN" "COMMENT" A[I,J] * FACTOR AND A[I1,J] - A[I,J];
                 F[I]:= R:= M * AA / D; G[I]:= S:= -M * BB / D;     
                 W:= A[I1,I]; X:= A[I,I1]; A[I1,I]:= Y:= X * S + W * R;       
                 A[I,I1]:= X:= X * R - W * S;   
                 AA:= A[I1,I1] - LAMBDA - X; BB:= -(MU + Y);        
                 "FOR" J:= I + 2 "STEP" 1 "UNTIL" N "DO"  
                 "BEGIN" W:= A[J,I]; X:= A[I,J];
                     A[J,I]:= Y:= X * S + W * R;
                     A[I,J]:= X:= X * R - W * S; A[J,I1]:= -Y;      
                     A[I1,J]:= A[I1,J] - X      
                 "END"      
             "END"
             "ELSE"         
             "BEGIN" "COMMENT" INTERCHANGE A[I1,J] AND    
                 A[I,J] - A[I1,J] * FACTOR;     
                 F[I]:= R:= AA / M; G[I]:= S:= BB / M;    
                 W:= A[I1,I1] - LAMBDA; AA:= A[I,I1] - R * W - S * MU;        
                 A[I,I1]:= W; BB:= A[I1,I] - S * W + R * MU;        
                 A[I1,I]:= -MU;       
                 "FOR" J:= I + 2 "STEP" 1 "UNTIL" N "DO"  
                 "BEGIN" W:= A[I1,J]; A[I1,J]:= A[I,J] - R * W;     
                     A[I,J]:= W;      
                     A[J,I1]:= A[J,I] - S * W; A[J,I]:= 0 
                 "END"      
             "END"
         "END"    
1SECTION 3.3.1.2.1            (JULY 1974)                        PAGE 20      
                                                                 ;  
        
        
         P[N]:= "TRUE"; D:= AA ** 2 + BB ** 2; "IF" D < MACHTOL ** 2
         "THEN" "BEGIN" AA:= MACHTOL; BB:= 0; D:= MACHTOL ** 2 "END";         
         A[N,N]:= D; F[N]:= AA; G[N]:= -BB;     
         "FOR" I:= 1 "STEP" 1 "UNTIL" N "DO"    
         "BEGIN" U[I]:= 1; V[I]:= 0 "END";      
         COUNT:= 0;         
      FORWARD: "IF" COUNT > MAX "THEN" "GOTO" OUTM;       
         "FOR" I:= 1 "STEP" 1 "UNTIL" N "DO"    
         "BEGIN" "IF" P[I] "THEN"     
              "BEGIN" W:= V[I]; V[I]:= G[I] * U[I] + F[I] * W;      
                 U[I]:= F[I] * U[I] - G[I] * W; "IF" I < N "THEN"   
                 "BEGIN" V[I + 1]:= V[I + 1] - V[I];      
                     U[I + 1]:= U[I + 1] - U[I] 
                 "END"      
             "END"
             "ELSE"         
             "BEGIN" AA:= U[I + 1]; BB:= V[I + 1];        
                 U[I + 1]:= U[I] - (F[I] * AA - G[I] * BB); U[I]:= AA;        
                 V[I + 1]:= V[I] - (G[I] * AA + F[I] * BB); V[I]:= BB         
             "END"
         "END" FORWARD;     
      BACKWARD: "FOR" I:= N "STEP" -1 "UNTIL" 1 "DO"      
         "BEGIN" I1:= I + 1;
             U[I]:= (U[I] - MATVEC(I1, N, I, A, U) + ("IF" P[I] "THEN"        
             TAMVEC(I1, N, I, A, V) "ELSE" A[I1,I] * V[I1])) / A[I,I];        
             V[I]:= (V[I] - MATVEC(I1, N, I, A, V) - ("IF" P[I] "THEN"        
             TAMVEC(I1, N, I, A, U) "ELSE" A[I1,I] * U[I1])) / A[I,I]         
         "END" BACKWARD;    
      NORMALISE: W:= 1 / SQRT(VECVEC(1, N, 0, U, U) +     
         VECVEC(1, N, 0, V, V));      
         "FOR" J:= 1 "STEP" 1 "UNTIL" N "DO"    
         "BEGIN" U[J]:= U[J] * W; V[J]:= V[J] * W "END";  
         COUNT:= COUNT + 1; "IF" W > TOL "THEN" "GOTO" FORWARD;     
      OUTM: EM[7]:= W; EM[9]:= COUNT  
     "END" COMVECHES;       
         "EOP"    
1SECTION 3.3.1.2.2            (JULY 1974)                         PAGE 1      
        
        
        
 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.
        
1SECTION 3.3.1.2.2            (JULY 1974)                         PAGE 2      
        
        
        
 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; 
        
     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.      
1SECTION 3.3.1.2.2            (JULY 1974)                         PAGE 3      
        
        
        
        
 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;      
        
     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;     
1SECTION 3.3.1.2.2            (JULY 1974)                         PAGE 4      
        
        
        
     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;        
        
1SECTION 3.3.1.2.2            (JULY 1974)                         PAGE 5      
        
        
        
     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.        
        
        
1SECTION 3.3.1.2.2            (JULY 1974)                         PAGE 6      
        
        
        
 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;      
        
     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.      
        
        
1SECTION 3.3.1.2.2            (JULY 1974)                         PAGE 7      
        
        
        
 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.        
        
        
1SECTION 3.3.1.2.2            (JULY 1974)                         PAGE 8      
        
        
        
 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;        
        
     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.   
        
        
1SECTION 3.3.1.2.2            (JULY 1974)                         PAGE 9      
        
        
        
 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;   
        
     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;        
1SECTION 3.3.1.2.2            (JULY 1974)                        PAGE 10      
        
        
        
             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.         
        
        
1SECTION 3.3.1.2.2            (JULY 1974)                        PAGE 11      
        
        
        
 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.
        
        
1SECTION 3.3.1.2.2            (JULY 1974)                        PAGE 12      
        
        
        
 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];      
         "INTEGER" "PROCEDURE" REAEIG3(A, N, EM, VAL, VEC);         
         "CODE" 34187;      
        
         "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 .  
1SECTION 3.3.1.2.2            (JULY 1974)                        PAGE 13      
        
        
        
     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];        
         "INTEGER" "PROCEDURE" COMEIGVAL(A, N, EM, RE, IM);         
         "CODE" 34192;      
        
         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:    
0"CODE" 34182;    
     "COMMENT" MCA 2412;    
     "INTEGER" "PROCEDURE" REAEIGVAL(A, N, EM, VAL); "VALUE" N;     
     "INTEGER" N; "ARRAY" A, EM, VAL; 
     "BEGIN" "INTEGER" I, J; "REAL" R;
         "ARRAY" D[1:N]; "INTEGER" "ARRAY" INT, INT0[1:N];
        
         "PROCEDURE" TFMREAHES(A, N, EM, INT); "CODE" 34170;        
         "PROCEDURE" EQILBR(A, N, EM, D, INT); "CODE" 34173;        
         "INTEGER" "PROCEDURE" REAVALQRI(A, N, EM, VAL); "CODE" 34180;        
        
         EQILBR(A, N, EM, D, INT0); TFMREAHES(A, N, EM, INT);       
         J:= REAEIGVAL:= REAVALQRI(A, N, EM, VAL);        
         "FOR" I:= J + 1 "STEP" 1 "UNTIL" N "DO"
         "FOR" J:= I + 1 "STEP" 1 "UNTIL" N "DO"
         "BEGIN" "IF" VAL[J] > VAL[I] "THEN"    
             "BEGIN" R:= VAL[I]; VAL[I]:= VAL[J]; VAL[J]:= R "END"  
         "END"    
     "END" REAEIGVAL;       
         "EOP"    
1SECTION 3.3.1.2.2            (JULY 1974)                        PAGE 14      
        
        
        
0"CODE" 34184;    
     "COMMENT" MCA 2414;    
     "INTEGER" "PROCEDURE" REAEIG1(A, N, EM, VAL, VEC); "VALUE" N;  
     "INTEGER" N; "ARRAY" A, EM, VAL, VEC;      
     "BEGIN" "INTEGER" I, K, MAX, J, L;         
         "REAL" RESIDU, R, MACHTOL;   
         "ARRAY" D, V[1:N], B[1:N,1:N];         
         "INTEGER" "ARRAY" INT, INT0[1:N];      
        
         "PROCEDURE" TFMREAHES(A, N, EM, INT); "CODE" 34170;        
         "PROCEDURE" BAKREAHES2(A, N, N1, N2, INT, VEC); "CODE" 34172;        
         "PROCEDURE" EQILBR(A, N, EM, D, INT); "CODE" 34173;        
         "PROCEDURE" BAKLBR(N, N1, N2, D, INT, VEC); "CODE" 34174;  
         "INTEGER" "PROCEDURE" REAVALQRI(A, N, EM, VAL); "CODE" 34180;        
         "PROCEDURE" REAVECHES(A, N, LAMBDA, EM, V); "CODE" 34181;  
         "PROCEDURE" REASCL(A, N, N1, N2); "CODE" 34183;  
        
         RESIDU:= 0; MAX:= 0; EQILBR(A, N, EM, D, INT0);  
         TFMREAHES(A, N, EM, INT);    
         "FOR" I:= 1 "STEP" 1 "UNTIL" N "DO"    
         "FOR" J:= ("IF" I = 1 "THEN" 1 "ELSE" I - 1)     
         "STEP" 1 "UNTIL" N "DO" B[I,J]:= A[I,J];         
         K:= REAEIG1:= REAVALQRI(B, N, EM, VAL);
         "FOR" I:= K + 1 "STEP" 1 "UNTIL" N "DO"
         "FOR" J:= I + 1 "STEP" 1 "UNTIL" N "DO"
         "BEGIN" "IF" VAL[J] > VAL[I] "THEN"    
             "BEGIN" R:= VAL[I]; VAL[I]:= VAL[J]; VAL[J]:= R "END"  
         "END";   
         MACHTOL:= EM[0] * EM[1];     
         "FOR" L:= K + 1 "STEP" 1 "UNTIL" N "DO"
         "BEGIN" "IF" L > 1 "THEN"    
             "BEGIN" "IF" VAL[L - 1] - VAL[L] < MACHTOL "THEN"      
                 VAL[L]:= VAL[L - 1] - MACHTOL  
             "END";         
             "FOR" I:= 1 "STEP" 1 "UNTIL" N "DO"
             "FOR" J:= ("IF" I = 1 "THEN" 1 "ELSE" I - 1) 
             "STEP" 1 "UNTIL" N "DO" B[I,J]:= A[I,J];     
             REAVECHES(B, N, VAL[L], EM, V);    
             "IF" EM[7] > RESIDU "THEN" RESIDU:= EM[7];   
             "IF" EM[9] > MAX "THEN" MAX:= EM[9];         
             "FOR" J:= 1 "STEP" 1 "UNTIL" N "DO" VEC[J,L]:= V[J]    
         "END";   
         EM[7]:= RESIDU; EM[9]:= MAX; 
         BAKREAHES2(A, N, K + 1, N, INT, VEC);  
         BAKLBR(N, K + 1, N, D, INT0, VEC);     
         REASCL(VEC, N, K + 1, N)     
     "END" REAEIG1;         
         "EOP"    
1SECTION 3.3.1.2.2            (JULY 1974)                        PAGE 15      
        
        
        
0"CODE" 34187;    
     "COMMENT" MCA 2417;    
     "INTEGER" "PROCEDURE" REAEIG3(A, N, EM, VAL, VEC); "VALUE" N;  
     "INTEGER" N; "ARRAY" A, EM, VAL, VEC;      
     "BEGIN" "INTEGER" I; "REAL" S;   
         "INTEGER" "ARRAY" INT, INT0[1:N]; "ARRAY" D[1:N];
        
         "PROCEDURE" TFMREAHES(A, N, EM, INT); "CODE" 34170;        
         "PROCEDURE" BAKREAHES2(A, N, N1, N2, INT, VEC); "CODE" 34172;        
         "PROCEDURE" EQILBR(A, N, EM, D, INT); "CODE" 34173;        
         "PROCEDURE" BAKLBR(N, N1, N2, D, INT, VEC); "CODE" 34174;  
         "PROCEDURE" REASCL(A, N, N1, N2); "CODE" 34183;  
         "INTEGER" "PROCEDURE" REAQRI(A, N, EM, VAL, VEC); "CODE" 34186;      
        
         EQILBR(A, N, EM, D, INT0); TFMREAHES(A, N, EM, INT);       
         I:= REAEIG3:= REAQRI(A, N, EM, VAL, VEC);        
         "IF" I = 0 "THEN"  
         "BEGIN" BAKREAHES2(A, N, 1, N, INT, VEC);        
             BAKLBR(N, 1, N, D, INT0, VEC); REASCL(VEC, N, 1, N)    
         "END"    
     "END" REAEIG3;         
         "EOP"    
0"CODE" 34192;    
     "COMMENT" MCA 2422;    
     "INTEGER" "PROCEDURE" COMEIGVAL(A, N, EM, RE, IM); "VALUE" N;  
     "INTEGER" N; "ARRAY" A, EM, RE, IM;        
     "BEGIN" "INTEGER" "ARRAY" INT, INT0[1:N];  
         "ARRAY" D[1:N];    
        
         "PROCEDURE" EQILBR(A, N, EM, D, INT); "CODE" 34173;        
         "PROCEDURE" TFMREAHES(A, N, EM, INT); "CODE" 34170;        
         "INTEGER" "PROCEDURE" COMVALQRI(A, N, EM, RE, IM);         
         "CODE" 34190;      
        
         EQILBR(A, N, EM, D, INT0); TFMREAHES(A, N, EM, INT);       
         COMEIGVAL:= COMVALQRI(A, N, EM, RE, IM)
     "END" COMEIGVAL;       
         "EOP"    
1SECTION 3.3.1.2.2            (JULY 1974)                        PAGE 16      
        
        
0"CODE" 34194;    
     "COMMENT" MCA 2424;    
     "INTEGER" "PROCEDURE" COMEIG1(A, N, EM, RE, IM, VEC);
     "VALUE" N; "INTEGER" N;
     "ARRAY" A, EM, RE, IM, VEC;      
     "BEGIN" "INTEGER" I, J, K, PJ, ITT;        
         "REAL" X, Y, MAX, NEPS;      
         "ARRAY" AB[1:N,1:N], D, U, V[1:N];     
         "INTEGER" "ARRAY" INT, INT0[1:N];      
        
         "PROCEDURE" TRANSFER;        
         "BEGIN" "INTEGER" I, J;      
             "FOR" I:= 1 "STEP" 1 "UNTIL" N "DO"
             "FOR" J:= ("IF" I = 1 "THEN" 1 "ELSE" I - 1) "STEP" 1  
             "UNTIL" N "DO" AB[I,J]:= A[I,J]    
         "END" TRANSFER;    
        
         "PROCEDURE" EQILBR(A, N, EM, D, INT); "CODE" 34173;        
         "PROCEDURE" TFMREAHES(A, N, EM, INT); "CODE" 34170;        
         "PROCEDURE" BAKREAHES2(A, N, N1, N2, INT, VEC); "CODE" 34172;        
         "PROCEDURE" BAKLBR(N, N1, N2, D, INT, VEC); "CODE" 34174;  
         "PROCEDURE" REAVECHES(A, N, LAMBDA, EM, V); "CODE" 34181;  
         "PROCEDURE" COMSCL(A, N, N1, N2, IM); "CODE" 34193;        
         "INTEGER" "PROCEDURE" COMVALQRI(A, N, EM, RE, IM);         
         "CODE" 34190;      
         "PROCEDURE" COMVECHES(A, N, LAMBDA, MU, EM, U, V);         
         "CODE" 34191;      
        
         EQILBR(A, N, EM, D, INT0); TFMREAHES(A, N, EM, INT); TRANSFER;       
         K:= COMEIG1:= COMVALQRI(AB, N, EM, RE, IM);      
         NEPS:= EM[0] * EM[1]; MAX:= 0; ITT:= 0;
         "FOR" I:= K + 1 "STEP" 1 "UNTIL" N "DO"
         "BEGIN" X:= RE[I]; Y:= IM[I]; PJ:= 0;  
          AGAIN: "FOR" J:= K + 1 "STEP" 1 "UNTIL" I - 1 "DO"        
             "BEGIN" "IF" ((X - RE[J]) ** 2 +   
                 (Y - IM[J]) ** 2 <= NEPS ** 2) "THEN"    
                 "BEGIN" "IF" PJ = J "THEN" NEPS:= EM[2] * EM[1]    
                     "ELSE" PJ:= J; X:= X + 2 * NEPS; "GOTO" AGAIN  
                 "END"      
             "END";         
             RE[I]:= X; TRANSFER; "IF" Y ^= 0 "THEN"      
             "BEGIN" COMVECHES(AB, N, RE[I], IM[I], EM, U, V);      
                 "FOR" J:= 1 "STEP" 1 "UNTIL" N "DO" VEC[J,I]:= U[J];         
                 I:= I + 1; RE[I]:= X 
             "END"
             "ELSE" REAVECHES(AB, N, X, EM, V); 
             "FOR" J:= 1 "STEP" 1 "UNTIL" N "DO" VEC[J,I]:= V[J];   
             "IF" EM[7] > MAX "THEN" MAX:= EM[7];         
             ITT:= "IF" ITT > EM[9] "THEN" ITT "ELSE" EM[9]         
         "END";   
         EM[7]:= MAX; EM[9]:= ITT; BAKREAHES2(A, N, K + 1, N, INT, VEC);      
         BAKLBR(N, K + 1, N, D, INT0, VEC); COMSCL(VEC, N, K + 1, N, IM)      
     "END" COMEIG1;         
         "EOP"    
1SECTION 3.3.2.1              (JULY 1974)                         PAGE 1      
        
        
        
 AUTHOR   : C.G. VAN DER LAAN.        
        
        
 CONTRIBUTORS : H.FIOLET, C.G. VAN DER LAAN.    
        
        
 INSTITUTE : MATHEMATICAL CENTRE.     
        
        
 RECEIVED: 730917.
        
        
 BRIEF DESCRIPTION:         
        
     THIS SECTION CONTAINS FOUR PROCEDURES FOR CALCULATING THE      
     EIGENVALUES OR THE EIGENVALUES AND EIGENVECTORS OF COMPLEX     
     HERMITIAN MATRICES.    
     EIGVALHRM CALCULATES THE EIGENVALUES OF A HERMITIAN MATRIX.    
     EIGHRM CALCULATES THE EIGENVALUES AND EIGENVECTORS OF A HERMITIAN        
     MATRIX.      
     QRIVALHRM CALCULATES THE EIGENVALUES OF A HERMITIAN MATRIX.    
     QRIHRM CALCULATES THE EIGENVALUES AND EIGENVECTORS OF A HERMITIAN        
     MATRIX.      
     WHEN A SMALL NUMBER OF EIGENVALUES OR EIGENVALUES AND EIGENVECTORS       
     IS REQUIRED, THE USE OF EIGVALHRM OR EIGHRM IS RECOMMANDED; WHEN         
     MORE THAN, SAY, 25 PERCENT OF THE EIGENSYSTEM IS REQUIRED. THE 
     PROCEDURES QRIVALHRM OR QRIHRM ARE TO BE USED.       
        
        
1SECTION 3.3.2.1              (JULY 1974)                         PAGE 2      
        
        
        
 SUBSECTION: EIGVALHRM.     
        
        
 CALLING SEQUENCE :         
        
     THE HEADING OF THE PROCEDURE READS :       
     "PROCEDURE" EIGVALHRM(A, N, NUMVAL, VAL, EM); "VALUE" N, NUMVAL;         
     "INTEGER" N, NUMVAL; "ARRAY" A, VAL, EM;   
        
     THE MEANING OF THE FORMAL PARAMETERS IS :  
     A:         <ARRAY IDENTIFIER>;   
                "ARRAY" A[1:N,1:N];   
                ENTRY: THE  REAL  PART  OF  THE  UPPER  TRIANGLE OF  THE      
                       HERMITIAN   MATRIX  MUST  BE  GIVEN IN THE  UPPER      
                       TRIANGULAR PART OF A (THE ELEMENTS A[I,J], I<=J);      
                       THE  IMAGINARY  PART OF THE STRICT LOWER TRIANGLE      
                       OF  THE   HERMITIAN  MATRIX  MUST BE GIVEN IN THE      
                       STRICT LOWER PART OF A (THE ELEMENTS A[I,J],I>J);      
                THE ELEMENTS AF A ARE ALTERED;  
     N:         <ARITHMETIC EXPRESSION>;        
                THE ORDER OF THE GIVEN MATRIX;  
     NUMVAL:    <ARITHMETIC EXPRESSION>;        
                EIGVALHRM  CALCULATES  THE LARGEST NUMVAL EIGENVALUES OF      
                THE HERMITIAN MATRIX; 
     VAL:       <ARRAY IDENTIFIER>;   
                "ARRAY" VAL[1:NUMVAL];
                EXIT:       
                IN ARRAY VAL THE LARGEST NUMVAL EIGENVALUES ARE     
                DELIVERED IN MONOTONICALLY NONINCREASING ORDER;     
     EM:        <ARRAY IDENTIFIER>;   
                "ARRAY" EM[0:3];      
                ENTRY:      
                EM[0]: THE MACHINE PRECISION;   
                EM[2]: THE RELATIVE TOLERANCE FOR THE EIGENVALUES;  
                       MORE PRECISELY: THE TOLERANCE FOR EACH EIGENVALUE      
                       LAMBDA, IS ABS(LAMBDA)*EM[2]+EM[1]*EM[0];    
                EXIT:       
                EM[1]: AN ESTIMATE OF A NORM OF THE ORIGINAL MATRIX;
                EM[3]: THE NUMBER OF ITERATIONS PERFORMED.
        
        
 PROCEDURES USED: 
        
     HSHHRMTRIVAL = CP34364,
     VALSYMTRI    = CP34151.
        
        
 REQUIRED CENTRAL MEMORY:   
     TWO AUXILIARY ARRAYS OF ORDER N AND N - 1 RESPECTIVELY ARE DECLARED      
        
        
 RUNNING TIME: PROPORTIONAL TO N CUBED.         
1SECTION 3.3.2.1              (JULY 1974)                         PAGE 3      
        
        
        
        
        
 LANGUAGE: ALGOL 60.        
        
        
 METHOD AND PERFORMANCE: SEE QRIHRM (THIS SECTION).       
        
        
 EXAMPLE OF USE:  SEE EIGHRM (THIS SECTION).    
        
        
 SUBSECTION: EIGHRM.        
        
        
 CALLING SEQUENCE :         
        
     THE HEADING OF THE PROCEDURE READS :       
     "PROCEDURE" EIGHRM(A, N, NUMVAL, VAL, VECR, VECI, EM);         
     "VALUE" N, NUMVAL; "INTEGER" N, NUMVAL;    
     "ARRAY" A, VAL, VECR, VECI, EM;  
        
     THE MEANING OF THE FORMAL PARAMETERS IS :  
     A:         <ARRAY IDENTIFIER>;   
                "ARRAY" A[1:N,1:N];   
                ENTRY: THE  REAL  PART  OF  THE  UPPER  TRIANGLE OF  THE      
                       HERMITIAN   MATRIX  MUST  BE  GIVEN IN THE  UPPER      
                       TRIANGULAR PART OF A (THE ELEMENTS A[I,J], I<=J);      
                       THE  IMAGINARY  PART OF THE STRICT LOWER TRIANGLE      
                       OF  THE   HERMITIAN  MATRIX  MUST BE GIVEN IN THE      
                       STRICT LOWER PART OF A (THE ELEMENTS A[I,J],I>J);      
                THE ELEMENTS AF A ARE ALTERED;  
     N:         <ARITHMETIC EXPRESSION>;        
                THE ORDER OF THE GIVEN MATRIX;  
     NUMVAL:    <ARITHMETIC EXPRESSION>;        
                EIGHRM  CALCULATES THE LARGEST NUMVAL EIGENVALUES OF THE      
                HERMITIAN MATRIX;     
     VAL:       <ARRAY IDENTIFIER>;   
                "ARRAY" VAL[1:NUMVAL];
                EXIT:       
                IN ARRAY VAL THE LARGEST NUMVAL EIGENVALUES ARE     
                DELIVERED IN MONOTONICALLY NONINCREASING ORDER;     
     VECR,VECI: <ARRAY IDENTIFIER>;   
                "ARRAY" VECR,VECI[1:N,1:NUMVAL];
                EXIT:       
                THE CALCULATED EIGENVECTORS;    
                THE  COMPLEX EIGENVECTOR WITH REAL PART VECR[1:N,I]  AND      
                IMAGINARY PART VECI[1:N,I] CORRESPONDS TO THE EIGENVALUE      
                VAL[I], I=1,...,NUMVAL;         
1SECTION 3.3.2.1              (JULY 1974)                         PAGE 4      
        
        
        
     EM:        <ARRAY IDENTIFIER>;   
                "ARRAY" EM[0:9];      
                ENTRY:      
                EM[0]: THE MACHINE PRECISION;   
                EM[2]: THE RELATIVE TOLERANCE FOR THE EIGENVALUES;  
                       MORE PRECISELY: THE TOLERANCE FOR EACH EIGENVALUE      
                       LAMBDA, IS ABS(LAMBDA)*EM[2]+EM[1]*EM[0];    
                EM[4]: THE ORTHOGONALIZATION PARAMETER (E.G. .01);  
                EM[6]: THE TOLERANCE FOR THE EIGENVECTORS;
                EM[8]: THE  MAXIMUM NUMBER OF INVERSE ITERATIONS ALLOWED      
                       FOR THE CALCULATION OF EACH EIGENVECTOR;     
                EXIT:       
                EM[1]: AN ESTIMATE OF A NORM OF THE ORIGINAL MATRIX;
                EM[3]: THE NUMBER OF ITERATIONS PERFORMED;
                EM[5]: THE  NUMBER  OF EIGENVECTORS INVOLVED IN THE LAST      
                       GRAM-SCHMIDT ORTHOGONALIZATION;    
                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      
                       EIGENVECTOR; IF, HOWEVER,  FOR  SOME   CALCULATED      
                       EIGENVECTOR, THE  EUCLIDEAN NORM OF THE  RESIDUES      
                       REMAINS     GREATER  THAN   EM[1]*EM[6],     THEN      
                       EM[9]:=EM[8]+1.
        
        
 PROCEDURES USED: 
        
     HSHHRMTRI = CP34363,   
     VALSYMTRI = CP34151,   
     VECSYMTRI = CP34152,   
     BAKHRMTRI = CP34365.   
        
        
 REQUIRED CENTRAL MEMORY:   
     THREE AUXILIARY ARRAYS OF ORDER N-1 AND TWO OF ORDER N ARE DECLARED      
        
        
 RUNNING TIME: PROPORTIONAL TO N CUBED.         
        
        
 LANGUAGE: ALGOL 60.        
        
        
 METHOD AND PERFORMANCE: SEE QRIHRM (THIS SECTION).       
        
        
1SECTION 3.3.2.1              (JULY 1974)                         PAGE 5      
        
        
        
 EXAMPLE OF USE:  
        
     LET  EIGHRM  CALCULATE THE LARGEST EIGENVALUE AND THE CORRESPONDING      
     EIGENVECTOR OF THE FOLLOWING MATRIX:       
     (SEE GREGORY AND KARNEY, CHAPTER 6, EXAMPLE 6.6)     
          3    1    0   +2I 
          1    3   -2I   0  
          0   +2I   1    1  
         -2I   0    1    1  
     THE EIGENVECTORS ARE NORMALIZED BY THE PROCEDURE SCLCOM (SEE   
     SECTION 1.2.11. ).     
        
     "BEGIN"      
     "COMMENT" GREGORY AND KARNEY,CHAPTER 6, EXAMPLE 6.6; 
     "PROCEDURE" SCLCOM(AR,AI,N,N1,N2);"CODE" 34360;      
     "PROCEDURE" INIMAT(LR,UR,LC,UC,A,X);"CODE" 31011;    
     "PROCEDURE" EIGHRM(A,N,NUMVAL,VAL,VECR,VECI,EM);"CODE" 34369;  
     "REAL" "ARRAY" A[1:4,1:4],VAL[1:1],VECR,VECI[1:4,1:1],EM[0:9]; 
     "INTEGER" I; 
     INIMAT(1,4,1,4,A,0);   
     A[1,1]:=A[2,2]:=3;     
     A[1,2]:=A[3,3]:=A[3,4]:=A[4,4]:=1;         
     A[3,2]:=2;A[4,1]:=-2;  
     EM[0]:=5"-14;EM[2]:="-12;        
     EM[4]:=.01;EM[6]:="-10;EM[8]:=5; 
     EIGHRM(A,4,1,VAL,VECR,VECI,EM);  
     SCLCOM(VECR,VECI,4,1,1);         
     OUTPUT(61,"(""("LARGEST EIGENVALUE: ")",N/")",VAL[1]);         
     OUTPUT(61,"(""("CORRESPONDING EIGENVECTOR:")",/")"); 
     "FOR" I:=1,2,3,4 "DO"  
     OUTPUT(61,"("+D.D,+D.DDD,"("*I")",/")",VECR[I,1],VECI[I,1]);   
     "FOR" I:=1,3,5,7,9 "DO"
     OUTPUT(61,"("/,"("EM[")",D,"("]: ")",+D.DDD"+DD")",I,EM[I]);   
     "END"        
        
     DELIVERS:    
        
     LARGEST EIGENVALUE: +4.8284271247462"000   
     CORRESPONDING EIGENVECTOR:       
     +1.0+0.000*I 
     +1.0+0.000*I 
     -0.0+0.414*I 
     +0.0-0.414*I 
        
     EM[1]: +6.000"+00      
     EM[3]: +1.800"+01      
     EM[5]: +1.000"+00      
     EM[7]: +5.303"-14      
     EM[9]: +1.000"+00      
        
        
1SECTION 3.3.2.1              (JULY 1974)                         PAGE 6      
        
        
        
 SUBSECTION: QRIVALHRM.     
        
        
 CALLING SEQUENCE :         
        
     THE HEADING OF THE PROCEDURE READS :       
     "INTEGER" "PROCEDURE" QRIVALHRM(A, N, VAL, EM); "VALUE" N;     
     "INTEGER" N; "ARRAY" A, VAL, EM; 
        
     THE MEANING OF THE FORMAL PARAMETERS IS:   
     A:         <ARRAY IDENTIFIER>;   
                "ARRAY" A[1:N,1:N];   
                ENTRY: THE  REAL  PART  OF  THE  UPPER  TRIANGLE OF  THE      
                       HERMITIAN   MATRIX  MUST  BE  GIVEN IN THE  UPPER      
                       TRIANGULAR PART OF A (THE ELEMENTS A[I,J], I<=J);      
                       THE  IMAGINARY  PART OF THE STRICT LOWER TRIANGLE      
                       OF  THE   HERMITIAN  MATRIX  MUST BE GIVEN IN THE      
                       STRICT LOWER PART OF A (THE ELEMENTS A[I,J],I>J);      
                THE ELEMENTS AF A ARE ALTERED;  
     N:         <ARITHMETIC EXPRESSION>;        
                THE ORDER OF THE GIVEN MATRIX;  
     VAL:       <ARRAY IDENTIFIER>;   
                "ARRAY" VAL[1:N];     
                EXIT:       
                THE CALCULATED EIGENVALUES;     
     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]: AN ESTIMATE OF A NORM OF THE ORIGINAL MATRIX;
                EM[3]: THE  MAXIMUM   ABSOLUTE VALUE OF THE   CODIAGONAL      
                       ELEMENTS NEGLECTED;      
                EM[5]: THE NUMBER OF ITERATIONS PERFORMED;
                       EM[5]:= EM[4]+1 IN THE CASE QRIVALHRM^=0;    
        
     QRIVALHRM:=0, PROVIDED  THE QR ITERATION IS COMPLETED  WITHIN EM[4]      
     ITERATIONS; OTHERWISE, QRIVALHRM:=THE NUMBER OF EIGENVALUES, K, NOT      
     CALCULATED AND ONLY THE LAST N-K ELEMENTS OF VAL ARE APPROXIMATE         
     EIGENVALUES OF THE ORIGINAL HERMITEAN MATRIX.        
        
        
 PROCEDURES USED: 
        
     HSHHRMTRIVAL = CP34364,
     QRIVALSYMTRI = CP34160.
        
        
 REQUIRED CENTRAL MEMORY:   
     TWO  AUXILIARY ARRAYS OF ORDER N  ARE DECLARED.      
        
 RUNNING TIME: PROPORTIONAL TO N CUBED.         
        
        
1SECTION 3.3.2.1              (JULY 1974)                         PAGE 7      
        
        
        
 LANGUAGE: ALGOL 60.        
        
 METHOD AND PERFORMANCE: SEE QRIHRM (THIS SECTION).       
        
        
 EXAMPLE OF USE: SEE QRIHRM (THIS SECTION).     
        
        
 SUBSECTION: QRIHRM.        
        
        
 CALLING SEQUENCE :         
        
     THE HEADING OF THE PROCEDURE READS :       
     "INTEGER" "PROCEDURE" QRIHRM(A, N, VAL, VR, VI, EM); "VALUE" N;
     "INTEGER" N; "ARRAY" A, VAL, VR, VI, EM;   
        
     THE MEANING OF THE FORMAL PARAMETERS IS:   
     A:         <ARRAY IDENTIFIER>;   
                "ARRAY" A[1:N,1:N];   
                ENTRY: THE  REAL  PART  OF  THE  UPPER  TRIANGLE OF  THE      
                       HERMITIAN   MATRIX  MUST  BE  GIVEN IN THE  UPPER      
                       TRIANGULAR PART OF A (THE ELEMENTS A[I,J], I<=J);      
                       THE  IMAGINARY  PART OF THE STRICT LOWER TRIANGLE      
                       OF  THE   HERMITIAN  MATRIX  MUST BE GIVEN IN THE      
                       STRICT LOWER PART OF A (THE ELEMENTS A[I,J],I>J);      
                THE ELEMENTS AF A ARE ALTERED;  
     N:         <ARITHMETIC EXPRESSION>;        
                THE ORDER OF THE GIVEN MATRIX;  
     VAL:       <ARRAY IDENTIFIER>;   
                "ARRAY" VAL[1:N];     
                EXIT:       
                THE CALCULATED EIGENVALUES;     
     VR,VI:     <ARRAY IDENTIFIER>;   
                "ARRAY" VR,VI[1:N,1:N];         
                EXIT:       
                THE CALCULATED EIGENVECTORS;    
                THE  COMPLEX  EIGENVECTOR WITH REAL PART VR[1:N,I]   AND      
                IMAGINARY  PART VI[1:N,I] CORRESPONDS TO THE  EIGENVALUE      
                VAL[I], I=1,...,N;    
     EM:        <ARRAY IDENTIFIER>;   
                "ARRAY" EM[0:5];      
                ENTRY:      
                EM[0]: THE MACHINE PRECISION;   
                EM[2]: THE RELATIVE TOLERANCE FOR THE QR ITERATION; 
                       (E.G. THE MACHINE PRECISION);      
                EM[4]: THE MAXIMUM ALLOWED NUMBER OF ITERATIONS;    
                       (E.G. 10 * N); 
                EXIT:       
                EM[1]: AN ESTIMATE OF A NORM OF THE ORIGINAL MATRIX;
                EM[3]: THE  MAXIMUM   ABSOLUTE VALUE OF THE   CODIAGONAL      
                       ELEMENTS NEGLECTED;      
                EM[5]: THE NUMBER OF ITERATIONS PERFORMED;
                       EM[5]:=EM[4]+1 IN THE CASE QRIHRM^=0;        
1SECTION 3.3.2.1              (JULY 1974)                         PAGE 8      
        
        
        
        
     QRIHRM:=0, PROVIDED  THE   PROCESS   IS   COMPLETED   WITHIN  EM[4]      
     ITERATIONS; OTHERWISE, QRIHRM:= THE  NUMBER OF EIGENVALUES, K, NOT       
     CALCULATED AND ONLY THE LAST N-K ELEMENTS OF  VAL ARE APPROXIMATE        
     EIGENVALUES AND THE COLUMNS OF THE ARRAYS VR,VI[1:N,N-K:N] ARE 
     APPROXIMATE EIGENVECTORS OF THE ORIGINAL HERMITEAN MATRIX .    
        
        
 PROCEDURES USED: 
        
     HSHHRMTRI = CP34363,   
     QRISYMTRI = CP34161,   
     BAKHRMTRI = CP34365.   
        
        
 REQUIRED CENTRAL MEMORY:   
     TWO AUXILIARY ARRAYS OF ORDER N - 1 AND TWO OF ORDER N ARE DECLARED      
        
        
 RUNNING TIME: PROPORTIONAL TO N CUBED.         
        
        
 LANGUAGE: ALGOL 60.        
        
        
 THE FOLLOWING HOLDS FOR THE FOUR PROCEDURES OF THIS SECTION:       
        
        
 METHOD AND PERFORMANCE:    
        
     FOR THE TRANSFORMATION OF THE GIVEN HERMITIAN MATRIX INTO A REAL         
     SYMMETRIC TRIDIAGONAL MATRIX, AND FOR THE CORRESPONDING BACK   
     TRANSFORMATION, PROCEDURES OF SECTION  3.2.1.2.2.1. ARE USED.  
     FOR THE CALCULATION OF THE EIGENVALUES AND EIGENVECTORS OF THE 
     RESULTING SYMMETRIC TRIDIAGONAL MATRIX, PROCEDURES OF SECTION  
     3.3.1.1.1. ARE USED.   
        
        
 EXAMPLE OF USE:  
        
     QRIHRM CALCULATES THE EIGENVALUES AND EIGENVECTORS OF THE      
     FOLLOWING MATRIX:      
     (SEE GREGORY AND KARNEY, CHAPTER 6, EXAMPLE 6.6)     
          3    1    0   +2I 
          1    3   -2I   0  
          0   +2I   1    1  
         -2I   0    1    1  
     THE EIGENVECTORS ARE NORMALIZED BY THE PROCEDURE SCLCOM (SEE   
     SECTION 1.2.11. ).     
     ONLY THE EIGENVECTORS CORRESPONDING TO VAL[2] AND VAL[3] ARE   
     PRINTED BY THE FOLLOWING PROGRAM:
        
1SECTION 3.3.2.1              (JULY 1974)                         PAGE 9      
        
        
     "BEGIN"      
     "COMMENT" GREGORY AND KARNEY,CHAPTER 6, EXAMPLE 6.6; 
     "INTEGER" "PROCEDURE" QRIHRM(A,N,VAL,VR,VI,EM);"CODE" 34371;   
     "PROCEDURE" INIMAT(LR,UR,LC,UC,A,X);"CODE" 31011;    
     "PROCEDURE" SCLCOM(AR,AI,N,N1,N2);"CODE" 34360;      
     "REAL" "ARRAY" A,VR,VI[1:4,1:4],VAL[1:4],EM[0:5];"INTEGER" I;  
     INIMAT(1,4,1,4,A,0);   
     A[1,1]:=A[2,2]:=3;     
     A[3,2]:=2;A[4,1]:=-2;  
     A[1,2]:=A[3,3]:=A[3,4]:=A[4,4]:=1;         
     EM[0]:=EM[2]:=5"-14;EM[4]:=20;   
     OUTPUT(61,"(""("QRIHRM: ")",D/")",QRIHRM(A,4,VAL,VR,VI,EM));   
     SCLCOM(VR,VI,4,2,3);   
     OUTPUT(61,"(""("EIGENVALUES: ")"")");      
     "FOR" I:=1,2,3,4 "DO" OUTPUT(61,"("/,"("VAL[")",D,"("]: ")",   
                                     +D.3DBB")",I,VAL[I]);
     OUTPUT(61,"("/,"("EIGENVECTORS CORRESPONDING TO")",/,
                  "(" VAL[2] , VAL[3] ")",/")");
     "FOR" I:=1,2,3,4 "DO"  
     OUTPUT(61,"("+D,+D,"("*I  ,  ")",+D,+D,"("*I")",/")",
               VR[I,2],VI[I,2],VR[I,3],VI[I,3]);
     "FOR" I:=1,3,5 "DO"    
      OUTPUT(61,"("/,"("EM[")",D,"("]:  ")",+D.DDD"+DD")",I,EM[I])  
     "END"        
        
   OUTPUT:        
     QRIHRM: 0    
     EIGENVALUES: 
     VAL[1]: +4.828         
     VAL[2]: +4.000         
     VAL[3]: -0.000         
     VAL[4]: -0.828         
     EIGENVECTORS CORRESPONDING TO    
      VAL[2] , VAL[3]       
     +1+0*I  ,  +0-1*I      
     -1+0*I  ,  +0+1*I      
     +0-1*I  ,  +1+0*I      
     +0-1*I  ,  +1+0*I      
        
     EM[1]:  +6.000"+00     
     EM[3]:  +3.804"-22     
     EM[5]:  +6.000"+00     
        
 SOURCE TEXT(S) : 
0"CODE" 34368;    
     "PROCEDURE" EIGVALHRM(A, N, NUMVAL, VAL, EM); "VALUE" N, NUMVAL;         
     "INTEGER" N, NUMVAL; "ARRAY" A, VAL, EM;   
     "BEGIN" "ARRAY" D[1:N], BB[1:N - 1];       
         "PROCEDURE" HSHHRMTRIVAL(A,N,D,BB,EM);"CODE" 34364;        
         "PROCEDURE" VALSYMTRI(D,BB,N,N1,N2,VAL,EM);"CODE" 34151;   
         HSHHRMTRIVAL(A, N, D, BB, EM);         
         VALSYMTRI(D, BB, N, 1, NUMVAL, VAL, EM)
     "END" EIGVALHRM;       
         "EOP"    
1SECTION 3.3.2.1              (JULY 1974)                        PAGE 10      
        
        
        
0"CODE" 34369;    
     "PROCEDURE" EIGHRM(A, N, NUMVAL, VAL, VECR, VECI, EM);         
     "VALUE" N, NUMVAL; "INTEGER" N, NUMVAL;    
     "ARRAY" A, VAL, VECR, VECI, EM;  
     "BEGIN" "ARRAY" BB, TR, TI[1:N - 1], D, B[1:N];      
         "PROCEDURE" HSHHRMTRI(A,N,D,B,BB,EM,TR,TI);"CODE" 34363;   
         "PROCEDURE" VALSYMTRI(D,BB,N,N1,N2,VAL,EM);"CODE" 34151;   
         "PROCEDURE" VECSYMTRI(D,B,N,N1,N2,VAL,VEC,EM);"CODE" 34152;
         "PROCEDURE" BAKHRMTRI(A,N,N1,N2,VECR,VECI,TR,TI);"CODE" 34365;       
         HSHHRMTRI(A, N, D, B, BB, EM, TR, TI); 
         VALSYMTRI(D, BB, N, 1, NUMVAL, VAL, EM); B[N]:= 0;         
         VECSYMTRI(D, B, N, 1, NUMVAL, VAL, VECR, EM);    
         BAKHRMTRI(A, N, 1, NUMVAL, VECR, VECI, TR, TI)   
     "END" EIGHRM;
         "EOP"    
0"CODE" 34370;    
     "INTEGER" "PROCEDURE" QRIVALHRM(A, N, VAL, EM); "VALUE" N;     
     "INTEGER" N; "ARRAY" A, VAL, EM; 
     "BEGIN" "ARRAY" B,BB[1:N];       
         "INTEGER" I;       
         "PROCEDURE" HSHHRMTRIVAL(A,N,D,BB,EM);"CODE" 34364;        
         "INTEGER" "PROCEDURE" QRIVALSYMTRI(D, BB, N, EM);"CODE" 34160;       
         HSHHRMTRIVAL(A, N, VAL, BB, EM); B[N]:=BB[N]:= 0;
         "FOR" I:=1 "STEP" 1 "UNTIL" N-1 "DO" B[I]:=SQRT(BB[I]);    
         QRIVALHRM:=QRIVALSYMTRI(VAL, BB, N, EM)
     "END" QRIVALHRM;       
         "EOP"    
0"CODE" 34371;    
     "INTEGER" "PROCEDURE" QRIHRM(A, N, VAL, VR, VI, EM); "VALUE" N;
     "INTEGER" N; "ARRAY" A, VAL, VR, VI, EM;   
     "BEGIN" "INTEGER" I, J;
         "ARRAY" B, BB[1:N], TR, TI[1:N - 1];   
         "PROCEDURE" HSHHRMTRI(A,N,D,B,BB,EM,TR,TI);"CODE" 34363;   
         "INTEGER" "PROCEDURE" QRISYMTRI(A,N,D,B,BB,EM);"CODE" 34161;         
         "PROCEDURE" BAKHRMTRI(A,N,N1,N2,VECR,VECI,TR,TI);"CODE" 34365;       
         HSHHRMTRI(A, N, VAL, B, BB, EM, TR, TI);         
         "FOR" I:= 1 "STEP" 1 "UNTIL" N "DO"    
         "BEGIN" VR[I,I]:= 1;         
             "FOR" J:= I + 1 "STEP" 1 "UNTIL" N "DO" VR[I,J]:= VR[J,I]:=      
             0    
         "END";   
         B[N]:= BB[N]:= 0;  
         I:= QRIHRM:= QRISYMTRI(VR, N, VAL, B, BB, EM);   
         BAKHRMTRI(A, N, I+1, N, VR, VI, TR, TI);         
     "END" QRIHRM;
         "EOP"