NUMAL Section 3.1.1.1.1.1.4
BEGIN SECTION : 3.1.1.1.1.1.4 (December, 1975)
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.
SUBSECTION: INV .
CALLING SEQUENCE:
THE HEADING OF THIS PROCEDURE IS:
"PROCEDURE" INV(A, N, P); "VALUE" N;
"INTEGER" N; "ARRAY" A; "INTEGER" "ARRAY" P;
"CODE" 34053;
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).
SUBSECTION: DECINV .
CALLING SEQUENCE:
THE HEADING OF THIS PROCEDURE IS:
"PROCEDURE" DECINV(A, N, AUX); "VALUE" N;
"INTEGER" N; "ARRAY" A, AUX;
"CODE" 34302;
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.
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" 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),/),/")");
"BEGIN" "INTEGER" I, J;
"FOR" I:= 1 "STEP" 1 "UNTIL" 4 "DO"
"FOR" J:= 1 "STEP" 1 "UNTIL" 4 "DO" INREAL(70, A[I,J])
"END" IN LIST;
AUX[2]:= "-14;
OUTPUT(71,"("/,"("INPUT:")",///")");
OUTLIST(71, LAYOUT, LIST);
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
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;
"CODE" 34235;
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).
SUBSECTION: GSSINV .
CALLING SEQUENCE:
THE HEADING OF THIS PROCEDURE IS:
"PROCEDURE" GSSINV(A, N, AUX); "VALUE" N;
"INTEGER" N; "ARRAY" A, AUX;
"CODE" 34236;
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.
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" 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),/),/")");
"BEGIN" "INTEGER" I, J;
"FOR" I:= 1 "STEP" 1 "UNTIL" 4 "DO"
"FOR" J:= 1 "STEP" 1 "UNTIL" 4 "DO" INREAL(70, A[I,J])
"END" IN LIST;
AUX[2]:= "-14; AUX[4]:= 8;
OUTPUT(71,"("/,"("INPUT:")",///")");
OUTLIST(71, LAYOUT, LIST);
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
SUBSECTION: GSSINVERB.
CALLING SEQUENCE:
THE HEADING OF THIS PROCEDURE IS:
"PROCEDURE" GSSINVERB(A, N, AUX); "VALUE" N;
"INTEGER" N; "ARRAY" A, AUX;
"CODE" 34244;
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.
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" 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),/),/")");
"BEGIN" "INTEGER" I, J;
"FOR" I:= 1 "STEP" 1 "UNTIL" 4 "DO"
"FOR" J:= 1 "STEP" 1 "UNTIL" 4 "DO" INREAL(70, A[I,J])
"END" IN LIST;
AUX[0]:= AUX[2]:= AUX[6]:= "-14;
OUTPUT(71,"("/,"("INPUT:")",///")");
OUTLIST(71, LAYOUT, LIST);
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"
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).
SOURCE TEXT(S):
"CODE" 34053;
"CODE" 34302;
"CODE" 34235;
"CODE" 34236;
"CODE" 34244;