1SECTION 3.3.2.2.1            (JULY 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 VALQRICOM AND QRICOM.     
     VALQRICOM CALCULATES THE EIGENVALUES OF A COMPLEX UPPER-HESSENBERG       
     MATRIX WITH A REAL SUBDIAGONAL.  
     QRICOM CALCULATES THE EIGENVECTORS AS WELL.
        
        
 KEYWORDS :       
        
     EIGENVALUES, 
     EIGENVECTORS,
     COMPLEX UPPER-HESSENBERG MATRIX, 
     QR-ITERATION.
        
        
1SECTION 3.3.2.2.1            (JULY 1974)                         PAGE 2      
        
        
        
 SUBSECTION: VALQRICOM.     
        
        
 CALLING SEQUENCE:
        
     THE HEADING OF THE PROCEDURE READS:        
     "INTEGER" "PROCEDURE" VALQRICOM(A1, A2, B, N, EM, VAL1, VAL2); 
     "VALUE" N; "INTEGER" N; "ARRAY" A1, A2, B, EM, VAL1, VAL2;     
        
     THE MEANING OF THE FORMAL PARAMETERS IS:   
     A1,A2:     <ARRAY IDENTIFIER>;   
                "ARRAY" A1,A2[1:N,1:N];         
                ENTRY:      
                THE REAL PART AND THE IMAGINARY PART OF THE UPPER   
                TRIANGLE OF THE UPPER-HESSENBERG MATRIX MUST BE GIVEN IN      
                THE CORRESPONDING PARTS OF THE ARRAYS A1 AND A2;    
                THE ELEMENTS IN THE UPPER TRIANGLE OF THE ARRAYS A1 AND       
                A2 ARE ALTERED;       
     B:         <ARRAY IDENTIFIER>;   
                "ARRAY" B[1:N-1];     
                ENTRY:      
                THE REAL SUBDIAGONAL OF THE UPPER-HESSENBERG MATRIX;
                THE ELEMENTS OF THE ARRAY B 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[1]: AN ESTIMATE OF THE NORM OF THE UPPER-HESSENBERG        
                       MATRIX (E.G. THE SUM OF THE INFINITY NORMS OF THE      
                       REAL AND IMAGINARY PARTS OF THE MATRIX);     
                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[3]: THE MAXIMUM ABSOLUTE VALUE OF THE SUBDIAGONAL
                       ELEMENTS NEGLECTED;      
                EM[5]: THE NUMBER OF ITERATIONS PERFORMED;
                       EM[5]:=EM[4]+1 IN THE CASE VALQRICOM^=0;     
     VAL1,VAL2: <ARRAY IDENTIFIER>;   
                "ARRAY" VAL1,VAL2[1:N];         
                EXIT:       
                THE REAL PART AND THE IMAGINARY PART OF THE CALCULATED        
                EIGENVALUES ARE DELIVERED IN THE ARRAYS VAL1 AND VAL2,        
                RESPECTIVELY;         
        
     VALQRICOM:=0, PROVIDED THE PROCESS IS COMPLETED WITHIN EM[4]   
     ITERATIONS; OTHERWISE, VALQRICOM:= THE NUMBER, K, OF EIGENVALUES         
     NOT CALCULATED AND ONLY THE LAST N-K ELEMENTS OF THE ARRAYS VAL1         
     AND VAL2 ARE APPROXIMATE EIGENVALUES OF THE UPPER-HESSENBERG   
     MATRIX.      
        
        
1SECTION 3.3.2.2.1            (JULY 1974)                         PAGE 3      
        
        
        
 PROCEDURES USED: 
        
     COMKWD    = CP34345,   
     ROTCOMROW = CP34358,   
     ROTCOMCOL = CP34357,   
     COMCOLCST = CP34352.   
        
 RUNNING TIME: PROPORTIONAL TO N**2 * NUMBER OF ITERATIONS.         
        
 LANGUAGE: ALGOL 60.        
        
 METHOD AND PERFORMANCE: SEE QRICOM (THIS SECTION).       
        
 EXAMPLE OF USE:  
     AS A FORMAL TEST OF THE PROCEDURE VALQRICOM THE ZEROS OF THE   
     POLYNOMIAL  X**4 + (4+2*I)* X**3 + (5+6*I)* X**2 + (2+6*I)* X + 2*I      
     ARE OBTAINED BY MEANS OF THE CALCULATION OF THE EIGENVALUES OF THE       
     FOLLOWING COMPANION MATRIX:      
     (SEE WILKINSON AND REINSCH, 1971, CONTRIBUTION II/15)
     -4-2*I  -5-6*I  -2-6*I    -2*I   
       1       0       0       0      
       0       1       0       0      
       0       0       1       0      
        
     "BEGIN"      
     "REAL" "ARRAY" A1,A2[1:4,1:4],B[1:3],EM[0:5],VAL1,VAL2[1:4];   
     "INTEGER" I; 
     "PROCEDURE" INIMAT(LR,UR,LC,UC,A,X);"CODE" 31011;    
     "INTEGER" "PROCEDURE" VALQRICOM(A1,A2,B,N,EM,VAL1,VAL2);       
     "CODE" 34372;
     INIMAT(1,4,1,4,A1,0);INIMAT(1,4,1,4,A2,0); 
     A1[1,1]:=-4;A1[1,2]:=-5;A1[1,3]:=A2[1,1]:=A2[1,4]:=-2;         
     A2[1,2]:=A2[1,3]:=-6;  
     B[1]:=B[2]:=B[3]:=1;   
     EM[0]:=5"-14;EM[1]:=27;EM[2]:="-12;EM[4]:=15;        
     OUTPUT(61,"(""("VALQRICOM: ")",D/")",      
            VALQRICOM(A1,A2,B,4,EM,VAL1,VAL2)); 
     OUTPUT(61,"(""("EIGENVALUES:")",/,"("REAL PART")",14B,         
                  "("IMAGINARY PART")",/")");   
     "FOR" I:=1,2,3,4 "DO" OUTPUT(61,"("N,N/")",VAL1[I],VAL2[I]);   
     OUTPUT(61,"(""("EM[3]: ")",D.D"+DD/,"("EM[5]: ")",3D")",       
                     EM[3],EM[5])     
     "END"        
        
   OUTPUT:        
     VALQRICOM: 0 
     EIGENVALUES: 
     REAL PART              IMAGINARY PART      
     -1.0000001920467"+000  -1.0000000831019"+000         
     -9.9999980795324"-001  -9.9999991689805"-001         
     -1.0000000047492"+000  +1.6824958523393"-007         
     -9.9999999525076"-001  -1.6824956397352"-007         
     EM[3]: 3.2"-14         
     EM[5]: 010   
        
        
1SECTION 3.3.2.2.1            (JULY 1974)                         PAGE 4      
        
        
        
 SUBSECTION: QRICOM.        
        
        
 CALLING SEQUENCE:
        
     THE HEADING OF THE PROCEDURE READS:        
     "INTEGER" "PROCEDURE" QRICOM(A1,A2,B,N,EM,VAL1,VAL2,VEC1,VEC2);
     "VALUE" N; "INTEGER" N;
     "ARRAY" A1, A2, B, EM, VAL1, VAL2, VEC1, VEC2;       
        
     THE MEANING OF THE FORMAL PARAMETERS IS:   
     A1,A2:     <ARRAY IDENTIFIER>;   
                "ARRAY" A1,A2[1:N,1:N];         
                ENTRY:      
                THE REAL PART AND THE IMAGINARY PART OF THE UPPER   
                TRIANGLE OF THE UPPER-HESSENBERG MATRIX MUST BE GIVEN IN      
                THE CORRESPONDING PARTS OF THE ARRAYS A1 AND A2;    
                THE ELEMENTS IN THE UPPER TRIANGLE OF THE ARRAYS A1 AND       
                A2 ARE ALTERED;       
     B:         <ARRAY IDENTIFIER>;   
                "ARRAY" B[1:N-1];     
                ENTRY:      
                THE REAL SUBDIAGONAL OF THE UPPER-HESSENBERG MATRIX;
                THE ELEMENTS OF THE ARRAY B 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[1]: AN ESTIMATE OF THE NORM OF THE UPPER-HESSENBERG        
                       MATRIX (E.G. THE SUM OF THE INFINITY NORMS OF THE      
                       REAL AND IMAGINARY PARTS OF THE MATRIX);     
                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[3]: THE MAXIMUM ABSOLUTE VALUE OF THE SUBDIAGONAL
                       ELEMENTS NEGLECTED;      
                EM[5] THE NUMBER OF ITERATIONS PERFORMED; 
                       EM[5]:=EM[4]+1 IN THE CASE QRICOM^=0;        
     VAL1,VAL2: <ARRAY IDENTIFIER>;   
                "ARRAY" VAL1,VAL2[1:N];         
                EXIT:       
                THE REAL PART AND THE IMAGINARY PART OF THE CALCULATED        
                EIGENVALUES ARE DELIVERED IN THE ARRAYS VAL1 AND VAL2,        
                RESPECTIVELY;         
     VEC1,VEC2: <ARRAY IDENTIFIER>;   
                "ARRAY" VEC1,VEC2[1:N,1:N];     
                EXIT:       
                THE EIGENVECTORS OF THE UPPER-HESSENBERG MATRIX;    
                THE EIGENVECTOR WITH REAL PART VEC1[1:N,J] AND IMAGINARY      
                PART VEC2[1:N,J] CORRESPONDS TO THE EIGENVALUE      
                VAL1[J] + VAL2[J] * I, J=1,...,N;         
        
1SECTION 3.3.2.2.1            (JULY 1974)                         PAGE 5      
        
        
        
     QRICOM:=0, PROVIDED THE PROCESS IS COMPLETED WITHIN EM[4]      
     ITERATIONS; OTHERWISE, QRICOM:= THE NUMBER, K, OF EIGENVALUES NOT        
     CALCULATED AND ONLY THE LAST N-K ELEMENTS OF THE ARRAYS VAL1 AND         
     VAL2 ARE APPROXIMATE EIGENVALUES OF THE UPPER-HESSENBERG MATRIX,         
     AND NO USEFUL EIGENVECTORS ARE DELIVERED.  
        
        
 PROCEDURES USED: 
        
     COMKWD    = CP34345,   
     ROTCOMROW = CP34358,   
     ROTCOMCOL = CP34357,   
     COMCOLCST = CP34352,   
     COMROWCST = CP34353,   
     MATVEC    = CP34011,   
     COMMATVEC = CP34354,   
     COMDIV    = CP34342.   
        
        
 REQUIRED CENTRAL MEMORY: TWO AUXILIARY ARRAYS OF ORDER N ARE DECLARED.       
        
        
 RUNNING TIME: PROPORTIONAL TO N**2 * NUMBER OF ITERATIONS.         
        
        
 LANGUAGE: ALGOL 60.        
        
        
 THE FOLLOWING HOLDS FOR BOTH PROCEDURES:       
        
        
 METHOD AND PERFORMANCE:    
        
     THE UPPER-HESSENBERG MATRIX IS TRANSFORMED BY MEANS OF FRANCIS' QR       
     ITERATION ( FRANCIS, 1961, AND WILKINSON, 1965 ) INTO A COMPLEX
     UPPER TRIANGULAR MATRIX. THE EIGENVALUES ARE THE DIAGONAL ELEMENTS       
     OF THE LATTER MATRIX. TO CALCULATE THE EIGENVECTORS WE FIRST SOLVE       
     THE RESULTING TRIANGULAR SYSTEM OF LINEAR EQUATIONS AND        
     SUBSEQUENTLY PERFORM THE CORRESPONDING BACKTRANSFORMATION.     
        
        
1SECTION 3.3.2.2.1            (JULY 1974)                         PAGE 6      
        
        
        
 EXAMPLE OF USE:  
        
     QRICOM CALCULATES THE EIGENVALUES AND EIGENVECTORS OF THE      
     FOLLOWING MATRIX:      
     (SEE WILKINSON AND REINSCH, 1971, CONTRIBUTION II/15)
     -4-2*I  -5-6*I  -2-6*I    -2*I   
       1       0       0       0      
       0       1       0       0      
       0       0       1       0      
        
     THE EIGENVECTORS ARE NORMALIZED BY THE PROCEDURE SCLCOM.       
     (SEE SECTION 1.2.11.). 
     ONLY THE EIGENVECTOR CORRESPONDING TO VAL1[1] + VAL2[1] * I IS 
     PRINTED BY THE FOLLOWING PROGRAM.
        
     "BEGIN"      
     "REAL" "ARRAY" A1,A2,VEC1,VEC2[1:4,1:4],B[1:3],      
                    EM[0:5],VAL1,VAL2[1:4];     
     "INTEGER" I; 
     "PROCEDURE" SCLCOM(AR,AI,N,N1,N2);"CODE" 34360;      
     "INTEGER" "PROCEDURE" QRICOM(A1,A2,B,N,EM,VAL1,VAL2,VEC1,VEC2);
     "CODE" 34373;
     "PROCEDURE" INIMAT(LR,UR,LC,UC,A,X);"CODE" 31011;    
     INIMAT(1,4,1,4,A1,0);INIMAT(1,4,1,4,A2,0); 
     A1[1,1]:=-4;A1[1,2]:=-5;A1[1,3]:=A2[1,1]:=A2[1,4]:=-2;         
     A2[1,2]:=A2[1,3]:=-6;  
     B[1]:=B[2]:=B[3]:=1;   
     EM[0]:=5"-14;EM[1]:=27;EM[2]:="-12;EM[4]:=15;        
     OUTPUT(61,"(""("QRICOM: ")",D/")",         
         QRICOM(A1,A2,B,4,EM,VAL1,VAL2,VEC1,VEC2));       
     OUTPUT(61,"(""("EIGENVALUES:")",/,"("REAL PART")",14B,         
                  "("IMAGINARY PART")",/")");   
     "FOR" I:=1,2,3,4 "DO" OUTPUT(61,"("N,N/")",VAL1[I],VAL2[I]);   
     SCLCOM(VEC1,VEC2,4,1,4);         
     OUTPUT(61,"(""("FIRST EIGENVECTOR:")",/,"("REAL PART")",14B,   
                  "("IMAGINARY PART")",/")");   
     "FOR" I:=1,2,3,4 "DO" OUTPUT(61,"("N,N/")",VEC1[I,1],VEC2[I,1]);         
     OUTPUT(61,"(""("EM[3]: ")",D.D"+DD/,"("EM[5]: ")",3D")",       
                     EM[3],EM[5])     
     "END"        
1SECTION 3.3.2.2.1            (JULY 1974)                         PAGE 7      
        
        
        
        
     OUTPUT:      
        
     QRICOM: 0    
     EIGENVALUES: 
     REAL PART              IMAGINARY PART      
     -9.9999980795324"-001  -9.9999991689805"-001         
     -190000001920467"+000  -1.0000000831019"+000         
     -1.0000000047492"+000  +1.6824958523393"-007         
     -9.9999999525076"-001  -1.6824956397352"-007         
     FIRST EIGENVECTOR:     
     REAL PART              IMAGINARY PART      
     +1.0000000000000"+000  -1.7763568394003"-015         
     -5.0000004155098"-001  +5.0000009602339"-001         
     -5.4472417687634"-008  -5.0000013757436"-001         
     +2.5000014403510"-001  +2.5000006232645"-001         
     EM[3]: 3.2"-14         
     EM[5]: 010   
        
        
 REFERENCES:      
        
     DEKKER, T.J. (1968),   
     ALGOL 60 PROCEDURES IN NUMERICAL ALGEBRA, PART 1,    
     MATH. CENTRE TRACTS 22, MATHEMATISCH CENTRUM;        
        
     DEKKER, T.J. AND W.HOFFMANN (1968),        
     ALGOL 60 PROCEDURES IN NUMERICAL ALGEBRA, PART 2,    
     MATH. CENTRE TRACTS 23, MATHEMATISCH CENTRUM;        
        
     FRANCIS, J.G.F. (1961),
     THE QR-TRANSFORMATION, PART 1 AND 2,       
     COMP.J. 4, P.265-271 AND P.332-345;        
        
     RUHE, A. (1966),       
     EIGENVALUES OF A COMPLEX MATRIX BY THE QR METHOD.    
     BIT, 6, P.350-358;     
        
     WILKINSON, J.H. (1965),
     THE ALGEBRAIC EIGENVALUE PROBLEM,
     CLARENDOM PRESS, OXFORD;         
        
        
1SECTION 3.3.2.2.1            (JULY 1974)                         PAGE 8      
        
        
        
 SOURCE TEXT(S) : 
0"CODE" 34372;    
     "INTEGER" "PROCEDURE" VALQRICOM(A1, A2, B, N, EM, VAL1, VAL2); 
     "VALUE" N; "INTEGER" N; "ARRAY" A1, A2, B, EM, VAL1, VAL2;     
     "BEGIN" "INTEGER" M, NM1, I, I1, Q, Q1, MAX, COUNT;  
         "REAL" R, Z1, Z2, DD1, DD2, CC, G1, G2, K1, K2, HC, A1NN,  
         A2NN, AIJ1, AIJ2, AI1I, KAPPA, NUI, MUI1, MUI2,  
         MUIM11, MUIM12, NUIM1, TOL;  
         "PROCEDURE" COMCOLCST(L,U,J,AR,AI,XR,XI);"CODE" 34352;     
         "PROCEDURE" ROTCOMCOL(L,U,I,J,AR,AI,CR,CI,S);"CODE" 34357; 
         "PROCEDURE" ROTCOMROW(L,U,I,J,AR,AI,CR,CI,S);"CODE" 34358; 
         "PROCEDURE" COMKWD(PR,PI,QR,QI,GR,GI,KR,KI);"CODE" 34345;  
         TOL:= EM[1] * EM[2]; MAX:= EM[4]; COUNT:= 0; R:= 0;        
         M:= N; "IF" N > 1 "THEN" HC:= B[N - 1];
     IN: NM1:= N - 1;       
         "FOR" I:= N, I - 1 "WHILE" ("IF" I >= 1 "THEN" ABS(B[I]) > TOL       
         "ELSE" "FALSE") "DO" Q:= I; "IF" Q > 1 "THEN"    
         "BEGIN" "IF" ABS(B[Q - 1]) > R "THEN" R:= ABS(B[Q - 1]) "END";       
         "IF" Q = N "THEN"  
         "BEGIN" VAL1[N]:= A1[N,N]; VAL2[N]:= A2[N,N]; N:= NM1;     
             "IF" N > 1 "THEN" HC:= B[N - 1];   
         "END"    
         "ELSE"   
         "BEGIN" DD1:= A1[N,N]; DD2:= A2[N,N]; CC:= B[NM1];         
             COMKWD((A1[NM1,NM1] - DD1) / 2, (A2[NM1,NM1] - DD2)    
             / 2, CC * A1[NM1,N], CC * A2[NM1,N], G1, G2, K1,       
             K2); "IF" Q = NM1 "THEN" 
             "BEGIN" VAL1[NM1]:= G1 + DD1; VAL2[NM1]:= G2 + DD2;    
                 VAL1[N]:= K1 + DD1; VAL2[N]:= K2 + DD2;  
                 N:= N - 2; "IF" N > 1 "THEN" HC:= B[N - 1];        
             "END"
             "ELSE"         
             "BEGIN" COUNT:= COUNT + 1;         
                 "IF" COUNT > MAX "THEN" "GOTO" OUT; Z1:= K1 + DD1; 
                 Z2:= K2 + DD2;       
                 "IF" ABS(CC) > ABS(HC) "THEN" Z1:= Z1 + ABS(CC);   
                 HC:= CC / 2; I:= Q1:= Q + 1;   
                 AIJ1:= A1[Q,Q] - Z1; AIJ2:= A2[Q,Q] - Z2;
                 AI1I:= B[Q];         
                 KAPPA:= SQRT(AIJ1 ** 2 + AIJ2 ** 2 + AI1I ** 2);   
                 MUI1:= AIJ1 / KAPPA; MUI2:= AIJ2 / KAPPA;
                 NUI:= AI1I / KAPPA; A1[Q,Q]:= KAPPA;     
                 A2[Q,Q]:= 0; A1[Q1,Q1]:= A1[Q1,Q1] - Z1; 
                 A2[Q1,Q1]:= A2[Q1,Q1] - Z2;    
                 ROTCOMROW(Q1, N, Q, Q1, A1, A2, MUI1, MUI2,        
                 NUI);      
                 ROTCOMCOL(Q, Q, Q, Q1, A1, A2, MUI1, - MUI2, -     
                 NUI); A1[Q,Q]:= A1[Q,Q] + Z1;  
                 A2[Q,Q]:= A2[Q,Q] + Z2;"COMMENT"         
1SECTION 3.3.2.2.1            (JULY 1974)                         PAGE 9      
                                                                  ; 
        
        
                 "FOR" I1:= Q1 + 1 "STEP" 1 "UNTIL" N "DO"
                 "BEGIN" AIJ1:= A1[I,I]; AIJ2:= A2[I,I];  
                     AI1I:= B[I];     
                     KAPPA:= SQRT(AIJ1 ** 2 + AIJ2 ** 2 + AI1I **   
                     2); MUIM11:= MUI1; MUIM12:= MUI2;    
                     NUIM1:= NUI; MUI1:= AIJ1 / KAPPA;    
                     MUI2:= AIJ2 / KAPPA; NUI:= AI1I / KAPPA;       
                     A1[I1,I1]:= A1[I1,I1] - Z1;
                     A2[I1,I1]:= A2[I1,I1] - Z2;
                     ROTCOMROW(I1, N, I, I1, A1, A2, MUI1,
                     MUI2, NUI); A1[I,I]:= MUIM11 * KAPPA;
                     A2[I,I]:= - MUIM12 * KAPPA;
                     B[I - 1]:= NUIM1 * KAPPA;  
        
                     ROTCOMCOL(Q, I, I, I1, A1, A2, MUI1, -         
                     MUI2, - NUI); A1[I,I]:= A1[I,I] + Z1;
                     A2[I,I]:= A2[I,I] + Z2; I:= I1;      
                 "END";     
                 AIJ1:= A1[N,N]; AIJ2:= A2[N,N];
                 KAPPA:= SQRT(AIJ1 ** 2 + AIJ2 ** 2);     
                 "IF" ("IF" KAPPA < TOL "THEN" "TRUE" "ELSE" AIJ2 ** 2        
                 <= EM[0] * AIJ1 ** 2) "THEN"   
                 "BEGIN" B[NM1]:= NUI * AIJ1;   
                     A1[N,N]:= AIJ1 * MUI1 + Z1;
                     A2[N,N]:= - AIJ1 * MUI2 + Z2         
                 "END"      
                 "ELSE"     
                 "BEGIN" B[NM1]:= NUI * KAPPA; A1NN:= MUI1 * KAPPA; 
                     A2NN:= - MUI2 * KAPPA; MUI1:= AIJ1 / KAPPA;    
                     MUI2:= AIJ2 / KAPPA;       
                     COMCOLCST(Q, NM1, N, A1, A2, MUI1, MUI2);      
                     A1[N,N]:= MUI1 * A1NN - MUI2 * A2NN + Z1;      
                     A2[N,N]:= MUI1 * A2NN + MUI2 * A1NN + Z2;      
                 "END";     
             "END"
         "END";   
         "IF" N > 0 "THEN" "GOTO" IN; 
     OUT: EM[3]:= R; EM[5]:= COUNT; VALQRICOM:= N;        
     "END" VALQRICOM;       
         "EOP"    
0"CODE" 34373;    
     "INTEGER" "PROCEDURE" QRICOM(A1, A2, B, N, EM, VAL1, VAL2, VEC1,         
     VEC2); "VALUE" N; "INTEGER" N;   
     "ARRAY" A1, A2, B, EM, VAL1, VAL2, VEC1, VEC2;       
     "BEGIN" "INTEGER" M, NM1, I, I1, J, Q, Q1, MAX, COUNT;         
         "REAL" R, Z1, Z2, DD1, DD2, CC, P1, P2, T1, T2, DELTA1,    
         DELTA2, MV1, MV2, H, H1, H2, G1, G2, K1, K2, HC, 
         AIJ12, AIJ22, A1NN, A2NN, AIJ1, AIJ2, AI1I, KAPPA,         
         NUI, MUI1, MUI2, MUIM11, MUIM12, NUIM1, TOL, MACHTOL;      
         "ARRAY" TF1, TF2[1:N];"COMMENT"        
1SECTION 3.3.2.2.1            (JULY 1974)                        PAGE 10      
                                                                 ;  
        
        
         "PROCEDURE" COMKWD(PR,PI,QR,QI,GR,GI,KR,KI);"CODE" 34345;  
         "PROCEDURE" ROTCOMROW(L,U,I,J,AR,AI,CR,CI,S);"CODE" 34358; 
         "PROCEDURE" ROTCOMCOL(L,U,I,J,AR,AI,CR,CI,S);"CODE" 34357; 
         "PROCEDURE" COMCOLCST(L,U,J,AR,AI,XR,XI);"CODE" 34352;     
         "PROCEDURE" COMROWCST(L,U,I,AR,AI,XR,XI);"CODE" 34353;     
         "REAL" "PROCEDURE" MATVEC(L,U,I,A,B);"CODE" 34011;         
         "PROCEDURE" COMMATVEC(L,U,I,AR,AI,BR,BI,RR,RI);"CODE" 34354;         
         "PROCEDURE" COMDIV(XR,XI,YR,YI,ZR,ZI);"CODE" 34342;        
         TOL:= EM[1] * EM[2]; MACHTOL:= EM[0] * EM[1];    
         MAX:= EM[4]; COUNT:= 0; R:= 0; M:= N;  
         "IF" N > 1 "THEN" HC:= B[N - 1];       
         "FOR" I:= 1 "STEP" 1 "UNTIL" N "DO"    
         "BEGIN" VEC1[I,I]:= 1; VEC2[I,I]:= 0;  
             "FOR" J:= I + 1 "STEP" 1 "UNTIL" N "DO" VEC1[I,J]:=    
             VEC1[J,I]:= VEC2[I,J]:= VEC2[J,I]:= 0        
         "END";   
     IN: NM1:= N - 1;       
         "FOR" I:= N, I - 1 "WHILE" ("IF" I >= 1 "THEN" ABS(B[I]) > TOL       
         "ELSE" "FALSE") "DO" Q:= I; "IF" Q > 1 "THEN"    
         "BEGIN" "IF" ABS(B[Q - 1]) > R "THEN" R:= ABS(B[Q - 1]) "END";       
         "IF" Q = N "THEN"  
         "BEGIN" VAL1[N]:= A1[N,N]; VAL2[N]:= A2[N,N]; N:= NM1;     
             "IF" N > 1 "THEN" HC:= B[N - 1];   
         "END"    
         "ELSE"   
         "BEGIN" DD1:= A1[N,N]; DD2:= A2[N,N]; CC:= B[NM1];         
             P1:= (A1[NM1,NM1] - DD1) * .5;     
             P2:= (A2[NM1,NM1] - DD2) * .5;     
             COMKWD(P1, P2, CC * A1[NM1,N], CC * A2[NM1,N], G1,     
             G2, K1, K2); "IF" Q = NM1 "THEN"   
             "BEGIN" A1[N,N]:= VAL1[N]:= G1 + DD1;        
                 A2[N,N]:= VAL2[N]:= G2 + DD2;  
                 A1[Q,Q]:= VAL1[Q]:= K1 + DD1;  
                 A2[Q,Q]:= VAL2[Q]:= K2 + DD2;  
                 KAPPA:= SQRT(K1 ** 2 + K2 ** 2 + CC ** 2);         
                 NUI:= CC / KAPPA; MUI1:= K1 / KAPPA;     
                 MUI2:= K2 / KAPPA; AIJ1:= A1[Q,N];       
                 AIJ2:= A2[Q,N]; H1:= MUI1 ** 2 - MUI2 ** 2;        
                 H2:= 2 * MUI1 * MUI2; H:= - NUI * 2;     
                 A1[Q,N]:= H * (P1 * MUI1 + P2 * MUI2) - NUI *      
                 NUI * CC + AIJ1 * H1 + AIJ2 * H2;        
                 A2[Q,N]:= H * (P2 * MUI1 - P1 * MUI2) + AIJ2 *     
                 H1 - AIJ1 * H2;      
                 ROTCOMROW(Q + 2, M, Q, N, A1, A2, MUI1, MUI2,      
                 NUI);      
                 ROTCOMCOL(1, Q - 1, Q, N, A1, A2, MUI1, -
                 MUI2, - NUI);        
                 ROTCOMCOL(1, M, Q, N, VEC1, VEC2, MUI1, -
                 MUI2, - NUI); N:= N - 2;       
                 "IF" N > 1 "THEN" HC:= B[N - 1]; B[Q]:= 0
             "END"
1SECTION 3.3.2.2.1            (JULY 1974)                        PAGE 11      
        
        
        
             "ELSE"         
             "BEGIN" COUNT:= COUNT + 1;         
                 "IF" COUNT > MAX "THEN" "GOTO" OUT; Z1:= K1 + DD1; 
                 Z2:= K2 + DD2;       
                 "IF" ABS(CC) > ABS(HC) "THEN" Z1:= Z1 + ABS(CC);   
                 HC:= CC / 2; Q1:= Q + 1; AIJ1:= A1[Q,Q] - Z1;      
                 AIJ2:= A2[Q,Q] - Z2; AI1I:= B[Q];        
                 KAPPA:= SQRT(AIJ1 ** 2 + AIJ2 ** 2 + AI1I ** 2);   
                 MUI1:= AIJ1 / KAPPA; MUI2:= AIJ2 / KAPPA;
                 NUI:= AI1I / KAPPA; A1[Q,Q]:= KAPPA;     
                 A2[Q,Q]:= 0; A1[Q1,Q1]:= A1[Q1,Q1] - Z1; 
                 A2[Q1,Q1]:= A2[Q1,Q1] - Z2;    
                 ROTCOMROW(Q1, M, Q, Q1, A1, A2, MUI1, MUI2,        
                 NUI);      
                 ROTCOMCOL(1, Q, Q, Q1, A1, A2, MUI1, - MUI2, -     
                 NUI); A1[Q,Q]:= A1[Q,Q] + Z1;  
                 A2[Q,Q]:= A2[Q,Q] + Z2;        
                 ROTCOMCOL(1, M, Q, Q1, VEC1, VEC2, MUI1, -         
                 MUI2, - NUI);        
                 "FOR" I:= Q1 "STEP" 1 "UNTIL" NM1 "DO"   
                 "BEGIN" I1:= I + 1; AIJ1:= A1[I,I]; AIJ2:= A2[I,I];
                     AI1I:= B[I];     
                     KAPPA:= SQRT(AIJ1 ** 2 + AIJ2 ** 2 + AI1I **   
                     2); MUIM11:= MUI1; MUIM12:= MUI2;    
                     NUIM1:= NUI; MUI1:= AIJ1 / KAPPA;    
                     MUI2:= AIJ2 / KAPPA; NUI:= AI1I / KAPPA;       
                     A1[I1,I1]:= A1[I1,I1] - Z1;
                     A2[I1,I1]:= A2[I1,I1] - Z2;
                     ROTCOMROW(I1, M, I, I1, A1, A2, MUI1,
                     MUI2, NUI); A1[I,I]:= MUIM11 * KAPPA;
                     A2[I,I]:= - MUIM12 * KAPPA;
                     B[I - 1]:= NUIM1 * KAPPA;  
                     ROTCOMCOL(1, I, I, I1, A1, A2, MUI1, -         
                     MUI2, - NUI); A1[I,I]:= A1[I,I] + Z1;
                     A2[I,I]:= A2[I,I] + Z2;    
                     ROTCOMCOL(1, M, I, I1, VEC1, VEC2, MUI1, -     
                     MUI2, - NUI);    
                 "END" "COMMENT"      
1SECTION 3.3.2.2.1            (JULY 1974)                        PAGE 12      
                                                                 ;  
        
        
                 AIJ1:= A1[N,N]; AIJ2:= A2[N,N]; AIJ12:= AIJ1 ** 2; 
                 AIJ22:= AIJ2 ** 2; KAPPA:= SQRT(AIJ12 + AIJ22);    
                 "IF" ("IF" KAPPA < TOL "THEN" "TRUE" "ELSE" AIJ22 <=         
                 EM[0] * AIJ12) "THEN"
                 "BEGIN" B[NM1]:= NUI * AIJ1;   
                     A1[N,N]:= AIJ1 * MUI1 + Z1;
                     A2[N,N]:= - AIJ1 * MUI2 + Z2         
                 "END"      
                 "ELSE"     
                 "BEGIN" B[NM1]:= NUI * KAPPA; A1NN:= MUI1 * KAPPA; 
                     A2NN:= - MUI2 * KAPPA; MUI1:= AIJ1 / KAPPA;    
                     MUI2:= AIJ2 / KAPPA;       
                     COMCOLCST(1, NM1, N, A1, A2, MUI1, MUI2);      
                     COMCOLCST(1, NM1, N, VEC1, VEC2, MUI1,         
                     MUI2); 
                     COMROWCST(N + 1, M, N, A1, A2, MUI1, -         
                     MUI2); 
                     COMCOLCST(N, M, N, VEC1, VEC2, MUI1, MUI2);    
                     A1[N,N]:= MUI1 * A1NN - MUI2 * A2NN + Z1;      
                     A2[N,N]:= MUI1 * A2NN + MUI2 * A1NN + Z2;      
                 "END";     
             "END";         
         "END";   
         "IF" N > 0 "THEN" "GOTO" IN; 
         "FOR" J:= M "STEP" - 1 "UNTIL" 2 "DO"  
         "BEGIN" TF1[J]:= 1; TF2[J]:= 0; T1:= A1[J,J]; T2:= A2[J,J];
             "FOR" I:= J - 1 "STEP" - 1 "UNTIL" 1 "DO"    
             "BEGIN" DELTA1:= T1 - A1[I,I]; DELTA2:= T2 - A2[I,I];  
                 COMMATVEC(I + 1, J, I, A1, A2, TF1, TF2, MV1,      
                 MV2);      
                 "IF" ABS(DELTA1) < MACHTOL "AND" ABS(DELTA2) <     
                 MACHTOL "THEN"       
                 "BEGIN" TF1[I]:= MV1 / MACHTOL;
                     TF2[I]:= MV2 / MACHTOL     
                 "END"      
                 "ELSE" COMDIV(MV1, MV2, DELTA1, DELTA2, TF1[I],    
                 TF2[I]);   
             "END";         
             "FOR" I:= 1 "STEP" 1 "UNTIL" M "DO" COMMATVEC(1, J, I, 
             VEC1, VEC2, TF1, TF2, VEC1[I,J], VEC2[I,J]); 
         "END";   
     OUT: EM[3]:= R; EM[5]:= COUNT; QRICOM:= N; 
     "END" QRICOM;
         "EOP"    
1SECTION 3.3.2.2.2            (JULY 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 EIGVALCOM AND EIGCOM.     
     EIGVALCOM CALCULATES THE EIGENVALUES OF A COMPLEX MATRIX AND   
     EIGCOM CALCULATES THE EIGENVECTORS AS WELL.
        
        
 KEYWORDS :       
        
     EIGENVALUES, 
     EIGENVECTORS,
     COMPLEX MATRICES,      
     EQUILIBRATION,         
     REDUCTION HESSENBERG FORM,       
     HOUSEHOLDER TRANSFORMATION,      
     QR-ITERATION.
        
        
1SECTION 3.3.2.2.2            (JULY 1974)                         PAGE 2      
        
        
        
 SUBSECTION: EIGVALCOM.     
        
        
 CALLING SEQUENCE:
        
     THE HEADING OF THE PROCEDURE READS:        
     "INTEGER" "PROCEDURE" EIGVALCOM(AR, AI, N, EM, VALR, VALI);    
     "VALUE" N; "INTEGER" N; "ARRAY" AR, AI, EM, VALR, VALI;        
        
     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 MUST       
                BE GIVEN IN THE ARRAYS AR AND AI, RESPECTIVELY;     
                THE ELEMENTS OF THE ARRAYS AR AND AI ARE ALTERED;   
     N:         <ARITHMETIC EXPRESSION>;        
                THE ORDER OF THE GIVEN MATRIX;  
     EM:        <ARRAY IDENTIFIER>;   
                "ARRAY" EM[0:7];      
                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 QR-ITERATIONS; 
                       (E.G. 10 * N); 
                EM[6]: THE MAXIMUM ALLOWED NUMBER OF ITERATIONS FOR 
                       EQUILIBRATING THE ORIGINAL MATRIX (E.G. N**2/2);       
                EXIT:       
                EM[1]: THE EUCLIDEAN NORM OF THE EQUILIBRATED MATRIX;         
                EM[3]: THE MAXIMUM ABSOLUTE VALUE OF THE SUBDIAGONAL
                       ELEMENTS NEGLECTED IN THE QR-ITERATION;      
                EM[5]: THE NUMBER OF QR-ITERATIONS PERFORMED;       
                       EM[5]:=EM[4]+1 IN THE CASE EIGVALCOM^=0;     
                EM[7]: THE NUMBER OF ITERATIONS PERFORMED FOR       
                       EQUILIBRATING THE ORIGINAL MATRIX; 
     VALR,VALI: <ARRAY IDENTIFIER>;   
                "ARRAY" VALR,VALI[1:N];         
                EXIT:       
                THE REAL PART AND THE IMAGINARY PART OF THE CALCULATED        
                EIGENVALUES ARE DELIVERED IN THE ARRAYS VALR AND VALI,        
                RESPECTIVELY;         
        
     EIGVALCOM:=0, PROVIDED THE QR-ITERATION IS COMPLETED WITHIN EM[4]        
     ITERATIONS; OTHERWISE, EIGVALCOM:= THE NUMBER, K, OF EIGENVALUES         
     NOT CALCULATED AND ONLY THE LAST N-K ELEMENTS OF THE ARRAYS VALR         
     AND VALI ARE APPROXIMATE EIGENVALUES OF THE ORIGINAL MATRIX.   
        
 PROCEDURES USED: 
        
     EQILBRCOM = CP34361,   
     COMEUCNRM = CP34359,   
     HSHCOMHES = CP34366,   
     VALQRICOM = CP34372.   
        
        
1SECTION 3.3.2.2.2            (JULY 1974)                         PAGE 3      
        
        
        
 REQUIRED CENTRAL MEMORY: FIVE  REAL ARRAYS OF ORDER N AND ONE INTEGER        
     ARRAY OF ORDER N ARE DECLARED.   
        
        
 RUNNING TIME: PROPORTIONAL TO N ** 2 * MAX(N,NUMBER OF ITERATIONS).
        
        
 LANGUAGE: ALGOL 60.        
        
        
 METHOD AND PERFORMANCE: SEE EIGCOM (THIS SECTION).       
        
        
 EXAMPLE OF USE: SEE EIGCOM (THIS SECTION).     
        
        
 SUBSECTION: EIGCOM.        
        
        
 CALLING SEQUENCE:
        
     THE HEADING OF THE PROCEDURE READS:        
     "INTEGER" "PROCEDURE" EIGCOM(AR, AI, N, EM, VALR, VALI, VR, VI);         
     "VALUE" N; "INTEGER" N; "ARRAY" AR, AI, EM, VALR, VALI, VR, VI;
        
     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 MUST       
                BE GIVEN IN THE ARRAYS AR AND AI, RESPECTIVELY;     
                THE ELEMENTS OF THE ARRAYS AR AND AI ARE ALTERED;   
     N:         <ARITHMETIC EXPRESSION>;        
                THE ORDER OF THE GIVEN MATRIX;  
     EM:        <ARRAY IDENTIFIER>;   
                "ARRAY" EM[0:7];      
                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 QR-ITERATIONS; 
                       (E.G. 10 * N); 
                EM[6]: THE MAXIMUM ALLOWED NUMBER OF ITERATIONS FOR 
                       EQUILIBRATING THE ORIGINAL MATRIX (E.G. N**2/2);       
                EXIT:       
                EM[1]: THE EUCLIDEAN NORM OF THE EQUILIBRATED MATRIX;         
                EM[3]: THE MAXIMUM ABSOLUTE VALUE OF THE SUBDIAGONAL
                       ELEMENTS NEGLECTED IN THE QR-ITERATION;      
                EM[5]: THE NUMBER OF QR-ITERATIONS PERFORMED;       
                       EM[5]:=EM[4]+1 IN THE CASE EIGCOM^=0;        
                EM[7]: THE NUMBER OF ITERATIONS PERFORMED FOR       
                       EQUILIBRATING THE ORIGINAL MATRIX; 
1SECTION 3.3.2.2.2            (JULY 1974)                         PAGE 4      
        
        
        
     VALR,VALI: <ARRAY IDENTIFIER>;   
                "ARRAY" VALR,VALI[1:N];         
                EXIT:       
                THE REAL PART AND THE IMAGINARY PART OF THE CALCULATED        
                EIGENVALUES ARE DELIVERED IN THE ARRAYS VALR AND VALI,        
                RESPECTIVELY;         
     VR,VI:     <ARRAY IDENTIFIER>;   
                "ARRAY" VR,VI[1:N,1:N];         
                EXIT:       
                THE EIGENVECTORS OF THE MATRIX; 
                THE NORMALIZED EIGENVECTOR WITH REAL PART VR[1:N,J] AND       
                IMAGINARY PART VI[1:N,J] CORRESPONDS TO THE EIGENVALUE        
                VALR[J] + VALI[J] * I, J=1,...,N;         
        
     EIGCOM:=0, PROVIDED THE QR-ITERATION IS COMPLETED WITHIN EM[4] 
     ITERATIONS; OTHERWISE, EIGCOM:= THE NUMBER, K, OF EIGENVALUES  
     NOT CALCULATED AND ONLY THE LAST N-K ELEMENTS OF THE ARRAYS VALR         
     AND VALI ARE APPROXIMATE EIGENVALUES OF THE ORIGINAL MATRIX AND NO       
     USEFUL EIGENVECTORS ARE DELIVERED.         
        
        
 PROCEDURES USED: 
        
     EQILBRCOM = CP34361,   
     COMEUCNRM = CP34359,   
     HSHCOMHES = CP34366,   
     QRICOM    = CP34373,   
     BAKCOMHES = CP34367,   
     BAKLBRCOM = CP34362,   
     SCLCOM    = CP34360.   
        
        
 REQUIRED CENTRAL MEMORY: FIVE  REAL ARRAYS OF ORDER N AND ONE INTEGER        
     ARRAY OF ORDER N ARE DECLARED.   
        
        
 RUNNING TIME: PROPORTIONAL TO N**2 * MAX(N, NUMBER OF ITERATIONS). 
        
        
 LANGUAGE : ALGOL 60.       
        
        
1SECTION 3.3.2.2.2            (JULY 1974)                         PAGE 5      
        
        
        
 THE FOLLOWING HOLDS FOR BOTH PROCEDURES:       
        
        
 METHOD AND PERFORMANCE:    
        
     FOR CALCULATING THE EIGENVALUES AND EIGENVECTORS OF A COMPLEX  
     MATRIX WE DISTINGUISH THE FOLLOWING STEPS: 
     1) THE MATRIX IS EQUILIBRATED (SEE ALSO SECTION 3.2.1.1.2.).   
     2) THE EQUILIBRATED MATRIX IS TRANSFORMED INTO HESSENBERG FORM BY        
        MEANS OF HOUSEHOLDER MATRICES (SEE ALSO SECTION 3.2.1.2.2.2.).        
     3) THE HESSENBERG MATRIX IS TRANSFORMED INTO AN UPPER TRIANGULAR         
        MATRIX BY MEANS OF QR-ITERATION WITH SHIFT OF ORIGIN AND    
        DEFLATION (SEE ALSO SECTION 3.3.2.2.1.).
     THE DIAGONAL ELEMENTS OF THE UPPER TRIANGULAR MATRIX ARE THE   
     EIGENVALUES OF THE ORIGINAL MATRIX.        
     THE EIGENVECTORS OF THE ORIGINAL MATRIX ARE OBTAINED BY CALCULATING      
     THE EIGENVECTORS OF THE UPPER TRIANGULAR MATRIX (3) FOLLOWED BY
     BACKTRANSFORMATIONS CORRESPONDING TO (2) AND (1).    
        
        
 REFERENCES:      
        
     DEKKER, T.J. (1968),   
     ALGOL 60 PROCEDURES IN NUMERICAL ALGEBRA, PART 1,    
     MATH. CENTRE TRACTS 22, MATHEMATISCH CENTRUM;        
        
     DEKKER, T.J. AND W.HOFFMANN (1968),        
     ALGOL 60 PROCEDURES IN NUMERICAL ALGEBRA, PART 2,    
     MATH. CENTRE TRACTS 23, MATHEMATISCH CENTRUM;        
        
     FRANCIS, J.G.F. (1961),
     THE QR-TRANSFORMATION, PART 1 AND 2,       
     COMP.J. 4, P.265-271 AND P.332-345;        
        
     MUELLER, D.J. (1966),  
     HOUSEHOLDER,S METHOD  FOR  COMPLEX MATRICES AND EIGENSYSTEMS   OF        
     HERMITIAN MATRICES,    
     NUMER.MATH., 8, P.72-92;         
        
     OSBORNE, E.E. (1960),  
     ON PRECONDITIONING OF MATRICES,  
     JACM., 7, P.338-354;   
        
     PARLETT, B.N. AND C.REINSCH (1969),        
     BALANCING A MATRIX FOR CALCULATION OF EIGENVALUES AND
     EIGENVECTORS,
     NUM. MATH., 13, P.293-304;       
        
     WILKINSON, J.H. (1965),
     THE ALGEBRAIC EIGENVALUE PROBLEM,
     CLARENDOM PRESS, OXFORD.         
        
        
1SECTION 3.3.2.2.2            (JULY 1974)                         PAGE 6      
        
        
        
 EXAMPLE OF USE:  
        
     EIGCOM CALCULATES THE EIGENVALUES AND THE EIGENVECTORS OF THE  
     FOLLOWING MATRIX:      
     (SEE WILKINSON AND REINSCH, 1971, CONTRIBUTION II/15)
     1+3*I  2+1*I  3+2*I  1+1*I       
     3+4*I  1+2*I  2+1*I  4+3*I       
     2+3*I  1+5*I  3+1*I  5+2*I       
     1+2*I  3+1*I  1+4*I  5+3*I       
        
     ONLY THE EIGENVECTOR CORRESPONDING TO VALR[1] + VALI[1] * I IS 
     PRINTED BY THE FOLLOWING PROGRAM.
        
     "BEGIN"      
     "INTEGER" "PROCEDURE" EIGCOM(AR,AI,N,EM,VALR,VALI,VR,VI);      
     "CODE" 34375;
     "REAL" "ARRAY" AR,AI,VR,VI[1:4,1:4],EM[0:7],VALR,VALI[1:4];    
     "INTEGER" I; 
     AR[1,1]:=AR[1,4]:=AR[2,2]:=AR[3,2]:=AR[4,1]:=AR[4,3]:=         
     AI[1,2]:=AI[1,4]:=AI[2,3]:=AI[3,3]:=AI[4,2]:=1;      
     AR[1,2]:=AR[2,3]:=AR[3,1]:=AI[1,3]:=AI[2,2]:=AI[3,4]:=AI[4,1]:=2;        
     AR[1,3]:=AR[2,1]:=AR[3,3]:=AR[4,2]:=       
     AI[1,1]:=AI[2,4]:=AI[3,1]:=AI[4,4]:=3;     
     AR[2,4]:=AI[2,1]:=AI[4,3]:=4;    
     AR[3,4]:=AR[4,4]:=AI[3,2]:=5;    
     EM[0]:=5"-14;EM[2]:="-12;EM[4]:=10;EM[6]:=10;        
     OUTPUT(61,"(""("EIGCOM: ")",D")",
               EIGCOM(AR,AI,4,EM,VALR,VALI,VR,VI));       
     OUTPUT(61,"("/,"("EIGENVALUES:")",/")");   
     "FOR" I:=1,2,3,4 "DO" OUTPUT(61,"("2(+D.4D),"(" * I")",/")",   
                             VALR[I],VALI[I]);  
     OUTPUT(61,"(""("FIRST EIGENVECTOR:")",/")");         
     "FOR" I:=1,2,3,4 "DO" OUTPUT(61,"("2(+D.4D),"(" * I")",/")",   
                             VR[I,1],VI[I,1]);  
     OUTPUT(61,"(""("EM[1]: ")",+DD.DD/,"("EM[3]: ")",+D.D"+DD/,    
       "("EM[5]: ")",+ZD/,"("EM[7]: ")",+ZD")",EM[1],EM[3],EM[5],EM[7]);      
     "END"        
        
     OUTPUT:      
        
     EIGCOM: 0    
     EIGENVALUES: 
     -3.3710-0.7705 * I     
     +9.7837+9.3225 * I     
     +1.3657-1.4011 * I     
     +2.2217+1.8490 * I     
     FIRST EIGENVECTOR:     
     -0.5061+0.5835 * I     
     +1.0000+0.0000 * I     
     +0.5183-0.7147 * I     
     -0.5535+0.0188 * I     
     EM[1]: +15.30
     EM[3]: +6.0"-12        
     EM[5]:  +7   
     EM[7]:  +4   
        
        
1SECTION 3.3.2.2.2            (JULY 1974)                         PAGE 7      
        
        
        
 SOURCE TEXT(S) : 
0"CODE" 34374;    
     "INTEGER" "PROCEDURE" EIGVALCOM(AR, AI, N, EM, VALR, VALI);    
     "VALUE" N; "INTEGER" N; "ARRAY" AR, AI, EM, VALR, VALI;        
     "BEGIN" "INTEGER" "ARRAY" INT[1:N];        
         "ARRAY" D, B, DEL, TR, TI[1:N];        
         "PROCEDURE" HSHCOMHES(AR,AI,N,EM,B,TR,TI,DEL);"CODE" 34366;
         "REAL" "PROCEDURE" COMEUCNRM(AR,AI,LW,N);"CODE" 34359;     
         "PROCEDURE" EQILBRCOM(A1,A2,N,EM,D,INT);"CODE" 34361;      
         "INTEGER" "PROCEDURE" VALQRICOM(A1,A2,B,N,EM,VAL1,VAL2);   
         "CODE" 34372;      
         EQILBRCOM(AR, AI, N, EM, D, INT);      
         EM[1]:= COMEUCNRM(AR, AI, N - 1, N);   
         HSHCOMHES(AR, AI, N, EM, B, TR, TI, DEL);        
         EIGVALCOM:= VALQRICOM(AR, AI, B, N, EM, VALR, VALI)        
     "END" EIGVALCOM;       
         "EOP"    
0"CODE" 34375;    
     "INTEGER" "PROCEDURE" EIGCOM(AR, AI, N, EM, VALR, VALI, VR, VI);         
     "VALUE" N; "INTEGER" N; "ARRAY" AR, AI, EM, VALR, VALI, VR, VI;
     "BEGIN" "INTEGER" I;   
         "INTEGER" "ARRAY" INT[1:N];  
         "ARRAY" D, B, DEL, TR, TI[1:N];        
         "PROCEDURE" EQILBRCOM(A1,A2,N,EM,D,INT);"CODE" 34361;      
         "REAL" "PROCEDURE" COMEUCNRM(AR,AI,LW,N);"CODE" 34359;     
         "PROCEDURE" HSHCOMHES(AR,AI,N,EM,B,TR,TI,DEL);"CODE" 34366;
         "INTEGER" "PROCEDURE" QRICOM(A1,A2,B,N,EM,VAL1,VAL2,VEC1,VEC2);      
         "CODE" 34373;      
         "PROCEDURE" BAKCOMHES(AR,AI,TR,TI,DEL,VR,VI,N,N1,N2);      
         "CODE" 34367;      
         "PROCEDURE" BAKLBRCOM(N,N1,N2,D,INT,VR,VI);"CODE" 34362;   
         "PROCEDURE" SCLCOM(AR,AI,N,N1,N2);"CODE" 34360;  
         EQILBRCOM(AR, AI, N, EM, D, INT);      
         EM[1]:= COMEUCNRM(AR, AI, N - 1, N);   
         HSHCOMHES(AR, AI, N, EM, B, TR, TI, DEL);        
         I:= EIGCOM:= QRICOM(AR, AI, B, N, EM, VALR, VALI, VR,      
         VI); "IF" I = 0 "THEN"       
         "BEGIN" BAKCOMHES(AR, AI, TR, TI, DEL, VR, VI, N, 1, N);   
             BAKLBRCOM(N, 1, N, D, INT, VR, VI);
             SCLCOM(VR, VI, N, 1, N)  
         "END"    
     "END" EIGCOM;
         "EOP"    
1SECTION : 3.5.1.1            (DECEMBER 1975)                     PAGE 1      
        
        
        
 AUTHORS : G.H.GOLUB AND C.REINSCH    
        
        
 CONTRIBUTOR : D.T.WINTER   
        
        
 INSTITUTE : MATHEMATICAL CENTRE      
        
        
 RECEIVED : 731217
        
        
 BRIEF DESCRIPTION :        
     THIS   SECTION   CONTAINS   TWO   PROCEDURES,    QRISNGVALBID   AND      
     QRISNGVALDECBID. BOTH PROCEDURES CALCULATE THE SINGULAR VALUES OF A      
     BIDIAGONAL MATRIX.  MOREOVER,  THE SECOND PROCEDURE  CALCULATES THE      
     THE  SINGULAR VALUES DECOMPOSITION  OF A FULL MATRIX  OF WHICH  THE      
     BIDIAGONAL AND THE PRE- AND POSTMULTIPLYING MATRICES, AS CALCULATED      
     BY HSHREABID ARE GIVEN.
        
        
 KEYWORDS :       
     SINGULAR VALUES        
     QR ITERATION 
     BIDIAGONAL MATRICES    
        
        
1SECTION : 3.5.1.1            (DECEMBER 1975)                     PAGE 2      
        
        
        
 SUBSECTION : QRISNGVALBID  
        
        
 CALLING SEQUENCE :         
     THE HEADING OF THE PROCEDURE IS :
     "INTEGER" "PROCEDURE" QRISNGVALBID(D, B, N, EM);     
     "VALUE" N; "INTEGER" N; "ARRAY" D, B, EM;  
        
        
     THE MEANING OF THE FORMAL PARAMETERS IS:   
     D:  <ARRAY IDENTIFIER>;
         "ARRAY" D[1:N];    
         ENTRY: THE DIAGONAL OF THE BIDIAGONAL MATRIX;    
         EXIT: THE SINGULAR VALUES;   
     B:  <ARRAY IDENTIFIER>;
         "ARRAY" B[1:N];    
         ENTRY: THE SUPER DIAGONAL OF THE BIDIAGONAL MATRIX,IN B[1:N-1];      
     N:  <ARITHMETIC EXPRESSION>;     
         THE LENGTH OF B AND D;       
     EM: <ARRAY IDENTIFIER>;
         "ARRAY" EM[1:7];   
         ENTRY: EM[1]: THE INFINITY NORM OF THE MATRIX;   
             EM[2]: THE RELATIVE PRECISION IN THE SINGULAR VALUES;  
             EM[4]: THE MAXIMAL NUMBER OF ITERATIONS TO BE PERFORMED;         
             EM[6]: THE MINIMAL NON-NEGLECTABLE SINGULAR VALUE;     
         EXIT: EM[3]: THE MAXIMAL NEGLECTED SUPERDIAGONAL ELEMENT;  
             EM[5]: THE NUMBER OF ITERATIONS PERFORMED;   
             EM[7]: THE NUMERICAL RANK OF THE MATRIX, I.E. THE NUMBER OF      
                    SINGULAR VALUES GREATER THAN OR EQUAL TO EM[6]. 
        
     MOREOVER :   
     QRISNGVALBID:=  THE NUMBER OF  SINGULAR VALUES  NOT FOUND,  I.E.  A      
         NUMBER NOT  EQUAL TO ZERO  IF THE NUMBER OF  ITERATIONS EXCEEDS      
         EM[4].   
        
        
 PROCEDURES USED : NONE     
        
 REQUIRED CENTRAL MEMORY : NO AUXILIARY ARRAYS ARE DECLARED         
        
        
 RUNNING TIME :   
     THE RUNNING TIME DEPENDS STRONGLY UPON THE PROPERTIES OF THE MATRIX      
        
        
 METHOD AND PERFORMANCE :   
     THE METHOD  IS DESCRIBED  IN DETAIL   IN [1].  THIS PROCEDURE  IS A      
     REWRITING  OF  PART  OF  THE  PROCEDURE   SVD  PUBLISHED  THERE  BY      
     G.H.GOLUB AND C.REINSCH.         
        
        
 LANGUAGE : ALGOL 60.       
        
        
1SECTION : 3.5.1.1            (DECEMBER 1975)                     PAGE 3      
        
        
        
 SUBSECTION : QRISNGVALDECBID         
        
        
 CALLING SEQUENCE :         
     THE HEADING OF THE PROCEDURE IS :
     "INTEGER" "PROCEDURE" QRISNGVALDECBID(D, B, M, N, U, V, EM);   
     "VALUE" M, N; "INTEGER" M, N; "ARRAY" D, B, U, V, EM;
        
        
     THE MEANING OF THE FORMAL PARAMETERS IS :  
     D:  <ARRAY IDENTIFIER>;
         "ARRAY" D[1:N];    
         ENTRY: THE DIAGONAL OF THE BIDIAGONAL MATRIX;    
         EXIT: THE SINGULAR VALUES;   
     B:  <ARRAY IDENTIFIER>;
         "ARRAY" B[1:N];    
         ENTRY: THE SUPER DIAGONAL OF THE BIDIAGONAL MATRIX,IN B[1:N-1];      
     M:  <ARITHMETIC EXPRESSION>;     
         THE NUMBER OF ROWS OF THE MATRIX U;    
     N:  <ARITHMETIC EXPRESSION>;     
         THE LENGTH OF  B  AND  D,  THE NUMBER OF COLUMNS OF  U  AND THE      
         NUMBER OF COLUMNS AND ROWS OF V;       
     U:  <ARRAY IDENTIFIER>;
         "ARRAY" U[1:M,1:N];
         ENTRY:  THE  PREMULTIPLYING  MATRIX  AS  PRODUCED  BY PRETFMMAT      
             (SECTION 3.2.2.1.1);     
         EXIT:  THE  PREMULTIPLYING  MATRIX  U  OF  THE  SINGULAR VALUES      
             DECOMPOSITION U * S * V';
     V:  <ARRAY IDENTIFIER>;
         "ARRAY" V[1:N,1:N];
         ENTRY:  THE TRANSPOSE OF THE POSTMULTIPLYING MATRIX AS PRODUCED      
             BY PSTTFMMAT (SECTION 3.2.2.1.1);  
         EXIT:  THE  TRANSPOSE  OF  THE POSTMULTIPLYING MATRIX  V OF THE      
             SINGULAR VALUES DECOMPOSITION;     
     EM: <ARRAY IDENTIFIER>;
         "ARRAY" EM[1:7];   
         ENTRY: EM[1]: THE INFINITY NORM OF THE MATRIX;   
             EM[2]: THE RELATIVE PRECISION IN THE SINGULAR VALUES;  
             EM[4]: THE MAXIMAL NUMBER OF ITERATIONS TO BE PERFORMED;         
             EM[6]: THE MINIMAL NON-NEGLECTABLE SINGULAR VALUE;     
         EXIT: EM[3]: THE MAXIMAL NEGLECTED SUPERDIAGONAL ELEMENT;  
             EM[5]: THE NUMBER OF ITERATIONS PERFORMED;   
             EM[7]: THE NUMERICAL RANK OF THE MATRIX, I.E. THE NUMBER OF      
                 SINGULAR VALUES GREATER THAN OR EQUAL TO EM[6].    
        
     MOREOVER :   
     QRISNGVALDECBID:= THE NUMBER OF SINGULAR VALUES NOT FOUND,  I.E.  A      
         NUMBER NOT  EQUAL TO ZERO  IF THE NUMBER OF  ITERATIONS EXCEEDS      
         EM[4].   
        
1SECTION : 3.5.1.1            (DECEMBER 1975)                     PAGE 4      
        
        
        
 PROCEDURES USED :
        
     ROTCOL = CP34040       
        
 REQUIRED CENTRAL MEMORY : NO AUXILIARY ARRAYS ARE DECLARED         
        
        
 RUNNING TIME :   
     THE RUNNING TIME DEPENDS STRONGLY UPON THE PROPERTIES OF THE MATRIX      
        
 METHOD AND PERFORMANCE :   
     THE METHOD  IS DESCRIBED  IN DETAIL   IN [1].  THIS PROCEDURE  IS A      
     REWRITING  OF  PART  OF  THE  PROCEDURE   SVD  PUBLISHED  THERE  BY      
     G.H.GOLUB AND C.REINSCH.         
        
        
 LANGUAGE : ALGOL 60        
        
        
 REFERENCES :     
     [1] WILKINSON, J.H. AND C.REINSCH
         HANDBOOK OF 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      
        
        
 SOURCE TEXT(S):  
0"CODE" 34270;    
 "INTEGER" "PROCEDURE" QRISNGVALBID(D, B, N, EM);         
 "VALUE" N; "INTEGER" N; "ARRAY" D, B, EM;      
 "BEGIN" "INTEGER" N1, K, K1, I, I1, COUNT, MAX, RNK;     
     "REAL" TOL, BMAX, Z, X, Y, G, H, F, C, S, MIN;       
     TOL:= EM[2] * EM[1]; COUNT:= 0; BMAX:= 0; MAX:= EM[4]; MIN:= EM[6];      
     RNK:= N;     
 IN: K:= N; N1:= N - 1;     
 NEXT: K:= K - 1; "IF" K > 0 "THEN"   
     "BEGIN" "IF" ABS(B[K]) >= TOL "THEN"       
         "BEGIN" "IF" ABS(D[K]) >= TOL "THEN" "GOTO" NEXT;
             C:= 0; S:= 1;  
             "FOR" I:= K "STEP" 1 "UNTIL" N1 "DO"         
             "BEGIN" F:= S * B[I]; B[I]:= C * B[I]; I1:= I + 1;     
                 "IF" ABS(F) < TOL "THEN" "GOTO" NEGLECT; 
                 G:= D[I1]; D[I1]:= H:= SQRT(F * F + G * G);        
                 C:= G / H; S:= - F / H         
             "END"
1SECTION : 3.5.1.1            (JULY 1974)                         PAGE 5      
                                                                  ; 
        
        
         NEGLECT: 
         "END"    
         "ELSE" "IF" ABS(B[K]) > BMAX "THEN" BMAX:= ABS(B[K])       
     "END";       
     "IF" K = N1 "THEN"     
     "BEGIN" "IF" D[N] < 0 "THEN" D[N]:= - D[N];
         "IF" D[N] <= MIN "THEN" RNK:= RNK - 1; N:= N1    
     "END"        
     "ELSE"       
     "BEGIN" COUNT:= COUNT + 1; "IF" COUNT > MAX "THEN" "GOTO" END; 
         K1:= K + 1; Z:= D[N]; X:= D[K1]; Y:= D[N1];      
         G:= "IF" N1 = 1 "THEN" 0 "ELSE" B[N1 - 1]; H:= B[N1];      
         F:= ((Y - Z) * (Y + Z) + (G - H) * (G + H)) / (2 * H * Y); 
         G:= SQRT(F * F + 1);         
         F:= ((X - Z) * (X + Z) + H * (Y / ("IF" F < 0 "THEN" F - G 
         "ELSE" F + G) - H)) / X; C:= S:= 1;    
         "FOR" I:= K1 + 1 "STEP" 1 "UNTIL" N "DO"         
         "BEGIN" I1:= I - 1; G:= B[I1]; Y:= D[I]; H:= S * G; G:= C * G;       
             Z:= SQRT(F * F + H * H); C:= F / Z; S:= H / Z;         
             "IF" I1 ^= K1 "THEN" B[I1 - 1]:= Z; F:= X * C + G * S; 
             G:= G * C - X * S; H:= Y * S; Y:= Y * C;     
             D[I1]:= Z:= SQRT(F * F + H * H); C:= F / Z; S:= H / Z; 
             F:= C * G + S * Y; X:= C * Y - S * G         
         "END";   
         B[N1]:= F; D[N]:= X
     "END";       
     "IF" N > 0 "THEN" "GOTO" IN;     
 END: EM[3]:= BMAX; EM[5]:= COUNT; EM[7]:= RNK; QRISNGVALBID:= N    
 "END" QRISNGVALBID;        
         "EOP"    
0"CODE" 34271;    
 "INTEGER" "PROCEDURE" QRISNGVALDECBID(D, B, M, N, U, V, EM);       
 "VALUE" M, N; "INTEGER" M, N; "ARRAY" D, B, U, V, EM;    
 "BEGIN" "INTEGER" N0, N1, K, K1, I, I1, COUNT, MAX, RNK; 
     "REAL" TOL, BMAX, Z, X, Y, G, H, F, C, S, MIN;       
        
     "PROCEDURE" ROTCOL(L, U, I, J, A, C, S);   
     "VALUE" L, U, I, J, C, S; "INTEGER" L, U, I, J;      
     "REAL" C, S; "ARRAY" A;
     "CODE" 34040;"COMMENT" 
1SECTION : 3.5.1.1            (JULY 1974)                         PAGE 6      
                                                                  ; 
        
        
        
     TOL:= EM[2] * EM[1]; COUNT:= 0; BMAX:= 0; MAX:= EM[4]; MIN:= EM[6];      
     RNK:= N0:= N;
 IN: K:= N; N1:= N - 1;     
 NEXT: K:= K - 1; "IF" K > 0 "THEN"   
     "BEGIN" "IF" ABS(B[K]) >= TOL "THEN"       
         "BEGIN" "IF" ABS(D[K]) >= TOL "THEN" "GOTO" NEXT;
             C:= 0; S:= 1;  
             "FOR" I:= K "STEP" 1 "UNTIL" N1 "DO"         
             "BEGIN" F:= S * B[I]; B[I]:= C * B[I]; I1:= I + 1;     
                 "IF" ABS(F) < TOL "THEN" "GOTO" NEGLECT; 
                 G:= D[I1]; D[I1]:= H:= SQRT(F * F + G * G);        
                 C:= G / H; S:= - F / H;        
                 ROTCOL(1, M, K, I1, U, C, S)   
             "END";         
         NEGLECT: 
         "END"    
         "ELSE" "IF" ABS(B[K]) > BMAX "THEN" BMAX:= ABS(B[K])       
     "END";       
     "IF" K = N1 "THEN"     
     "BEGIN" "IF" D[N] < 0 "THEN"     
         "BEGIN" D[N]:= - D[N];       
             "FOR" I:= 1 "STEP" 1 "UNTIL" N0 "DO" V[I,N]:= - V[I,N] 
         "END";   
         "IF" D[N] <= MIN "THEN" RNK:= RNK - 1; N:= N1    
     "END"        
     "ELSE"       
     "BEGIN" COUNT:= COUNT + 1; "IF" COUNT > MAX "THEN" "GOTO" END; 
         K1:= K + 1; Z:= D[N]; X:= D[K1]; Y:= D[N1];      
         G:= "IF" N1 = 1 "THEN" 0 "ELSE" B[N1 - 1]; H:= B[N1];      
         F:= ((Y - Z) * (Y + Z) + (G - H) * (G + H)) / (2 * H * Y); 
         G:= SQRT(F * F + 1);         
         F:= ((X - Z) * (X + Z) + H * (Y / ("IF" F < 0 "THEN" F - G 
         "ELSE" F + G) - H)) / X; C:= S:= 1;    
         "FOR" I:= K1 + 1 "STEP" 1 "UNTIL" N "DO"         
         "BEGIN" I1:= I - 1; G:= B[I1]; Y:= D[I]; H:= S * G; G:= C * G;       
             Z:= SQRT(F * F + H * H); C:= F / Z; S:= H / Z;         
             "IF" I1 ^= K1 "THEN" B[I1 - 1]:= Z; F:= X * C + G * S; 
             G:= G * C - X * S; H:= Y * S; Y:= Y * C;     
             ROTCOL(1, N0, I1, I, V, C, S);     
             D[I1]:= Z:= SQRT(F * F + H * H); C:= F / Z; S:= H / Z; 
             F:= C * G + S * Y; X:= C * Y - S * G;        
             ROTCOL(1, M, I1, I, U, C, S)       
         "END";   
         B[N1]:= F; D[N]:= X
     "END";       
     "IF" N > 0 "THEN" "GOTO" IN;     
 END: EM[3]:= BMAX; EM[5]:= COUNT; EM[7]:= RNK; QRISNGVALDECBID:= N 
 "END" QRISNGVALDECBID;     
         "EOP"    
1SECTION : 3.5.1.2            (JULY 1974)                         PAGE 1      
        
        
        
 AUTHOR : D.T.WINTER        
        
        
 INSTITUTE : MATHEMATICAL CENTRE      
        
        
 RECEIVED : 731217
        
        
 BRIEF DESCRIPTION :        
     THIS SECTION CONTAINS  TWO PROCEDURES,  QRISNGVAL AND QRISNGVALDEC.      
     QRISNGVAL  CALCULATES  THE SINGULAR VALUES  OF A GIVEN MATRIX. 
     QRISNGVALDEC  CALCULATES  THE  SINGULAR  VALUES  DECOMPOSITION 
     U * S * V',  WITH U AND V ORTHOGONAL AND S POSITIVE DIAGONAL.  
        
        
 KEYWORDS :       
     SINGULAR VALUES        
     QR ITERATION 
        
        
1SECTION : 3.5.1.2            (DECEMBER 1975)                     PAGE 2      
        
        
        
 SUBSECTION : QRISNGVAL     
        
        
 CALLING SEQUENCE :         
     THE HEADING OF THE PROCEDURE IS :
     "INTEGER" "PROCEDURE" QRISNGVAL(A, M, N, VAL, EM);   
     "VALUE" M, N; "INTEGER" M, N; "ARRAY" A, VAL, EM;    
        
     THE MEANING OF THE FORMAL PARAMETERS IS :  
     A:  <ARRAY IDENTIFIER>;
         "ARRAY" A[1:M,1:N];
         ENTRY: THE INPUT MATRIX;     
         EXIT: DATA CONCERNING THE TRANSFORMATION TO BIDIAGONAL FORM;         
     M:  <ARITHMETIC EXPRESSION>;     
         THE NUMBER OF ROWS OF A;     
     N:  <ARITHMETIC EXPRESSION>;     
         THE NUMBER OF COLUMNS OF A, N SHOULD SATISFY N <= M;       
     VAL: <ARRAY IDENTIFIER>;         
         "ARRAY" VAL[1:N];  
         EXIT: THE SINGULAR VALUES;   
     EM: <ARRAY IDENTIFIER>;
         "ARRAY" EM[0:7];   
         ENTRY: EM[0]: THE MACHINE PRECISION;   
             EM[2]: THE RELATIVE PRECISION IN THE SINGULAR VALUES;  
             EM[4]: THE MAXIMAL NUMBER OF ITERATIONS TO BE PERFORMED;         
             EM[6]: THE MINIMAL NON-NEGLECTABLE SINGULAR VALUE;     
         EXIT: EM[1]: THE INFINITY NORM OF THE MATRIX;    
             EM[3]: THE MAXIMAL NEGLECTED SUPERDIAGONAL ELEMENT;    
             EM[5]: THE NUMBER OF ITERATIONS PERFORMED;   
             EM[7]: THE NUMERICAL RANK OF THE MATRIX, I.E. THE NUMBER OF      
                 SINGULAR VALUES GREATER THAN OR EQUAL TO EM[6].    
        
     MOREOVER :   
     QRISNGVAL:= THE NUMBER OF SINGULAR VALUES NOT FOUND,  I.E. A NUMBER      
         NOT EQUAL TO ZERO IF THE NUMBER OF ITERATIONS EXCEEDS EM[4].         
        
 PROCEDURES USED :
     HSHREABID    = CP34260 
     QRISNGVALBID = CP34270 
        
 REQUIRED CENTRAL MEMORY : AN AUXILIARY ARRAY OF N REALS IS DECLARED
        
 RUNNING TIME :   
     THE RUNNING TIME DEPENDS UPON THE PROPERTIES OF THE MATRIX, HOWEVER      
     THE PROCESS OF BIDIAGONALIZATION DOMINATES,  AND ITS RUNNING TIME        
     IS PROPORTIONAL TO  (M + N) * N * N        
        
        
 METHOD AND PERFORMANCE :   
     THE MATRIX IS FIRST TRANSFORMED TO BIDIAGONAL FORM BY THE PROCEDURE      
     HSHREABID  (SECTION 3.2.2.1.1) ,  AND THEN  THE SINGULAR VALUES ARE      
     CALCULATED BY QRISNGVALBID (SECTION 3.5.1.1).        
        
 LANGUAGE : ALGOL 60        
1SECTION : 3.5.1.2            (DECEMBER 1975)                     PAGE 3      
        
        
        
 SUBSECTION : QRISNGVALDEC  
        
        
 CALLING SEQUENCE :         
     THE HEADING OF THE PROCEDURE IS :
     "INTEGER" "PROCEDURE" QRISNGVALDEC(A, M, N, VAL, V, EM);       
     "VALUE" M, N; "INTEGER" M, N; "ARRAY" A, VAL, V, EM; 
        
        
     THE MEANING OF THE FORMAL PARAMETERS IS:   
     A:  <ARRAY IDENTIFIER>;
         "ARRAY" A[1:M,1:N];
         ENTRY: THE GIVEN MATRIX;     
         EXIT:  THE  MATRIX  U  IN  THE  SINGULAR  VALUES  DECOMPOSITION      
             U * S * V';    
     M:  <ARITHMETIC EXPRESSION>;     
         THE NUMBER OF ROWS OF A;     
     N:  <ARITHMETIC EXPRESSION>;     
         THE NUMBER OF COLUMNS OF A, N SHOULD SATISFY N <= M;       
     VAL: <ARRAY IDENTIFIER>;         
         "ARRAY" VAL[1:N];  
         EXIT: THE SINGULAR VALUES;   
     V:  <ARRAY IDENTIFIER>;
         "ARRAY" V[1:N,1:N];
         EXIT: THE TRANSPOSE OF MATRIX V IN THE SINGULAR VALUES     
             DECOMPOSITION; 
     EM: <ARRAY IDENTIFIER>;
         "ARRAY" EM[0:7];   
         ENTRY: EM[0]: THE MACHINE PRECISION;   
             EM[2]: THE RELATIVE PRECISION IN THE SINGULAR VALUES;  
             EM[4]: THE MAXIMAL NUMBER OF ITERATIONS TO BE PERFORMED;         
             EM[6]: THE MINIMAL NON-NEGLECTABLE SINGULAR VALUE;     
         EXIT: EM[1]: THE INFINITY NORM OF THE MATRIX;    
             EM[3]: THE MAXIMAL NEGLECTED SUPER DIAGONAL ELEMENT;   
             EM[5]: THE NUMBER OF ITERATIONS PERFORMED;   
             EM[7]: THE NUMERICAL RANK OF THE MATRIX, I.E. THE NUMBER OF      
                 SINGULAR VALUES GREATER THAN OR EQUAL TO EM[6].    
        
     MOREOVER :   
     QRISNGVALDEC:=  THE NUMBER OF  SINGULAR VALUES  NOT FOUND,  I.E.  A      
         NUMBER NOT  EQUAL TO ZERO  IF THE NUMBER OF  ITERATIONS EXCEEDS      
         EM[4].   
        
 PROCEDURES USED :
     HSHREABID       = CP34260        
     PSTTFMMAT       = CP34261        
     PRETFMMAT       = CP34262        
     QRISNGVALDECBID = CP34271        
        
        
 REQUIRED CENTRAL MEMORY : AN AUXILIARY ARRAY OF N ELEMENTS IS DECLARED       
        
        
1SECTION : 3.5.1.2            (DECEMBER 1975)                     PAGE 4      
        
        
        
 RUNNING TIME :   
     THE RUNNING TIME DEPENDS UPON THE PROPERTIES OF THE MATRIX, HOWEVER      
     THE PROCESS OF BIDIAGONALIZATION DOMINATES,  AND ITS RUNNING TIME        
     IS PROPORTIONAL TO (M + N) * N * N         
        
 METHOD AND PERFORMANCE:    
     THE MATRIX IS FIRST TRANSFORMED TO BIDIAGONAL FORM BY THE PROCEDURE      
     HSHREABID (SECTION 3.2.2.1.1), THE TWO TRANSFORMING MATRICES ARE         
     CALCULATED BY THE PROCEDURES PSTTFMMAT AND PRETFMMAT (SECTIONS 
     3.2.2.1.2 AND 3.2.2.1.3 RESPECTIVELY), AND FINALLY THE SINGULAR
     VALUES DECOMPOSITION IS CALCULATED BY QRISNGVALDECBID (SECTION 
     3.5.1.1).    
        
 LANGUAGE : ALGOL 60        
        
 REFERENCES :     
         WILKINSON, J.H. AND C.REINSCH
         HANDBOOK OF AUTOMATIC COMPUTATION, VOL. 2        
         LINEAR ALGEBRA     
         HEIDELBERG (1971)  
        
 EXAMPLE OF USE : 
     AS THE PROCEDURE QRISNGVALDEC  CALCULATES THE  SINGULAR VALUES OF A      
     MATRIX IN EXACTLY THE SAME WAY AS  QRISNGVAL, WE GIVE  HER ONLY  AN      
     EXAMPLE OF  USE OF  THE PROCEDURE  QRISNGVALDEC.  FIRST  WE GIVE  A      
     PROGRAM, AND THEN THE RESULTS OF THIS PROGRAM:       
        
     "BEGIN" "ARRAY" A[1:6,1:5], V[1:5,1:5], VAL[1:5], EM[0:7];     
         "INTEGER" I, J;    
         "INTEGER" "PROCEDURE" QRISNGVALDEC(A, M, N, VAL, V, EM);   
         "VALUE" M, N; "INTEGER" M, N; "ARRAY" A, VAL, V, EM;       
         "CODE" 34273;      
        
         "FOR" I:= 1 "STEP" 1 "UNTIL" 6 "DO"    
         "FOR" J:= 1 "STEP" 1 "UNTIL" 5 "DO"    
         A[I,J]:= 1 / (I + J - 1);    
         EM[0]:= "-14; EM[2]:= "-12; EM[4]:= 25; EM[6]:= "-10;      
         I:= QRISNGVALDEC(A, 6, 5, VAL, V, EM); 
         OUTPUT(61, "("3B, "("NUMBER SINGULAR VALUES NOT FOUND : ")",         
         3ZD, /, 3B, "("INFINITY NORM : ")", N, /, 3B,    
         "("MAX NEGLECTED SUBDIAGONAL ELEMENT : ")", N, /, 3B,      
         "("NUMBER ITERATIONS : ")", 3ZD, /, 3B,
         "("NUMERICAL RANK : ")", 3ZD, /")", I, EM[1], EM[3], EM[5],
         EM[7]);  
         OUTPUT(61, "("/, 3B, "("SINGULAR VALUES : ")", /")");      
         "FOR" I:= 1 "STEP" 1 "UNTIL" 5 "DO"    
         OUTPUT(61, "("/, 3B, N")", VAL[I]);    
         OUTPUT(61, "("/, /, 3B, "("MATRIX U, FIRST 3 COLUMNS")", /")");      
         "FOR" I:= 1 "STEP" 1 "UNTIL" 6 "DO"    
         OUTPUT(61, "("/, 3B, 3(N)")", A[I,1], A[I,2], A[I,3]);     
         OUTPUT(61, "("/, /, 13B, "("LAST 2 COLUMNS")", /")");      
         "FOR" I:= 1 "STEP" 1 "UNTIL" 6 "DO"    
         OUTPUT(61, "("/, 13B, 2(N)")", A[I,4], A[I,5])   
     "END"        
1SECTION : 3.5.1.2            (JULY 1974)                         PAGE 5      
        
        
        
        
     NUMBER SINGULAR VALUES NOT FOUND :    0    
     INFINITY NORM : +2.2833333333334"+000      
     MAX NEGLECTED SUBDIAGONAL ELEMENT : +5.7786437871158"-014      
     NUMBER ITERATIONS :    5         
     NUMERICAL RANK :    5  
        
     SINGULAR VALUES :      
        
     +1.5921172587262"+000  
     +2.2449595426097"-001  
     +1.3610556101029"-002  
     +4.3245382038374"-004  
     +6.4001947134260"-006  
        
     MATRIX U, FIRST 3 COLUMNS        
        
     -7.5497918208386"-001  +6.1011090790645"-001  -2.3287173869184"-001      
     -4.3909273679284"-001  -2.2602102994174"-001  +7.0245315582712"-001      
     -3.1703146681544"-001  -3.7306964696148"-001  +2.1607293656979"-001      
     -2.4999458583084"-001  -3.9557817833576"-001  -1.4665595223684"-001      
     -2.0704999076883"-001  -3.8483260608872"-001  -3.6803786187007"-001      
     -1.7699734614538"-001  -3.6458192866515"-001  -4.9868122801331"-001      
        
               LAST 2 COLUMNS         
        
               +5.8625326935176"-002  -1.0184205426735"-002         
               -4.8169088124009"-001  +1.7189132301455"-001         
               +5.4982292571999"-001  -5.9788920283495"-001         
               +4.0633053815463"-001  +4.5989617524697"-001         
               -6.1755991033503"-002  +4.3029765325422"-001         
               -5.4158416488948"-001  -4.6499203623570"-001         
        
        
1SECTION : 3.5.1.2            (JULY 1974)                         PAGE 6      
        
        
        
 SOURCE TEXT(S):  
0"CODE" 34272;    
 "INTEGER" "PROCEDURE" QRISNGVAL(A, M, N, VAL, EM);       
 "VALUE" M, N; "INTEGER" M, N; "ARRAY" A, VAL, EM;        
 "BEGIN" "ARRAY" B[1:N];    
        
     "PROCEDURE" HSHREABID(A, M, N, D, B, EM);  
     "VALUE" M, N; "INTEGER" M, N; "ARRAY" D, B, EM;      
     "CODE" 34260;
        
     "INTEGER" "PROCEDURE" QRISNGVALBID(D, B, N, EM);     
     "VALUE" N; "INTEGER" N; "ARRAY" D, B, EM;  
     "CODE" 34270;
        
     HSHREABID(A, M, N, VAL, B, EM);  
     QRISNGVAL:= QRISNGVALBID(VAL, B, N, EM)    
 "END" QRISNGVAL; 
         "EOP"    
0"CODE" 34273;    
 "INTEGER" "PROCEDURE" QRISNGVALDEC(A, M, N, VAL, V, EM); 
 "VALUE" M, N; "INTEGER" M, N; "ARRAY" A, VAL, V, EM;     
 "BEGIN" "ARRAY" B[1:N];    
        
     "PROCEDURE" HSHREABID(A, M, N, D, B, EM);  
     "VALUE" M, N; "INTEGER" M, N; "ARRAY" A, D, B, EM;   
     "CODE" 34260;
        
     "PROCEDURE" PSTTFMMAT(A, N, V, B);         
     "VALUE" N; "INTEGER" N; "ARRAY" A, V, B;   
     "CODE" 34261;
        
     "PROCEDURE" PRETFMMAT(A, M, N, D);         
     "VALUE" M, N; "INTEGER" M, N; "ARRAY" A, D;
     "CODE" 34262;
        
     "INTEGER" "PROCEDURE" QRISNGVALDECBID(D, B, M, N, U, V, EM);   
     "VALUE" M, N; "INTEGER" M, N; "ARRAY" D, B, U, V, EM;
     "CODE" 34271;
        
     HSHREABID(A, M, N, VAL, B, EM);  
     PSTTFMMAT(A, N, V, B);  PRETFMMAT(A, M, N, VAL);     
     QRISNGVALDEC:= QRISNGVALDECBID(VAL, B, M, N, A, V, EM)         
 "END" QRISNGVALDEC;        
         "EOP"    
1SECTION : 3.6.3              (JULY 1974)                         PAGE 1      
        
        
        
 AUTHOR: C.G. VAN DER LAAN. 
        
        
 INSTITUTE: MATHEMATICAL CENTRE.      
        
        
 RECEIVED: 730815.
        
        
 BRIEF DESCRIPTION:         
        
     COMKWD CALCULATES THE ROOTS OF A QUADRATIC EQUATION WITH COMPLEX         
     COEFFICIENTS.
        
        
 KEYWORDS:        
        
     ZEROS,QUADRATIC EQUATION,POLYNOMIAL EQUATION,COMPLEX COEFFICIENTS.       
        
        
 CALLING SEQUENCE:
        
     THE HEADING OF THE PROCEDURE READS:        
     "PROCEDURE"COMKWD(PR,PI,QR,QI,GR,GI,KR,KI);
     "VALUE"PR,PI,QR,QI;"REAL"PR,PI,QR,QI,GR,GI,KR,KI;    
        
     THE MEANING OF THE FORMAL PARAMETERS IS:   
     PR,PI,QR,QI:<ARITHMETIC EXPRESSION>;       
                 ENTRY:PR,QR ARE THE REAL PARTS AND PI,QI ARE THE   
                       IMAGINARY PARTS OF THE COEFFICIENTS OF THE   
                       QUADRATIC EQUATION:      
                       X**2-2*(PR+I*PI)*X-(QR+I*QI)=0;    
     GR,GI,KR,KI:<VARIABLE>;
                 EXIT:THE REAL PARTS AND THE IMAGINARY PARTS OF THE 
                      DINOMIAL ARE DELIVERED IN GR,KR AND GI,KI,    
                      RESPECTIVELY;MOREOVER,THE MODULUS OF GR+I*GI IS         
                      GREATER OR EQUAL THE MODULUS OF KR+I*KI.      
        
        
 PROCEDURES USED: 
        
     COMMUL=CP34341;        
     COMDIV=CP34342;        
     COMSQRT=CP34343.       
        
        
 LANGUAGE: ALGOL 60.        
        
1SECTION : 3.6.3              (JULY 1974)                         PAGE 2      
        
        
        
        
 EXAMPLE OF USE:  
        
 "BEGIN""REAL"GR,GI,KR,KI;  
 "PROCEDURE" COMKWD(PR,PI,QR,QI,GR,GI,KR,KI);   
 "CODE"34345;     
 COMKWD(-.1,.3,.11,.02,GR,GI,KR,KI);  
 OUTPUT(61,"(""("X**2-2(-.1+.3*I)*X-( .11+.02*I) HAS ROOTS")",/,    
     -D.DD,+D.DD,"("*I")",/,
     -D.DD,+D.DD,"("*I")"")",GR,GI,KR,KI)       
 "END"  
        
 X**2-2(-.1+.3*I)*X-( .11+.02*I) HAS ROOTS      
 -0.30+0.40*I     
  0.10+0.20*I     
        
        
 SOURCE TEXT(S):  
0"CODE"34345;     
 "PROCEDURE" COMKWD(PR,PI,QR,QI,GR,GI,KR,KI);   
 "VALUE" PR,PI,QR,QI;"REAL" PR,PI,QR,QI,GR,GI,KR,KI;      
 "BEGIN"
 "PROCEDURE"COMMUL (AR,AI,BR,BI,RR,RI);         
 "CODE"34341;     
 "PROCEDURE"COMDIV(XR,XI,YR,YI,ZR,ZI);
 "CODE"34342;     
 "PROCEDURE"COMSQRT(AR,AI,PR,PI);     
 "CODE"34343;     
 "IF" QR=0 & QI = 0 "THEN"  
 "BEGIN" KR:=KI:=0 ;GR:=PR*2;GI:=PI*2 "END" "ELSE"        
 "IF" PR=0 & PI= 0 "THEN"   
 "BEGIN" COMSQRT(QR,QI,GR,GI);KR:=-GR;KI:=-GI "END" "ELSE"
 "BEGIN" "REAL" HR,HI;      
 "IF" ABS(PR) > 1 "OR" ABS(PI) >1 "THEN" "BEGIN"
 COMDIV(QR,QI,PR,PI,HR,HI); 
 COMDIV(HR,HI,PR,PI,HR,HI); 
 COMSQRT(1+HR,HI,HR,HI);    
 COMMUL(PR,PI,HR+1,HI,GR,GI);         
 "END" "ELSE" "BEGIN" COMSQRT(QR+(PR+PI)*(PR-PI),QI+ PR*PI*2,HR,HI);
        "IF" PR * HR + PI * HI > 0 "THEN"       
        "BEGIN" GR := PR + HR;GI := PI + HI "END" "ELSE"  
        "BEGIN" GR := PR - HR;GI:= PI - HI "END";         
 "END"; 
 COMDIV(-QR,-QI,GR,GI,KR,KI);         
 "END"  
 "END" COMKWD;    
         "EOP"    
1SECTION : 4.1                (JULY 1974)                         PAGE 1      
        
        
        
 AUTHOR     :    J.W. DANIEL.         
        
        
 REVISOR    :    J. KOK.    
        
        
 INSTITUTE  :    MATHEMATICAL CENTRE. 
        
        
 RECEIVED   :    730528 (EULER).      
                 730917 (SUMPOSSERIES).         
        
        
 BRIEF DESCRIPTION  :       
        
     THIS  SECTION  CONTAINS  TWO  PROCEDURES  FOR  THE  SUMMATION  OF        
     CONVERGENT INFINITE SERIES:      
        
     EULER PERFORMS THE  SUMMATION OF AN ALTERNATING SERIES.        
        
     SUMPOSSERIES   PERFORMS  THE  SUMMATION OF A CONVERGENT SERIES WITH      
     POSITIVE MONOTONICALLY  DECREASING TERMS USING THE VAN WIJNGAARDEN       
     TRANSFORMATION OF  THE SERIES TO AN ALTERNATING SERIES.        
        
 KEYWORDS   :     
        
     SUMMATION,   
     SERIES,      
     VAN WIJNGAARDEN TRANSFORMATION.  
        
        
 SUBSECTION : EULER.        
        
 CALLING SEQUENCE   :       
        
     THE HEADING OF THE PROCEDURE IS    :       
     "REAL""PROCEDURE" EULER(AI, I, EPS, TIM);  
     "VALUE" EPS, TIM; "INTEGER" I, TIM; "REAL" AI, EPS;  
        
     EULER  :    DELIVERS THE COMPUTED SUM OF THE INFINITE SERIES   
                 A[I], I:= 0,1,...    .         
        
     THE MEANING OF THE FORMAL PARAMETERS IS:   
     AI     :    <ARITHMETIC EXPRESSION>;       
                 THE SUMMAND,         
                 THIS EXPRESSION WILL BE DEPENDENT ON THE JENSEN    
                 PARAMETER I;         
                 AI IS THE I-TH TERM OF THE SERIES (I >= 0).        
     I      :    <VARIABLE>;
                 JENSEN PARAMETER.    
     EPS,TIM:    <ARITHMETIC EXPRESSION>;       
                 THE SUMMATION IS CONTINUED UNTIL TIM SUCCESSIVE TERMS        
                 OF THE TRANSFORMED SERIES ARE IN ABSOLUTE VALUE LESS         
                 THAN EPS.  
        
        
        
1SECTION : 4.1                (DECEMBER 1975)                     PAGE 2      
        
        
        
 PROCEDURES USED    :    NONE.        
        
        
        
 REQUIRED CENTRAL MEMORY    :         
        
     EXECUTION FIELD LENGTH     :    25.        
        
        
 LANGUAGE   :    ALGOL 60.  
        
        
 METHOD AND PERFORMANCE :   
        
     EULER PERFORMS THE SUMMATION OF AN ALTERNATING SEQUENCE BY APPLYING      
     EULER'S TRANSFORMATION. BY THIS TRANSFORMATION THE SEQUENCE OF 
     TERMS IS REPLACED BY THE SEQUENCE OF MEANS OF TWO SUCCESSIVE TERMS.      
     IF NECESSARY THE NEW SEQUENCE IS AGAIN TRANSFORMED BY EULER'S  
     TRANSFORMATION. THE SUMMATION STOPS WHEN TIM SUCCESSIVE TERMS OF         
     THE (ONCE OR SEVERAL TIMES TRANSFORMED) SEQUENCE ARE IN ABSOLUTE         
     VALUE LESS THAN EPS.   
        
        
 REFERENCES :     
        
     P.NAUR, ED.  :  REVISED REPORT ON THE ALGORITHMIC LANGUAGE     
     ALGOL 60. COPENHAGEN (1964).     
        
        
 EXAMPLE OF USE     :       
        
     THE PROGRAM :
        
     "BEGIN""INTEGER" K;    
         "REAL""PROCEDURE" EULER(A, B, C, D); "CODE" 32010;         
        
         OUTPUT(61, "("+.8D"+2D")",   
             EULER((- 1) ** K / (K + 1) ** 2, K, "- 6, 100))        
     "END"        
        
     DELIVERS   : 
        
           +.82246703"+00.  
        
        
        
1SECTION : 4.1                (JULY 1974)                         PAGE 3      
        
        
        
 SUBSECTION : SUMPOSSERIES. 
        
 CALLING SEQUENCE   :       
        
     THE HEADING OF THE PROCEDURE IS    :       
        
     "REAL""PROCEDURE" SUMPOSSERIES(AI, I, MAXADDUP, MAXZERO, MAXRECURS,      
                                    MACHEXP, TIM);        
     "VALUE" MAXADDUP, MAXZERO, MAXRECURS, MACHEXP, TIM;  
     "REAL" AI, I, MAXZERO; "INTEGER" MAXADDUP, MAXRECURS, MACHEXP, TIM;      
        
     SUMPOSSERIES : DELIVERS THE COMPUTED SUM OF THE INFINITE SERIES
                    A[I] , I:= 1,2,....     .   
        
     THE MEANING OF THE FORMAL PARAMETERS IS:   
     AI :        <ARITHMETIC EXPRESSION>;       
                 THE SUMMAND,         
                 THIS EXPRESSION SHOULD BE DEPENDENT ON THE JENSEN  
                 PARAMETER I;         
                 AI IS THE I-TH TERM OF THE SERIES (I >= 1).        
     I :         <VARIABLE>;
                 JENSEN PARAMETER.    
     MAXADDUP :  <ARITHMETIC EXPRESSION>;       
                 UPPER LIMIT FOR THE NUMBER OF STRAIGHTFORWARD ADDITIONS      
                 .
     MAXZERO, TIM :         
                 <ARITHMETIC EXPRESSION>;       
                 TOLERANCES EPS AND TIM NEEDED FOR A CALL OF THE    
                 PROCEDURE EULER (THIS SECTION). MAXZERO IS ALSO USED         
                 AS A TOLERANCE FOR MAXADDUP STRAIGHTFORWARD ADDITIONS.       
     MAXRECURS : <ARITHMETIC EXPRESSION>;       
                 UPPER LIMIT FOR THE RECURSION DEPTH OF THE         
                 VAN WIJNGAARDEN TRANSFORMATIONS.         
     MACHEXP :   <ARITHMETIC EXPRESSION>;       
                 IN ORDER TO AVOID OVERFLOW AND EVALUATION OF THOSE 
                 TERMS WHICH CAN BE NEGLECTED, MACHEXP HAS TO BE THE
                 LARGEST ADMISSIBLE VALUE FOR WHICH TERMS WITH INDEX
                 K * (2 ** MACHEXP) CAN BE COMPUTED (K IS SMALL).   
                 OTHERWISE OVERFLOW MIGHT OCCUR IN COMPUTING A VALUE FOR      
                 THE JENSEN PARAMETER I, WHICH CAN BE AN UNUSUALLY  
                 HIGH POWER OF 2.     
        
        
        
 PROCEDURES USED    :    EULER = CP32010.       
        
        
 REQUIRED CENTRAL MEMORY    :         
        
     EXECUTION FIELD LENGTH :  ABOUT 1000 * RECURSION DEPTH.        
        
        
 LANGUAGE   :    ALGOL 60.  
        
        
1SECTION : 4.1                (JULY 1974)                         PAGE 4      
        
        
        
 METHOD AND PERFORMANCE :   
        
     WHEN THE TERMS AI WITH INDICES   
     MAXADDUP + 1, ... , MAXADDUP + TIM ARE ALL LESS THAN MAXZERO,  
     CONVERGENCE IS ASSUMED AND SUMPOSSERIES DELIVERS THE SUM OF THE
     SERIES BY STRAIGHTFORWARD ADDITION UNTIL TIM SERIAL TERMS ARE LESS       
     THAN MAXZERO. OTHERWISE THE VAN WIJNGAARDEN TRANSFORMATION IS  
     APPLIED, YIELDING AN ALTERNATING SERIES WHICH IS SUMMED UP WITH
     EULER'S METHOD. SINCE THE TERMS OF THIS ALTERNATING SERIES ARE 
     THEMSELVES INFINITE SERIES WITH POSITIVE TERMS, THE HERE DESCRIBED       
     PROCESS  IS RECURSIVELY CALLED FOR THE SUMMATION OF EACH       
     TERM THAT IS WANTED BY EULER'S METHOD.     
     HOWEVER, ONLY STRAIGHTFORWARD ADDITION IS APPLIED IF THE ALLOWED         
     RECURSION LEVEL (SPECIFIED BY MAXRECURS) HAS BEEN REACHED.     
     IN THE RECURSION THE PROCESS ASKS FOR TERMS AI WITH INDICES OF 
     THE TYPE J * (2 ** K), IN WHICH K CAN BE VERY LARGE. IN ORDER  
     TO AVOID OVERFLOW AN UPPER BOUND FOR K MUST BE GIVEN IN MACHEXP.         
     IF K EXCEEDS THIS BOUND THE CORRESPONDING TERM IS TAKEN TO BE ZERO.      
        
        
        
 REFERENCES :     
        
     [1] DANIEL, J.W.   :   
         SUMMATION OF A SERIES OF POSITIVE TERMS BY CONDENSATION    
         TRANSFORMATIONS. MATH. OF COMP. V.23, P.91-96 (1969).      
     [2] WIJNGAARDEN, A. VAN  :       
         COURSE SCIENTIFIC COMPUTING B, PROCESS ANALYSIS (DUTCH)    
         MATHEMATISCH CENTRUM CR-18 (1965).     
        
        
 EXAMPLE OF USE : 
        
     THE PROGRAM    :       
        
     "BEGIN""COMMENT" 730808, EXAMPLE OF THE USE OF SUMPOSSERIES;   
        "REAL""PROCEDURE" SUMPOSSERIES(A, B, C, D, E, F, G); "CODE" 
        32020;    
        
        "INTEGER" I;        
        OUTPUT(61, "("/, +.12D"+DD")",
        SUMPOSSERIES(1 / I ** 2, I, 100, "- 7, 8, 1068, 10))        
     "END"        
        
     DELIVERS   : 
        
         +.164493406604"+01 
         NUMBER OF TERMS USED   :  462,         
         RECURSION DEPTH        :    1.         
        
        
1SECTION : 4.1                (JULY 1974)                         PAGE 5      
        
        
        
 SOURCE TEXT(S) : 
0"CODE" 32010;    
 "REAL""PROCEDURE" EULER(AI, I, EPS, TIM);      
 "VALUE" EPS, TIM; "INTEGER" I, TIM; "REAL" AI, EPS;      
 "BEGIN""INTEGER" K, N, T; "REAL" MN, MP, DS, SUM; "ARRAY" M[0:15]; 
     N:= T:= 0; I:= 0; M[0]:= AI; SUM:= M[0] / 2;         
 NEXT TERM: I:= I + 1; MN:= AI;       
     "FOR" K:= 0 "STEP" 1 "UNTIL" N "DO"        
     "BEGIN" MP:= (MN + M[K]) / 2; M[K]:= MN; MN:= MP "END";        
     "IF" ABS(MN) < ABS(M[N]) & N < 15 "THEN"   
     "BEGIN" DS:= MN / 2; N:= N + 1; M[N]:= MN "END" "ELSE" DS:= MN;
     SUM:= SUM + DS; T:= "IF" ABS(DS) < EPS "THEN" T + 1 "ELSE" 0;  
     "IF" T < TIM "THEN" "GO TO" NEXT TERM;     
     EULER:= SUM  
 "END" EULER;     
         "EOP"    
0"CODE" 32020;    
 "REAL" "PROCEDURE" SUMPOSSERIES(AI, I, MAXADDUP, MAXZERO,MAXRECURS,
                                 MACHEXP, TIM); 
 "VALUE" MAXADDUP, MAXZERO, MAXRECURS, MACHEXP, TIM;      
 "REAL" AI, I, MAXZERO; "INTEGER" MAXADDUP, MAXRECURS, MACHEXP, TIM;
 "BEGIN" "INTEGER" RECURS, VL, VL2, VL4;        
        "REAL" "PROCEDURE" EULER(AI, I, EPS, TIM); "CODE" 32010;    
        
        "REAL" "PROCEDURE" SUMUP(AI, I); "REAL" AI, I;    
        "BEGIN" "INTEGER" J; "REAL" SUM, NEXTTERM;        
           I:= MAXADDUP + 1; J:= 1;   
        CHECK ADD: "IF" AI <= MAXZERO "THEN"    
           "BEGIN""IF" J < TIM "THEN" 
              "BEGIN" J:= J + 1; I:= I + 1; "GO TO" CHECK ADD "END" 
           "END""ELSE"      
           "IF" RECURS ^= MAXRECURS "THEN""GO TO" TRANSFORMSERIES;  
           SUM:= 0; I:= 0; J:= 0;     
        ADD LOOP: I:= I + 1; NEXTTERM:= AI;     
           J:= "IF" NEXTTERM <= MAXZERO "THEN" J + 1 "ELSE" 0;      
           SUM:= SUM + NEXTTERM;      
           "IF" J < TIM "THEN""GO TO" ADD LOOP; 
           SUMUP:= SUM; "GO TO" GOTSUM;         
        TRANSFORMSERIES:    
           "BEGIN""BOOLEAN" JODD; "INTEGER" J2; "ARRAY" V[1:VL];    
        
              "REAL""PROCEDURE" BJK(J, K); "VALUE" J, K; "REAL" K;  
              "INTEGER" J;  
              "BEGIN""REAL" COEFF;    
                 "IF" K > MACHEXP "THEN" BJK:= 0 "ELSE"   
                 "BEGIN" COEFF:= 2 ** (K - 1); I:= J * COEFF;       
                    BJK:= COEFF * AI  
                 "END"      
              "END" BJK;
1SECTION : 4.1                (JULY 1974)                         PAGE 6      
                                                                  ; 
        
        
        
              "REAL""PROCEDURE" VJ(J); "VALUE" J; "INTEGER" J;      
              "BEGIN""REAL" TEMP, K;  
                 "IF" JODD "THEN"     
                 "BEGIN" JODD:= "FALSE"; RECURS:= RECURS + 1;       
                    TEMP:= VJ:= SUMUP(BJK(J, K), K);      
                    RECURS:= RECURS - 1;        
                    "IF" J <= VL "THEN" V[J]:= TEMP "ELSE"
                    "IF" J <= VL2 "THEN" V[J - VL]:= TEMP 
                 "END""ELSE"
                 "BEGIN" JODD:= "TRUE"; "IF" J > VL4 "THEN"         
                    "BEGIN" RECURS:= RECURS + 1;
                       VJ:= - SUMUP(BJK(J, K), K); RECURS:= RECURS - 1        
                    "END""ELSE"       
                    "BEGIN" J2:= J2 + 1; I:= J2;
                       "IF" J > VL2 "THEN" VJ:= - (V[J2 - VL] - AI) / 2       
                       "ELSE"         
                       "BEGIN" TEMP:= V["IF" J <= VL "THEN" J "ELSE"
                          J - VL]:= (V[J2] - AI) / 2; VJ:= - TEMP   
                       "END"
                    "END"   
                 "END"      
              "END" VJ;     
        
              J2:= 0;       
              JODD:= "TRUE"; SUMUP:= EULER(VJ(J + 1), J, MAXZERO, TIM)        
           "END" TRANSFORMSERIES;     
        GOTSUM:   
        "END" SUMUP;        
        
        RECURS:= 0; VL:= 1000; VL2:= 2 * VL; VL4:= 2 * VL2;         
        SUMPOSSERIES:= SUMUP(AI, I)   
     "END" SUMPOSSERIES;    
         "EOP"    
1SECTION : 4.2.1              (JULY 1974)                         PAGE 1
  
  
  
 SECTION 4.2.1 CONTAINS TWO ALTERNATIVE PROCEDURES FOR THE COMPUTATION
 OF A DEFINITE INTEGRAL.
  
 A.  THE PROCEDURE QADRAT USES HIGH ORDER INTEGRATION RULES (UP TO 16-TH
     ORDER) AND IS APPROPRIATE FOR THE EVALUATION OVER A FINITE
     INTERVAL.
  
 B.  THE PROCEDURE  INTEGRAL  USES A  5-TH  ORDER METHOD AND CAN ALSO BE
     USED TO CALCULATE THE  INTEGRAL OVER A NUMBER OF CONSECUTIVE
     INTERVALS.  MOREOVER THE PROCEDURE CAN BE USED FOR THE COMPUTATION
     OF THE DEFINITE INTEGRAL  OVER AN INFINITE INTERVAL.
  
 FOR A COMPARISION  OF A  NUMBER OF  PROCEDURES  THAT  EVALUATE DEFINITE
 INTEGRALS : SEE REF[2].
  
 REFERENCES :
  
      [1] T.J.DEKKER AND C.J.ROOTHART.
          INTRODUCTION TO NUMERICAL ANALYSIS. (DUTCH).
          MATH. CENTRE  REPORT CR 244/74, AMSTERDAM.
      [2] C.J.ROOTHART AND H. FIOLET.
          QUDRATURE PROCEDURES.
          MATH. CENTRE  REPORT MR 137/72, AMSTERDAM.
  
  
1SECTION : 4.2.1.A            (JULY 1974)                         PAGE 1
  
  
  
 AUTHORS: C.J.ROOTHART.
  
  
 CONTRIBUTORS: P.W.HEMKER.
  
  
 INSTITUTE: MATHEMATICAL CENTRE.
  
  
 RECEIVED: 730530.
  
  
 BRIEF DESCRIPTION:
  
     QADRAT  COMPUTES  THE  DEFINITE  INTEGRAL  OF  A  FUNCTION  OF  ONE
     VARIABLE  OVER  A  FINITE  INTERVAL.
  
  
 KEYWORDS:
  
     INTEGRATION,
     QUADRATURE,
     DEFINITE INTEGRAL.
  
  
 CALLING SEQUENCE:
  
     THE HEADING OF THE PROCEDURE READS:
     "REAL" "PROCEDURE" QADRAT (X, A, B, FX, E);
     "VALUE" A, B; "REAL" X, A, B, FX; "ARRAY" E;
  
     QADRAT: DELIVERS THE  COMPUTED VALUE OF THE  DEFINITE INTEGRAL FROM
             A TO B  OF THE FUNCTION  F(X);
  
     THE MEANING OF THE FORMAL PARAMETERS IS:
     X:      <VARIABLE>;
             INTEGRATION VARIABLE; X CAN BE USED AS A JENSEN-PARAMETER
             DURING THE EVALUATIONS OF FX;
     A,B:    <ARITHMETIC EXPRESSION>;
             (A,B) DENOTES THE INTERVAL OF INTEGRATION;
             B < A  IS ALLOWED;
     FX:     <ARITHMETIC EXPRESSION>;
             FX  DENOTES  THE INTEGRAND  F(X). THIS  EXPRESSION  WILL BE
             DEPENDENT ON THE JENSEN-PARAMETER X.
     E:      <ARRAY IDENTIFIER>;
             "ARRAY" E[1:3];
             ENTRY: E[1]: THE RELATIVE ACCURACY REQUIRED;
                    E[2]: THE ABSOLUTE ACCURACY REQUIRED;
             EXIT:  E[3]: THE  NUMBER  OF  ELEMENTARY INTEGRATIONS  WITH
                    H < ABS(B-A) * E[1].
  
  
1SECTION : 4.2.1.A            (JULY 1974)                         PAGE 2
  
  
  
 PROCEDURES USED: NONE.
  
  
 REQUIRED CENTRAL MEMORY:
  
     EXECUTION FIELD LENGTH: CIRCA 16 + 9 * RECURSION DEPTH.
  
  
 RUNNING TIME: DEPENDS STRONGLY  ON THE DEFINITE  INTEGRAL  TO  COMPUTE.
  
  
 LANGUAGE: ALGOL 60.
  
  
 METHOD AND PERFORMANCE: SEE REF[1].
  
  
  
 REFERENCES:
  
     [1].C.J.ROOTHART AND H.FIOLET.
         QUADRATURE PROCEDURES.
         MATH.CENTRE, AMSTERDAM. REPORT MR 137/72.
  
  
 EXAMPLE OF USE:
  
     "BEGIN"
         "REAL" "PROCEDURE" QADRAT(X,A,B,FX,E); "CODE" 32070;
         "ARRAY" E[1:3]; "REAL" T,Q;
  
         E[1]:= E[2]:= "-9;
         Q:= QADRAT(T,0,3.141592653589  ,SIN(T),E);
         OUTPUT(61,"("/,+.15D"+3D,3B,3ZD,/")",Q,E[3]);
     "END"
  
     DELIVERS:
  
     +.200000000079740"+001      0
  
  
1SECTION : 4.2.1.A            (NOVEMBER 1976)                     PAGE 3
  
  
  
 SOURCE TEXT(S):
0"CODE" 32070;
    "REAL" "PROCEDURE" QADRAT(X, A, B, FX, E);
    "VALUE" A, B; "REAL" X, A, B, FX; "ARRAY" E;
    "BEGIN" "REAL" F0, F2, F3, F5, F6, F7, F9,
         F14, V, W, HMIN, HMAX, RE, AE;
  
       "REAL" "PROCEDURE" LINT(X0, XN, F0, F2, F3, F5, F6, F7, F9, F14);
       "REAL" X0, XN, F0, F2, F3, F5, F6, F7, F9, F14;
       "BEGIN" "REAL" H, XM, F1, F4, F8, F10, F11, F12, F13;
          XM:= (X0 + XN) / 2; H:= (XN - X0) / 32; X:= XM + 4 * H;
          F8:= FX; X:= XN - 4 * H; F11:= FX; X:= XN - 2 * H; F12:= FX;
          V:= 0.330580178199226 * F7 + 0.173485115707338 * (F6 + F8) +
          0.321105426559972*(F5 + F9) + 0.135007708341042 * (F3 + F11)
          + 0.165714514228223*(F2 + F12) + 0.393971460638127"- 1 * (F0
          + F14); X:= X0 + H; F1:= FX; X:= XN - H; F13:= FX;
          W:= 0.260652434656970 * F7 + 0.239063286684765 * (F6 + F8) +
          0.263062635477467*(F5 + F9) + 0.218681931383057 * (F3 + F11)
          + 0.275789764664284"- 1 * (F2 + F12) + 0.105575010053846* (F1
          + F13) + 0.157119426059518"- 1 * (F0 + F14);
          "IF" ABS(H) < HMIN "THEN" E[3]:= E[3] + 1;
          "IF" ABS(V - W) < ABS(W) * RE + AE "OR" ABS(H) < HMIN
          "THEN" LINT:= H * W "ELSE"
          "BEGIN" X:= X0 + 6 * H; F4:= FX; X:= XN - 6 * H; F10:= FX;
             V:= 0.245673430093324* F7 + 0.255786258286921* (F6 + F8) +
             0.228526063690406*(F5 + F9) + 0.500557131525460"- 1*(F4 +
             F10) + 0.177946487736780*(F3 + F11)+0.584014599347449"- 1
             * (F2 + F12) + 0.874830942871331"- 1 * (F1 + F13) +
             0.189642078648079"- 1 * (F0 + F14);
             LINT:= "IF" ABS(V - W) < ABS(V) * RE + AE "THEN" H * V
                                                       "ELSE"
             LINT(X0, XM, F0, F1, F2, F3, F4, F5, F6, F7) - LINT(XN,
             XM, F14, F13, F12, F11, F10, F9, F8, F7)
          "END"
       "END" LINT;
  
       HMAX:= (B - A) / 16; "IF" HMAX=0 "THEN"
       "BEGIN" QADRAT:= 0; "GOTO" RETURN "END";
       RE:= E[1]; AE:= 2 * E[2] / ABS(B - A); E[3]:= 0;
       HMIN:= ABS(B - A) * RE; X:= A; F0:= FX;
       X:= A + HMAX; F2:= FX; X:= A + 2 * HMAX; F3:= FX;
       X:= A + 4 * HMAX; F5:= FX; X:= A + 6 * HMAX; F6:= FX;
       X:= A + 8 * HMAX; F7:= FX; X:= B - 4 * HMAX; F9:= FX; X:= B;
       F14:= FX;
       QADRAT:= LINT(A, B, F0, F2, F3, F5, F6, F7, F9, F14) * 16;
       RETURN:
    "END" QADRAT;
         "EOP"
1SECTION : 4.2.1.B            (JULY 1974)                         PAGE 1
  
  
  
 AUTHOR: H.N.GLORIE.
  
  
 CONTRIBUTOR: H.FIOLET.
  
  
 INSTITUTE: MATHEMATICAL CENTRE.
  
  
 RECEIVED: 730606.
  
  
 BRIEF DESCRIPTION:
  
     INTEGRAL  CALCULATES  THE DEFINITE INTEGRAL OF A  FUNCTION  OF  ONE
     VARIABLE, OVER  A  FINITE OR INFINITE INTERVAL OR OVER A NUMBER  OF
     CONSECUTIVE INTERVALS.
  
  
 KEYWORDS:
  
     DEFINITE INTEGRAL,
     INFINITE INTERVAL,
     SIMPSON RULE,
     RICHARDSON CORRECTION.
  
  
1SECTION : 4.2.1.B            (JULY 1974)                         PAGE 2
  
  
  
 CALLING SEQUENCE:
  
     THE HEADING OF THE PROCEDURE READS:
     "REAL" "PROCEDURE" INTEGRAL(X,A,B,FX,E,UA,UB);
     "VALUE" A,B;"REAL" X,A,B,FX;"ARRAY" E;"BOOLEAN" UA,UB;
  
     INTEGRAL:
  
         DELIVERS THE  COMPUTED  VALUE  OF THE DEFINITE INTEGRAL OF  THE
         FUNCTION  FROM A TO B; AFTER SUCCESSIVE CALLS OF THE PROCEDURE,
         THE  INTEGRAL  OVER  THE TOTAL INTERVAL IS  DELIVERED, I.E. THE
         VALUE  OF  A IN  THE LAST CALL WITH UA="TRUE" IS  THE
         STARTING POINT OF THE INTERVAL.
  
     THE MEANING OF THE FORMAL PARAMETERS IS:
     X:      <VARIABLE>;
             INTEGRATION VARIABLE;X CAN BE USED AS JENSEN-PARAMETER
             FOR  FX.
     A,B:    <ARITHMETIC EXPRESSION>;
             (A,B) DENOTES THE INTERVAL OF INTEGRATION;B<A IS ALLOWED;
     FX:     <ARITHMETIC EXPRESSION>;
             THE INTEGRAND F(X);
     E:      <ARRAY IDENTIFIER>;
             "ARRAY" E[1:6];
             ENTRY : E[1]:THE RELATIVE ACCURACY REQUIRED;
                     E[2]:THE ABSOLUTE ACCURACY REQUIRED;
             EXIT:   E[3]:THE NUMBER OF SKIPPED INTEGRATION STEPS;
                     E[4]:THE VALUE OF THE INTEGRAL FROM A TO B;
                     E[5]:"IF" UB "THEN" B "ELSE" 0;
                     E[6]:"IF" UB "THEN" F(B) "ELSE" 0;
     UA:     <BOOLEAN EXPRESSION>;
             DETERMINES THE STARTING POINT OF THE INTEGRATION;
             STARTING POINT:="IF" UA "THEN" A "ELSE" E[5];
     UB:     <BOOLEAN EXPRESSION>;
             DETERMINES THE FINAL POINT OF THE INTEGRATION;
             FINAL POINT:="IF" UB "THEN" B "ELSE"
                       "IF" B>A "THEN" +INFINITY "ELSE" -INFINITY;
             IN  THE  CASE UB="FALSE" , THE VALUE OF B IS STILL RELEVANT
             (SEE METHOD AND PERFORMANCE).
  
  
 PROCEDURES USED: NONE.
  
  
 REQUIRED CENTRAL MEMORY:
  
     EXECUTION FIELD LENGTH: CIRCA 16 + 5 * RECURSION LEVEL.
  
  
 RUNNING TIME: DEPENDS STRONGLY ON THE INTEGRAL TO COMPUTE.
  
  
1SECTION : 4.2.1.B            (MARCH 1977)                        PAGE 3
  
  
  
 LANGUAGE: ALGOL 60.
  
  
 METHOD AND PERFORMANCE:
  
     INTEGRAL USES THE  SUBPROCEDURE  QAD  FOR THE  CALCULATION  OF  THE
     DEFINITE INTEGRAL OVER A  FINITE INTERVAL. THIS IS DONE BY MEANS OF
     SIMPSON 'S RULE WITH RICHARDSON CORRECTION.IF THE FOURTH DERIVATIVE
     IS TOO LARGE (AND THUS THE CORRECTION TERM) , THE TOTAL INTERVAL IS
     SPLIT INTO TWO EQUAL PARTS AND  THE  INTEGRATION PROCESS IS INVOKED
     RECURSIVELY. THIS IS DONE IN  SUCH A  WAY  THAT THE TOTAL AMOUNT OF
     RICHARDSON  CORRECTIONS  IS  SLIGHTLY SMALLER  THAN  OR  EQUAL  TO
       E[1] * ABS ( THE INTEGRAL FROM A TO B OF (FX) ) + E[2].
     THE INTEGRATION OVER AN INFINITE INTERVAL REQUIRES TWO CALLS OF THE
     PROCEDURE QAD. IN THE FIRST CALL QAD COMPUTES THE DEFINITE INTEGRAL
     FROM  A  TO  B .  IN  THE  INTERVAL  FROM  B  TO  + OR -  INFINITY
     THE   INTEGRAND  IS  TRANSFORMED  BY   MEANS  OF  THE  SUBSTITUTION
     Z=1/(X+1-B) TO THE DEFINITE INTEGRAL OF F(B-1+1/Z)/Z**2 FROM  0  TO
     1. FOR  THE  INTEGRATION  OF  A  DEFINITE  INTEGRAL  OVER A  FINITE
     INTERVAL  THE  USE OF  QADRAT IS RECOMMENDED ,ESPECIALLY  WHEN HIGH
     ACCURACY IS REQUIRED.
  
 REFERENCES:
  
     [1] T.J.DEKKER AND C.J.ROOTHART.
         INTRODUCTION TO NUMERICAL ANALYSIS.(DUTCH)
         MATH.CENTR. REPORT CR 24/71, AMSTERDAM.
     [2] C.J.ROOTHART AND H.FIOLET.
         QUADRATURE PROCEDURES.
         MATH.CENTR. REPORT MR 137/72, AMSTERDAM.
  
 EXAMPLE OF USE:
  
     "BEGIN"
         "REAL" "PROCEDURE" INTEGRAL(X,A,B,FX,E,UA,UB);"CODE" 32051;
         "ARRAY" E[1:6];"REAL" A,B,X;"BOOLEAN" UA,UB;
         UA:="TRUE";E[1]:=E[2]:="-14;
         "FOR" B:=2,4,20,100 "DO"
         "BEGIN" UB:=B<50;
             A:=INTEGRAL(X,-1,-B,10/X**2,E,UA,UB);
             OUTPUT(61,"("N,B+DDB,N,2(B+DDDB),/")",
                    A,E[3],E[4],E[5],E[6]);
             UA:="FALSE"
         "END"
     "END"
  
     DELIVERS:
     -4.9999999999999"+000   +00 -4.9999999999999   -002  +003
     -7.4999999999998"+000   +00 -7.4999999999998   -004  +001
     -9.5000000000000"+000   +00 -9.5000000000000   -020  +000
     -9.9999999999998"+000   +01 -9.9999999999998   +000  +000   .
  
  
1SECTION : 4.2.1.B            (MARCH 1977)                        PAGE 4
  
  
  
 SOURCE TEXT(S):
0"CODE" 32051;
     "REAL" "PROCEDURE" INTEGRAL(X, A, B, FX, E, UA, UB);
     "VALUE" A,B;"REAL" X, A, B, FX; "ARRAY" E; "BOOLEAN" UA, UB;
     "BEGIN"
         "REAL" "PROCEDURE" TRANSF;
         "BEGIN" Z:= 1 / X; X:= Z + B1; TRANSF:= FX * Z * Z "END";
         "REAL" "PROCEDURE" QAD(FX); "REAL" FX;
         "BEGIN" "REAL" T, V, SUM, HMIN;
             "PROCEDURE" INT;
             "BEGIN" "REAL" X3, X4, F3, F4, H;
                 X4:= X2; X2:= X1; F4:= F2; F2:= F1;
             ANEW: X:= X1:= (X0 + X2) * .5; F1:= FX;
                 X:= X3:= (X2 + X4) * .5; F3:= FX; H:= X4 - X0;
                 V:= (4 * (F1 + F3) +2  * F2 + F0 + F4) * 15;
                 T:= 6 * F2 -4  * (F1 + F3) + F0 + F4;
                 "IF" ABS(T) < ABS(V) * RE + AE "THEN"
                 SUM:=SUM + (V - T) * H "ELSE"
                 "IF" ABS(H) < HMIN "THEN" E[3]:= E[3] +1
                 "ELSE"
                 "BEGIN" INT; X2:= X3; F2:= F3; "GOTO" ANEW "END";
                 X0:= X4; F0:= F4
             "END" INT;
  
             HMIN:= ABS(X0 - X2) * RE; X:= X1:= (X0 + X2) * .5;
             F1:=FX;SUM:= 0; INT; QAD:= SUM / 180
         "END" QAD;
         "REAL" X0, X1, X2, F0, F1, F2, RE, AE, B1, Z;
         RE:= E[1]; "IF" UB "THEN" AE:= E[2] * 180 / ABS(B - A)
         "ELSE" AE:= E[2] * 90 / ABS(B - A); "IF" UA "THEN"
         "BEGIN" E[3]:= E[4]:= 0; X:= X0:= A; F0:= FX "END"
         "ELSE"
         "BEGIN" X:= X0:= A:= E[5]; F0:= E[6] "END";
         E[5]:= X:= X2:= B; E[6]:= F2:= FX; E[4]:= E[4] + QAD(FX);
         "IF"  ^UB "THEN"
         "BEGIN" "IF" A < B "THEN"
             "BEGIN" B1:= B -1 ; X0:= 1 "END"
             "ELSE"
             "BEGIN" B1:= B +1 ; X0:= -1  "END";
             F0:= E[6]; E[5]:= X2:= 0; E[6]:= F2:= 0;
             AE:= E[2] * 90;
             E[4]:= E[4] - QAD(TRANSF)
         "END";
         INTEGRAL:= E[4]
     "END" INTEGRAL;
           "EOP"
1SECTION: 0.0                 (JANUARY 1976)                      PAGE 0      
1SECTION: 0.0                 (JANUARY 1976)                      PAGE 0      
1SECTION : 5.1.2.2.1          (DECEMBER 1975)                     PAGE 1      
        
        
        
 AUTHOR: J.C.P.BUS. 
        
        
 INSTITUTE: MATHEMATICAL CENTRE.      
        
        
 RECEIVED: 730620.
        
 BRIEF DESCRIPTION: 
     THIS SECTION CONTAINS FOUR PROCEDURES, LINEMIN, RNK1UPD, DAVUPD AND      
     FLEUPD, THAT ARE AUXILIARY PROCEDURES FOR THE PROCEDURES RNK1MIN 
     AND FLEMIN (SECTION 5.1.2.2.2).  
        
        
 KEYWORDS:
     AUXILIARY PROCEDURE.   
        
        
 SUBSECTION: LINEMIN.       
        
 CALLING SEQUENCE:
        
     THE HEADING OF THIS AUXILIARY PROCEDURE IS:
     "PROCEDURE" LINEMIN(N, X, D, ND, ALFA, G, FUNCT, F0, F1, DF0, DF1,       
     EVLMAX, STRONGSEARCH, IN);       
     "VALUE" N, ND, F0, DF0, STRONGSEARCH;      
     "INTEGER" N, EVLMAX; "BOOLEAN" STRONGSEARCH; 
     "REAL" ND, ALFA, F0, F1, DF0, DF1; 
     "ARRAY" X, D, G, IN; "REAL" "PROCEDURE" FUNCT;"CODE" 34210;    
        
     THE MEANING OF THE FORMAL PARAMETERS IS:   
     N:      <ARITHMETIC EXPRESSION>; 
             THE NUMBER OF VARIABLES OF THE GIVEN FUNCTION F;       
     X:      <ARRAY IDENTIFIER>;      
             "ARRAY" X[1 : N];
             ENTRY:  A VECTOR  X0, SUCH THAT  F  IS DECREASING IN X0, IN      
                     THE DIRECTION GIVEN BY D;  
             EXIT:   THE  CALCULATED  APPROXIMATION  OF  THE  VECTOR FOR      
                     WHICH F IS MINIMAL ON THE LINE DEFINED BY:     
                     X0 + ALFA * D, (ALFA > 0); 
     D:      <ARRAY IDENTIFIER>;      
             "ARRAY" D[1 : N];
             ENTRY:  THE DIRECTION  OF THE LINE  ON WHICH  F  HAS  TO BE      
                     MINIMIZED;       
     ND:     <ARITHMETIC EXPRESSION>; 
             ENTRY:  THE EUCLIDEAN NORM OF THE VECTOR GIVEN IN D[1 : N];      
     ALFA:   <VARIABLE>;    
             THE INDEPENDENT VARIABLE, THAT DEFINES  THE POSITION ON THE      
             LINE ON WHICH F HAS TO BE MINIMIZED; 
             THIS LINE IS DEFINED BY X0 + ALFA * D, (ALFA > 0);     
             ENTRY:  AN ESTIMATE ALFA0 OF THE VALUE FOR WHICH       
                     H(ALFA) = F(X0 + ALFA * D), (ALFA > 0), IS MINIMAL;      
             EXIT:   THE CALCULATED APPROXIMATION ALFAM OF THE VALUE FOR      
                     WHICH H(ALFA) IS MINIMAL;  
1SECTION : 5.1.2.2.1          (FEBRUARY 1979)                     PAGE 2      
        
        
        
     G:      <ARRAY IDENTIFIER>;      
             "ARRAY" G[1 : N];
             EXIT:   THE GRADIENT OF  F  AT THE CALCULATED APPROXIMATION      
                     OF THE MINIMUM;  
     FUNCT:  <PROCEDURE IDENTIFIER>;  
             THE HEADING OF THIS PROCEDURE SHOULD BE:     
             "REAL" "PROCEDURE" FUNCT(N, X, G); "VALUE" N;
             "INTEGER" N; "ARRAY" X, G; 
             A CALL OF FUNCT SHOULD EFFECTUATE :
             1:  FUNCT:= F(X);
             2:  THE VALUE OF G[I], (I = 1, ..., N),  BECOMES  THE VALUE      
                 OF THE I - TH COMPONENT OF THE GRADIENT OF F AT X; 
     F0:     <ARITHMETIC EXPRESSION>; 
             ENTRY:  THE VALUE OF H(0), (SEE ALFA);       
     F1:     <VARIABLE>;    
             ENTRY:  THE VALUE OF H(ALFA0);     
             EXIT:   THE VALUE OF H(ALFAM), (SEE ALFA);   
     DF0:    <ARITHMETIC EXPRESSION>; 
             ENTRY:  THE VALUE OF THE DERIVATIVE OF H AT ALFA = 0;  
     DF1:    <VARIABLE>;    
             ENTRY:  THE VALUE OF THE DERIVATIVE OF H AT ALFA = ALFA0;
             EXIT:   THE VALUE OF THE DERIVATIVE OF H AT ALFA = ALFAM;
     EVLMAX: <VARIABLE>;    
             ENTRY:  THE MAXIMUM ALLOWED NUMBER OF CALLS OF FUNCT;  
             EXIT:   THE NUMBER OF TIMES FUNCT HAS BEEN CALLED;     
     STRONGSEARCH:
             <BOOLEAN EXPRESSION>;    
             IF  THE VALUE OF STRONGSEARCH IS  TRUE,  THEN  THE  PROCESS      
             MAKES USE OF TWO STOPPING CRITERIA:
             A:  THE NUMBER OF TIMES FUNCT HAS BEEN CALLED EXCEEDS THE
                 GIVEN VALUE OF EVLMAX; 
             B:  AN INTERVAL IS FOUND  WITH LENGTH  LESS  THAN TWO TIMES      
                 THE PRESCRIBED PRECISION,ON WICH A MINIMUM IS EXPECTED;      
             IF THE VALUE OF STRONGSEARCH  IS FALSE,  THE PROCESS  MAKES      
             ALSO USE OF A THIRD STOPPING CRITERION :     
             C:  MU <= (H(ALFAK) - H(ALFA0)) / (ALFAK * DF0) <= 1 - MU,       
                 WHEREBY ALFAK IS THE CURRENT ITERATE AND MU A      
                 PRESCRIBED CONSTANT; 
     IN:     <ARRAY IDENTIFIER>;      
             ENTRY: 
             "ARRAY" IN[1:3]; 
             IN[1]:  THE  RELATIVE  PRECISION,  EPSR,  NECESSARY FOR THE      
                     STOPPING CRITERION B, (SEE STRONGSEARCH);      
             IN[2]:  THE  ABSOLUTE  PRECISION,  EPSA,  NECESSARY FOR THE      
                     STOPPING CRITERION B, (SEE STRONGSEARCH);      
             THE PRESCRIBED PRECISION, EPS, AT ALFA = ALFAK IS GIVEN BY:      
             EPS = NORM ( X0 + ALPHA * D ) * EPSR + EPSA, WHERE     
             NORM ( . ) DENOTES THE EUCLIDEAN NORM.       
             IN[3]: THE PARAMETER MU NECESSARY FOR STOPPING CRITERION C;      
                    THIS  PARAMETER  MUST  SATISFY:   0 < MU < 0.5 ; IN       
                    PRACTICE,A CHOICE OF MU = 0.0001 IS ADVISED.    
        
        
1SECTION : 5.1.2.2.1          (FEBRUARY 1979)                     PAGE 3      
        
        
        
 DATA AND RESULTS:
        
     LINEMIN   CALCULATES   AN    APPROXIMATION   OF   A  MINIMUM  OF  A      
     HIGHER - DIMENSIONAL FUNCTION ON A GIVEN LINE;       
     THE QUANTITY DF0 MUST SATIFY: DF0 < 0;     
     IF MOREOVER DF1 > 0,  THEN  THE PROCEDURE WILL YIELD A RESULT  THAT      
     SATISFIES ONE OF THE CHOSEN STOPPING CRITERIA, (SEE STRONGSEARCH),       
     OTHERWISE WE CAN NOT GUARANTEE SUCH A RESULT.
        
        
 PROCEDURES USED: 
        
     VECVEC = CP34010,      
     ELMVEC = CP34020,      
     DUPVEC = CP31030.      
        
        
 REQUIRED CENTRAL MEMORY:   
     N WORDS.     
        
        
        
 METHOD AND PERFORMANCE:    
        
     AN APPROXIMATION TO THE MINIMUM  ON  THE  GIVEN LINE  IS CALCULATED      
     WITH CUBIC INTERPOLATION ([2]);THE STOPPING CRITERION USED WHEN THE      
     VALUE OF STRONGSEARCH IS FALSE  IS DESCRIBED  IN  [3]  AND  [4];  A      
     DETAILED DESCRIPTION OF THIS PROCEDURE IS GIVEN IN [1].
        
        
 REFERENCES:      
     [1] BUS, J. C. P.      
         MINIMIZATION OF FUNCTIONS OF SEVERAL VARIABLES (DUTCH).    
         MATHEMATICAL CENTRE, AMSTERDAM, NR 29/72 (1972). 
     [2] DAVIDON, W. C.     
         VARIABLE METRIC METHOD FOR MINIMIZATION. 
         ARGONNE NAT. LAB. REPORT, ANL 5990 (1959).       
     [3] FLETCHER, R.       
         A NEW APPROACH TO VARIABLE METRIC ALGORITHMS.    
         COMP. J. 6, (1963), P.163 - 168.       
     [4] GOLDSTEIN, A. A. AND PRICE, J. F.      
         AN EFFECTIVE ALGORITHM FOR MINIMIZATION. 
         NUMER. MATH. 10, (1967), P.184 - 189.  
        
        
1SECTION : 5.1.2.2.1          (FEBRUARY 1979)                     PAGE 4      
        
        
        
 SUBSECTION: RNK1UPD.       
        
        
 CALLING SEQUENCE:
        
     THE HEADING OF THIS AUXILIARY PROCEDURE IS:
     "PROCEDURE" RNK1UPD(H, N, V, C); "VALUE" N, C;       
     "INTEGER" N; "REAL" C; "ARRAY" H, V;       
     "CODE" 34211;
        
     THE MEANING OF THE FORMAL PARAMETERS IS:   
     N:      <ARITHMETIC EXPRESSION>; 
             THE ORDER OF THE  SYMMETRIC  MATRIX, WHOSE UPPERTRIANGLE IS      
             STORED COLUMNWISE IN THE ONE - DIMENSIONAL ARRAY H;    
     C:      <ARITHMETIC EXPRESSION>; 
             SEE V; 
     V:      <ARRAY IDENTIFIER>;      
             "ARRAY" V[1 : N];
             THE GIVEN MATRIX IS UPDATED (ANOTHER MATRIX IS ADDED TO IT)      
             WITH A SYMMETRIC MATRIX , U, OF RANK ONE, DEFINED BY:  
                        U[I,J] = C * V[I] * V[J]; 
     H:      <ARRAY IDENTIFIER>;      
             "ARRAY" H[1 : N * (N + 1) // 2];   
             ENTRY:  THE  UPPERTRIANGLE   (STORED  COLUMNWISE, I.E. : 
                     A[I,J] = H[(J-1)*J//2+I], 1 <= I <= J <= N)    
                     OF THE  SYMMETRIC MATRIX THAT HAS TO BE UPDATED; 
             EXIT:   THE  UPPERTRIANGLE   (STORED  COLUMNWISE)   OF  THE      
                     UPDATED MATRIX.  
        
        
 PROCEDURES USED: 
        
     ELMVEC = CP34020.      
        
        
 REQUIRED CENTRAL MEMORY:   
        
     NO AUXILIARY ARRAYS ARE DECLARED IN RNK1UPD. 
        
        
1SECTION : 5.1.2.2.1          (FEBRUARY 1979)                     PAGE 5      
        
        
        
 SUBSECTION: DAVUPD.
        
        
 CALLING SEQUENCE:
        
     THE HEADING OF THIS AUXILIARY PROCEDURE IS:
     "PROCEDURE" DAVUPD(H, N, V, W, C1, C2); "VALUE" N, C1, C2;     
     "INTEGER" N; "REAL" C1, C2; "ARRAY" H, V, W; 
     "CODE" 34212;
        
     THE MEANING OF THE FORMAL PARAMETERS IS:   
     N:      <ARITHMETIC EXPRESSION>; 
             THE ORDER  OF  THE  SYMMETRIC MATRIX WHOSE UPPERTRIANGLE IS      
             STORED COLUMNWISE IN THE ONE - DIMENSIONAL ARRAY H;    
     C1:     <ARITHMETIC EXPRESSION>; 
             SEE W; 
     C2:     <ARITHMETIC EXPRESSION>; 
             SEE W; 
     V:      <ARRAY IDENTIFIER>;      
             "ARRAY" V[1 : N];
             SEE W; 
     W:      <ARRAY IDENTIFIER>;      
             "ARRAY" W[1 : N];
             THE GIVEN MATRIX  IS UPDATED  WITH  A SYMMETRIC MATRIX U OF      
             RANK TWO, DEFINED BY:    
             U[I,J] = C1 * V[I] * V[J] - C2 * W[I] * W[J];
     H:      <ARRAY IDENTIFIER>;      
             "ARRAY" H[1 : N * (N + 1) // 2];   
             ENTRY:  THE UPPERTRIANGLE (STORED COLUMNWISE, I.E. :   
                     A[I,J] = H[(J - 1) * J // 2 + I], 1 <= I <= J <= N)      
                     OF THE MATRIX  THAT HAS TO BE UPDATED; 
             EXIT:   THE  UPPERTRIANGLE   (STORED  COLUMNWISE)   OF  THE      
                     UPDATED MATRIX.  
        
        
 PROCEDURES USED: NONE.     
        
        
 REQUIRED CENTRAL MEMORY:   
        
     NO AUXILIARY ARRAYS ARE DECLARED IN DAVUPD.
        
        
1SECTION : 5.1.2.2.1          (FEBRUARY 1979)                     PAGE 6      
        
        
        
 SUBSECTION: FLEUPD.
        
        
 CALLING SEQUENCE:
        
     THE HEADING OF THIS AUXILIARY PROCEDURE IS:
     "PROCEDURE" FLEUPD(H, N, V, W, C1, C2); "VALUE" N, C1, C2;     
     "INTEGER" N; "REAL" C1, C2; "ARRAY" H, V, W; 
     "CODE" 34213;
        
     THE MEANING OF THE FORMAL PARAMETERS IS:   
     N:      <ARITHMETIC EXPRESSION>; 
             THE ORDER OF THE SYMMETRIC MATRIX  WHOSE  UPPERTRIANGLE  IS      
             STORED COLUMNWISE IN THE ONE - DIMENSIONAL ARRAY H;    
     C1:     <ARITHMETIC EXPRESSION>; 
             SEE W; 
     C2:     <ARITHMETIC EXPRESSION>; 
             SEE W; 
     V:      <ARRAY IDENTIFIER>;      
             "ARRAY" V[1 : N];
             SEE W; 
     W:      <ARRAY IDENTIFIER>;      
             "ARRAY" W[1 : N];
             THE GIVEN MATRIX IS UPDATED  WITH A SYMMETRIC MATRIX  U  OF      
             RANK TWO, DEFINED BY:    
             U[I,J]= C2 * V[I] * V[J] - C1 *(V[I] * W[J] + W[I] * V[J]);      
     H:      <ARRAY IDENTIFIER>;      
             "ARRAY" H[1 : N * (N + 1) // 2];   
             ENTRY:  THE UPPERTRIANGLE (STORED COLUMNWISE, I.E. :   
                     A[I,J] = H[(J - 1) * J // 2 + I], 1 <= I <= J <= N)      
                     OF THE MATRIX THAT HAS TO BE UPDATED;
             EXIT:   THE  UPPERTRIANGLE   (STORED  COLUMNWISE)   OF  THE      
                     UPDATED MATRIX.  
        
        
 PROCEDURE USED: NONE.      
        
        
 REQUIRED CENTRAL MEMORY:   
        
     NO AUXILIARY ARRAYS ARE DECLARED IN FLEUPD.
1SECTION : 5.1.2.2.1          (DECEMBER 1975)                     PAGE 7      
        
        
        
 SOURCE TEXT(S):  
0"CODE" 34210;    
     "PROCEDURE" LINEMIN(N, X, D, ND, ALFA, G, FUNCT, F0, F1, DF0, DF1,       
     EVLMAX, STRONGSEARCH, IN); "VALUE" N, ND, F0, DF0, STRONGSEARCH; 
     "INTEGER" N, EVLMAX; "BOOLEAN" STRONGSEARCH; 
     "REAL" ND, ALFA, F0, F1, DF0, DF1; 
     "ARRAY" X, D, G, IN;   
     "REAL" "PROCEDURE" FUNCT;
     "BEGIN" "INTEGER" I, EVL;
         "BOOLEAN" NOTININT;
         "REAL" F,OLDF,DF,OLDDF,MU,ALFA0,Q,W,Y,Z,RELTOL,ABSTOL      
         ,EPS, AID; 
         "ARRAY" X0[1:N];   
         "REAL" "PROCEDURE" VECVEC(L, U, SHIFT, A, B); "CODE" 34010;
         "PROCEDURE" ELMVEC(L, U, SHIFT, A, B, X); "CODE" 34020;    
         "PROCEDURE" DUPVEC(L, U, SHIFT, A, B); "CODE" 31030;       
        
         RELTOL:= IN[1]; ABSTOL:= IN[2]; MU:= IN[3]; EVL:= 0;       
         ALFA0:= 0; OLDF:= F0; OLDDF:= DF0; Y:= ALFA; NOTININT:= "TRUE";      
         DUPVEC(1, N, 0, X0, X);      
         EPS:= (SQRT(VECVEC(1, N, 0, X, X)) * RELTOL + ABSTOL) / ND;
         Q:= (F1 - F0) / (ALFA * DF0);
     INT: "IF" NOTININT "THEN" NOTININT:= DF1 < 0 "AND" Q > MU;     
         AID:= ALFA; "IF" DF1 >= 0 "THEN"       
         "BEGIN" Z:= 3 * (OLDF - F1) / ALFA + OLDDF + DF1;
             W:= SQRT(Z ** 2 - OLDDF * DF1);    
             ALFA:= ALFA * (1 - (DF1 + W - Z) / (DF1 - OLDDF + W * 2));       
             "IF" ALFA < EPS "THEN" ALFA:= EPS "ELSE"     
             "IF" AID - ALFA < EPS "THEN" ALFA:= AID - EPS
         "END" CUBIC INTERPOLATION    
         "ELSE" "IF" NOTININT "THEN"  
         "BEGIN" ALFA0:= ALFA:= Y; OLDDF:= DF1; OLDF:= F1 "END"     
         "ELSE" ALFA:= 0.5 * ALFA; Y:= ALFA + ALFA0;      
         DUPVEC(1, N, 0, X, X0); ELMVEC(1, N, 0, X, D, Y);
         EPS:= (SQRT(VECVEC(1, N, 0, X, X)) * RELTOL + ABSTOL) / ND;
         F:= FUNCT(N, X, G); EVL:= EVL +1 ; DF:= VECVEC(1, N, 0, D, G);       
         Q:= (F - F0) / (Y * DF0);    
         "IF" ("IF" NOTININT "OR" STRONGSEARCH "THEN" "TRUE" "ELSE" 
         Q < MU "OR" Q > 1 - MU) "AND" EVL < EVLMAX "THEN"
         "BEGIN" "IF" NOTININT "OR" DF > 0 "OR" Q < MU "THEN"       
             "BEGIN" DF1:= DF; F1:= F "END"     
             "ELSE" 
             "BEGIN" ALFA0:= Y; ALFA:= AID - ALFA; OLDDF:= DF; OLDF:= F       
             "END"; 
             "IF" ALFA > EPS * 2 "THEN" "GOTO" INT
         "END";   
         ALFA:= Y; EVLMAX:= EVL; DF1:= DF; F1:= F 
     "END" LINEMIN; 
         "EOP"    
1SECTION : 5.1.2.2.1          (DECEMBER 1975)                     PAGE 8      
        
        
0"CODE" 34211;    
    "PROCEDURE" RNK1UPD(H, N, V, C);"VALUE" N, C; "INTEGER" N;      
    "REAL" C;"ARRAY" H, V;  
     "BEGIN" "INTEGER" J, K;
         "PROCEDURE" ELMVEC(L, U, SHIFT, A, B, X); "CODE" 34020;    
         K:= 0;   
         "FOR" J:= 1, J + K "WHILE" K < N "DO"  
         "BEGIN" K:= K +1 ; 
             ELMVEC(J, J + K - 1, 1 - J, H, V, V[K] * C)  
         "END"    
     "END" RNK1UPD; 
         "EOP"    
0"CODE" 34212;    
     "PROCEDURE" DAVUPD(H, N, V, W, C1, C2); "VALUE" N, C1, C2;     
     "INTEGER" N; "REAL" C1, C2; "ARRAY" H, V, W; 
     "BEGIN" "INTEGER" I, J, K;       
         "REAL" VK, WK;     
         K:= 0;   
         "FOR" J:= 1, J + K "WHILE" K < N "DO"  
         "BEGIN" K:= K +1 ; VK:= V[K] * C1; WK:= W[K] * C2; 
             "FOR" I:= 0 "STEP" 1 "UNTIL" K -1  "DO"      
             H[I + J]:= H[I + J] + V[I + 1] * VK - W[I + 1] * WK    
         "END"    
     "END" DAVUPD;
         "EOP"    
0"CODE" 34213;    
     "PROCEDURE" FLEUPD(H, N, V, W, C1, C2); "VALUE" N, C1, C2;     
     "INTEGER" N; "REAL" C1, C2; "ARRAY" H, V, W; 
     "BEGIN" "INTEGER" I, J, K;       
         "REAL" VK, WK;     
         K:= 0; "FOR" J:= 1, J + K "WHILE" K < N "DO"     
         "BEGIN" K:= K +1; VK:= - W[K] * C1 + V[K] * C2; WK:= V[K] * C1;      
             "FOR" I:= 0 "STEP" 1 "UNTIL" K - 1  "DO"     
             H[I + J]:= H[I + J] + V[I + 1] * VK -W[I + 1] * WK     
         "END"    
     "END" FLEUPD;
         "EOP"