NUMAL Section 3.1.1.1.1.1.1

BEGIN SECTION : 3.1.1.1.1.1.1 (February, 1979)

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.


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.

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).


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;

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).


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).


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).

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).


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.

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;

    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).

SOURCE TEXT(S):
"CODE" 34300;

"CODE" 34231;
"CODE" 34240;
"CODE" 34241;
"CODE" 34242;
"CODE" 34252;