1SECTION: 2.3                 (MAY 1974)                          PAGE 1
 
 
 
 AUTHOR   : H.FIOLET
 
 
 INSTITUTE: MATHEMATICAL CENTRE.
 
 
 RECEIVED: 731105.
 
 
 BRIEF DESCRIPTION:
 
     JFRAC CALCULATES A TERMINATING CONTINUED FRACTION.
 
 
 KEYWORDS:
 
     CONTINUED FRACTION,
     TERMINATING CONTINUED FRACTION.
 
 
 CALLING SEQUENCE:
     THE HEADING OF THE PROCEDURE READS:
     "REAL" "PROCEDURE" JFRAC(N,A,B);
     "VALUE" N;"INTEGER" N;"ARRAY" A,B;
     "CODE" 35083;
 
     JFRAC DELIVERS THE VALUE OF THE TERMINATING CONTINUED FRACTION:
     B[0]+A[1]/(B[1]+A[2]/(B[2]+A[3]/(B[3]+ . . . + A[N]/B[N])))...))
 
     THE MEANING OF THE FORMAL PARAMETERS IS:
     N:         <ARITHMETIC EXPRESSION>;
                THE UPPER INDEX OF THE ARRAYS A AND B;
     A,B:       <ARRAY IDENTIFIER>;
                "ARRAY" A[1:N];
                "ARRAY" B[0:N];
                THE ELEMENTS OF THE CONTINUED FRACTION:
                B[0]+A[1]/(B[1]+A[2]/(B[2]+A[3]/(B[3]+ . . . +
                + A[N]/B[N])))...)).
 
 
 PROCEDURES USED: NONE.
 
 
 RUNNING TIME: PROPORTIONAL TO N.
 
 
 
 
1SECTION: 2.3                 (MAY 1974)                          PAGE 2
 
 
 
 EXAMPLE OF USE:
 
     "BEGIN"
     "REAL" "PROCEDURE" JFRAC(N,A,B);"CODE" 35083;
     "REAL" "ARRAY" P[1:10],Q[0:10];
     "INTEGER" I;
     "FOR" I:=1 "STEP" 1 "UNTIL" 10 "DO"
     "BEGIN" P[I]:=1;Q[I]:=2 "END";
     Q[0]:=1;
     "FOR" I:=7 "STEP" 1 "UNTIL" 10 "DO"
     OUTPUT(61,"("N/")",JFRAC(I,P,Q))
     "END"
 
     DELIVERS:
 
     +1.4142156862745"+000
     +1.4142131979695"+000
     +1.4142136248949"+000
     +1.4142135516461"+000     .
 
 
 SOURCE TEXT:
 
 "CODE" 35083;
     "REAL" "PROCEDURE" JFRAC(N,A,B);
     "VALUE" N;"INTEGER" N;"ARRAY" A,B;
     "BEGIN" "REAL" D;"INTEGER" I;
         D:=0;
         "FOR" I:=N "STEP" -1 "UNTIL" 1 "DO" D:=A[I]/(B[I]+D);
         JFRAC:=D+B[0]
     "END" JFRAC;
         "EOP"
1SECTION : 2.4.1             (DECEMBER 1978)                     PAGE 1
 
 
 
 AUTHOR:       C.G. VAN DER LAAN.
 
 
 CONTRIBUTORS: C.G. VAN DER LAAN, M. VOORINTHOLT.
 
 
 INSTITUTE:    REKENCENTRUM DER RIJKSUNIVERSITEIT GRONINGEN.
 
 
 RECEIVED:     780601.
 
 
 BRIEF DESCRIPTION:
 
      WE CONSIDER THE REPRESENTATIONS
                                           N
              POWER SUM                :  SUM A[K]*X**K,
                                          K=0
                                           N
              CHEBYSHEV SUM            :  SUM A[K]*T[K](X),
                                          K=0
                                           N
              SHIFTED CHEBYSHEV SUM    :  SUM A[K]*T'[K](X),
                                          K=0
                                           N           K-1
              NEWTON SUM               :  SUM (A[K] * PROD (X-X[J])).
                                          K=0          J=0
 
      THE SHIFTED CHEBYSHEV POLYNOMIAL T'[N] IS
      DEFINED BY  T'[N](X) = T[N](2*X-1).
 
 
      THIS SECTION CONTAINS THE TRANSFORMATIONS:
 
      PROCEDURE NAME :  TRANSFORMATION
      ---------------:------------------------------------------------
         POLCHS      :  POWER SUM INTO CHEBYSHEV SUM
         CHSPOL      :  CHEBYSHEV SUM INTO POWER SUM
         POLSHTCHS   :  POWER SUM INTO SHIFTED CHEBYSHEV SUM
         SHTCHSPOL   :  SHIFTED CHEBYSHEV SUM INTO POWER SUM
         GRNNEW      :  POWER SUM INTO NEWTON SUM
         NEWGRN      :  NEWTON SUM INTO POWER SUM
         LINTFMPOL   :  POWER SUM IN X INTO POWER SUM IN Y, X=P*Y+Q
 
 
 KEYWORDS:
 
      TRANSFORMATION OF POLYNOMIAL REPRESENTATION.
1SECTION : 2.4.1             (DECEMBER 1978)                     PAGE 2
 
 
 
 RUNNING TIME:  PROPORTIONAL TO THE SQUARED DEGREE OF THE POLYNOMIAL.
 
 METHOD AND PERFORMANCE:
 
      ALTHOUGH THE TRANSFORMATION OF REPRESENTATIONS OF POLYNOMIALS
      COULD HAVE BEEN OBTAINED BY FAST EVALUATION AND FAST INTERPOLATION
      WE IMPLEMENTED THE ALGORITHM OF HAMMING (1973, 474,475), BECAUSE
      OF ITS SIMPLE APPEARANCE. AN EXPLANATION OF THE HAMMING ALGORITHM
      IS GIVEN IN VAN DER LAAN (1977,224-229).
 
 
 REFERENCES:
 
      HAMMING, R.W. (1973):
         NUMERICAL METHODS FOR SCIENTISTS AND ENGINEERS.
         MCGRAW-HILL.
      LAAN, C.G. VAN DER (1977):
         APPROXIMATIE VAN FUNCTIES EN DATA.
         IN: RIELE, H.J.J. TE (ED.);
             COLLOQUIUM NUMERIEKE PROGRAMMATUUR, DEEL 2,
             MC SYLLABUS 29.2, MATHEMATISCH CENTRUM AMSTERDAM.
 
 
 SUBSECTION: POLCHS.
 
 
 CALLING SEQUENCE:
 
      THE DECLARATION OF THE PROCEDURE IN THE CALLING PROGRAM READS:
      "PROCEDURE" POLCHS(N,A);
      "VALUE" N; "INTEGER" N; "ARRAY" A;
      "CODE" 31051;
 
      THE MEANING OF THE FORMAL PARAMETERS IS:
      N:      <ARITHMETIC EXPRESSION>;
              ENTRY: THE DEGREE OF THE POLYNOMIAL;
      A:      <ARRAY IDENTIFIER>;
              "ARRAY" A[0:N];
              ENTRY: THE COEFFICIENTS OF THE POWER SUM;
              EXIT:  THE COEFFICIENTS OF THE CHEBYSHEV SUM;
 
 
 PROCEDURES USED: NONE.
 
 
 EXAMPLE OF USE: SEE NEXT SUBSECTION.
1SECTION : 2.4.1             (DECEMBER 1978)                     PAGE 3
 
 
 
 SUBSECTION: CHSPOL.
 
 
 CALLING SEQUENCE:
 
      THE DECLARATION OF THE PROCEDURE IN THE CALLING PROGRAM READS:
      "PROCEDURE" CHSPOL(N,A);
      "VALUE" N; "INTEGER" N; "ARRAY" A;
      "CODE" 31052;
 
      THE MEANING OF THE FORMAL PARAMETERS IS:
      N:      <ARITHMETIC EXPRESSION>;
              ENTRY: THE DEGREE OF THE POLYNOMIAL;
      A:      <ARRAY IDENTIFIER>;
              "ARRAY" A[0:N];
              ENTRY: THE COEFFICIENTS OF THE CHEBYSHEV SUM;
              EXIT:  THE COEFFICIENTS OF THE POWER SUM;
 
 
 PROCEDURES USED: NONE.
 
 
 EXAMPLE OF USE:
 
      AS AN EXAMPLE WE TRANSFORMED
              THE POWER SUM:  1 + 2*X + 3*X**2
      INTO
              ITS CHEBYSHEV SUM;
      AS A CHECK WE TRANSFORMED THE LATTER REPRESENTATION
      BACK INTO THE ORIGINAL POWER SUM.
 
      "BEGIN"
          "PROCEDURE" POLCHS(N,A);
          "VALUE" N; "INTEGER" N; "ARRAY" A;
          "CODE" 31051;
          "PROCEDURE" CHSPOL(N,A);
          "VALUE" N; "INTEGER" N; "ARRAY" A;
          "CODE" 31052;
          "ARRAY" A[0:2];
          A[0]:=1; A[1]:=2; A[2]:=3;
          OUTPUT(61,"("/,16B,"("A[0]")",4B,"("A[1]")",4B,"("A[2]")"")");
          OUTPUT(61,"("/,"(" INPUT")",5B,"(":")",3(2B,+ZD.DD)")",A);
          POLCHS(2,A);
          OUTPUT(61,"("/,"(" POLCHS")",4B,"(":")",3(2B,+ZD.DD)")",A);
          CHSPOL(2,A);
          OUTPUT(61,"("/,"(" CHSPOL")",4B,"(":")",3(2B,+ZD.DD)")",A);
      "END" PROGRAM;
1SECTION : 2.4.1             (DECEMBER 1978)                     PAGE 4
 
 
 
                     A[0]    A[1]    A[2]
      INPUT     :   +1.00   +2.00   +3.00
      POLCHS    :   +2.50   +2.00   +1.50
      CHSPOL    :   +1.00   +2.00   +3.00
 
 
 
 SUBSECTION: POLSHTCHS.
 
 
 CALLING SEQUENCE:
 
      THE DECLARATION OF THE PROCEDURE IN THE CALLING PROGRAM READS:
      "PROCEDURE" POLSHTCHS(N,A);
      "VALUE" N; "INTEGER" N; "ARRAY" A;
      "CODE" 31053;
 
      THE MEANING OF THE FORMAL PARAMETERS IS:
      N:      <ARITHMETIC EXPRESSION>;
              ENTRY: THE DEGREE OF THE POLYNOMIAL;
      A:      <ARRAY IDENTIFIER>;
              "ARRAY" A[0:N];
              ENTRY: THE COEFFICIENTS OF THE POWER SUM;
              EXIT:  THE COEFFICIENTS OF THE SHIFTED CHEBYSHEV SUM;
 
 
 PROCEDURES USED: LINTFMPOL = CP31250,
                  POLCHS    = CP31051.
 
 
 EXAMPLE OF USE: SEE NEXT SUBSECTION.
 
 
 SUBSECTION: SHTCHSPOL.
 
 
 CALLING SEQUENCE:
 
      THE DECLARATION OF THE PROCEDURE IN THE CALLING PROGRAM READS:
      "PROCEDURE" SHTCHSPOL(N,A);
      "VALUE" N; "INTEGER" N; "ARRAY" A;
      "CODE" 31054;
 
      THE MEANING OF THE FORMAL PARAMETERS IS:
      N:      <ARITHMETIC EXPRESSION>;
              ENTRY: THE DEGREE OF THE POLYNOMIAL;
      A:      <ARRAY IDENTIFIER>;
              "ARRAY" A[0:N];
              ENTRY: THE COEFFICIENTS OF THE SHIFTED CHEBYSHEV SUM;
              EXIT:  THE COEFFICIENTS OF THE POWER SUM.
1SECTION : 2.4.1             (DECEMBER 1978)                     PAGE 5
 
 
 
 PROCEDURES USED: LINTFMPOL = CP31250,
                  CHSPOL    = CP31052.
 
 EXAMPLE OF USE:
 
      AS AN EXAMPLE WE TRANSFORMED
              THE POWER SUM:  1 + 2*X + 3*X**2
      INTO
              ITS SHIFTED CHEBYSHEV SUM;
      AS A CHECK WE TRANSFORMED THE LATTER REPRESENTATION
      BACK INTO THE ORIGINAL POWER SUM.
 
      "BEGIN"
          "PROCEDURE" POLSHTCHS(N,A);
          "VALUE" N; "INTEGER" N; "ARRAY" A;
          "CODE" 31053;
          "PROCEDURE" SHTCHSPOL(N,A);
          "VALUE" N; "INTEGER" N; "ARRAY" A;
          "CODE" 31054;
          "ARRAY" A[0:2];
          A[0]:=1; A[1]:=2; A[2]:=3;
          OUTPUT(61,"("/,16B,"("A[0]")",4B,"("A[1]")",4B,"("A[2]")"")");
          OUTPUT(61,"("/,"(" INPUT")",5B,"(":")",3(2B,+ZD.DD)")",A);
          POLSHTCHS(2,A);
          OUTPUT(61,"("/,"(" POLSHTCHS")",B,"(":")",3(2B,+ZD.DD)")",A);
          SHTCHSPOL(2,A);
          OUTPUT(61,"("/,"(" SHTCHSPOL")",B,"(":")",3(2B,+ZD.DD)")",A);
      "END" PROGRAM;
 
                     A[0]    A[1]    A[2]
      INPUT     :   +1.00   +2.00   +3.00
      POLSHTCHS :   +3.13   +2.50   +0.38
      SHTCHSPOL :   +1.00   +2.00   +3.00
 
 
 SUBSECTION: GRNNEW.
 
 
 CALLING SEQUENCE:
 
      THE DECLARATION OF THE PROCEDURE IN THE CALLING PROGRAM READS:
      "PROCEDURE" GRNNEW(N,X,A);
      "VALUE" N; "INTEGER" N; "ARRAY" X,A;
      "CODE" 31055;
1SECTION : 2.4.1             (DECEMBER 1978)                     PAGE 6
 
 
 
      THE MEANING OF THE FORMAL PARAMETERS IS:
      N:      <ARITHMETIC EXPRESSION>;
              ENTRY: THE DEGREE OF THE POLYNOMIAL;
      X:      <ARRAY IDENTIFIER>;
              "ARRAY" X[0:N-1];
              ENTRY: THE INTERPOLATION POINTS;
      A:      <ARRAY IDENTIFIER>;
              ENTRY: THE COEFFICIENTS OF THE POWER SUM;
              EXIT:  THE COEFFICIENTS OF THE NEWTON SUM;
 
 
 PROCEDURES USED: NONE.
 
 
 EXAMPLE OF USE: SEE NEXT SUBSECTION.
 
 
 
 SUBSECTION: NEWGRN.
 
 
 CALLING SEQUENCE:
 
      THE DECLARATION OF THE PROCEDURE IN THE CALLING PROGRAM READS:
      "PROCEDURE" NEWGRN(N,X,A);
      "VALUE" N; "INTEGER" N; "ARRAY" X,A;
      "CODE" 31050;
 
      THE MEANING OF THE FORMAL PARAMETERS IS:
      N:      <ARITHMETIC EXPRESSION>;
              ENTRY: THE DEGREE OF THE POLYNOMIAL;
      X:      <ARRAY IDENTIFIER>;
              "ARRAY" X[0:N-1];
              ENTRY: THE INTERPOLATION POINTS;
      A:      <ARRAY IDENTIFIER>;
              "ARRAY" A[0:N];
              ENTRY: THE COEFFICIENTS OF THE NEWTON SUM;
              EXIT:  THE COEFFICIENTS OF THE POWER SUM;
 
 
 PROCEDURES USED : ELMVEC    = CP34020.
 
 
 EXAMPLE OF USE:
 
      AS AN EXAMPLE WE TRANSFORMED
              THE POWER SUM:  1 + 2*X + 3*X**2
      INTO
              ITS NEWTON SUM WITH
              INTERPOLATION POINTS: X[0]:=1.0, X[1]:=2.0;
      AS A CHECK WE TRANSFORMED THE LATTER REPRESENTATION
      BACK INTO THE ORIGINAL POWER SUM.
1SECTION : 2.4.1             (DECEMBER 1978)                     PAGE 7
 
 
 
      "BEGIN"
          "PROCEDURE" GRNNEW(N,X,A);
          "VALUE" N; "INTEGER" N; "ARRAY" X,A;
          "CODE" 31055;
          "PROCEDURE" NEWGRN(N,X,A);
          "VALUE" N; "INTEGER" N; "ARRAY" X,A;
          "CODE" 31050;
          "ARRAY" X[0:1], A[0:2];
          A[0]:=1; A[1]:=2; A[2]:=3;
          X[0]:=1; X[1]:=2;
          OUTPUT(61,"("/,16B,"("A[0]")",4B,"("A[1]")",4B,"("A[2]")"")");
          OUTPUT(61,"("/,"(" INPUT")",5B,"(":")",3(2B,+ZD.DD)")",A);
          GRNNEW(2,X,A);
          OUTPUT(61,"("/,"(" GRNNEW")",4B,"(":")",3(2B,+ZD.DD)")",A);
          NEWGRN(2,X,A);
          OUTPUT(61,"("/,"(" NEWGRN")",4B,"(":")",3(2B,+ZD.DD)")",A);
      "END" PROGRAM;
 
                     A[0]    A[1]    A[2]
      INPUT     :   +1.00   +2.00   +3.00
      GRNNEW    :   +6.00  +11.00   +3.00
      NEWGRN    :   +1.00   +2.00   +3.00
 
 
 
 SUBSECTION: LINTFMPOL.
 
 
 CALLING SEQUENCE:
 
      THE DECLARATION OF THE PROCEDURE IN THE CALLING PROGRAM READS:
      "PROCEDURE" LINTFMPOL(P,Q,N,A); "VALUE" N,P,Q;
      "INTEGER" N; "REAL" P,Q; "ARRAY" A;
      "CODE" 31250;
 
      THE MEANING OF THE FORMAL PARAMETERS IS:
      N:      <ARITHMETIC EXPRESSION>;
              ENTRY: THE DEGREE OF THE POLYNOMIAL;
      P,Q:    <ARITHMETIC EXPRESSION>;
              ENTRY: DEFINING THE LINEAR TRANSFORMATION OF
                     THE INDEPENDENT VARIABLE X=P*Y+Q;
                     (P=0 GIVES THE VALUE OF THE POLYNOMIAL
                      WITH ARGUMENT Q.)
      A:      <ARRAY IDENTIFIER>;
              "ARRAY" A[0:N];
              ENTRY: THE COEFFICIENTS OF THE POWER SUM IN X;
              EXIT:  THE COEFFICIENTS OF THE POWER SUM IN Y;
 
 
 PROCEDURES USED: NORDERPOL  = CP31242.
1SECTION : 2.4.1             (DECEMBER 1978)                     PAGE 8
 
 
 
 EXAMPLE OF USE:
 
      AS AN EXAMPLE WE TRANSFORMED
              THE POWER SUM:  1 + 2*X + 3*X**2
      INTO
              ITS POWER SUM IN Y WITH
              X = 2*Y + 3;
      AS A CHECK WE TRANSFORMED THE LATTER REPRESENTATION
      BACK INTO THE ORIGINAL POWER SUM.
 
      "BEGIN"
          "PROCEDURE" LINTFMPOL(P,Q,N,A);
          "VALUE" N,P,Q; "INTEGER" N;
          "REAL" P,Q; "ARRAY" A;
          "CODE" 31250;
          "ARRAY" A[0:2];
          A[0]:=1; A[1]:=2; A[2]:=3;
          OUTPUT(61,"("/,16B,"("A[0]")",4B,"("A[1]")",4B,"("A[2]")"")");
          OUTPUT(61,"("/,"(" INPUT")",5B,"(":")",3(2B,+ZD.DD)")",A);
          LINTFMPOL(2,3,2,A);
          OUTPUT(61,"("/,"(" LINTFMPOL")",B,"(":")",3(2B,+ZD.DD)")",A);
          LINTFMPOL(1/2,-3/2,2,A);
          OUTPUT(61,"("/,"(" LINTFMPOL")",B,"(":")",3(2B,+ZD.DD)")",A);
      "END" PROGRAM;
 
                     A[0]    A[1]    A[2]
      INPUT     :   +1.00   +2.00   +3.00
      LINTFMPOL :  +34.00  +40.00  +12.00   (POWER SUM IN Y)
      LINTFMPOL :   +1.00   +2.00   +3.00   (POWER SUM IN X)
1SECTION : 2.4.1             (DECEMBER 1978)                     PAGE 9
 
 
 
 SOURCE TEXTS:
 
 "CODE" 31051;
 "PROCEDURE" POLCHS(N,A);
 "VALUE" N; "INTEGER" N; "ARRAY" A;
 "IF" N>1 "THEN"
 "BEGIN"
    "COMMENT"  SCALING;
    "INTEGER" K,L,TWOPOW;
    TWOPOW:=2;
    "FOR" K:=1 "STEP" 1 "UNTIL" N-2 "DO"
    "BEGIN"
       A[K]:=A[K]/TWOPOW; TWOPOW:=TWOPOW*2;
    "END";
    A[N-1]:=2*A[N-1]/TWOPOW;
    A[N]:=A[N]/TWOPOW;
    A[N-2]:=A[N-2]+A[N];
    "COMMENT" N<=2 READY;
    "FOR" K:=N-2 "STEP" -1 "UNTIL" 1 "DO"
    "BEGIN"
       A[K-1]:=A[K-1]+A[K+1]; A[K]:=A[K]*2 + A[K+2];
       "FOR" L:=K+1 "STEP" 1 "UNTIL" N-2 "DO"
       A[L]:=A[L]+A[L+2];
    "END";
 "END" POLCHS;
         "EOP"
 
 
 
 "CODE" 31052;
 "PROCEDURE" CHSPOL(N,A);
 "VALUE" N; "INTEGER" N; "ARRAY" A;
 "IF" N>1 "THEN"
 "BEGIN"
     "INTEGER" K,L,TWOPOW;
     "FOR" K:=0 "STEP" 1 "UNTIL" N-2 "DO"
     "BEGIN"
         "FOR" L:=N-2 "STEP" -1 "UNTIL" K "DO"
         A[L]:=A[L]-A[L+2]; A[K+1]:=A[K+1]/2;
     "END";
     TWOPOW:=2;
     "FOR" K:=1 "STEP" 1 "UNTIL" N-2 "DO"
     "BEGIN"
         A[K]:=A[K]*TWOPOW; TWOPOW:=TWOPOW*2;
     "END";
     A[N-1]:=TWOPOW*A[N-1];
     A[N]:=TWOPOW*A[N];
 "END" CHSPOL;
         "EOP"
1SECTION : 2.4.1             (DECEMBER 1978)                     PAGE 10
 
 
 
 "CODE" 31053;
 "PROCEDURE" POLSHTCHS(N,A);
 "VALUE" N; "INTEGER" N; "ARRAY" A;
 "BEGIN"
     "PROCEDURE" LINTFMPOL(P,Q,N,A);
     "VALUE" N,P,Q; "INTEGER" N; "REAL" P,Q; "ARRAY" A;
     "CODE"31250;
     "PROCEDURE" POLCHS(N,A);
     "VALUE" N; "INTEGER" N; "ARRAY" A;
     "CODE"31051;
     LINTFMPOL(.5,.5,N,A);
     POLCHS(N,A);
 "END" POLSHTCHS;
         "EOP"
 
 
 
 "CODE" 31054;
 "PROCEDURE" SHTCHSPOL(N,A);
 "VALUE" N; "INTEGER" N; "ARRAY" A;
 "BEGIN"
     "PROCEDURE" CHSPOL(N,A);
     "VALUE" N; "INTEGER" N; "ARRAY" A;
     "CODE"31052;
     "PROCEDURE" LINTFMPOL(P,Q,N,A);
     "VALUE" N,P,Q; "INTEGER" N; "REAL" P,Q; "ARRAY" A;
     "CODE"31250;
     CHSPOL(N,A);
     LINTFMPOL(2,-1,N,A);
 "END" SHTCHSPOL;
         "EOP"
 
 
 
 "CODE" 31055;
 "PROCEDURE" GRNNEW(N,X,A);
 "VALUE" N; "INTEGER" N; "ARRAY" X,A;
 "BEGIN"
     "PROCEDURE" ELMCEV(L,U,SHIFT,A,B,X);
     "VALUE" L,U,SHIFT,X; "INTEGER" L,U,SHIFT;
     "REAL" X; "ARRAY" A,B;
     "FOR" L:=L "STEP" -1 "UNTIL" U "DO" A[L]:=A[L]+B[L+SHIFT]*X;
     "INTEGER" K;
     "FOR" K:=N-1 "STEP" -1 "UNTIL" 0 "DO"
     ELMCEV(N-1,N-1-K,1,A,A,X[N-1-K]);
 "END" GRNNEW;
         "EOP"
1SECTION : 2.4.1             (DECEMBER 1978)                     PAGE 11
 
 
 
 "CODE" 31050;
 "PROCEDURE" NEWGRN(N,X,A);
 "VALUE" N; "INTEGER" N; "ARRAY" X,A;
 "BEGIN"
     "PROCEDURE" ELMVEC(L,U,SHIFT,A,B,X);
     "VALUE" L,U,SHIFT,X; "INTEGER" L,U,SHIFT;
     "REAL" X; "ARRAY" A,B;
     "CODE" 34020;
     "INTEGER" K;
     "FOR" K:=N-1 "STEP" -1 "UNTIL" 0 "DO"
     ELMVEC(K,N-1,1,A,A,-X[K]);
 "END" NEWGRN;
         "EOP"
 
 
 
 "CODE" 31250;
 "PROCEDURE" LINTFMPOL(P,Q,N,A);
 "VALUE" N,P,Q; "INTEGER" N; "REAL" P,Q; "ARRAY" A;
 "BEGIN"
     "PROCEDURE" NORDERPOL(N,K,X,A);
     "VALUE" N,K,X; "INTEGER" N,K; "REAL" X; "ARRAY" A;
     "CODE" 31242;
     "INTEGER" K;
     "REAL" PPOWER;
     NORDERPOL(N,N,Q,A);
     PPOWER:=P;
     "FOR" K:=1 "STEP" 1 "UNTIL" N "DO"
     "BEGIN"
         A[K]:=PPOWER*A[K];
         PPOWER:=P*PPOWER;
     "END";
 "END" LINTFMPOL;
         "EOP"
1SECTION: 3.1.1.1.1.1.1       (FEBRUARY 1979)                     PAGE 1
 
 
 
 AUTHORS: J. C. P. BUS AND T. J. DEKKER.
 
 
 CONTRIBUTOR: J.C.P. BUS AND P. A. BEENTJES.
 
 
 INSTITUTE: MATHEMATICAL CENTRE.
 
 
 RECEIVED: 730830.
 
 
 BRIEF DESCRIPTION:
     THIS SECTION  CONTAINS  SIX  PROCEDURES:
     A: DEC PERFORMS A TRIANGULAR DECOMPOSITION WITH PARTIAL PIVOTING;
     B: GSSELM PERFORMS A TRIANGULAR DECOMPOSITION WITH A COMBINATION OF
        PARTIAL AND COMPLETE PIVOTING;
     C: ONENRMINV DELIVERS THE 1-NORM OF THE INVERSE OF A MATRIX WHOSE
        TRIANGULARLY DECOMPOSED FORM HAS BEEN DELIVERED BY DEC OR GSSELM;
     D: ERBELM CALCULATES A ROUGH UPPERBOUND FOR THE SOLUTION OF A LINEAR
        SYSTEM WHOSE MATRIX IS TRIANGULARLY DECOMPOSED BY GSSELM;
     E: GSSERB PERFORMS A TRIANGULAR DECOMPOSTION OF THE MATRIX OF A
        LINEAR SYSTEM AND CALCULATES AN UPPERBOUND FOR THE RELATIVE ERROR
        OF THE SOLUTION OF THAT SYSTEM;
     F: GSSNRI PERFORMS A TRIANGULAR DECOMPOSITION AND CALCULATES THE
        1-NORM OF THE INVERSE MATRIX;
 
     THE METHOD USED IN DEC AND GSSELM YIELDS A LOWER-TRIANGULAR MATRIX
     L AND A UNIT UPPER-TRIANGULAR MATRIX U SUCH THAT THE PRODUCT
     LU EQUALS THE GIVEN MATRIX WITH PERMUTED ROWS (DEC) OR ROWS AND
     COLUMNS (GSSELM); IN DEC, ONLY PARTIAL PIVOTING IS USED ([3],
     [4, P.115], [5, P.201]); THE PIVOTING STRATEGY IN GSSELM IS
     A COMBINATION OF PARTIAL AND COMPLETE PIVOTING ([2], [1]); IN THIS
     STRATEGY THE PROCESS WILL SWITCH TO COMPLETE PIVOTING IF PARTIAL
     PIVOTING MIGHT NOT YIELD STABLE RESULTS; SO IN GSSELM THE
     EFFICIENCY OF PARTIAL PIVOTING IS COMBINED WITH THE STABILITY OF
     COMPLETE PIVOTING;
     SINCE, IN EXCEPTIONAL CASES, PARTIAL PIVOTING MAY YIELD USELESS
     RESULTS, EVEN FOR WELL-CONDITIONED MATRICES, THE USER IS ADVISED TO
     USE GSSELM; HOWEVER, IF THE NUMBER OF VARIABLES IS SMALL RELATIVE
     TO THE NUMBER OF BINARY DIGITS IN THE MANTISSA ( 48 FOR THE CYBER ),
     THEN DEC MAY ALSO BE USED;
 
 
 KEYWORDS:
 
     LU DECOMPOSITION,
     TRIANGULAR DECOMPOSITION,
     GAUSSIAN ELIMINATION.
 
 
 
1SECTION: 3.1.1.1.1.1.1       (FEBRUARY 1979)                     PAGE 2
 
 
 
 SUBSECTION: DEC.
 
 
 CALLING SEQUENCE:
 
     THE HEADING OF THIS PROCEDURE IS:
     "PROCEDURE" DEC(A, N, AUX, P); "VALUE" N;
     "INTEGER" N; "ARRAY" A, AUX; "INTEGER" "ARRAY" P;"CODE" 34300;
 
     THE MEANING OF THE FORMAL PARAMETERS IS:
     A:      <ARRAY IDENTIFIER>;
             "ARRAY" A[1:N,1:N];
             ENTRY:THE MATRIX;
             EXIT:THE  CALCULATED  LOWER-TRIANGULAR  MATRIX  AND  UNIT
             UPPERTRIANGULAR  MATRIX  WITH ITS  UNIT  DIAGONAL  OMITTED;
     N:      <ARITHMETIC EXPRESSION>;
             THE ORDER OF THE MATRIX;
     AUX:    <ARRAY IDENTIFIER>;
             "ARRAY" AUX[1:3];
             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; SEE
                     METHOD AND PERFORMANCE;
             EXIT:
             AUX[1]: IF  R  IS THE NUMBER OF ELIMINATION STEPS PERFORMED
                     (SEE  AUX[3]),   THEN   AUX[1]  EQUALS  1   IF  THE
                     DETERMINANT OF THE PRINCIPAL SUBMATRIX OF  ORDER  R
                     IS POSITIVE, ELSE AUX[1] = -1;
             AUX[3]: THE  NUMBER  OF  ELIMINATION  STEPS  PERFORMED;  IF
                     AUX[3] < N THEN THE PROCESS IS BROKEN OFF BECAUSE
                     THE  SELECTED  PIVOT  IS  TOO  SMALL  RELATIVE TO
                     THE MAXIMUM OF THE EUCLIDEAN NORMS OF THE ROWS OF
                     THE GIVEN MATRIX;
     P:      <ARRAY IDENTIFIER>;
             "INTEGER""ARRAY" P[1:N];
             EXIT:THE PIVOTAL INDICES.
 
1SECTION: 3.1.1.1.1.1.1       (FEBRUARY 1979)                     PAGE 3
 
 
 
 PROCEDURES USED:
 
     MATMAT = CP34013,
     MATTAM = CP34015,
     ICHROW = CP34032.
 
 
 REQUIRED CENTRAL MEMORY:
 
     A REAL ARRAY OF ORDER N IS DECLARED.
 
 
 RUNNING TIME: PROPORTIONAL TO N ** 3.
 
 
 METHOD AND PERFORMANCE:
 
     THE METHOD USED IN DEC IS TRIANGULAR DECOMPOSITION WITH STABILIZING
     ROW INTERCHANGES, ALSO CALLED "PARTIAL PIVOTING"; SEE ALSO [3,P.19]
     AND [5,P.201].
 
 
 EXAMPLE OF USE: SEE DECSOL (SECTION 3.1.1.1.1.1.3).
 
 
 
1SECTION: 3.1.1.1.1.1.1       (FEBRUARY 1979)                     PAGE 4
 
 
 
 SUBSECTION: GSSELM   .
 
 
 CALLING SEQUENCE:
 
     THE HEADING OF THIS PROCEDURE IS:
     "PROCEDURE" GSSELM(A, N, AUX, RI, CI); "VALUE" N;
     "INTEGER" N; "ARRAY" A, AUX; "INTEGER" "ARRAY" RI, CI;"CODE" 34231;
 
     THE MEANING OF THE FORMAL PARAMETERS IS:
     A:      <ARRAY IDENTIFIER>;
             "ARRAY" A[1:N,1:N];
             ENTRY:THE N-TH ORDER MATRIX;
             EXIT:THE  CALCULATED  LOWER-TRIANGULAR  MATRIX  AND  UNIT
             UPPERTRIANGULAR  MATRIX  WITH ITS UNIT  DIAGONAL  OMITTED;
     N:      <ARITHMETIC EXPRESSION>;
             THE ORDER OF THE MATRIX;
     AUX:    <ARRAY IDENTIFIER>;
             "ARRAY" AUX[1:7];
             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; SEE
                     METHOD AND PERFORMANCE;
             AUX[4]: A  VALUE  WHICH  IS USED  FOR CONTROLLING PIVOTING;
                     USUALLY, AUX[4] = 8 WILL GIVE GOOD RESULTS; SEE
                     METHOD AND PERFORMANCE;
             EXIT:
             AUX[1]: IF  R  IS THE NUMBER OF ELIMINATION STEPS PERFORMED
                     (SEE  AUX[3]),   THEN   AUX[1]  EQUALS   1  IF  THE
                     DETERMINANT OF THE PRINCIPAL SUBMATRIX  OF ORDER  R
                     IS POSITIVE, ELSE AUX[1] = -1;
             AUX[3]: THE  NUMBER  OF  ELIMINATION  STEPS  PERFORMED;  IF
                     AUX[3] < N THEN THE PROCESS HAS BEEN BROKEN OFF,
                     BECAUSE THE SELECTED PIVOT IS TOO SMALL RELATIVE TO
                     THE MAXIMUM OF THE MODULI OF ELEMENTS OF THE GIVEN
                     MATRIX;
             AUX[5]: THE  MODULUS  OF  AN  ELEMENT  WHICH  IS OF MAXIMUM
                     ABSOLUTE VALUE FOR THE MATRIX WHICH HAD BEEN GIVEN
                     IN ARRAY A;
             AUX[7]: AN UPPER BOUND FOR THE GROWTH (I. E. THE MODULUS OF
                     AN ELEMENT  WHICH IS OF MAXIMUM ABSOLUTE VALUE  FOR
                     THE MATRICES OCCURRING DURING ELIMINATION);
     RI:     <ARRAY IDENTIFIER>;
             "INTEGER""ARRAY" RI[1:N];
             EXIT: THE PIVOTAL ROW-INDICES;
     CI:     <ARRAY IDENTIFIER>;
             "INTEGER""ARRAY" CI[1:N];
             EXIT: THE PIVOTAL COLUMN-INDICES;
 
 
1SECTION: 3.1.1.1.1.1.1       (DECEMBER 1979)                     PAGE 5
 
 
 
 PROCEDURES USED:
 
     ROWCST      = CP31132,
     ELMROW      = CP34024,
     MAXELMROW   = CP34025,
     ICHCOL      = CP34031,
     ICHROW      = CP34032,
     ABSMAXMAT   = CP31069.
 
 REQUIRED CENTRAL MEMORY: NO EXTRA ARRAYS ARE DECLARED.
 
 RUNNING TIME: PROPORTIONAL TO N ** 3.
 
 METHOD AND PERFORMANCE:
     THE PROCESS OF GAUSSIAN ELIMINATION IS PERFORMED IN AT MOST N STEPS
     , WHERE N DENOTES THE ORDER OF THE MATRIX; PARTIAL PIVOTING WILL BE
     USED AS LONG AS THE CALCULATED UPPER BOUND FOR THE GROWTH ([2],
     [1]), IS LESS THAN A CRITICAL VALUE THAT EQUALS AUX[4] * N TIMES
     THE MODULUS OF AN ELEMENT WHICH IS OF MAXIMUM ABSOLUTE VALUE FOR
     THE GIVEN MATRIX;
     IN THE PARTIAL PIVOTING STRATEGY, THAT ELEMENT IS CHOSEN AS
     PIVOT IN THE K-TH STEP, WHOSE ABSOLUTE VALUE IS MAXIMAL
     FOR THE K-TH COLUMN OF THE LOWER-TRIANGULAR MATRIX L; HOWEVER, IF
     THE UPPER BOUND FOR THE GROWTH EXCEEDS THIS CRITICAL VALUE IN THE
     K-TH STEP, THEN A PIVOT IS SELECTED IN THE J-TH STEP (J = K,..., N)
     , IN SUCH A WAY, THAT ITS ABSOLUTE VALUE IS MAXIMAL FOR THE
     REMAINING SUBMATRIX OF ORDER N - K + 1 (COMPLETE PIVOTING);
     SINCE  IN PRACTICE,  IF WE CHOOSE  AUX[4] PROPERLY, THE UPPER BOUND
     FOR THE GROWTH RARELY EXCEEDS  THIS CRITICAL VALUE  ([2], [4]),  WE
     WILL  USUALLY  TAKE  ADVANTAGE OF THE GREATER SPEED OF PARTIAL
     PIVOTING (ORDER N - K + 1 IN THE K-TH STEP), WHILE IN A FEW
     DOUBTFUL CASES NUMERICAL DIFFICULTIES WILL BE RECOGNIZED AND THE
     PROCESS WILL SWITCH TO COMPLETE PIVOTING (ORDER
     (N - K + 1) ** 2 IN THE K-TH STEP);
     USING GSSELM,THE UPPER BOUND FOR THE RELATIVE ERROR IN THE SOLUTION
     OF A LINEAR SYSTEM ([4], [5]), WILL BE AT MOST AUX[4] * N TIMES THE
     UPPER BOUND USING GAUSSIAN ELIMINATION WITH COMPLETE PIVOTING ONLY;
     USUALLY, HOWEVER, THIS WILL BE A CRUDE OVERESTIMATE;
     THE CHOICE AUX[4] < 1 / N WILL RESULT IN COMPLETE PIVOTING ONLY
     , WHILE  PARTIAL PIVOTING WILL BE USED IN EVERY STEP IF WE CHOOSE
     AUX[4] > (2 ** (N - 1)) / N; USUALLY, AUX[4] = 8 WILL GIVE
     GOOD RESULTS ([2], [1]);
     THE PROCESS WILL ALSO SWITCH TO COMPLETE PIVOTING IF THE MODULUS OF
     THE PIVOT OBTAINED WITH PARTIAL PIVOTING IS LESS THAN A CERTAIN
     TOLERANCE, WHICH EQUALS THE GIVEN RELATIVE TOLERANCE AUX[2] TIMES
     THE MODULUS OF AN ELEMENT WHICH IS OF MAXIMUM ABSOLUTE VALUE FOR
     THE GIVEN MATRIX; IF ALL ELEMENTS IN THE REMAINING SUBMATRIX ARE
     SMALLER IN ABSOLUTE VALUE THAN THIS TOLERANCE THEN THE PROCESS IS
     BROKEN OFF AND THE PREVIOUS STEPNUMBER IS DELIVERED IN AUX[3];
     IN  CONTRAST  WITH  THE  METHOD  USED  IN  DEC  (THIS SECTION),  NO
     EQUILIBRATING IS DONE  IN THIS PIVOTING STRATEGY;  THE USER HIMSELF
     HAS TO TAKE CARE FOR A REASONABLE SCALING OF THE MATRIX ELEMENTS.
 
 
 EXAMPLE OF USE: SEE GSSSOL (SECTION 3.1.1.1.1.1.3).
1SECTION: 3.1.1.1.1.1.1       (FEBRUARY 1979)                     PAGE 6
 
 
 
 SUBSECTION: ONENRMINV.
 
 
 CALLING SEQUENCE:
 
     THE HEADING OF THIS PROCEDURE IS
     "REAL" "PROCEDURE" ONENRMINV(A, N); "VALUE" N;
     "INTEGER" N; "ARRAY" A; "CODE" 34240;
 
     ONENRMINV:= THE  1-NORM  OF THE CALCULATED  INVERSE  OF THE MATRIX,
             WHOSE TRIANGULARLY DECOMPOSED FORM IS GIVEN IN ARRAY A;
 
     THE MEANING OF THE FORMAL PARAMETERS IS:
     A:      <ARRAY IDENTIFIER>;
             "ARRAY" A[1:N,1:N];
             ENTRY:
             THE TRIANGULARLY DECOMPOSED FORM OF A MATRIX,  AS DELIVERED
             BY  GSSELM  OR  DEC  (THIS SECTION);
     N:      <ARITHMETIC EXPRESSION>;
             THE ORDER  OF  THE MATRIX,  WHOSE  TRIANGULARLY  DECOMPOSED
             FORM HAS BEEN GIVEN IN ARRAY A.
 
 
 PROCEDURES USED:
 
     MATVEC  = CP34011.
 
 
 REQUIRED CENTRAL MEMORY:
 
     ONE REAL ARRAY OF ORDER N IS DECLARED.
 
 
 RUNNING TIME:   PROPORTIONAL TO N ** 3.
 
 
 METHOD AND PERFORMANCE:
 
     THE INVERSE OF THE MATRIX  WHOSE  TRIANGULARLY DECOMPOSED FORM,  AS
     DELIVERED BY GSSELM OR DEC, HAS BEEN GIVEN IN ARRAY A, IS CALCULATED
     WITH FORWARD AND BACK SUBSTITUTION ([3],[4],[5]); ONLY THE 1-NORM OF
     THIS INVERSE IS DELIVERED BY  ONENRMINV;  THE ELEMENTS  OF ARRAY  A
     REMAIN UNALTERED.
 
 
 EXAMPLE OF USE: SEE GSSSOLERB (SECTION 3.1.1.1.1.1.3).
 
 
 
1SECTION: 3.1.1.1.1.1.1       (FEBRUARY 1979)                     PAGE 7
 
 
 
 SUBSECTION: ERBELM.
 
 
 CALLING SEQUENCE:
 
     THE HEADING OF THIS PROCEDURE IS:
     "PROCEDURE" ERBELM(N, AUX, NRMINV); "VALUE" N, NRMINV;
     "INTEGER" N; "REAL" NRMINV; "ARRAY" AUX; "CODE" 34241;
 
     THE MEANING OF THE FORMAL PARAMETERS IS:
     N:      <ARITHMETIC EXPRESSION>;
             THE ORDER OF THE LINEAR SYSTEM IN CONSIDERATION;
     AUX:    <ARRAY IDENTIFIER>;
             "ARRAY" AUX[0:11];
             ENTRY:
             AUX[0]: THE MACHINE PRECISION;
             AUX[5]: THE MODULUS  OF  AN ELEMENT,  WHICH  IS  OF MAXIMUM
                     ABSOLUTE VALUE FOR THE MATRIX OF THE LINEAR SYSTEM;
                     THIS VALUE IS DELIVERED  BY GSSELM  IN AUX[5] (THIS
                     SECTION);
             AUX[6]: AN  UPPER BOUND  FOR  THE  RELATIVE  ERROR  IN  THE
                     ELEMENTS OF THE MATRIX OF THE LINEAR SYSTEM;
             AUX[7]: AN  UPPER BOUND  FOR  THE  GROWTH  DURING  GAUSSIAN
                     ELIMINATION;  THIS VALUE IS DELIVERED IN  AUX[7] BY
                     GSSELM (THIS SECTION);
             EXIT:
             AUX[9]: THE VALUE OF NRMINV;
             AUX[11]: A ROUGH UPPER BOUND  FOR THE RELATIVE ERROR IN THE
                     SOLUTION   OF  A  LINEAR  SYSTEM    WHEN   GAUSSIAN
                     ELIMINATION  IS USED  FOR  THE CALCULATION  OF THIS
                     SOLUTION; IF NO USE CAN BE MADE OF THE FORMULA  FOR
                     THE  ERROR BOUND  (SEE:  METHOD  AND  PERFORMANCE),
                     BECAUSE OF A VERY BAD CONDITION OF THE MATRIX, THEN
                     AUX[11]:= -1;
     NRMINV: <ARITHMETIC EXPRESSION>;
             THE 1-NORM  OF  THE INVERSE  OF THE MATRIX  OF  THE  LINEAR
             SYSTEM MUST BE GIVEN IN NRMINV; THIS VALUE MAY BE OBTAINED
             BY ONENRMINV (THIS SECTION).
1SECTION: 3.1.1.1.1.1.1       (FEBRUARY 1979)                     PAGE 8
 
 
 
 PROCEDURES USED: NONE.
 
 
 REQUIRED CENTRAL MEMORY: NO EXTRA ARAYS ARE DECLARED.
 
 
 METHOD AND PERFORMANCE:
 
     WHEN CALLED AFTER GSSELM, ERBELM WILL CALCULATE A ROUGH UPPER BOUND
     FOR THE RELATIVE ERROR IN THE SOLUTION OF THE LINEAR SYSTEM,  WHOSE
     MATRIX HAS BEEN DECOMPOSED INTO TRIANGULAR FORM BY GSSELM (THIS
     SECTION), BY  ([3], [4], [5]):
             NORM(DX) / NORM(X) <= P / (1 - P),
     WHERE : P = Q * NORM(C) / (1 - Q * NORM(C)),
             Q = G * (.75 * N ** 3 + 4.5 * N ** 2) * EPS + EPSA,
             C IS THE CALCULATED INVERSE OF THE MATRIX,
             G   THE  UPPER  BOUND  FOR   THE  GROWTH   DURING  GAUSSIAN
             ELIMINATION, AS DELIVERED BY GSSELM (THIS SECTION),
             N THE ORDER OF THE MATRIX,
             EPSA  AN UPPER BOUND  FOR  THE RELATIVE ERROR IN THE MATRIX
             ELEMENTS,
             EPS THE MACHINE PRECISION AND
             NORM(.) DENOTES THE 1-NORM.
 
     THIS PROCEDURE IS USED IN E.G. GSSERB (THIS SECTION) AND GSSINVERB
     (SECTION 3.1.1.1.1.1.4)
 
 
 EXAMPLE OF USE: SEE GSSSOLERB (SECTION 3.1.1.1.1.1.3).
 
 
 
1SECTION: 3.1.1.1.1.1.1       (FEBRUARY 1979)                     PAGE 9
 
 
 
 SUBSECTION: GSSERB    .
 
 CALLING SEQUENCE:
 
     THE HEADING OF THIS PROCEDURE IS:
     "PROCEDURE" GSSERB(A, N, AUX, RI, CI); "VALUE" N;
     "INTEGER" N; "ARRAY" A, AUX; "INTEGER" "ARRAY" RI, CI; "CODE" 34242;
 
     THE MEANING OF THE FORMAL PARAMETERS IS:
     A:      <ARRAY IDENTIFIER>;
             "ARRAY" A[1:N,1:N];
             ENTRY:THE MATRIX TO BE DECOMPOSED;
             EXIT: THE CALCULATED LOWER-TRIANGULAR MATRIX AND UNIT UPPER
             TRIANGULAR MATRIX,WITH ITS UNIT DIAGONAL OMITTED;
     N:      <ARITHMETIC EXPRESSION>;
             THE ORDER OF THE MATRIX;
     AUX:    <ARRAY IDENTIFIER>;
             "ARRAY" AUX[0:11];
             ENTRY:(SEE ALSO GSSELM IN THIS SECTION);
             AUX[0]: THE MACHINE PRECISION;
             AUX[2]: A RELATIVE TOLERANCE;
             AUX[4]: A VALUE WHICH IS USED FOR CONTROLLING PIVOTING;
             AUX[6]: AN UPPER BOUND  FOR  THE RELATIVE PRECISION  OF THE
                     MATRIX ELEMENTS;
             EXIT:
             AUX[1]: IF  R IS THE NUMBER OF ELIMINATION STEPS PERFORMED,
                     AUX[1] EQUALS 1 IF THE DETERMINANT OF THE PRINCIPAL
                     SUBMATRIX OF ORDER R IS POSITIVE, ELSE AUX[1] = -1;
             AUX[3]: THE NUMBER OF ELIMINATION STEPS PERFORMED;
             AUX[5]: THE  MODULUS  OF  AN ELEMENT,  WHICH  IS OF MAXIMUM
                     ABSOLUTE VALUE FOR THE MATRIX WHICH HAS BEEN GIVEN
                     IN ARRAY A;
             AUX[7]: AN UPPER BOUND FOR THE GROWTH;
             AUX[9]: IF AUX[3] = N, THEN AUX[9] WILL EQUAL THE 1-NORM OF
                     THE INVERSE MATRIX, ELSE AUX[9] WILL BE UNDEFINED;
             AUX[11]: IF AUX[3] = N, THEN THE VALUE OF AUX[11] WILL BE A
                     ROUGH  UPPER BOUND  FOR  THE RELATIVE ERROR  IN THE
                     SOLUTION OF LINEAR SYSTEMS  WITH A MATRIX  AS GIVEN
                     IN ARRAY A, ELSE AUX[11] WILL BE UNDEFINED;  IF  NO
                     USE CAN BE MADE OF THE FORMULA FOR  THE ERROR BOUND
                     AS GIVEN ABOVE (SUBSECTION  ERBELM),  BECAUSE OF  A
                     VERY BAD CONDITION OF THE MATRIX, THEN AUX[11]:=-1;
     RI:     <ARRAY IDENTIFIER>;
             "INTEGER""ARRAY" RI[1:N];
             EXIT: THE PIVOTAL ROW-INDICES.
     CI:     <ARRAY IDENTIFIER>;
             "INTEGER""ARRAY" CI[1:N];
             EXIT: THE PIVOTAL COLUMN-INDICES.
 
 PROCEDURES USED:
 
     GSSELM      = CP34231,
     ONENRMINV   = CP34240,
     ERBELM      = CP34241.
1SECTION: 3.1.1.1.1.1.1       (FEBRUARY 1979)                    PAGE 10
 
 
 
 REQUIRED CENTRAL MEMORY: NO EXTRA ARRAYS ARE DECLARED.
 
 
 RUNNING TIME: PROPORTIONAL TO N ** 3.
 
 
 METHOD AND PERFORMANCE:
 
     GSSERB  USES  GSSELM  (THIS  SECTION)  TO  PERFORM  THE  TRIANGULAR
     DECOMPOSITION  OF THE MATRIX  GIVEN  IN  ARRAY  A  AND  ERBELM  AND
     ONENRMINV  (THIS  SECTION) TO  CALCULATE  AN  UPPER BOUND  FOR  THE
     RELATIVE ERROR IN THE SOLUTION OF LINEAR SYSTEMS  WITH A MATRIX  AS
     GIVEN IN ARRAY  A; IF  AUX[3] < N,  THEN  THE EFFECT OF  GSSERB  IS
     MERELY THAT OF GSSELM.
 
 EXAMPLE OF USE: SEE GSSSOLERB (SECTION 3.1.1.1.1.1.3).
 
 
 
 SUBSECTION: GSSNRI   .
 
 
 CALLING SEQUENCE:
 
     THE HEADING OF THIS PROCEDURE IS:
     "PROCEDURE" GSSNRI(A, N, AUX, RI, CI); "VALUE" N;
     "INTEGER" N; "ARRAY" A, AUX; "INTEGER" "ARRAY" RI, CI; "CODE" 34252;
 
     THE MEANING OF THE FORMAL PARAMETERS IS:
     A:      <ARRAY IDENTIFIER>;
             "ARRAY" A[1:N,1:N];
             ENTRY:THE MATRIX TO BE DECOMPOSED;
             EXIT: THE CALCULATED LOWER-TRIANGULAR MATRIX AND UNIT UPPER
             TRIANGULAR MATRIX,WITH ITS DIAGONAL OMITTED;
     N:      <ARITHMETIC EXPRESSION>;
             THE ORDER OF THE MATRIX;
     AUX:    <ARRAY IDENTIFIER>;
             "ARRAY" AUX[1:9];
             ENTRY:(SEE ALSO GSSELM IN THIS SECTION);
             AUX[2]: A RELATIVE TOLERANCE;
             AUX[4]: A VALUE USED FOR CONTROLLING PIVOTING;
             EXIT:
             AUX[1]: IF R  IS THE NUMBER OF ELIMINATION STEPS PERFORMED,
                     THEN AUX[1]  EQUALS  1  IF THE DETERMINANT  OF  THE
                     PRINCIPAL SUBMATRIX OF ORDER  R  IS POSITIVE,  ELSE
                     AUX[1] = -1;
             AUX[3]: THE NUMBER OF ELIMINATION STEPS PERFORMED;
             AUX[5]: THE  MODULUS  OF  AN  ELEMENT,  WHICH IS OF MAXIMUM
                     ABSOLUTE VALUE FOR THE MATRIX WHICH HAD BEEN GIVEN
                     IN ARRAY A;
             AUX[7]: AN UPPER BOUND FOR THE GROWTH;
             AUX[9]: IF AUX[3] = N, THEN AUX[9] WIL EQUAL  THE 1-NORM OF
                     THE INVERSE MATRIX, ELSE AUX[9] WILL BE UNDEFINED;
1SECTION: 3.1.1.1.1.1.1       (FEBRUARY 1979)                    PAGE 11
 
 
 
     RI:     <ARRAY IDENTIFIER>;
             "INTEGER""ARRAY" RI[1:N];
             EXIT: THE PIVOTAL ROW INDICES.
     CI:     <ARRAY IDENTIFIER>;
             "INTEGER""ARRAY" CI[1:N];
             EXIT:THE PIVOTAL COLUMN INDICES.
 
 
 PROCEDURES USED:
 
     GSSELM      = CP34231,
     ONENRMINV   = CP34240.
 
 
 REQUIRED CENTRAL MEMORY: NO EXTRA ARRAYS ARE DECLARED.
 
 
 RUNNING TIME: PROPORTIONAL TO N ** 3.
 
 
 METHOD AND PERFORMANCE:
 
     GSSNRI  USES  GSSELM  (THIS  SECTION)  TO  PERFORM  THE  TRIANGULAR
     DECOMPOSITION OF THE MATRIX GIVEN IN ARRAY  A  AND  ONENRMINV (THIS
     SECTION)  TO  CALCULATE  THE  1-NORM  OF  THE  INVERSE  MATRIX;  IF
     AUX[3] < N,  THEN THE EFFECT OF  GSSNRI  IS  MERELY THAT OF  GSSELM
     (THIS SECTION).
 
 
 EXAMPLE OF USE: SEE GSSITISOLERB (SECTION 3.1.1.1.1.1.5).
 
 
 
 REFERENCES:
 
     [1] BUS, J. C. P.
         LINEAR SYSTEMS  WITH CALCULATION OF ERROR BOUNDS  AND ITERATIVE
         REFINEMENT (DUTCH).
         MATHEMATICAL CENTRE, AMSTERDAM, LR 3. 4. 19 (1972).
     [2] BUSINGER, P. A.
         MONITORING THE NUMERICAL STABILITY OF GAUSSIAN ELIMINATION.
         NUMER. MATH. 16, 360-361 (1971).
     [3] DEKKER, T. J.
         ALGOL 60 PROCEDURES IN NUMERICAL ALGEBRA, PART 1.
         MATHEMATICAL CENTRE, AMSTERDAM, TRACT 22 (1968).
     [4] WILKINSON, J. H.
         ROUNDING ERRORS IN ALGEBRAIC PROCESSES.
         LONDON (1963).
     [5] WILKINSON, J. H.
         THE ALGEBRAIC EIGENVALUE PROBLEM.
         OXFORD (1965).
 
 
1SECTION: 3.1.1.1.1.1.1       (MAY 1974)                         PAGE 12
 
 
 
 SOURCE TEXT(S):
 
 "CODE" 34300;
     "PROCEDURE" DEC(A, N, AUX, P); "VALUE" N; "INTEGER" N;
     "ARRAY" A, AUX; "INTEGER" "ARRAY" P;
     "BEGIN" "INTEGER" I, K, K1, PK, D;
         "REAL" R, S, EPS;
         "ARRAY" V[1:N];
         "REAL" "PROCEDURE" MATMAT(L, U, I, J, A, B); "CODE" 34013;
         "REAL" "PROCEDURE" MATTAM(L, U, I, J, A, B); "CODE" 34015;
         "PROCEDURE" ICHROW(L, U, I, J, A); "CODE" 34032;
         R:= -1;
         "FOR" I:= 1 "STEP" 1 "UNTIL" N "DO"
         "BEGIN" S:= SQRT(MATTAM(1,N,I,I,A,A));
             "IF" S > R "THEN" R:= S; V[I]:= 1/S
         "END";
         EPS:= AUX[2] * R; D:= 1;
         "FOR" K:= 1 "STEP" 1 "UNTIL" N "DO"
         "BEGIN" R:= - 1; K1:= K - 1;
             "FOR" I:= K "STEP" 1 "UNTIL" N "DO"
             "BEGIN" A[I,K]:= A[I,K] - MATMAT(1, K1, I, K, A, A);
                 S:= ABS(A[I,K]) * V[I]; "IF" S > R "THEN"
                 "BEGIN" R:= S; PK:= I "END"
             "END" LOWER;
             P[K]:= PK; V[PK]:= V[K]; S:= A[PK,K];
             "IF" ABS(S) < EPS "THEN" "GOTO" END;
             "IF" S < 0 "THEN" D:= -D; "IF" PK ^= K "THEN"
             "BEGIN" D:= - D; ICHROW(1, N, K, PK, A) "END";
             "FOR" I:= K + 1 "STEP" 1 "UNTIL" N "DO"
             A[K,I]:= (A[K,I] - MATMAT(1, K1, K, I, A, A)) / S
         "END" LU;
         K:= N + 1;
      END: AUX[1]:= D; AUX[3]:= K - 1
     "END" DEC;
         "EOP"
 
 "CODE" 34231;
     "PROCEDURE" GSSELM(A, N, AUX, RI, CI); "VALUE" N; "INTEGER" N;
     "ARRAY" A, AUX;
     "INTEGER" "ARRAY" RI, CI;
     "BEGIN" "INTEGER" I, J, P, Q, R, R1, JPIV, RANK, SIGNDET;
         "REAL" CRIT, PIVOT, RGROW, MAX, AID, MAX1, EPS;
         "BOOLEAN" PARTIAL;
         "PROCEDURE" ELMROW(L, U, I, J, A, B, X); "CODE" 34024;
         "INTEGER" "PROCEDURE" MAXELMROW(L, U, I, J, A, B, X);
         "CODE" 34025;
         "PROCEDURE" ICHCOL(L, U, I, J, A); "CODE" 34031;
         "PROCEDURE" ICHROW(L, U, I, J, A); "CODE" 34032;
                                                               "COMMENT"
1SECTION: 3.1.1.1.1.1.1       (DECEMBER 1979)                    PAGE 13
                                                                 ;
 
 
         "PROCEDURE" ROWCST(L, U, I, A, X); "CODE" 31132;
         "REAL" "PROCEDURE" ABSMAXMAT(LR, UR, LC, UC, I, J, A);
         "CODE" 31069;
         AUX[5]:= RGROW:= ABSMAXMAT(1, N, 1, N, I, J, A);
         CRIT:= N * RGROW * AUX[4]; EPS:= RGROW * AUX[2]; MAX:= 0;
         RANK:= N; SIGNDET:= 1; PARTIAL:= RGROW ^= 0;
         "FOR" Q:= 1 "STEP" 1 "UNTIL" N "DO" "IF" Q ^= J "THEN"
         "BEGIN" AID:= ABS(A[I,Q]);
             "IF" AID > MAX "THEN" MAX:= AID
         "END";
         RGROW:= RGROW + MAX;
         "FOR" R:= 1 "STEP" 1 "UNTIL" N "DO"
         "BEGIN" R1:= R + 1; "IF" I ^= R "THEN"
             "BEGIN" SIGNDET:= - SIGNDET; ICHROW(1, N, R, I, A) "END";
             "IF" J ^= R "THEN"
             "BEGIN" SIGNDET:= - SIGNDET; ICHCOL(1, N, R, J, A) "END";
             RI[R]:= I; CI[R]:= J; PIVOT:= A[R,R];
             "IF" PIVOT < 0 "THEN" SIGNDET:= - SIGNDET;
             "IF" PARTIAL "THEN"
             "BEGIN" MAX:= MAX1:= 0; J:= R1;
                 ROWCST(R1, N, R, A, 1 / PIVOT);
                 "FOR" P:= R1 "STEP" 1 "UNTIL" N "DO"
                 "BEGIN" ELMROW(R1, N, P, R, A, A, - A[P,R]);
                     AID:= ABS(A[P,R1]); "IF" MAX < AID "THEN"
                     "BEGIN" MAX:= AID; I:= P "END";
                  "END";
                 "FOR" Q:= R1 + 1 "STEP" 1 "UNTIL" N "DO"
                 "BEGIN" AID:= ABS(A[I,Q]);
                     "IF" MAX1 < AID "THEN" MAX1:= AID
                 "END";
                 AID:= RGROW; RGROW:= RGROW + MAX1;
                 "IF" RGROW > CRIT "OR" MAX < EPS "THEN"
                 "BEGIN" PARTIAL:= "FALSE"; RGROW:= AID;
                     MAX:= ABSMAXMAT(R1, N, R1, N, I, J, A)
                 "END"
             "END" PARTIAL PIVOTINGSTEP
             "ELSE"
             "BEGIN" "IF" MAX <= EPS "THEN"
                 "BEGIN" RANK:= R - 1;
                     "IF" PIVOT < 0 "THEN" SIGNDET:= - SIGNDET;"GOTO"OUT
                 "END";
                 MAX:= - 1;
                 ROWCST(R1, N, R, A, 1 / PIVOT);
                 "FOR" P:= R1 "STEP" 1 "UNTIL" N "DO"
                 "BEGIN" JPIV:= MAXELMROW(R1, N, P, R, A, A, - A[P,R]);
                     AID:= ABS(A[P,JPIV]); "IF" MAX < AID "THEN"
                     "BEGIN" MAX:= AID; I:= P; J:= JPIV "END"
                 "END";
                 "IF" RGROW < MAX "THEN" RGROW:= MAX
             "END" COMPLETE PIVOTINGSTEP
         "END" ELIMINATIONSTEP;
     OUT: AUX[1]:= SIGNDET; AUX[3]:= RANK; AUX[7]:= RGROW
     "END" GSSELM;
         "EOP"
1SECTION: 3.1.1.1.1.1.1       (MAY 1974)                         PAGE 14
 
 
 
 "CODE" 34240;
     "REAL" "PROCEDURE" ONENRMINV(A, N); "VALUE" N; "INTEGER" N;
     "ARRAY" A;
     "BEGIN" "INTEGER" I, J;
         "REAL" NORM, MAX, AID;
         "ARRAY" Y[1:N];
         "REAL" "PROCEDURE" MATVEC(L, U, I, A, B); "CODE" 34011;
         NORM:= 0;
         "FOR" J:= 1 "STEP" 1 "UNTIL" N "DO"
         "BEGIN" "FOR" I:= 1 "STEP" 1 "UNTIL" N "DO" Y[I]:= "IF" I < J
             "THEN" 0 "ELSE" "IF" I = J "THEN" 1 / A[I,I] "ELSE"
             - MATVEC(J, I - 1, I, A, Y) / A[I,I];
             MAX:= 0;
             "FOR" I:= N "STEP" - 1 "UNTIL" 1 "DO"
             "BEGIN" AID:= Y[I]:= Y[I] - MATVEC(I + 1, N, I, A, Y);
                 MAX:= MAX + ABS(AID)
             "END";
             "IF" NORM < MAX "THEN" NORM:= MAX
         "END";
         ONENRMINV:= NORM
     "END" ONENRMINV;
         "EOP"
 "CODE" 34241;
     "PROCEDURE" ERBELM(N, AUX, NRMINV); "VALUE" N, NRMINV;
     "INTEGER" N; "REAL" NRMINV;
     "ARRAY" AUX;
     "BEGIN" "REAL" AID, EPS;
         EPS:= AUX[0]; AID:= (1.06 * EPS * (.75 * N + 4.5) * N ** 2
         * AUX[7] + AUX[5] * AUX[6]) * NRMINV;
         AUX[11]:= "IF" 2 * AID >= (1 - EPS) "THEN" - 1 "ELSE"
         AID / (1 - 2 * AID); AUX[9]:= NRMINV
     "END" ERBELM;
         "EOP"
 "CODE" 34242;
     "PROCEDURE" GSSERB(A, N, AUX, RI, CI); "VALUE" N; "INTEGER" N;
     "ARRAY" A, AUX;
     "INTEGER" "ARRAY" RI, CI;
     "BEGIN" "PROCEDURE" GSSELM(A, N, AUX, RI, CI); "CODE" 34231;
         "REAL" "PROCEDURE" ONENRMINV(A, N); "CODE" 34240;
         "PROCEDURE" ERBELM(N, AUX, NRMINV); "CODE" 34241;
         GSSELM(A, N, AUX, RI, CI);
         "IF" AUX[3] = N "THEN" ERBELM(N, AUX, ONENRMINV(A, N))
     "END" GSSERB;
         "EOP"
 "CODE" 34252;
     "PROCEDURE" GSSNRI(A, N, AUX, RI, CI); "VALUE" N; "INTEGER" N;
     "ARRAY" A, AUX;
     "INTEGER" "ARRAY" RI, CI;
     "BEGIN" "PROCEDURE" GSSELM(A, N, AUX, RI, CI); "CODE" 34231;
         "REAL" "PROCEDURE" ONENRMINV(A, N); "CODE" 34240;
         GSSELM(A, N, AUX, RI, CI);
         "IF" AUX[3] = N "THEN" AUX[9]:= ONENRMINV(A, N)
     "END" GSSNRI;
         "EOP"
1SECTION: 3.1.1.1.1.1.2       (MAY 1974)                          PAGE 1      
        
        
        
 AUTHOR: J. C. P. BUS.      
        
        
 CONTRIBUTOR: J.C.P. BUS AND P. A. BEENTJES.    
        
        
 INSTITUTE: MATHEMATICAL CENTRE.      
        
        
 RECEIVED: 730831.
        
        
 BRIEF DESCRIPTION:         
     THIS SECTION  CONTAINS  A  PROCEDURE FOR CALCULATING THE       
     DETERMINANT OF A TRIANGULAR DECOMPOSED MATRIX;       
        
        
 KEYWORDS:        
        
     DETERMINANT. 
        
        
 CALLING SEQUENCE:
        
     THE HEADING OF THIS PROCEDURE IS:
     "REAL" "PROCEDURE" DETERM(A, N, SIGN); "VALUE" N, SIGN;        
     "INTEGER" N, SIGN; "ARRAY" A;    
        
     DETERM: DELIVERS THE CALCULATED VALUE OF THE DETERMINANT OF THE
             MATRIX;        
        
     THE MEANING OF THE FORMAL PARAMETERS IS:   
     A:      <ARRAY IDENTIFIER>;      
             "ARRAY" A[1:N, 1:N];     
             ENTRY:  THE DIAGONAL ELEMENTS OF THE         
                     LOWER-TRIANGULAR MATRIX L, OBTAINED BY TRIANGULAR        
                     DECOMPOSITION OF THE MATRIX, HAS TO BE GIVEN IN
                     A[I,I], I= 1, ..., N;      
     N:      <ARITHMETIC EXPRESSION>; 
             THE  ORDER  OF  THE  MATRIX,  WHOSE  DETERMINANT  HAS TO BE      
             CALCULATED;    
     SIGN:   <ARITHMETIC EXPRESSION>; 
             ENTRY:  IF THE DETERMINANT OF THE MATRIX  IS POSITIVE  THEN      
                     THE VALUE OF SIGN SHOULD BE +1, ELSE -1; THIS VALUE      
                     IS DELIVERED BY GSSELM OR DEC IN  AUX[1],  (SECTION      
                     3.1.1.1.1.1.1).  
        
        
 PROCEDURES USED: NONE.     
1SECTION: 3.1.1.1.1.1.2       (MAY 1974)                          PAGE 2      
        
        
        
 RUNNING TIME: PROPORTIONAL TO N.     
        
        
 LANGUAGE:   ALGOL 60.      
        
        
 METHOD AND PERFORMANCE:    
     A LOWER-TRIANGULAR MATRIX L HAS TO BE      
     GIVEN, SUCH THAT FOR SOME UNIT UPPER-TRIANGULAR MATRIX U THE   
     PRODUCT LU EQUALS THE MATRIX (WITH PERMUTED ROWS AND COLUMNS); 
     THE SIGN OF THE DETERMINANT ALSO HAS TO BE GIVEN; THESE DATA ARE         
     DELIVERED IN THE MATRIX AND AUX[1] BY THE PROCEDURES GSSELM OR DEC       
     (SECTION 3.1.1.1.1.1.1) AND THE PROCEDURES GSSERB, GSSNRI (SECTION       
     3.1.1.1.1.1.1), DECSOL, GSSSOL, GSSSOLERB (SECTION 3.1.1.1.1.1.3),       
     GSSITISOL AND GSSITISOLERB (SECTION 3.1.1.1.1.1.5), WHICH MAKE USE       
     OF GSSELM OR DEC.      
     THE  CALCULATION   OF  THE  DETERMINANT  IS  DONE  STRAIGTH  ON  BY      
     CALCULATING   THE  PRODUCT   OF   THE  DIAGONAL  ELEMENTS   OF  THE      
     LOWER-TRIANGULAR MATRIX GIVEN IN ARRAY A;  THE USER IS WARNED, THAT      
     OVERFLOW MAY OCCUR IF THE ORDER OF THE MATRIX IS LARGE.        
        
        
 EXAMPLE OF USE:  
        
     THE DETERMINANT  OF THE FOURTH ORDER SEGMENT  OF THE HILBERT MATRIX      
     MAY BE OBTAINED BY THE FOLLOWING PROGRAM:  
        
     "BEGIN" "INTEGER" I, J; "REAL" D; "INTEGER" "ARRAY" RI, CI[1:4];         
         "ARRAY" A[1:4, 1:4], AUX[1:7];         
         "PROCEDURE" GSSELM(A, N, AUX, RI, CI); "CODE" 34231;       
         "REAL" "PROCEDURE" DETERM(A, N, SIGN); "CODE" 34303;       
         "FOR" I:= 1 "STEP" 1 "UNTIL" 4 "DO"    
         "FOR" J:= 1 "STEP" 1 "UNTIL" 4 "DO" A[I, J]:= 1 / (I + J - 1);       
         AUX[2]:= "-14; AUX[4]:= 8;   
         GSSELM(A, 4, AUX, RI, CI);   
         D:= "IF" AUX[3] = 4 "THEN" DETERM(A, 4, AUX[1]) "ELSE" 0;  
         OUTPUT(71, "(""("DETERMINANT =")"B+.15D"+3D")", D)         
     "END"        
        
     RESULT:      
        
     DETERMINANT = +.165343915345370"-006       
        
        
 SOURCE TEXT(S):  
 "CODE" 34303;    
     "REAL" "PROCEDURE" DETERM(A, N, SIGN); "VALUE" N, SIGN;        
     "INTEGER" N, SIGN; "ARRAY" A;    
     "BEGIN" "INTEGER" I; "REAL" DET; 
         DET:= 1; 
         "FOR" I:= 1 "STEP" 1 "UNTIL" N "DO" DET:= A[I, I] * DET;   
         DETERM:= SIGN * ABS(DET)     
     "END" DETERM;
         "EOP"    
1SECTION: 3.1.1.1.1.1.3       (MAY 1974)                          PAGE 1      
        
        
        
 AUTHORS: J. C. P. BUS AND T. J. DEKKER.        
        
        
 CONTRIBUTOR: J.C.P. BUS AND P. A. BEENTJES.    
        
        
 INSTITUTE: MATHEMATICAL CENTRE.      
        
        
 RECEIVED: 730915.
        
        
 BRIEF DESCRIPTION:         
        
     THIS SECTION CONTAINS FIVE PROCEDURES:     
     SOL SOLVES THE LINEAR SYSTEM WHOSE MATRIX HAS BEEN TRIANGULARLY
     DECOMPOSED BY DEC;     
     DECSOL SOLVES A LINEAR SYSTEM WHOSE ORDER IS SMALL RELATIVE TO THE       
     NUMBER OF BINARY DIGITS IN THE NUMBER REPRSENTATION; 
     SOLELM SOLVES A LINEAR SYSTEM WHOSE MATRIX HAS BEEN TRIANGULARLY         
     DECOMPOSED BY GSSELM OR GSSERB(SECTION 3.1.1.1.1.1.1.).        
     GSSSOL SOLVES A LINEAR SYSTEM;   
     GSSSOLERB  SOLVES  A  LINEAR  SYSTEM  AND  CALCULATES  A  ROUGH
     UPPERBOUND FOR THE RELATIVE ERROR IN THE CALCULATED SOLUTION;  
        
     THE
     DIFFERENCE BETWEEN  DECSOL ON THE ONE SIDE AND GSSSOL AND GSSSOLERB      
     ON  THE  OTHER  SIDE  LIES  IN  THE  METHOD  USED   FOR  TRIANGULAR      
     DECOMPOSITION,  PARTICULARLY IN THE PIVOTING STRATEGY;  DECSOL USES      
     DEC,  GSSSOL  AND  GSSSOLERB  USE  GSSELM TO PERFORM THE TRIANGULAR      
     DECOMPOSITION (SECTION 3.1.1.1.1.1.1); SINCE, IN EXCEPTIONAL CASES,      
     DEC  MAY YIELD USELESS RESULTS,  ONE IS ADVISED  TO USE  GSSSOL  OR      
     GSSSOLERB; HOWEVER, IF THE ORDER OF THE LINEAR SYSTEM IS VERY SMALL      
     RELATIVE  TO   THE  NUMBER   OF   BINARY  DIGITS   IN   THE  NUMBER      
     REPRESENTATION, THEN DECSOL ALSO MAY BE USED.        
        
        
 KEYWORDS:        
        
     ALGEBRAIC EQUATIONS,   
     LINEAR SYSTEMS.        
        
        
        
1SECTION: 3.1.1.1.1.1.3       (MAY 1974)                          PAGE 2      
        
        
        
 SUBSECTION: SOL .
        
        
 CALLING SEQUENCE:
        
     THE HEADING OF THIS PROCEDURE IS:
     "PROCEDURE" SOL(A, N, P, B); "VALUE" N;    
     "INTEGER" N; "ARRAY" A, B; "INTEGER" "ARRAY" P;      
        
     THE MEANING OF THE FORMAL PARAMETERS IS:   
     A:      <ARRAY IDENTIFIER>;      
             "ARRAY" A[1:N, 1:N];     
             ENTRY:  THE TRIANGULARLY  DECOMPOSED FORM  OF THE MATRIX OF      
                     THE  LINEAR  SYSTEM  AS  PRODUCED BY  DEC  (SECTION      
                     3.1.1.1.1.1.1);  
     N:      <ARITHMETIC EXPRESSION>; 
             THE ORDER OF THE MATRIX; 
     P:      <ARRAY IDENTIFIER>;      
             "INTEGER""ARRAY" P[1:N]; 
             ENTRY:THE PIVOTAL INDICES, AS PRODUCED BY DEC.         
     B:      <ARRAY IDENTIFIER>;      
             "ARRAY" B[1:N];
             ENTRY:  THE RIGHT-HAND SIDE OF THE LINEAR SYSTEM;      
             EXIT:   THE SOLUTION OF THE LINEAR SYSTEM.   
        
        
 PROCEDURES USED: 
        
     MATVEC = CP34011.      
        
        
 RUNNING TIME: PROPORTIONAL TO N ** 2.
        
        
 LANGUAGE:   ALGOL 60.      
        
        
 METHOD AND PERFORMANCE:    
        
     SOL SHOULD BE CALLED AFTER  DEC (SECTION 3.1.1.1.1.1.1)  AND SOLVES      
     THE LINEAR SYSTEM WITH A MATRIX, WHOSE TRIANGULARLY DECOMPOSED FORM      
     AS PRODUCED BY  DEC  IS GIVEN IN ARRAY  A, AND A RIGHT-HAND SIDE AS      
     GIVEN IN ARRAY B; SOL LEAVES A AND P UNALTERED,  SO, AFTER ONE CALL      
     OF DEC, SEVERAL CALLS OF SOL MAY FOLLOW FOR SOLVING SEVERAL SYSTEMS      
     HAVING THE SAME MATRIX BUT DIFFERENT RIGHT-HAND SIDES.         
        
        
 EXAMPLE OF USE: SEE DECSOL (THIS SECTION).     
1SECTION: 3.1.1.1.1.1.3       (MAY 1974)                          PAGE 3      
        
        
        
 SUBSECTION: DECSOL   .     
        
        
 CALLING SEQUENCE:
        
     THE HEADING OF THIS PROCEDURE IS:
     "PROCEDURE" DECSOL(A, N, AUX, B); "VALUE" N;         
     "INTEGER" N; "ARRAY" A, AUX, B;  
        
     THE MEANING OF THE FORMAL PARAMETERS IS:   
     A:      <ARRAY IDENTIFIER>;      
             "ARRAY" A[1:N, 1:N];     
             ENTRY:THE N-TH ORDER MATRIX;       
             EXIT: THE  CALCULATED  LOWER-TRIANGULAR  MATRIX  AND  UNIT       
             UPPERTRIANGULAR  MATRIX  WITH ITS  UNIT  DIAGONAL  OMITTED;      
     N:      <ARITHMETIC EXPRESSION>; 
             THE ORDER OF THE MATRIX; 
     AUX:    <ARRAY IDENTIFIER>;      
             "ARRAY" AUX[1:3];        
             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[1]: IF  R  IS THE NUMBER OF ELIMINATION STEPS PERFORMED      
                     (SEE  AUX[3]),   THEN   AUX[1]  EQUALS   1  IF  THE      
                     DETERMINANT OF THE PRINCIPAL SUBMATRIX  OF ORDER  R      
                     IS POSITIVE, ELSE AUX[1] = -1;       
             AUX[3]: THE  NUMBER  OF  ELIMINATION  STEPS  PERFORMED;  IF      
                     AUX[3] < N  THEN  THE PROCESS  IS TERMINATED AND NO      
                     SOLUTION WILL BE CALCULATED;         
     B:      <ARRAY IDENTIFIER>;      
             "ARRAY" B[1:N];
             ENTRY:  THE RIGHT-HAND SIDE OF THE LINEAR SYSTEM;      
             EXIT:   IF AUX[3] = N,  THEN THE CALCULATED SOLUTION OF THE      
                     LINEAR SYSTEM IS OVERWRITTEN ON B,  ELSE  B REMAINS      
                     UNALTERED.       
        
        
 PROCEDURES USED: 
        
     DEC = CP34300,         
     SOL = CP34051.         
        
        
 REQUIRED CENTRAL MEMORY:   
        
     EXECUTION FIELD LENGTH: DECSOL DECLARES AN AUXILIARY ARRAY  OF TYPE      
                             INTEGER AND ORDER N.         
        
 RUNNING TIME: PROPORTIONAL TO N ** 3.
        
 LANGUAGE:   ALGOL 60.      
1SECTION: 3.1.1.1.1.1.3       (MAY 1974)                          PAGE 4      
        
        
        
 METHOD AND PERFORMANCE:    
        
        
     DECSOL  USES  DEC  TO PERFORM  THE TRIANGULAR DECOMPOSITION  OF THE      
     MATRIX  AND  SOL  TO CALCULATE  THE SOLUTION  WITH FORWARD AND BACK      
     SUBSTITUTION;  SINCE  DECSOL  MAY YIELD USELESS RESULTS,  EVEN  FOR      
     WELL-CONDITIONED MATRICES (SEE DEC, SECTION 3.1.1.1.1.1.1),  DECSOL      
     SHOULD ONLY BE USED IF THE ORDER OF THE MATRIX IS SMALL RELATIVE TO      
     THE NUMBER  OF  BINARY  DIGITS  IN  THE  NUMBER REPRESENTATION;  IF      
     AUX[3] < N, THEN THE EFFECT OF DECSOL IS MERELY THAT OF DEC.   
        
        
 EXAMPLE OF USE:  
        
     LET A BE THE FOURTH ORDER SEGMENT OF THE HILBERT MATRIX AND B THE        
     THIRD COLUMN OF A, THEN THE SOLUTION OF THE LINEAR SYSTEM      
     AX = B IS GIVEN BY THE THIRD UNIT VECTOR AND MAY BE CALCULATED BY        
     THE FOLLOWING PROGRAM: 
        
     "BEGIN" "INTEGER" I, J;
         "ARRAY" A[1:4, 1:4], B[1:4], AUX[1:3]; 
         "PROCEDURE" DECSOL(A, N, AUX, B); "CODE" 34301;  
         "PROCEDURE" LIST(ITEM); "PROCEDURE" ITEM;        
         "BEGIN" "INTEGER" I;         
             "FOR" I:= 1 "STEP" 1 "UNTIL" 4 "DO" ITEM(B[I]);        
             "FOR" I:= 1 "STEP" 2 "UNTIL" 3 "DO" ITEM(AUX[I])       
         "END" LIST;        
         "PROCEDURE" LAYOUT;
         FORMAT("("*, "("SOLUTION:")"B+.15D"+3D,/,3(10B+.15D"+3D,/),
         "("SIGN(DET) = ")"+D,/,"("NUMBER OF ELIMINATIONSTEPS = ")" 
         +D")");  
        
         "FOR" I:= 1 "STEP" 1 "UNTIL" 4 "DO"    
         "BEGIN" "FOR" J:= 1 "STEP" 1 "UNTIL" 4 "DO"      
             A[I,J]:= 1 / (I + J - 1); B[I]:= A[I,3]      
         "END";   
         AUX[2]:= "-14;     
         DECSOL(A, 4, AUX, B);        
         OUTLIST(71, LAYOUT, LIST)    
     "END"        
        
     RESULTS:     
        
     SOLUTION: +.000000000000000"+000 
               +.000000000000000"+000 
               +.100000000000000"+001 
               +.000000000000000"+000 
     SIGN(DET) = +1         
     NUMBER OF ELIMINATIONSTEPS = +4  
        
        
        
1SECTION: 3.1.1.1.1.1.3       (MAY 1974)                          PAGE 5      
        
        
 SUBSECTION: SOLELM  .      
        
        
 CALLING SEQUENCE:
        
     THE HEADING OF THIS PROCEDURE IS:
     "PROCEDURE" SOLELM(A, N, RI, CI, B); "VALUE" N;      
     "INTEGER" N; "ARRAY" A, B; "INTEGER" "ARRAY" RI, CI; 
        
     THE MEANING OF THE FORMAL PARAMETERS IS:   
     A:      <ARRAY IDENTIFIER>;      
             "ARRAY" A[1:N, 1:N];     
             ENTRY:  THE TRIANGULARLY  DECOMPOSED FORM  OF THE MATRIX OF      
                     THE LINEAR SYSTEM AS  PRODUCED BY  GSSELM  (SECTION      
                     3.1.1.1.1.1.1);  
     N:      <ARITHMETIC EXPRESSION>; 
             THE ORDER OF THE MATRIX; 
     RI:     <ARRAY IDENTIFIER>;      
             "INTEGER""ARRAY" RI[1:N];
             ENTRY:THE PIVOTAL ROW INDICES,  AS PRODUCED BY  GSSELM;
     CI:     <ARRAY IDENTIFIER>;      
             "INTEGER""ARRAY" CI[1:N];
             ENTRY:THE PIVOTAL COLUMN INDICES, AS PRODUCED BY  GSSELM;        
     B:      <ARRAY IDENTIFIER>;      
             "ARRAY" B[1:N];
             ENTRY:  THE RIGHT-HAND SIDE OF THE LINEAR SYSTEM;      
             EXIT:   THE SOLUTION OF THE LINEAR SYSTEM.   
        
        
 PROCEDURES USED: 
        
     SOL = CP34051.         
        
        
 RUNNING TIME: PROPORTIONAL TO N ** 2.
        
        
 LANGUAGE:   ALGOL 60.      
        
        
 METHOD AND PERFORMANCE:    
        
     SOLELM   SHOULD  BE  CALLED  AFTER   GSSELM   OR   GSSERB  (SECTION      
     3.1.1.1.1.1.1) AND SOLVES THE LINEAR SYSTEM WITH THE MATRIX,  WHOSE      
     TRIANGULARLY  DECOMPOSED  FORM  AS PRODUCED BY  GSSELM  IS GIVEN IN      
     ARRAY A, AND A RIGHT-HAND SIDE AS GIVEN IN ARRAY B;  SOLELM  LEAVES      
     A,  RI AND CI  UNALTERED,  SO, AFTER ONE CALL OF GSSELM  OR GSSERB,      
     SEVERAL CALLS OF  SOLELM  MAY FOLLOW  FOR  SOLVING  SEVERAL SYSTEMS      
     HAVING THE SAME MATRIX  BUT DIFFERENT RIGHT-HAND SIDES.        
        
 EXAMPLE OF USE: SEE GSSSOL OR GSSSOLERB (THIS SECTION).  
1SECTION: 3.1.1.1.1.1.3       (MAY 1974)                          PAGE 6      
        
        
        
 SUBSECTION: GSSSOL   .     
        
        
 CALLING SEQUENCE:
        
     THE HEADING OF THIS PROCEDURE IS:
     "PROCEDURE" GSSSOL(A, N, AUX, B); "VALUE" N;         
     "INTEGER" N; "ARRAY" A, AUX, B;  
        
     THE MEANING OF THE FORMAL PARAMETERS IS:   
     A:      <ARRAY IDENTIFIER>;      
             "ARRAY" A[1:N, 1:N];     
             ENTRY:THE N-TH ORDER MATRIX;       
             EXIT: THE  CALCULATED  LOWER-TRIANGULAR  MATRIX  AND  UNIT       
             UPPERTRIANGULAR  MATRIX  WITH ITS  UNIT  DIAGONAL  OMITTED;      
     N:      <ARITHMETIC EXPRESSION>; 
             THE ORDER OF THE MATRIX; 
     AUX:    <ARRAY IDENTIFIER>;      
             "ARRAY" AUX[1:7];        
             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;     
             AUX[4]: A VALUE WHICH IS USED FOR CONTROLLING PIVOTING (SEE      
                     GSSELM, SECTION 3.1.1.1.1.1.1);      
             EXIT:
             AUX[1]: IF  R  IS THE NUMBER OF ELIMINATION STEPS PERFORMED      
                     (SEE  AUX[3]),   THEN   AUX[1]  EQUALS   1  IF  THE      
                     DETERMINANT OF THE PRINCIPAL SUBMATRIX  OF ORDER  R      
                     IS POSITIVE, ELSE AUX[1] = -1;       
             AUX[3]: THE  NUMBER  OF  ELIMINATION  STEPS  PERFORMED;  IF      
                     AUX[3] < N  THEN  THE PROCESS  IS TERMINATED AND NO      
                     SOLUTION WILL HAVE BEEN CALCULATED;  
             AUX[5]: THE  MODULUS  OF  AN ELEMENT  WHICH  IS  OF MAXIMUM      
                     ABSOLUTE VALUE FOR THE MATRIX GIVEN IN ARRAY A;
             AUX[7]: AN UPPER BOUND FOR THE GROWTH (SEE GSSELM,  SECTION      
                     3.1.1.1.1.1.1);  
     B:      <ARRAY IDENTIFIER>;      
             "ARRAY" B[1:N];
             ENTRY:  THE RIGHT-HAND SIDE OF THE LINEAR SYSTEM;      
             EXIT:   IF AUX[3] = N,  THEN THE CALCULATED SOLUTION OF THE      
                     LINEAR SYSTEM IS OVERWRITTEN ON B,  ELSE  B REMAINS      
                     UNALTERED.       
        
 PROCEDURES USED: 
        
     SOLELM  = CP34061,     
     GSSELM  = CP34231.     
        
 REQUIRED CENTRAL MEMORY:   
        
     EXECUTION FIELD LENGTH: GSSSOL  DECLARES  TWO  AUXILIARY  ARRAYS OF      
                             TYPE INTEGER AND ORDER N.    
1SECTION: 3.1.1.1.1.1.3       (MAY 1974)                          PAGE 7      
        
        
        
 RUNNING TIME: PROPORTIONAL TO N ** 3.
        
        
 LANGUAGE:   ALGOL 60.      
        
        
 METHOD AND PERFORMANCE:    
        
     GSSSOL USES  GSSELM (SECTION 3.1.1.1.1.1.1) TO PERFORM A TRIANGULAR      
     DECOMPOSITION OF THE MATRIX AND  SOLELM (THIS SECTION) TO CALCULATE      
     THE SOLUTION  OF THE GIVEN LINEAR SYSTEM;  IF AUX[3] < N,  THEN THE      
     EFFECT OF GSSSOL IS MERELY THAT OF GSSELM. 
        
        
 EXAMPLE OF USE:  
        
     LET A BE THE FOURTH ORDER SEGMENT OF THE HILBERT MATRIX AND B THE        
     THIRD COLUMN OF A, THEN THE SOLUTION OF THE LINEAR SYSTEM      
     AX = B IS GIVEN BY THE THIRD UNIT VECTOR AND MAY BE CALCULATED BY        
     THE FOLLOWING PROGRAM: 
        
     "BEGIN" "INTEGER" I, J;
         "ARRAY" A[1:4, 1:4], B[1:4], AUX[1:7]; 
         "PROCEDURE" GSSSOL(A, N, AUX, B); "CODE" 34232;  
         "PROCEDURE" LIST(ITEM); "PROCEDURE" ITEM;        
         "BEGIN" "INTEGER" I;         
             "FOR" I:= 1 "STEP" 1 "UNTIL" 4 "DO" ITEM(B[I]);        
             "FOR" I:= 1 "STEP" 2 "UNTIL" 7 "DO" ITEM(AUX[I])       
         "END" LIST;        
         "PROCEDURE" LAYOUT;
         FORMAT("("*, "("SOLUTION:")"B+.15D"+3D,/,3(10B+.15D"+3D,/),
         "("SIGN(DET) = ")"+D,/,"("NUMBER OF ELIMINATIONSTEPS = ")" 
         +D,/,"("MAX(ABS(A[I,J]))= ")"+.15D"+3D,/,        
         "("UPPER BOUND GROWTH: ")"+.15D"+3D")");         
        
         "FOR" I:= 1 "STEP" 1 "UNTIL" 4 "DO"    
         "BEGIN" "FOR" J:= 1 "STEP" 1 "UNTIL" 4 "DO"      
             A[I,J]:= 1 / (I + J - 1); B[I]:= A[I,3]      
         "END";   
         AUX[2]:= "-14; AUX[4]:= 8;   
         GSSSOL(A, 4, AUX, B);        
         OUTLIST(71, LAYOUT, LIST)    
     "END"        
        
     RESULTS:     
        
     SOLUTION: +.888178419700120"-014 
               -.497379915032070"-013 
               +.100000000000010"+001 
               +.000000000000000"+000 
     SIGN(DET) = +1         
     NUMBER OF ELIMINATIONSTEPS = +4  
     MAX(ABS(A[I,J]))= +.100000000000000"+001   
     UPPER BOUND GROWTH: +.159619047619050"+001 
1SECTION: 3.1.1.1.1.1.3       (MAY 1974)                          PAGE 8      
        
        
        
 SUBSECTION: GSSSOLERB.     
        
 CALLING SEQUENCE:
     THE HEADING OF THIS PROCEDURE IS:
     "PROCEDURE" GSSSOLERB(A, N, AUX, B); "VALUE" N;      
     "INTEGER" N; "ARRAY" A, AUX, B;  
        
     THE MEANING OF THE FORMAL PARAMETERS IS:   
     A:      <ARRAY IDENTIFIER>;      
             "ARRAY" A[1:N, 1:N];     
             ENTRY:THE N-TH ORDER MATRIX;       
             EXIT: THE  CALCULATED  LOWER-TRIANGULAR  MATRIX  AND  UNIT       
             UPPERTRIANGULAR  MATRIX  WITH ITS  UNIT  DIAGONAL  OMITTED;      
     N:      <ARITHMETIC EXPRESSION>; 
             THE ORDER OF THE MATRIX; 
     AUX:    <ARRAY IDENTIFIER>;      
             "ARRAY" AUX[0:11];       
             ENTRY:         
             AUX[0]: THE MACHINE PRECISION;     
             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;     
             AUX[4]: A VALUE WHICH IS USED FOR CONTROLLING PIVOTING (SEE      
                     GSSELM, SECTION 3.1.1.1.1.1.1);      
             AUX[6]: AN UPPER BOUND  FOR  THE RELATIVE PRECISION  OF THE      
                     GIVEN MATRIX ELEMENTS;     
             EXIT:
             AUX[1]: IF  R  IS THE NUMBER OF ELIMINATION STEPS PERFORMED      
                     (SEE  AUX[3]),   THEN   AUX[1]  EQUALS   1  IF  THE      
                     DETERMINANT OF THE PRINCIPAL SUBMATRIX  OF ORDER  R      
                     IS POSITIVE, ELSE AUX[1] = -1;       
             AUX[3]: THE  NUMBER  OF  ELIMINATION  STEPS  PERFORMED;  IF      
                     AUX[3] < N  THEN  THE PROCESS  IS TERMINATED AND NO      
                     SOLUTION OR ERROR BOUND WILL HAVE BEEN CALCULATED;       
             AUX[5]: THE  MODULUS  OF  AN ELEMENT  WHICH  IS  OF MAXIMUM      
                     ABSOLUTE VALUE FOR THE MATRIX GIVEN IN ARRAY A;
             AUX[7]: AN UPPER BOUND FOR THE GROWTH (SEE GSSELM,  SECTION      
                     3.1.1.1.1.1.1);  
             AUX[9]: IF AUX[3] = N, THEN AUX[9] WILL EQUAL THE 1-NORM OF      
                     THE INVERSE MATRIX, ELSE AUX[9] WILL BE UNDEFINED;       
             AUX[11]: IF AUX[3] = N THEN THE VALUE OF AUX[11]  WILL BE A      
                     ROUGH  UPPER BOUND  FOR  THE RELATIVE ERROR  IN THE      
                     CALCULATED  SOLUTION  OF  THE  GIVEN LINEAR SYSTEM,      
                     ELSE  AUX[11]  WILL BE UNDEFINED;  IF NO USE CAN BE      
                     MADE OF THE FORMULA FOR THE ERROR BOUND AS GIVEN IN      
                     SECTION 3.1.1.1.1.1.1 (SUBSECTION ERBELM),  BECAUSE      
                     OF  A  VERY  BAD  CONDITION  OF  THE  MATRIX,  THEN      
                     AUX[11]:= -1;    
     B:      <ARRAY IDENTIFIER>;      
             "ARRAY" B[1:N];
             ENTRY:  THE RIGHT-HAND SIDE OF THE LINEAR SYSTEM;      
             EXIT:   IF AUX[3] = N,  THEN THE CALCULATED SOLUTION OF THE      
                     LINEAR SYSTEM, ELSE  B REMAINS UNALTERED.      
1SECTION: 3.1.1.1.1.1.3       (MAY 1974)                          PAGE 9      
        
        
        
 PROCEDURES USED: 
        
     SOLELM  = CP34061,     
     GSSERB  = CP34242.     
        
 REQUIRED CENTRAL MEMORY:   
        
     EXECUTION FIELD LENGTH: GSSSOLERB  DECLARES TWO AUXILIARY ARRAYS OF      
                             TYPE INTEGER AND ORDER N.    
        
 RUNNING TIME: PROPORTIONAL TO N ** 3.
        
 LANGUAGE:   ALGOL 60.      
        
 METHOD AND PERFORMANCE:    
        
     GSSSOLERB  USES  GSSERB  (SECTION  3.1.1.1.1.1.1)  TO  PERFORM  THE      
     TRIANGULAR DECOMPOSITION OF THE MATRIX  AND  TO CALCULATE  AN UPPER      
     BOUND  FOR THE RELATIVE ERROR  IN  THE CALCULATED SOLUTION  OF  THE      
     GIVEN LINEAR SYSTEM, AND  SOLELM  (THIS SECTION)  TO CALCULATE THIS      
     SOLUTION;  IF AUX[3] < N,  THEN THE EFFECT OF  GSSSOLERB  IS MERELY      
     THAT OF GSSELM (SECTION 3.1.1.1.1.1.1).    
        
 EXAMPLE OF USE:  
     LET A BE THE FOURTH ORDER SEGMENT OF THE HILBERT MATRIX AND B THE        
     THIRD COLUMN OF A, THEN THE SOLUTION OF THE LINEAR SYSTEM      
     AX = B IS GIVEN BY THE THIRD UNIT VECTOR AND THIS SOLUTION, AS 
     WELL AS AN UPPER BOUND FOR THE RELATIVE ERROR IN THE CALCULATED ONE      
     , MAY BE OBTAINED BY THE FOLLOWING PROGRAM:
        
     "BEGIN" "INTEGER" I, J;
         "ARRAY" A[1:4, 1:4], B[1:4], AUX[0:11];
         "PROCEDURE" GSSSOLERB(A, N, AUX, B); "CODE" 34243;         
         "PROCEDURE" LIST(ITEM); "PROCEDURE" ITEM;        
         "BEGIN" "INTEGER" I;         
             "FOR" I:= 1 "STEP" 1 "UNTIL" 4 "DO" ITEM(B[I]);        
             "FOR" I:= 1 "STEP" 2 "UNTIL" 11 "DO" ITEM(AUX[I])      
         "END" LIST;        
         "PROCEDURE" LAYOUT;
         FORMAT("("*, "("SOLUTION:")"B+.15D"+3D,/,3(10B+.15D"+3D,/),
         "("SIGN(DET) = ")"+D,/,"("NUMBER OF ELIMINATIONSTEPS = ")" 
         +D,/,"("MAX(ABS(A[I,J]))= ")"+.15D"+3D,/,        
         "("UPPER BOUND GROWTH: ")"+.15D"+3D,/, 
         "("1-NORM OF THE INVERSE MATRIX:")"B+.15D"+3D,/, 
         "("UPPER BOUND REL. ERR. IN THE CALC. SOL.")"    
         B+.15D"+3D")");    
         "FOR" I:= 1 "STEP" 1 "UNTIL" 4 "DO"    
         "BEGIN" "FOR" J:= 1 "STEP" 1 "UNTIL" 4 "DO"      
             A[I,J]:= 1 / (I + J - 1); B[I]:= A[I,3]      
         "END";   
         AUX[0]:= AUX[2]:= "-14; AUX[4]:= 8; AUX[6]:= "-14;         
         GSSSOLERB(A, 4, AUX, B);     
         OUTLIST(71, LAYOUT, LIST)    
     "END"        
1SECTION: 3.1.1.1.1.1.3       (MAY 1974)                         PAGE 10      
        
        
        
     RESULTS:     
        
     SOLUTION: +.888178419700120"-014 
               -.497379915032070"-013 
               +.100000000000010"+001 
               +.000000000000000"+000 
     SIGN(DET) = +1         
     NUMBER OF ELIMINATIONSTEPS = +4  
     MAX(ABS(A[I,J]))= +.100000000000000"+001   
     UPPER BOUND GROWTH: +.159619047619050"+001 
     1-NORM OF THE INVERSE MATRIX: +.136199999998790"+005 
     UPPER BOUND REL. ERR. IN THE CALC. SOL. +.277896269157090"-007 
        
        
        
 REFERENCES:      
        
     [1] BUS, J. C. P.      
         LINEAR SYSTEMS  WITH CALCULATION OF ERROR BOUNDS  AND ITERATIVE      
         REFINEMENT (DUTCH).
         MATHEMATICAL CENTRE, AMSTERDAM, LR 3. 4. 19 (1972).        
     [2] DEKKER, T. J.      
         ALGOL 60 PROCEDURES IN NUMERICAL ALGEBRA, PART 1.
         MATHEMATICAL CENTRE, AMSTERDAM, TRACT 22 (1968). 
        
 SOURCE TEXT(S):  
        
 "CODE" 34051;    
     "PROCEDURE" SOL(A, N, P, B); "VALUE" N; "INTEGER" N; "ARRAY" A, B;       
     "INTEGER" "ARRAY" P;   
     "BEGIN" "INTEGER" K, PK;         
         "REAL" R;
         "REAL" "PROCEDURE" MATVEC(L, U, I, A, B); "CODE" 34011;    
         "FOR" K:= 1 "STEP" 1 "UNTIL" N "DO"    
         "BEGIN" R:= B[K]; PK:= P[K]; 
             B[K]:= (B[PK] - MATVEC(1, K - 1, K, A, B)) / A[K,K];   
             "IF" PK ^= K "THEN" B[PK]:= R      
         "END";   
         "FOR" K:= N "STEP" - 1 "UNTIL" 1 "DO"  
         B[K]:= B[K] - MATVEC(K + 1, N, K, A, B)
     "END" SOL;   
         "EOP"    
1SECTION: 3.1.1.1.1.1.3       (MAY 1974)                         PAGE 11      
        
        
        
0"CODE" 34301;    
     "PROCEDURE" DECSOL(A, N, AUX, B); "VALUE" N; "INTEGER" N;      
     "ARRAY" A, AUX, B;     
     "BEGIN" "INTEGER" "ARRAY" P[1:N];
         "PROCEDURE" SOL(A, N, P, B); "CODE" 34051;       
         "PROCEDURE" DEC(A, N, AUX, P); "CODE" 34300;     
         DEC(A, N, AUX, P); 
         "IF" AUX[3] = N "THEN" SOL(A, N, P, B) 
     "END" DECSOL;
         "EOP"    
0"CODE" 34061;    
     "PROCEDURE" SOLELM(A, N, RI, CI, B); "VALUE" N; "INTEGER" N;   
     "ARRAY" A, B;
     "INTEGER" "ARRAY" RI, CI;        
     "BEGIN" "INTEGER" R, CIR;        
         "REAL" W;
         "PROCEDURE" SOL(A, N, P, B); "CODE" 34051;       
         SOL(A, N, RI, B);  
         "FOR" R:= N "STEP" - 1 "UNTIL" 1 "DO"  
         "BEGIN" CIR:= CI[R]; "IF" CIR ^= R "THEN"        
             "BEGIN" W:= B[R]; B[R]:= B[CIR]; B[CIR]:= W "END"      
         "END"    
     "END" SOLELM;
         "EOP"    
 "CODE" 34232;    
     "PROCEDURE" GSSSOL(A, N, AUX, B); "VALUE" N; "INTEGER" N;      
     "ARRAY" A, AUX, B;     
     "BEGIN" "INTEGER" "ARRAY" RI, CI[1:N];     
         "PROCEDURE" SOLELM(A, N, RI, CI, B); "CODE" 34061;         
         "PROCEDURE" GSSELM(A, N, AUX, RI, CI); "CODE" 34231;       
         GSSELM(A, N, AUX, RI, CI);   
         "IF" AUX[3] = N "THEN" SOLELM(A, N, RI, CI, B)   
     "END" GSSSOL;
         "EOP"    
0"CODE" 34243;    
     "PROCEDURE" GSSSOLERB(A, N, AUX, B); "VALUE" N; "INTEGER" N;   
     "ARRAY" A, AUX, B;     
     "BEGIN" "INTEGER" "ARRAY" RI, CI[1:N];     
         "PROCEDURE" SOLELM(A, N, RI, CI, B); "CODE" 34061;         
         "PROCEDURE" GSSERB(A, N, AUX, RI, CI); "CODE" 34242;       
         GSSERB(A, N, AUX, RI, CI);   
         "IF" AUX[3] = N "THEN" SOLELM(A, N, RI, CI, B)   
     "END" GSSSOLERB;       
         "EOP"    
1SECTION: 3.1.1.1.1.1.4       (DECEMBER 1975)                     PAGE 1      
        
        
        
 AUTHORS: J. C. P. BUS AND T. J. DEKKER.        
        
        
 CONTRIBUTOR: J.C.P. BUS AND P.A. BEENTJES.     
        
        
 INSTITUTE: MATHEMATICAL CENTRE.      
        
        
 RECEIVED: 730920.
        
        
 BRIEF DESCRIPTION:         
        
     THIS SECTION CONTAINS  FIVE  PROCEDURES  FOR INVERSION OF MATRICES:      
     INV CALCULATES THE INVERSE OF A MATRIX THAT HAS BEEN TRIANGULARLY        
     DECOMPOSED BY DEC;     
     DECINV CALCULATES THE INVERSE OF A MATRIX WHOSE ORDER IS SMALL 
     RELATIVE  TO  THE  NUMBER  OF  BINARY  DIGITS  IN  THE  NUMBER 
     REPRESENTATION;        
     INV1 CALCULATES THE INVERSE OF A MATRIX THAT HAS BEEN TRIANGULARLY       
     DECOMPOSED BY GSSELM OR GSSERB.THE 1-NORM OF THE INVERSE MATRIX
     MIGHT ALSO BE CALCULATED         
     GSSINV CALCULATES THE INVERSE OF A MATRIX; 
     GSSINVERB CALCULATES THE INVERSE OF A MATRIX AND ITS 1-NORM.   
     A ROUGH UPPERBOUND FOR THE RELATIVE ERROR IN THE CALCULATED INVERSE      
     MATRIX IS ALSO GIVEN;  
        
     THE DIFFERENCE         
     BETWEEN  DECINV  ON THE ONE SIDE AND  GSSINV  AND  GSSINVERB ON THE      
     OTHER SIDE  LIES IN THE METHOD  USED  FOR TRIANGULAR DECOMPOSITION,      
     PARTICULARLY IN THE PIVOTING STRATEGY;  DECINV USES DEC, GSSINV AND      
     GSSINVERB  USE  GSSELM TO PERFORM THE TRIANGULAR DECOMPOSITION; THE      
     USER  IS  ADVISED  TO  USE   GSSINV  OR   GSSINVERB   (SEE  SECTION      
     3.1.1.1.1.1.1).        
        
        
 KEYWORDS:        
        
     MATRIX INVERSION.      
        
        
        
1SECTION: 3.1.1.1.1.1.4       (MAY 1974)                          PAGE 2      
        
        
        
 SUBSECTION: INV .
        
        
 CALLING SEQUENCE:
        
     THE HEADING OF THIS PROCEDURE IS:
     "PROCEDURE" INV(A, N, P); "VALUE" N;       
     "INTEGER" N; "ARRAY" A; "INTEGER" "ARRAY" P;         
        
     THE MEANING OF THE FORMAL PARAMETERS IS:   
     A:      <ARRAY IDENTIFIER>;      
             "ARRAY" A[1:N, 1:N];     
             ENTRY:  THE TRIANGULARLY  DECOMPOSED FORM  OF THE MATRIX AS      
                     PRODUCED BY  DEC  (SECTION 3.1.1.1.1.1.1);     
             EXIT:   THE CALCULATED INVERSE MATRIX;       
     N:      <ARITHMETIC EXPRESSION>; 
             THE ORDER OF THE MATRIX; 
     P:      <ARRAY IDENTIFIER>;      
             "INTEGER""ARRAY P[1:N];  
             ENTRY:THE PIVOTAL INDICES, AS PRODUCED BY DEC;         
        
        
 PROCEDURES USED: 
        
     MATMAT     = CP34013,  
     ICHCOL     = CP34031,  
     DUPCOLVEC  = CP31034.  
        
        
 REQUIRED CENTRAL MEMORY:   
        
     EXECUTION FIELD LENGTH: INV  DECLARES  AN AUXILIARY ARRAY  OF  TYPE      
                             REAL AND ORDER N.  
        
        
 RUNNING TIME: PROPORTIONAL TO N ** 3.
        
        
 LANGUAGE:   ALGOL 60.      
        
        
 METHOD AND PERFORMANCE:    
        
     INV  SHOULD  BE  CALLED  AFTER  DEC  (SECTION  3.1.1.1.1.1.1)   AND      
     CALCULATES THE INVERSE OF THE MATRIX, WHOSE TRIANGULARLY DECOMPOSED      
     FORM AS PRODUCED BY  DEC IS GIVEN IN ARRAY A; THE INVERSE MATRIX IS      
     OVERWRITTEN ON A.      
        
        
 EXAMPLE OF USE: SEE DECINV (THIS SECTION).     
        
        
        
1SECTION: 3.1.1.1.1.1.4       (DECEMBER 1975)                     PAGE 3      
        
        
        
 SUBSECTION: DECINV   .     
        
        
 CALLING SEQUENCE:
        
     THE HEADING OF THIS PROCEDURE IS:
     "PROCEDURE" DECINV(A, N, AUX); "VALUE" N;  
     "INTEGER" N; "ARRAY" A, AUX;     
        
     THE MEANING OF THE FORMAL PARAMETERS IS:   
     A:      <ARRAY IDENTIFIER>;      
             "ARRAY" A[1:N, 1:N];     
             ENTRY:  THE MATRIX, WHOSE INVERSE HAS TO BE CALCULATED;
             EXIT:   IF  AUX[3] = N, THEN THE CALCULATED INVERSE MATRIX;      
     N:      <ARITHMETIC EXPRESSION>; 
             THE ORDER OF THE MATRIX; 
     AUX:    <ARRAY IDENTIFIER>;      
             "ARRAY" AUX[1:3];        
             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[1]: IF  R  IS THE NUMBER OF ELIMINATION STEPS PERFORMED      
                     (SEE  AUX[3]),   THEN   AUX[1]  EQUALS   1  IF  THE      
                     DETERMINANT OF THE PRINCIPAL SUBMATRIX  OF ORDER  R      
                     IS POSITIVE, ELSE AUX[1] = -1;       
             AUX[3]: THE  NUMBER  OF  ELIMINATION  STEPS  PERFORMED;  IF      
                     AUX[3] < N, THEN  THE PROCESS  IS TERMINATED AND NO      
                     INVERSE WILL HAVE BEEN CALCULATED.   
        
        
 PROCEDURES USED: 
        
     DEC = CP34300,         
     INV = CP34053.         
        
        
 REQUIRED CENTRAL MEMORY:   
        
     EXECUTION FIELD LENGTH: DECINV DECLARES AN AUXILIARY ARRAY  OF TYPE      
                             INTEGER AND ORDER N.         
        
        
 RUNNING TIME: PROPORTIONAL TO N ** 3.
        
        
 LANGUAGE:   ALGOL 60.      
        
        
1SECTION: 3.1.1.1.1.1.4       (MAY 1974)                          PAGE 4      
        
        
        
 METHOD AND PERFORMANCE:    
        
     DECINV  USES  DEC (SECTION 3.1.1.1.1.1.1) TO PERFORM THE TRIANGULAR      
     DECOMPOSITION OF A MATRIX AND INV TO CALCULATE ITS INVERSE;  DECINV      
     SHOULD ONLY BE USED IF THE ORDER OF THE MATRIX IS SMALL RELATIVE TO      
     THE NUMBER OF BINARY DIGITS IN THE NUMBER REPRESENTATION (SEE DEC);      
     IF AUX[3] < N, THEN THE EFFECT OF DECINV IS MERELY THAT OF DEC.
        
        
 EXAMPLE OF USE:  
        
     THE FOLLOWING PROGRAM CALCULATES  THE INVERSE  OF  THE INPUT MATRIX      
     AND PRINTS THE RESULTS:
        
     "BEGIN"      
         "ARRAY" A[1:4, 1:4], AUX[1:3];         
         "PROCEDURE" DECINV(A, N, AUX); "CODE" 34302;     
         "PROCEDURE" LIST(ITEM); "PROCEDURE" ITEM;        
         "BEGIN" "INTEGER" I, J;      
             "FOR" I:= 1 "STEP" 1 "UNTIL" 4 "DO"
             "FOR" J:= 1 "STEP" 1 "UNTIL" 4 "DO" ITEM(A[I,J])       
         "END" LIST;        
         "PROCEDURE" LAYOUT;
         FORMAT("("4(4B,4(B+ZDB),/),/")");      
        
         INLIST(70, LAYOUT, LIST); AUX[2]:= "-14;         
         OUTPUT(71,"("/,"("CALCULATED INVERSE:")",/")");  
         DECINV(A, 4, AUX); 
         OUTLIST(71, LAYOUT, LIST);   
         OUTPUT(71, "(""("AUX[1]=")"B+D,/,"("AUX[3]=")"B+D")",      
         AUX[1], AUX[3])    
     "END"        
        
     INPUT:       
        
      + 4  + 2  + 4  + 1    
      +30  +20  +45  +12    
      +20  +15  +36  +10    
      +35  +28  +70  +20    
        
        
     RESULTS:     
        
     CALCULATED INVERSE:    
       +4   -2   +4   -1    
      -30  +20  -45  +12    
      +20  -15  +36  -10    
      -35  +28  -70  +20    
        
     AUX[1]= +1   
     AUX[3]= +4   
        
        
        
1SECTION: 3.1.1.1.1.1.4       (MAY 1974)                          PAGE 5      
        
        
        
 SUBSECTION: INV1  .        
        
        
 CALLING SEQUENCE:
        
     THE HEADING OF THIS PROCEDURE IS:
     "REAL" "PROCEDURE" INV1(A, N, RI, CI, WITHNORM);     
     "VALUE" N, WITHNORM; "INTEGER" N; "BOOLEAN" WITHNORM; "ARRAY" A;         
     "INTEGER" "ARRAY" RI, CI;        
        
     INV1:   IF THE VALUE OF  WITHNORM  IS TRUE, THEN THE VALUE OF  INV1      
             WILL EQUAL THE  1-NORM  OF THE CALCULATED  INVERSE  MATRIX,      
             ELSE INV1:= 0; 
        
     THE MEANING OF THE FORMAL PARAMETERS IS:   
     A:      <ARRAY IDENTIFIER>;      
             "ARRAY" A[1:N, 1:N];     
             ENTRY:  THE TRIANGULARLY  DECOMPOSED FORM  OF THE MATRIX AS      
                     PRODUCED BY  GSSELM  (SECTION 3.1.1.1.1.1.1);  
             EXIT:   THE CALCULATED INVERSE MATRIX;       
     N:      <ARITHMETIC EXPRESSION>; 
             THE ORDER OF THE MATRIX; 
     RI:     <ARRAY IDENTIFIER>;      
             "INTEGER""ARRAY" RI[1:N];
             ENTRY:THE PIVOTAL ROW INDICES,  AS PRODUCED BY  GSSELM;
     CI:     <ARRAY IDENTIFIER>;      
             "INTEGER""ARRAY" CI[1:N];
             ENTRY:THE PIVOTAL COLUMN INDICES, AS PRODUCED BY  GSSELM;        
     WITHNORM: <BOOLEAN EXPRESSION>;  
             IF THE VALUE OF  WITHNORM  IS TRUE,  THEN THE 1-NORM OF THE      
             INVERSE MATRIX  WILL  BE CALCULATED  AND  ASSIGNED TO INV1,      
             ELSE INV1:= 0; 
        
 PROCEDURES USED: 
        
     ICHROW = CP34032,      
     INV    = CP34053.      
        
 RUNNING TIME: PROPORTIONAL TO N ** 3.
        
 LANGUAGE:   ALGOL 60.      
        
        
 METHOD AND PERFORMANCE:    
        
     INV1   SHOULD  BE  CALLED  AFTER   GSSELM   OR   GSSERB    (SECTION      
     3.1.1.1.1.1.1),WHICH DELIVERS THE TRIANGULARLY DECOMPOSED FORM OF A      
     MATRIX;  INV1  CALCULATES  THE INVERSE MATRIX  AND  ALSO ITS 1-NORM      
     MIGHT BE CALCULATED; THE INVERSE MATRIX IS OVERWRITTEN ON A.   
        
 EXAMPLE OF USE: SEE GSSINV AND GSSINVERB (THIS SECTION). 
1SECTION: 3.1.1.1.1.1.4       (DECEMBER 1975)                     PAGE 6      
        
        
        
 SUBSECTION: GSSINV   .     
        
 CALLING SEQUENCE:
        
     THE HEADING OF THIS PROCEDURE IS:
     "PROCEDURE" GSSINV(A, N, AUX); "VALUE" N;  
     "INTEGER" N; "ARRAY" A, AUX;     
        
     THE MEANING OF THE FORMAL PARAMETERS IS:   
     A:      <ARRAY IDENTIFIER>;      
             "ARRAY" A[1:N, 1:N];     
             ENTRY:  THE MATRIX, WHOSE INVERSE HAS TO BE CALCULATED;
             EXIT:   IF  AUX[3] = N, THEN THE CALCULATED INVERSE MATRIX;      
     N:      <ARITHMETIC EXPRESSION>; 
             THE ORDER OF THE MATRIX; 
     AUX:    <ARRAY IDENTIFIER>;      
             "ARRAY" AUX[1:9];        
             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;     
             AUX[4]: A VALUE WHICH IS USED FOR CONTROLLING PIVOTING (SEE      
                     GSSELM, SECTION 3.1.1.1.1.1.1);      
             EXIT:
             AUX[1]: IF  R  IS THE NUMBER OF ELIMINATION STEPS PERFORMED      
                     (SEE  AUX[3]),   THEN   AUX[1]  EQUALS   1  IF  THE      
                     DETERMINANT OF THE PRINCIPAL SUBMATRIX  OF ORDER  R      
                     IS POSITIVE, ELSE AUX[1] = -1;       
             AUX[3]: THE  NUMBER  OF  ELIMINATION  STEPS  PERFORMED;  IF      
                     AUX[3] < N  THEN  THE PROCESS  IS TERMINATED AND NO      
                     SOLUTION WILL BE CALCULATED;         
             AUX[5]: THE  MODULUS  OF  AN ELEMENT  WHICH  IS  OF MAXIMUM      
                     ABSOLUTE VALUE FOR THE MATRIX GIVEN IN ARRAY A;
             AUX[7]: AN UPPER BOUND FOR THE GROWTH (SEE GSSELM,  SECTION      
                     3.1.1.1.1.1.1);  
             AUX[9]: IF AUX[3] = N, THEN AUX[9] WILL EQUAL THE 1-NORM OF      
                     THE CALCULATED INVERSE MATRIX, ELSE AUX[9]  WILL BE      
                     UNDEFINED.       
        
 PROCEDURES USED: 
        
     INV1    = CP34235,     
     GSSELM  = CP34231.     
        
 REQUIRED CENTRAL MEMORY:   
        
     EXECUTION FIELD LENGTH: GSSINV  DECLARES  TWO  AUXILIARY  ARRAYS OF      
                             TYPE INTEGER AND ORDER N.    
        
 RUNNING TIME: PROPORTIONAL TO N ** 3.
        
 LANGUAGE:   ALGOL 60.      
1SECTION: 3.1.1.1.1.1.4       (MAY 1974)                          PAGE 7      
        
        
        
 METHOD AND PERFORMANCE:    
        
     GSSINV USES  GSSELM (SECTION 3.1.1.1.1.1.1) TO PERFORM A TRIANGULAR      
     DECOMPOSITION OF THE MATRIX AND  INV1 (THIS  SECTION)  TO CALCULATE      
     THE INVERSE MATRIX;  IF AUX[3] = N,  THEN THE EFFECT OF  GSSINV  IS      
     MERELY THAT OF GSSELM. 
        
        
 EXAMPLE OF USE:  
        
     THE FOLLOWING PROGRAM CALCULATES  THE INVERSE  OF  THE INPUT MATRIX      
     AND PRINTS THE RESULTS:
        
     "BEGIN"      
         "ARRAY" A[1:4, 1:4], AUX[1:9];         
         "PROCEDURE" GSSINV(A, N, AUX); "CODE" 34236;     
         "PROCEDURE" LIST(ITEM); "PROCEDURE" ITEM;        
         "BEGIN" "INTEGER" I, J;      
             "FOR" I:= 1 "STEP" 1 "UNTIL" 4 "DO"
             "FOR" J:= 1 "STEP" 1 "UNTIL" 4 "DO" ITEM(A[I,J])       
         "END" LIST;        
         "PROCEDURE" LAYOUT;
         FORMAT("("4(4B,4(B+ZDB),/),/")");      
        
         INLIST(70, LAYOUT, LIST); AUX[2]:= "-14; AUX[4]:= 8;       
         GSSINV(A, 4, AUX); 
         OUTPUT(71,"("/,"("CALCULATED INVERSE:")",/")");  
         OUTLIST(71, LAYOUT, LIST);   
         OUTPUT(71, "("4B"("AUX ELEMENTS:")",/,2(4B+D,/), 
         3(4B+.15D"+3D,/)")", AUX[1], AUX[3], AUX[5], AUX[7], AUX[9])         
     "END"        
        
     INPUT:       
        
      + 4  + 2  + 4  + 1    
      +30  +20  +45  +12    
      +20  +15  +36  +10    
      +35  +28  +70  +20    
        
        
     RESULTS:     
        
     CALCULATED INVERSE:    
       +4   -2   +4   -1    
      -30  +20  -45  +12    
      +20  -15  +36  -10    
      -35  +28  -70  +20    
        
     AUX ELEMENTS:
     +1 
     +4 
     +.700000000000000"+002 
     +.112528571428570"+003 
     +.154999999999730"+003 
1SECTION: 3.1.1.1.1.1.4       (DECEMBER 1975)                     PAGE 8      
        
        
        
 SUBSECTION: GSSINVERB.     
        
        
 CALLING SEQUENCE:
        
     THE HEADING OF THIS PROCEDURE IS:
     "PROCEDURE" GSSINVERB(A, N, AUX); "VALUE" N;         
     "INTEGER" N; "ARRAY" A, AUX;     
        
     THE MEANING OF THE FORMAL PARAMETERS IS:   
     A:      <ARRAY IDENTIFIER>;      
             "ARRAY" A[1:N, 1:N];     
             ENTRY:  THE MATRIX, WHOSE INVERSE HAS TO BE CALCULATED;
             EXIT:   IF AUX[3] = N,  THEN THE CALCULATED INVERSE MATRIX;      
     N:      <ARITHMETIC EXPRESSION>; 
             THE ORDER OF THE MATRIX; 
     AUX:    <ARRAY IDENTIFIER>;      
             "ARRAY" AUX[0:11];       
             ENTRY:         
             AUX[0]: THE MACHINE PRECISION;     
             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;     
             AUX[4]: A VALUE WHICH IS USED FOR CONTROLLING PIVOTING (SEE      
                     GSSELM, SECTION 3.1.1.1.1.1.1);      
             AUX[6]: AN UPPER BOUND  FOR  THE RELATIVE PRECISION  OF THE      
                     GIVEN MATRIX ELEMENTS;     
             EXIT:
             AUX[1]: IF  R  IS THE NUMBER OF ELIMINATION STEPS PERFORMED      
                     (SEE  AUX[3]),   THEN   AUX[1]  EQUALS   1  IF  THE      
                     DETERMINANT OF THE PRINCIPAL SUBMATRIX  OF ORDER  R      
                     IS POSITIVE, ELSE AUX[1] = -1;       
             AUX[3]: THE  NUMBER  OF  ELIMINATION  STEPS  PERFORMED;  IF      
                     AUX[3] < N  THEN  THE PROCESS  IS TERMINATED AND NO      
                     SOLUTION WILL HAVE BEEN CALCULATED;  
             AUX[5]: THE  MODULUS  OF  AN ELEMENT  WHICH  IS  OF MAXIMUM      
                     ABSOLUTE VALUE FOR THE MATRIX GIVEN IN ARRAY A;
             AUX[7]: AN UPPER BOUND FOR THE GROWTH (SEE GSSELM,  SECTION      
                     3.1.1.1.1.1.1);  
             AUX[9]: IF AUX[3] = N, THEN AUX[9] WILL EQUAL THE 1-NORM OF      
                     THE INVERSE MATRIX, ELSE AUX[9] WILL BE UNDEFINED;       
             AUX[11]: IF AUX[3] = N THEN THE VALUE OF AUX[11]  WILL BE A      
                     ROUGH  UPPER BOUND  FOR  THE RELATIVE ERROR  IN THE      
                     CALCULATED  INVERSE MATRIX,  ELSE  AUX[11]  WILL BE      
                     BE UNDEFINED;  IF NO USE CAN BE MADE OF THE FORMULA      
                     FOR  THE  ERROR  BOUND    AS   GIVEN   IN   SECTION      
                     3.1.1.1.1.1.1  (SUBSECTION  ERBELM),   BECAUSE OF A      
                     VERY BAD CONDITION OF THE MATRIX, THEN AUX[11]:=-1.      
        
        
1SECTION: 3.1.1.1.1.1.4       (MAY 1974)                          PAGE 9      
        
        
        
 PROCEDURES USED: 
        
     INV1    = CP34235,     
     GSSELM  = CP34231,     
     ERBELM  = CP34241.     
        
 REQUIRED CENTRAL MEMORY:   
        
     EXECUTION FIELD LENGTH: GSSINVERB  DECLARES TWO AUXILIARY ARRAYS OF      
                             TYPE INTEGER AND ORDER N.    
        
        
 RUNNING TIME: PROPORTIONAL TO N ** 3.
        
        
 LANGUAGE:   ALGOL 60.      
        
        
 METHOD AND PERFORMANCE:    
        
     GSSINVERB  USES  GSSELM  (SECTION  3.1.1.1.1.1.1)  TO  PERFORM  THE      
     TRIANGULAR  DECOMPOSITION  OF THE MATRIX,  INV1  (THIS SECTION)  TO      
     CALCULATE THE INVERSE MATRIX AND ITS  1-NORM  AND  ERBELM  (SECTION      
     3.1.1.1.1.1.1) TO CALCULATE AN UPPER BOUND  FOR  THE RELATIVE ERROR      
     IN  THE  CALCULATED  INVERSE;  IF  AUX[3] < N,  THEN  THE EFFECT OF      
     GSSINVERB IS MERELY THAT OF GSSELM.        
        
        
 EXAMPLE OF USE:  
        
     THE FOLLOWING PROGRAM CALCULATES  THE INVERSE  OF  THE INPUT MATRIX      
     WITH AN UPPER BOUND  FOR THE RELATIVE ERROR IN IT  AND  PRINTS  THE      
     RESULTS:     
        
     "BEGIN"      
         "ARRAY" A[1:4, 1:4], AUX[0:11];        
         "PROCEDURE" GSSINVERB(A, N, AUX); "CODE" 34244;  
         "PROCEDURE" LIST(ITEM); "PROCEDURE" ITEM;        
         "BEGIN" "INTEGER" I, J;      
             "FOR" I:= 1 "STEP" 1 "UNTIL" 4 "DO"
             "FOR" J:= 1 "STEP" 1 "UNTIL" 4 "DO" ITEM(A[I,J])       
         "END" LIST;        
         "PROCEDURE" LAYOUT;
         FORMAT("("4(4B,4(B+ZDB),/),/")");      
        
         INLIST(70, LAYOUT, LIST); AUX[0]:= AUX[2]:= AUX[6]:= "-14; 
         AUX[4]:= 8; GSSINVERB(A, 4, AUX);      
         OUTPUT(71,"("/,"("CALCULATED INVERSE:")",/")");  
         OUTLIST(71, LAYOUT, LIST);   
         OUTPUT(71, "("4B"("AUX ELEMENTS:")",/,2(4B+D,/), 
         4(4B+.15D"+3D,/)")", AUX[1], AUX[3], AUX[5], AUX[7], AUX[9],         
         AUX[11]) 
     "END"        
        
1SECTION: 3.1.1.1.1.1.4       (MAY 1974)                         PAGE 10      
        
        
        
     INPUT:       
        
      + 4  + 2  + 4  + 1    
      +30  +20  +45  +12    
      +20  +15  +36  +10    
      +35  +28  +70  +20    
        
        
     RESULTS:     
        
     CALCULATED INVERSE:    
       +4   -2   +4   -1    
      -30  +20  -45  +12    
      +20  -15  +36  -10    
      -35  +28  -70  +20    
        
     AUX ELEMENTS:
     +1 
     +4 
     +.700000000000000"+002 
     +.112528571428570"+003 
     +.154999999999730"+003 
     +.222946341369190"-007 
        
        
        
 REFERENCES:      
        
     [1] BUS, J. C. P.      
         LINEAR SYSTEMS  WITH CALCULATION OF ERROR BOUNDS  AND ITERATIVE      
         REFINEMENT (DUTCH).
         MATHEMATICAL CENTRE, AMSTERDAM, LR 3. 4. 19 (1972).        
     [2] DEKKER, T. J.      
         ALGOL 60 PROCEDURES IN NUMERICAL ALGEBRA, PART 1.
         MATHEMATICAL CENTRE, AMSTERDAM, TRACT 22 (1968). 
        
        
        
1SECTION: 3.1.1.1.1.1.4       (MAY 1974)                         PAGE 11      
        
        
        
 SOURCE TEXT(S):  
        
 "CODE" 34053;    
     "PROCEDURE" INV(A, N, P); "VALUE" N; "INTEGER" N; "ARRAY" A;   
     "INTEGER" "ARRAY" P;   
     "BEGIN" "INTEGER" J, K, K1;      
         "REAL" R;
         "ARRAY" V[1:N];    
         "REAL" "PROCEDURE" MATMAT(L, U, I, J, A, B); "CODE" 34013; 
         "PROCEDURE" ICHCOL(L, U, I, J, A); "CODE" 34031; 
         "PROCEDURE" DUPCOLVEC(L, U, J, A, B); "CODE" 31034;        
         "FOR" K:= N "STEP" - 1 "UNTIL" 1 "DO"  
         "BEGIN" K1:= K + 1;
             "FOR" J:= N "STEP" - 1 "UNTIL" K1 "DO"       
             "BEGIN" A[J,K1]:= V[J];  
                 V[J]:= - MATMAT(K1, N, K, J, A, A)       
             "END";         
             R:= A[K,K];    
             "FOR" J:= N "STEP" - 1 "UNTIL" K1 "DO"       
             "BEGIN" A[K,J]:= V[J];   
                 V[J]:= - MATMAT(K1, N, J, K, A, A) / R   
             "END";         
             V[K]:= (1 - MATMAT(K1, N, K, K, A, A)) / R   
         "END";   
         DUPCOLVEC(1, N, 1, A, V);    
         "FOR" K:= N - 1 "STEP" - 1 "UNTIL" 1 "DO"        
         "BEGIN" K1:= P[K]; "IF" K1 ^= K "THEN" 
             ICHCOL(1, N, K, K1, A)   
         "END"    
     "END" INV;   
         "EOP"    
0"CODE" 34302;    
     "PROCEDURE" DECINV(A, N, AUX); "VALUE" N; "INTEGER" N;         
     "ARRAY" A, AUX;        
     "BEGIN" "INTEGER" "ARRAY" P[1:N];
         "PROCEDURE" DEC(A, N, AUX, P); "CODE" 34300;     
         "PROCEDURE" INV(A, N, P); "CODE" 34053;
         DEC(A, N, AUX, P); "IF" AUX[3] = N "THEN" INV(A, N, P)     
     "END" DECINV;
         "EOP"    
1SECTION: 3.1.1.1.1.1.4       (MAY 1974)                         PAGE 12      
        
        
        
 "CODE" 34235;    
     "REAL" "PROCEDURE" INV1(A, N, RI, CI, WITHNORM);     
     "VALUE" N, WITHNORM; "INTEGER" N; "BOOLEAN" WITHNORM;
     "ARRAY" A; "INTEGER" "ARRAY" RI, CI;       
     "BEGIN" "INTEGER" L, K, K1;      
         "REAL" AID, NRMINV;
         "PROCEDURE" ICHROW(L, U, I, J, A); "CODE" 34032; 
         "PROCEDURE" INV(A, N, P); "CODE" 34053;
         INV(A, N, RI); NRMINV:= 0; "IF" WITHNORM "THEN"  
         "FOR" L:= 1 "STEP" 1 "UNTIL" N "DO"    
         NRMINV:= NRMINV + ABS(A[L,N]);         
         "FOR" K:= N - 1 "STEP" - 1 "UNTIL" 1 "DO"        
         "BEGIN" "IF" WITHNORM "THEN" 
             "BEGIN" AID:= 0;         
                 "FOR" L:= 1 "STEP" 1 "UNTIL" N "DO"      
                 AID:= AID + ABS(A[L,K]);       
                 "IF" NRMINV < AID "THEN" NRMINV:= AID    
             "END";         
             K1:= CI[K]; "IF" K1 ^= K "THEN" ICHROW(1, N, K, K1, A) 
         "END";   
         INV1:= NRMINV      
     "END" INV1;  
         "EOP"    
0"CODE" 34236;    
     "PROCEDURE" GSSINV(A, N, AUX); "VALUE" N; "INTEGER" N;         
     "ARRAY" A, AUX;        
     "BEGIN" "INTEGER" "ARRAY" RI, CI[1:N];     
         "PROCEDURE" GSSELM(A, N, AUX, RI, CI); "CODE" 34231;       
         "REAL" "PROCEDURE" INV1(A, N, RI, CI, WITHNORM); "CODE" 34235;       
         GSSELM(A, N, AUX, RI, CI);   
         "IF" AUX[3] = N "THEN" AUX[9]:= INV1(A, N, RI, CI, "TRUE") 
     "END" GSSINV;
         "EOP"    
0"CODE" 34244;    
     "PROCEDURE" GSSINVERB(A, N, AUX); "VALUE" N; "INTEGER" N;      
     "ARRAY" A, AUX;        
     "BEGIN" "INTEGER" "ARRAY" RI, CI[1:N];     
         "PROCEDURE" GSSELM(A, N, AUX, RI, CI); "CODE" 34231;       
         "REAL" "PROCEDURE" INV1(A, N, RI, CI, WITHNORM); "CODE" 34235;       
         "PROCEDURE" ERBELM(N, AUX, NRMINV); "CODE" 34241;
         GSSELM(A, N, AUX, RI, CI);   
         "IF" AUX[3] = N "THEN"       
         ERBELM(N, AUX, INV1(A, N, RI, CI, "TRUE"))       
     "END" GSSINVERB;       
         "EOP"    
1SECTION: 3.1.1.1.1.1.5       (FEBRUARY 1979)                     PAGE 1      
        
        
        
 AUTHOR: J. C. P. BUS.      
        
        
 CONTRIBUTOR: J.C.P. BUS AND  P. A. BEENTJES.   
        
        
 INSTITUTE: MATHEMATICAL CENTRE.      
        
        
 RECEIVED: 731008.
        
        
 BRIEF DESCRIPTION: 
        
     THIS   SECTION  CONTAINS   FOUR   PROCEDURES  FOR   CALCULATING  AN      
     ITERATIVELY  IMPROVED  SOLUTION  OF  A SYSTEM  OF LINEAR EQUATIONS:      
     ITISOL SOLVES A LINEAR SYSTEM WHOSE MATRIX HAS BEEN TRIANGULARLY 
     DECOMPOSED BY GSSELM OR GSSERB. THIS SOLUTION THUS OBTAINED    
     NUMERICALLY, IS IMPROVED ITERATIVELY;      
     GSSITISOL SOLVES A LINEAR SYSTEM AND THIS SOLUTION THUS OBTAINED 
     NUMERICALLY, IS IMPROVED ITERATIVELY;      
     ITISOLERB SOLVES A LINEAR SYSTEM WHOSE MATRIX HAS BEEN TRIANGULARLY      
     DECOMPOSED BY GSSNRI.THIS SOLUTION IS IMPROVED ITERATIVELY.MOREOVER      
     A REALISTIC  UPPERBOUND FOR THE RELATIVE ERROR IN THE SOLUTION IS
     CALCULATED.  
     GSSITISOLERB SOLVES A LINEAR SYSTEM.THIS SOLUTION IS IMPROVED  
     ITERATIVELY AND A REALISTIC UPPERBOUND FOR THE RELATIVE ERROR IN 
     THE SOLUTION IS CALCULATED;      
        
        
 KEYWORDS:
        
     ALGEBRAIC EQUATIONS,   
     LINEAR SYSTEMS,
     ITERATIVE REFINEMENT.  
        
        
        
1SECTION: 3.1.1.1.1.1.5       (FEBRUARY 1979)                     PAGE 2      
        
        
        
 SUBSECTION: ITISOL  .      
        
 CALLING SEQUENCE:
        
     THE HEADING OF THIS PROCEDURE IS:
     "PROCEDURE" ITISOL(A, LU, N, AUX, RI, CI, B); "VALUE" N;       
     "INTEGER" N; "ARRAY" A, LU, AUX, B; "INTEGER" "ARRAY" RI, CI;  
     "CODE" 34250;
        
     THE MEANING OF THE FORMAL PARAMETERS IS:   
     A:      <ARRAY IDENTIFIER>;      
             "ARRAY" A[1:N,1:N];      
             ENTRY:THE MATRIX OF THE LINEAR SYSTEM;       
     LU:     <ARRAY IDENTIFIER>;      
             "ARRAY" LU[1:N, 1:N];    
             ENTRY:THE TRIANGULARLY DECOMPOSED FORM OF THE MATRIX GIVEN       
             IN A, AS DELIVERED BY GSSELM  (SECTION 3.1.1.1.1.1.1); 
     N:      <ARITHMETIC EXPRESSION>; 
             THE ORDER OF THE MATRIX; 
     AUX:    <ARRAY IDENTIFIER>;      
             "ARRAY" AUX[10:13];      
             ENTRY: 
             AUX[10]: A RELATIVE TOLERANCE  FOR THE SOLUTION VECTOR;  IF      
                     THE  1-NORM  OF THE VECTOR  OF CORRECTIONS  TO  THE      
                     SOLUTION,  DIVIDED BY  THE 1-NORM OF THE CALCULATED      
                     SOLUTION, IS SMALLER THAN AUX[10], THEN THE PROCESS      
                     WILL STOP; THE USER SHOULD NOT CHOOSE  THE VALUE OF      
                     AUX[10] SMALLER THAN THE RELATIVE PRECISION  OF THE      
                     ELEMENTS OF THE MATRIX  AND THE RIGHT-HAND SIDE  OF      
                     THE LINEAR SYSTEM; 
             AUX[12]: THE MAXIMUM NUMBER  OF ITERATIONS  ALLOWED FOR THE      
                     REFINEMENT  OF  THE  SOLUTION;  IF  THE  NUMBER  OF      
                     ITERATIONS EXCEEDS THE VALUE OF  AUX[12],  THEN THE      
                     PROCESS  WILL  BE BROKEN OFF;  USUALLY  AUX[12] = 5      
                     WILL GIVE GOOD RESULTS;    
             EXIT:
             AUX[11]: THE 1-NORM  OF  THE VECTOR  OF  CORRECTIONS TO THE      
                     SOLUTION IN THE LAST ITERATION STEP, DIVIDED BY THE      
                     1-NORM OF THE CALCULATED SOLUTION;   
                     IF  AUX[11] > AUX[10],  THEN  THE PROCESS  HAS BEEN      
                     BROKEN OFF,   BECAUSE  THE  NUMBER  OF  ITERATIONS       
                     EXCEEDED THE VALUE GIVEN IN AUX[12]; 
             AUX[13]: THE 1-NORM OF THE RESIDUAL VECTOR  (SEE METHOD AND      
                     PERFORMANCE;     
     RI:     <ARRAY IDENTIFIER>;      
             "INTEGER""ARRAY" RI[1:N];
             ENTRY:THE PIVOTAL ROW-INDICES,  AS PRODUCED BY  GSSELM;
     CI:     <ARRAY IDENTIFIER>;      
             "INTEGER""ARRAY" CI[1:N];
             ENTRY:THE PIVOTAL COLUMN-INDICES, AS PRODUCED BY  GSSELM;
     B:      <ARRAY IDENTIFIER>;      
             "ARRAY" B[1:N];
             ENTRY:  THE RIGHT-HAND SIDE OF THE LINEAR SYSTEM;      
             EXIT:   THE CALCULATED  SOLUTION OF THE LINEAR SYSTEM. 
        
1SECTION: 3.1.1.1.1.1.5       (FEBRUARY 1979)                     PAGE 3      
        
        
        
 PROCEDURES USED: 
        
     SOLELM      = CP34061, 
     INIVEC      = CP31010, 
     DUPVEC      = CP31030, 
     LNGMATVEC   = CP34411. 
        
 REQUIRED CENTRAL MEMORY:   
        
     TWO REAL ARRAYS, BOTH OF ORDER N, ARE DECLARED.      
        
 RUNNING TIME: THE NUMBER OF  ARITHMETICAL OPERATIONS  IN EACH ITERATION      
               STEP IS PROPORTIONAL TO N ** 2.  
        
 METHOD AND PERFORMANCE:    
        
     ITISOL   SHOULD  BE  CALLED  AFTER   GSSELM   OR   GSSERB  (SECTION      
     3.1.1.1.1.1.1) AND SOLVES THE LINEAR SYSTEM WITH A MATRIX  GIVEN 
     IN  ARRAY  A, AND A RIGHT-HAND SIDE  GIVEN IN  ARRAY  B; ONCE A
     SOLUTION IS CALCULATED WITH  SOLELM  (SECTION 3.1.1.1.1.1.3),  THIS      
     SOLUTION WILL BE REFINED ITERATIVELY  UNTIL THE CALCULATED RELATIVE      
     CORRECTION  TO THIS SOLUTION  WILL BE LESS THAN  A PRESCRIBED VALUE      
     (SEE AUX[10]); 
     EACH ITERATION  OF THE REFINEMENT PROCESS CONSISTS OF THE FOLLOWING      
     THREE STEPS (SEE [1], [2], [3]): 
     1   CALCULATE, IN DOUBLE PRECISION, THE RESIDUAL VECTOR R,  DEFINED      
         BY:             R = AX - B,  
         WHERE   X  DENOTES  THE  SOLUTION,  OBTAINED  IN  THE  PREVIOUS      
         ITERATION,  B  THE RIGHT-HAND SIDE OF THE LINEAR SYSTEM,  GIVEN      
         IN B[1:N], AND A THE MATRIX GIVEN IN A[1:N, 1:N];
     2   CALCULATE THE SOLUTION C, SAY, OF THE LINEAR SYSTEM: AC = R, 
         WITH THE AID OF THE TRIANGULARLY DECOMPOSED MATRIX  AS GIVEN IN      
         LU[1:N, 1:N];      
     3   CALCULATE THE NEW SOLUTION: XNEW = X - C;
     CONDITION OF THE MATRIX IS NOT TOO BAD,  THEN  THE PRECISION OF THE      
     CALCULATED SOLUTION WILL BE OF THE ORDER OF THE PRECISION ASKED FOR      
     IN  AUX[10];  HOWEVER,  IF THE CONDITION OF THE MATRIX IS VERY BAD,      
     THEN  THIS PROCESS WILL POSSIBLY  NOT  CONVERGE OR,  IN EXCEPTIONAL      
     CASES,  CONVERGE TO  A USELESS RESULT;  IF  THE USER  WANTS TO MAKE      
     CERTAIN  ABOUT  THE PRECISION  OF  THE CALCULATED SOLUTION, THEN HE      
     HAS TO USE  ITISOLERB  (THIS SECTION),  WHICH NEEDS THE CALCULATION      
     (OF ORDER  N ** 3) OF THE INVERSE MATRIX  TO GET AN UPPER BOUND FOR      
     THE CONDITION NUMBER  AND  WHICH  GIVES A REALISTIC UPPER BOUND FOR      
     THE RELATIVE ERROR IN THE CALCULATED SOLUTION;       
     ITISOL  LEAVES  A,  LU,  RI  AND CI UNALTERED, SO AFTER ONE CALL OF      
     GSSELM SEVERAL CALLS OF ITISOL MAY FOLLOW TO CALCULATE THE SOLUTION      
     OF SEVERAL LINEAR SYSTEMS WITH THE SAME MATRIX BUT DIFFERENT RIGHT-      
     HAND SIDES.  
        
 EXAMPLE OF USE: SEE GSSITISOL (THIS SECTION).  
1SECTION: 3.1.1.1.1.1.5       (FEBRUARY 1979)                     PAGE 4      
        
        
        
 SUBSECTION: GSSITISOL   .  
        
        
 CALLING SEQUENCE:
        
     THE HEADING OF THIS PROCEDURE IS:
     "PROCEDURE" GSSITISOL(A, N, AUX, B); "VALUE" N;      
     "INTEGER" N; "ARRAY" A, AUX, B; "CODE" 34251;
        
     THE MEANING OF THE FORMAL PARAMETERS IS:   
     A:      <ARRAY IDENTIFIER>;      
             "ARRAY" A[1:N,1:N];      
             ENTRY:THE N-TH ORDER MATRIX;       
             EXIT: THE  CALCULATED  LOWER-TRIANGULAR  MATRIX  AND  UNIT       
             UPPERTRIANGULAR  MATRIX  WITH  ITS  UNIT  DIAGONAL OMITTED;      
     N:      <ARITHMETIC EXPRESSION>; 
             THE ORDER OF THE MATRIX; 
     AUX:    <ARRAY IDENTIFIER>;      
             "ARRAY" AUX[1:13];       
             ENTRY: 
             AUX[2]: A RELATIVE TOLERANCE FOR THE PROCESS OF TRIANGULAR       
                     DECOMPOSITION; A REASONABLE CHOICE FOR THIS VALUE IS     
                     AN ESTIMATE OF THE RELATIVE PRECISION OF THE MATRIX      
                     ELEMENTS; HOWEVER,IT SHOULD NOT BE CHOSEN SMALLER
                     THANN THE MACHINE PRECISION; 
             AUX[4]: A VALUE WHICH IS USED FOR CONTROLLING PIVOTING (SEE      
                     GSSELM, SECTION 3.1.1.1.1.1.1);      
             AUX[10]: A RELATIVE TOLERANCE  FOR THE SOLUTION VECTOR;  IF      
                     THE  1-NORM  OF THE VECTOR  OF CORRECTIONS  TO  THE      
                     SOLUTION,  DIVIDED BY  THE 1-NORM OF THE CALCULATED      
                     SOLUTION, IS SMALLER THAN AUX[10], THEN THE PROCESS      
                     WILL STOP; THE USER SHOULD NOT CHOOSE  THE VALUE OF      
                     AUX[10] SMALLER THAN THE RELATIVE PRECISION  OF THE      
                     ELEMENTS OF THE MATRIX  AND THE RIGHT-HAND SIDE  OF      
                     THE LINEAR SYSTEM; 
             AUX[12]: THE MAXIMUM NUMBER  OF ITERATIONS  ALLOWED FOR THE      
                     REFINEMENT  OF  THE  SOLUTION;  IF  THE  NUMBER  OF      
                     ITERATIONS EXCEEDS THE VALUE OF  AUX[12],  THEN THE      
                     PROCESS  WILL  BE BROKEN OFF;  USUALLY  AUX[12] = 5      
                     WILL GIVE GOOD RESULTS;    
             EXIT:
             AUX[1]: IF  R  IS THE NUMBER OF ELIMINATION STEPS PERFORMED      
                     (SEE  AUX[3]),   THEN   AUX[1]  EQUALS   1  IF  THE      
                     DETERMINANT OF THE PRINCIPAL SUBMATRIX  OF ORDER  R      
                     IS POSITIVE, ELSE AUX[1] = -1;       
             AUX[3]: THE  NUMBER  OF  ELIMINATION  STEPS  PERFORMED;  IF      
                     AUX[3] < N THEN THE PROCESS HAS BEEN BROKEN OFF AND      
                     NO SOLUTION WILL HAVE BEEN CALCULATED; 
             AUX[5]: THE  MODULUS  OF  AN ELEMENT  WHICH  IS  OF MAXIMUM      
                     ABSOLUTE VALUE FOR THE MATRIX WHICH HAD BEEN GIVEN       
                     IN ARRAY A;      
1SECTION: 3.1.1.1.1.1.5       (FEBRUARY 1979)                     PAGE 5      
        
        
        
             AUX[7]: AN UPPER BOUND FOR THE GROWTH (SEE GSSELM,  SECTION      
                     3.1.1.1.1.1.1);  
             AUX[11]: IF AUX[3] < N,  THEN  AUX[11]  WILL  BE UNDEFINED,      
                     ELSE  AUX[11]  EQUALS  THE 1-NORM  OF THE VECTOR OF      
                     CORRECTIONS  TO  THE  SOLUTION  IN  THE  LAST STEP,      
                     DIVIDED BY THE 1-NORM OF THE CALCULATED SOLUTION;
                     IF  AUX[11] > AUX[10],  THEN  THE PROCESS  HAS BEEN      
                     BROKEN OFF,   BECAUSE   THE  NUMBER  OF  ITERATIONS      
                     EXCEEDED THE VALUE GIVEN IN AUX[12]; 
             AUX[13]: IF AUX[3] = N,  THEN  THE VALUE  OF  AUX[13]  WILL      
                     EQUAL THE 1-NORM OF THE RESIDUAL VECTOR (SEE ITISOL      
                     IN THIS SECTION), ELSE AUX[13] WILL BE UNDEFINED;
     B:      <ARRAY IDENTIFIER>;      
             "ARRAY" B[1:N];
             ENTRY:  THE RIGHT-HAND SIDE OF THE LINEAR SYSTEM;      
             EXIT:   IF AUX[3] = N,  THEN THE CALCULATED SOLUTION OF THE      
                     LINEAR SYSTEM IS OVERWRITTEN ON B,  ELSE  B REMAINS      
                     UNALTERED.       
        
        
 PROCEDURES USED: 
        
     DUPMAT  = CP31035,     
     GSSELM  = CP34231,     
     ITISOL  = CP34250.     
        
        
 REQUIRED CENTRAL MEMORY:   
        
     THREE REAL ARRAYS, TWO OF ORDER N AND ONE OF ORDER N ** 2, ARE 
     DECLARED. FURTHERMORE, TWO INTEGER ARRAYS OF ORDER N ARE USED. 
        
        
 RUNNING TIME: PROPORTIONAL TO N ** 3.
        
        
 METHOD AND PERFORMANCE:    
        
     GSSITISOL   USES   GSSELM  (SECTION  3.1.1.1.1.1.1)  TO  PERFORM  A      
     TRIANGULAR DECOMPOSITION OF THE MATRIX AND ITISOL (THIS SECTION) TO      
     CALCULATE  AN  ITERATIVELY  REFINED  SOLUTION  OF  THE GIVEN LINEAR      
     SYSTEM; IF AUX[3] < N, THEN THE EFFECT OF GSSITISOL  IS MERELY THAT      
     OF  GSSELM;  IF THE CONDITION OF THE MATRIX IS VERY BAD,  THEN,  IN      
     EXCEPTIONAL CASES,  THE CALCULATED SOLUTION  MAY  BE  USELESS  (SEE      
     ITISOL, IN THIS SECTION).
1SECTION: 3.1.1.1.1.1.5       (FEBRUARY 1979)                     PAGE 6      
        
        
        
 EXAMPLE OF USE:  
        
     LET A BE THE FOURTH ORDER SEGMENT OF THE HILBERT MATRIX,       
     MULTIPLIED WITH 840 TO GET INTEGER ELEMENTS, AND B THE THIRD COLUMN      
     OF A, THEN THE SOLUTION OF THE LINEAR SYSTEM AX = B IS GIVEN BY
     THE THIRD UNIT VECTOR AND MAY BE OBTAINED BY THE FOLLOWING     
     PROGRAM:     
        
     "BEGIN" "INTEGER" I, J;
         "ARRAY" A[1:4, 1:4], B[1:4], AUX[1:13];
         "PROCEDURE" GSSITISOL(A, N, AUX, B); "CODE" 34251; 
         "PROCEDURE" LIST(ITEM); "PROCEDURE" ITEM;
         "BEGIN" "INTEGER" I; 
             "FOR" I:= 1 "STEP" 1 "UNTIL" 4 "DO" ITEM(B[I]);
             "FOR" I:= 1 "STEP" 2 "UNTIL" 7 "DO" ITEM(AUX[I]);      
             ITEM(AUX[11]); ITEM(AUX[13])       
         "END" LIST;
         "PROCEDURE" LAYOUT;
         FORMAT("("*, "("SOLUTION:")"B+.15D"+3D,/,3(10B+.15D"+3D,/),
         "("SIGN(DET) = ")"+D,/,"("NUMBER OF ELIMINATIONSTEPS = ")" 
         +D,/,"("MAX(ABS(A[I,J]))= ")"+.15D"+3D,/,
         "("UPPER BOUND GROWTH: ")"+.15D"+3D,/, 
         "("NORM LAST CORRECTION VECTOR: ")"+.15D"+3D,/,  
         "("NORM RESIDUAL VECTOR: ")"+.15D"+3D")");       
        
         "FOR" I:= 1 "STEP" 1 "UNTIL" 4 "DO"    
         "BEGIN" "FOR" J:= 1 "STEP" 1 "UNTIL" 4 "DO"      
             A[I,J]:= 840 // (I + J - 1); B[I]:= A[I,3]   
         "END";   
         AUX[2]:= "-14; AUX[4]:= 8; AUX[10]:= "-14; AUX[12]:= 5;    
         GSSITISOL(A, 4, AUX, B);     
         OUTLIST(71, LAYOUT, LIST)    
     "END"
        
     RESULTS:     
        
     SOLUTION: +.000000000000000"+000 
               +.000000000000000"+000 
               +.100000000000000"+001 
               +.000000000000000"+000 
     SIGN(DET) = +1 
     NUMBER OF ELIMINATIONSTEPS = +4  
     MAX(ABS(A[I,J]))= +.840000000000000"+003   
     UPPER BOUND GROWTH: +.134080000000000"+004 
     NORM LAST CORRECTION VECTOR: +.000000000000000"+000  
     NORM RESIDUAL VECTOR: +.000000000000000"+000 
        
        
        
1SECTION: 3.1.1.1.1.1.5       (FEBRUARY 1979)                     PAGE 7      
        
        
        
 SUBSECTION: ITISOLERB.     
        
        
 CALLING SEQUENCE:
        
     THE HEADING OF THIS PROCEDURE IS:
     "PROCEDURE" ITISOLERB(A, LU, N, AUX, RI, CI, B); "VALUE" N;    
     "INTEGER" N; "ARRAY" A, LU, AUX, B; "INTEGER" "ARRAY" RI, CI;  
     "CODE" 34253;
        
     THE MEANING OF THE FORMAL PARAMETERS IS:   
     A:      <ARRAY IDENTIFIER>;      
             "ARRAY"[1:N,1:N];
             ENTRY:THE  MATRIX  OF  THE  LINEAR  SYSTEM;  
     LU:     <ARRAY IDENTIFIER>;      
             "ARRAY" LU[1:N,1:N];     
             ENTRY:THE TRIANGULARLY DECOMPOSED FORM OF THE MATRIX GIVEN       
             IN A AS DELIVERED BY GSSNRI  (SECTION 3.1.1.1.1.1.1);  
     N:      <ARITHMETIC EXPRESSION>; 
             THE ORDER OF THE MATRIX; 
     AUX:    <ARRAY IDENTIFIER>;      
             "ARRAY" AUX[0:13];       
             ENTRY: 
             AUX[0]: THE MACHINE PRECISION;     
             AUX[5]: THE  MODULUS  OF  AN  ELEMENT,  WHICH IS OF MAXIMUM      
                     ABSOLUTE VALUE FOR THE MATRIX OF THE LINEAR SYSTEM;      
                     THIS  VALUE  IS  DELIVERED   BY   GSSNRI   (SECTION      
                     3.1.1.1.1.1.1) IN AUX[5];  
             AUX[6]: AN  UPPER  BOUND  FOR  THE  RELATIVE  ERROR  IN THE      
                     ELEMENTS OF THE MATRIX OF THE LINEAR SYSTEM;   
             AUX[7]: AN  UPPER  BOUND  FOR  THE GROWTH  DURING  GAUSSIAN      
                     ELIMINATION; THIS VALUE IS DELIVERED BY  GSSNRI  IN      
                     AUX[7];
             AUX[8]: AN  UPPER  BOUND  FOR  THE  RELATIVE  ERROR  IN THE      
                     ELEMENTS  OF  THE  RIGHT-HAND  SIDE  OF  THE LINEAR      
                     SYSTEM;
             AUX[9]: THE 1-NORM  OF  THE INVERSE MATRIX;  THIS  VALUE IS      
                     DELIVERED BY GSSNRI IN AUX[9];       
             AUX[10]: A RELATIVE TOLERANCE  FOR THE SOLUTION VECTOR;  IF      
                     THE  1-NORM  OF THE VECTOR  OF CORRECTIONS  TO  THE      
                     SOLUTION,  DIVIDED BY  THE 1-NORM OF THE CALCULATED      
                     SOLUTION, IS SMALLER THAN AUX[10], THEN THE PROCESS      
                     WILL STOP; THE USER SHOULD NOT CHOOSE  THE VALUE OF      
                     AUX[10] SMALLER THAN THE RELATIVE PRECISION  OF THE      
                     ELEMENTS OF THE MATRIX  AND THE RIGHT-HAND SIDE  OF      
                     THE LINEAR SYSTEM, GIVEN IN AUX[6] AND AUX[8]; 
             AUX[12]: THE MAXIMUM NUMBER  OF ITERATIONS  ALLOWED FOR THE      
                     REFINEMENT  OF  THE  SOLUTION;  IF  THE  NUMBER  OF      
                     ITERATIONS EXCEEDS THE VALUE OF  AUX[12],  THEN THE      
                     PROCESS  WILL  BE BROKEN OFF;  USUALLY  AUX[12] = 5      
                     WILL GIVE GOOD RESULTS;    
1SECTION: 3.1.1.1.1.1.5       (FEBRUARY 1979)                     PAGE 8      
        
        
        
             EXIT:
             AUX[11]:A REALISTIC UPPERBOUND FOR THE RELATIVE ERROR IN 
                     THE CALCULATED SOLUTION; HOWEVER, IF NO USE CAN BE       
                     MADE OF THE ERROR-FORMULA, THEN AUX[11] := -1; 
             AUX[13]:THE 1-NORM OF THE RESIDUAL VECTOR  (SEE METHOD AND       
                     PERFORMANCE);    
     RI:     <ARRAY IDENTIFIER>;      
             "INTEGER""ARRAY" RI[1:N];
             ENTRY:THE PIVOTAL ROW-INDICES,  AS PRODUCED BY  GSSNRI;
     CI:     <ARRAY IDENTIFIER>;      
             "INTEGER""ARRAY" CI[1:N];
             EXIT:THE PIVOTAL COLUMN-INDICES, AS PRODUCED BY  GSSNRI; 
     B:      <ARRAY IDENTIFIER>;      
             "ARRAY" B[1:N];
             ENTRY:  THE RIGHT-HAND SIDE OF THE LINEAR SYSTEM;      
             EXIT:   THE SOLUTION OF THE LINEAR SYSTEM.   
        
 PROCEDURES USED: 
        
     ITISOL  = CP34250.     
        
 REQUIRED CENTRAL MEMORY:TWO REAL AARRYS, BOTH OF ORDER N, ARE DECLARED.      
        
 RUNNING TIME: THE NUMBER OF  ARITHMETICAL OPERATIONS  IN EACH ITERATION      
               STEP OF THE REFINEMENT PROCESS IS PROPORTIONAL TO N ** 2.      
        
 METHOD AND PERFORMANCE:    
        
     ITISOLERB  SHOULD BE CALLED  AFTER  GSSNRI (SECTION 3.1.1.1.1.1.1),      
     WHICH DELIVERS  THE TRIANGULARLY DECOMPOSED FORM OF THE MATRIX  AND      
     THE PROPER VALUES  FOR THE ODD ELEMENTS  OF  ARRAY  AUX;  ITISOLERB      
     CALCULATES, WITH THE USE OF  ITISOL (THIS SECTION),  AN ITERATIVELY      
     IMPROVED SOLUTION OF THE LINEAR SYSTEM  WITH A MATRIX  AS  GIVEN IN      
     ARRAY A AND A RIGHT-HAND SIDE AS GIVEN IN  B;  MOREOVER,  ITISOLERB      
     CALCULATES  A REALISTIC UPPER BOUND  FOR THE RELATIVE ERROR  IN THE      
     CALCULATED  SOLUTION, BY (SEE [1], [2]):   
         NORM(DX) / NORM(X) <= P / (1 - P),     
     WHERE : P = ( NORM(R) / NORM(X) + DB / NORM(X) + DA )
         * NORM(C) / (1 - Q * NORM(C) ) 
         FOR Q SEE SECTION 3.1.1.1.1.1.1 (SUBSECTION ERBELM),       
         R IS THE RESIDUAL VECTOR (SEE ITISOL IN THIS SECTION),     
         X IS THE CALCULATED SOLUTION,
         DB IS THE UPPER BOUND FOR THE RELATIVE ERROR  IN THE RIGHT-HAND      
         SIDE ELEMENTS,     
         DA IS THE UPPER BOUND FOR THE RELATIVE ERROR IN THE MATRIX 
         ELEMENTS, C IS THE CALCULATED INVERSE MATRIX,    
         AND THE 1-NORM OF A VECTOR OR A MATRIX IS DENOTED BY: NORM(.)
1SECTION: 3.1.1.1.1.1.5       (FEBRUARY 1979)                     PAGE 9      
        
        
        
     IF 1 - P < AUX[0], THEN THE VALUE -1 IS DELIVERED IN AUX[11];  
     ITISOLERB  LEAVES A, LU, RI AND CI UNALTERED,  SO AFTER ONE CALL OF      
     GSSNRI  SEVERAL CALLS OF  ITISOLERB  MAY FOLLOW,  TO CALCULATE  THE      
     SOLUTION  OF  SEVERAL  LINEAR  SYSTEMS  WITH  THE SAME  MATRIX  BUT      
     DIFFERENT RIGHT-HAND SIDES.      
        
 EXAMPLE OF USE: SEE GSSITISOLERB (THIS SECTION). 
        
        
        
 SUBSECTION: GSSITISOLERB.  
        
        
 CALLING SEQUENCE:
        
     THE HEADING OF THIS PROCEDURE IS:
     "PROCEDURE" GSSITISOLERB(A, N, AUX, B); "VALUE" N;   
     "INTEGER" N; "ARRAY" A, AUX, B; "CODE" 34254;
        
     THE MEANING OF THE FORMAL PARAMETERS IS:   
     A:      <ARRAY IDENTIFIER>;      
             "ARRAY" A[1:N,1:N];      
             ENTRY:THE N-TH ORDER MATRIX;       
             EXIT: THE  CALCULATED  LOWER-TRIANGULAR  MATRIX  AND  UNIT       
             UPPER-TRIANGULAR  MATRIX WITH ITS  UNIT  DIAGONAL  OMITTED;      
     N:      <ARITHMETIC EXPRESSION>; 
             THE ORDER OF THE MATRIX; 
     AUX:    <ARRAY IDENTIFIER>;      
             "ARRAY" AUX[0:13];       
             ENTRY: 
             AUX[0]: THE MACHINE PRECISION;     
             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;     
             AUX[4]: A VALUE WHICH IS USED FOR CONTROLLING PIVOTING (SEE      
                     GSSELM, SECTION 3.1.1.1.1.1.1);      
             AUX[6]: AN UPPER BOUND FOR THE RELATIVE ERROR IN THE MATRIX      
                     ELEMENTS OF THE LINEAR SYSTEM;       
             AUX[8]: AN  UPPER  BOUND  FOR  THE  RELATIVE  ERROR  IN THE      
                     ELEMENTS OF THE RIGHT-HAND SIDE;     
             AUX[10]: A RELATIVE TOLERANCE  FOR THE SOLUTION VECTOR;  IF      
                     THE  1-NORM  OF THE VECTOR  OF CORRECTIONS  TO  THE      
                     SOLUTION,  DIVIDED BY  THE 1-NORM OF THE CALCULATED      
                     SOLUTION, IS SMALLER THAN AUX[10], THEN THE PROCESS      
                     WILL STOP; THE USER SHOULD NOT CHOOSE  THE VALUE OF      
                     AUX[10] SMALLER THAN THE RELATIVE PRECISION  OF THE      
                     ELEMENTS OF THE MATRIX  AND THE RIGHT-HAND SIDE  OF      
                     THE LINEAR SYSTEM (AUX[10] >= AUX[2]); 
             AUX[12]: THE MAXIMUM NUMBER  OF ITERATIONS  ALLOWED FOR THE      
                     REFINEMENT  OF  THE  SOLUTION;  IF  THE  NUMBER  OF      
                     ITERATIONS EXCEEDS THE VALUE OF  AUX[12],  THEN THE      
                     PROCESS  WILL  BE BROKEN OFF;  USUALLY  AUX[12] = 5      
                     WILL GIVE GOOD RESULTS;    
1SECTION: 3.1.1.1.1.1.5       (FEBRUARY 1979)                    PAGE 10      
        
        
        
             EXIT:
             AUX[1]: IF  R  IS THE NUMBER OF ELIMINATION STEPS PERFORMED      
                     (SEE  AUX[3]),   THEN   AUX[1]  EQUALS   1  IF  THE      
                     DETERMINANT OF THE PRINCIPAL SUBMATRIX  OF ORDER  R      
                     IS POSITIVE, ELSE AUX[1] = -1;       
             AUX[3]: THE  NUMBER  OF  ELIMINATION  STEPS  PERFORMED;  IF      
                     AUX[3] < N  THEN  THE PROCESS  HAS BEEN BROKEN OFF       
                     AND NO SOLUTION WILL HAVE BEEN CALCULATED;     
             AUX[5]: THE  MODULUS  OF  AN ELEMENT  WHICH  IS  OF MAXIMUM      
                     ABSOLUTE VALUE FOR THE MATRIX WHICH HAD BEEN GIVEN       
                     IN ARRAY A;      
             AUX[7]: AN UPPER BOUND FOR THE GROWTH (SEE GSSELM,  SECTION      
                     3.1.1.1.1.1.1);  
             AUX[9]: IF AUX[3] = N THEN AUX[9] EQUALS THE 1-NORM  OF THE      
                     CALCULATED  INVERSE  MATRIX,  ELSE  AUX[9]  WILL BE      
                     UNDEFINED;       
             AUX[11]: IF AUX[3] < N,  THEN  AUX[11]  WILL  BE UNDEFINED,      
                     ELSE THE VALUE OF AUX[11]  EQUALS A REALISTIC UPPER      
                     BOUND  FOR  THE RELATIVE ERROR  IN  THE  CALCULATED      
                     SOLUTION;  HOWEVER,  IF NO USE  CAN BE MADE  OF THE      
                     ERROR FORMULA (SEE ITISOLERB IN THIS SECTION), THEN      
                     AUX[11]:= -1;    
             AUX[13]: IF AUX[3] = N,  THEN  AUX[13] EQUALS THE 1-NORM OF      
                     THE RESIDUAL VECTOR  (SEE  ITISOL IN THIS SECTION),      
                     ELSE AUX[13] WILL BE UNDEFINED;      
     B:      <ARRAY IDENTIFIER>;      
             "ARRAY" B[1:N];
             ENTRY:  THE RIGHT-HAND SIDE OF THE LINEAR SYSTEM;      
             EXIT:   IF AUX[3] = N,  THEN THE CALCULATED SOLUTION OF THE      
                     LINEAR SYSTEM IS OVERWRITTEN ON B,  ELSE  B REMAINS      
                     UNALTERED.       
        
 PROCEDURES USED: 
        
     DUPMAT      = CP31035, 
     GSSNRI      = CP34252, 
     ITISOLERB   = CP34253. 
        
 REQUIRED CENTRAL MEMORY:   
        
     THREE REAL ARRAYS, TWO OF ORDER N AND ONE OF ORDER N ** 2 , ARE
     USED. FURTHERMORE, TWO INTEGER ARRAYS OF ORDER N ARE USED.     
        
 RUNNING TIME: PROPORTIONAL TO N ** 3.
        
 METHOD AND PERFORMANCE:    
        
     GSSITISOLERB USES  GSSNRI  (SECTION  3.1.1.1.1.1.1)  TO  PERFORM  A      
     TRIANGULAR DECOMPOSITION OF THE MATRIX AND ITISOLERB (THIS SECTION)      
     TO CALCULATE AN ITERATIVELY REFINED  SOLUTION  OF  THE GIVEN LINEAR      
     SYSTEM  AND  A REALISTIC UPPER BOUND FOR THE RELATIVE ERROR IN THIS      
     SOLUTION; IF AUX[3] < N, THEN THE EFFECT OF  GSSITISOLERB IS MERELY      
     THAT OF GSSELM (SECTION 3.1.1.1.1.1.1).    
1SECTION: 3.1.1.1.1.1.5       (FEBRUARY 1979)                    PAGE 11      
        
        
        
 EXAMPLE OF USE:  
        
     LET A BE THE FOURTH ORDER SEGMENT OF THE HILBERT MATRIX,       
     MULTIPLIED WITH 840 TO GET INTEGER ELEMENTS, AND B THE THIRD COLUMN      
     OF A, THEN THE SOLUTION OF THE LINEAR SYSTEM AX = B IS GIVEN BY
     THE THIRD UNIT VECTOR  AND  THIS SOLUTION,  AS WELL AS  A REALISTIC      
     UPPER BOUND FOR THE RELATIVE ERROR IN IT,  MAY BE OBTAINED BY THE
     FOLLOWING PROGRAM:     
        
     "BEGIN" "INTEGER" I, J;
         "ARRAY" A[1:4, 1:4], B[1:4], AUX[0:13];
         "PROCEDURE" GSSITISOLERB(A, N, AUX, B); "CODE" 34254;      
         "PROCEDURE" LIST(ITEM); "PROCEDURE" ITEM;
         "BEGIN" "INTEGER" I; 
             "FOR" I:= 1 "STEP" 1 "UNTIL" 4 "DO" ITEM(B[I]);
             "FOR" I:= 1 "STEP" 2 "UNTIL" 13 "DO" ITEM(AUX[I])      
         "END" LIST;
         "PROCEDURE" LAYOUT;
         FORMAT("("*, "("SOLUTION:")"B+.15D"+3D,/,3(10B+.15D"+3D,/),
         "("SIGN(DET) = ")"+D,/,"("NUMBER OF ELIMINATIONSTEPS = ")" 
         +D,/,"("MAX(ABS(A[I,J]))= ")"+.15D"+3D,/,
         "("UPPER BOUND GROWTH: ")"+.15D"+3D,/, 
         "("NORM CALCULATED INVERSE MATRIX: ")"+.15D"+3D,/, 
         "("UPPER BOUND FOR THE RELATIVE ERROR: ")"+.15D"+3D,/,     
         "("NORM RESIDUAL VECTOR: ")"+.15D"+3D")");       
        
         "FOR" I:= 1 "STEP" 1 "UNTIL" 4 "DO"    
         "BEGIN" "FOR" J:= 1 "STEP" 1 "UNTIL" 4 "DO"      
             A[I,J]:= 840 // (I + J - 1); B[I]:= A[I,3]   
         "END";   
         AUX[0]:= AUX[2]:= "-14; AUX[4]:= 8; AUX[6]:= AUX[8]:= 0;   
         AUX[10]:= "-14; AUX[12]:= 5; 
         GSSITISOLERB(A, 4, AUX, B);  
         OUTLIST(71, LAYOUT, LIST)    
     "END"
        
     RESULTS:     
        
     SOLUTION: +.000000000000000"+000 
               +.000000000000000"+000 
               +.100000000000000"+001 
               +.000000000000000"+000 
     SIGN(DET) = +1 
     NUMBER OF ELIMINATIONSTEPS = +4  
     MAX(ABS(A[I,J]))= +.840000000000000"+003   
     UPPER BOUND GROWTH: +.134080000000000"+004 
     NORM CALCULATED INVERSE MATRIX: +.162142857143540"+002 
     UPPER BOUND FOR THE RELATIVE ERROR: +.000000000000000"+000     
     NORM RESIDUAL VECTOR: +.000000000000000"+000 
1SECTION: 3.1.1.1.1.1.5       (MAY 1974)                         PAGE 12      
        
        
        
 REFERENCES:      
        
     [1] BUS, J. C. P.      
         LINEAR SYSTEMS  WITH CALCULATION  OF ERROR BOUNDS AND ITERATIVE      
         REFINEMENT (DUTCH).
         MATHEMATICAL CENTRE, AMSTERDAM, LR 3.4.19, (1972). 
     [2] DEKKER, T. J.      
         NUMERICAL ALGEBRA (DUTCH).   
         MATHEMATICAL CENTRE, AMSTERDAM, SYLLABUS 12, (1971).       
     [3] WILKINSON, J. H.   
         THE ALGEBRAIC EIGENVALUE PROBLEM.      
         OXFORD (1965).     
        
        
        
 SOURCE TEXT(S):  
        
 "CODE" 34250;    
     "PROCEDURE" ITISOL(A, LU, N, AUX, RI, CI, B); "VALUE" N;       
     "INTEGER" N; 
     "ARRAY" A, LU, AUX, B; "INTEGER" "ARRAY" RI, CI;     
     "BEGIN" "INTEGER" I, ITER, MAXITER;
         "REAL" MAXERX, ERX, NRMRES, NRMSOL, R, RR;       
         "ARRAY" RES, SOL[1:N];       
         "PROCEDURE" SOLELM(A, N, RI, CI, B); "CODE" 34061; 
         "PROCEDURE" INIVEC(L, U, A, X); "CODE" 31010;    
         "PROCEDURE" DUPVEC(L, U, SHIFT, A, B); "CODE" 31030;       
         "PROCEDURE" LNGMATVEC(A, B, C, D, E, F, G, H, I);
         "CODE" 34411;      
         MAXERX:= ERX:= AUX[10]; MAXITER:= AUX[12];       
         INIVEC(1, N, SOL, 0); DUPVEC(1, N, 0, RES, B);   
         "FOR" ITER:= 1, ITER + 1 "WHILE" ITER <= MAXITER & 
         MAXERX < ERX "DO"  
         "BEGIN" SOLELM(LU, N, RI, CI, RES); ERX:= NRMSOL:= NRMRES:= 0;       
             "FOR" I:= 1 "STEP" 1 "UNTIL" N "DO"
             "BEGIN" R:= RES[I]; ERX:= ERX + ABS(R); RR:= SOL[I] + R; 
                 SOL[I]:= RR; NRMSOL:= NRMSOL + ABS(RR)   
             "END"; 
             ERX:= ERX / NRMSOL;      
             "FOR" I:= 1 "STEP" 1 "UNTIL" N "DO"
             "BEGIN" LNGMATVEC(1, N, I, A, SOL, - B[I], 0, R, RR);  
                 R:= - (R + RR); RES[I]:= R; NRMRES:= NRMRES + ABS(R) 
             "END"
         "END" ITERATION;   
         DUPVEC(1, N, 0, B, SOL);     
         AUX[11]:= ERX; AUX[13]:= NRMRES
     "END" ITISOL;
         "EOP"    
1SECTION: 3.1.1.1.1.1.5       (MAY 1974)                         PAGE 13      
        
        
        
0"CODE" 34251;    
     "PROCEDURE" GSSITISOL(A, N, AUX, B); "VALUE" N; "INTEGER" N;   
     "ARRAY" A, AUX, B;     
     "BEGIN" "INTEGER" I, J;
         "ARRAY" AA[1:N,1:N]; 
         "INTEGER" "ARRAY" RI, CI[1:N]; 
         "PROCEDURE" GSSELM(A, N, AUX, RI, CI); "CODE" 34231;       
         "PROCEDURE" ITISOL(A, LU, N, AUX, RI, CI, B); "CODE" 34250;
         "PROCEDURE" DUPMAT(L, U, I, J, A, B); "CODE" 31035;
         DUPMAT(1, N, 1, N, AA, A);   
         GSSELM(A, N, AUX, RI, CI);   
         "IF" AUX[3] = N "THEN" ITISOL(AA, A, N, AUX, RI, CI, B)    
     "END" GSSITISOL;       
         "EOP"    
        
 "CODE" 34253;    
     "PROCEDURE" ITISOLERB(A, LU, N, AUX, RI, CI, B); "VALUE" N;    
     "INTEGER" N; 
     "ARRAY" A, LU, AUX, B; "INTEGER" "ARRAY" RI, CI;     
     "BEGIN" "INTEGER" I;   
         "REAL" NRMSOL, NRMINV, NRMB, ALFA, TOLA, EPS;    
         "PROCEDURE" ITISOL(A, LU, N, AUX, RI, CI, B); "CODE" 34250;
         EPS:= AUX[0];      
         NRMINV:= AUX[9]; TOLA:= AUX[5] * AUX[6]; NRMB:= NRMSOL:= 0;
         "FOR" I:= 1 "STEP" 1 "UNTIL" N "DO" NRMB:= NRMB + ABS(B[I]); 
         ITISOL(A, LU, N, AUX, RI, CI, B);      
         "FOR" I:= 1 "STEP" 1 "UNTIL" N "DO"    
         NRMSOL:= NRMSOL + ABS(B[I]); 
         ALFA:= 1 - (1.06 * EPS * AUX[7] * (.75 * N + 4.5) * N ** 2 
                + TOLA) * NRMINV;     
         "IF" ALFA < EPS "THEN" AUX[11]:= - 1 "ELSE"      
         "BEGIN" ALFA:= ((AUX[13] + AUX[8] * NRMB) / NRMSOL + TOLA) * 
                        NRMINV / ALFA;
             AUX[11]:= "IF" 1 - ALFA < EPS "THEN" - 1 "ELSE"
                       ALFA / (1 - ALFA)
         "END"    
     "END" ITISOLERB;       
         "EOP"    
0"CODE" 34254;    
     "PROCEDURE" GSSITISOLERB(A, N, AUX, B); "VALUE" N; "INTEGER" N;
     "ARRAY" A, AUX, B;     
     "BEGIN" "INTEGER" I, J;
         "ARRAY" AA[1:N,1:N]; 
         "INTEGER" "ARRAY" RI, CI[1:N]; 
         "PROCEDURE" GSSNRI(A, N, AUX, RI, CI); "CODE" 34252;       
         "PROCEDURE" ITISOLERB(A, LU, N, AUX, RI, CI, B); "CODE" 34253;       
         "PROCEDURE" DUPMAT(L, U, I, J, A, B); "CODE" 31035;
         DUPMAT(1, N, 1, N, AA, A);   
         GSSNRI(A, N, AUX, RI, CI);   
         "IF" AUX[3] = N "THEN" ITISOLERB(AA, A, N, AUX, RI, CI, B) 
     "END" GSSITISOLERB;    
         "EOP"    
1SECTION: 3.1.1.1.1.2.1       (MAY 1974)                          PAGE 1      
        
        
        
 AUTHOR: T.J. DEKKER.       
        
        
 CONTRIBUTORS: S.P.N. VAN KAMPEN, J. KOK.       
        
        
 INSTITUTE: MATHEMATICAL CENTRE.      
        
 RECEIVED: 731015.
        
 BRIEF DESCRIPTION:         
        
     THIS SECTION CONTAINS TWO PROCEDURES:      
     A) CHLDEC2 CALCULATES THE CHOLESKY  DECOMPOSITION OF A POSITIVE
     DEFINITE SYMMETRIC MATRIX WH0SE UPPER TRIANGLE IS GIVEN IN A TWO-        
     DIMENSIONAL ARRAY;     
     B) CHLDEC1 CALCULATES THE CHOLESKY  DECOMPOSITION OF A POSITIVE
     DEFINITE SYMMETRIC MATRIX WHOSE UPPER TRIANGLE IS GIVEN COLUMNWISE       
     IN A ONE-DIMENSIONAL ARRAY.      
        
 KEYWORDS:        
        
     LINEAR EQUATIONS,      
     POSITIVE DEFINITE SYMMETRIC MATRIX,        
     CHOLESKY DECOMPOSITION.
        
        
 SUBSECTION: CHLDEC2.       
        
 CALLING SEQUENCE:
        
     THE HEADING OF THE PROCEDURE IS: 
     "PROCEDURE" CHLDEC2(A, N, AUX); "VALUE" N; "INTEGER" N;        
     "ARRAY" A, AUX;        
        
     THE MEANING OF THE FORMAL PARAMETERS IS:   
     A:      <ARRAY IDENTIFIER>;      
             "ARRAY" A[1:N,1:N];      
             ENTRY: THE UPPER  TRIANGLE OF THE  POSITIVE DEFINITE MATRIX      
                    MUST BE GIVEN IN THE UPPER-TRIANGULAR PART OF A (THE      
                    ELEMENTS A[I,J], I <= J);   
             EXIT:  THE  CHOLESKY  DECOMPOSITION  OF   THE  MATRIX  IS        
                    DELIVERED IN THE UPPER TRIANGLE OF A; 
     N:      <ARITHMETIC EXPRESSION>; 
             THE ORDER OF THE MATRIX; 
     AUX:    <ARRAY IDENTIFIER>;      
             "ARRAY" AUX[2:3];        
             ENTRY: AUX[2]: A  RELATIVE TOLERANCE  USED TO  CONTROL  THE      
                    CALCULATION OF THE DIAGONAL ELEMENTS; 
             NORMAL EXIT:   AUX[3]:= N;         
             ABNORMAL EXIT: IF THE  DECOMPOSITION CANNOT BE  CARRIED OUT      
                    BECAUSE  THE  MATRIX  IS (NUMERICALLY) NOT  POSITIVE      
                    DEFINITE,  AUX[3]:= K - 1, WHERE K IS THE LAST STAGE      
                    NUMBER. 
1SECTION: 3.1.1.1.1.2.1       (MAY 1974)                          PAGE 2      
        
        
        
 PROCEDURES USED: TAMMAT = CP34014.   
        
        
 RUNNING TIME: ROUGHLY PROPORTIONAL TO N CUBED. 
        
        
 LANGUAGE:   ALGOL 60.      
        
        
 METHOD AND PERFORMANCE:    
        
     CHLDEC2 PERFORMS THE CHOLESKY DECOMPOSITION OF A SYMMETRIC     
     POSITIVE DEFINITE MATRIX.        
     THE METHOD USED IS CHOLESKY'S SQUARE ROOT METHOD WITHOUT       
     PIVOTING (SEE REF[1] AND [2]). IF THE GIVEN SYMMETRIC MATRIX IS
     POSITIVE DEFINITE, THE METHOD YIELDS AN UPPER-TRIANGULAR MATRIX U        
     SUCH THAT U'U EQUALS THE GIVEN MATRIX.     
     THE PROCESS IS TERMINATED AT STAGE K, IF THE K-TH DIAGONAL ELEMENT       
     OF THE GIVEN MATRIX MINUS THE SUM OF THE SQUARED ELEMENTS OF THE         
     K-TH COLUMN OF U IS LESS THAN A TOLERANCE TIMES THE MAXIMUM    
     DIAGONAL ELEMENT OF THE GIVEN MATRIX. IN THIS CASE THE MATRIX, 
     POSSIBLY MODIFIED BY ROUNDING ERRORS, IS NOT POSITIVE DEFINITE.
        
        
 REFERENCES:      
        
     [1]. T.J. DEKKER.      
          ALGOL 60 PROCEDURES IN NUMERICAL ALGEBRA, PART 1.         
          MC TRACT 22, 1968, MATH. CENTR., AMSTERDAM.     
        
     [2]. J.H. WILKINSON.   
          THE ALGEBRAIC EIGENVALUE PROBLEM.     
          CLARENDON PRESS, OXFORD, 1965.        
        
        
 EXAMPLE OF USE:  
        
     SEE EXAMPLE OF USE OF CHLINV2, SECTION 3.1.1.1.1.2.4.
        
        
1SECTION: 3.1.1.1.1.2.1       (MAY 1974)                          PAGE 3      
        
        
        
 SUBSECTION: CHLDEC1.       
        
        
 CALLING SEQUENCE:
        
     THE HEADING OF THE PROCEDURE IS: 
     "PROCEDURE" CHLDEC1(A, N, AUX); "VALUE" N; "INTEGER" N;        
     "ARRAY" A, AUX;        
        
     THE MEANING OF THE FORMAL PARAMETERS IS:   
     A:      <ARRAY IDENTIFIER>;      
             "ARRAY" A[1 : (N + 1) * N // 2];   
             ENTRY: THE UPPER-TRIANGULAR PART OF THE  POSITIVE  DEFINITE      
                    SYMMETRIC MATRIX MUST BE GIVEN COLUMNWISE IN ARRAY A      
                    (THE (I,J)-TH ELEMENT OF THE MATRIX MUST BE GIVEN IN      
                    A[(J - 1) * J // 2 + I] FOR 1 <= I <= J <= N);  
             EXIT:  THE  CHOLESKY  DECOMPOSITION  OF  THE   MATRIX  IS        
                    DELIVERED COLUMNWISE IN A.  
     N:      <ARITHMETIC EXPRESSION>; 
             THE ORDER OF THE MATRIX; 
     AUX:    <ARRAY IDENTIFIER>;      
             "ARRAY" AUX[2:3];        
             ENTRY: AUX[2]: A  RELATIVE TOLERANCE  USED TO  CONTROL  THE      
                    CALCULATION OF THE DIAGONAL ELEMENTS; 
             NORMAL EXIT:   AUX[3]:= N;         
             ABNORMAL EXIT: IF THE  DECOMPOSITION CANNOT BE  CARRIED OUT      
                    BECAUSE  THE  MATRIX  IS (NUMERICALLY) NOT  POSITIVE      
                    DEFINITE,  AUX[3]:= K - 1, WHERE K IS THE LAST STAGE      
                    NUMBER. 
        
        
 PROCEDURES USED: VECVEC = CP34010.   
        
        
 RUNNING TIME: ROUGHLY PROPORTIONAL TO N CUBED. 
        
        
 LANGUAGE: ALGOL 60.        
        
        
 METHOD AND PERFORMANCE:    
        
     CHLDEC1 PERFORMS THE CHOLESKY DECOMPOSITION OF A SYMMETRIC     
     POSITIVE DEFINITE MATRIX, WHOSE UPPER TRIANGLE IS STORED IN A ONE-       
     DIMENSIONAL ARRAY, BY CHOLESKY'S SQUARE ROOT METHOD WITHOUT    
     PIVOTING.    
     SEE ALSO METHOD AND PERFORMANCE OF CHLDEC2, (THIS SECTION).    
        
        
 EXAMPLE OF USE:  
        
     SEE EXAMPLE OF USE OF CHLINV1, SECTION 3.1.1.1.1.2.4.
1SECTION: 3.1.1.1.1.2.1       (DECEMBER 1975)                     PAGE 4      
        
        
        
 SOURCE TEXT(S) : 
        
 "CODE" 34310;    
     "PROCEDURE" CHLDEC2(A, N, AUX); "VALUE" N; "INTEGER" N;        
     "ARRAY" A, AUX;        
     "BEGIN" "INTEGER" K, J; "REAL" R, EPSNORM; 
         "REAL" "PROCEDURE" TAMMAT(L, U, I, J, A, B); "CODE" 34014; 
        
         R:= 0;   
         "FOR" K:= 1 "STEP" 1 "UNTIL" N "DO"    
         "IF" A[K,K] > R "THEN" R:= A[K,K];     
         EPSNORM:= AUX[2] * R;        
         "FOR" K:= 1 "STEP" 1 "UNTIL" N "DO"    
         "BEGIN" R:= A[K,K] - TAMMAT(1, K - 1, K, K, A, A);         
             "IF" R <= EPSNORM "THEN" 
             "BEGIN" AUX[3]:= K - 1; "GOTO" END "END";    
             A[K,K]:= R:= SQRT(R);    
             "FOR" J:= K + 1 "STEP" 1 "UNTIL" N "DO"      
             A[K,J]:= (A[K,J] - TAMMAT(1, K - 1, J, K, A, A)) / R   
         "END";   
         AUX[3]:= N;        
     END:         
     "END" CHLDEC2;         
         "EOP"    
        
 "CODE" 34311;    
     "PROCEDURE" CHLDEC1(A, N, AUX); "VALUE" N; "INTEGER" N;        
     "ARRAY" A, AUX;        
     "BEGIN" "INTEGER" J, K, KK, KJ, LOW, UP; "REAL" R, EPSNORM;    
         "REAL" "PROCEDURE" VECVEC(L, U, SHIFT, A, B); "CODE" 34010;
        
         R:= 0; KK:= 0;     
         "FOR" K:= 1 "STEP" 1 "UNTIL" N "DO"    
         "BEGIN" KK:= KK + K; "IF" A[KK] > R "THEN" R:= A[KK] "END";
         EPSNORM:= AUX[2] * R; KK:= 0;
         "FOR" K:= 1 "STEP" 1 "UNTIL" N "DO"    
         "BEGIN" KK:= KK + K; LOW:= KK - K + 1; UP:= KK - 1;        
             R:= A[KK] - VECVEC(LOW, UP, 0, A, A);        
             "IF" R <= EPSNORM "THEN" 
             "BEGIN" AUX[3]:= K - 1; "GOTO" END "END";    
             A[KK]:= R:= SQRT(R); KJ:= KK + K;  
             "FOR" J:= K + 1 "STEP" 1 "UNTIL" N "DO"      
             "BEGIN" A[KJ]:= (A[KJ] - 
                     VECVEC(LOW, UP, KJ - KK, A, A)) / R; 
                 KJ:= KJ + J
             "END"
         "END";   
         AUX[3]:= N;        
     END:         
     "END" CHLDEC1;         
         "EOP"    
1SECTION: 3.1.1.1.1.2.2       (FEBRUARY 1979)                     PAGE 1      
        
        
        
 CONTRIBUTORS: S.P.N. VAN KAMPEN, J. KOK.       
        
        
 INSTITUTE: MATHEMATICAL CENTRE.      
        
        
 RECEIVED: 731015.
        
        
 BRIEF DESCRIPTION: 
        
     THIS SECTION CONTAINS TWO PROCEDURES:      
     CHLDETERM2, FOR THE CALCULATION OF THE DETERMINANT OF A
     SYMMETRIC POSITIVE DEFINITE MATRIX WHOSE CHOLESKY MATRIX IS  GIVEN       
     IN THE UPPER TRIANGLE OF A TWO-DIMENSIONAL ARRAY;    
     CHLDETERM1, FOR THE CALCULATION OF THE DETERMINANT OF A
     SYMMETRIC POSITIVE DEFINITE MATRIX WHOSE CHOLESKY MATRIX IS  GIVEN       
     COLUMNWISE IN A ONE-DIMENSIONAL ARRAY.     
        
        
 KEYWORDS:
        
     DETERMINANT, 
     POSITIVE DEFINITE SYMMETRIC MATRIX,
     CHOLESKY DECOMPOSITION.
        
        
 SUBSECTION: CHLDETERM2.    
        
        
 CALLING SEQUENCE:
        
     THE HEADING OF THE PROCEDURE IS: 
     "REAL" "PROCEDURE" CHLDETERM2(A, N); "VALUE" N; "INTEGER" N;   
     "ARRAY" A; "CODE" 34312; 
        
     CHLDETERM2  := THE  DETERMINANT  OF THE  SYMMETRIC  POSITIVE   
     DEFINITE MATRIX OF WHICH THE CHOLESKY MATRIX IS STORED IN A;   
        
     THE MEANING OF THE FORMAL PARAMETERS IS:   
     A:      <ARRAY IDENTIFIER>;      
             "ARRAY" A[1:N,1:N];      
             ENTRY: THE  UPPER-TRIANGULAR PART OF THE  CHOLESKY MATRIX
                    AS  PRODUCED BY CHLDEC2 (SECTION  3.1.1.1.1.2.1) OR       
                    CHLDECSOL2 (SECTION  3.1.1.1.1.2.3)  MUST  BE GIVEN       
                    IN THE UPPER TRIANGLE OF A; 
             EXIT:  THE CONTENTS OF A ARE NOT CHANGED;    
     N:      <ARITHMETIC EXPRESSION>; 
             THE ORDER OF THE MATRIX. 
1SECTION: 3.1.1.1.1.2.2       (FEBRUARY 1979)                     PAGE 2      
        
        
        
 PROCEDURES USED: NONE.     
        
        
 REQUIRED CENTRAL MEMORY: NO EXTRA ARRAYS ARE DECLARED.   
        
        
 RUNNING TIME: PROPORTIONAL TO N.     
        
        
 METHOD AND PERFORMANCE:    
        
     THE PROCEDURE CHLDETERM2  SHOULD BE CALLED AFTER A  SUCCESSFUL CALL      
     OF CHLDEC2 OR CHLDECSOL2, I.E. IF AUX[3] = N;
     CHLDETERM2  SHOULD  NOT  BE CALLED  IF OVERFLOW IS TO BE  EXPECTED.      
        
        
 EXAMPLE OF USE:  
        
     SEE EXAMPLE OF USE OF CHLDECINV2 (SECTION 3.1.1.1.1.2.4)       
        
        
 SUBSECTION: CHLDETERM1.    
        
        
 CALLING SEQUENCE:
        
     THE HEADING OF THE PROCEDURE IS: 
     "REAL" "PROCEDURE" CHLDETERM1(A, N); "VALUE" N; "INTEGER" N;   
     "ARRAY" A; "CODE" 34313; 
        
     CHLDETERM1  := THE  DETERMINANT  OF THE  SYMMETRIC  POSITIVE   
     DEFINITE MATRIX OF WHICH THE CHOLESKY MATRIX IS STORED IN A;   
        
     THE MEANING OF THE FORMAL PARAMETERS IS:   
     A:      <ARRAY IDENTIFIER>;      
             "ARRAY" A[1 : (N + 1) * N // 2];   
             ENTRY: THE  UPPER-TRIANGULAR PART OF THE  CHOLESKY MATRIX
                    AS  PRODUCED BY CHLDEC1 (SECTION  3.1.1.1.1.2.1) OR       
                    CHLDECSOL1 (SECTION  3.1.1.1.1.2.3)  MUST  BE GIVEN       
                    COLUMNWISE IN ARRAY A;      
             EXIT:  THE CONTENTS OF A ARE NOT CHANGED;    
     N:      <ARITHMETIC EXPRESSION>; 
             THE ORDER OF THE MATRIX. 
1SECTION: 3.1.1.1.1.2.2       (FEBRUARY 1979)                     PAGE 3      
        
        
        
 PROCEDURES USED: NONE.     
        
        
 REQUIRED CENTRAL MEMORY: NO EXTRA ARRAYS ARE DECLARED.   
        
        
 RUNNING TIME: PROPORTIONAL TO N.     
        
        
 METHOD AND PERFORMANCE:    
        
     THE PROCEDURE CHLDETERM1  SHOULD BE CALLED AFTER A  SUCCESSFUL CALL      
     OF CHLDEC1 OR CHLDECSOL1, I.E. IF AUX[3] = N;
     CHLDETERM1  SHOULD  NOT  BE  CALLED  IF OVERFLOW IS TO BE EXPECTED.      
        
        
        
 EXAMPLE OF USE:  
        
     SEE EXAMPLE OF USE OF CHLDECINV1, SECTION 3.1.1.1.1.2.4.       
        
        
 SOURCE TEXT(S) : 
        
 "CODE" 34312;    
     "REAL" "PROCEDURE" CHLDETERM2(A, N); "VALUE" N; "INTEGER" N;   
     "ARRAY" A;   
     "BEGIN" "INTEGER" K; "REAL" D;   
         D:= 1;   
         "FOR" K:= 1 "STEP" 1 "UNTIL" N "DO" D:= A[K,K] * D;
         CHLDETERM2:= D * D 
     "END" CHLDETERM2;      
         "EOP"    
        
 "CODE" 34313;    
     "REAL" "PROCEDURE" CHLDETERM1(A, N); "VALUE" N; "INTEGER" N;   
     "ARRAY" A;   
     "BEGIN" "INTEGER" K, KK; "REAL" D; 
         D:= 1; KK:= 0;     
         "FOR" K:= 1 "STEP" 1 "UNTIL" N "DO"    
         "BEGIN" KK:= KK + K; D:= A[KK] * D "END";
         CHLDETERM1:= D * D 
     "END" CHLDETERM1;      
         "EOP"    
1SECTION: 3.1.1.1.1.2.3       (MAY 1974)                          PAGE 1      
        
        
        
 AUTHOR: T.J. DEKKER.       
        
        
 CONTRIBUTORS: S.P.N. VAN KAMPEN, J. KOK.       
        
        
 INSTITUTE: MATHEMATICAL CENTRE.      
        
        
 RECEIVED: 731015.
        
        
 BRIEF DESCRIPTION:         
        
     THIS SECTION CONTAINS FOUR PROCEDURES:     
     A) CHLSOL2, FOR THE SOLUTION OF A SYSTEM OF LINEAR EQUATIONS IF THE      
     COEFFICIENT MATRIX HAS BEEN DECOMPOSED BY A CALL OF THE PROCEDURE        
     CHLDEC2, SECTION 3.1.1.1.1.2.1., OR CHLDECSOL2;      
     B) CHLSOL1, FOR THE SOLUTION OF A SYSTEM OF LINEAR EQUATIONS IF THE      
     COEFFICIENT MATRIX HAS BEEN DECOMPOSED BY A CALL OF THE PROCEDURE        
     CHLDEC1, SECTION 3.1.1.1.1.2.1., OR CHLDECSOL1;      
     C) CHLDECSOL2, FOR THE SOLUTION OF A SYSTEM OF LINEAR EQUATIONS BY       
     CHOLESKY'S SQUARE ROOT METHOD;   
     THE COEFFICIENT MATRIX HAS TO BE SYMMETRIC POSITIVE DEFINITE AND         
     MUST BE GIVEN IN THE UPPER TRIANGLE OF A TWO-DIMENSIONAL ARRAY;
     D) CHLDECSOL1, FOR THE SOLUTION OF A SYSTEM OF LINEAR EQUATIONS BY       
     CHOLESKY'S SQUARE ROOT METHOD;   
     THE COEFFICIENT MATRIX HAS TO BE SYMMETRIC POSITIVE DEFINITE AND         
     MUST BE GIVEN COLUMNWISE IN A ONE-DIMENSIONAL ARRAY. 
        
        
 KEYWORDS:        
        
     LINEAR EQUATIONS,      
     POSITIVE DEFINITE SYMMETRIC MATRIX,        
     CHOLESKY DECOMPOSITION.
        
        
1SECTION: 3.1.1.1.1.2.3       (MAY 1974)                          PAGE 2      
        
        
        
 SUBSECTION: CHLSOL2.       
        
 CALLING SEQUENCE:
        
     THE HEADING OF THE PROCEDURE IS: 
     "PROCEDURE" CHLSOL2(A, N, B); "VALUE" N; "INTEGER" N; "ARRAY" A, B;      
        
     THE MEANING OF THE FORMAL PARAMETERS IS:   
     A:      <ARRAY IDENTIFIER>;      
             "ARRAY" A[1:N,1:N];      
             ENTRY: THE  UPPER-TRIANGULAR PART OF THE  CHOLESKY MATRIX        
                    AS  PRODUCED BY CHLDEC2, SECTION  3.1.1.1.1.2.1., OR      
                    CHLDECSOL2  (THIS  SECTION), MUST  BE  GIVEN IN  THE      
                    UPPER TRIANGLE OF A;        
             EXIT:  THE CONTENTS OF A ARE NOT CHANGED;    
     N:      <ARITHMETIC EXPRESSION>; 
             THE ORDER OF THE MATRIX; 
     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: 
        
     MATVEC = CP34011,      
     TAMVEC = CP34012.      
        
 RUNNING TIME: ROUGHLY PROPORTIONAL TO N SQUARED.         
        
 LANGUAGE: ALGOL 60.        
        
 METHOD AND PERFORMANCE:    
        
     THE PROCEDURE CHLSOL2 CALCULATES THE SOLUTION OF A SYSTEM OF LINEAR      
     EQUATIONS, PROVIDED THAT THE COEFFICIENT MATRIX HAS BEEN DECOMPOSED      
     BY A SUCCESSFUL CALL OF CHLDEC2 OR CHLDECSOL2;       
     THE SOLUTION IS OBTAINED BY CARRYING OUT THE FORWARD AND BACK  
     SUBSTITUTION WITH THE CHOLESKY MATRIX AND THE RIGHT HAND SIDE. 
     THE RIGHT HAND SIDE IS OVERWRITTEN BY THE SOLUTION BUT THE     
     ELEMENTS OF THE CHOLESKY MATRIX ARE NOT CHANGED, THUS SEVERAL  
     SYSTEMS OF LINEAR EQUATIONS WITH THE SAME COEFFICIENT MATRIX BUT         
     DIFFERENT RIGHT HAND SIDES CAN BE SOLVED BY SUCCESSIVE CALLS OF
     CHLSOL2. SEE ALSO REF [1].       
        
 REFERENCES:      
        
     [1]. T.J. DEKKER.      
          ALGOL 60 PROCEDURES IN NUMERICAL ALGEBRA, PART 1.         
          MC TRACT 22, 1968, MATH. CENTR., AMSTERDAM.     
        
 EXAMPLE OF USE:  
        
     SEE EXAMPLE OF USE OF CHLINV2, SECTION 3.1.1.1.1.2.4.
1SECTION: 3.1.1.1.1.2.3       (DECEMBER 1975)                     PAGE 3      
        
        
        
 SUBSECTION: CHLSOL1.       
        
        
 CALLING SEQUENCE:
        
     THE HEADING OF THE PROCEDURE IS: 
     "PROCEDURE" CHLSOL1(A, N, B); "VALUE" N; "INTEGER" N; "ARRAY" A, B;      
        
     THE MEANING OF THE FORMAL PARAMETERS IS:   
     A:      <ARRAY IDENTIFIER>;      
             "ARRAY" A[1 : (N + 1) * N // 2];   
             ENTRY: THE  UPPER-TRIANGULAR PART OF THE  CHOLESKY MATRIX        
                    AS  PRODUCED BY CHLDEC1, SECTION  3.1.1.1.1.2.1., OR      
                    CHLDECSOL1 (THIS SECTION), MUST  BE GIVEN COLUMNWISE      
                    IN ARRAY A;       
             EXIT:  THE CONTENTS OF A ARE NOT CHANGED;    
     N:      <ARITHMETIC EXPRESSION>; 
             THE ORDER OF THE MATRIX; 
     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,      
     SEQVEC = CP34016.      
        
        
 RUNNING TIME: ROUGHLY PROPORTIONAL TO N SQUARED.         
        
        
 LANGUAGE: ALGOL 60.        
        
        
 METHOD AND PERFORMANCE:    
        
     THE PROCEDURE CHLSOL1 CALCULATES THE SOLUTION OF A SYSTEM OF LINEAR      
     EQUATIONS, PROVIDED THAT THE COEFFICIENT MATRIX HAS BEEN DECOMPOSED      
     BY A SUCCESSFUL CALL OF CHLDEC1 OR CHLDECSOL1;       
     SEVERAL SYSTEMS WITH THE SAME COEFFICIENT MATRIX BUT DIFFERENT 
     RIGHT HAND SIDES CAN BE SOLVED BY SUCCESSIVE CALLS OF CHLSOL1. 
     SEE ALSO METHOD AND PERFORMANCE OF CHLSOL2 (THIS SECTION).     
        
        
 EXAMPLE OF USE:  
        
     SEE EXAMPLE OF USE OF CHLINV1, SECTION 3.1.1.1.1.2.4.
        
        
1SECTION: 3.1.1.1.1.2.3       (DECEMBER 1975)                     PAGE 4      
        
        
        
 SUBSECTION: CHLDECSOL2.    
        
        
 CALLING SEQUENCE:
        
     THE HEADING OF THE PROCEDURE IS: 
     "PROCEDURE" CHLDECSOL2(A, N, AUX, B); "VALUE" N; "INTEGER" N;  
     "ARRAY" A, AUX, B;     
        
     THE MEANING OF THE FORMAL PARAMETERS IS:   
     A:      <ARRAY IDENTIFIER>;      
             "ARRAY" A[1:N,1:N];      
             ENTRY: THE UPPER  TRIANGLE OF THE  POSITIVE DEFINITE MATRIX      
                    MUST BE GIVEN IN THE UPPER-TRIANGULAR PART OF A (THE      
                    ELEMENTS A[I,J], I <= J);   
             EXIT:  THE  CHOLESKY  DECOMPOSITION   OF  THE  MATRIX  IS        
                    DELIVERED IN THE UPPER TRIANGLE OF A; 
     N:      <ARITHMETIC EXPRESSION>; 
             THE ORDER OF THE MATRIX; 
     AUX:    <ARRAY IDENTIFIER>;      
             "ARRAY" AUX[2:3];        
             ENTRY: AUX[2]: A  RELATIVE TOLERANCE  USED TO  CONTROL  THE      
                    CALCULATION OF THE DIAGONAL ELEMENTS; (SEE METHOD         
                    AND PERFORMANCE OF CHLDEC2, SECTION 3.1.1.1.1.2.1);       
             NORMAL EXIT:   AUX[3]:= N;         
             ABNORMAL EXIT: IF THE  DECOMPOSITION CANNOT BE  CARRIED OUT      
                    BECAUSE  THE  MATRIX  IS (NUMERICALLY) NOT  POSITIVE      
                    DEFINITE,  AUX[3]:= K - 1, WHERE K IS THE LAST STAGE      
                    NUMBER. 
     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: 
     CHLDEC2 = CP34310,     
     CHLSOL2 = CP34390.     
        
 RUNNING TIME: ROUGHLY PROPORTIONAL TO N CUBED. 
        
 LANGUAGE: ALGOL 60.        
        
 METHOD AND PERFORMANCE:    
        
     THE PROCEDURE CHLDECSOL2 SOLVES A SYSTEM OF LINEAR EQUATIONS WITH        
     A SYMMETRIC POSITIVE DEFINITE COEFFICIENT MATRIX BY CALLING    
     CHLDEC2, SECTION 3.1.1.1.1.2.1., AND, IF THIS CALL WAS         
     SUCCESSFUL, CHLSOL2 (THIS SECTION).        
     SEE ALSO CHLDEC2, SECTION 3.1.1.1.1.2.1.   
        
 EXAMPLE OF USE:  
        
     SEE EXAMPLE OF USE OF CHLDECINV2, SECTION 3.1.1.1.1.2.4.       
1SECTION: 3.1.1.1.1.2.3       (DECEMBER 1975)                     PAGE 5      
        
        
        
 SUBSECTION: CHLDECSOL1.    
        
 CALLING SEQUENCE:
        
     THE HEADING OF THE PROCEDURE IS: 
     "PROCEDURE" CHLDECSOL1(A, N, AUX, B); "VALUE" N; "INTEGER" N;  
     "ARRAY" A, AUX, B;     
        
     THE MEANING OF THE FORMAL PARAMETERS IS:   
     A:      <ARRAY IDENTIFIER>;      
             "ARRAY" A[1 : (N + 1) * N // 2];   
             ENTRY: THE UPPER-TRIANGULAR PART OF THE  POSITIVE  DEFINITE      
                    SYMMETRIC MATRIX MUST BE GIVEN COLUMNWISE IN ARRAY A      
                    (THE (I,J)-TH ELEMENT OF THE MATRIX MUST BE GIVEN IN      
                    A[(J - 1) * J // 2 + I] FOR 1 <= I <= J <= N);  
             EXIT:  THE  CHOLESKY  DECOMPOSITION   OF  THE  MATRIX  IS        
                    DELIVERED COLUMNWISE IN A.  
     N:      <ARITHMETIC EXPRESSION>; 
             THE ORDER OF THE MATRIX; 
     AUX:    <ARRAY IDENTIFIER>;      
             "ARRAY" AUX[2:3];        
             ENTRY: AUX[2]: A  RELATIVE TOLERANCE  USED TO  CONTROL  THE      
                    CALCULATION OF THE DIAGONAL ELEMENTS; 
             NORMAL EXIT:   AUX[3]:= N;         
             ABNORMAL EXIT: IF THE  DECOMPOSITION CANNOT BE  CARRIED OUT      
                    BECAUSE  THE  MATRIX  IS (NUMERICALLY) NOT  POSITIVE      
                    DEFINITE,  AUX[3]:= K - 1, WHERE K IS THE LAST STAGE      
                    NUMBER. 
     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: 
        
     CHLDEC1 = CP34311,     
     CHLSOL1 = CP34391.     
        
 RUNNING TIME: ROUGHLY PROPORTIONAL TO N CUBED. 
        
 LANGUAGE: ALGOL 60.        
        
 METHOD AND PERFORMANCE:    
     THE PROCEDURE CHLDECSOL1 SOLVES A SYSTEM OF LINEAR EQUATIONS WITH        
     A SYMMETRIC POSITIVE DEFINITE COEFFICIENT MATRIX BY CALLING    
     CHLDEC1, SECTION 3.1.1.1.1.2.1., AND, IF THIS CALL WAS         
     SUCCESSFUL, CHLSOL1 (THIS SECTION).        
     THE UPPER TRIANGLE OF THE COEFFICIENT MATRIX MUST BE STORED COLUMN-      
     WISE IN A ONE-DIMENSIONAL ARRAY. 
     SEE ALSO CHLDEC1, SECTION 3.1.1.1.1.2.1.   
        
 EXAMPLE OF USE:  
     SEE EXAMPLE OF USE OF CHLDECINV1, SECTION 3.1.1.1.1.2.4.       
1SECTION: 3.1.1.1.1.2.3       (MAY 1974)                          PAGE 6      
        
        
        
 SOURCE TEXT(S) : 
        
 "CODE" 34390;    
     "PROCEDURE" CHLSOL2(A, N, B); "VALUE" N; "INTEGER" N; "ARRAY" A, B;      
     "BEGIN" "INTEGER" I;   
         "REAL" "PROCEDURE" MATVEC(L, U, I, A, B); "CODE" 34011;    
         "REAL" "PROCEDURE" TAMVEC(L, U, I, A, B); "CODE" 34012;    
        
         "FOR" I:= 1 "STEP" 1 "UNTIL" N "DO"    
         B[I]:= (B[I] - TAMVEC(1, I - 1, I, A, B)) / A[I,I];        
         "FOR" I:= N "STEP" - 1 "UNTIL" 1 "DO"  
         B[I]:= (B[I] - MATVEC(I + 1, N, I, A, B)) / A[I,I]         
     "END" CHLSOL2;         
         "EOP"    
 "CODE" 34391;    
     "PROCEDURE" CHLSOL1(A, N, B); "VALUE" N; "INTEGER" N; "ARRAY" A, B;      
     "BEGIN" "INTEGER" I, II;         
         "REAL" "PROCEDURE" VECVEC(L, U, SHIFT, A, B); "CODE" 34010;
         "REAL" "PROCEDURE" SEQVEC(L, U, I1, SHIFT, A, B); "CODE" 34016;      
        
         II:= 0;  
         "FOR" I:= 1 "STEP" 1 "UNTIL" N "DO"    
         "BEGIN" II:= II + I;         
             B[I]:= (B[I] - VECVEC(1, I - 1, II - I, B, A)) / A[II] 
         "END";   
         "FOR" I:= N "STEP" - 1 "UNTIL" 1 "DO"  
         "BEGIN" B[I]:= (B[I] -       
                 SEQVEC(I + 1, N, II + I, 0, A, B)) / A[II];        
             II:= II - I    
         "END"    
     "END" CHLSOL1;         
         "EOP"    
 "CODE" 34392;    
     "PROCEDURE" CHLDECSOL2(A, N, AUX, B); "VALUE" N; "INTEGER" N;  
     "ARRAY" A, AUX, B;     
     "BEGIN"      
         "PROCEDURE" CHLDEC2(A, N, AUX); "CODE" 34310;    
         "PROCEDURE" CHLSOL2(A, N, B); "CODE" 34390;      
        
         CHLDEC2(A, N, AUX);
         "IF" AUX[3] = N "THEN" CHLSOL2(A, N, B)
     "END" CHLDECSOL2;      
         "EOP"    
 "CODE" 34393;    
     "PROCEDURE" CHLDECSOL1(A, N, AUX, B); "VALUE" N; "INTEGER" N;  
     "ARRAY" A, AUX, B;     
     "BEGIN"      
         "PROCEDURE" CHLDEC1(A, N, AUX); "CODE" 34311;    
         "PROCEDURE" CHLSOL1(A, N, B); "CODE" 34391;      
        
         CHLDEC1(A, N, AUX);
         "IF" AUX[3] = N "THEN" CHLSOL1(A, N, B)
     "END" CHLDECSOL1;      
         "EOP"