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;