1SECTION: 3.1.1.1.1.2.4 (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) CHLINV2, FOR THE INVERSION OF A SYMMETRIC POSITIVE DEFINITE MATRIX, IF THE MATRIX HAS BEEN DECOMPOSED BY A CALL OF THE PROCEDURE CHLDEC2, SECTION 3.1.1.1.1.2.1., OR CHLDECSOL2, SECTION 3.1.1.1.1.2.3.; B) CHLINV1, FOR THE INVERSION OF A SYMMETRIC POSITIVE DEFINITE MATRIX, IF THE MATRIX HAS BEEN DECOMPOSED BY A CALL OF THE PROCEDURE CHLDEC1, SECTION 3.1.1.1.1.2.1., OR CHLDECSOL1, SECTION 3.1.1.1.1.2.3.; C) CHLDECINV2, FOR THE INVERSION OF A MATRIX 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) CHLDECINV1, FOR THE INVERSION OF A MATRIX 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: MATRIX INVERSION, POSITIVE DEFINITE SYMMETRIC MATRIX, CHOLESKY DECOMPOSITION. 1SECTION: 3.1.1.1.1.2.4 (MAY 1974) PAGE 2 SUBSECTION: CHLINV2. CALLING SEQUENCE: THE HEADING OF THE PROCEDURE IS: "PROCEDURE" CHLINV2(A, N); "VALUE" N; "INTEGER" N; "ARRAY" A; THE MEANING OF THE FORMAL PARAMETERS IS: A: ; "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 UPPER-TRIANGULAR PART OF THE INVERSE MATRIX IS DELIVERED IN THE UPPER TRIANGLE OF A; N: ; THE ORDER OF THE MATRIX. PROCEDURES USED: MATVEC = CP34011, TAMVEC = CP34012, DUPVECROW = CP31031. REQUIRED CENTRAL MEMORY: EXECUTION FIELD LENGTH: N. RUNNING TIME: ROUGHLY PROPORTIONAL TO N CUBED. LANGUAGE: ALGOL 60. METHOD AND PERFORMANCE: THE PROCEDURE CHLINV2 CALCULATES THE INVERSE OF A MATRIX, PROVIDED THAT THE MATRIX HAS BEEN DECOMPOSED BY A SUCCESSFUL CALL OF CHLDEC2 OR CHLDECSOL2; THE INVERSE, X, OF U'U, WHERE U IS THE CHOLESKY MATRIX, IS OBTAINED FROM THE CONDITIONS THAT X BE SYMMETRIC AND UX BE A LOWER-TRIANGULAR MATRIX WHOSE MAIN DIAGONAL ELEMENTS ARE THE RECIPROCALS OF THE DIAGONAL ELEMENTS OF U. HEREWITH THE UPPER- TRIANGULAR ELEMENTS OF X ARE CALCULATED BY BACK SUBSTITUTION. THE UPPER TRIANGLE OF THE INVERSE MATRIX IS DELIVERED IN THE UPPER TRIANGLE OF THE GIVEN ARRAY. SEE ALSO REF[1]. REFERENCES: [1]. T.J. DEKKER. ALGOL 60 PROCEDURES IN NUMERICAL ALGEBRA, PART 1. MC TRACT 22, 1968, MATH. CENTR., AMSTERDAM. 1SECTION: 3.1.1.1.1.2.4 (MAY 1974) PAGE 3 EXAMPLE OF USE: THE SYMMETRIC POSITIVE DEFINITE COEFFICIENT MATRIX (THE PASCAL MATRIX OF ORDER 4) OF THE SYSTEM OF EQUATIONS X1 + X2 + X3 + X4 = 2 X1 + 2 * X2 + 3 * X3 + 4 * X4 = 4 X1 + 3 * X2 + 6 * X3 + 10 * X4 = 8 X1 + 4 * X2 + 10 * X3 + 20 * X4 = 16 IS STORED IN THE TWO-DIMENSIONAL ARRAY PASCAL2. THE INVERSE OF THE COEFFICIENT MATRIX AND THE SOLUTION OF THE LINEAR SYSTEM ARE CALCULATED BY THE FOLLOWING PROGRAM: "BEGIN" "COMMENT" TEST CHLDEC2, CHLSOL2 AND CHLINV2; "INTEGER" I, J; "ARRAY" PASCAL2[1:4,1:4], B[1:4], AUX[2:3]; "PROCEDURE" CHLDEC2(A, N, AUX); "CODE" 34310; "PROCEDURE" CHLSOL2(A, N, B); "CODE" 34390; "PROCEDURE" CHLINV2(A, N); "CODE" 34400; "FOR" J:= 1 "STEP" 1 "UNTIL" 4 "DO" "BEGIN" PASCAL2[1,J]:= 1; "FOR" I:= 2 "STEP" 1 "UNTIL" J "DO" PASCAL2[I,J]:= "IF" I = J "THEN" PASCAL2[I-1,J] * 2 "ELSE" PASCAL2[I,J-1] + PASCAL2[I-1,J]; B[J]:= 2 ** J "END"; AUX[2]:= "-11; CHLDEC2(PASCAL2, 4, AUX); "IF" AUX[3] = 4 "THEN" "BEGIN" CHLSOL2(PASCAL2, 4, B); CHLINV2(PASCAL2, 4) "END" "ELSE" OUTPUT(61, "(""("MATRIX NOT POSITIVE DEFINITE")", /")"); OUTPUT(61, "("4B")"); OUTPUT(61, "(""("SOLUTION WITH CHLDEC2 AND CHLSOL2:")", /")"); "FOR" I:= 1 "STEP" 1 "UNTIL" 4 "DO" OUTPUT(61, "("4B+D.5D")", B[I]); OUTPUT(61, "("//, 4B")"); OUTPUT(61, "(""("INVERSE MATRIX WITH CHLINV2:")", /, 4B")"); "FOR" I:= 1 "STEP" 1 "UNTIL" 4 "DO" "BEGIN" "FOR" J:= 1 "STEP" 1 "UNTIL" 4 "DO" "IF" J < I "THEN" OUTPUT(61, "("12B")") "ELSE" OUTPUT(61, "("+ZD.5D3B")", PASCAL2[I,J]); OUTPUT(61, "("/, 4B")") "END" "END" THIS PROGRAM DELIVERS: SOLUTION WITH CHLDEC2 AND CHLSOL2: +0.00000 +4.00000 -4.00000 +2.00000 INVERSE MATRIX WITH CHLINV2: +4.00000 -6.00000 +4.00000 -1.00000 +14.00000 -11.00000 +3.00000 +10.00000 -3.00000 +1.00000 1SECTION: 3.1.1.1.1.2.4 (MAY 1974) PAGE 4 SUBSECTION: CHLINV1. CALLING SEQUENCE: THE HEADING OF THE PROCEDURE IS: "PROCEDURE" CHLINV1(A, N); "VALUE" N; "INTEGER" N; "ARRAY" A; THE MEANING OF THE FORMAL PARAMETERS IS: A: ; "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 UPPER-TRIANGULAR PART OF THE INVERSE MATRIX IS DELIVERED COLUMNWISE IN ARRAY A; N: ; THE ORDER OF THE MATRIX. PROCEDURES USED: SEQVEC = CP34016, SYMMATVEC = CP34018. REQUIRED CENTRAL MEMORY: EXECUTION FIELD LENGTH: N. RUNNING TIME: ROUGHLY PROPORTIONAL TO N CUBED. LANGUAGE: ALGOL 60. METHOD AND PERFORMANCE: THE PROCEDURE CHLINV1 CALCULATES THE INVERSE OF A MATRIX, PROVIDED THAT THE MATRIX HAS BEEN DECOMPOSED BY A SUCCESSFUL CALL OF CHLDEC1 OR CHLDECSOL1; THE UPPER TRIANGLE OF THE INVERSE MATRIX IS DELIVERED COLUMNWISE IN THE ONE-DIMENSIONAL ARRAY. SEE ALSO METHOD AND PERFORMANCE OF CHLINV2 (THIS SECTION). EXAMPLE OF USE: THE SYMMETRIC POSITIVE DEFINITE COEFFICIENT MATRIX (THE PASCAL MATRIX OF ORDER 4) OF THE SYSTEM OF EQUATIONS X1 + X2 + X3 + X4 = 2 X1 + 2 * X2 + 3 * X3 + 4 * X4 = 4 X1 + 3 * X2 + 6 * X3 + 10 * X4 = 8 X1 + 4 * X2 + 10 * X3 + 20 * X4 = 16 IS STORED IN THE ONE-DIMENSIONAL ARRAY PASCAL1. THE INVERSE OF THE COEFFICIENT MATRIX AND THE SOLUTION OF THE LINEAR SYSTEM ARE CALCULATED BY THE FOLLOWING PROGRAM: 1SECTION: 3.1.1.1.1.2.4 (MAY 1974) PAGE 5 "BEGIN" "COMMENT" TEST CHLDEC1, CHLSOL1 AND CHLINV1; "INTEGER" I, J, JJ; "ARRAY" PASCAL1[1:(4 + 1) * 4 // 2], B[1:4], AUX[2:3]; "PROCEDURE" CHLDEC1(A, N, AUX); "CODE" 34311; "PROCEDURE" CHLSOL1(A, N, B); "CODE" 34391; "PROCEDURE" CHLINV1(A, N); "CODE" 34401; JJ:= 1; "FOR" J:= 1 "STEP" 1 "UNTIL" 4 "DO" "BEGIN" PASCAL1[JJ]:= 1; "FOR" I:= 2 "STEP" 1 "UNTIL" J "DO" PASCAL1[JJ + I - 1]:= "IF" I = J "THEN" PASCAL1[JJ + I - 2] * 2 "ELSE" PASCAL1[JJ + I - 2] + PASCAL1[JJ + I - J]; B[J]:= 2 ** J; JJ:= JJ + J "END"; AUX[2]:= "-11; CHLDEC1(PASCAL1, 4, AUX); "IF" AUX[3] = 4 "THEN" "BEGIN" CHLSOL1(PASCAL1, 4, B); CHLINV1(PASCAL1, 4) "END" "ELSE" OUTPUT(61, "(""("MATRIX NOT POSITIVE DEFINITE")", /")"); OUTPUT(61,"("4B")"); OUTPUT(61, "(""("SOLUTION WITH CHLDEC1 AND CHLSOL1:")", /")"); "FOR" I:= 1 "STEP" 1 "UNTIL" 4 "DO" OUTPUT(61, "("4B+D.5D")", B[I]); OUTPUT(61, "("2/, 4B")"); OUTPUT(61, "(""("INVERSE MATRIX WITH CHLINV1:")", /, 4B")"); "FOR" I:= 1 "STEP" 1 "UNTIL" 4 "DO" "BEGIN" "FOR" J:= 1 "STEP" 1 "UNTIL" 4 "DO" "IF" J < I "THEN" OUTPUT(61, "("12B")") "ELSE" OUTPUT(61, "("+ZD.5D3B")", PASCAL1[(J - 1) * J // 2 + I]); OUTPUT(61, "("/, 4B")") "END" "END" THIS PROGRAM DELIVERS: SOLUTION WITH CHLDEC1 AND CHLSOL1: +0.00000 +4.00000 -4.00000 +2.00000 INVERSE MATRIX WITH CHLINV1: +4.00000 -6.00000 +4.00000 -1.00000 +14.00000 -11.00000 +3.00000 +10.00000 -3.00000 +1.00000 1SECTION: 3.1.1.1.1.2.4 (MAY 1974) PAGE 6 SUBSECTION: CHLDECINV2. CALLING SEQUENCE: THE HEADING OF THE PROCEDURE IS: "PROCEDURE" CHLDECINV2(A, N, AUX); "VALUE" N; "INTEGER" N; "ARRAY" A, AUX; THE MEANING OF THE FORMAL PARAMETERS IS: A: ; "ARRAY" A[1:N,1:N]; ENTRY: THE UPPER TRIANGLE OF THE POSITIVE DEFINITE MATRIX MUST BE GIVEN IN THE UPPER TRIANGLE OF A (THE ELEMENTS A[I,J], I <= J); EXIT: THE UPPER-TRIANGULAR PART OF THE INVERSE MATRIX IS DELIVERED IN THE UPPER TRIANGLE OF A. N: ; THE ORDER OF THE MATRIX; AUX: ; "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: CHLDEC2 = CP34310, CHLINV2 = CP34400. RUNNING TIME: ROUGHLY PROPORTIONAL TO N CUBED. LANGUAGE: ALGOL 60. METHOD AND PERFORMANCE: THE PROCEDURE CHLDECINV2 CALCULATES THE INVERSE OF A SYMMETRIC POSITIVE DEFINITE MATRIX BY CALLING CHLDEC2 AND, IF THIS CALL WAS SUCCESSFUL, CHLINV2. THE UPPER TRIANGLE OF THE INVERSE MATRIX IS DELIVERED IN THE UPPER TRIANGLE OF THE GIVEN ARRAY. SEE ALSO METHOD AND PERFORMANCE OF CHLINV2 (THIS SECTION) AND CHLDEC2, SECTION 3.1.1.1.1.2.1. 1SECTION: 3.1.1.1.1.2.4 (MAY 1974) PAGE 7 EXAMPLE OF USE: THE SYMMETRIC POSITIVE DEFINITE COEFFICIENT MATRIX (THE PASCAL MATRIX OF ORDER 4) OF THE SYSTEM OF EQUATIONS X1 + X2 + X3 + X4 = 2 X1 + 2 * X2 + 3 * X3 + 4 * X4 = 4 X1 + 3 * X2 + 6 * X3 + 10 * X4 = 8 X1 + 4 * X2 + 10 * X3 + 20 * X4 = 16 IS STORED IN THE TWO-DIMENSIONAL ARRAY PASCAL2. THE DETERMINANT AND THE INVERSE OF THE COEFFICIENT MATRIX AND THE SOLUTION OF THE LINEAR SYSTEM ARE CALCULATED BY THE FOLLOWING PROGRAM: "BEGIN" "COMMENT" TEST CHLDECSOL2, CHLDETERM2 AND CHLDECINV2; "INTEGER" I, J; "ARRAY" PASCAL2[1:4,1:4], B[1:4], AUX[2:3]; "REAL" DETERMINANT; "PROCEDURE" CHLDECSOL2(A, N, AUX, B); "CODE" 34392; "REAL" "PROCEDURE" CHLDETERM2(A, N); "CODE" 34312; "PROCEDURE" CHLDECINV2(A, N, AUX); "CODE" 34402; "FOR" J:= 1 "STEP" 1 "UNTIL" 4 "DO" "BEGIN" PASCAL2[1,J]:= 1; "FOR" I:= 2 "STEP" 1 "UNTIL" J "DO" PASCAL2[I,J]:= "IF" I = J "THEN" PASCAL2[I-1,J] * 2 "ELSE" PASCAL2[I,J-1] + PASCAL2[I-1,J]; B[J]:= 2 ** J "END"; AUX[2]:= "-11; CHLDECSOL2(PASCAL2, 4, AUX, B); "IF" AUX[3] = 4 "THEN" DETERMINANT:= CHLDETERM2(PASCAL2, 4) "ELSE" OUTPUT(61, "(""("MATRIX NOT POSITIVE DEFINITE")", /")"); OUTPUT(61, "("4B")"); OUTPUT(61, "(""("SOLUTION WITH CHLDECSOL2:")", /")"); "FOR" I:= 1 "STEP" 1 "UNTIL" 4 "DO" OUTPUT(61, "("4B+D.5D")", B[I]); OUTPUT(61, "("//, 4B, "("DETERMINANT WITH CHLDETERM2: ")", +D.5D, 2/, 4B")", DETERMINANT); "FOR" J:= 1 "STEP" 1 "UNTIL" 4 "DO" "BEGIN" PASCAL2[1,J]:= 1; "FOR" I:= 2 "STEP" 1 "UNTIL" J "DO" PASCAL2[I,J]:= "IF" I = J "THEN" PASCAL2[I-1,J] * 2 "ELSE" PASCAL2[I,J-1] + PASCAL2[I-1,J] "END"; CHLDECINV2(PASCAL2, 4, AUX); OUTPUT(61, "(""("INVERSE MATRIX WITH CHLDECINV2:")", /, 4B")"); "FOR" I:= 1 "STEP" 1 "UNTIL" 4 "DO" "BEGIN" "FOR" J:= 1 "STEP" 1 "UNTIL" 4 "DO" "IF" J < I "THEN" OUTPUT(61, "("12B")") "ELSE" OUTPUT(61, "("+ZD.5D3B")", PASCAL2[I,J]); OUTPUT(61, "("/, 4B")") "END" "END" 1SECTION: 3.1.1.1.1.2.4 (DECEMBER 1975) PAGE 8 THIS PROGRAM DELIVERS: SOLUTION WITH CHLDECSOL2: +0.00000 +4.00000 -4.00000 +2.00000 DETERMINANT WITH CHLDETERM2: +1.00000 INVERSE MATRIX WITH CHLDECINV2: +4.00000 -6.00000 +4.00000 -1.00000 +14.00000 -11.00000 +3.00000 +10.00000 -3.00000 +1.00000 SUBSECTION: CHLDECINV1. CALLING SEQUENCE: THE HEADING OF THE PROCEDURE IS: "PROCEDURE" CHLDECINV1(A, N, AUX); "VALUE" N; "INTEGER" N; "ARRAY" A, AUX; THE MEANING OF THE FORMAL PARAMETERS IS: A: ; "ARRAY" A[1:(N + 1) * N // 2]; ENTRY: THE UPPER-TRIANGULAR PART OF THE SYMMETRIC POSITIVE DEFINITE 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 UPPER-TRIANGULAR PART OF THE INVERSE MATRIX IS DELIVERED COLUMNWISE IN ARRAY A; N: ; THE ORDER OF THE MATRIX; AUX: ; "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. 1SECTION: 3.1.1.1.1.2.4 (MAY 1974) PAGE 9 PROCEDURES USED: CHLDEC1 = CP34311, CHLINV1 = CP34401. RUNNING TIME: ROUGHLY PROPORTIONAL TO N CUBED. LANGUAGE: ALGOL 60. METHOD AND PERFORMANCE: THE PROCEDURE CHLDECINV1 CALCULATES THE INVERSE OF A SYMMETRIC POSITIVE DEFINITE MATRIX BY CALLING CHLDEC1 AND, IF THIS CALL WAS SUCCESSFUL, CHLINV1. THE UPPER TRIANGLE OF THE INVERSE MATRIX IS DELIVERED COLUMNWISE IN THE GIVEN ONE-DIMENSIONAL ARRAY. SEE ALSO METHOD AND PERFORMANCE OF CHLINV2, (THIS SECTION) AND CHLDEC1, SECTION 3.1.1.1.1.2.1. EXAMPLE OF USE: THE SYMMETRIC POSITIVE DEFINITE COEFFICIENT MATRIX (THE PASCAL MATRIX OF ORDER 4) OF THE SYSTEM OF EQUATIONS X1 + X2 + X3 + X4 = 2 X1 + 2 * X2 + 3 * X3 + 4 * X4 = 4 X1 + 3 * X2 + 6 * X3 + 10 * X4 = 8 X1 + 4 * X2 + 10 * X3 + 20 * X4 = 16 IS STORED IN THE ONE-DIMENSIONAL ARRAY PASCAL1. THE DETERMINANT AND THE INVERSE OF THE COEFFICIENT MATRIX AND THE SOLUTION OF THE LINEAR SYSTEM ARE CALCULATED BY THE FOLLOWING PROGRAM: 1SECTION: 3.1.1.1.1.2.4 (DECEMBER 1975) PAGE 10 "BEGIN" "COMMENT" TEST CHLDECSOL1, CHLDETERM1 AND CHLDECINV1; "INTEGER" I, J, JJ; "ARRAY" PASCAL1[1:(4 + 1) * 4 // 2], B[1:4], AUX[2:3]; "REAL" DETERMINANT; "PROCEDURE" CHLDECSOL1(A, N, AUX, B); "CODE" 34393; "REAL" "PROCEDURE" CHLDETERM1(A, N); "CODE" 34313; "PROCEDURE" CHLDECINV1(A, N, AUX); "CODE" 34403; JJ:= 1; "FOR" J:= 1 "STEP" 1 "UNTIL" 4 "DO" "BEGIN" PASCAL1[JJ]:= 1; "FOR" I:= 2 "STEP" 1 "UNTIL" J "DO" PASCAL1[JJ + I - 1]:= "IF" I = J "THEN" PASCAL1[JJ + I - 2] * 2 "ELSE" PASCAL1[JJ + I - 2] + PASCAL1[JJ + I - J]; B[J]:= 2 ** J; JJ:= JJ + J "END"; AUX[2]:= "-11; CHLDECSOL1(PASCAL1, 4, AUX, B); "IF" AUX[3] = 4 "THEN" DETERMINANT:= CHLDETERM1(PASCAL1, 4) "ELSE" OUTPUT(61, "(""("MATRIX NOT POSITIVE DEFINITE")", /")"); OUTPUT(61, "("4B")"); OUTPUT(61, "(""("SOLUTION WITH CHLDECSOL1:")", /")"); "FOR" I:= 1 "STEP" 1 "UNTIL" 4 "DO" OUTPUT(61, "("4B+D.5D")", B[I]); OUTPUT(61, "("//, 4B, "("DETERMINANT WITH CHLDETERM1: ")", +D.5D, 2/, 4B")", DETERMINANT); JJ:= 1; "FOR" J:= 1 "STEP" 1 "UNTIL" 4 "DO" "BEGIN" PASCAL1[JJ]:= 1; "FOR" I:= 2 "STEP" 1 "UNTIL" J "DO" PASCAL1[JJ + I - 1]:= "IF" I = J "THEN" PASCAL1[JJ + I - 2] * 2 "ELSE" PASCAL1[JJ + I - 2] + PASCAL1[JJ + I - J]; JJ:= JJ + J "END"; CHLDECINV1(PASCAL1, 4, AUX); OUTPUT(61, "(""("INVERSE MATRIX WITH CHLDECINV1:")", /, 4B")"); "FOR" I:= 1 "STEP" 1 "UNTIL" 4 "DO" "BEGIN" "FOR" J:= 1 "STEP" 1 "UNTIL" 4 "DO" "IF" J < I "THEN" OUTPUT(61, "("12B")") "ELSE" OUTPUT(61, "("+ZD.5D3B")", PASCAL1[(J - 1) * J // 2 + I]); OUTPUT(61, "("/, 4B")") "END" "END" 1SECTION: 3.1.1.1.1.2.4 (MAY 1974) PAGE 11 THIS PROGRAM DELIVERS: SOLUTION WITH CHLDECSOL1: +0.00000 +4.00000 -4.00000 +2.00000 DETERMINANT WITH CHLDETERM1: +1.00000 INVERSE MATRIX WITH CHLDECINV1: +4.00000 -6.00000 +4.00000 -1.00000 +14.00000 -11.00000 +3.00000 +10.00000 -3.00000 +1.00000 0SOURCE TEXT(S) : "CODE" 34400; "PROCEDURE" CHLINV2(A, N); "VALUE" N; "INTEGER" N; "ARRAY" A; "BEGIN" "REAL" R; "INTEGER" I, J, I1; "ARRAY" U[1:N]; "PROCEDURE" DUPVECROW(L, U, I, A, B); "CODE" 31031; "REAL" "PROCEDURE" MATVEC(L, U, I, A, B); "CODE" 34011; "REAL" "PROCEDURE" TAMVEC(L, U, I, A, B); "CODE" 34012; "FOR" I:= N "STEP" - 1 "UNTIL" 1 "DO" "BEGIN" R:= 1 / A[I,I]; I1:= I + 1; DUPVECROW(I1, N, I, U, A); "FOR" J:= N "STEP" - 1 "UNTIL" I1 "DO" A[I,J]:= - (TAMVEC(I1, J, J, A, U) + MATVEC(J + 1, N, J, A, U)) * R; A[I,I]:= (R - MATVEC(I1, N, I, A, U)) * R "END" "END" CHLINV2; "EOP" "CODE" 34401; "PROCEDURE" CHLINV1(A, N); "VALUE" N; "INTEGER" N; "ARRAY" A; "BEGIN" "INTEGER" I, II, I1, J, IJ; "REAL" R; "ARRAY" U[1:N]; "REAL" "PROCEDURE" SEQVEC(L, U, I1, SHIFT, A, B); "CODE" 34016; "REAL" "PROCEDURE" SYMMATVEC(L, U, I, A, B); "CODE" 34018; II:= (N + 1) * N // 2; "FOR" I:= N "STEP" - 1 "UNTIL" 1 "DO" "BEGIN" R:= 1 / A[II]; I1:= I + 1; IJ:= II + I; "FOR" J:= I1 "STEP" 1 "UNTIL" N "DO" "BEGIN" U[J]:= A[IJ]; IJ:= IJ + J "END"; "FOR" J:= N "STEP" - 1 "UNTIL" I1 "DO" "BEGIN" IJ:= IJ - J; A[IJ]:= -SYMMATVEC(I1, N, J, A, U) * R "END"; A[II]:= (R - SEQVEC(I1, N, II + I, 0, A, U)) * R; II:= II - I "END" "END" CHLINV1; "EOP" 1SECTION: 3.1.1.1.1.2.4 (MAY 1974) PAGE 12 "CODE" 34402; "PROCEDURE" CHLDECINV2(A, N, AUX); "VALUE" N; "INTEGER" N; "ARRAY" A, AUX; "BEGIN" "PROCEDURE" CHLDEC2(A, N, AUX); "CODE" 34310; "PROCEDURE" CHLINV2(A, N); "CODE" 34400; CHLDEC2(A, N, AUX); "IF" AUX[3] = N "THEN" CHLINV2(A, N) "END" CHLDECINV2; "EOP" "CODE" 34403; "PROCEDURE" CHLDECINV1(A, N, AUX); "VALUE" N; "INTEGER" N; "ARRAY" A, AUX; "BEGIN" "PROCEDURE" CHLDEC1(A, N, AUX); "CODE" 34311; "PROCEDURE" CHLINV1(A, N); "CODE" 34401; CHLDEC1(A, N, AUX); "IF" AUX[3] = N "THEN" CHLINV1(A, N) "END" CHLDECINV1; "EOP" 1SECTION 3.1.1.2.1.1 (MAY 1974) PAGE 1 AUTHOR : T.J. DEKKER. CONTRIBUTOR : J. KOK. INSTITUTE : MATHEMATICAL CENTRE. RECEIVED : 731015. BRIEF DESCRIPTION : THIS SECTION CONTAINS TWO PROCEDURES : A) LSQORTDEC, FOR THE HOUSEHOLDER TRIANGULARIZATION WITH COLUMN INTERCHANGES OF THE COEFFICIENT MATRIX OF A LINEAR LEAST SQUARES PROBLEM; B) LSQDGLINV, FOR THE CALCULATION OF THE DIAGONAL ELEMENTS OF THE INVERSE OF M'M, WHERE M IS THE COEFFICIENT MATRIX OF A LINEAR LEAST SQUARES PROBLEM. KEY WORDS : LINEAR LEAST SQUARES PROBLEM, HOUSEHOLDER TRIANGULARIZATION. SUBSECTION : LSQORTDEC. CALLING SEQUENCE: THE HEADING OF THE PROCEDURE IS : "PROCEDURE" LSQORTDEC(A, N, M, AUX, AID, CI); "VALUE" N, M; "INTEGER" N, M; "INTEGER""ARRAY" CI; "ARRAY" A, AUX, AID; THE MEANING OF THE FORMAL PARAMETERS IS : A : ; "ARRAY" A[1 : N,1 : M]; ENTRY :THE COEFFICIENT MATRIX OF THE LINEAR LEAST SQUARES PROBLEM; EXIT : IN THE UPPER TRIANGLE OF A (THE ELEMENTS A[I,J] WITH I < J) THE SUPERDIAGONAL ELEMENTS OF THE UPPER-TRIANGULAR MATRIX, PRODUCED BY THE HOUSEHOLDER TRANSFORMATION; IN THE OTHER PART OF THE COLUMNS OF A THE SIGNIFICANT ELEMENTS OF THE GENERATING VECTORS OF THE HOUSEHOLDER MATRICES USED FOR THE HOUSEHOLDER TRIANGULARIZATION; 1SECTION 3.1.1.2.1.1 (MAY 1974) PAGE 2 N : ; NUMBER OF ROWS OF THE MATRIX; M : ; NUMBER OF COLUMNS OF THE MATRIX (N >= M); AUX : ; "ARRAY" AUX[2 : 5]; ENTRY : AUX[2] CONTAINS A RELATIVE TOLERANCE USED FOR CALCULATING THE DIAGONAL ELEMENTS OF THE UPPER-TRIANGULAR MATRIX; EXIT : AUX[3] DELIVERS THE NUMBER OF THE DIAGONAL ELEMENTS OF THE UPPER-TRIANGULAR MATRIX WHICH ARE FOUND NOT NEGLIGIBLE; NORMAL EXIT AUX[3] = M; AUX[5] := THE MAXIMUM OF THE EUCLIDEAN NORMS OF THE COLUMNS OF THE GIVEN MATRIX; AID : ; "ARRAY" AID[1 : M]; NORMAL EXIT (AUX[3] = M) : AID CONTAINS THE DIAGONAL ELEMENTS OF THE UPPER-TRIANGULAR MATRIX PRODUCED BY THE HOUSEHOLDER TRIANGULARIZATION; CI : ; "INTEGER""ARRAY" CI[1 : M]; EXIT : CI CONTAINS THE PIVOTAL INDICES OF THE INTERCHANGES OF THE COLUMNS OF THE GIVEN MATRIX. PROCEDURES USED : TAMMAT = CP34014, ELMCOL = CP34023, ICHCOL = CP34031. REQUIRED CENTRAL MEMORY : EXECUTION FIELD LENGTH : AN ARRAY OF M ELEMENTS IS DECLARED. RUNNING TIME: (C1 * M + C2) * M * (N - M / 3); THE CONSTANTS C1 AND C2 DEPEND ON THE ARITHMETIC OF THE COMPUTER. LANGUAGE : ALGOL 60. 1SECTION 3.1.1.2.1.1 (MAY 1974) PAGE 3 METHOD AND PERFORMANCE : THE PROCEDURE LSQORTDEC IS A MODIFICATION OF THE PROCEDURE LSQDEC DUE TO T.J. DEKKER (SEE REF [1]), WHERE A DERIVATION IS GIVEN OF A SET OF PROCEDURES BY P. BUSINGER AND G.H. GOLUB (SEE REF [2]). THE METHOD IS HOUSEHOLDER TRIANGULARIZATION WITH COLUMN INTERCHANGES. LET M DENOTE THE GIVEN MATRIX. LSQORTDEC PRODUCES AN N-TH ORDER ORTHOGONAL MATRIX Q AND AN N * M UPPER-TRIANGULAR MATRIX R SUCH THAT R EQUALS QM WITH PERMUTED COLUMNS, Q IS THE PRODUCT OF AT MOST M HOUSEHOLDER MATRICES WHICH ARE REPRESENTED BY THEIR GENERATING VECTORS. M IS REDUCED TO R IN AT MOST M STAGES : AT THE K-TH STAGE THE K-TH COLUMN OF THE (ALREADY MODIFIED) MATRIX IS INTERCHANGED WITH THE COLUMN OF MAXIMUM EUCLIDEAN NORM (THE PIVOTAL COLUMN); THEN THE MATRIX IS MULTIPLIED WITH A HOUSEHOLDER MATRIX SUCH, THAT THE SUBDIAGONAL ELEMENTS OF THE K-TH COLUMN BECOME ZERO, WHILE THE FIRST K - 1 COLUMNS REMAIN UNCHANGED. THE PROCESS TERMINATES PREMATURELY, IF AT SOME STAGE THE EUCLIDEAN NORM OF THE PIVOTAL COLUMN IS LESS THAN SOME TOLERANCE, VIZ. A GIVEN TOLERANCE (AUX[2]) TIMES THE MAXIMUM OF THE EUCLIDEAN NORMS OF THE COLUMNS OF THE GIVEN MATRIX. LSQORTDEC DELIVERS THE SIGNIFICANT ELEMENTS OF THE GENERATING VECTOR OF THE K-TH HOUSEHOLDER MATRIX (THE FIRST K - 1 ELEMENTS OF THIS VECTOR BEING ZERO) IN THE LOWER TRIANGLE PART OF THE K-TH COLUMN OF THE ARRAY A (A[I,K] FOR I >= K). OF THE RESULTING UPPER-TRIANGULAR MATRIX THE DIAGONAL ELEMENTS ARE DELIVERED SEPARATELY IN AN ARRAY AID, AND THE REMAINING ELEMENTS IN THE SUPER-TRIANGULAR PART OF THE ARRAY A. FOR THE SOLUTION OF LEAST SQUARES PROBLEMS, ONLY CALLS WITH N >= M ARE USEFUL. REFERENCES : [1] DEKKER, T.J. : ALGOL 60 PROCEDURES IN NUMERICAL ALGEBRA, PART 1, MC TRACT 22, 1968, MATHEMATISCH CENTRUM, AMSTERDAM. [2] BUSINGER, P. AND G.H. GOLUB : LINEAR LEAST SQUARES SOLUTION BY HOUSEHOLDER TRANSFORMATIONS, NUM. MATH. 7 (1965), PP. 269 - 276. EXAMPLE OF USE : SEE EXAMPLE OF USE OF LSQSOL. 1SECTION 3.1.1.2.1.1 (MAY 1974) PAGE 4 SUBSECTION : LSQDGLINV. CALLING SEQUENCE : THE HEADING OF THE PROCEDURE IS : "PROCEDURE" LSQDGLINV(A, M, AID, CI, DIAG); "VALUE" M; "INTEGER" M; "INTEGER""ARRAY" CI; "ARRAY" A, AID, DIAG; THE MEANING OF THE FORMAL PARAMETERS IS : A, M, AID, CI : SEE CALLING SEQUENCE OF LSQORTDEC (THIS SECTION); THE CONTENTS OF A, AID AND CI SHOULD BE PRODUCED BY A SUCCESSFUL CALL OF LSQORTDEC (AUX[3] = M) . DIAG : ; "ARRAY" DIAG[1 : M]; EXIT : THE DIAGONAL ELEMENTS OF THE INVERSE OF M'M WHERE M IS THE MATRIX OF THE LINEAR LEAST SQUARES PROBLEM. PROCEDURES USED : VECVEC = CP34010, TAMVEC = CP34012. RUNNING TIME : (C3 * M + C4) * M * M; THE CONSTANTS C3 AND C4 DEPEND ON THE ARITHMETIC OF THE COMPUTER. LANGUAGE : ALGOL 60. METHOD AND PERFORMANCE : LSQDGLINV SHOULD BE CALLED AFTER A SUCCESSFUL CALL OF LSQORTDEC, I.E. IF AUX[3] = M. LSQDGLINV CALCULATES THE DIAGONAL ELEMENTS OF THE INVERSE OF M'M, WHERE M IS THE MATRIX OF A LINEAR LEAST SQUARES PROBLEM. THESE VALUES CAN BE USED FOR THE COMPUTATION OF THE STANDARD DEVIATIONS OF LEAST SQUARES SOLUTIONS. EXAMPLE OF USE : SEE EXAMPLE OF USE OF LSQSOL. 1SECTION 3.1.1.2.1.1 (MAY 1974) PAGE 5 SOURCE TEXT(S) : "CODE" 34134; "PROCEDURE" LSQORTDEC(A, N, M, AUX, AID, CI); "VALUE" N, M; "INTEGER" N, M; "ARRAY" A, AUX, AID; "INTEGER" "ARRAY" CI; "BEGIN" "INTEGER" J, K, KPIV; "REAL" BETA, SIGMA, NORM, W, EPS, AKK, AIDK; "ARRAY" SUM[1:M]; "REAL""PROCEDURE" TAMMAT(L, U, I, J, A, B);"CODE" 34014; "PROCEDURE" ELMCOL(L, U, I, J, A, B, X);"CODE" 34023; "PROCEDURE" ICHCOL(L, U, I, J, A);"CODE" 34031; NORM:= 0; AUX[3]:= M; "FOR" K:= 1 "STEP" 1 "UNTIL" M "DO" "BEGIN" W:= SUM[K]:= TAMMAT(1, N, K, K, A, A); "IF" W > NORM "THEN" NORM:= W "END"; W:= AUX[5]:= SQRT(NORM); EPS:= AUX[2] * W; "FOR" K:= 1 "STEP" 1 "UNTIL" M "DO" "BEGIN" SIGMA:= SUM[K]; KPIV:= K; "FOR" J:= K + 1 "STEP" 1 "UNTIL" M "DO" "IF" SUM[J] > SIGMA "THEN" "BEGIN" SIGMA:= SUM[J]; KPIV:= J "END"; "IF" KPIV ^= K "THEN" "BEGIN" SUM[KPIV]:= SUM[K]; ICHCOL(1, N, K, KPIV, A) "END"; CI[K]:= KPIV; AKK:= A[K,K]; SIGMA:= TAMMAT(K, N, K, K, A, A); W:= SQRT(SIGMA); AIDK:= AID[K]:= "IF" AKK < 0 "THEN" W "ELSE" - W; "IF" W < EPS "THEN" "BEGIN" AUX[3]:= K - 1; "GO TO" ENDDEC "END"; BETA:= 1 / (SIGMA - AKK * AIDK); A[K,K]:= AKK - AIDK; "FOR" J:= K + 1 "STEP" 1 "UNTIL" M "DO" "BEGIN" ELMCOL(K, N, J, K, A, A, - BETA * TAMMAT(K, N, K, J, A, A)); SUM[J]:= SUM[J] - A[K,J] ** 2 "END" "END" FOR K; ENDDEC: "END" LSQORTDEC; "EOP" 1SECTION 3.1.1.2.1.1 (MAY 1974) PAGE 6 "CODE" 34132; "PROCEDURE" LSQDGLINV(A, M, AID, CI, DIAG); "VALUE" M; "INTEGER" M; "ARRAY" A, AID, DIAG; "INTEGER" "ARRAY" CI; "BEGIN" "INTEGER" J, K, CIK; "REAL" W; "REAL""PROCEDURE" VECVEC(L, U, S, A, B);"CODE" 34010; "REAL""PROCEDURE" TAMVEC(L, U, I, A, B);"CODE" 34012; "FOR" K:= 1 "STEP" 1 "UNTIL" M "DO" "BEGIN" DIAG[K]:= 1 / AID[K]; "FOR" J:= K + 1 "STEP" 1 "UNTIL" M "DO" DIAG[J]:= - TAMVEC(K, J - 1, J, A, DIAG) / AID[J]; DIAG[K]:= VECVEC(K, M, 0, DIAG, DIAG) "END"; "FOR" K:= M "STEP" - 1 "UNTIL" 1 "DO" "BEGIN" CIK:= CI[K]; "IF" CIK ^= K "THEN" "BEGIN" W:= DIAG[K]; DIAG[K]:= DIAG[CIK]; DIAG[CIK]:= W "END" "END" "END" LSQDGLINV; "EOP" 1SECTION 3.1.1.2.1.2 (MAY 1974) PAGE 1 AUTHOR : T.J. DEKKER. CONTRIBUTOR : J. KOK. INSTITUTE : MATHEMATICAL CENTRE. RECEIVED : 731015. BRIEF DESCRIPTION : THIS SECTION CONTAINS TWO PROCEDURES : A) LSQSOL, FOR THE SOLUTION OF A LINEAR LEAST SQUARES PROBLEM IF THE COEFFICIENT MATRIX HAS BEEN DECOMPOSED BY LSQORTDEC (SECTION 3.1.1.2.1.1.); B) LSQORTDECSOL, FOR THE SOLUTION OF A LINEAR LEAST SQUARES PROBLEM BY HOUSEHOLDER TRIANGULARIZATION WITH COLUMN INTERCHANGES AND FOR THE CALCULATION OF THE DIAGONAL OF THE INVERSE OF M'M, WHERE M IS THE COEFFICIENT MATRIX. KEY WORDS : LINEAR LEAST SQUARES PROBLEM, HOUSEHOLDER TRIANGULARIZATION. SUBSECTION : LSQSOL. CALLING SEQUENCE: THE HEADING OF THE PROCEDURE IS : "PROCEDURE" LSQSOL(A, N, M, AID, CI, B); "VALUE" N, M; "INTEGER" N, M; "INTEGER""ARRAY" CI; "ARRAY" A, AID, B; THE MEANING OF THE FORMAL PARAMETERS IS : A, N, M, AID, CI : SEE CALLING SEQUENCE OF LSQORTDEC (SECTION 3.1.1.2.1.1.); THE CONTENTS OF THE ARRAYS A, AID AND CI SHOULD BE PRODUCED BY A SUCCESSFUL CALL OF LSQORTDEC, I.E. IF AUX[3] = M; B : ; "ARRAY" B[1 : N]; ENTRY : B CONTAINS THE RIGHT HAND SIDE OF A LINEAR LEAST SQUARES PROBLEM; EXIT : B[1 : M] CONTAINS THE SOLUTION OF THE PROBLEM; B[M + 1 : N] CONTAINS A VECTOR WITH EUCLIDEAN LENGTH EQUAL TO THE EUCLIDEAN LENGTH OF THE RESIDUE VECTOR. 1SECTION 3.1.1.2.1.2 (MAY 1974) PAGE 2 PROCEDURES USED : MATVEC = CP34011, TAMVEC = CP34012, ELMVECCOL = CP34021. RUNNING TIME: (C5 * M + C6) * N; THE CONSTANTS C5 AND C6 DEPEND UPON THE ELEMENTARY ARITHMETIC OF THE COMPUTER. LANGUAGE : ALGOL 60. METHOD AND PERFORMANCE : LSQSOL SHOULD BE CALLED AFTER A SUCCESSFUL CALL OF LSQORTDEC (SECTION 3.1.1.2.1.1.), I.E. IF AUX[3] = M. LSQSOL YIELDS THE LEAST SQUARES SOLUTION OF THE OVERDETERMINED SYSTEM WITH THE DECOMPOSED COEFFICIENT MATRIX IN ARRAY A AND THE RIGHT HAND SIDE IN ARRAY B. FIRST THE ORTHOGONAL TRANSFORMATION WITH THE HOUSEHOLDER MATRICES IS PERFORMED ON THE RIGHT HAND SIDE. NEXT THE SYSTEM OF THE FIRST M EQUATIONS AND WITH AN UPPER-TRIANGULAR COEFFICIENT MATRIX IS SOLVED BY BACK SUBSTITUTION, YIELDING A SOLUTION WITH M PERMUTED COMPONENTS DUE TO THE COLUMN INTERCHANGES OF THE TRIANGULARIZATION. FINALLY THE ORDER OF THE M COMPONENTS IS RESTORED. SEE ALSO METHOD AND PERFORMANCE OF LSQORTDEC (SECTION 3.1.1.2.1.1.). THE LEAST SQUARES SOLUTIONS OF SEVERAL OVERDETERMINED SYSTEMS WITH THE SAME COEFFICIENT MATRIX CAN BE SOLVED BY SUCCESSIVE CALLS OF LSQSOL WITH DIFFERENT RIGHT HAND SIDES. EXAMPLE OF USE : THE NEXT PROGRAM SOLVES THE SYSTEM - 2 * X1 + X2 = 0 - X1 + X2 = 1 X1 + X2 = 2 2 * X1 + X2 = 2 X1 + 2 * X2 = 3 1SECTION 3.1.1.2.1.2 (MAY 1974) PAGE 3 "BEGIN""COMMENT" 730912, TEST LSQORTDEC, LSQSOL, LSQDGLINV; "ARRAY" A, C[1 : 5,1 : 2], B, X[1 : 5], DIAG, AID[1 : 2], AUX[2 : 5]; "INTEGER""ARRAY" PIV[1 : 2]; "INTEGER" I, J; "REAL" H; "REAL""PROCEDURE" VECVEC(L, U, S, A, B);"CODE" 34010; "PROCEDURE" LSQORTDEC(A, N, M, AUX, AID, CI);"CODE" 34134; "PROCEDURE" LSQSOL(A, N, M, AID, CI, B);"CODE" 34131; "PROCEDURE" LSQDGLINV(A, M, AID, CI, DIAG);"CODE" 34132; "REAL""PROCEDURE" SUM(I, A, B, X); "VALUE" A, B; "INTEGER" I, A, B; "REAL" X; "BEGIN""REAL" S; S:= 0; "FOR" I:= A "STEP" 1 "UNTIL" B "DO" S:= S + X; SUM:= S "END" SUM; AUX[2]:= "-12; I:= J:= 1; "FOR" H:= - 2, - 1, 1, 2, 1, 1, 1, 1, 1, 2 "DO" "BEGIN" A[I,J]:= C[I,J]:= H; "IF" I < 5 "THEN" I:= I + 1 "ELSE" "BEGIN" I:= 1; J:= J + 1 "END" "END"; "FOR" H:= 0, 1, 2, 2, 3 "DO" "BEGIN" B[I]:= X[I]:= H; I:= I + 1 "END"; LSQORTDEC(A, 5, 2, AUX, AID, PIV); "IF" AUX[3] = 2 "THEN" "BEGIN" LSQSOL(A, 5, 2, AID, PIV, X); LSQDGLINV(A, 2, AID, PIV, DIAG); OUTPUT(61, "("/, "("AUX[2, 3, 5] = ")" +.4D"+DD5B, 3ZD5B, +.4D"+DD/, "("LSQ SOLUTION :")", 2(2B+.8D"+DD), / "("RESIDUE (DELIVERED) :")" +.8D"+DD/, "("RESIDUE (CHECKED) :")" +.8D"+DD/, "("DIAGONAL OF INVERSE M'M :")", 2(2B+.8D"+DD)")", AUX[2], AUX[3], AUX[5], X[1], X[2], SQRT(VECVEC(3, 5, 0, X, X)), SQRT(SUM(I, 1, 5, (B[I] - C[I,1] * X[1] - C[I,2] * X[2]) ** 2)), DIAG[1], DIAG[2]) "END" "END" "EOP" END OF PROGRAM DELIVERS : AUX[2, 3, 5] = +.1000"-11 2 +.3317"+01 LSQ SOLUTION : +.50000000"+00 +.12500000"+01 RESIDUE (DELIVERED) :+.50000000"+00 RESIDUE (CHECKED) :+.50000000"+00 DIAGONAL OF INVERSE M'M : +.95238095"-01 +.13095238"+00 1SECTION 3.1.1.2.1.2 (MAY 1974) PAGE 4 SUBSECTION : LSQORTDECSOL. CALLING SEQUENCE : THE HEADING OF THE PROCEDURE IS : "PROCEDURE" LSQORTDECSOL(A, N, M, AUX, DIAG, B); "VALUE" N, M; "INTEGER" N, M; "ARRAY" A, AUX, DIAG, B; THE MEANING OF THE FORMAL PARAMETERS IS : A : ; "ARRAY" A[1 : N,1 : M]; ENTRY : A CONTAINS THE COEFFICIENT MATRIX OF THE LINEAR LEAST SQUARES PROBLEM; EXIT : IN THE UPPER TRIANGLE OF A (THE ELEMENTS A[I,J] WITH I < J) THE SUPERDIAGONAL ELEMENTS OF THE UPPER-TRIANGULAR MATRIX, PRODUCED BY THE HOUSEHOLDER TRANSFORMATION; IN THE OTHER PART OF THE COLUMNS OF A THE SIGNIFICANT ELEMENTS OF THE GENERATING VECTORS OF THE HOUSEHOLDER MATRICES USED FOR THE HOUSEHOLDER TRIANGULARIZATION; N : ; NUMBER OF ROWS OF THE MATRIX; M : ; NUMBER OF COLUMNS OF THE MATRIX (N >= M); AUX : ; "ARRAY" AUX[2 : 5]; ENTRY : AUX[2] CONTAINS A RELATIVE TOLERANCE USED FOR CALCULATING THE DIAGONAL ELEMENTS OF THE UPPER-TRIANGULAR MATRIX; EXIT : AUX[3] DELIVERS THE NUMBER OF THE DIAGONAL ELEMENTS OF THE UPPER-TRIANGULAR MATRIX WHICH ARE FOUND NOT NEGLIGIBLE; NORMAL EXIT AUX[3] = M; AUX[5] := THE MAXIMUM OF THE EUCLIDEAN NORMS OF THE COLUMNS OF THE GIVEN MATRIX; DIAG : ; "ARRAY" DIAG[1 : M]; EXIT : THE DIAGONAL ELEMENTS OF THE INVERSE OF M'M WHERE M IS THE MATRIX OF THE LINEAR LEAST SQUARES PROBLEM; B : ; "ARRAY" B[1 : N]; ENTRY : B CONTAINS THE RIGHT HAND SIDE OF A LINEAR LEAST SQUARES PROBLEM; EXIT : B[1 : M] CONTAINS THE SOLUTION OF THE PROBLEM; B[M + 1 : N] CONTAINS A VECTOR WITH EUCLIDEAN LENGTH EQUAL TO THE EUCLIDEAN LENGTH OF THE RESIDUE VECTOR. PROCEDURES USED : LSQORTDEC = CP34134, LSQDGLINV = CP34132, LSQSOL = CP34131. 1SECTION 3.1.1.2.1.2 (MAY 1974) PAGE 5 REQUIRED CENTRAL MEMORY : EXECUTION FIELD LENGTH : AN INTEGER ARRAY AND A REAL ARRAY, BOTH OF M ELEMENTS, ARE DECLARED. RUNNING TIME : ROUGHLY PROPORTIONAL TO N * M ** 2. LANGUAGE : ALGOL 60. METHOD AND PERFORMANCE : LSQORTDECSOL SOLVES AN OVERDETERMINED SYSTEM OF N LINEAR EQUATIONS IN M UNKNOWNS BY CALLING LSQORTDEC AND, IF THIS CALL WAS SUCCESSFUL LSQDGLINV AND LSQSOL. LSQORTDECSOL DELIVERS THE LEAST SQUARES SOLUTION AND THE DIAGONAL OF THE INVERSE OF M'M, WHERE M IS THE COEFFICIENT MATRIX OF THE SYSTEM. SEE SECTION 3.1.1.2.1.1., AND LSQSOL (THIS SECTION). EXAMPLE OF USE : THE PROGRAM "BEGIN""COMMENT" 730914, TEST LSQORTDECSOL; "ARRAY" A, C[1 : 5,1 : 2], B, X[1 : 5], DIAG[1 : 2], AUX[2 : 5]; "INTEGER" I, J; "REAL" H; "REAL""PROCEDURE" VECVEC(L, U, S, A, B);"CODE" 34010; "PROCEDURE" LSQORTDECSOL(A, N, M, AUX, DIAG, B);"CODE" 34135; "REAL""PROCEDURE" SUM(I, A, B, X); "VALUE" A, B; "INTEGER" I, A, B; "REAL" X; "BEGIN""REAL" S; S:= 0; "FOR" I:= A "STEP" 1 "UNTIL" B "DO" S:= S + X; SUM:= S "END" SUM; AUX[2]:= "-12; I:= J:= 1; "FOR" H:= - 2, - 1, 1, 2, 1, 1, 1, 1, 1, 2 "DO" "BEGIN" A[I,J]:= C[I,J]:= H; "IF" I < 5 "THEN" I:= I + 1 "ELSE" "BEGIN" I:= 1; J:= J + 1 "END" "END"; "FOR" H:= 0, 1, 2, 2, 3 "DO" "BEGIN" B[I]:= X[I]:= H; I:= I + 1 "END"; LSQORTDECSOL(A, 5, 2, AUX, DIAG, X); "IF" AUX[3] = 2 "THEN" OUTPUT(61, "("/, "("AUX[2, 3, 5] = ")" +.4D"+DD5B, 3ZD5B, +.4D"+DD/, "("LSQ SOLUTION :")", 2(2B+.8D"+DD), / "("RESIDUE (DELIVERED) :")" +.8D"+DD/, "("RESIDUE (CHECKED) :")" +.8D"+DD/, "("DIAGONAL OF INVERSE M'M :")", 2(2B+.8D"+DD)")", AUX[2], AUX[3], AUX[5], X[1], X[2], SQRT(VECVEC(3, 5, 0, X, X)), SQRT(SUM(I, 1, 5, (B[I] - C[I,1] * X[1] - C[I,2] * X[2]) ** 2)), DIAG[1], DIAG[2]) "END" "EOP" END OF PROGRAM 1SECTION 3.1.1.2.1.2 (MAY 1974) PAGE 6 WHICH SOLVES THE PROBLEM OF THE EXAMPLE OF USE OF LSQSOL, DELIVERS : AUX[2, 3, 5] = +.1000"-11 2 +.3317"+01 LSQ SOLUTION : +.50000000"+00 +.12500000"+01 RESIDUE (DELIVERED) :+.50000000"+00 RESIDUE (CHECKED) :+.50000000"+00 DIAGONAL OF INVERSE M'M : +.95238095"-01 +.13095238"+00 SOURCE TEXT(S) : 0"CODE" 34131; "PROCEDURE" LSQSOL(A, N, M, AID, CI, B); "VALUE" N, M; "INTEGER" N, M; "ARRAY" A, AID, B; "INTEGER" "ARRAY" CI; "BEGIN" "INTEGER" K, CIK; "REAL" W; "REAL""PROCEDURE" MATVEC(L, U, I, A, B);"CODE" 34011; "REAL""PROCEDURE" TAMVEC(L, U, I, A, B);"CODE" 34012; "PROCEDURE" ELMVECCOL(L, U, I, A, B, X);"CODE" 34021; "FOR" K:= 1 "STEP" 1 "UNTIL" M "DO" ELMVECCOL(K, N, K, B, A, TAMVEC(K, N, K, A, B) / (AID[K] * A[K,K])); "FOR" K:= M "STEP" - 1 "UNTIL" 1 "DO" B[K]:= (B[K] - MATVEC (K + 1, M, K, A, B)) / AID[K]; "FOR" K:= M "STEP" - 1 "UNTIL" 1 "DO" "BEGIN" CIK:= CI[K]; "IF" CIK ^= K "THEN" "BEGIN" W:= B[K]; B[K]:= B[CIK]; B[CIK]:= W "END" "END" "END" LSQSOL; "EOP" 0"CODE" 34135; "PROCEDURE" LSQORTDECSOL(A, N, M, AUX, DIAG, B); "VALUE" N, M; "INTEGER" N, M; "ARRAY" A, AUX, DIAG, B; "BEGIN" "ARRAY" AID[1:M]; "INTEGER" "ARRAY" CI[1:M]; "PROCEDURE" LSQORTDEC(A, N, M, AUX, AID, CI);"CODE" 34134; "PROCEDURE" LSQDGLINV(A, M, AID, CI, DIAG);"CODE" 34132; "PROCEDURE" LSQSOL(A, N, M, AID, CI, B);"CODE" 34131; LSQORTDEC(A, N, M, AUX, AID, CI); "IF" AUX[3] = M "THEN" "BEGIN" LSQDGLINV(A, M, AID, CI, DIAG); LSQSOL(A, N, M, AID, CI, B) "END" "END" LSQORTDECSOL; "EOP" 1SECTION 3.1.1.3.1.1 (FEBRUARY 1979) PAGE 1 AUTHOR : D.T.WINTER INSTITUTE : MATHEMATICAL CENTRE RECEIVED : 731217 BRIEF DESCRIPTION : THIS SECTION CONTAINS 2 PROCEDURES FOR THE SOLUTION OF AN OVERDETERMINED SYSTEM OF LINEAR EQUATIONS: SOLSVDOVR SOLVES AN OVERDETERMINED SYSTEM OF LINEAR EQUATIONS , MULTIPLYING THE RIGHT-HAND SIDE BY THE PSEUDO-INVERSE OF THE GIVEN MATRIX; THE SINGULAR VALUES DECOMPOSITION SHOULD BE AVAILABLE. SOLOVR CALCULATES THE SINGULAR VALUES DECOMPOSITION AND SOLVES AN OVERDETERMINED SYSTEM OF LINEAR EQUATIONS. KEYWORDS : BEST LEAST-SQUARES SOLUTION SINGULAR VALUES PSEUDO-INVERSE SUBSECTION : SOLSVDOVR CALLING SEQUENCE : THE HEADING OF THE PROCEDURE IS : "PROCEDURE" SOLSVDOVR(U, VAL, V, M, N, X, EM); "VALUE" M,N; "INTEGER" M,N; "ARRAY" U, VAL, V, X, EM; "CODE" 34280; THE MEANING OF THE FORMAL PARAMETERS IS : U: ; "ARRAY" U[1:M,1:N]; ENTRY:THE MATRIX U IN THE SINGULAR VALUES DECOMPOSITION U*S*V'. VAL: ; "ARRAY" VAL[1:N]; ENTRY:THE SINGULAR VALUES. V: ; "ARRAY" V[1:N,1:N]; ENTRY:THE MATRIX V IN THE SINGULAR VALUES DECOMPOSITION. N: ; THE NUMBER OF UNKNOWNS. M: ; THE LENGTH OF THE RIGHT-HAND SIDE VECTOR, N SHOULD SATISFY N <= M. X: ; "ARRAY" X[1:M]; ENTRY: THE RIGHT-HAND SIDE VECTOR; EXIT: THE SOLUTION VECTOR IN X[1:N]. EM: ; "ARRAY" EM[6:6]; ENTRY: EM[6]: THE MINIMAL NON-NEGLECTABLE SINGULAR VALUE. 1SECTION 3.1.1.3.1.1 (FEBRUARY 1979) PAGE 2 PROCEDURES USED : MATVEC = CP34011 TAMVEC = CP34012 REQUIRED CENTRAL MEMORY : AN AUXILIARY ARRAY OF N REALS IS DECLARED. RUNNING TIME : ROUGHLY PROPORTIONAL TO (M + N) * N METHOD AND PERFORMANCE : THE SOLUTION IS FOUND IN THREE STEPS : 1. U' * X = X1 IS CALCULATED, 2. VAL+ * X1 = X2 IS CALCULATED, HERE VAL+ DENOTES THE DIAGONAL MATRIX OBTAINED FROM VAL BY SETTING VAL+[I,I] = 1/VAL[I] IF VAL[I] GREATER THAN OR EQUAL TO EM[6], AND 0 OTHERWISE, 3. THE SOLUTION V * X2 IS CALCULATED. SUBSECTION : SOLOVR CALLING SEQUENCE : THE HEADING OF THE PROCEDURE IS : "INTEGER" "PROCEDURE" SOLOVR(A, M, N, X, EM); "VALUE" M, N; "INTEGER" M, N; "ARRAY" A, X, EM; "CODE" 34281; SOLOVR:= THE NUMBER OF SINGULAR VALUES NOT FOUND, I.E. ZERO IF ALL SINGULAR VALUES ARE CALCULATED. THE MEANING OF THE FORMAL PAREMETERS IS : A: ; "ARRAY" A[1:M,1:N]; ENTRY: THE MATRIX OF THE SYSTEM; M: ; THE NUMBER OF ROWS OF A; N: ; THE NUMBER OF COLUMNS OF A, N <= M; X: ; "ARRAY" X[1:M]; ENTRY: THE RIGHT-HAND SIDE VECTOR; EXIT: THE SOLUTION VECTOR IN X[1:N]; EM: ; "ARRAY" EM[0:7]; ENTRY: EM[0]: THE MACHINE PRECISION; EM[2]: THE RELATIVE PRECISION OF THE SINGULAR VALUES; EM[4]: THE MAXIMAL NUMBER OF ITERATIONS TO BE PERFORMED IN THE SINGULAR VALUES DECOMPOSITION; EM[6]: THE MINIMAL NON-NEGLECTABLE SINGULAR VALUE; EXIT:EM[1]: THE INFINITY NORM OF THE MATRIX; EM[3]: THE MAXIMAL NEGLECTED SUPERDIAGONAL ELEMENT; EM[5]: THE NUMBER OF ITERATIONS PERFORMED IN THE SINGULAR VALUES DECOMPOSITION; EM[7]: THE NUMERICAL RANK OF THE MATRIX, I.E. THE NUMBER OF SINGULAR VALUES GREATER THAN OR EQUAL TO EM[6]. 1SECTION 3.1.1.3.1.1 (FEBRUARY 1979) PAGE 3 PROCEDURES USED : QRISNGVALDEC = CP34273 SOLSVDOVR = CP34280 REQUIRED CENTRAL MEMORY : AUXILIARY ARRAYS ARE DECLARED TO A TOTAL OF (N + 2) * N REALS RUNNING TIME : ROUGHLY PROPORTIONAL TO (M + N) * N * N METHOD AND PERFORMANCE : THE SOLUTION IS FOUND IN TWO STEPS : 1. THE SINGULAR VALUES DECOMPOSITION IS CALCULATED BY MEANS OF THE PROCEDURE QRISNGVALDEC (SECTION 3.5.1.2); 2. THE SOLUTION IS CALCULATED BY MEANS OF THE PROCEDURE SOLSVDOVR, (THIS SECTION); REFERENCES : WILKINSON, J.H. AND C.REINSCH HANDBOOK OF AUTOMATIC COMPUTATION, VOL. 2 (CONTRIBUTION I-10) LINEAR ALGEBRA HEIDELBERG (1971) EXAMPLE OF USE : FIRST A PROGRAM IS GIVEN, AND THEN THE RESULTS OF THIS PROGRAM : "BEGIN" "ARRAY" A[1:8,1:5], B[1:8], EM[0:7]; "INTEGER" I; "INTEGER" "PROCEDURE" SOLOVR(A, M, N, X, EM); "VALUE" M, N; "INTEGER" M, N; "ARRAY" A, X, EM; "CODE" 34281; A[1,1]:=22; A[1,2]:= A[2,3]:=10; A[1,3]:= A[7,1]:= A[8,5]:=2; A[1,4]:= A[3,5]:=3; A[1,5]:= A[2,2]:=7; A[2,1]:=14; A[2,5]:=8; A[2,4]:= A[8,3]:=0; A[3,1]:= A[3,3]:= A[6,5]:=-1; A[3,2]:=13; A[3,4]:=-11; A[4,1]:=-3; A[4,2]:= A[4,4]:= A[5,4]:= A[8,4]:=-2; A[4,3]:=13; A[4,5]:= A[5,5]:= A[8,1]:=4; A[5,1]:= A[6,1]:=9; A[5,2]:=8; A[5,3]:= A[6,2]:= A[7,5]:=1; A[6,3]:=-7; A[6,4]:= A[7,4]:= A[8,2]:=5; A[7,2]:=-6; A[7,3]:=6; B[1]:=-1; B[2]:=2; B[3]:= B[7]:=1; B[4]:=4; B[5]:= B[8]:=0; B[6]:=-3; EM[0]:="-14; EM[2]:="-12; EM[4]:=80; EM[6]:="-10; I:= SOLOVR(A, 8, 5, B, EM); OUTPUT(61, "("4B, "("NUMBER SINGULAR VALUES NOT FOUND : ")", 3ZD,/, 4B, "("NORM : ")", N,/, 4B, "("MAX NEGL SUBD ELEM : ")", N,/, 4B, "("NUMBER ITERATIONS : ")", 3ZD, /, 4B, "("RANK : ")", 3ZD, /")", I, EM[1], EM[3], EM[5], EM[7]); OUTPUT(61, "("/, 4B, "("SOLUTION VECTOR")",/,/, 5(4B, N, /)")", B[1], B[2], B[3], B[4], B[5]) "END" 1SECTION 3.1.1.3.1.1 (MAY 1974) PAGE 4 NUMBER SINGULAR VALUES NOT FOUND : 0 NORM : +4.4000000000000"+001 MAX NEGL SUBD ELEM : +4.3977072741076"-014 NUMBER ITERATIONS : 6 RANK : 3 SOLUTION VECTOR -8.3333333333334"-002 +1.0989227456287"-015 +2.5000000000000"-001 -8.3333333333332"-002 +8.3333333333334"-002 SOURCE TEXT(S): 0"CODE" 34280; "PROCEDURE" SOLSVDOVR(U, VAL, V, M, N, X, EM); "VALUE" M, N; "INTEGER" M, N; "ARRAY" U, VAL, V, X, EM; "BEGIN" "INTEGER" I; "REAL" MIN; "ARRAY" X1[1:N]; "REAL" "PROCEDURE" MATVEC(L, U, I, A, B); "VALUE" L, U, I; "INTEGER" L, U, I; "ARRAY" A, B; "CODE" 34011; "REAL" "PROCEDURE" TAMVEC(L, U, I, A, B); "VALUE" L, U, I; "INTEGER" L, U, I; "ARRAY" A, B; "CODE" 34012; MIN:= EM[6]; "FOR" I:= 1 "STEP" 1 "UNTIL" N "DO" X1[I]:= "IF" VAL[I] <= MIN "THEN" 0 "ELSE" TAMVEC(1, M, I, U, X) / VAL[I]; "FOR" I:= 1 "STEP" 1 "UNTIL" N "DO" X[I]:= MATVEC(1, N, I, V, X1) "END" SOLSVDOVR; "EOP" 0"CODE" 34281; "INTEGER" "PROCEDURE" SOLOVR(A, M, N, X, EM); "VALUE" M, N; "INTEGER" M, N; "ARRAY" A, X, EM; "BEGIN" "INTEGER" I; "ARRAY" VAL[1:N], V[1:N,1:N]; "INTEGER" "PROCEDURE" QRISNGVALDEC(A, M, N, VAL, V, EM); "VALUE" M, N; "INTEGER" M, N; "ARRAY" A, VAL, V, EM; "CODE" 34273; "PROCEDURE" SOLSVDOVR(U, VAL, V, M, N, X, EM); "VALUE" M, N; "INTEGER" M, N; "ARRAY" U, VAL, V, X, EM; "CODE" 34280; SOLOVR:= I:= QRISNGVALDEC(A, M, N, VAL, V, EM); "IF" I = 0 "THEN" SOLSVDOVR(A, VAL, V, M, N, X, EM) "END" SOLOVR; "EOP" 1SECTION 3.1.1.3.1.2 (DECEMBER 1975) PAGE 1 AUTHOR : D.T.WINTER INSTITUTE : MATHEMATICAL CENTRE RECEIVED : 731217 BRIEF DESCRIPTION : THIS SECTION CONTAINS 2 PROCEDURES FOR THE SOLUTION OF AN UNDERDETERMINED SYSTEM OF LINEAR EQUATIONS: SOLUND EXPECTS AS INPUT THE MATRIX OF THE SYSTEM OF EQUATIONS, CALCULATES THE SINGULAR VALUES DECOMPOSITION BY MEANS OF THE PROCEDURE QRISNGVALDEC, AND SOLVES THE SYSTEM BY MEANS OF THE PROCEDURE SOLSVDUND. SOLSVDUND ASSUMES THAT THE MATRIX IS ALREADY DECOMPOSED AND SOLVES THE SYSTEM OF EQUATIONS, MULTIPLYING THE RIGHT-HAND SIDE BY THE PSEUDO-INVERSE OF THE GIVEN MATRIX. KEYWORDS : BEST LEAST-SQUARES SOLUTION SINGULAR VALUES PSEUDO-INVERSE 1SECTION 3.1.1.3.1.2 (DECEMBER 1975) PAGE 2 SUBSECTION : SOLSVDUND CALLING SEQUENCE : THE HEADING OF THE PROCEDURE IS : "PROCEDURE" SOLSVDUND(U, VAL, V, M, N, X, EM); "VALUE" M, N; "INTEGER" M, N; "ARRAY" U, VAL, V, X, EM; THE MEANING OF THE FORMAL PARAMETERS IS : U: ; "ARRAY" U[1:M,1:N]; ENTRY:THE MATRIX U IN THE SINGULAR VALUES DECOMPOSITION V*S*U'. VAL: ; "ARRAY" VAL[1:N]; ENTRY:THE SINGULAR VALUES; V: ; "ARRAY" V[1:N,1:N]; ENTRY:THE MATRIX V IN THE SINGULAR VALUES DECOMPOSITION. N: ; THE LENGTH OF THE RIGHT-HAND SIDE VECTOR; M: ; THE NUMBER OF UNKNOWNS, N SHOULD SATISFY N <= M; X: ; "ARRAY" X[1:M]; ENTRY: THE RIGHT-HAND SIDE VECTOR IN X[1:N]; EXIT: THE SOLUTION VECTOR IN X[1:M]; EM: ; "ARRAY" EM[6:6]; ENTRY: EM[6]: THE MINIMAL NON-NEGLECTABLE SINGULAR VALUE. PROCEDURES USED : MATVEC = CP34011 TAMVEC = CP34012 REQUIRED CENTRAL MEMORY : AN AUXILIARY ARRAY OF N REALS IS DECLARED RUNNING TIME : ROUGHLY PROPORTIONAL TO (M + N) * N METHOD AND PERFORMANCE : THE SOLUTION IS FOUND IN THREE STEPS : 1. V' * X = X1 IS CALCULATED, 2. VAL+ * X1 = X2 IS CALCULATED, HERE VAL+ DENOTES THE DIAGONAL MATRIX OBTAINED FROM VAL BY SETTING VAL+[I,I] = 1/VAL[I] IF VAL[I] GREATER THAN OR EQUAL TO EM[6], AND 0 OTHERWISE, 3. THE SOLUTION U * X2 IS CALCULATED. LANGUAGE : ALGOL 60 1SECTION 3.1.1.3.1.2 (DECEMBER 1975) PAGE 3 SUBSECTION : SOLUND CALLING SEQUENCE : THE HEADING OF THE PROCEDURE IS : "INTEGER" "PROCEDURE" SOLUND(A, M, N, X, EM); "VALUE" M, N; "INTEGER" M, N; "ARRAY" A, X, EM; SOLUND:= THE NUMBER OF SINGULAR VALUES NOT FOUND, I.E. ZERO IF ALL SINGULAR VALUES ARE CALCULATED. THE MEANING OF THE FORMAL PAREMETERS IS : A: ; "ARRAY" A[1:M,1:N]; ENTRY: THE TRANSPOSE OF THE MATRIX; M: ; THE NUMBER OF ROWS OF A; N: ; THE NUMBER OF COLUMNS OF A, N <= M; X: ; "ARRAY" X[1:M]; ENTRY: THE RIGHT-HAND SIDE VECTOR IN X[1:N]; EXIT: THE SOLUTION VECTOR. EM: ; "ARRAY" EM[0:7]; ENTRY: EM[0]: THE MACHINE PRECISION; EM[2]: THE RELATIVE PRECISION FOR THE SINGULAR VALUES; EM[4]: THE MAXIMAL NUMBER OF ITERATIONS TO BE PERFORMED IN THE SINGULAR VALUES DECOMPOSITION; EM[6]: THE MINIMAL NON-NEGLECTABLE SINGULAR VALUE; EXIT: EM[1]: THE INFINITY NORM OF THE MATRIX; EM[3]: THE MAXIMAL NEGLECTED SUPERDIAGONAL ELEMENT; EM[5]: THE NUMBER OF ITERATIONS PERFORMED IN THE SINGULAR VALUES DECOMPOSITION; EM[7]: THE NUMERICAL RANK OF THE MATRIX, I.E. THE NUMBER OF SINGULAR VALUES GREATER THAN OR EQUAL TO EM[6]. PROCEDURES USED : QRISNGVALDEC = CP34273 SOLSVDUND = CP34282 REQUIRED CENTRAL MEMORY : AUXILIARY ARRAYS ARE DECLARED TO A TOTAL OF (N + 1) * N REALS METHOD AND PERFORMANCE : THE SOLUTION IS FOUND IN TWO STEPS : 1. THE SINGULAR VALUES DECOMPOSITION IS CALCULATED BY MEANS OF THE PROCEDURE QRISNGVALDEC; 2. THE SOLUTION IS CALCULATED BY MEANS OF THE PROCEDURE SOLSVDUND. 1SECTION 3.1.1.3.1.2 (MAY 1974) PAGE 4 RUNNING TIME : ROUGHLY PROPORTIONAL TO (M + N) * N * N LANGUAGE : ALGOL-60 REFERENCES : WILKINSON, J.H. AND C.REINSCH HANDBOOK OF AUTOMATIC COMPUTATION, VOL. 2 LINEAR ALGEBRA HEIDELBERG (1971) EXAMPLE OF USE : FIRST WE GIVE A PROGRAM, AND THAN THE RESULTS OF THIS PROGRAM : "BEGIN" "ARRAY" A[1:8,1:5], B[1:8], EM[0:7]; "INTEGER" I; "INTEGER" "PROCEDURE" SOLUND(A, M, N, X, EM); "VALUE" M, N; "INTEGER" M, N; "ARRAY" A, X, EM; "CODE" 34283; A[1,1]:=22; A[1,2]:= A[2,3]:=10; A[1,3]:= A[7,1]:= A[8,5]:=2; A[1,4]:= A[3,5]:=3; A[1,5]:= A[2,2]:=7; A[2,1]:=14; A[2,5]:=8; A[2,4]:= A[8,3]:=0; A[3,1]:= A[3,3]:= A[6,5]:=-1; A[3,2]:=13; A[3,4]:=-11; A[4,1]:=-3; A[4,2]:= A[4,4]:= A[5,4]:= A[8,4]:=-2; A[4,3]:=13; A[4,5]:= A[5,5]:= A[8,1]:=4; A[5,1]:= A[6,1]:=9; A[5,2]:=8; A[5,3]:= A[6,2]:= A[7,5]:=1; A[6,3]:=-7; A[6,4]:= A[7,4]:= A[8,2]:=5; A[7,2]:=-6; A[7,3]:=6; B[1]:=-1; B[2]:=2; B[3]:=1; B[4]:=4; B[5]:=0; EM[0]:="-14; EM[2]:="-12; EM[4]:=80; EM[6]:="-10; I:= SOLUND(A, 8, 5, B, EM); OUTPUT(61, "("4B, "("NUMBER SINGULAR VALUES NOT FOUND : ")", 3ZD,/, 4B, "("NORM : ")",N, /, 4B, "("MAX NEGL SUBD ELEM : ")", N,/, 4B, "("NUMBER ITERATIONS : ")", 3ZD, /, 4B, "("RANK : ")", 3ZD, /")", I, EM[1], EM[3], EM[5], EM[7]); OUTPUT(61, "("/, 4B, "("SOLUTION VECTOR")", /,/, 8(4B, N,/)")", B[1], B[2], B[3], B[4], B[5], B[6], B[7], B[8]) "END" NUMBER SINGULAR VALUES NOT FOUND : 0 NORM : +4.4000000000000"+001 MAX NEGL SUBD ELEM : +4.3977072741076"-014 NUMBER ITERATIONS : 6 RANK : 3 SOLUTION VECTOR +1.6410256410255"-002 +1.4807692307694"-002 -4.8397435897438"-002 +1.0000000000002"-002 -6.7948717948740"-003 +1.1602564102565"-002 +2.9999999999996"-002 -8.3974358974328"-003 1SECTION 3.1.1.3.1.2 (MAY 1974) PAGE 5 SOURCE TEXT(S): 0"CODE" 34282; "PROCEDURE" SOLSVDUND(U, VAL, V, M, N, X, EM); "VALUE" M, N; "INTEGER" M, N; "ARRAY" U, VAL, V, X, EM; "BEGIN" "INTEGER" I; "REAL" MIN; "ARRAY" X1[1:N]; "REAL" "PROCEDURE" MATVEC(L, U, I, A, B); "VALUE" L, U, I; "INTEGER" L, U, I; "ARRAY" A, B; "CODE" 34011; "REAL" "PROCEDURE" TAMVEC(L, U, I, A, B); "VALUE" L, U, I; "INTEGER" L, U, I; "ARRAY" A, B; "CODE" 34012; MIN:= EM[6]; "FOR" I:= 1 "STEP" 1 "UNTIL" N "DO" X1[I]:= "IF" VAL[I] <= MIN "THEN" 0 "ELSE" TAMVEC(1, N, I, V, X) / VAL[I]; "FOR" I:= 1 "STEP" 1 "UNTIL" M "DO" X[I]:= MATVEC(1, N, I, U, X1) "END" SOLSVDUND; "EOP" 0"CODE" 34283; "INTEGER" "PROCEDURE" SOLUND(A, M, N, X, EM); "VALUE" M, N; "INTEGER" M, N; "ARRAY" A, X, EM; "BEGIN" "INTEGER" I; "ARRAY" VAL[1:N], V[1:N,1:N]; "INTEGER" "PROCEDURE" QRISNGVALDEC(A, M, N, VAL, V, EM); "VALUE" M, N; "INTEGER" M, N; "ARRAY" A, VAL, V, EM; "CODE" 34273; "PROCEDURE" SOLSVDUND(U, VAL, V, M, N, X, EM); "VALUE" M, N; "INTEGER" M, N; "ARRAY" U, VAL, V, X, EM; "CODE" 34282; SOLUND:= I:= QRISNGVALDEC(A, M, N, VAL, V, EM); "IF" I = 0 "THEN" SOLSVDUND(A, VAL, V, M, N, X, EM) "END" SOLUND; "EOP" 1SECTION 3.1.1.3.1.3 (DECEMBER 1975) PAGE 1 AUTHOR : D.T.WINTER INSTITUTE : MATHEMATICAL CENTRE RECEIVED : 731217 BRIEF DESCRIPTION : THIS SECTION CONTAINS 2 PROCEDURES FOR THE CALCULATION OF THE HOMOGENEOUS EQUATIONS A * X = 0 AND X' * A = 0, WHERE A DENOTES A MATRIX, AND X A VECTOR: HOMSOLSVD ASSUMES THAT THE SINGULAR VALUES DECOMPOSITION OF A HAS BEEN GIVEN. HOMSOL FIRST CALCULATES THE SINGULAR VALUES DECOMPOSITION BY MEANS OF THE PROCEDURE QRISNGVALDEC. KEYWORDS : HOMOGENEOUS SOLUTION SINGULAR VALUES SUBSECTION : HOMSOLSVD CALLING SEQUENCE : THE HEADING OF THE PROCEDURE IS : "PROCEDURE" HOMSOLSVD(U, VAL, V, M, N); "VALUE" M, N; "INTEGER" M, N; "ARRAY" U, VAL, V; THE MEANING OF THE FORMAL PARAMETERS IS : U: ; "ARRAY" U[1:M,1:N]; ENTRY:THE MATRIX U IN THE SINGULAR VALUES DECOMPOSITION U*S*V'. EXIT: THE COLUMNS OF U THAT CORRESPOND TO THE ELEMENTS OF VAL WITH A VALUE SMALLER THAN SOME SMALL CONSTANT MAY BE SEEN AS THE SOLUTIONS OF X' * A = 0; VAL: ; "ARRAY" VAL[1:N]; ENTRY:THE SINGULAR VALUES; EXIT :THE ARRAY WILL BE ORDERED IN SUCH A WAY THAT VAL[I] < VAL[J] IF J < I; V: ; "ARRAY" V[1:N,1:N]; ENTRY:THE MATRIX V IN THE SINGULAR VALUES DECOMPOSITION; EXIT:THE COLUMNS OF V THAT CORRESPOND TO THE ELEMENTS OF VAL THAT ARE SMALLER THAN SOME SMALL CONSTANT MAY BE SEEN AS THE SOLUTIONS OF THE EQUATION A * X = 0; M: ; THE NUMBER OF ROWS OF U; N: ; THE NUMBER OF COLUMNS OF U; 1SECTION 3.1.1.3.1.3 (DECEMBER 1975) PAGE 2 PROCEDURES USED : ICHCOL = CP34031 RUNNING TIME : PROPORTIONAL TO N ' 2 METHOD AND PERFORMANCE : THE PROCEDURE DOES NOTHING MORE THAN A SIMPLE SORTING PROCESS ON THE ELEMENTS OF THE ARRAY VAL, AT THE SAME TIME THE COLUMNS OF U AND V ARE INTERCHANGED, ACCORDING TO THE INTERCHANGING OF THE ELEMENTS VAL. LANGUAGE : ALGOL 60 SUBSECTION : HOMSOL CALLING SEQUENCE : THE HEADING OF THE PROCEDURE IS : "INTEGER" "PROCEDURE" HOMSOL(A, M, N, V, EM); "VALUE" M, N; "INTEGER" M, N; "ARRAY" A, V, EM; HOMSOL:= THE NUMBER OF SINGULAR VALUES NOT FOUND, I.E. ZERO IF ALL SINGULAR VALUES ARE CALCULATED. THE MEANING OF THE FORMAL PAREMETERS IS : A: ; "ARRAY" A[1:M,1:N]; ENTRY: THE MATRIX; EXIT: THE COLUMNS OF A THAT CORRESPOND TO THE ELEMENTS OF VAL THAT ARE SMALLER THAN SOME SMALL CONSTANT, MAY BE SEEN AS THE SOLUTIONS OF THE EQUATION X' * A = 0. M: ; THE NUMBER OF ROWS OF A. N: ; THE NUMBER OF COLUMNS OF A. 1SECTION 3.1.1.3.1.3 (DECEMBER 1975) PAGE 3 V: ; "ARRAY" V[1:N,1:N]; EXIT: THE COLUMNS OF V THAT CORRESPOND TO ELEMENTS OF VAL SMALLER THAN SOME SMALL CONSTANT MAY BE SEEN AS THE SOLUTIONS OF THE EQUATION A * X = 0. EM: ; "ARRAY" EM[0:7]; ENTRY: EM[0]: THE MACHINE PRECISION; EM[2]: THE RELATIVE PRECISION FOR THE SINGULAR VALUES; EM[4]: THE MAXIMAL NUMBER OF ITERATIONS TO BE PERFORMED IN THE SINGULAR VALUES DECOMPOSITION. EM[6]: THE MINIMAL NON-NEGLECTABLE SINGULAR VALUE; EXIT: EM[1]: THE INFINITY NORM OF THE MATRIX; EM[3]: THE MAXIMAL NEGLECTED SUPERDIAGONAL ELEMENT; EM[5]: THE NUMBER OF ITERATIONS PERFORMED IN THE SINGULAR VALUES DECOMPOSITION; EM[7]: THE NUMERICAL RANK OF THE MATRIX, I.E. THE NUMBER OF SINGULAR VALUES GREATER THAN OR EQUAL TO EM[6]. PROCEDURES USED : QRISNGVALDEC = CP34273 HOMSOLSVD = CP34284 METHOD AND PERFORMANCE : THE SOLUTION IS FOUND IN TWO STEPS : 1. THE SINGULAR VALUES DECOMPOSITION IS CALCULATED BY MEANS OF THE PROCEDURE QRISNGVALDEC; 2. THE SINGULAR VALUES ARE ORDERED BY MEANS OF THE PROCEDURE HOMSOLSVD. RUNNING TIME : ROUGHLY PROPORTIONAL TO (M + N) * N * N LANGUAGE : ALGOL 60 REFERENCES : WILKINSON, J.H. AND C.REINSCH HANDBOOK OF AUTOMATIC COMPUTATION, VOL. 2 LINEAR ALGEBRA HEIDELBERG (1971) 1SECTION 3.1.1.3.1.3 (NOVEMBER 1976) PAGE 4 EXAMPLE OF USE : FIRST WE GIVE A PROGRAM, AND THAN THE RESULTS OF THIS PROGRAM : "BEGIN" "ARRAY" A[1:8,1:5], V[1:5,1:5], EM[0:7]; "INTEGER" I, J; A[1,1]:=22; A[1,2]:= A[2,3]:=10; A[1,3]:= A[7,1]:= A[8,5]:=2; A[1,4]:= A[3,5]:=3; A[1,5]:= A[2,2]:=7; A[2,1]:=14; A[2,5]:=8; A[2,4]:= A[8,3]:=0; A[3,1]:= A[3,3]:= A[6,5]:=-1; A[3,2]:=13; A[3,4]:=-11; A[4,1]:=-3; A[4,2]:= A[4,4]:= A[5,4]:= A[8,4]:=-2; A[4,3]:=13; A[4,5]:= A[5,5]:= A[8,1]:=4; A[5,1]:= A[6,1]:=9; A[5,2]:=8; A[5,3]:= A[6,2]:= A[7,5]:=1; A[6,3]:=-7; A[6,4]:= A[7,4]:= A[8,2]:=5; A[7,2]:=-6; A[7,3]:=6; EM[0]:="-14; EM[2]:="-12; EM[4]:=80; EM[6]:="-10; I:= HOMSOL(A, 8, 5, V, EM); OUTPUT(61, "("4B, "("NUMBER SINGULAR VALUES NOT FOUND : ")", 3ZD,/, 4B, "("NORM : ")", N,/, 4B, "("MAX NEGL SUBD ELEM : ")", N,/, 4B, "("NUMBER ITERATIONS : ")", 3ZD, /, 4B, "("RANK : ")", 3ZD, /")", I, EM[1], EM[3], EM[5], EM[7]); "FOR" J:= EM[7] + 1 "STEP" 1 "UNTIL" 5 "DO" OUTPUT(61, "("/, 4B, "("COLUMN NUMBER : ")", D, 5(/ ,4B, 2(N)), 3(/, 4B, N), /")", J, A[1,J], V[1,J], A[2,J], V[2,J], A[3,J], V[3,J], A[4,J], V[4,J], A[5,J], V[5,J], A[6,J], A[7,J], A[8,J]) "END" NUMBER SINGULAR VALUES NOT FOUND : 0 NORM : +4.4000000000000"+001 MAX NEGL SUBD ELEM : +4.3977072741076"-014 NUMBER ITERATIONS : 6 RANK : 3 COLUMN NUMBER : 4 +3.4708599800002"-001 -4.1909548511171"-001 -6.0723369623011"-001 +4.4050912303713"-001 +1.2207461910546"-001 -5.2004549247434"-002 +6.1882574433898"-001 +6.7605914021670"-001 -4.6344371870996"-003 +4.1297730284731"-001 +3.3409859838125"-001 -3.3528410857408"-002 -1.3547246422274"-002 COLUMN NUMBER : 5 -2.5533109413182"-001 +0.0000000000000"+000 -1.7359809248754"-001 -4.1854806384909"-001 -2.2081225414163"-001 -3.4879005320758"-001 +4.1165471593410"-002 -2.4415303724531"-001 +9.2044247057656"-001 +8.0221712237742"-001 -2.8895953996492"-002 +6.1327596621994"-002 -4.9058079025100"-002 1SECTION 3.1.1.3.1.3 (MAY 1974) PAGE 5 SOURCE TEXT(S): "CODE" 34284; "PROCEDURE" HOMSOLSVD(U, VAL, V, M, N); "VALUE" M, N; "INTEGER" M, N; "ARRAY" U, VAL, V; "BEGIN" "INTEGER" I, J; "REAL" X; "PROCEDURE" ICHCOL(L, U, I, J, A); "VALUE" L, U, I, J; "INTEGER" L, U, I, J; "ARRAY" A; "CODE" 34031; "FOR" I:= N "STEP" - 1 "UNTIL" 2 "DO" "FOR" J:= I - 1 "STEP" - 1 "UNTIL" 1 "DO" "IF" VAL[I] > VAL[J] "THEN" "BEGIN" X:= VAL[I]; VAL[I]:= VAL[J]; VAL[J]:= X; ICHCOL(1, M, I, J, U); ICHCOL(1, N, I, J, V) "END" "END" HOMSOLSVD; "EOP" 0"CODE" 34285; "INTEGER" "PROCEDURE" HOMSOL(A, M, N, V, EM); "VALUE" M, N; "INTEGER" M, N; "ARRAY" A, V, EM; "BEGIN" "INTEGER" I; "ARRAY" VAL[1:N]; "INTEGER" "PROCEDURE" QRISNGVALDEC(A, M, N, VAL, V, EM); "VALUE" M, N; "INTEGER" M, N; "ARRAY" A, VAL, V, EM; "CODE" 34273; "PROCEDURE" HOMSOLSVD(U, VAL, V, M, N); "VALUE" M, N; "INTEGER" M, N; "ARRAY" U, VAL, V; "CODE" 34284; HOMSOL:= I:= QRISNGVALDEC(A, M, N, VAL, V, EM); "IF" I = 0 "THEN" HOMSOLSVD(A, VAL, V, M, N) "END" HOMSOL; "EOP" 1SECTION 3.1.1.3.1.4 (DECEMBER 1979) PAGE 1 AUTHOR : D.T.WINTER INSTITUTE : MATHEMATICAL CENTRE RECEIVED : 731217 BRIEF DESCRIPTION : THIS SECTION CONTAINS TWO PROCEDURES FOR THE CALCULATION OF THE PSEUDO-INVERSE OF A MATRIX: PSDINVSVD ASSUMES THAT THE MATRIX IS GIVEN AS SINGULAR VALUES DECOMPOSITION. PSDINV FIRST CALCULATES THIS DECOMPOSITION. KEYWORDS : PSEUDO-INVERSE SINGULAR VALUES SUBSECTION : PSDINVSVD CALLING SEQUENCE : THE HEADING OF THE PROCEDURE IS : "PROCEDURE" PSDINVSVD((U, VAL, V, M, N, EM); "VALUE" M, N; "INTEGER" M, N; "ARRAY" U, VAL, V, EM; "CODE" 34286; THE MEANING OF THE FORMAL PARAMETERS IS : U: ; "ARRAY" U[1:M,1:N]; ENTRY: THE MATRIX U IN THE SINGULAR VALUES DECOMPOSITION U * S * V'; EXIT: THE TRANSPOSE OF THE PSEUDO-INVERSE. VAL: ; "ARRAY" VAL[1:N]; THE SINGULAR VALUES. V: ; "ARRAY" V[1:N,1:N]; THE MATRIX V IN THE SINGULAR VALUES DECOMPOSITION U * S * V'. M: ; THE NUMBER OF ROWS OF U. N: ; THE NUMBER OF COLUMNS OF V. EM: ; "ARRAY" EM[6:6]; ENTRY: EM[6]:THE MINIMAL NON-NEGLECTABLE SINGULAR VALUE. 1SECTION 3.1.1.3.1.4 (DECEMBER 1979) PAGE 2 PROCEDURES USED : MATVEC = CP34011 REQUIRED CENTRAL MEMORY : AN AUXILIARY ARRAY OF N REALS IS DECLARED RUNNING TIME : ROUGHLY PROPORTIONAL TO M * N * N METHOD AND PERFORMANCE : THE PSEUDO-INVERSE IS CALCULATED IN TWO STEPS : 1. THE MATRIX X = VAL+ * U' IS CALCULATED, WHERE VAL+ DENOTES THE DIAGONAL MATRIX OBTAINED FROM VAL BY PUTTING VAL+[I,I] = 1/VAL[I] IF VAL[I] GREATER THAN OR EQUAL TO EM[6], AND VAL+[I,I] = 0 OTHERWISE. 2. THE PSEUDO INVERSE (V * X) IS CALCULATED. SUBSECTION : PSDINV CALLING SEQUENCE : THE HEADING OF THE PROCEDURE IS : "INTEGER" "PROCEDURE" PSDINV(A, M, N, EM); "VALUE" M, N; "INTEGER" M, N; "ARRAY" A, EM; "CODE" 34287; PSDINV:= THE NUMBER OF SINGULAR VALUES NOT FOUND, I.E. ZERO IF ALL SINGULAR VALUES ARE CALCULATED. THE MEANING OF THE FORMAL PARAMETERS IS : A ; "ARRAY" A[1:M,1:N]; ENTRY : THE GIVEN MATRIX; EXIT : THE TRANSPOSE OF THE PSEUDO-INVERSE; M: ; THE NUMBER OF ROWS OF A; N: ; THE NUMBER OF COLUMNS OF A, N<= M; EM: ; "ARRAY" EM[0:7]; ENTRY: EM[0]: THE MACHINE PRECISION; EM[2]: THE RELATIVE PRECISION FOR THE SINGULAR VALUES; EM[4]: THE MAXIMAL NUMBER OF ITERATIONS TO BE PERFORMED; EM[6]: THE MINIMAL NON-NEGLECTABLE SINGULAR VALUE; EXIT: EM[1]: THE INFINITY NORM OF THE MATRIX; EM[3]: THE MAXIMAL NEGLECTED SUPERDIAGONAL ELEMENT; EM[5]: THE NUMBER OF ITERATIONS PERFORMED IN THE SINGULAR VALUES DECOMPOSITION; EM[7]: THE NUMERICAL RANK OF THE MATRIX, I.E. THE NUMBER OF SINGULAR VALUES GREATER THAN OR EQUAL TO EM[6]; 1SECTION 3.1.1.3.1.4 (DECEMBER 1979) PAGE 3 PROCEDURES USED : QRISNGVALDEC = CP34273 PSDINVSVD = CP34286 REQUIRED CENTRAL MEMORY : AUXILIARY ARRAY'S ARE DECLARED TO A TOTAL OF (N + 1) * N REALS RUNNING TIME : ROUGHLY PROPORTIONAL TO (2M + N) * N * N METHOD AND PERFORMANCE : FIRST THE SINGULAR VALUES DECOMPOSITION IS CALCULATED, AND THAN THE PSEUDO-INVERSE IS CALCULATED BY PSDINVSVD. REFERENCES : WILKINSON, J.H. AND C.REINSCH HANDBOOK OF AUTOMATIC COMPUTATION, VOL.2 LINEAR ALGEBRA HEIDELBERG (1971) 1SECTION 3.1.1.3.1.4 (DECEMBER 1979) PAGE 4 EXAMPLE OF USE : FIRST WE GIVE A PROGRAM, AND THEN THE RESULTS OF THIS PROGRAM : "BEGIN" "ARRAY" A[1:8,1:5], EM[0:7]; "INTEGER" I, J; "INTEGER" "PROCEDURE" PSDINV(A, M, N, EM); "VALUE" M, N; "INTEGER" M, N; "ARRAY" A, EM; "CODE" 34287; A[1,1]:=22; A[1,2]:= A[2,3]:=10; A[1,3]:= A[7,1]:= A[8,5]:=2; A[1,4]:= A[3,5]:=3; A[1,5]:= A[2,2]:=7; A[2,1]:=14; A[2,5]:=8; A[2,4]:= A[8,3]:=0; A[3,1]:= A[3,3]:= A[6,5]:=-1; A[3,2]:=13; A[3,4]:=-11; A[4,1]:=-3; A[4,2]:= A[4,4]:= A[5,4]:= A[8,4]:=-2; A[4,3]:=13; A[4,5]:= A[5,5]:= A[8,1]:=4; A[5,1]:= A[6,1]:=9; A[5,2]:=8; A[5,3]:= A[6,2]:= A[7,5]:=1; A[6,3]:=-7; A[6,4]:= A[7,4]:= A[8,2]:=5; A[7,2]:=-6; A[7,3]:=6; EM[0]:="-14; EM[2]:="-12; EM[4]:=80; EM[6]:="-10; I:= PSDINV(A, 8, 5, EM); OUTPUT(61, "("4B, "("NUMBER SINGULAR VALUES NOT FOUND : ")", 3ZD,/, 4B, "("NORM : ")", N,/, 4B, "("MAX NEGL SUBD ELEM : ")", N,/, 4B, "("NUMBER ITERATIONS : ")", 3ZD, /, 4B, "("RANK : ")", 3ZD, /")", I, EM[1], EM[3], EM[5], EM[7]); OUTPUT(61, "("/, 4B, "("TRANSPOSE OF PSEUDO-INVERSE")", /, 4B, "("FIRST THREE COLUMNS")", /, /, 8(4B, 3(N), /), /, /, 4B, "("LAST TWO COLUMNS")", /, /, 8(15B, 2(N), /)")", A[1,1], A[1,2], A[1,3], A[2,1], A[2,2], A[2,3], A[3,1], A[3,2], A[3,3], A[4,1], A[4,2], A[4,3], A[5,1], A[5,2], A[5,3], A[6,1], A[6,2], A[6,3], A[7,1], A[7,2], A[7,3], A[8,1], A[8,2], A[8,3], A[1,4], A[1,5], A[2,4], A[2,5], A[3,4], A[3,5], A[4,4], A[4,5], A[5,4], A[5,5], A[6,4], A[6,5], A[7,4], A[7,5], A[8,4], A[8,5]) "END" NUMBER SINGULAR VALUES NOT FOUND : 0 NORM : +4.4000000000000"+001 MAX NEGL SUBD ELEM : +4.3977072741076"-014 NUMBER ITERATIONS : 6 RANK : 3 TRANSPOSE OF PSEUDO-INVERSE FIRST THREE COLUMNS +2.1129807692308"-002 +4.6153846153850"-003 -2.1073717948727"-003 +9.3108974358974"-003 +2.2115384615376"-003 +2.0528846153848"-002 -1.1097756410256"-002 +2.7403846153848"-002 -3.8862179487199"-003 -7.9166666666667"-003 -5.0000000000007"-003 +3.3750000000001"-002 +5.5128205128205"-003 +9.8076923076935"-003 -8.9743589743826"-004 +1.4318910256410"-002 -2.5961538461548"-003 -2.0136217948716"-002 +4.8958333333335"-003 -1.4999999999998"-002 +1.5312499999996"-002 +1.5064102564102"-003 +7.4038461538447"-003 -1.6987179487147"-003 1SECTION 3.1.1.3.1.4 (MAY 1974) PAGE 5 LAST TWO COLUMNS +7.6041666666662"-003 +3.8060897435894"-003 -2.0833333333295"-004 +1.0016025641026"-002 -2.7604166666667"-002 +4.2067307692303"-003 -5.4166666666662"-003 +1.0416666666667"-002 -5.0000000000005"-003 +3.2051282051275"-003 +1.2812500000000"-002 -6.2099358974354"-003 +1.2395833333332"-002 +2.6041666666656"-003 -4.9999999999993"-003 +1.6025641025649"-003 SOURCE TEXT(S): 0"CODE" 34286; "PROCEDURE" PSDINVSVD(U, VAL, V, M, N, EM); "VALUE" M, N; "INTEGER" M, N; "ARRAY" U, VAL, V, EM; "BEGIN" "INTEGER" I, J; "REAL" MIN, VALI; "ARRAY" X[1:N]; "REAL" "PROCEDURE" MATVEC(L, U, I, A, B); "VALUE" L, U, I; "INTEGER" L, U, I; "ARRAY" A, B; "CODE" 34011; MIN:= EM[6]; "FOR" I:= 1 "STEP" 1 "UNTIL" N "DO" "IF" VAL[I] > MIN "THEN" "BEGIN" VALI:= 1 / VAL[I]; "FOR" J:= 1 "STEP" 1 "UNTIL" M "DO" U[J,I]:= U[J,I] * VALI "END" "ELSE" "FOR" J:= 1 "STEP" 1 "UNTIL" M "DO" U[J,I]:= 0; "FOR" I:= 1 "STEP" 1 "UNTIL" M "DO" "BEGIN" "FOR" J:= 1 "STEP" 1 "UNTIL" N "DO" X[J]:= U[I,J]; "FOR" J:= 1 "STEP" 1 "UNTIL" N "DO" U[I,J]:= MATVEC(1, N, J, V, X) "END" "END" PSDINVSVD; "EOP" 0"CODE" 34287; "INTEGER" "PROCEDURE" PSDINV(A, M, N, EM); "VALUE" M, N; "INTEGER" M, N; "ARRAY" A, EM; "BEGIN" "INTEGER" I; "ARRAY" VAL[1:N], V[1:N,1:N]; "INTEGER" "PROCEDURE" QRISNGVALDEC(A, M, N, VAL, V, EM); "VALUE" M, N; "INTEGER" M, N; "ARRAY" A, VAL, V, EM; "CODE" 34273; "PROCEDURE" PSDINVSVD(U, VAL, V, M, N, EM); "VALUE" M, N; "INTEGER" M, N; "ARRAY" U, VAL, V, EM; "CODE" 34286; PSDINV:= I:= QRISNGVALDEC(A, M, N, VAL, V, EM); "IF" I = 0 "THEN" PSDINVSVD(A, VAL, V, M, N, EM) "END" PSDINV; "EOP" 1SECTION 3.1.2.1.1.1.1.1 (DECEMBER 1975) PAGE 1 AUTHOR : T.J. DEKKER. REVISOR : J. KOK. INSTITUTE : MATHEMATICAL CENTRE. RECEIVED : 730903. BRIEF DESCRIPTION : THIS SECTION CONTAINS THE PROCEDURE DECBND FOR THE DECOMPOSITION OF A BAND MATRIX BY GAUSSIAN ELIMINATION WITH STABILIZING ROW INTERCHANGES (PARTIAL PIVOTING). KEY WORDS : LINEAR EQUATIONS, PARTIAL PIVOTING, GAUSSIAN ELIMINATION, BAND MATRIX. CALLING SEQUENCE : THE HEADING OF THE PROCEDURE IS : "PROCEDURE" DECBND(A, N, LW, RW, AUX, M, P); "VALUE" N, LW, RW; "INTEGER" N, LW, RW; "INTEGER""ARRAY" P; "ARRAY" A, M, AUX; THE MEANING OF THE FORMAL PARAMETERS IS : A : ; "ARRAY" A[1 : (LW + RW) * (N - 1) + N]; ENTRY : A CONTAINS ROWWISE THE BAND ELEMENTS OF THE BAND MATRIX IN SUCH A WAY THAT THE (I,J)-TH ELEMENT OF THE MATRIX IS GIVEN IN A[( LW + RW ) * (I - 1) + J], I=1,...,N AND J=MAX(1,I-LW),...,MIN(N,I+RW). THE VALUES OF THE REMAINING ELEMENTS OF A ARE IRRELEVANT. EXIT : THE BAND ELEMENTS OF THE GAUSSIAN ELIMINATED MATRIX, WHICH IS AN UPPERTRIANGULAR BAND MATRIX WITH (LW + RW) CODIAGONALS, ARE ROWWISE DELIVERED IN A AS FOLLOWS: THE (I,J)-TH ELEMENT OF U IS A [ ( LW + RW ) * (I - J) + J ],I=1,...,N AND J=I,...,MIN(N,I + LW + RW). N : ; ORDER OF THE BAND MATRIX; LW : ; NUMBER OF LEFT CODIAGONALS OF A; 1SECTION 3.1.2.1.1.1.1.1 (JUNE 1974) PAGE 2 RW : ; NUMBER OF RIGHT CODIAGONALS OF A; AUX : ; "ARRAY" AUX[1 : 5]; ENTRY :AUX[2] = EPS IS A RELATIVE TOLERANCE TO CONTROL THE ELIMINATION; THE PROCESS IS DISCONTINUED IF (EPS > PIVOT[I] / EUCLIDEAN NORM OF I-TH ROW) IN THE I-TH ELIMINATION STEP; NORMAL EXIT : AUX[1] = SIGN OF THE DETERMINANT OF THE MATRIX (+1 OR -1); AUX[3] = N; AUX[5] = MINIMUM ABSOLUTE VALUE OF PIVOT[I] / EUCLIDEAN NORM OF THE I-TH ROW; ABNORMAL EXIT : IF THE ELIMINATION CANNOT BE CARRIED OUT, I.E. IF TEMP (THE QUANTITY ABS(PIVOT[I] / EUCLIDEAN NORM OF THE I-TH ROW)) IS TOO SMALL IN THE I-TH ELIMINATION STEP : AUX[3] = I - 1, AUX[5] = TEMP; M : ; "ARRAY" M[1 : LW * (N - 2) + 1]; EXIT : THE GAUSSIAN MULTIPLIERS OF ALL ELIMINATIONS IN SUCH A WAY THAT THE I-TH MULTIPLIER OF THE J-TH STEP IS M [ LW * (J - 1) + I - J ]. P : ; "INTEGER""ARRAY" P[1 : N]; EXIT : THE PIVOTAL INDICES. PROCEDURES USED : VECVEC = CP34010, ELMVEC = CP34020, ICHVEC = CP34030. REQUIRED CENTRAL MEMORY : EXECUTION FIELD LENGTH : A REAL ARRAY OF N ELEMENTS IS DECLARED. RUNNING TIME : (C1 * LW + C2) * (LW + RW + 1) * N; THE CONSTANTS C1 AND C2 DEPEND UPON THE ARITHMETIC OF THE COMPUTER. 1SECTION 3.1.2.1.1.1.1.1 (JUNE 1974) PAGE 3 LANGUAGE : ALGOL 60. METHOD AND PERFORMANCE : DECBND PERFORMS THE DECOMPOSITION OF A MATRIX WHOSE NON-ZERO ELEMENTS ARE IN BAND FORM, AND WHOSE BAND ELEMENTS ARE STORED ROWWISE IN A ONE-DIMENSIONAL ARRAY. THE METHOD USED IS GAUSSIAN ELIMINATION WITH STABILIZING ROW INTERCHANGES (PARTIAL PIVOTING). THE GAUSSIAN ELIMINATION IS PERFORMED IN N STEPS. IN THE K-TH STEP, K = 1, ... , N, A PIVOT IS SELECTED IN THE K-TH COLUMN OF THE REMAINING SUBMATRIX OF ORDER N - K + 1 (THIS COLUMN CONTAINS AT MOST LW + 1 NON-ZERO ELEMENTS); THEN THE PIVOTAL ROW IS INTERCHANGED WITH THE K-TH ROW; SUBSEQUENTLY THE K-TH UNKNOWN IS ELIMINATED IN THE LAST N - K ROWS (ONLY THE FIRST LW OF THESE LAST ROWS ARE INVOLVED HERE). THE PIVOT IS SELECTED IN SUCH A WAY THAT ITS ABSOLUTE VALUE DIVIDED BY THE EUCLIDEAN NORM OF THE CORRESPONDING ROW OF THE MATRIX IS MAXIMAL. THUS, THE MATRIX IS EQUILIBRATED IN THIS PIVOTING STRATEGY SUCH, THAT THE ROWS EFFECTIVELY OBTAIN UNIT EUCLIDEAN NORM. THE PROCEDURE DELIVERS THE BAND ELEMENTS OF THE ELIMINATED MATRIX (WHICH IS AN UPPER TRIANGULAR MATRIX WITH LW + RW SUPERDIAGONALS) AND THE GAUSSIAN MULTIPLIERS FOR EACH ELIMINATION . THE ELIMINATION CANNOT BE CARRIED OUT IF THE ABSOLUTE VALUE OF THE PIVOT IS LESS THAN A GIVEN RELATIVE TOLERANCE (AUX[2]) TIMES THE EUCLIDEAN NORM OF THE CORRESPONDING ROW OF THE MATRIX. THEN THE PREVIOUS STEP NUMBER OF THE ELIMINATION IS DELIVERED (IN AUX[3], WHICH ELSE TAKES THE VALUE N). SEE ALSO REF [1], SECTION 212. REFERENCE : [1] DEKKER, T.J. : ALGOL 60 PROCEDURES IN NUMERICAL ALGEBRA, PART 1, MC TRACT 22, 1968, MATHEMATISCH CENTRUM, AMSTERDAM. EXAMPLE OF USE : SEE EXAMPLE OF USE OF SOLBND. 1SECTION 3.1.2.1.1.1.1.1 (JUNE 1974) PAGE 4 SOURCE TEXT(S) : 0"CODE" 34320; "PROCEDURE" DECBND(A, N, LW, RW, AUX, M, P); "VALUE" N, LW, RW; "INTEGER" N, LW, RW; "INTEGER" "ARRAY" P; "ARRAY" A, M, AUX; "BEGIN" "INTEGER" I, J, K, KK, KK1, PK, MK, IK, LW1, F, Q, W, W1, W2, NRW, IW, SDET; "REAL" R, S, EPS, MIN; "ARRAY" V[1:N]; "REAL""PROCEDURE" VECVEC(A, B, C, D, E); "CODE" 34010; "PROCEDURE" ELMVEC(A, B, C, D, E, F); "CODE" 34020; "PROCEDURE" ICHVEC(A, B, C, D); "CODE" 34030; F:= LW; W1:= LW + RW; W:= W1 + 1; W2:= W - 2; IW:= 0; SDET:= 1; NRW:= N - RW; LW1:= LW + 1; Q:= LW - 1; "FOR" I:= 2 "STEP" 1 "UNTIL" LW "DO" "BEGIN" Q:= Q - 1; IW:= IW + W1; "FOR" J:= IW - Q "STEP" 1 "UNTIL" IW "DO" A[J]:= 0 "END"; IW:= - W2; Q:= - LW; "FOR" I:= 1 "STEP" 1 "UNTIL" N "DO" "BEGIN" IW:= IW + W; "IF" I <= LW1 "THEN" IW:= IW - 1; Q:= Q + W; "IF" I > NRW "THEN" Q:= Q - 1; V[I]:= SQRT(VECVEC(IW, Q, 0, A, A)) "END"; EPS:= AUX[2]; MIN:= 1; KK:= - W1; MK:= - LW; "IF" F > NRW "THEN" W2:= W2 + NRW - F; "FOR" K:= 1 "STEP" 1 "UNTIL" N "DO" "BEGIN""IF" F < N "THEN" F:= F + 1; IK:= KK:= KK + W; MK:= MK + LW; S:= ABS(A[KK]) / V[K]; PK:= K; KK1:= KK + 1; "FOR" I:= K + 1 "STEP" 1 "UNTIL" F "DO" "BEGIN" IK:= IK + W1; M[MK + I - K]:= R:= A[IK]; A[IK]:= 0; R:= ABS(R) / V[I]; "IF" R > S "THEN" "BEGIN" S:= R; PK:= I "END" "END"; "IF" S < MIN "THEN" MIN:= S; "IF" S < EPS "THEN" "BEGIN" AUX[3]:= K - 1; AUX[5]:= S; "GO TO" END "END"; "IF" K + W2 >= N "THEN" W2:= W2 - 1; P[K]:= PK; "IF" PK ^= K "THEN" "BEGIN" V[PK]:= V[K]; PK:= PK - K; ICHVEC(KK1, KK1 + W2, PK * W1, A); SDET:= - SDET; R:= M[MK + PK]; M[MK + PK]:= A[KK]; A[KK]:= R "END""ELSE" R:= A[KK]; "IF" R < 0 "THEN" SDET:= - SDET; IW:= KK1; LW1:= F - K + MK; "FOR" I:= MK + 1 "STEP" 1 "UNTIL" LW1 "DO" "BEGIN" M[I]:= S:= M[I] / R; IW:= IW + W1; ELMVEC(IW, IW + W2, KK1 - IW, A, A, - S) "END" "END"; AUX[3]:= N; AUX[5]:= MIN; END: AUX[1]:= SDET "END" DECBND; "EOP" 1SECTION 3.1.2.1.1.1.1.2 (JUNE 1974) PAGE 1 CONTRIBUTOR : J. KOK. INSTITUTE : MATHEMATICAL CENTRE. RECEIVED : 730903. BRIEF DESCRIPTION : THIS SECTION CONTAINS THE PROCEDURE DETERMBND FOR THE CALCULATION OF THE DETERMINANT OF A BAND MATRIX. KEY WORDS : DETERMINANT, BAND MATRIX. CALLING SEQUENCE : THE HEADING OF THE PROCEDURE IS : "REAL""PROCEDURE" DETERMBND(A, N, LW, RW, SGNDET); "VALUE" N, LW, RW, SGNDET; "INTEGER" N, LW, RW, SGNDET; "ARRAY" A; DETERMBND DELIVERS THE DETERMINANT OF THE MATRIX. THE MEANING OF THE FORMAL PARAMETERS IS : A, N, LW, RW : SEE CALLING SEQUENCE OF DECBND (SECTION 3.1.2.1.1.1.1.1.); ENTRY : THE CONTENTS OF A ARE AS PRODUCED BY DECBND OR DECSOLBND (SECTION 3.1.2.1.1.1.1.3.); SGNDET : ; ENTRY : THE SIGN OF THE DETERMINANT AS DELIVERED IN AUX[1] BY DECBND, IF THE ELIMINATION BY DECBND WAS SUCCESSFUL. 1SECTION 3.1.2.1.1.1.1.2 (DECEMBER 1975) PAGE 2 PROCEDURES USED : NONE. RUNNING TIME : PROPORTIONAL TO N. LANGUAGE : ALGOL 60. METHOD AND PERFORMANCE : DETERMBND CAN BE CALLED AFTER DECBND OR DECSOLBND ONLY IF THE GAUSSIAN ELIMINATION WAS SUCCESSFUL, I.E. IF AUX[3] = N. THE FUNCTION VALUE OF DETERMBND IS THE DETERMINANT OF THE GAUSSIAN ELIMINATED UPPER TRIANGULAR MATRIX PROVIDED WITH THE CORRECT SIGN THAT IS DELIVERED BY DECBND OR DECSOLBND IN AUX[1]. DETERMBND SHOULD NOT BE CALLED WHEN OVERFLOW CAN BE EXPECTED. EXAMPLE OF USE : SEE EXAMPLES OF USE OF SOLBND AND DECSOLBND. SOURCE TEXT(S) : 0"CODE" 34321; "REAL""PROCEDURE" DETERMBND(A, N, LW, RW, SGNDET); "VALUE" N, LW, RW, SGNDET; "INTEGER" N, LW, RW, SGNDET; "ARRAY" A; "BEGIN""INTEGER" I, L; "REAL" P; L:= 1; P:= 1; LW:= LW + RW + 1; "FOR" I:= 1 "STEP" 1 "UNTIL" N "DO" "BEGIN" P:= A[L] * P; L:= L + LW "END"; DETERMBND:= ABS(P) * SGNDET "END" DETERMBND; "EOP" 1SECTION 3.1.2.1.1.1.1.3 (JUNE 1974) PAGE 1 AUTHOR : T.J. DEKKER. REVISOR : J. KOK. INSTITUTE : MATHEMATICAL CENTRE. RECEIVED : 730903. BRIEF DESCRIPTION : THIS SECTION CONTAINS TWO PROCEDURES. A) SOLBND, FOR THE SOLUTION OF ONE OR MORE SYSTEMS OF LINEAR EQUATIONS WITH THE SAME COEFFICIENT MATRIX, IF THIS MATRIX HAS BEEN DECOMPOSED BY A CALL OF THE PROCEDURE DECBND (SECTION 3.1.2.1.1.1.1.1.). B) DECSOLBND, FOR THE SOLUTION OF ONE SYSTEM OF LINEAR EQUATIONS BY GAUSSIAN ELIMINATION WITH STABILIZING ROW INTERCHANGES (PARTIAL PIVOTING) IF THE COEFFICIENT MATRIX IS IN BAND FORM AND IS STORED ROWWISE IN A ONE-DIMENSIONAL ARRAY. KEY WORDS : LINEAR EQUATIONS, PARTIAL PIVOTING, GAUSSIAN ELIMINATION, BAND MATRIX. SUBSECTION : SOLBND. CALLING SEQUENCE : THE HEADING OF THE PROCEDURE IS : "PROCEDURE" SOLBND(A, N, LW, RW, M, P, B); "VALUE" N, LW, RW; "INTEGER" N, LW, RW; "INTEGER""ARRAY" P; "ARRAY" A, M, B; THE MEANING OF THE FORMAL PARAMETERS IS : A, N, LW, RW, M, P : SEE CALLING SEQUENCE OF DECBND, ENTRY : THE CONTENTS OF THE ARRAYS A, M, P ARE AS PRODUCED BY DECBND; B : ; "ARRAY" B[1 : N]; ENTRY : THE RIGHT HAND SIDE OF THE SYSTEM OF LINEAR EQUATIONS; EXIT : THE SOLUTION OF THE SYSTEM. 1SECTION 3.1.2.1.1.1.1.3 (JUNE 1974) PAGE 2 PROCEDURES USED : VECVEC = CP34010, ELMVEC = CP34020. RUNNING TIME : (C3 * LW + C4 * RW + C5) * N; THE CONSTANTS C3, C4 AND C5 DEPEND UPON THE ARITHMETIC OF THE COMPUTER. LANGUAGE : ALGOL 60. METHOD AND PERFORMANCE : SOLBND CALCULATES THE SOLUTION OF A SYSTEM OF LINEAR EQUATIONS, PROVIDED THAT THE MATRIX WAS DECOMPOSED BY A SUCCESSFUL CALL OF DECBND (SECTION 3.1.2.1.1.1.1.1.). THE SOLUTION OF THE LINEAR SYSTEM IS OBTAINED BY CARRYING OUT THE ELIMINATIONS, FOR WHICH THE GAUSSIAN MULTIPLIERS ARE SAVED, ON THE RIGHT HAND SIDE, AND BY SOLVING THE NEW SYSTEM WITH THE UPPER TRIANGULAR BAND MATRIX , AS PRODUCED BY DECBND, BY BACK SUBSTITUTION. THE SOLUTIONS OF SEVERAL SYSTEMS WITH THE SAME COEFFICIENT MATRIX CAN BE OBTAINED BY SUCCESSIVE CALLS OF SOLBND. EXAMPLE OF USE : THE FOLLOWING PROGRAM SOLVES THE SYSTEM OF SIMULTANEOUS EQUATIONS 2 * X1 - X2 = 1 - X1 + 2 * X2 - X3 = 0 - X2 + 2 * X3 - X4 = 0 - X3 + 2 * X4 - X5 = 0 - X4 + 2 * X5 = 1 "BEGIN""COMMENT" 730822, TEST DECBND, SOLBND AND DETERMBND; "PROCEDURE" DECBND(A, N, LW, RW, AUX, M, P); "CODE" 34320; "PROCEDURE" SOLBND(A, N, LW, RW, M, P, B); "CODE" 34071; "REAL""PROCEDURE" DETERMBND(A, N, LW, RW, SGNDET); "CODE" 34321; "INTEGER" I; "INTEGER""ARRAY" ROWIND[1 : 5]; "ARRAY" BAND[1 : 13], MULT[1 : 4], RIGHT, AUX[1 : 5]; "FOR" I:= 1 "STEP" 1 "UNTIL" 13 "DO" BAND[I]:= "IF" (I + 1) // 3 * 3 < I "THEN" 2 "ELSE" - 1; RIGHT[1]:= RIGHT[5]:= 1; "FOR" I:= 2, 3, 4 "DO" RIGHT[I]:= 0; AUX[2]:= "- 12; DECBND(BAND, 5, 1, 1, AUX, MULT, ROWIND); "IF" AUX[3] = 5 "THEN" "BEGIN" SOLBND(BAND, 5, 1, 1, MULT, ROWIND, RIGHT); OUTPUT(61, "("5(+2Z.4D2B), /"("DETERMINANT IS ")" +.8D"+DD ")", (RIGHT[I], I:= 1 : 5), DETERMBND(BAND, 5, 1, 1, AUX[1])) "END" "END" DELIVERS : +1.0000 +1.0000 +1.0000 +1.0000 +1.0000 DETERMINANT IS +.60000000"+01 1SECTION 3.1.2.1.1.1.1.3 (DECEMBER 1975) PAGE 3 SUBSECTION : DECSOLBND. CALLING SEQUENCE : THE HEADING OF THE PROCEDURE IS : "PROCEDURE" DECSOLBND(A, N, LW, RW, AUX, B); "VALUE" N, LW, RW; "INTEGER" N, LW, RW; "ARRAY" A, AUX, B; THE MEANING OF THE FORMAL PARAMETERS IS : A, N, LW, RW, AUX : SEE DECBND (SECTION : 3.1.2.1.1.1.1.1); B : SEE SOLBND (THIS SECTION). PROCEDURES USED : VECVEC = CP34010, ELMVEC = CP34020, ICHVEC = CP34030. REQUIRED CENTRAL MEMORY : EXECUTION FIELD LENGTH : A REAL ARRAY OF N ELEMENTS AND A REAL ARRAY OF LW + 1 ELEMENTS ARE DECLARED. RUNNING TIME : (C1 * LW + C6) * (LW + RW + 1) * N; THE CONSTANTS C1 AND C6 DEPEND UPON THE ARITHMETIC OF THE COMPUTER. LANGUAGE : ALGOL 60. 1SECTION 3.1.2.1.1.1.1.3 (JUNE 1974) PAGE 4 METHOD AND PERFORMANCE : DECSOLBND PERFORMS GAUSSIAN ELIMINATION IN THE SAME WAY AS DECBND , MEANWHILE ALSO CARRYING OUT THE ELIMINATION WITH THE GIVEN RIGHT HAND SIDE. THE SOLUTION OF THE ELIMINATED SYSTEM IS OBTAINED BY BACK SUBSTITUTION. EXAMPLE OF USE : THE PROGRAM "BEGIN""COMMENT" 730822, TEST DECSOLBND AND DETERMBND; "PROCEDURE" DECSOLBND(A, N, LW, RW, AUX, B); "CODE" 34322; "REAL""PROCEDURE" DETERMBND(A, N, LW, RW, SGNDET); "CODE" 34321; "INTEGER" I; "ARRAY" BAND[1 : 13], RIGHT, AUX[1 : 5]; "FOR" I:= 1 "STEP" 1 "UNTIL" 13 "DO" BAND[I]:= "IF" (I + 1) // 3 * 3 < I "THEN" 2 "ELSE" - 1; RIGHT[1]:= RIGHT[5]:= 1; "FOR" I:= 2, 3, 4 "DO" RIGHT[I]:= 0; AUX[2]:= "- 12; DECSOLBND(BAND, 5, 1, 1, AUX, RIGHT); "IF" AUX[3] = 5 "THEN" "BEGIN" OUTPUT(61, "("5(+2Z.4D2B), /"("DETERMINANT IS ")" +.8D"+DD ")", (RIGHT[I], I:= 1 : 5), DETERMBND(BAND, 5, 1, 1, AUX[1])) "END" "END" WHICH SOLVES THE SAME PROBLEM AS THE PROGRAM IN THE EXAMPLE OF USE OF SOLBND, DELIVERS : +1.0000 +1.0000 +1.0000 +1.0000 +1.0000 DETERMINANT IS +.60000000"+01 1SECTION 3.1.2.1.1.1.1.3 (JUNE 1974) PAGE 5 SOURCE TEXT(S) : 0"CODE" 34071; "PROCEDURE" SOLBND(A, N, LW, RW, M, P, B); "VALUE" N, LW, RW; "INTEGER" N, LW, RW; "INTEGER" "ARRAY" P; "ARRAY" A, B, M; "BEGIN" "INTEGER" F, I, K, KK, W, W1, W2, SHIFT; "REAL" S; "REAL""PROCEDURE" VECVEC(A, B, C, D, E); "CODE" 34010; "PROCEDURE" ELMVEC(A, B, C, D, E, F); "CODE" 34020; F:= LW; SHIFT:= - LW; W1:= LW - 1; "FOR" K:= 1 "STEP" 1 "UNTIL" N "DO" "BEGIN""IF" F < N "THEN" F:= F + 1; SHIFT:= SHIFT + W1; I:=P[K]; S:= B[I]; "IF" I ^= K "THEN" "BEGIN" B[I]:= B[K]; B[K]:= S "END"; ELMVEC(K + 1, F, SHIFT, B, M, - S) "END"; W1:= LW + RW; W:= W1 + 1; KK:= (N + 1) * W - W1; W2:= - 1; SHIFT:= N * W1; "FOR" K:= N "STEP" - 1 "UNTIL" 1 "DO" "BEGIN" KK:= KK - W; SHIFT:= SHIFT - W1; "IF" W2 < W1 "THEN" W2:= W2 + 1; B[K]:= (B[K] - VECVEC(K + 1, K + W2, SHIFT, B, A)) / A[KK] "END" "END" SOLBND; "EOP" 0"CODE" 34322; "PROCEDURE" DECSOLBND(A, N, LW, RW, AUX, B); "VALUE" N, LW, RW; "INTEGER" N, LW, RW; "ARRAY" A, B, AUX; "BEGIN""INTEGER" I, J, K, KK, KK1, PK, IK, LW1, F, Q, W, W1, W2,IW, NRW, SHIFT, SDET; "REAL" R, S, EPS, MIN; "ARRAY" M[0:LW], V[1:N]; "REAL""PROCEDURE" VECVEC(A, B, C, D, E); "CODE" 34010; "PROCEDURE" ELMVEC(A, B, C, D, E, F); "CODE" 34020; "PROCEDURE" ICHVEC(A, B, C, D); "CODE" 34030; F:= LW; SDET:= 1; W1:= LW + RW; W:= W1 + 1; W2:= W - 2; IW:= 0; NRW:= N - RW; LW1:= LW + 1; Q:= LW - 1; "FOR" I:= 2 "STEP" 1 "UNTIL" LW "DO" "BEGIN" Q:= Q - 1; IW:= IW + W1; "FOR" J:= IW - Q "STEP" 1 "UNTIL" IW "DO" A[J]:= 0 "END"; "COMMENT" 1SECTION 3.1.2.1.1.1.1.3 (JUNE 1974) PAGE 6 ; IW:= - W2; Q:= - LW; "FOR" I:= 1 "STEP" 1 "UNTIL" N "DO" "BEGIN" IW:= IW + W; "IF" I <= LW1 "THEN" IW:= IW - 1; Q:= Q + W; "IF" I > NRW "THEN" Q:= Q - 1; V[I]:= SQRT(VECVEC(IW, Q, 0, A, A)) "END"; EPS:= AUX[2]; MIN:= 1; KK:= - W1; "IF" F > NRW "THEN" W2:= W2 + NRW - F; "FOR" K:= 1 "STEP" 1 "UNTIL" N "DO" "BEGIN""IF" F < N "THEN" F:= F + 1; IK:= KK:= KK + W; S:= ABS(A[KK]) / V[K]; PK:= K; KK1:= KK + 1; "FOR" I:= K + 1 "STEP" 1 "UNTIL" F "DO" "BEGIN" IK:= IK + W1; M[I - K]:= R:= A[IK]; A[IK]:= 0; R:= ABS(R) / V[I]; "IF" R > S "THEN" "BEGIN" S:= R; PK:= I "END" "END"; "IF" S < MIN "THEN" MIN:= S; "IF" S < EPS "THEN" "BEGIN" AUX[3]:= K - 1; AUX[5]:= S; "GO TO" END "END"; "IF" K + W2 >= N "THEN" W2:= W2 - 1; "IF" PK ^= K "THEN" "BEGIN" V[PK]:= V[K]; PK:= PK - K; ICHVEC(KK1, KK1 + W2, PK * W1, A); SDET:= - SDET; R:= B[K]; B[K]:= B[PK + K]; B[PK + K]:= R; R:= M[PK]; M[PK]:= A[KK]; A[KK]:= R "END" "ELSE" R:= A[KK]; IW:= KK1; LW1:= F - K; "IF" R < 0 "THEN" SDET:= - SDET; "FOR" I:= 1 "STEP" 1 "UNTIL" LW1 "DO" "BEGIN" M[I]:= S:= M[I] / R; IW:= IW + W1; ELMVEC(IW, IW + W2, KK1 - IW, A, A, - S); B[K + I]:= B[K + I] - B[K] * S "END" "END"; AUX[3]:= N; AUX[5]:= MIN; KK:= (N + 1) * W - W1; W2:= - 1; SHIFT:= N * W1; "FOR" K:= N "STEP" - 1 "UNTIL" 1 "DO" "BEGIN" KK:= KK - W; SHIFT:= SHIFT - W1; "IF" W2 < W1 "THEN" W2:= W2 + 1; B[K]:= (B[K] - VECVEC(K + 1, K + W2, SHIFT, B, A)) / A[KK] "END"; END: AUX[1]:= SDET "END" DECSOLBND; "EOP"