NUMAL Section 3.1.1.1.1.1.5

BEGIN SECTION : 3.1.1.1.1.1.5 (February, 1979)

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.


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.

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


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;

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

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" 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


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;

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

    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;

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

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" 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

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;
"CODE" 34251;

"CODE" 34253;
"CODE" 34254;