1SECTION 3.1.2.1.1.1.2.1      (JUNE 1974)                         PAGE 1      
        
        
        
 AUTHOR: W. HOFFMANN.       
        
        
 CONTRIBUTOR: J. C. P. BUS. 
        
        
 INSTITUTE: MATHEMATICAL CENTRE.      
        
        
 RECEIVED: 731210.
        
        
 BRIEF DESCRIPTION:         
        
     THIS SECTION CONTAINS TWO  PREPARATORY  PROCEDURES  FOR THE    
     SOLUTION  OF  SYSTEMS  OF  LINEAR  ALGEBRAIC  EQUATIONS  WITH A
     TRIDIAGONAL  MATRIX;   
     DECTRI  PERFORMS  A  TRIANGULAR  DECOMPOSITION  OF  A  TRIDIAGONAL       
     MATRIX.      
     DECTRIPIV  PERFORMS A TRIANGULAR DECOMPOSITION OF A TRIDIAGONAL
     MATRIX, USING PARTIAL PIVOTING TO STABILIZE THE PROCESS.       
        
        
 KEYWORDS:        
        
     LU DECOMPOSITION,      
     TRIANGULAR DECOMPOSITION,        
     TRIDIAGONAL MATRIX.    
        
        
        
 SUBSECTION: DECTRI.        
        
        
 CALLING SEQUENCE:
        
     THE HEADING OF THIS PROCEDURE IS:
     "PROCEDURE" DECTRI(SUB, DIAG, SUPER, N, AUX);        
     "VALUE" N; "INTEGER" N; "ARRAY" SUB, DIAG, SUPER, AUX;         
        
     THE MEANING OF THE FORMAL PARAMETERS IS:   
     SUB:    <ARRAY IDENTIFIER>;      
             "ARRAY" SUB[1: N - 1];   
             ENTRY:  THE  SUBDIAGONAL   OF  THE  GIVEN  MATRIX  T,  SAY;      
                     T[I + 1, I]  SHOULD   BE  GIVEN  IN  SUB[I], I = 1,      
                     ..., N - 1;      
             EXIT:   SUPPOSE L DENOTES THE LOWER-BIDIAGONAL MATRIX, SUCH      
                     THAT  LU = T,  FOR SOME UPPER-BIDIAGONAL MATRIX  U,      
                     WITH UNIT DIAGONAL ELEMENTS, THEN L[I + 1, I]  WILL      
                     BE DELIVERED IN SUB[I], I = 1, ..., AUX[3] - 1;
1SECTION 3.1.2.1.1.1.2.1      (JUNE 1974)                         PAGE 2      
        
        
        
     DIAG:   <ARRAY IDENTIFIER>;      
             "ARRAY" DIAG[1: N];      
             ENTRY:  THE DIAGONAL OF T;         
             EXIT:   L[I, I] WILL BE DELIVERED IN  DIAG[I],  I = 1, ...,      
                     AUX[3];
     SUPER:  <ARRAY IDENTIFIER>;      
             "ARRAY" SUPER[1: N - 1]; 
             ENTRY:  THE SUPERDIAGONAL OF T; T[I, I + 1] SHOULD BE GIVEN      
                     IN SUPER[I], I = 1, ..., N - 1;      
             EXIT:   U[I, I + 1]  WILL BE DELIVERED IN  SUPER[I], I = 1,      
                     ..., AUX[3] - 1; 
     N:      <ARITHMETICAL EXPRESSION>;         
             THE ORDER OF THE MATRIX; 
     AUX:    <ARRAY IDENTIFIER>;      
             "ARRAY" AUX[2:5];        
             ENTRY :        
             AUX[2]: A RELATIVE TOLERANCE;  A REASONABLE CHOICE FOR THIS      
                     VALUE IS AN ESTIMATE OF  THE RELATIVE PRECISION  OF      
                     THE  MATRIX ELEMENTS,  HOWEVER, IT  SHOULD  NOT  BE      
                     CHOSEN SMALLER THAN THE MACHINE PRECISION;     
             EXIT:
             AUX[3]: THE NUMBER OF ELIMINATION STEPS PERFORMED;     
             AUX[5]: IF AUX[3] = N, THEN AUX[5] WILL EQUAL THE INFINITY-      
                     NORM OF THE MATRIX,  ELSE  AUX[5]  IS SET  EQUAL TO      
                     THE  VALUE  OF  THAT  ELEMENT   WHICH   CAUSES  THE      
                     BREAKDOWN OF THE DECOMPOSITION.      
        
        
 PROCEDURES USED: NONE.     
        
        
 LANGUAGE:   ALGOL 60.      
        
        
 METHOD AND PERFORMANCE:    
        
     THE METHOD USED IN DECTRI YIELDS A LOWER-BIDIAGONAL MATRIX L  AND A      
     UNIT UPPER-BIDIAGONAL MATRIX U, SUCH THAT THE PRODUCT LU EQUALS THE      
     GIVEN TRIDIAGONAL MATRIX;  THE PROCESS  IS TERMINATED  IN THE  K-TH      
     STEP, IF THE MODULUS OF THE K-TH DIAGONAL ELEMENT IS SMALLER THAN A      
     CERTAIN SMALL VALUE,  WHICH IS GIVEN BY  AUX[2]  MULTIPLIED BY  THE      
     1-NORM OF THE K-TH ROW; IN THIS CASE AUX[3] WILL BE GIVEN THE VALUE      
     K - 1  AND  AUX[5]  WILL BE GIVEN  THE VALUE  OF  THE K-TH DIAGONAL      
     ELEMENT.     
        
        
 EXAMPLE OF USE: SEE DECSOLTRI (SECTION 3.1.2.1.1.1.2.3). 
1SECTION 3.1.2.1.1.1.2.1      (DECEMBER 1975)                     PAGE 3      
        
        
        
 SUBSECTION: DECTRIPIV.     
        
 CALLING SEQUENCE:
        
     THE HEADING OF THIS PROCEDURE IS:
     "PROCEDURE" DECTRIPIV(SUB, DIAG, SUPER, N, AID, AUX, PIV);     
     "VALUE" N; "INTEGER" N; "ARRAY" SUB, DIAG, SUPER, AID, AUX;    
     "BOOLEAN" "ARRAY" PIV; 
        
     THE MEANING OF THE FORMAL PARAMETERS IS:   
     SUB:    <ARRAY IDENTIFIER>;      
             "ARRAY" SUB[1: N - 1];   
             ENTRY:  THE  SUBDIAGONAL   OF  THE  GIVEN  MATRIX  T,  SAY;      
                     T[I+1,I] SHOULD BE GIVEN IN SUB[I],I > 1,...,N - 2       
             EXIT:   LET  T'  DENOTE THE MATRIX  T  WITH  PERMUTED ROWS;      
                     SUPPOSE L DENOTES THE LOWER-BIDIAGONAL MATRIX, SUCH      
                     THAT LU = T', FOR SOME UNIT UPPER-TRIANGULAR MATRIX      
                     U, THEN  L[I + 1, I]  WILL BE DELIVERED IN  SUB[I],      
                     I = 1, ..., AUX[3] - 1;   NOTE  THAT   U   HAS  TWO      
                     CODIAGONALS, BECAUSE OF THE PARTIAL PIVOTING DURING      
                     THE DECOMPOSITION;         
     DIAG:   <ARRAY IDENTIFIER>;      
             "ARRAY" DIAG[1: N];      
             ENTRY:  THE DIAGONAL OF T;         
             EXIT:   L[I,I] WILL BE DELIVERED IN DIAG[I],I=1,...,AUX[3];      
     SUPER:  <ARRAY IDENTIFIER>;      
             "ARRAY" SUPER[1: N - 1]; 
             ENTRY:  THE SUPERDIAGONAL OF T; T[I, I + 1] SHOULD BE GIVEN      
                     IN SUPER[I], I = 1, ..., N - 1;      
             EXIT:   U[I, I + 1]  WILL BE DELIVERED IN  SUPER[I], I = 1,      
                     ..., AUX[3] - 1; 
     N:      <ARITHMETICAL EXPRESSION>;         
             THE ORDER OF THE MATRIX; 
     AID:    <ARRAY IDENTIFIER>;      
             "ARRAY" AID[1: N - 2];   
             EXIT:U[I,I+2] WILL BE DELIVERED IN AID[I],I=1,...,AUX[3]-2;      
     AUX:    <ARRAY IDENTIFIER>;      
             "ARRAY" AUX[2:5];        
             ENTRY:         
             AUX[2]: A RELATIVE TOLERANCE;  A REASONABLE CHOICE FOR THIS      
                     VALUE IS AN ESTIMATE OF  THE RELATIVE PRECISION  OF      
                     THE  MATRIX ELEMENTS,  HOWEVER, IT  SHOULD  NOT  BE      
                     CHOSEN SMALLER THAN THE MACHINE PRECISION;     
             EXIT:
             AUX[3]: THE NUMBER OF ELIMINATION STEPS PERFORMED;     
             AUX[5]: IF AUX[3] = N, THEN AUX[5] WILL EQUAL THE INFINITY-      
                     NORM OF THE MATRIX,  ELSE  AUX[5]  IS SET  EQUAL TO      
                     THE  VALUE  OF  THAT  ELEMENT   WHICH   CAUSES  THE      
                     BREAKDOWN OF THE DECOMPOSITION.      
     PIV:    <ARRAY IDENTIFIER>;      
             "BOOLEAN""ARRAY" PIV[1 : N - 1];   
             THE VALUE OF PIV[I] WILL BE TRUE IF THE I-TH AND (I + 1)-TH      
             ROW ARE INTERCHANGED, I = 1, ..., MIN(AUX[3], N - 1),  ELSE      
             PIV[I] WILL BE FALSE.    
1SECTION 3.1.2.1.1.1.2.1      (JUNE 1974)                         PAGE 4      
        
        
        
 PROCEDURES USED: NONE.     
        
 LANGUAGE:   ALGOL 60.      
        
 METHOD AND PERFORMANCE:    
        
     THE METHOD USED IN DECTRIPIV YIELDS A LOWER-BIDIAGONAL MATRIX L AND      
     A UNIT UPPER-TRIANGULAR MATRIX  U  WITH TWO CODIAGONALS,  SUCH THAT      
     THE PRODUCT  LU  EQUALS THE GIVEN TRIDIAGONAL MATRIX  WITH PERMUTED      
     ROWS; PARTIAL PIVOTING IS USED DURING THE TRIANGULAR DECOMPOSITION,      
     I.E. THAT ELEMENT OF THE K-TH COLUMN OF L IS CHOSEN AS PIVOT IN THE      
     K-TH STEP, WHOSE MODULUS DIVIDED BY THE 1-NORM OF THE CORRESPONDING      
     ROW OF THE GIVEN MATRIX IS MAXIMAL;  THE PROCESS  IS TERMINATED  IN      
     THE K-TH STEP,  IF THE MODULUS OF  THE K-TH  PIVOT ELEMENT  IS LESS      
     THAN A CERTAIN SMALL VALUE, WHICH IS GIVEN BY AUX[2]  MULTIPLIED BY      
     THE 1-NORM OF THE  CORRESPONDING ROW;  IN THIS CASE  AUX[3] WILL BE      
     GIVEN THE VALUE  K - 1, AND  AUX[5]  WILL BE GIVEN THE VALUE OF THE      
     K-TH PIVOT ELEMENT.    
        
 EXAMPLE OF USE: SEE SOLTRIPIV (SECTION 3.1.2.1.1.1.2.3). 
        
        
        
 SOURCE TEXTS:    
0"CODE" 34423;    
   "PROCEDURE" DECTRI(SUB, DIAG, SUPER, N, AUX);
   "VALUE" N; "INTEGER" N; "ARRAY" SUB, DIAG, SUPER, AUX; 
   "BEGIN" "INTEGER" I, N1; 
      "REAL" D, R, S, U, NORM, NORM1, TOL;      
      TOL:= AUX[2]; D:= DIAG[1]; R:= SUPER[1];  
      NORM:= NORM1:= ABS(D) + ABS(R); 
      "IF" ABS(D) <= NORM1 * TOL "THEN"         
      "BEGIN" AUX[3]:= 0; AUX[5]:= D; "GOTO" EXIT "END";  
      U:= SUPER[1]:= R / D; S:= SUB[1]; N1:= N - 1;       
      "FOR" I:= 2 "STEP" 1 "UNTIL" N1 "DO"      
      "BEGIN" D:= DIAG[I]; R:= SUPER[I];        
         NORM1:= ABS(S) + ABS(D) + ABS(R);      
         D:= DIAG[I]:= D - U * S;     
         "IF" ABS(D) <= NORM1 * TOL "THEN"      
         "BEGIN" AUX[3]:= I - 1; AUX[5]:= D; "GOTO" EXIT "END";     
         U:= SUPER[I]:= R / D; S:= SUB[I];      
         "IF" NORM1 > NORM "THEN" NORM:= NORM1  
      "END";      
      D:= DIAG[N]; NORM1:= ABS(D) + ABS(S);     
      D:= DIAG[N]:= D - U * S;        
      "IF" ABS(D) <= NORM1 * TOL "THEN"         
      "BEGIN" AUX[3]:= N1; AUX[5]:= D; "GOTO" EXIT "END"; 
      "IF" NORM1 > NORM "THEN" NORM:= NORM1;    
      AUX[3]:= N; AUX[5]:= NORM;      
  EXIT: 
  "END" DECTRI;   
         "EOP"    
1SECTION 3.1.2.1.1.1.2.1      (JUNE 1974)                         PAGE 5      
        
        
0"CODE" 34426;    
   "PROCEDURE" DECTRIPIV(SUB, DIAG, SUPER, N, AID, AUX, PIV);       
   "VALUE" N; "INTEGER" N; "ARRAY" SUB, DIAG, SUPER, AID, AUX;      
   "BOOLEAN" "ARRAY" PIV;   
   "BEGIN" "INTEGER" I, I1, N1, N2;   
      "REAL" D, R, S, U, T, Q, V, W, NORM, NORM1, NORM2, TOL;       
      TOL:= AUX[2]; D:= DIAG[1]; R:= SUPER[1];  
      NORM:= NORM2:= ABS(D) + ABS(R); N2:= N - 2;         
      "FOR" I:= 1 "STEP" 1 "UNTIL" N2 "DO"      
      "BEGIN" I1:= I + 1; S:= SUB[I]; T:= DIAG[I1]; Q:= SUPER[I1];  
         NORM1:= NORM2; NORM2:= ABS(S) + ABS(T) + ABS(Q); 
         "IF" NORM2 > NORM "THEN" NORM:= NORM2; 
         "IF" ABS(D) * NORM2 < ABS(S) * NORM1 "THEN"      
         "BEGIN" "IF" ABS(S) <= TOL * NORM2 "THEN"        
            "BEGIN" AUX[3]:= I - 1; AUX[5]:= S; "GOTO" EXIT "END";  
            DIAG[I]:= S; U:= SUPER[I]:= T / S;  
            V:= AID[I]:= Q / S; SUB[I]:= D;     
            W:= SUPER[I1]:= -V * D; D:= DIAG[I1]:= R - U * D;       
            R:= W; NORM2:= NORM1; PIV[I]:= "TRUE"         
         "END" "ELSE"       
         "BEGIN" "IF" ABS(D) <= TOL * NORM1 "THEN"        
            "BEGIN" AUX[3]:= I - 1; AUX[5]:= D; "GOTO" EXIT "END";  
            U:= SUPER[I]:= R / D; D:= DIAG[I1]:= T - U * S;         
            AID[I]:= 0; PIV[I]:= "FALSE"; R:= Q 
         "END"    
      "END";      
      N1:= N - 1; S:= SUB[N1]; T:= DIAG[N]; NORM1:= NORM2;
      NORM2:= ABS(S) + ABS(T); "IF" NORM2 > NORM "THEN" NORM:= NORM2;         
      "IF" ABS(D) * NORM2 < ABS(S) * NORM1 "THEN"         
      "BEGIN" "IF" ABS(S) <= TOL * NORM2 "THEN" 
         "BEGIN" AUX[3]:= N2; AUX[5]:= S; "GOTO" EXIT "END";        
         DIAG[N1]:= S; U:= SUPER[N1]:= T / S; SUB[N1]:= D;
         D:= DIAG[N]:= R - U * D; NORM2:= NORM1; PIV[N1]:= "TRUE"   
      "END" "ELSE"
      "BEGIN" "IF" ABS(D) <= TOL * NORM1 "THEN" 
         "BEGIN" AUX[3]:= N2; AUX[5]:= D; "GOTO" EXIT "END";        
         U:= SUPER[N1]:= R / D; D:= DIAG[N]:= T - U * S;  
         PIV[N1]:= "FALSE"  
      "END";      
      "IF" ABS(D) <= TOL * NORM2 "THEN"         
      "BEGIN" AUX[3]:= N1; AUX[5]:= D; "GOTO" EXIT "END"; 
      AUX[3]:= N; AUX[5]:= NORM;      
  EXIT: 
  "END" DECTRIPIV;
         "EOP"    
1SECTION 3.1.2.1.1.1.2.3      (JUNE 1974)                         PAGE 1      
        
        
        
 AUTHOR: W. HOFFMANN.       
        
        
 CONTRIBUTOR: J. C. P. BUS. 
        
        
 INSTITUTE: MATHEMATICAL CENTRE.      
        
        
 RECEIVED: 731210.
        
        
 BRIEF DESCRIPTION:         
        
     THIS SECTION CONTAINS FOUR PROCEDURES      
     FOR SOLVING A SYSTEM OF LINEAR EQUATIONS  WITH A     
     TRIDIAGONAL MATRIX;    
     SOLTRI CALCULATES A SOLUTION BY FORWARD AND BACK SUBSTITUTION IF         
     THE TRIANGULAR DECOMPOSED FORM AS DELIVERED BY DECTRI IS GIVEN.
     DECSOLTRI  PERFORMS  THE TRIANGULAR DECOMPOSITION OF THE GIVEN 
     MATRIX ( NOT USING ANY PIVOTING STRATEGY DURING THE PROCESS ) AND        
     CALCULATES THE SOLUTION BY FORWARD AND BACK SUBSTITUTION.      
     SOLTRIPIV CALCULATES A SOLUTION BY FORWARD AND BACK SUBTITUTION,IF       
     THE TRIANGULAR DECOMPOSED FORM AS DELIVERED BY DECTRIPIV IS GIVEN.       
     DECSOLTRIPIV PERFORMS THE TRIANGULAR DECOMPOSITION OF THE GIVEN
     MATRIX ( USING PARTIAL PIVOTING ) AND CALCULATES THE SOLUTION BY         
     FORWARD AND BACK SUBSTITUTION.   
        
        
 KEYWORDS:        
        
     ALGEBRAIC EQUATIONS,   
     LINEAR SYSTEMS,        
     TRIDIAGONAL MATRIX,    
     FORWARD AND BACK SUBSTITUTION.   
        
        
        
1SECTION 3.1.2.1.1.1.2.3      (JUNE 1974)                         PAGE 2      
        
        
        
 SUBSECTION: SOLTRI.        
        
        
 CALLING SEQUENCE:
        
     THE HEADING OF THIS PROCEDURE IS:
     "PROCEDURE" SOLTRI(SUB, DIAG, SUPER, N, B);
     "VALUE" N; "INTEGER" N; "ARRAY" SUB, DIAG, SUPER, B; 
        
     THE MEANING OF THE FORMAL PARAMETERS IS:   
     SUB:    <ARRAY IDENTIFIER>;      
             "ARRAY" SUB[1, N - 1];   
             ENTRY  :  THE   SUBDIAGONAL  OF  THE         
             LOWER-BIDIAGONAL MATRIX,  AS DELIVERED BY  DECTRI  (SECTION      
             3.1.2.1.1.1.2.1);        
     DIAG:   <ARRAY IDENTIFIER>;      
             "ARRAY" DIAG[1:N];       
             ENTRY :  THE DIAGONAL OF THE LOWER-
             BIDIAGONAL MATRIX, AS DELIVERED BY DECTRI;   
     SUPER:  <ARRAY IDENTIFIER>;      
             "ARRAY" SUPER[1: N - 1]; 
             ENTRY : THE SUPERDIAGONAL OF THE   
             UPPER-BIDIAGONAL MATRIX AS DELIVERED BY DECTRI;        
     N:      <ARITHMETICAL EXPRESSION>;         
             THE ORDER OF THE MATRIX; 
     B:      <ARRAY IDENTIFIER>;      
             "ARRAY" B[1:N];
             ENTRY:  THE RIGHT-HAND SIDE OF THE LINEAR SYSTEM;      
             EXIT:   THE CALCULATED SOLUTION OF THE LINEAR SYSTEM.  
        
        
 PROCEDURES USED: NONE.     
        
        
 LANGUAGE:   ALGOL 60.      
        
        
 METHOD AND PERFORMANCE:    
        
     SOLTRI   CALCULATES   THE  SOLUTION  OF  A  LINEAR  SYSTEM  WITH  A      
     TRIDIAGONAL  MATRIX,   WITH  FORWARD  AND  BACK  SUBSTITUTION;  THE      
     TRIANGULARLY DECOMPOSED FORM OF THE MATRIX,  AS PRODUCED BY  DECTRI      
     (SECTION  3.1.2.1.1.1.2.1),  SHOULD BE GIVEN;  ONE CALL  OF  DECTRI      
     FOLLOWED BY SEVERAL CALLS OF  SOLTRI  MAY BE USED  TO SOLVE SEVERAL      
     LINEAR SYSTEMS  HAVING  THE SAME TRIDIAGONAL MATRIX,  BUT DIFFERENT      
     RIGHT-HAND SIDES.      
        
        
 EXAMPLE OF USE: SEE DECSOLTRI (THIS SECTION).  
        
        
        
1SECTION 3.1.2.1.1.1.2.3      (JUNE 1974)                         PAGE 3      
        
        
        
 SUBSECTION: DECSOLTRI.     
        
        
 CALLING SEQUENCE:
        
     THE HEADING OF THIS PROCEDURE IS:
     "PROCEDURE" DECSOLTRI(SUB, DIAG, SUPER, N, AUX, B);  
     "VALUE" N; "INTEGER" N; "ARRAY" SUB, DIAG, SUPER, AUX, B;      
        
     THE MEANING OF THE FORMAL PARAMETERS IS:   
     SUB:    <ARRAY IDENTIFIER>;      
             "ARRAY" SUB[1: N - 1];   
             ENTRY:  THE  SUBDIAGONAL   OF  THE  GIVEN  MATRIX  T,  SAY;      
                     T[I + 1, I]  SHOULD   BE  GIVEN  IN  SUB[I], I = 1,      
                     ..., N - 1;      
             EXIT:   SUPPOSE L DENOTES THE LOWER-BIDIAGONAL MATRIX, SUCH      
                     THAT  LU = T,  FOR SOME UPPER-BIDIAGONAL MATRIX  U,      
                     WITH UNIT DIAGONAL ELEMENTS, THEN L[I + 1, I]  WILL      
                     BE DELIVERED IN SUB[I], I = 1, ..., AUX[3] - 1;
     DIAG:   <ARRAY IDENTIFIER>;      
             "ARRAY" DIAG[1: N];      
             ENTRY:  THE DIAGONAL OF T;         
             EXIT:   L[I, I] WILL BE DELIVERED IN  DIAG[I],  I = 1, ...,      
                     AUX[3];
     SUPER:  <ARRAY IDENTIFIER>;      
             "ARRAY" SUPER[1: N - 1]; 
             ENTRY:  THE SUPERDIAGONAL OF T; T[I, I + 1] SHOULD BE GIVEN      
                     IN SUPER[I], I = 1, ..., N - 1;      
             EXIT:   U[I, I + 1]  WILL BE DELIVERED IN  SUPER[I], I = 1,      
                     ..., AUX[3] - 1; 
     N:      <ARITHMETICAL EXPRESSION>;         
             THE ORDER OF THE MATRIX; 
     AUX:    <ARRAY IDENTIFIER>;      
             "ARRAY" AUX[2:5];        
             ENTRY :        
             AUX[2]: A RELATIVE TOLERANCE;  A REASONABLE CHOICE FOR THIS      
                     VALUE IS AN ESTIMATE OF  THE RELATIVE PRECISION  OF      
                     THE  MATRIX ELEMENTS,  HOWEVER, IT  SHOULD  NOT  BE      
                     CHOSEN SMALLER THAN THE MACHINE PRECISION;     
             EXIT :         
             AUX[3]: THE NUMBER OF ELIMINATION STEPS PERFORMED;     
             AUX[5]: IF AUX[3] = N, THEN AUX[5] WILL EQUAL THE INFINITY-      
                     NORM OF THE MATRIX  (SEE  SECTION  3.1.2.1.1.1.2.1,      
                     SUBSECTION DECTRI);        
     B:      <ARRAY IDENTIFIER>;      
             "ARRAY" B[1:N];
             ENTRY:  THE RIGHT-HAND SIDE OF THE LINEAR SYSTEM;      
             EXIT:   IF AUX[3] = N,  THEN  THE  SOLUTION  OF  THE LINEAR      
                     SYSTEM  IS  OVERWRITTEN  ON  B,  ELSE   B   REMAINS      
                     UNALTERED.       
1SECTION 3.1.2.1.1.1.2.3      (JUNE 1974)                         PAGE 4      
        
        
        
 PROCEDURES USED: 
        
     DECTRI = CP34423,      
     SOLTRI = CP34424.      
        
        
 LANGUAGE:   ALGOL 60.      
        
        
 METHOD AND PERFORMANCE:    
        
     DECSOLTRI  CALCULATES  THE  SOLUTION  OF  A  LINEAR  SYSTEM  WITH A      
     TRIDIAGONAL MATRIX; THE TRIANGULAR DECOMPOSITION IS DONE BY CALLING      
     DECTRI    (SECTION  3.1.2.1.1.1.2.1)   AND  THE  FORWARD  AND  BACK      
     SUBSTITUTION BY CALLING SOLTRI (THIS SECTION);  IF AUX[3] < N, THEN      
     THE EFFECT OF DECSOLTRI IS MERELY THAT OF DECTRI.    
        
        
 EXAMPLE OF USE:  
        
     LET  T  BE A TRIDIAGONAL MATRIX  WITH SUBDIAGONAL AND SUPERDIAGONAL      
     ELEMENTS I * 2 AND I RESPECTIVELY (I = 1, ..., N - 1), AND DIAGONAL      
     ELEMENTS I + 10 (I = 1, ..., N); LET  B  BE THE SECOND COLUMN OF T;      
     THEN  THE SOLUTION  OF THE LINEAR SYSTEM   TX = B   IS GIVEN BY THE      
     SECOND UNIT VECTOR;  BY THE FOLLOWING PROGRAM  WE  MAY  SOLVE  THIS      
     SYSTEM AND PRINT THE ERROR IN THE CALCULATED SOLUTION.         
        
     "BEGIN"      
         "PROCEDURE" DECSOLTRI(L, D, U, N, A, B); "CODE" 34425;     
         "REAL" "PROCEDURE" VECVEC(L, U, S, A, B); "CODE" 34010;    
         "INTEGER" I;       
         "ARRAY" D, SUB, SUPER, B[1:30], AUX[2:5];        
         "FOR" I:= 1 "STEP" 1 "UNTIL" 30 "DO"   
         "BEGIN" SUB[I]:= I * 2; SUPER[I]:= I; D[I]:= I + 10;       
             B[I]:= 0       
         "END"; B[1]:= 1; B[2]:= 12; B[3]:= 4;  
         AUX[2]:= "-14;     
         DECSOLTRI(SUB, D, SUPER, 30, AUX, B);  
         OUTPUT(71, "("/,"("AUX[3] AND AUX[5]:")",2(/,N)")",        
         AUX[3], AUX[5]);   
         B[2]:= B[2] - 1;   
         OUTPUT(71, "("//"("ERROR IN THE SOLUTION: ")", N")",       
         SQRT(VECVEC(1, 30, 0, B, B)))
     "END"        
        
     RESULTS:     
        
     AUX[3] AND AUX[5]:     
     +3.0000000000000"+001  
     +1.2400000000000"+002  
        
     ERROR IN THE SOLUTION: +0.0000000000000"+000         
1SECTION 3.1.2.1.1.1.2.3      (JUNE 1974)                         PAGE 5      
        
        
        
 SUBSECTION: SOLTRIPIV.     
        
        
 CALLING SEQUENCE:
        
     THE HEADING OF THIS PROCEDURE IS:
     "PROCEDURE" SOLTRIPIV(SUB, DIAG, SUPER, N, AID, PIV, B);       
     "VALUE" N; "INTEGER" N; "ARRAY" SUB, DIAG, SUPER, AID, B;      
     "BOOLEAN" "ARRAY" PIV; 
        
     THE MEANING OF THE FORMAL PARAMETERS IS:   
     SUB:    <ARRAY IDENTIFIER>;      
             "ARRAY" SUB[1, N - 1];   
             ENTRY : THE SUBDIAGONAL OF THE     
             LOWER-BIDIAGONAL MATRIX, AS DELIVERED BY DECTRIPIV (SECTION      
             3.1.2.1.1.1.2.1);        
     DIAG:   <ARRAY IDENTIFIER>;      
             "ARRAY" DIAG[1:N];       
             ENTRY: THE DIAGONAL OF THE LOWER-  
             BIDIAGONAL MATRIX, AS DELIVERED BY DECTRIPIV;
     SUPER:  <ARRAY IDENTIFIER>;      
             "ARRAY" SUPER[1: N - 1]; 
             ENTRY : THE FIRST CODIAGONAL OF    
             THE UPPER-TRIANGULAR MATRIX AS DELIVERED BY DECTRIPIV; 
     N:      <ARITHMETICAL EXPRESSION>;         
             THE ORDER OF THE MATRIX; 
     AID:    <ARRAY IDENTIFIER>;      
             "ARRAY" AID[1: N - 2];   
             ENTRY:  THE SECOND CODIAGONAL OF   
             THE UPPER-TRIANGULAR MATRIX AS DELIVERED BY DECTRIPIV; 
     PIV:    <ARRAY IDENTIFIER>;      
             "BOOLEAN""ARRAY" PIV[1: N-1];      
             ENTRY:  THE PIVOT-       
             INFORMATION AS DELIVERED BY DECTRIPIV;       
     B:      <ARRAY IDENTIFIER>;      
             "ARRAY" B[1:N];
             ENTRY:  THE RIGHT-HAND SIDE OF THE LINEAR SYSTEM;      
             EXIT:   THE CALCULATED SOLUTION OF THE LINEAR SYSTEM.  
        
        
 PROCEDURES USED: NONE.     
        
 LANGUAGE:   ALGOL 60.      
        
        
 METHOD AND PERFORMANCE:    
        
     SOLTRIPIV  CALCULATES  THE SOLUTION  OF  A  LINEAR  SYSTEM  WITH  A      
     TRIDIAGONAL  MATRIX,   WITH  FORWARD  AND  BACK  SUBSTITUTION;  THE      
     TRIANGULARLY DECOMPOSED FORM OF THE MATRIX,AS PRODUCED BY DECTRIPIV      
     (SECTION 3.1.2.1.1.1.2.1), SHOULD BE GIVEN;  ONE CALL OF  DECTRIPIV      
     FOLLOWED BY SEVERAL CALLS OF SOLTRIPIV MAY BE USED TO SOLVE SEVERAL      
     LINEAR SYSTEMS  HAVING  THE SAME TRIDIAGONAL MATRIX,  BUT DIFFERENT      
     RIGHT-HAND SIDES.      
1SECTION 3.1.2.1.1.1.2.3      (JUNE 1974)                         PAGE 6      
        
        
        
 EXAMPLE OF USE:  
        
     LET  T  BE THE MATRIX  AS GIVEN IN THE EXAMPLE OF USE OF  DECSOLTRI      
     (THIS SECTION)  AND LET B1 AND B2 BE THE SECOND AND THIRD COLUMN OF      
     T, THEN  THE SOLUTIONS OF THE LINEAR SYSTEMS  TX = B1  AND  TX = B2      
     ARE GIVEN BY THE SECOND AND THIRD UNIT VECTOR  RESPECTIVELY; IN THE      
     FOLLOWING PROGRAM  THESE SYSTEMS  ARE SOLVED  AND THE ERRORS IN THE      
     CALCULATED SOLUTIONS ARE PRINTED.
        
     "BEGIN"      
         "PROCEDURE" DECTRIPIV(L, D, U, N, A, AX, P); "CODE" 34426; 
         "PROCEDURE" SOLTRIPIV(L, D, U, N, A, P, B); "CODE" 34427;  
         "REAL" "PROCEDURE" VECVEC(L, U, S, A, B); "CODE" 34010;    
         "INTEGER" I;       
         "ARRAY" D, SUB, SUPER, AID, B1, B2[1:30], AUX[2:5];        
         "BOOLEAN" "ARRAY" PIV[1:29]; 
         "FOR" I:= 1 "STEP" 1 "UNTIL" 30 "DO"   
         "BEGIN" SUB[I]:= I * 2; SUPER[I]:= I; D[I]:= I + 10;       
             B1[I]:= B2[I]:= 0        
         "END"; B1[1]:= 1; B1[2]:= 12; B1[3]:= 4;         
         B2[2]:= 2; B2[3]:= 13; B2[4]:= 6;      
         AUX[2]:= "-14;     
         DECTRIPIV(SUB, D, SUPER, 30, AID, AUX, PIV);     
         SOLTRIPIV(SUB, D, SUPER, 30, AID, PIV, B1);      
         SOLTRIPIV(SUB, D, SUPER, 30, AID, PIV, B2);      
         OUTPUT(71, "("/,"("AUX[3] AND AUX[5]:")",2(/,N)")",        
         AUX[3], AUX[5]);   
         B1[2]:= B1[2] - 1; B2[3]:= B2[3] - 1;  
         OUTPUT(71, "("//"("ERROR IN B1: ")",N,"("ERROR IN B2: ")",N")",      
         SQRT(VECVEC(1, 30, 0, B1, B1)), SQRT(VECVEC(1, 30, 0, B2, B2)))      
     "END"        
        
     RESULTS:     
        
     AUX[3] AND AUX[5]:     
     +3.0000000000000"+001  
     +1.2400000000000"+002  
        
     ERROR IN B1: +0.0000000000000"+000         
     ERROR IN B2: +0.0000000000000"+000         
1SECTION 3.1.2.1.1.1.2.3      (JUNE 1974)                         PAGE 7      
        
        
        
 SUBSECTION: DECSOLTRIPIV.  
        
        
 CALLING SEQUENCE:
        
     THE HEADING OF THIS PROCEDURE IS:
     "PROCEDURE" DECSOLTRIPIV(SUB, DIAG, SUPER, N, AUX, B);         
     "VALUE" N; "INTEGER" N; "ARRAY" SUB, DIAG, SUPER, AUX, B;      
        
     THE MEANING OF THE FORMAL PARAMETERS IS:   
     SUB:    <ARRAY IDENTIFIER>;      
             "ARRAY" SUB[1: N - 1];   
             ENTRY:  THE  SUBDIAGONAL   OF  THE  GIVEN  MATRIX  T,  SAY;      
                     T[I + 1, I]  SHOULD   BE  GIVEN  IN  SUB[I], I = 1,      
                     ..., N - 1;      
             EXIT:   THE ELEMENTS OF SUB WILL BE OVERWRITTEN;       
     DIAG:   <ARRAY IDENTIFIER>;      
             "ARRAY" DIAG[1: N];      
             ENTRY:  THE DIAGONAL OF T;         
             EXIT:   THE ELEMENTS OF DIAG WILL BE OVERWRITTEN;      
     SUPER:  <ARRAY IDENTIFIER>;      
             "ARRAY" SUPER[1: N - 1]; 
             ENTRY:  THE SUPERDIAGONAL OF T; T[I, I + 1] SHOULD BE GIVEN      
                     IN SUPER[I], I = 1, ..., N - 1;      
             EXIT:   THE ELEMENTS OF SUPER WILL BE OVERWRITTEN;     
     N:      <ARITHMETICAL EXPRESSION>;         
             THE ORDER OF THE MATRIX; 
     AUX:    <ARRAY IDENTIFIER>;      
             "ARRAY" AUX[2:5];        
             ENTRY:         
             AUX[2]: A RELATIVE TOLERANCE;  A REASONABLE CHOICE FOR THIS      
                     VALUE IS AN ESTIMATE OF  THE RELATIVE PRECISION  OF      
                     THE  MATRIX ELEMENTS,  HOWEVER, IT  SHOULD  NOT  BE      
                     CHOSEN SMALLER THAN THE MACHINE PRECISION;     
             EXIT:
             AUX[3]: THE NUMBER OF ELIMINATION STEPS PERFORMED;     
             AUX[5]: IF AUX[3] = N, THEN AUX[5] WILL EQUAL THE INFINITY-      
                     NORM OF THE MATRIX  (SEE SECTION  3.1.2.1.1.1.2.1.,      
                     SUBSECTION DECTRIPIV);     
     B:      <ARRAY IDENTIFIER>;      
             "ARRAY" B[1 : N];        
             ENTRY:  THE RIGHT-HAND SIDE OF THE LINEAR SYSTEM;      
             EXIT:   IF  AUX[3] = N,  THEN  THE SOLUTION  OF  THE LINEAR      
                     SYSTEM WILL BE OVERWRITTEN ON B, ELSE B WILL REMAIN      
                     UNALTERED.       
1SECTION 3.1.2.1.1.1.2.3      (DECEMBER 1975)                     PAGE 8      
        
        
        
 PROCEDURES USED: NONE.     
        
        
 REQUIRED CENTRAL MEMORY:   
        
     EXECUTION FIELD LENGTH: ONE  AUXILIARY ARRAY  OF  TYPE BOOLEAN  AND      
                             ORDER N IS DECLARED IN DECSOLTRIPIV;   
        
        
 LANGUAGE:   ALGOL 60.      
        
        
 METHOD AND PERFORMANCE:    
        
     ONE CALL OF  DECSOLTRIPIV  IS EQUIVALENT WITH CALLING CONSECUTIVELY      
     DECTRIPIV  (SECTION 3.1.2.1.1.1.2.1) AND  SOLTRIPIV (THIS SECTION);      
     HOWEVER, DECSOLTRIPIV DOES NOT MAKE USE OF DECTRIPIV AND SOLTRIPIV,      
     TO SAVE MEMORY SPACE AND TIME; THIS IS ONLY  TRUE  IN THE CASE THAT      
     LINEAR SYSTEMS WITH  DIFFERENT MATRICES HAVE TO BE SOLVED; IF  
     AUX[3] < N THEN        
     DECSOLTRIPIV  IS TERMINATED PREMATURELY (SEE  DECTRIPIV  IN SECTION      
     3.1.2.1.1.1.2.1).      
        
        
 EXAMPLE OF USE:  
        
     THE SAME LINEAR SYSTEM AS GIVEN IN THE EXAMPLE OF USE OF  DECSOLTRI      
     MAY BE SOLVED WITH DECSOLTRIPIV BY THE FOLLOWING PROGRAM:      
        
     "BEGIN"      
         "PROCEDURE" DECSOLTRIPIV(L, D, U, N, A, B); "CODE" 34428;  
         "REAL" "PROCEDURE" VECVEC(L, U, S, A, B); "CODE" 34010;    
         "INTEGER" I;       
         "ARRAY" D, SUB, SUPER, B[1:30], AUX[2:5];        
         "FOR" I:= 1 "STEP" 1 "UNTIL" 30 "DO"   
         "BEGIN" SUB[I]:= I * 2; SUPER[I]:= I; D[I]:= I + 10;       
             B[I]:= 0       
         "END"; B[1]:= 1; B[2]:= 12; B[3]:= 4;  
         AUX[2]:= "-14;     
         DECSOLTRIPIV(SUB, D, SUPER, 30, AUX, B);         
         OUTPUT(71, "("/,"("AUX[3] AND AUX[5]:")",2(/,N)")",        
         AUX[3], AUX[5]);   
         B[2]:= B[2] - 1;   
         OUTPUT(71, "("//"("ERROR IN THE SOLUTION: ")", N")",       
         SQRT(VECVEC(1, 30, 0, B, B)))
     "END"        
        
     RESULTS:     
        
     AUX[3] AND AUX[5]:     
     +3.0000000000000"+001  
     +1.2400000000000"+002  
        
     ERROR IN THE SOLUTION: +0.0000000000000"+000         
1SECTION 3.1.2.1.1.1.2.3      (JUNE 1974)                         PAGE 9      
        
        
        
 SOURCE TEXTS:    
0"CODE" 34424;    
   "PROCEDURE" SOLTRI(SUB, DIAG, SUPER, N, B);  
   "VALUE" N; "INTEGER" N; "ARRAY" SUB, DIAG, SUPER, B;   
   "BEGIN" "INTEGER" I;     
      "REAL" R;   
      R:= B[1]:= B[1] / DIAG[1];      
      "FOR" I:= 2 "STEP" 1 "UNTIL" N "DO"       
      R:= B[I]:= (B[I] - SUB[I - 1] * R) / DIAG[I];       
      "FOR" I:= N - 1 "STEP" -1 "UNTIL" 1 "DO"  
      R:= B[I] := B[I] - SUPER[I] * R 
   "END" SOLTRI;  
         "EOP"    
        
0"CODE" 34425;    
  "PROCEDURE" DECSOLTRI(SUB, DIAG, SUPER, N, AUX, B);     
  "VALUE" N; "INTEGER" N; "ARRAY" SUB, DIAG, SUPER, AUX, B;         
  "BEGIN" "PROCEDURE" DECTRI(SUB, DIAG, SUPER, N, AUX); "CODE" 34423;         
          "PROCEDURE" SOLTRI( SUB, DIAG, SUPER, N, B);  "CODE" 34424;         
          DECTRI(SUB, DIAG, SUPER, N, AUX); "IF" AUX[3]= N "THEN"   
          SOLTRI(SUB, DIAG, SUPER, N, B)        
  "END" DECSOLTRI;
         "EOP"    
        
0"CODE" 34427;    
   "PROCEDURE" SOLTRIPIV(SUB, DIAG, SUPER, N, AID, PIV, B);         
   "VALUE" N; "INTEGER" N; "ARRAY" SUB, DIAG, SUPER, AID, B;        
   "BOOLEAN" "ARRAY" PIV;   
   "BEGIN" "INTEGER" I, N1; 
      "REAL" BI, BI1, R, S, T;        
      N1:= N - 1; 
      "FOR" I:= 1 "STEP" 1 "UNTIL" N1 "DO"      
      "BEGIN" "IF" PIV[I] "THEN"      
         "BEGIN" BI:= B[I+1]; BI1:= B[I] "END"  
         "ELSE"   
         "BEGIN" BI:= B[I]; BI1:= B[I+1] "END"; 
         R:= B[I]:= BI / DIAG[I];     
         B[I+1]:= BI1 - SUB[I] * R    
      "END";      
      R:= B[N]:= B[N] / DIAG[N];      
      T:= B[N1]:= B[N1] - SUPER[N1] * R;        
      "FOR" I:= N - 2 "STEP" -1 "UNTIL" 1 "DO"  
      "BEGIN" S:= R; R:= T; T:= B[I]:= B[I] - SUPER[I] * R -        
         ("IF" PIV[I] "THEN" AID[I] * S "ELSE" 0)         
      "END"       
  "END" SOLTRIPIV;
         "EOP"    
        
1SECTION 3.1.2.1.1.1.2.3      (JUNE 1974)                        PAGE 10      
        
        
0"CODE" 34428;    
   "PROCEDURE" DECSOLTRIPIV(SUB, DIAG, SUPER, N, AUX, B); 
   "VALUE" N; "INTEGER" N; "ARRAY" SUB, DIAG, SUPER, AUX, B;        
   "BEGIN" "INTEGER" I, I1, N1, N2;   
      "REAL" D, R, S, U, T, Q, V, W, NORM, NORM1, NORM2, TOL,       
      BI, BI1, BI2;         
      "BOOLEAN" "ARRAY" PIV[1:N];     
      TOL:= AUX[2]; D:= DIAG[1]; R:= SUPER[1]; BI:= B[1]; 
      NORM:= NORM2:= ABS(D) + ABS(R); N2:= N - 2;         
      "FOR" I:= 1 "STEP" 1 "UNTIL" N2 "DO"      
      "BEGIN" I1:= I + 1; S:= SUB[I]; T:= DIAG[I1]; Q:= SUPER[I1];  
         BI1:= B[I1];       
         NORM1:= NORM2; NORM2:= ABS(S) + ABS(T) + ABS(Q); 
         "IF" NORM2 > NORM "THEN" NORM:= NORM2; 
         "IF" ABS(D) * NORM2 < ABS(S) * NORM1 "THEN"      
         "BEGIN" "IF" ABS(S) <= TOL * NORM2 "THEN"        
            "BEGIN" AUX[3]:= I - 1; AUX[5]:= S; "GOTO" EXIT "END";  
            U:= SUPER[I]:= T / S; BI1:= B[I]:= BI1 / S;   
            BI:= BI - BI1 * D; V:= SUB[I]:= Q / S;        
            W:= SUPER[I1]:= -V * D; D:= DIAG[I1]:= R - U * D;       
            R:= W; NORM2:= NORM1; PIV[I]:= "TRUE"         
         "END" "ELSE"       
         "BEGIN" "IF" ABS(D) <= TOL * NORM1 "THEN"        
            "BEGIN" AUX[3]:= I - 1; AUX[5]:= D; "GOTO" EXIT "END";  
            U:= SUPER[I]:= R / D; BI:= B[I]:= BI / D;     
            BI:= BI1 - BI * S; D:= DIAG[I1]:= T - U * S;  
            PIV[I]:= "FALSE"; R:= Q   
         "END"    
      "END";      
      N1:= N - 1; S:= SUB[N1]; T:= DIAG[N]; NORM1:= NORM2; BI1:= B[N];        
      NORM2:= ABS(S) + ABS(T); "IF" NORM2 > NORM "THEN" NORM:= NORM2;         
      "IF" ABS(D) * NORM2 < ABS(S) * NORM1 "THEN"         
      "BEGIN" "IF" ABS(S) <= TOL * NORM2 "THEN" 
         "BEGIN" AUX[3]:= N2; AUX[5]:= S; "GOTO" EXIT "END";        
         U:= SUPER[N1]:= T / S; BI1:= B[N1]:= BI1 / S;    
         BI:= BI - BI1 * D; D:= R - U * D; NORM2:= NORM1  
      "END" "ELSE"
      "BEGIN" "IF" ABS(D) <= TOL * NORM1 "THEN" 
         "BEGIN" AUX[3]:= N2; AUX[5]:= D; "GOTO" EXIT "END";        
         U:= SUPER[N1]:= R / D; BI:= B[N1]:= BI / D;      
         BI:= BI1 - BI * S; D:= T - U * S       
      "END";      
      "IF" ABS(D) <= TOL * NORM2 "THEN"         
      "BEGIN" AUX[3]:= N1; AUX[5]:= D; "GOTO" EXIT "END"; 
      AUX[3]:= N; AUX[5]:= NORM;      
      BI1:= B[N]:= BI / D; BI:= B[N1]:= B[N1] - SUPER[N1] * BI1;    
      "FOR" I:= N - 2 "STEP" -1 "UNTIL" 1 "DO"  
      "BEGIN" BI2:= BI1; BI1:= BI;    
          BI:= B[I]:= B[I] - SUPER[I] * BI1 -   
          ("IF" PIV[I] "THEN" SUB[I] * BI2 "ELSE" 0)      
      "END";      
  EXIT: 
  "END" DECSOLTRIPIV;       
         "EOP"    
1SECTION 3.1.2.1.1.2.1.1      (DECEMBER 1975)                     PAGE 1      
        
        
        
 AUTHOR      :    T.J. DEKKER.        
        
        
 CONTRIBUTOR :    J. KOK.   
        
        
 INSTITUTE   :    MATHEMATICAL CENTRE.
        
        
 RECEIVED    :    731001.   
        
        
 BRIEF DESCRIPTION  :       
     THIS SECTION CONTAINS THE PROCEDURE CHLDECBND        
     FOR THE CHOLESKY DECOMPOSITION OF A SYMMETRIC POSITIVE DEFINITE
     BAND MATRIX. 
        
 KEYWORDS   :     
     LINEAR EQUATIONS,      
     CHOLESKY DECOMPOSITION,
     SYMMETRIC POSITIVE DEFINITE MATRIX,        
     BAND MATRIX. 
        
 CALLING SEQUENCE:
        
     THE HEADING OF THE PROCEDURE IS :
     "PROCEDURE" CHLDECBND(A, N, W, AUX); "VALUE" N, W; "INTEGER" N, W;       
     "ARRAY" A, AUX;        
        
     THE MEANING OF THE FORMAL PARAMETERS IS :  
         A    :  <ARRAY IDENTIFIER>;  
                 "ARRAY" A[1 : W * (N - 1) + N];
                 ENTRY  : A CONTAINS COLUMNWISE (I.E. THE(I,J)-TH   
                 ELEMENT OF THE MATRIX IS A[(J-1)*W+I], J=1,..,N,   
                 I=MAX(1,J-W),..,J) THE UPPER-TRIANGULAR BAND       
                 ELEMENTS OF THE SYMMETRIC BAND MATRIX;   
                 EXIT   : THE BAND ELEMENTS OF THE CHOLESKY         
                 MATRIX, WHICH IS AN UPPER-TRIANGULAR BAND MATRIX WITH        
                 W SUPERDIAGONALS, ARE DELIVERED COLUMNWISE IN A;   
         N    :  <ARITHMETIC EXPRESSION>;       
                 ORDER OF THE BAND MATRIX;      
         W    :  <ARITHMETIC EXPRESSION>;       
                 NUMBER OF  SUPERDIAGONALS OF THE MATRIX; 
         AUX  :  <ARRAY IDENTIFIER>;  
                 "ARRAY" AUX[2 : 3];  
                 ENTRY  : AUX[2] IS A RELATIVE TOLERANCE TO CONTROL THE       
                  CALCULATION OF THE DIAGONAL ELEMENTS OF THE       
                  CHOLESKY MATRIX (SEE METHOD AND PERFORMANCE);     
                 NORMAL EXIT  :       
                  AUX[3] = N;         
                 ABNORMAL EXIT :      
                  AUX[3] = K - 1, WHERE K IS THE INDEX OF THE DIAGONAL        
                  ELEMENT OF THE CHOLESKY MATRIX THAT CANNOT BE     
                  CALCULATED.         
1SECTION 3.1.2.1.1.2.1.1      (JUNE 1974)                         PAGE 2      
        
        
        
 PROCEDURES USED :
        
     VECVEC = CP34010.      
        
        
 RUNNING TIME :   
        
     (C1 * W + C2) * W * N; 
     THE CONSTANTS C1 AND C2 DEPEND UPON THE    
     ARITHMETIC OF THE COMPUTER.      
        
        
 LANGUAGE    :    ALGOL 60. 
        
        
 METHOD AND PERFORMANCE :   
        
     CHLDECBND PERFORMS THE CHOLESKY DECOMPOSITION OF A SYMMETRIC   
     POSITIVE DEFINITE MATRIX, WHOSE NON-ZERO ELEMENTS ARE IN BAND FORM,      
     AND WHOSE UPPER-TRIANGULAR BAND ELEMENTS ARE STORED COLUMNWISE IN        
     IONAL ARRAY. 
     THE METHOD USED IS CHOLESKY'S SQUARE ROOT METHOD. IF THE GIVEN 
     MATRIX IS POSITIVE DEFINITE, THEN THIS METHOD YIELDS AN UPPER- 
     TRIANGULAR BAND MATRIX, THE CHOLESKY MATRIX. THE NUMBER OF     
     NON-ZERO SUPERDIAGONALS OF THE GIVEN MATRIX AND ITS CHOLESKY   
     MATRIX ARE EQUAL. THE PROCESS IS COMPLETED IN N STAGES, AT EACH
     STAGE PRODUCING A ROW OF THE CHOLESKY MATRIX. HOWEVER, THE     
     PROCESS IS DISCONTINUED IF AT SOME STAGE, SAY K, THE K-TH DIAGONAL       
     ELEMENT OF THE GIVEN MATRIX MINUS THE SUM OF SQUARES OF THE    
     SUPERDIAGONAL ELEMENTS OF THE K-TH COLUMN OF THE CHOLESKY MATRIX         
     (THE SQUARE ROOT OF THIS QUANTITY BEING THE K-TH DIAGONAL ELEMENT        
     OF THE CHOLESKY MATRIX) IS EITHER NEGATIVE OR LESS THAN A GIVEN
     RELATIVE TOLERANCE (AUX[2]) TIMES THE MAXIMAL DIAGONAL ELEMENT 
     THE GIVEN MATRIX. IN THIS CASE THE GIVEN MATRIX, POSSIBLY MODIFIED       
     BY ROUNDING ERRORS, IS NOT POSITIVE DEFINITE. THIS IS INDICATED IN       
     AUX[3], BY WHICH THE VALUE K - 1 IS DELIVERED. IF THE
     DECOMPOSITION IS CARRIED OUT FULLY, AUX[3] BECOMES N.
     THE PROCEDURE DELIVERS THE BAND ELEMENTS OF THE CHOLESKY MATRIX.         
     SEE ALSO REF [1], SECTION 222.   
        
        
 REFERENCE   :    
        
     [1] DEKKER, T.J. :     
         ALGOL 60 PROCEDURES IN NUMERICAL ALGEBRA, PART 1,
         MC TRACT 22, 1968, MATHEMATISCH CENTRUM, AMSTERDAM.        
        
        
 EXAMPLE OF USE : 
        
     SEE EXAMPLE OF USE OF CHLSOLBND. 
1SECTION 3.1.2.1.1.2.1.1      (JUNE 1974)                         PAGE 3      
        
        
        
 SOURCE TEXT(S) : 
0"CODE" 34330;    
 "PROCEDURE" CHLDECBND(A, N, W, AUX); "VALUE" N, W; "INTEGER" N, W; 
 "ARRAY" A, AUX;  
 "BEGIN" "INTEGER" J, K, JMAX, KK, KJ, W1, START;         
         "REAL" R, EPS, MAX;
         "REAL""PROCEDURE" VECVEC(L, U, S, A, B); "CODE" 34010;     
         MAX:= 0; KK:= - W; W1:= W + 1;         
         "FOR" J:= 1 "STEP" 1 "UNTIL" N "DO"    
         "BEGIN" KK:= KK + W1; "IF" A[KK] > MAX "THEN" MAX:= A[KK]"END";      
         JMAX:= W; W1:= W + 1; KK:= - W; EPS:= AUX[2] * MAX;        
         "FOR" K:= 1 "STEP" 1 "UNTIL" N "DO"    
         "BEGIN""IF" K + W > N "THEN" JMAX:= JMAX - 1; KK:= KK + W1;
             START:= KK - K + 1;      
             R:= A[KK] - VECVEC("IF" K <= W1 "THEN" START "ELSE" KK - W,      
             KK - 1, 0, A, A); "IF" R <= EPS "THEN"       
             "BEGIN" AUX[3]:= K - 1; "GO TO" END "END";   
             A[KK]:= R:= SQRT(R); KJ:= KK;      
             "FOR" J:= 1 "STEP" 1 "UNTIL" JMAX "DO"       
             "BEGIN" KJ:= KJ + W;     
                 A[KJ]:= (A[KJ] - VECVEC("IF" K + J <= W1 "THEN" START        
                 "ELSE" KK - W + J, KK - 1, KJ - KK, A, A)) / R     
             "END"
         "END";   
         AUX[3]:= N;        
     END:         
     "END" CHLDECBND;       
         "EOP"    
1SECTION 3.1.2.1.1.2.1.2      (JUNE 1974)                         PAGE 1      
        
        
        
 CONTRIBUTOR :    J. KOK.   
        
        
 INSTITUTE   :    MATHEMATICAL CENTRE.
        
        
 RECEIVED    :    731001.   
        
        
 BRIEF DESCRIPTION  :       
        
     THIS SECTION CONTAINS ONE PROCEDURE,       
     CHLDETERMBND, FOR THE CALCULATION OF THE DETERMINANT OF A SYMMETRIC      
     POSITIVE DEFINITE BAND MATRIX.   
        
        
 KEY WORDS  :     
        
     DETERMINANT, 
     SYMMETRIC POSITIVE DEFINITE MATRIX,        
     BAND MATRIX. 
        
        
 CALLING SEQUENCE  :        
        
     THE HEADING OF THE PROCEDURE IS :
        
     "REAL""PROCEDURE" CHLDETERMBND(A, N, W); "VALUE" N, W;         
     "INTEGER" N, W; "ARRAY" A;       
        
     CHLDETERMBND DELIVERS THE DETERMINANT OF THE SYMMETRIC POSITIVE
     DEFINITE BAND MATRIX WHOSE CHOLESKY MATRIX IS STORED IN A.     
        
     THE MEANING OF THE FORMAL PARAMETERS IS :  
        
     A, N, W :   SEE CALLING SEQUENCE OF CHLDECBND        
                 (SECTION 3.1.2.1.1.2.1.1.);    
                 THE CONTENTS OF A ARE AS PRODUCED BY CHLDECBND OR  
                 CHLDECSOLBND (SECTION 3.1.2.1.1.2.1.3.). 
        
        
 PROCEDURES USED :    NONE. 
        
        
 RUNNING TIME :    PROPORTIONAL TO N. 
        
        
 LANGUAGE    :    ALGOL 60. 
        
        
1SECTION 3.1.2.1.1.2.1.2      (JUNE 1974)                         PAGE 2      
        
        
        
 METHOD AND PERFORMANCE :   
        
     CHLDETERMBND CAN BE CALLED AFTER CHLDECBND OR CHLDECSOLBND ONLY IF       
     THE CHOLESKY DECOMPOSITION WAS SUCCESSFUL, I.E. IF   
     AUX[3] = N.  
     THE FUNCTION VALUE OF CHLDETERMBND IS THE SQUARE OF THE DETERMINANT      
     OF THE CHOLESKY MATRIX.
     CHLDETERMBND SHOULD NOT BE CALLED WHEN OVERFLOW CAN BE EXPECTED.         
        
        
 EXAMPLE OF USE : 
        
     SEE EXAMPLES OF USE OF CHLSOLBND AND CHLDECSOLBND.   
        
        
        
        
        
        
 SOURCE TEXT(S) : 
0"CODE" 34331;    
 "REAL""PROCEDURE" CHLDETERMBND(A, N, W); "VALUE" N, W; "INTEGER" N,W;        
 "ARRAY" A;       
 "BEGIN""INTEGER" J, KK, W1; "REAL" P;
         W1:= W + 1; KK:= - W; P:= 1; 
         "FOR" J:= 1 "STEP" 1 "UNTIL" N "DO"    
         "BEGIN" KK:= KK + W1; P:= A[KK] * P "END";       
         CHLDETERMBND:= P * P         
     "END" CHLDETERMBND;    
         "EOP"    
1SECTION 3.1.2.1.1.2.1.3      (DECEMBER 1975)                     PAGE 1      
        
        
        
 AUTHOR      :    T.J. DEKKER.        
        
        
 CONTRIBUTOR :    J. KOK.   
        
        
 INSTITUTE   :    MATHEMATICAL CENTRE.
        
        
 RECEIVED    :    731001.   
        
 BRIEF DESCRIPTION  :       
     THIS SECTION CONTAINS TWO PROCEDURES.      
     A) CHLSOLBND, FOR THE SOLUTION OF ONE OR MORE SYSTEMS OF LINEAR
     EQUATIONS WITH THE SAME COEFFICIENT MATRIX, WHICH IS SYMMETRIC,
     POSITIVE DEFINITE AND IN BANDFORM, PROVIDED THAT THIS
     MATRIX HAS BEEN DECOMPOSED BY A CALL OF THE PROCEDURE CHLDECBND
     (SECTION 3.1.2.1.1.2.1.1.).      
     B) CHLDECSOLBND, FOR THE SOLUTION OF ONE SYSTEM OF LINEAR      
     EQUATIONS BY CHOLESKY'S SQUARE ROOT METHOD, PROVIDED THAT THE  
     SYMMETRIC POSITIVE DEFINITE COEFFICIENT MATRIX IS IN BAND FORM 
     AND IS STORED COLUMNWISE IN A ONE-DIMENSIONAL ARRAY. 
        
        
 KEYWORDS  :      
     LINEAR EQUATIONS,      
     CHOLESKY DECOMPOSITION,
     SYMMETRIC POSITIVE DEFINITE MATRIX,        
     BAND MATRIX. 
        
        
 SUBSECTION :    CHLSOLBND. 
        
 CALLING SEQUENCE  :        
        
        
     THE HEADING OF THE PROCEDURE IS :
     "PROCEDURE" CHLSOLBND(A, N, W, B); "VALUE" N, W; "INTEGER" N, W;         
     "ARRAY" A, B;
        
     THE MEANING OF THE FORMAL PARAMETERS IS :  
         A, N, W :   SEE CALLING SEQUENCE OF CHLDECBND,   
                 THE CONTENTS OF THE ARRAY A ARE AS PRODUCED BY     
                 CHLDECBND; 
         B  :    <ARRAY IDENTIFIER>;  
                 "ARRAY" B[1 : N];    
                 ENTRY : THE RIGHT HAND SIDE OF THE SYSTEM OF LINEAR
                 EQUATIONS; 
                 EXIT  :  THE SOLUTION OF THE SYSTEM.     
        
 PROCEDURES USED :
        
     VECVEC = CP34010,      
     SCAPRD1= CP34017.      
        
        
1SECTION 3.1.2.1.1.2.1.3      (JUNE 1974)                         PAGE 2      
        
        
        
 RUNNING TIME:    
        
     (C3 * W + C4) * N;     
     THE CONSTANTS C3 AND C4 DEPEND UPON THE    
     ARITHMETIC OF THE COMPUTER.      
        
        
 LANGUAGE    :    ALGOL 60. 
        
        
 METHOD AND PERFORMANCE :   
     THE PROCEDURE CHLSOLBND CALCULATES THE SOLUTION OF A SYSTEM OF 
     LINEAR EQUATIONS, PROVIDED THAT THE COEFFICIENT MATRIX WAS     
     DECOMPOSED BY A SUCCESSFUL CALL OF CHLDECBND (SECTION
     3.1.2.1.1.2.1.1.).     
     THE SOLUTION OF THE LINEAR SYSTEM IS OBTAINED BY CARRYING OUT THE        
     FORWARD AND BACK SUBSTITUTION WITH THE CHOLESKY MATRIX AND THE 
     RIGHT HAND SIDE. THE LATTER IS OVERWRITTEN BY THE SOLUTION.    
     THE SOLUTIONS OF SEVERAL SYSTEMS WITH THE SAME COEFFICIENT MATRIX        
     CAN BE OBTAINED BY SUCCESSIVE CALLS OF CHLSOLBND.    
        
        
 EXAMPLE OF USE : 
        
     THE FOLLOWING PROGRAM SOLVES THE SYSTEM OF SIMULTANEOUS EQUATIONS        
        
         2 * X1   -   X2                             = 1  
         -   X1 + 2 * X2   -   X3                    = 0  
                  -   X2 + 2 * X3   -   X4           = 0  
                           -   X3 + 2 * X4   -   X5  = 0  
                                    -   X4 + 2 * X5  = 1  
        
     "BEGIN""COMMENT" 730829, TEST CHLDECBND, CHLSOLBND AND         
        CHLDETERMBND;       
        "PROCEDURE" CHLDECBND(A, N, W, AUX); "CODE" 34330;
        "PROCEDURE" CHLSOLBND(A, N, W, B); "CODE" 34332;  
        "REAL""PROCEDURE" CHLDETERMBND(A, N, W); "CODE" 34331;      
        "INTEGER" I;        
        "ARRAY" SYMBAND[1 : 9], RIGHT[1 : 5], AUX[2 : 3]; 
        "FOR" I:= 1 "STEP" 1 "UNTIL" 9 "DO"     
        SYMBAND[I]:= "IF" I // 2 * 2 < I "THEN" 2 "ELSE" - 1;       
        RIGHT[1]:= RIGHT[5]:= 1;      
        "FOR" I:= 2, 3, 4 "DO" RIGHT[I]:= 0; AUX[2]:= "- 12;        
        CHLDECBND(SYMBAND, 5, 1, AUX);
        "IF" AUX[3] = 5 "THEN"        
        "BEGIN" CHLSOLBND(SYMBAND, 5, 1, RIGHT);
           OUTPUT(61, "("5(+2Z.4D2B), /"("DETERMINANT IS  ")" +.8D"+DD        
           ")", (RIGHT[I], I:= 1 : 5), CHLDETERMBND(SYMBAND, 5, 1)) 
        "END"     
     "END".       
        
     THIS PROGRAM DELIVERS: 
      +1.0000   +1.0000   +1.0000   +1.0000   +1.0000     
     DETERMINANT IS  +.60000000"+01   
1SECTION 3.1.2.1.1.2.1.3      (JUNE 1974)                         PAGE 3      
        
        
        
 SUBSECTION :    CHLDECSOLBND.        
        
        
 CALLING SEQUENCE :         
        
     THE HEADING OF THE PROCEDURE IS :
        
        
     "PROCEDURE" CHLDECSOLBND(A, N, W, AUX, B); "VALUE" N, W;       
     "INTEGER" N, W; "ARRAY" A, AUX, B;         
        
     THE MEANING OF THE FORMAL PARAMETERS IS :  
        
         A, N, W, AUX    : SEE CHLDECBND;       
         B               : SEE CHLSOLBND.       
        
        
 PROCEDURES USED :
        
     CHLDECBND = CP34330,   
     CHLSOLBND = CP34332.   
        
        
 RUNNING TIME:    
        
     (C1 * W + C5) * W * N; 
     THE CONSTANTS C1 AND C5 DEPEND UPON THE ELEMENTARY   
     ARITHMETIC OF THE COMPUTER.      
        
        
 LANGUAGE    :    ALGOL 60. 
        
        
 METHOD AND PERFORMANCE :   
        
     CHLDECSOLBND SOLVES A SYSTEM OF LINEAR EQUATIONS BY CALLING    
     CHLDECBND AND, PROVIDED THAT THE CHOLESKY DECOMPOSITION WAS    
     SUCCESSFUL, CHLSOLBND. THE COEFFICIENT MATRIX OF THIS SYSTEM HAS TO      
     BE A SYMMETRIC POSITIVE DEFINITE BAND MATRIX WHOSE UPPER-TRIANGULAR      
     BAND ELEMENTS ARE STORED COLUMNWISE IN IONAL ARRAY.  
        
        
 EXAMPLE OF USE : 
        
     THE PROGRAM  
        
     "BEGIN""COMMENT" 730829, TEST CHLDECSOLBND AND CHLDETERMBND;   
        "PROCEDURE" CHLDECSOLBND(A, N, W, AUX, B); "CODE" 34333;    
        "REAL""PROCEDURE" CHLDETERMBND(A, N, W); "CODE" 34331;      
        "INTEGER" I;        
        "ARRAY" SYMBAND[1 : 9], RIGHT[1 : 5], AUX[2 : 3]; 
1SECTION 3.1.2.1.1.2.1.3      (JUNE 1974)                         PAGE 4      
        
        
        
        "FOR" I:= 1 "STEP" 1 "UNTIL" 9 "DO"     
        SYMBAND[I]:= "IF" I // 2 * 2 < I "THEN" 2 "ELSE" - 1;       
        RIGHT[1]:= RIGHT[5]:= 1;      
        "FOR" I:= 2, 3, 4 "DO" RIGHT[I]:= 0; AUX[2]:= "- 12;        
        CHLDECSOLBND(SYMBAND, 5, 1, AUX, RIGHT);
        "IF" AUX[3] = 5 "THEN"        
        "BEGIN"   
           OUTPUT(61, "("5(+2Z.4D2B), /"("DETERMINANT IS  ")" +.8D"+DD        
           ")", (RIGHT[I], I:= 1 : 5), CHLDETERMBND(SYMBAND, 5, 1)) 
        "END"     
     "END"        
        
     WHICH SOLVES THE SAME PROBLEM AS THE PROGRAM IN THE EXAMPLE  OF USE      
     OF CHLSOLBND, DELIVERS:
        
      +1.0000   +1.0000   +1.0000   +1.0000   +1.0000     
     DETERMINANT IS  +.60000000"+01   
        
 SOURCE TEXT(S) : 
0"CODE" 34332;    
 "PROCEDURE" CHLSOLBND(A, N, W, B); "VALUE" N, W; "INTEGER" N, W;   
 "ARRAY" A, B;    
 "BEGIN" "INTEGER" I, K, IMAX, KK, W1;
         "REAL""PROCEDURE" VECVEC(L, U, S, A, B); "CODE" 34010;     
         "REAL""PROCEDURE" SCAPRD1(LA, SA, LB, SB, N, A, B);        
         "CODE" 34017;      
         KK:= - W; W1:= W + 1;        
         "FOR" K:= 1 "STEP" 1 "UNTIL" N "DO"    
         "BEGIN" KK:= KK + W1;        
             B[K]:= (B[K] - VECVEC("IF" K <= W1 "THEN" 1 "ELSE" K - W,        
             K - 1, KK - K, B, A)) / A[KK]      
         "END";   
         IMAX:= - 1;        
         "FOR" K:= N "STEP" - 1 "UNTIL" 1 "DO"  
         "BEGIN""IF" IMAX < W "THEN" IMAX:= IMAX + 1;     
             B[K]:= (B[K] - SCAPRD1(KK + W, W, K + 1, 1, IMAX, A, B)) /       
             A[KK]; KK:= KK - W1      
         "END"    
     "END" CHLSOLBND;       
         "EOP"    
0"CODE" 34333;    
 "PROCEDURE" CHLDECSOLBND(A, N, W, AUX, B); "VALUE" N, W; 
 "INTEGER" N, W; "ARRAY" A, AUX, B;   
 "BEGIN""PROCEDURE" CHLDECBND(A, N, W, AUX); "CODE" 34330;
         "PROCEDURE" CHLSOLBND(A, N, W, B); "CODE" 34332; 
         CHLDECBND(A, N, W, AUX);     
         "IF" AUX[3] = N "THEN" CHLSOLBND(A, N, W, B)     
     "END" CHLDECSOLBND;    
         "EOP"    
1SECTION 3.1.2.1.1.2.2.1      (JUNE 1974)                         PAGE 1      
        
        
        
 AUTHOR: W. HOFFMANN.       
        
 CONTRIBUTOR: J. C. P. BUS. 
        
 INSTITUTE: MATHEMATICAL CENTRE.      
        
 RECEIVED: 731215.
        
 BRIEF DESCRIPTION:         
     THIS  SECTION  CONTAINS  A  PROCEDURE,  DECSYMTRI,   TO  PERFORM  A      
     TRIANGULAR DECOMPOSITION OF A SYMMETRIC TRIDIAGONAL MATRIX.    
        
 KEYWORDS:        
        
     LU DECOMPOSITION,      
     TRIANGULAR DECOMPOSITION,        
     SYMMETRIC TRIDIAGONAL MATRIX.    
        
 CALLING SEQUENCE:
        
     THE HEADING OF THIS PROCEDURE IS:
     "PROCEDURE" DECSYMTRI(DIAG, CO, N, AUX);   
     "VALUE" N; "INTEGER" N; "ARRAY" DIAG, CO, AUX;       
        
     THE MEANING OF THE FORMAL PARAMETERS IS:   
     DIAG:   <ARRAY IDENTIFIER>;      
             "ARRAY" DIAG[1: N];      
             ENTRY:  THE DIAGONAL OF THE GIVEN MATRIX T, SAY;       
             EXIT:   SUPPOSE U DENOTES THE UNIT UPPER-BIDIAGONAL MATRIX,      
                     SUCH THAT  U'DU = T   FOR SOME  DIAGONAL MATRIX  D,      
                     WHERE U' DENOTES THE TRANSPOSED MATRIX; THEN D[I,I]      
                     WILL BE DELIVERED IN DIAG[I], I = 1, ..., AUX[3];        
     CO:     <ARRAY IDENTIFIER>;      
             "ARRAY" CO[1: N - 1];    
             ENTRY:  THE CODIAGONAL OF T; T[I, I + 1] SHOULD BE GIVEN IN      
                     CO[I], I = 1, ..., N - 1;  
             EXIT:   U[I, I + 1] WILL BE DELIVERED IN CO[I], I = 1, ...,      
                     AUX[3] - 1;      
     N:      <ARITHMETICAL EXPRESSION>;         
             THE ORDER OF THE MATRIX; 
     AUX:    <ARRAY IDENTIFIER>;      
             "ARRAY" AUX[2:5];        
             ENTRY:         
             AUX[2]: A RELATIVE TOLERANCE;  A REASONABLE CHOICE FOR THIS      
                     VALUE IS AN ESTIMATE OF  THE RELATIVE PRECISION  OF      
                     THE  MATRIX ELEMENTS,  HOWEVER, IT  SHOULD  NOT  BE      
                     CHOSEN SMALLER THAN THE MACHINE PRECISION;     
             EXIT:
             AUX[3]: THE NUMBER OF ELIMINATION STEPS PERFORMED;     
             AUX[5]: IF AUX[3] = N, THEN AUX[5] WILL EQUAL THE INFINITY-      
                     NORM OF THE MATRIX,  ELSE  AUX[5]  IS SET  EQUAL TO      
                     THE  VALUE  OF  THAT  ELEMENT   WHICH   CAUSES  THE      
                     BREAKDOWN OF THE DECOMPOSITION.      
1SECTION 3.1.2.1.1.2.2.1      (JUNE 1974)                         PAGE 2      
        
        
        
 PROCEDURES USED: NONE.     
        
        
 LANGUAGE:   ALGOL 60.      
        
        
 METHOD AND PERFORMANCE:    
        
     THE METHOD USED IN DECSYMTRI  YIELDS A UNIT UPPER-BIDIAGONAL MATRIX      
     U AND A DIAGONAL MATRIX  D, SUCH THAT THE PRODUCT  U'DU  EQUALS THE      
     GIVEN  SYMMETRIC TRIDIAGONAL MATRIX;  THE PROCESS  IS TERMINATED IN      
     THE K-TH STEP  IF  THE MODULUS  OF  THE  K-TH  DIAGONAL ELEMENT  IS      
     SMALLER THAN  A CERTAIN  SMALL  VALUE,  WHICH  IS GIVEN  BY  AUX[2]      
     MULTIPLIED BY THE 1-NORM OF THE K-TH ROW;  IN THIS CASE AUX[3] WILL      
     BE GIVEN THE VALUE  K - 1  AND  AUX[5]  WILL BE GIVEN  THE VALUE OF      
     THE K-TH DIAGONAL ELEMENT.       
        
        
 EXAMPLE OF USE: SEE DECSOLSYMTRI (SECTION 3.1.2.1.1.2.2.3).        
        
        
 SOURCE TEXT:     
0"CODE" 34420;    
       "PROCEDURE" DECSYMTRI(DIAG, CO, N, AUX); "VALUE" N; "INTEGER" N;       
       "ARRAY" DIAG, CO, AUX;         
       "BEGIN"  "INTEGER" I, N1;      
          "REAL" D, R, S, U, TOL, NORM, NORMR;  
          TOL:= AUX[2]; D:= DIAG[1]; R:= CO[1]; 
          NORM:= NORMR:= ABS(D) + ABS(R);       
          "IF" ABS(D) <= NORMR * TOL "THEN"     
          "BEGIN" AUX[3]:= 0; AUX[5]:= D; "GOTO" EXIT "END";        
          U:= CO[1]:= R / D; N1:= N - 1;        
          "FOR" I:= 2 "STEP" 1 "UNTIL" N1 "DO"  
          "BEGIN" S:= R; R:= CO[I]; D:= DIAG[I];
             NORMR:= ABS(S) + ABS(D) + ABS(R);  
             D:= DIAG[I]:= D - U * S; 
             "IF" ABS(D) <= NORMR * TOL "THEN"  
             "BEGIN" AUX[3]:= I - 1; AUX[5]:= D; "GOTO" EXIT "END"; 
             U:= CO[I]:= R / D; "IF" NORMR > NORM "THEN" NORM:= NORMR         
          "END";  
          D:= DIAG[N]; NORMR:= ABS(D) + ABS(R); 
          D:= DIAG[N]:= D - U * R;    
          "IF" ABS(D) <= NORMR * TOL "THEN"     
          "BEGIN" AUX[3]:= N1; AUX[5]:= D; "GOTO" EXIT "END";       
          "IF" NORMR > NORM "THEN" NORM:= NORMR;
          AUX[3]:= N; AUX[5]:= NORM;  
       EXIT:      
       "END" DECSYMTRI;     
         "EOP"    
1SECTION 3.1.2.1.1.2.2.3      (FEBRUARY 1979)                     PAGE 1      
        
        
        
 AUTHOR: W. HOFFMANN.       
        
        
 CONTRIBUTOR: J. C. P. BUS. 
        
        
 INSTITUTE: MATHEMATICAL CENTRE.      
        
        
 RECEIVED: 731215.
        
        
 BRIEF DESCRIPTION: 
        
     THIS SECTION CONTAINS TWO PROCEDURES       
     FOR  SOLVING  A  SYSTEM  OF  LINEAR  EQUATIONS   WITH  A  SYMMETRIC      
     TRIDIAGONAL  MATRIX;   
     SOLSYMTRI  CALCULATES  A  SOLUTION   IF  THE 
     TRIANGULARLY DECOMPOSED FORM,  AS DELIVERED BY  DECSYMTRI  (SECTION      
     3.1.2.1.1.2.2.1),  IS GIVEN;     
     DECSOLSYMTRI  PERFORMS THE TRIANGULAR      
     DECOMPOSITION  AS WELL AS  THE  FORWARD  AND  BACK SUBSTITUTION  TO      
     CALCULATE THE SOLUTION  OF THE GIVEN LINEAR SYSTEM.  
        
        
 KEYWORDS:
        
     ALGEBRAIC EQUATIONS,   
     LINEAR SYSTEMS,
     SYMMETRIC TRIDIAGONAL MATRIX,    
     FORWARD AND BACK SUBSTITUTION.   
        
        
        
1SECTION 3.1.2.1.1.2.2.3      (FEBRUARY 1979)                     PAGE 2      
        
        
        
 SUBSECTION: SOLSYMTRI.     
        
        
 CALLING SEQUENCE:
        
     THE HEADING OF THIS PROCEDURE IS:
     "PROCEDURE" SOLSYMTRI(DIAG, CO, N, B);     
     "VALUE" N; "INTEGER" N; "ARRAY" DIAG, CO, B; "CODE" 34421;     
        
     THE MEANING OF THE FORMAL PARAMETERS IS:   
     DIAG:   <ARRAY IDENTIFIER>;      
             "ARRAY" DIAG[1:N];       
             ENTRY:  THE  DIAGONAL  MATRIX,  AS 
             DELIVERED BY DECSYMTRI (SECTION 3.1.2.1.1.2.2.1);      
     CO:     <ARRAY IDENTIFIER>;      
             "ARRAY" CO[1: N - 1];    
             ENTRY:  THE CODIAGONAL OF THE UNIT 
             UPPER-BIDIAGONAL MATRIX AS DELIVERED BY DECSYMTRI;     
     N:      <ARITHMETIC EXPRESSION>; 
             THE ORDER OF THE MATRIX; 
     B:      <ARRAY IDENTIFIER>;      
             "ARRAY" B[1:N];
             ENTRY:  THE RIGHT HAND SIDE OF THE LINEAR SYSTEM;      
             EXIT:   THE CALCULATED SOLUTION OF THE LINEAR SYSTEM.  
        
        
 PROCEDURES USED: NONE.     
        
        
 METHOD AND PERFORMANCE:    
        
     SOLSYMTRI  CALCULATES  THE  SOLUTION  OF  A  LINEAR  SYSTEM  WITH A      
     SYMMETRIC TRIDIAGONAL MATRIX, WITH FORWARD  AND  BACK SUBSTITUTION;      
     THE TRIANGULARLY  DECOMPOSED  FORM  OF  THE MATRIX,  AS PRODUCED BY      
     DECSYMTRI (SECTION  3.1.2.1.1.2.2.1), SHOULD BE GIVEN;  ONE CALL OF      
     DECSYMTRI FOLLOWED BY SEVERAL CALLS OF  SOLSYMTRI  MAY BE  USED  TO      
     SOLVE SEVERAL LINEAR SYSTEMS  HAVING THE SAME SYMMETRIC TRIDIAGONAL      
     MATRIX, BUT DIFFERENT RIGHT HAND SIDES.    
        
        
 EXAMPLE OF USE: SEE DECSOLSYMTRI (THIS SECTION). 
        
        
        
1SECTION 3.1.2.1.1.2.2.3      (FEBRUARY 1979)                     PAGE 3      
        
        
        
 SUBSECTION: DECSOLSYMTRI.  
        
        
 CALLING SEQUENCE:
        
     THE HEADING OF THIS PROCEDURE IS:
     "PROCEDURE" DECSOLSYMTRI(DIAG, CO, N, AUX, B);       
     "VALUE" N; "INTEGER" N; "ARRAY" DIAG, CO, AUX, B; "CODE" 34422;
        
     THE MEANING OF THE FORMAL PARAMETERS IS:   
     DIAG:   <ARRAY IDENTIFIER>;      
             "ARRAY" DIAG[1: N];      
             ENTRY:  THE DIAGONAL OF THE GIVEN MATRIX T, SAY;       
             EXIT:   SUPPOSE U DENOTES THE UNIT UPPER-BIDIAGONAL MATRIX,      
                     SUCH THAT  U'DU = T   FOR SOME  DIAGONAL MATRIX  D,      
                     WHERE U' DENOTES THE TRANSPOSED MATRIX; THEN D[I,I]      
                     WILL BE DELIVERED IN DIAG[I], I = 1, ..., AUX[3];
     CO:     <ARRAY IDENTIFIER>;      
             "ARRAY" CO[1: N - 1];    
             ENTRY:  THE CODIAGONAL OF T; T[I, I + 1] SHOULD BE GIVEN IN      
                     CO[I], I = 1, ..., N - 1;  
             EXIT:   U[I, I + 1] WILL BE DELIVERED IN CO[I], I = 1, ...,      
                     AUX[3] - 1;      
     N:      <ARITHMETIC EXPRESSION>; 
             THE ORDER OF THE MATRIX; 
     AUX:    <ARRAY IDENTIFIER>;      
             "ARRAY" AUX[2:5];
             ENTRY: 
             AUX[2]: A RELATIVE TOLERANCE;  A REASONABLE CHOICE FOR THIS      
                     VALUE IS AN ESTIMATE OF  THE RELATIVE PRECISION  OF      
                     THE  MATRIX ELEMENTS,  HOWEVER, IT  SHOULD  NOT  BE      
                     CHOSEN SMALLER THAN THE MACHINE PRECISION;     
             EXIT:
             AUX[3]: THE NUMBER OF ELIMINATION STEPS PERFORMED;     
             AUX[5]: IF AUX[3] = N, THEN AUX[5] WILL EQUAL THE INFINITY-      
                     NORM OF THE MATRIX,  ELSE  AUX[5]  IS SET  EQUAL TO      
                     THE  VALUE  OF  THAT  ELEMENT   WHICH   CAUSES  THE      
                     BREAKDOWN OF THE DECOMPOSITION;      
     B:      <ARRAY IDENTIFIER>;      
             "ARRAY" B[1:N];
             ENTRY:  THE RIGHT HAND SIDE OF THE LINEAR SYSTEM;      
             EXIT:   IF  AUX[3] = N  THEN  THE  SOLUTION  OF  THE LINEAR      
                     SYSTEM  IS  OVERWRITTEN  ON   B,  ELSE   B  REMAINS      
                     UNALTERED.       
        
        
1SECTION 3.1.2.1.1.2.2.3      (FEBRUARY 1979)                     PAGE 4      
        
        
        
 PROCEDURES USED: 
        
     DECSYMTRI = CP34420,   
     SOLSYMTRI = CP34421.   
        
        
 METHOD AND PERFORMANCE:    
        
     DECSOLSYMTRI  CALCULATES  THE SOLUTION  OF  A LINEAR SYSTEM  WITH A      
     SYMMETRIC TRIDIAGONAL MATRIX;  THE TRIANGULAR DECOMPOSITION IS DONE      
     BY CALLING DECSYMTRI (SECTION 3.1.2.1.1.2.2.1)  AND THE FORWARD AND      
     BACK SUBSTITUTION BY CALLING SOLSYMTRI (THIS SECTION); IF AUX[3]<N,      
     THEN THE EFFECT OF DECSOLSYMTRI IS MERELY THAT OF DECSYMTRI.   
        
        
 EXAMPLE OF USE:  
        
     LET T BE A SYMMETRIC TRIDIAGONAL MATRIX OF ORDER 100  WITH DIAGONAL      
     ELEMENTS I (I = 1, ..., 100) AND CODIAGONAL ELEMENTS I * 2  (I = 1,      
     ..., 99); LET THE RIGHT HAND SIDE  B  BE GIVEN BY THE SECOND COLUMN      
     OF T; THEN THE SOLUTION OF THE LINEAR SYSTEM TX = B IS GIVEN BY THE      
     SECOND UNIT VECTOR;  BY  THE FOLLOWING PROGRAM  WE  MAY SOLVE  THIS      
     SYSTEM AND PRINT THE ERROR IN THE CALCULATED SOLUTION. 
        
     "BEGIN"      
         "PROCEDURE" DECSOLSYMTRI(D, C, N, A, B); "CODE" 34422;     
         "REAL" "PROCEDURE" VECVEC(L, U, S, A, B); "CODE" 34010;    
         "INTEGER" I; "ARRAY" D, CO, B[1:100], AUX[2:5];  
         "FOR" I:= 1 "STEP" 1 "UNTIL" 100 "DO"  
         "BEGIN" D[I]:= I; CO[I]:= I * 2; B[I]:= 0 "END"; 
         B[1]:= B[2]:= 2; B[3]:= 4;   
         AUX[2]:= "-14;     
         DECSOLSYMTRI(D, CO, 100, AUX, B);      
         B[2]:= B[2] - 1;   
         OUTPUT(71, "("/,"("    AUX[3] AND AUX[5]:")",2(/4B,N)")",  
         AUX[3], AUX[5]);   
         OUTPUT(71, "("//,"("    ERROR IN THE SOLUTION:")",N,/")",  
         SQRT(VECVEC(1, 100, 0, B, B))) 
     "END"
        
     RESULTS:     
        
     AUX[3] AND AUX[5]:     
     +1.0000000000000"+002  
     +4.9300000000000"+002  
        
     ERROR IN THE SOLUTION:+0.0000000000000"+000
1SECTION 3.1.2.1.1.2.2.3      (JUNE 1974)                         PAGE 5      
        
        
        
 SOURCE TEXTS:    
0"CODE" 34421;    
       "PROCEDURE" SOLSYMTRI(DIAG, CO, N, B); "VALUE" N; "INTEGER" N; 
       "ARRAY" DIAG, CO, B; 
       "BEGIN" "INTEGER" I; 
          "REAL" R, S;      
          R:= B[1]; B[1]:= R / DIAG[1]; 
          "FOR" I:= 2 "STEP" 1 "UNTIL" N "DO"   
          "BEGIN" R:= B[I] - CO[I-1] * R; B[I]:= R / DIAG[I] "END"; 
          S:= B[N]; 
          "FOR" I:= N - 1 "STEP" -1 "UNTIL" 1 "DO"
          S:= B[I]:= B[I] - CO[I] * S 
       "END" SOLSYMTRI;     
         "EOP"    
        
0"CODE" 34422;    
       "PROCEDURE" DECSOLSYMTRI(DIAG, CO, N, AUX, B); "VALUE" N;    
       "INTEGER" N; "ARRAY" DIAG, CO, AUX, B;   
       "BEGIN" "PROCEDURE" DECSYMTRI(DIAG, CO, N, AUX); "CODE" 34420; 
               "PROCEDURE" SOLSYMTRI(DIAG, CO, N, B);   "CODE" 34421; 
               DECSYMTRI(DIAG, CO, N, AUX); "IF" AUX[3] = N "THEN"  
               SOLSYMTRI(DIAG, CO, N, B)
       "END" DECSOLSYMTRI;  
         "EOP"    
1SECTION 3.1.2.2.1            (DECEMBER 1975)                     PAGE 1
 
 
 
 AUTHOR: P.W.HEMKER.
 
 
 INSTITUTE: MATHEMATICAL CENTRE.
 
 
 RECEIVED: 730615.
 
 
 BRIEF DESCRIPTION:
 
     CONJ GRAD  SOLVES A  LINEAR  SYSTEM OF  EQUATIONS BY THE  METHOD OF
     CONJUGATE GRADIENTS. THE SYSTEM  HAS TO BE  SYMMETRIC AND  POSITIVE
     DEFINITE.
 
 
 CALLING SEQUENCE:
 
     THE HEADING OF THE PROCEDURE READS:
     "PROCEDURE" CONJ GRAD (MATVEC, X, R, L, N, GO ON, ITERATE, NORM2);
     "VALUE" L, N; "BOOLEAN" GO ON; "INTEGER" L, N, ITERATE;
     "REAL" NORM2; "ARRAY" X, R; "PROCEDURE" MATVEC;
 
     THE MEANING OF THE FORMAL PARAMETERS IS:
     MATVEC: <PROCEDURE IDENTIFIER>;
             THE HEADING OF THIS PROCEDURE READS:
             "PROCEDURE" MATVEC( P, Q); "ARRAY" P, Q;
             THIS PROCEDURE DEFINES THE MATRIX A (THE COEFFICIENT MATRIX
             OF THE SYSTEM) AS FOLLOWS :
             AT EACH CALL MATVEC DELIVERS IN Q THE MATRIX-VECTOR PRODUCT
             AP; P AND Q ARE ONE - DIMENSIONAL ARRAYS:
             "ARRAY" P,Q[L:N];
     X:      <ARRAY IDENTIFIER>;
             "ARRAY" X[L:N];
             ENTRY: AN INITIAL APPROXIMATION TO THE SOLUTION X;
             EXIT:  THE SOLUTION;
     R:      <ARRAY IDENTIFIER>;
             "ARRAY" R[L:N];
             ENTRY: THE RIGHT-HAND SIDE OF THE SYSTEM;
             EXIT:  THE RESIDUE B - AX, COMPUTED RECURSIVELY;
     L,N:    <ARITHMETIC EXPRESSION>;
             L AND N ARE RESPECTIVELY THE  LOWER AND  UPPER BOUND OF THE
             ARRAYS  X,R,P,Q;
     GO ON:  <BOOLEAN EXPRESSION>;
             GO ON INDICATES THE CONTINUATION OF THE PROCESS.
             THIS EXPRESSION MAY DEPEND ON THE JENSEN PARAMETERS ITERATE
             AND NORM2. WITH THIS BOOLEAN EXPRESSION THE USER CONTROLS
             THE  CONTINUATION OF THE  PROCESS. IF  GO ON:= "FALSE"  THE
             ITERATION PROCESS IS STOPPED.
     ITERATE:<IDENTIFIER>;
             DELIVERS THE NUMBER OF ITERATION STEPS ALREADY PERFORMED;
     NORM2:  <IDENTIFIER>;
             DELIVERS THE SQUARED EUCLIDIC NORM OF THE RESIDUE,COMPUTED
             RECURSIVELY
1SECTION 3.1.2.2.1            (JUNE 1974)                         PAGE 2
 
 
 
 PROCEDURES USED:
 
     VECVEC = CP34010 ,
     ELMVEC = CP34020 .
 
 
 REQUIRED CENTRAL MEMORY:
 
     EXECUTION FIELD LENGTH: 7 + 2 * ( N - L + 1 ).
 
 
 RUNNING TIME:
 
     THE  RUNNING TIME IS PROPORTIONAL TO THE  NUMBER OF ITERATION STEPS
     PERFORMED. EACH  ITERATION STEP  REQUIRES  ONE  EVALUATION  OF  THE
     PROCEDURE MATVEC, THE  EVALUATION OF TWO SCALAR - VECTOR - PRODUCTS
     AND ONE VECTOR - VECTOR - PRODUCT.
 
 
 
 LANGUAGE: ALGOL 60.
 
 
 METHOD AND PERFORMANCE: SEE REF[1].
 
 
 REFERENCES:
 
     [1].J.K.REID.
         ON THE METHOD OF  CONJUGATE GRADIENTS FOR THE SOLUTION OF LARGE
         SPARSE SYSTEMS OF LINEAR EQUATIONS.
         IN:LARGE SPARSE SETS OF LINEAR EQUATIONS (J.K.REID ED)1971.
 
 
1SECTION 3.1.2.2.1            (DECEMBER 1979)                     PAGE 3
 
 
 
 EXAMPLE OF USE:
 
     "BEGIN"
         "PROCEDURE" CONJ GRAD(A,X,R,L,N,GO ON,IT,NO);"CODE" 34220;
         "ARRAY" X,B[0:12];
         "INTEGER" IT,I;
         "REAL" NO;
         "PROCEDURE" A(X,B);
         "ARRAY" X,B;
         "BEGIN" B[0]:=2*X[0]-X[1];
             "FOR" I:=1 "STEP" 1 "UNTIL" 11 "DO"
             B[I]:=-X[I-1]+2*X[I]-X[I+1];
             B[12]:=2*X[12]-X[11]
         "END" A;
         "FOR" I:=0 "STEP" 1 "UNTIL" 12 "DO" X[I]:=B[I]:=0;
         B[0]:=1;B[12]:=4;
         CONJ GRAD(A,X,B,0,12,IT<20 "AND" NO>"-10,IT,NO);
         OUTPUT(61,"(""("IT= ")",B3D5B,"("NO= ")",N,//,10B,
         "("X")",20B,"("R")",//")",IT,NO);
         "FOR" I:=0 "STEP" 1 "UNTIL" 12 "DO"
         OUTPUT(61,"("N,5B,N,/")",X[I],B[I])
     "END"
 
     DELIVERS:
     IT= 013     NO= +3.3424581859911"-027
 
                X                    R
 
     +1.2142857142857"+000     -7.1054273576010"-015
     +1.4285714285715"+000     +1.5151278924296"-014
     +1.6428571428572"+000     -1.3184703260130"-014
     +1.8571428571429"+000     +1.6718441615946"-014
     +2.0714285714286"+000     -1.5514524667596"-014
     +2.2857142857144"+000     +2.2130179956186"-014
     +2.5000000000001"+000     -2.2524167805437"-014
     +2.7142857142858"+000     +2.0834049529361"-014
     +2.9285714285715"+000     -1.8674557504802"-014
     +3.1428571428572"+000     +1.9163204503355"-014
     +3.3571428571429"+000     -1.2366043539824"-014
     +3.5714285714286"+000     +8.2548347242718"-015
     +3.7857142857143"+000     +4.4408920985006"-016  .
 
 
1SECTION 3.1.2.2.1            (JUNE 1974)                         PAGE 4
 
 
 
 SOURCE TEXT(S):
0"CODE" 34220;
     "PROCEDURE" CONJ GRAD( MATVEC, X, R, L, N, GO ON, ITERATE, NORM2);
     "VALUE" L, N; "PROCEDURE" MATVEC; "ARRAY" X, R; "BOOLEAN" GO ON;
     "INTEGER" L, N, ITERATE; "REAL" NORM2;
     "BEGIN" "ARRAY" P, AP[ L: N];
         "INTEGER" I;
         "REAL"  A, B, PRR, RRP;
         "REAL" "PROCEDURE" VECVEC( A, B, C, D, E); "CODE" 34010;
         "PROCEDURE" ELMVEC( A, B, C, D, E, F); "CODE" 34020;
         "FOR" ITERATE:= 0, ITERATE + 1 "WHILE" GO ON "DO"
         "BEGIN" "IF" ITERATE = 0 "THEN"
             "BEGIN" MATVEC( X, P);
                 "FOR" I:= L "STEP" 1 "UNTIL" N "DO"
                 P[ I]:= R[ I]:= R[ I] - P[ I];
                 PRR:= VECVEC( L, N, 0, R, R)
             "END" "ELSE"
             "BEGIN" B:= RRP / PRR; PRR:= RRP;
                 "FOR" I:= L "STEP" 1 "UNTIL" N "DO"
                 P[ I]:= R[ I] + B * P[ I]
             "END";
             MATVEC( P, AP);
             A:= PRR / VECVEC( L, N, 0, P, AP);
             ELMVEC( L, N, 0, X, P, A);
             ELMVEC( L, N, 0, R, AP, -A);
             NORM2:= RRP:= VECVEC( L, N, 0, R, R)
         "END"
     "END" CONJ GRAD;
         "EOP"
1SECTION:3.2.1.1.1            (JUNE 1974)                         PAGE 1      
        
        
        
 AUTHORS  :      T.J.DEKKER AND W.HOFFMANN.     
        
        
 CONTRIBUTORS: W.HOFFMANN, J.G.VERWER.
        
        
 INSTITUTE:    MATHEMATICAL CENTRE.   
        
        
 RECEIVED:     731022.      
        
        
 BRIEF DESCRIPTION:         
        
     THIS SECTION  CONTAINS TWO PROCEDURES.     
     A) EQILBR  EQUILIBRATES A MATRIX BY MEANS OF A DIAGONAL  SIMILARITY      
     TRANSFORMATION,        
     B) BAKLBR  PERFORMS  THE  CORRESPONDING  BACK TRANSFORMATION ON THE      
     COLUMNS OF A MATRIX AND SHOULD BE CALLED AFTER EQILBR.         
        
 KEYWORDS:        
        
     SIMILARITY TRANSFORMATION,       
     EQUILIBRATION.         
        
        
 SUBSECTION: EQILBR.        
        
        
 CALLING SEQUENCE:
     THE HEADING OF THE PROCEDURE IS: 
     "PROCEDURE" EQILBR(A, N, EM, D, INT); "VALUE" N;     
     "INTEGER" N; "ARRAY" A, EM, D; "INTEGER" "ARRAY" INT;
        
     THE MEANING OF THE FORMAL PARAMETERS IS:   
     N:      <ARITHMETIC EXPRESSION>; 
             THE ORDER OF THE GIVEN MATRIX;     
     A:      <ARRAY IDENTIFIER>;      
             "ARRAY"A[1:N,1:N];       
             ENTRY:  THE MATRIX TO BE EQUILIBRATED;       
             EXIT:   THE EQUILIBRATED MATRIX;   
     EM:     <ARRAY IDENTIFIER>;      
             "ARRAY"EM[0:0];
             ENTRY:  EM[0], THE MACHINE PRECISION;        
     D:      <ARRAY IDENTIFIER>;      
             "ARRAY"D[1:N]; 
             EXIT:   THE  MAIN  DIAGONAL  OF  THE TRANSFORMING  DIAGONAL      
                     MATRIX;
     INT:    <ARRAY IDENTIFIER>;      
             "INTEGER""ARRAY" INT[1:N];         
             EXIT:   INFORMATION DEFINING THE POSSIBLE INTERCHANGING  OF      
                     SOME ROWS AND THE CORRESPONDING COLUMNS;       
        
        
1SECTION:3.2.1.1.1            (JUNE 1974)                         PAGE 2      
        
        
        
 PROCEDURES USED: 
     TAMMAT          =      CP34014,  
     MATTAM          =      CP34015,  
     ICHCOL          =      CP34031,  
     ICHROW          =      CP34032.  
        
        
 RUNNING TIME: ROUGHLY PROPORTIONAL TO N SQUARED.         
        
 LANGUAGE:   ALGOL 60.      
        
        
        
 METHOD AND PERFORMANCE:    
     THE    MATRIX  IS EQUILIBRATED  BY  MEANS  OF OSBORNE'S  DIAGONAL        
     SIMILARITY    TRANSFORMATION    POSSIBLY   WITH   INTERCHANGES [2].      
     THE TRANSFORMING DIAGONAL MATRIX AND THE  EQUILIBRATED  MATRIX  ARE      
     CALCULATED ITERATIVELY:
     IN EACH STEP A CERTAIN COLUMN OF  THE  MATRIX IS MULTIPLIED BY, AND      
     THE CORRESPONDING ROW DIVIDED BY, A FACTOR WHICH IS CHOSEN IN  SUCH      
     A WAY THAT THE  CONSIDERED  COLUMN AND ROW OBTAIN  ROUGHLY THE SAME      
     EUCLIDEAN  NORM (IN  FACT, THE  FACTOR  IS  ROUNDED  TO THE NEAREST      
     INTEGRAL  POWER  OF  2, IN  ORDER TO PREVENT ROUNDING  ERRORS); THE      
     COLUMNS AND ROWS ARE HANDLED IN CYCLIC ORDER. IF  THE  MATRIX  DOES      
     NOT CONTAIN COLUMNS OR ROWS  WHOSE  OFF-DIAGONAL  ELEMENTS ARE 0 OR      
     NEARLY 0, THEN THE PROCESS (WITH UNROUNDED FACTORS) CONVERGES,  AND      
     IN  PRACTICE  A  FEW  STEPS  ARE  NEEDED  TO  OBTAIN  A  REASONABLY      
     EQUILIBRATED MATRIX [2].         
     IF ALL OFF-DIAGONAL  ELEMENTS OF SOME CONSIDERED COLUMN (ROW) ARE 0      
     OR NEARLY 0, THEN THIS  COLUMN (ROW) IS INTERCHANGED WITH THE FIRST      
     NONZERO COLUMN (LAST NONZERO ROW)  OF THE  MATRIX, AND, IN ORDER TO      
     HAVE A SIMILARITY TRANSFORMATION, THE CORRESPONDING  ROWS (COLUMNS)      
     ARE ALSO INTERCHANGED;  THEN  FOR THE  FURTHER  EQUILIBRATION,  THE      
     SUBMATRIX IS CONSIDERED WHICH DOES NOT  CONTAIN  SUCH  ZERO COLUMNS      
     AND ROWS AND THE CORRESPONDING ROWS AND COLUMNS. THE  EQUILIBRATION      
     PROCESS IS CONTINUED UNTIL, IN A WHOLE CYCLE NO FACTOR > 2 OR < 0.5      
     AND NO ZERO COLUMN OR ROW IS FOUND, OR UNTIL  (N + 1) * N ** 2 ROWS      
     AND COLUMNS HAVE BEEN CONSIDERED.
        
        
1SECTION:3.2.1.1.1            (JUNE 1974)                         PAGE 3      
        
        
        
 SUBSECTION: BAKLBR.        
        
 CALLING SEQUENCE:
     THE HEADING OF THE PROCEDURE IS: 
     "PROCEDURE" BAKLBR(N, N1, N2, D, INT, VEC); "VALUE" N, N1, N2; 
     "INTEGER" N, N1, N2; "ARRAY" D, VEC; "INTEGER" "ARRAY" INT;    
        
     THE MEANING OF THE FORMAL PARAMETERS IS:   
     N:      <ARITHMETIC EXPRESSION>; 
             THE LENGTH OF THE VECTORS TO BE TRANSFORMED; 
     N1, N2: <ARITHMETIC EXPRESSION>; 
             THE  SERIAL   NUMBERS  OF  THE  FIRST AND LAST VECTOR TO BE      
             TRANSFORMED;   
     VEC:    <ARRAY IDENTIFIER>;      
             "ARRAY"VEC[1:N,N1:N2];   
             ENTRY:  THE   N2 - N1 + 1   VECTORS   OF   LENGTH N  TO  BE      
                     TRANSFORMED;     
             EXIT:   THE N2 - N1 + 1 VECTORS OF LENGTH N RESULTING  FROM      
                     THE BACK TRANSFORMATION;   
     D:      <ARRAY IDENTIFIER>;      
             "ARRAY"D[1:N]; 
             ENTRY:  THE MAIN DIAGONAL  OF  THE  TRANSFORMING   DIAGONAL      
                     MATRIX OF ORDER N, AS PRODUCED BY EQILBR;      
     INT:    <ARRAY IDENTIFIER>;      
             "INTEGER""ARRAY" INT[1:N];         
             ENTRY:  INFORMATION DEFINING THE POSSIBLE INTERCHANGING  OF      
                     SOME ROWS AND COLUMNS, AS PRODUCED BY EQILBR.  
        
 PROCEDURES USED: 
     ICHROW          =      CP34032.  
        
 RUNNING TIME: ROUGHLY PROPORTIONAL TO (N2 - N1 + 1) * N. 
        
 LANGUAGE:   ALGOL 60.      
        
 METHOD AND PERFORMANCE:    
     THE  BACK  TRANSFORMATION, WHICH  CORRESPONDS  WITH  THE   DIAGONAL      
     SIMILARITY   TRANSFORMATION  AS  PERFORMED  BY  EQILBR,  TRANSFORMS      
     A  VECTOR X INTO  A  VECTOR  DX  AND  PERFORMS  THE   CORRESPONDING      
     INTERCHANGES. THE MATRIX D IS THE DIAGONAL MATRIX OF  THE  DIAGONAL      
     SIMILARITY TRANSFORMATION.       
        
 REFERENCES:      
     [1]      DEKKER, T. J. AND HOFFMANN, W.    
              ALGOL 60 PROCEDURES IN NUMERICAL ALGEBRA, PART 2,     
              MATHEMATICAL CENTRE TRACTS 23,    
              MATHEMATISCH CENTRUM, AMSTERDAM, 1968;      
     [2]      OSBORNE, E. E, ON PRECONDITIONING OF MATRICES,        
              J. ACM 7(1960) 338-354. 
        
 EXAMPLES OF USE: 
     EXAMPLES OF USE OF EQILBR AND BAKLBR CAN BE FOUND IN THE PROCEDURES      
     FOR  CALCULATING  EIGENVALUES  AND  EIGENVECTORS  AS  DESCRIBED  IN      
     SECTION 3.3.1.2.2.     
1SECTION:3.2.1.1.1            (JUNE 1974)                         PAGE 4      
        
        
        
 SOURCE TEXT(S) : 
0"CODE" 34173;    
     "COMMENT" MCA 2405;    
     "PROCEDURE" EQILBR(A, N, EM, D, INT); "VALUE" N; "INTEGER" N;  
      "ARRAY" A, EM, D; "INTEGER" "ARRAY" INT;  
     "BEGIN" "INTEGER" I, IM, I1, P, Q, J, T, COUNT, EXPONENT, NI;  
         "REAL" C, R, EPS, OMEGA, FACTOR;       
        
         "PROCEDURE" MOVE(K); "VALUE" K; "INTEGER" K;     
         "BEGIN" "REAL" DI; 
             NI:= Q - P; T:= T + 1; "IF" K ^= I "THEN"    
             "BEGIN" ICHCOL(1, N, K, I, A); ICHROW(1, N, K, I, A);  
                 DI:= D[I]; D[I]:= D[K]; D[K]:= DI        
             "END"
         "END" MOVE;        
        
         "REAL" "PROCEDURE" TAMMAT(L, U, I, J, A, B); "CODE" 34014; 
         "REAL" "PROCEDURE" MATTAM(L, U, I, J, A, B); "CODE" 34015; 
         "PROCEDURE" ICHCOL(L, U, I, J, A); "CODE" 34031; 
         "PROCEDURE" ICHROW(L, U, I, J, A); "CODE" 34032; 
        
         FACTOR:= 1 / (2 * LN(2)); "COMMENT" MORE GENERALLY: LN(BASE);        
         EPS:= EM[0]; OMEGA:= 1 / EPS; T:= P:= 1; Q:= NI:= I:= N;   
         COUNT:= (N + 1) * N // 2;    
         "FOR" J:= 1 "STEP" 1 "UNTIL" N "DO"    
         "BEGIN" D[J]:= 1; INT[J]:= 0 "END";    
         "FOR" I:= "IF" I < Q "THEN" I + 1 "ELSE" P       
         "WHILE" COUNT > 0 "AND" NI > 0 "DO"    
         "BEGIN" COUNT:= COUNT - 1; IM:= I - 1; I1:= I + 1;         
             C:= SQRT(TAMMAT(P, IM, I, I, A, A) +         
             TAMMAT(I1, Q, I, I, A, A));        
             R:= SQRT(MATTAM(P, IM, I, I, A, A) +         
             MATTAM(I1, Q, I, I, A, A));        
             "IF" C * OMEGA <= R * EPS "THEN"   
             "BEGIN" INT[T]:= I; MOVE(P); P:= P + 1 "END" 
             "ELSE" "IF" R * OMEGA <= C * EPS "THEN"      
             "BEGIN" INT[T]:= -I; MOVE(Q); Q:= Q - 1 "END"
             "ELSE"         
             "BEGIN" EXPONENT:= LN(R / C) * FACTOR;       
                 "IF" ABS(EXPONENT) > 1 "THEN"  
                 "BEGIN" NI:= Q - P; C:= 2 ** EXPONENT; R:= 1 / C;  
                     D[I]:= D[I] * C; 
                     "FOR" J:= 1 "STEP" 1 "UNTIL" IM,     
                     I1 "STEP" 1 "UNTIL" N "DO" 
                     "BEGIN" A[J,I]:= A[J,I] * C;         
                         A[I,J]:= A[I,J] * R    
                     "END"  
                 "END" "ELSE" NI:= NI - 1       
             "END"
         "END"    
     "END" EQILBR;
         "EOP"    
1SECTION:3.2.1.1.1            (JUNE 1974)                         PAGE 5      
        
        
0"CODE" 34174;    
     "COMMENT" MCA 2406;    
     "PROCEDURE" BAKLBR(N, N1, N2, D, INT, VEC); "VALUE" N, N1, N2; 
     "INTEGER" N, N1, N2; "ARRAY" D, VEC; "INTEGER" "ARRAY" INT;    
     "BEGIN" "INTEGER" I, J, K, P, Q; "REAL" DI;
        
         "PROCEDURE" ICHROW(L, U, I, J, A); "CODE" 34032; 
        
         P:= 1; Q:= N;      
         "FOR" I:= 1 "STEP" 1 "UNTIL" N "DO"    
         "BEGIN" DI:= D[I]; "IF" DI ^= 1 "THEN" 
             "FOR" J:= N1 "STEP" 1 "UNTIL" N2 "DO" VEC[I,J]:=       
             VEC[I,J] * DI; K:= INT[I];         
             "IF" K > 0 "THEN" P:= P + 1 "ELSE" 
             "IF" K < 0 "THEN" Q:= Q - 1        
         "END";   
         "FOR" I:= P - 1 + N - Q "STEP" -1 "UNTIL" 1 "DO" 
         "BEGIN" K:= INT[I]; "IF" K > 0 "THEN"  
             "BEGIN" P:= P - 1; "IF" K ^= P "THEN"        
                 ICHROW(N1, N2, K, P, VEC)      
             "END" "ELSE"   
             "BEGIN" Q:= Q + 1; "IF" -K ^= Q "THEN"       
                 ICHROW(N1, N2, -K, Q, VEC)     
             "END"
         "END"    
     "END" BAKLBR;
         "EOP"    
1SECTION:3.2.1.1.2            (JUNE 1974)                         PAGE 1      
        
        
        
 AUTHOR   : C.G. VAN DER LAAN.        
        
        
 CONTRIBUTORS : H.FIOLET, C.G. VAN DER LAAN.    
        
        
 INSTITUTE : MATHEMATICAL CENTRE.     
        
        
 RECEIVED: 731008.
        
        
 BRIEF DESCRIPTION :        
        
     THIS SECTION CONTAINS THE PROCEDURES EQILBRCOM AND BAKLBRCOM.  
     EQILBRCOM EQUILIBRATES A GIVEN MATRIX.     
     BAKLBRCOM TRANSFORMS THE EIGENVECTORS OF THE EQUILIBRATED MATRIX         
     INTO THE EIGENVECTORS OF THE ORIGINAL MATRIX.        
        
        
 KEYWORDS :       
        
     COMPLEX MATRIX,        
     EIGENVECTORS,
     EQUILIBRATION.         
        
        
 SUBSECTION: EQILBRCOM.     
        
        
 CALLING SEQUENCE:
        
     THE HEADING OF THE PROCEDURE READS:        
     "PROCEDURE" EQILBRCOM(A1, A2, N, EM, D, INT); "VALUE" N;       
     "INTEGER" N; "ARRAY" A1, A2, EM, D; "INTEGER" "ARRAY" INT;     
        
     THE MEANING OF THE FORMAL PARAMETERS IS:   
     A1,A2::    <ARRAY IDENTIFIER>;   
                "ARRAY" A1,A2[1:N,1:N];         
                ENTRY:      
                THE REAL PART AND IMAGINARY PART OF THE MATRIX TO  BE         
                EQUILIBRATED  MUST BE GIVEN IN THE ARRAYS A1 AND A2,
                RESPECTIVELY;         
                EXIT:       
                THE REAL PART AND THE IMAGINARY PART OF THE EQUILIBRATED      
                MATRIX ARE DELIVERED IN THE ARRAYS A1 AND A2,       
                RESPECTIVELY;         
     N:         <ARITHMETIC EXPRESSION>;        
                THE ORDER OF THE GIVEN MATRIX;  
     EM:        <ARRAY IDENTIFIER>;   
                "ARRAY"EM[0:7];       
                ENTRY:      
                EM[0]: THE MACHINE PRECISION;   
                EM[6]: THE MAXIMUM ALLOWED NUMBER OF ITERATIONS;    
                EXIT:       
                EM[7]: THE NUMBER OF ITERATIONS PERFORMED;
1SECTION:3.2.1.1.2            (JUNE 1974)                         PAGE 2      
        
        
        
     D:         <ARRAY IDENTIFIER>:   
                "ARRAY" D[1:N];       
                EXIT:       
                THE SCALING FACTORS OF THE DIAGONAL SIMILARITY      
                TRANSFORMATION;       
     INT:       <ARRAY IDENTIFIER>;   
                "INTEGER""ARRAY" INT[1:N];      
                EXIT:       
                INFORMATION CONCERNING THE POSSIBLE INTERCHANGING OF
                SOME ROWS AND CORRESPONDING COLUMNS.      
        
        
 PROCEDURES USED: 
        
     ICHCOL = CP34031,      
     ICHROW = CP34332,      
     TAMMAT = CP34014,      
     MATTAM = CP34015.      
        
        
 RUNNING TIME: PROPORTIONAL TO N * NUMBER OF ITERATIONS.  
        
        
 LANGUAGE: ALGOL 60.        
        
        
 METHOD AND PERFORMANCE: SEE BAKLBRCOM.         
        
        
 SUBSECTION: BAKLBRCOM.     
        
        
 CALLING SEQUENCE:
        
     THE HEADING OF THE PROCEDURE READS:        
     "PROCEDURE" BAKLBRCOM(N, N1, N2, D, INT, VR, VI);    
     "VALUE" N, N1, N2; "INTEGER" N, N1, N2; "ARRAY" D, VR, VI;     
     "INTEGER" "ARRAY" INT; 
        
     THE MEANING OF THE FORMAL PARAMETERS IS:   
     N:         <ARITHMETIC EXPRESSION>;        
                THE ORDER OF THE MATRIX OF WHICH THE EIGENVECTORS ARE         
                CALCULATED; 
     N1,N2:     <ARITHMETIC EXPRESSION>;        
                THE EIGENVECTORS CORRESPONDING TO THE EIGENVALUES WITH        
                INDICES N1,...,N2 ARE TO BE TRANSFORMED;  
     D:         <ARRAY IDENTIFIER>:   
                "ARRAY" D[1:N];       
                ENTRY: THE SCALING FACTORS OF THE DIAGONAL SIMILARITY         
                TRANSFORMATION AS DELIVERED BY EQILBRCOM; 
     INT:       <ARRAY IDENTIFIER>;   
                "INTEGER""ARRAY" INT[1:N];      
                ENTRY: INFORMATION DEFINING THE INTERCHANGING OF SOME         
                ROWS AND COLUMNS, AS DELIVERED BY EQILBRCOM;        
1SECTION:3.2.1.1.2            (DECEMBER 1975)                     PAGE 3      
        
        
        
     VR,VI:     <ARRAY IDENTIFIER>;   
                "ARRAY" VR,VI[1:N,N1:N2];       
                ENTRY:      
                THE BACK TRANSFORMATION IS PERFORMED ON THE EIGENVECTORS      
                WITH THE REAL PARTS GIVEN IN ARRAY VR AND THE IMAGINARY       
                PARTS GIVEN IN ARRAY VI;        
                EXIT:       
                THE REAL PARTS AND IMAGINARY PARTS OF THE RESULTING 
                EIGENVECTORS ARE DELIVERED IN THE COLUMNS OF THE ARRAYS       
                VR AND VI, RESPECTIVELY.        
        
        
 PROCEDURES USED: BAKLBR = CP34174.   
        
        
 RUNNING TIME: ROUGHLY PROPORTIONAL TO N * (N2-N1).       
        
        
 LANGUAGE: ALGOL 60.        
        
        
 THE FOLLOWING HOLDS FOR BOTH PROCEDURES:       
        
        
 METHOD AND PERFORMANCE:    
        
     A MATRIX M IS SAID TO BE EQUILIBRATED, WHEN THE DIAGONAL ELEMENTS        
     OF M''M - MM'' ARE ZERO, WHERE '' STANDS FOR CONJUGATING       
     AND TRANSPOSING.  IN  EQILBRCOM  THE  MATRIX  M  IS  EQUILIBRATED        
     BY MEANS OF OSBORNE'S DIAGONAL SIMILARITY TRANSFORMATION WITH  
     POSSIBLE INTERCHANGES (OSBORNE, 1960).     
     BAKLBRCOM PERFORMS THE CORRESPONDING BACK TRANSFORMATION.      
     LET THE EIGENVECTORS OF THE EQUILIBRATED MATRIX BE GIVEN IN THE
     COLUMNS OF MATRIX V. THE EIGENVECTORS OF THE ORIGINAL MATRIX ARE         
     OBTAINED BY MULTIPLYING (OR POSSIBLE INTERCHANGING) THE ROWS OF THE      
     MATRIX V WITH THE SCALING FACTORS. AS THE SCALING FACTORS ARE REAL       
     QUANTITIES, THE TRANSFORMATION IS PERFORMED BY CALLING THE     
     PROCEDURE BAKLBR FOR BOTH VR AND VI (DEKKER AND HOFFMANN, 1968).         
        
        
 REFERENCES:      
        
     DEKKER, T.J. AND W.HOFFMANN (1968),        
     ALGOL 60 PROCEDURES IN NUMERICAL ALGEBRA, PART 2,    
     MATH. CENTRE TRACTS 23, MATHEMATISCH CENTRUM;        
        
     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;       
1SECTION:3.2.1.1.2            (JUNE 1974)                         PAGE 4      
        
        
        
 EXAMPLE OF USE:  
        
     BAKLBRCOM IS USED IN THE PROCEDURE EIGCOM (SEE SECTION 3.3.2.2.2.).      
     AS A FORMAL TEST OF THE PROCEDURE EQILBRCOM, THE FOLLOWING MATRIX        
     WAS USED:    
       1     0  1024*I      
       0     1    0         
     I/1024  0    2         
        
     "BEGIN" "INTEGER" I,J; 
     "REAL" "ARRAY" A1,A2[1:3,1:3],EM[0:7],D,INT[1:3];    
     "PROCEDURE" INIMAT(LR,UR,LC,UC,A,X);"CODE" 31011;    
     "PROCEDURE" EQILBRCOM(A1,A2,N,EM,D,INT);"CODE" 34361;
     EM[0]:=5"-14;EM[6]:=10;
     INIMAT(1,3,1,3,A1,0);INIMAT(1,3,1,3,A2,0); 
     A1[1,1]:=A1[2,2]:=1;A1[3,3]:=2;  
     A2[1,3]:=2**10;A2[3,1]:=2**(-10);
     EQILBRCOM(A1,A2,3,EM,D,INT);     
     OUTPUT(61,"(""("EQUILIBRATED MATRIX:")",/")");       
     OUTPUT(61,"("3(D2B),/,2(D2B),"("I")",/,D2B,"("I")",2BD/")",    
              A1[1,1],A1[1,2],A1[1,3],A1[2,1],A1[2,2],A1[3,1],A1[3,3]);       
     OUTPUT(61,"("/,"("EM[7]:")",5BD/")",EM[7]);
     OUTPUT(61,"(""("D[1:3]: ")",3(3ZD2B),/")",D[1],D[2],D[3]);     
     OUTPUT(61,"(""("INT[1:3]: ")",BD,3B,2BD,4BD")",      
                  INT[1],INT[2],INT[3])         
     "END"        
        
     OUTPUT:      
     EQUILIBRATED MATRIX:   
     1  0  0      
     0  1  I      
     0  I  2      
        
     EM[7]:     4 
     D[1:3]:    1  1024     1         
     INT[1:3]:  2     0     0         
        
        
 SOURCE TEXT(S) : 
0"CODE" 34361;    
     "PROCEDURE" EQILBRCOM(A1, A2, N, EM, D, INT); "VALUE" N;       
     "INTEGER" N; "ARRAY" A1, A2, EM, D; "INTEGER" "ARRAY" INT;     
     "BEGIN" "INTEGER" I, P, Q, J, T, COUNT, EXPONENT, NI, IM, I1;  
         "REAL" C, R, EPS;  
         "PROCEDURE" ICHCOL(L,U,I,J,A);"CODE" 34031;      
         "PROCEDURE" ICHROW(L,U,I,J,A);"CODE" 34032;      
         "REAL" "PROCEDURE" TAMMAT(L,U,I,J,A,B);"CODE" 34014;       
         "REAL" "PROCEDURE" MATTAM(L,U,I,J,A,B);"CODE" 34015;       
                                                               "COMMENT"      
1SECTION:3.2.1.1.2            (JUNE 1974)                         PAGE 5      
                                                                  ; 
        
        
         "PROCEDURE" MOVE(K); "VALUE" K; "INTEGER" K;     
         "BEGIN" "REAL" DI; 
             NI:= Q - P; T:= T + 1; "IF" K ^= I "THEN"    
             "BEGIN" ICHCOL(1, N, K, I, A1); ICHROW(1, N, K, I, A1);
                 ICHCOL(1, N, K, I, A2); ICHROW(1, N, K, I, A2);    
                 DI:= D[I]; D[I]:= D[K]; D[K]:= DI        
             "END"
         "END" MOVE;        
         EPS:= EM[0] ** 4; T:= P:= 1; Q:= NI:= I:= N;     
         COUNT:= EM[6];     
         "FOR" J:= 1 "STEP" 1 "UNTIL" N "DO"    
         "BEGIN" D[J]:= 1; INT[J]:= 0 "END";    
         "FOR" I:= "IF" I < Q "THEN" I + 1 "ELSE" P "WHILE" COUNT > 0         
         "AND" NI > 0 "DO"  
         "BEGIN" COUNT:= COUNT - 1; IM:= I - 1; I1:= I + 1;         
             C:= TAMMAT(P, IM, I, I, A1, A1) + TAMMAT(I1, Q, I,     
             I, A1, A1) + TAMMAT(P, IM, I, I, A2, A2) +   
             TAMMAT(I1, Q, I, I, A2, A2);       
             R:= MATTAM(P, IM, I, I, A1, A1) + MATTAM(I1, Q, I,     
             I, A1, A1) + MATTAM(P, IM, I, I, A2, A2) +   
             MATTAM(I1, Q, I, I, A2, A2); "IF" C / EPS <= R "THEN"  
             "BEGIN" INT[T]:= I; MOVE(P); P:= P + 1 "END" 
             "ELSE" "IF" R / EPS <= C "THEN"    
             "BEGIN" INT[T]:= - I; MOVE(Q); Q:= Q - 1 "END"         
             "ELSE"         
             "BEGIN" EXPONENT:= LN(R / C) * 0.36067;      
                 "IF" ABS(EXPONENT) > 1 "THEN"  
                 "BEGIN" NI:= Q - P; C:= 2 ** EXPONENT;   
                     D[I]:= D[I] * C; 
                     "FOR" J:= 1 "STEP" 1 "UNTIL" IM, I1 "STEP" 1   
                     "UNTIL" N "DO"   
                     "BEGIN" A1[J,I]:= A1[J,I] * C;       
                         A1[I,J]:= A1[I,J] / C; 
                         A2[J,I]:= A2[J,I] * C; 
                         A2[I,J]:= A2[I,J] / C  
                     "END"  
                 "END"      
                 "ELSE" NI:= NI - 1   
             "END"
         "END";   
         EM[7]:= EM[6] - COUNT        
     "END" EQILBRCOM;       
         "EOP"    
0"CODE" 34362;    
     "PROCEDURE" BAKLBRCOM(N, N1, N2, D, INT, VR, VI);    
     "VALUE" N, N1, N2; "INTEGER" N, N1, N2; "ARRAY" D, VR, VI;     
     "INTEGER" "ARRAY" INT; 
     "BEGIN"      
         "PROCEDURE" BAKLBR(N,N1,N2,D,INT,VEC);"CODE" 34174;        
         BAKLBR(N, N1, N2, D, INT, VR);         
         BAKLBR(N, N1, N2, D, INT, VI)
     "END" BAKLBRCOM;       
         "EOP"