1SECTION : 0.0 (JULY 1979) PAGE 0 1SECTION : 0.0 (JULY 1979) PAGE 0 1SECTION : 0.0 (JULY 1979) PAGE 0 1SECTION : 0.0 (JULY 1979) PAGE 0 1SECTION : 1.5.2 (JANUARY 1976) PAGE 1 AUTHORS: J. WOLLESWINKEL AND D. WINTER. CONTRIBUTOR: J. C. P. BUS. INSTITUTE: MATHEMATICAL CENTRE. RECEIVED: 750501. BRIEF DESCRIPTION: THIS SECTION CONTAINS FOURTEEN PROCEDURES FOR CALCULATING, WITH DUOBLE LENGTH ARITHMETIC, THE (SCALAR) PRODUCTS OF SINGLE LENGTH VECTORS AND MATRICES. LNGVECVEC: CALCULATES THE SCALAR PRODUCT OF TWO VECTORS GIVEN IN ONE-DIMENSIONAL ARRAYS; LNGMATVEC: CALCULATES THE SCALAR PRODUCT OF A VECTOR GIVEN IN A ONE-DIMENSIONAL ARRAY AND A ROW OF A MATRIX GIVEN IN A TWO DIMENSIONAL ARRAY; LNGTAMVEC: CALCULATES THE SCALAR PRODUCT OF A VECTOR GIVEN IN A ONE-DIMENSIONAL ARRAY AND A COLUMN OF A MATRIX GIVEN IN A TWO-DIMENSIONAL ARRAY; LNGMATMAT: CALCULATES THE SCALAR PRODUCT OF A ROW OF A MATRIX AND A COLUMN OF ANOTHER MATRIX, WHICH ARE BOTH GIVEN IN TWO-DIMENSIONAL ARRAYS; LNGTAMMAT: CALCULATES THE SCALAR PRODUCT OF COLUMNS OF MATRICES, WHICH ARE BOTH GIVEN IN TWO-DIMENSIONAL ARRAYS; LNGMATTAM: CALCULATES THE SCALAR PRODUCT OF ROWS OF MATRICES, WHICH ARE BOTH GIVEN IN TWO-DIMENSIONAL ARRAYS; LNGSEQVEC: CALCULATES THE SCALAR PRODUCT OF TWO VECTORS GIVEN IN ONE-DIMENSIONAL ARRAYS, WHERE THE MUTUAL SPACINGS BETWEEN THE INDICES OF THE FIRST VECTOR CHANGE LINEARLY; LNGSCAPRD1: CALCULATES THE SCALAR PRODUCT OF TWO VECTORS GIVEN IN ONE-DIMENSIONAL ARRAYS, WHERE THE SPACINGS OF BOTH VECTORS ARE CONSTANT; LNGSYMMATVEC: CALCULATES THE SCALAR PRODUCT OF A VECTOR GIVEN IN A ONE-DIMENSIONAL ARRAY AND A ROW OF A SYMMETRIC MATRIX, WHOSE UPPER TRIANGLE IS STORED COLUMN-WISE IN A ONE- DIMENSIONAL ARRAY; THE ABOVE PROCEDURES HAVE THE POSSIBILITY OF ADDING A DOUBLE LENGTH SCALAR TO THE CALCULATED SCALAR PRODUCT; 1SECTION : 1.5.2 (JANUARY 1976) PAGE 2 LNGFULMATVEC: CALCULATES THE VECTOR A * B, WHERE A IS A GIVEN MATRIX AND B IS A VECTOR; LNGFULTAMVEC: CALCULATES THE VECTOR A' * B, WHERE A' IS THE TRANSPOSED OF THE MATRIX A AND B IS A VECTOR; LNGFULSYMMATVEC: CALCULATES THE VECTOR A * B, WHERE A IS A SYMMETRIC MATRIX, WHOSE UPPERTRIANGLE IS STORED COLUMNWISE IN A ONE-DIMENSIONAL ARRAY, AND B IS A VECTOR; LNGRESVEC: CALCULATES THE RESIDUAL VECTOR A * B + X * C, WHERE A IS A GIVEN MATRIX, B AND C ARE VECTORS AND X IS A SCALAR; LNGSYMRESVEC: CALCULATES THE RESIDUAL VECTOR A * B + X * C, WHERE A IS A SYMMETRIC MATRIX, WHOSE UPPERTRIANGLE IS STORED COLUMNWISE IN A ONE-DIMENSIONAL ARRAY, B AND C ARE VECTORS AND, X IS A SCALAR. IN THIS SECTION (X, XX) DENOTES A DOUBLE-LENGTH NUMBER WITH HEAD X AND TAIL XX (SEE METHOD AND PERFORMANCE). KEYWORDS: ELEMENTARY OPERATIONS, VECTOR OPERATIONS, SCALAR PRODUCTS, DOUBLE-LENGTH ARITHMETIC. 1SECTION : 1.5.2 (JANUARY 1976) PAGE 3 SUBSECTION: LNGVECVEC. CALLING SEQUENCE: THE HEADING OF THE PROCEDURE READS: "PROCEDURE" LNGVECVEC(L, U, SHIFT, A, B, C, CC, D, DD); "VALUE" L, U, SHIFT, C, CC; "INTEGER" L, U, SHIFT; "REAL" C, CC, D, DD; "ARRAY" A, B; "CODE" 34410; THE MEANING OF THE FORMAL PARAMETERS IS: L: ; THE LOWER BOUND OF THE RUNNING SUBSCRIPT; U: ; THE UPPER BOUND OF THE RUNNING SUBSCRIPT; SHIFT: ; THE INDEX-SHIFTING PARAMETER OF THE VECTOR GIVEN IN B; A: ; ONE OF THE VECTORS SHOULD BE GIVEN IN ARRAY A[L:U]; B: ; THE OTHER VECTOR SHOULD BE GIVEN IN ARRAY B[L + SHIFT : U + SHIFT]; C, CC: ; THE HEAD AND TAIL OF THE DOUBLE-LENGTH SCALAR THAT HAS TO BE ADDED TO THE CALCULATED SCALAR PRODUCT; IF CC IS NOT A TAIL TO C THEN AN ERROR MESSAGE WILL BE PRINTED (SEE METHOD AND PERFORMANCE); D, DD: ; EXIT: THE HEAD AND TAIL OF THE CALCULATED DOUBLE-LENGTH RESULT. DATA AND RESULTS: (D, DD):= (C, CC) + THE SCALAR PRODUCT OF THE VECTORS GIVEN IN THE ARRAYS A[L:U] AND B[L + SHIFT : U + SHIFT]. LANGUAGE: COMPASS. 1SECTION : 1.5.2 (JANUARY 1976) PAGE 4 SUBSECTION: LNGMATVEC. CALLING SEQUENCE: THE HEADING OF THE PROCEDURE READS: "PROCEDURE" LNGMATVEC(L, U, I, A, B, C, CC, D, DD); "VALUE" L, U, I, C, CC; "INTEGER" L, U, I; "REAL" C, CC, D, DD; "ARRAY" A, B; "CODE" 34411; THE MEANING OF THE FORMAL PARAMETERS IS: L: ; THE LOWER BOUND OF THE RUNNING SUBSCRIPT; U: ; THE UPPER BOUND OF THE RUNNING SUBSCRIPT; I: ; THE INDEX OF THE ROW-VECTOR GIVEN IN ARRAY A; A: ; THE ROW-VECTOR SHOULD BE GIVEN IN ARRAY A[I:I, L:U]; B: ; THE VECTOR SHOULD BE GIVEN IN ARRAY B[L:U]; C, CC: ; THE HEAD AND TAIL OF THE DOUBLE-LENGTH SCALAR THAT HAS TO BE ADDED TO THE CALCULATED SCALAR PRODUCT; IF CC IS NOT A TAIL TO C THEN AN ERROR MESSAGE WILL BE PRINTED (SEE METHOD AND PERFORMANCE); D, DD: ; EXIT: THE HEAD AND TAIL OF THE CALCULATED DOUBLE-LENGTH RESULT. DATA AND RESULTS: (D, DD):= (C, CC) + THE SCALAR PRODUCT OF THE VECTORS GIVEN IN THE ARRAYS A[I:I, L:U] AND B[L:U]. LANGUAGE: COMPASS. 1SECTION : 1.5.2 (JANUARY 1976) PAGE 5 SUBSECTION: LNGTAMVEC. CALLING SEQUENCE: THE HEADING OF THE PROCEDURE READS: "PROCEDURE" LNGTAMVEC(L, U, I, A, B, C, CC, D, DD); "VALUE" L, U, I, C, CC; "INTEGER" L, U, I; "REAL" C, CC, D, DD; "ARRAY" A, B; "CODE" 34412; THE MEANING OF THE FORMAL PARAMETERS IS: L: ; THE LOWER BOUND OF THE RUNNING SUBSCRIPT; U: ; THE UPPER BOUND OF THE RUNNING SUBSCRIPT; I: ; THE INDEX OF THE COLUMN-VECTOR GIVEN IN ARRAY A; A: ; THE COLUMN-VECTOR SHOULD BE GIVEN IN ARRAY A[L:U, I:I]; B: ; THE VECTOR SHOULD BE GIVEN IN ARRAY B[L:U]; C, CC: ; THE HEAD AND TAIL OF THE DOUBLE-LENGTH SCALAR THAT HAS TO BE ADDED TO THE CALCULATED SCALAR PRODUCT; IF CC IS NOT A TAIL TO C THEN AN ERROR MESSAGE WILL BE PRINTED (SEE METHOD AND PERFORMANCE); D, DD: ; EXIT: THE HEAD AND TAIL OF THE CALCULATED DOUBLE-LENGTH RESULT. DATA AND RESULTS: (D, DD):= (C, CC) + THE SCALAR PRODUCT OF THE VECTORS GIVEN IN THE ARRAYS A[L:U, I:I] AND B[L:U]. LANGUAGE: COMPASS. 1SECTION : 1.5.2 (JANUARY 1976) PAGE 6 SUBSECTION: LNGMATMAT. CALLING SEQUENCE: THE HEADING OF THE PROCEDURE READS: "PROCEDURE" LNGMATMAT(L, U, I, J, A, B, C, CC, D, DD); "VALUE" L, U, I, J, C, CC; "INTEGER" L, U, I, J; "REAL" C, CC, D, DD; "ARRAY" A, B; "CODE" 34413; THE MEANING OF THE FORMAL PARAMETERS IS: L: ; THE LOWER BOUND OF THE RUNNING SUBSCRIPT; U: ; THE UPPER BOUND OF THE RUNNING SUBSCRIPT; I: ; THE INDEX OF THE ROW-VECTOR GIVEN IN ARRAY A; J: ; THE INDEX OF THE COLUMN-VECTOR GIVEN IN ARRAY B; A: ; THE ROW-VECTOR SHOULD BE GIVEN IN ARRAY A[I:I, L:U]; B: ; THE COLUMN-VECTOR SHOULD BE GIVEN IN ARRAY B[L:U, J:J]; C, CC: ; THE HEAD AND TAIL OF THE DOUBLE-LENGTH SCALAR THAT HAS TO BE ADDED TO THE CALCULATED SCALAR PRODUCT; IF CC IS NOT A TAIL TO C THEN AN ERROR MESSAGE WILL BE PRINTED (SEE METHOD AND PERFORMANCE); D, DD: ; EXIT: THE HEAD AND TAIL OF THE CALCULATED DOUBLE-LENGTH RESULT. DATA AND RESULTS: (D, DD):= (C, CC) + THE SCALAR PRODUCT OF THE VECTORS GIVEN IN THE ARRAYS A[I:I, L:U] AND B[L:U, J:J]. LANGUAGE: COMPASS. 1SECTION : 1.5.2 (JANUARY 1976) PAGE 7 SUBSECTION: LNGTAMMAT. CALLING SEQUENCE: THE HEADING OF THE PROCEDURE READS: "PROCEDURE" LNGTAMMAT(L, U, I, J, A, B, C, CC, D, DD); "VALUE" L, U, I, J, C, CC; "INTEGER" L, U, I, J; "REAL" C, CC, D, DD; "ARRAY" A, B; "CODE" 34414; THE MEANING OF THE FORMAL PARAMETERS IS: L: ; THE LOWER BOUND OF THE RUNNING SUBSCRIPT; U: ; THE UPPER BOUND OF THE RUNNING SUBSCRIPT; I: ; THE INDEX OF THE COLUMN-VECTOR GIVEN IN ARRAY A; J: ; THE INDEX OF THE COLUMN-VECTOR GIVEN IN ARRAY B; A: ; ONE OF THE COLUMN-VECTORS SHOULD BE GIVEN IN ARRAY A[L:U, I:I]; B: ; THE OTHER COLUMN-VECTOR SHOULD BE GIVEN IN ARRAY B[L:U, J:J]; C, CC: ; THE HEAD AND TAIL OF THE DOUBLE-LENGTH SCALAR THAT HAS TO BE ADDED TO THE CALCULATED SCALAR PRODUCT; IF CC IS NOT A TAIL TO C THEN AN ERROR MESSAGE WILL BE PRINTED (SEE METHOD AND PERFORMANCE); D, DD: ; EXIT: THE HEAD AND TAIL OF THE CALCULATED DOUBLE-LENGTH RESULT. DATA AND RESULTS: (D, DD):= (C, CC) + THE SCALAR PRODUCT OF THE VECTORS GIVEN IN THE ARRAYS A[L:U, I:I] AND B[L:U, J:J]. LANGUAGE: COMPASS. 1SECTION : 1.5.2 (JANUARY 1976) PAGE 8 SUBSECTION: LNGMATTAM. CALLING SEQUENCE: THE HEADING OF THE PROCEDURE READS: "PROCEDURE" LNGMATTAM(L, U, I, J, A, B, C, CC, D, DD); "VALUE" L, U, I, J, C, CC; "INTEGER" L, U, I, J; "REAL" C, CC, D, DD; "ARRAY" A, B; "CODE" 34415; THE MEANING OF THE FORMAL PARAMETERS IS: L: ; THE LOWER BOUND OF THE RUNNING SUBSCRIPT; U: ; THE UPPER BOUND OF THE RUNNING SUBSCRIPT; I: ; THE INDEX OF THE ROW-VECTOR GIVEN IN ARRAY A; J: ; THE INDEX OF THE ROW-VECTOR GIVEN IN ARRAY B; A: ; ONE OF THE ROW-VECTORS SHOULD BE GIVEN IN ARRAY A[I:I,L:U]; B: ; THE OTHER ROW-VECTOR SHOULD BE GIVEN IN ARRAY B[J:J, L:U]; C, CC: ; THE HEAD AND TAIL OF THE DOUBLE-LENGTH SCALAR THAT HAS TO BE ADDED TO THE CALCULATED SCALAR PRODUCT; IF CC IS NOT A TAIL TO C THEN AN ERROR MESSAGE WILL BE PRINTED (SEE METHOD AND PERFORMANCE); D, DD: ; EXIT: THE HEAD AND TAIL OF THE CALCULATED DOUBLE-LENGTH RESULT. DATA AND RESULTS: (D, DD):= (C, CC) + THE SCALAR PRODUCT OF THE VECTORS GIVEN IN THE ARRAYS A[I:I, L:U] AND B[J:J, L:U]. LANGUAGE: COMPASS. 1SECTION : 1.5.2 (JANUARY 1976) PAGE 9 SUBSECTION: LNGSEQVEC. CALLING SEQUENCE: THE HEADING OF THE PROCEDURE READS: "PROCEDURE" LNGSEQVEC(L, U, IL, SHIFT, A, B, C, CC, D, DD); "VALUE" L, U, IL, SHIFT, C, CC; "INTEGER" L, U, IL, SHIFT; "REAL" C, CC, D, DD; "ARRAY" A, B; "CODE" 34416; THE MEANING OF THE FORMAL PARAMETERS IS: L: ; THE LOWER BOUND OF THE RUNNING SUBSCRIPT; U: ; THE UPPER BOUND OF THE RUNNING SUBSCRIPT; IL: ; THE FIRST INDEX OF THE VECTOR GIVEN IN ARRAY A; SHIFT: ; THE INDEX-SHIFTING PARAMETER OF THE VECTOR GIVEN IN B; A: ; ONE OF THE VECTORS SHOULD BE GIVEN IN ARRAY A[P:Q], WHERE P = MIN(IL + (J + L - 1) * (J - 1) // 2; J = L, ..., U) AND Q = MAX(IL + (J + L - 1) * (J - 1) // 2; J = L, ..., U); B: ; THE OTHER VECTOR SHOULD BE GIVEN IN ARRAY B[L + SHIFT : U + SHIFT]; C, CC: ; THE HEAD AND TAIL OF THE DOUBLE-LENGTH SCALAR THAT HAS TO BE ADDED TO THE CALCULATED SCALAR PRODUCT; IF CC IS NOT A TAIL TO C THEN AN ERROR MESSAGE WILL BE PRINTED (SEE METHOD AND PERFORMANCE); D, DD: ; EXIT: THE HEAD AND TAIL OF THE CALCULATED DOUBLE-LENGTH RESULT. DATA AND RESULTS: (D, DD):= (C, CC) + THE SCALAR PRODUCT OF THE VECTORS GIVEN IN THE ARRAYS A[P:Q] (SEE THE MEANING OF PARAMETER A ) AND B[L + SHIFT: U + SHIFT], WHERE THE ELEMENTS OF THE FIRST VECTOR ARE A[IL + (J + L - 1) * (J - 1) // 2], J = L, ..., U. LANGUAGE: COMPASS. 1SECTION : 1.5.2 (JANUARY 1976) PAGE 10 SUBSECTION: LNGSCAPRD1. CALLING SEQUENCE: THE HEADING OF THE PROCEDURE READS: "PROCEDURE" LNGSCAPRD1(LA, SA, LB, SB, N, A, B, C, CC, D, DD); "VALUE" LA, SA, LB, SB, N, C, CC; "INTEGER" LA, SA, LB, SB, N; "REAL" C, CC, D, DD; "ARRAY" A, B; "CODE" 34417; THE MEANING OF THE FORMAL PARAMETERS IS: LA: ; THE FIRST INDEX OF THE VECTOR GIVEN IN ARRAY A; SA: ; THE SPACING OF THE VECTOR GIVEN IN ARRAY A; LB: ; THE FIRST INDEX OF THE VECTOR GIVEN IN ARRAY B; SB: ; THE SPACING OF THE VECTOR GIVEN IN ARRAY B; N: ; THE UPPER BOUND OF THE RUNNING SUBSCRIPT; IF N < 1, THEN ON EXIT THE RESULT (D, DD) WILL SATISFY: (D, DD) = (C, CC); A: ; ONE OF THE VECTORS SHOULD BE GIVEN IN ARRAY A[MIN(LA, LA + (N - 1) * SA) : MAX(LA, LA + (N - 1) * SA)]; B: ; THE OTHER VECTOR SHOULD BE GIVEN IN ARRAY B[MIN(LB, LB + (N - 1) * SB) : MAX(LB, LB + (N - 1) * SB)]; C, CC: ; THE HEAD AND TAIL OF THE DOUBLE-LENGTH SCALAR THAT HAS TO BE ADDED TO THE CALCULATED SCALAR PRODUCT; IF CC IS NOT A TAIL TO C THEN AN ERROR MESSAGE WILL BE PRINTED (SEE METHOD AND PERFORMANCE); D, DD: ; EXIT: THE HEAD AND TAIL OF THE CALCULATED DOUBLE-LENGTH RESULT. DATA AND RESULTS: (D, DD):= (C, CC) + THE SCALAR PRODUCT OF THE VECTORS GIVEN IN THE ARRAYS A[MIN(LA, LA + (N - 1) * SA):MAX(LA, LA + (N - 1) * SA)] AND B[MIN(LB, LB + (N - 1) * SB) : MAX(LB, LB + (N - 1) * SB)], WHERE THE ELEMENTS OF THE VECTORS ARE A[LA + (J - 1) * SA] AND B[LB + (J - 1) * SB], FOR J = 1, .., N. LANGUAGE: COMPASS. 1SECTION : 1.5.2 (JANUARY 1976) PAGE 11 SUBSECTION: LNGSYMMATVEC. CALLING SEQUENCE: THE HEADING OF THE PROCEDURE READS: "PROCEDURE" LNGSYMMATVEC(L, U, I, A, B, C, CC, D, DD); "VALUE" L, U, I, C, CC; "INTEGER" L, U, I; "REAL" C, CC, D, DD; "ARRAY" A, B; "CODE" 34418; THE MEANING OF THE FORMAL PARAMETERS IS: L: ; THE LOWER BOUND OF THE RUNNING SUBSCRIPT; U: ; THE UPPER BOUND OF THE RUNNING SUBSCRIPT; I: ; THE INDEX OF THE ROW OF THE SYMMETRIC MATRIX, WHOSE UPPER TRIANGLE IS STORED COLUMN-WISE IN THE ONE-DIMENSIONAL ARRAY A; A: ; THE ROW OF THE SYMMETRIC MATRIX SHOULD BE GIVEN IN ARRAY A[P:Q], WHERE, IF I > L THEN P = (I - 1) * I // 2 + L ELSE P = (L - 1) * L // 2 + I, AND IF I > U THEN Q = (I - 1) * I // 2 + U ELSE Q = (U - 1) * U // 2 + I; B: ; THE VECTOR SHOULD BE GIVEN IN ARRAY B[L:U]; C, CC: ; THE HEAD AND TAIL OF THE DOUBLE-LENGTH SCALAR THAT HAS TO BE ADDED TO THE CALCULATED SCALAR PRODUCT; IF CC IS NOT A TAIL TO C THEN AN ERROR MESSAGE WILL BE PRINTED (SEE METHOD AND PERFORMANCE); D, DD: ; EXIT: THE HEAD AND TAIL OF THE CALCULATED DOUBLE-LENGTH RESULT. PROCEDURES USED: DPMUL = CP31103; LNGADD = CP31105. DATA AND RESULTS: (D, DD):= (C, CC) + THE SCALAR PRODUCT OF THE VECTORS GIVEN IN THE ARRAYS A[P:Q] (SEE THE MEANING OF PARAMETER A) AND B[L:U], WHERE THE ELEMENTS OF THE FIRST VECTOR ARE: IF L < I THEN A[(I - 1) * I // 2 + J], J = L, ..., MIN(U, I - 1) AND A[(J - 1) * J // 2 + I], J = I, ..., U, RESPECTIVELY, OTHERWISE A[(J - 1) * J // 2 + I], J = L, ..., U; THE VALUES OF L, U, I SHOULD BE POSITIVE. LANGUAGE: ALGOL 60. 1SECTION : 1.5.2 (JANUARY 1976) PAGE 12 SUBSECTION: LNGFULMATVEC. CALLING SEQUENCE: THE HEADING OF THE PROCEDURE READS: "PROCEDURE" LNGFULMATVEC(LR, UR, LC, UC, A, B, C); "VALUE" LR, UR, LC, UC, B; "INTEGER" LR, UR, LC, UC; "ARRAY" A, B, C; "CODE" 31505; THE MEANING OF THE FORMAL PARAMETERS IS: LR, UR: ; LOWER AND UPPER BOUND OF THE ROW-INDEX; LC, UC: ; LOWER AND UPPER BOUND OF THE COLUMN-INDEX; A: ; "ARRAY" A[LR:UR,LC:UC]; THE MATRIX; B: ; "ARRAY" B[LC:UC]; THE VECTOR; C: ; "ARRAY" C[LR:UR]; THE RESULT A * B, CALCULATED WITH DOUBLE LENGTH ARITHMETIC, IS DELIVERED IN C. LANGUAGE: COMPASS. METHOD AND PERFORMANCE: ELEMENTARY. 1SECTION : 1.5.2 (JANUARY 1976) PAGE 13 SUBSECTION: LNGFULTAMVEC. CALLING SEQUENCE: THE HEADING OF THE PROCEDURE READS: "PROCEDURE" LNGFULTAMVEC(LR, UR, LC, UC, A, B, C); "VALUE" LR, UR, LC, UC, B; "INTEGER" LR, UR, LC, UC; "ARRAY" A, B, C; "CODE" 31506; THE MEANING OF THE FORMAL PARAMETERS IS: LR, UR: ; LOWER AND UPPER BOUND OF THE ROW-INDEX; LC, UC: ; LOWER AND UPPER BOUND OF THE COLUMN-INDEX; A: ; "ARRAY" A[LR:UR,LC:UC]; THE MATRIX; B: ; "ARRAY" B[LR:UR]; THE VECTOR; C: ; "ARRAY" C[LC:UC]; THE RESULT A' * B, CALCULATED WITH DOUBLE LENGTH ARITHMETIC, IS DELIVERED IN C; HERE A' DENOTES THE TRANSPOSED OF THE MATRIX A. LANGUAGE: COMPASS. METHOD AND PERFORMANCE: ELEMENTARY. 1SECTION : 1.5.2 (JANUARY 1976) PAGE 14 SUBSECTION: LNGFULSYMMATVEC. CALLING SEQUENCE: THE HEADING OF THE PROCEDURE READS: "PROCEDURE" LNGFULSYMMATVEC(LR, UR, LC, UC, A, B, C); "VALUE" LR, UR, LC, UC, B; "INTEGER" LR, UR, LC, UC; "ARRAY" A, B, C; "CODE" 31507; THE MEANING OF THE FORMAL PARAMETERS IS: LR, UR: ; LOWER AND UPPER BOUND OF THE ROW-INDEX; LR >= 1; LC, UC: ; LOWER AND UPPER BOUND OF THE COLUMN-INDEX; LC >= 1; A: ; "ARRAY" A[L:U], WHERE: L = MIN(LR * (LR - 1) // 2 + LC, LC * (LC - 1) // 2 + LR), U = MAX(UR * (UR - 1) // 2 + UC, UC * (UC - 1) // 2 + UR) AND THE (I,J)-TH ELEMENT OF THE SYMMETRIC MATRIX SHOULD BE GIVEN IN A[J * (J - 1) // 2 + I]; B: ; "ARRAY" B[LC:UC]; THE VECTOR; C: ; "ARRAY" C[LR:UR]; THE RESULT A * B, CALCULATED WITH DOUBLE LENGTH ARITHMETIC, IS DELIVERED IN C. PROCEDURES USED: LNGSYMMATVEC = CP34418. LANGUAGE: ALGOL 60. METHOD AND PERFORMANCE: ELEMENTARY. 1SECTION : 1.5.2 (JANUARY 1976) PAGE 15 SUBSECTION: LNGRESVEC. CALLING SEQUENCE: THE HEADING OF THE PROCEDURE READS: "PROCEDURE" LNGRESVEC(LR, UR, LC, UC, A, B, C, X); "VALUE" LR, UR, LC, UC, B, X; "INTEGER" LR, UR, LC, UC; "REAL" X; "ARRAY" A, B, C; "CODE" 31508; THE MEANING OF THE FORMAL PARAMETERS IS: LR, UR: ; LOWER AND UPPER BOUND OF THE ROW-INDEX; LC, UC: ; LOWER AND UPPER BOUND OF THE COLUMN-INDEX; A: ; "ARRAY" A[LR:UR,LC:UC]; THE MATRIX; B: ; "ARRAY" B[LC:UC]; THE VECTOR; X: ; THE VALUE OF THE MULTIPLYING SCALAR; C: ; "ARRAY" C[LR:UR]; THE RESULT A * B + X * C, CALCULATED WITH DOUBLE LENGTH ARITHMETIC, IS OVERWRITTEN ON C. LANGUAGE: COMPASS. METHOD AND PERFORMANCE: ELEMENTARY. 1SECTION : 1.5.2 (JANUARY 1976) PAGE 16 SUBSECTION: LNGSYMRESVEC. CALLING SEQUENCE: THE HEADING OF THE PROCEDURE READS: "PROCEDURE" LNGSYMRESVEC(LR, UR, LC, UC, A, B, C, X); "VALUE" LR, UR, LC, UC, B, X; "INTEGER" LR, UR, LC, UC; "REAL" X; "ARRAY" A, B, C; "CODE" 31509; THE MEANING OF THE FORMAL PARAMETERS IS: LR, UR: ; LOWER AND UPPER BOUND OF THE ROW-INDEX; LR >= 1; LC, UC: ; LOWER AND UPPER BOUND OF THE COLUMN-INDEX; LC >= 1; A: ; "ARRAY" A[L:U], WHERE: L = MIN(LR * (LR - 1) // 2 + LC, LC * (LC - 1) // 2 + LR), U = MAX(UR * (UR - 1) // 2 + UC, UC * (UC - 1) // 2 + UR) AND THE (I,J)-TH ELEMENT OF THE SYMMETRIC MATRIX SHOULD BE GIVEN IN A[J * (J - 1) // 2 + I]; B: ; "ARRAY" B[LC:UC]; THE VECTOR; X: ; THE VALUE OF THE MULTIPLYING SCALAR; C: ; "ARRAY" C[LR:UR]; THE RESULT A * B + X * C, CALCULATED WITH DOUBLE LENGTH ARITHMETIC, IS OVERWRITTEN ON C. PROCEDURES USED: DPMUL = CP31103; LNGSYMMATVEC = CP34418. LANGUAGE: ALGOL 60. METHOD AND PERFORMANCE: ELEMENTARY. 1SECTION : 1.5.2 (JANUARY 1976) PAGE 17 METHOD AND PERFORMANCE: ALL PROCEDURES GIVEN IN THIS SECTION MAKE USE OF THE DOUBLE-LENGTH ARITHMETIC OPERATIONS AVAILABLE IN HARDWARE ON THE CYBER 73. LET (X, XX) DENOTE A DOUBLE-LENGTH NUMBER, THEN WE WILL SAY THAT XX IS A TAIL TO X, WHEN THE FOLLOWING CONDITIONS HOLD: X = FL1(X + XX), (X, XX) = FL2(X + XX), WHERE FL1(. + .) AND FL2(. + .) DENOTE SINGLE-LENGTH RESPECTIVELY DOUBLE-LENGTH ADDITION WITH TRUNCATING OF THE RESULT TO 48 AND 96 BITS RESPECTIVELY. WHEN A PROCEDURE DELIVERS A DOUBLE LENGTH RESULT IN D AND DD, THEN THESE RESULTS ARE SUCH THAT DD IS A TAIL TO D; WHEN ONE SHOULD PROVIDE AN INITIALIZING DOUBLE LENGTH SCALAR IN C AND CC, THEN CC SHOULD BE A TAIL TO C, OTHERWISE THE FOLLOWING ERROR MESSAGE WILL BE PRINTED: DP PARAMETER TAIL ERROR AND EXECUTION OF THE PROGRAM WILL TERMINATE IN THE USUAL WAY. NOTE THAT CC = 0 IS A TAIL TO C FOR ALL VALUES OF C. FURTHERMORE, IT SEEMS WORTHWHILE TO NOTE THAT THE ARRAY B MUST BE A VALUE PARAMETER IN ALGOL 60, HOWEVER, IN THE COMPASS ROUTINE THE DUPLICATION OF THIS ARRAY IS ONLY DONE IF NECESSARY, I.E. IF THE ACTUAL PARAMETERS B AND C ARE THE SAME. SOURCE TEXTS: THE PROCEDURES IN THIS SECTION ARE WRITTEN IN COMPASS, EXCEPT FOR LNGSYMMATVEC, LNGFULSYMMATVEC AND LNGSYMRESVEC. WE GIVE EQUIVALENT ALGOL 60 TEXTS OF THE COMPASS ROUTINES. FOR THE COMPASS TEXT SEE APPENDIX C, SECTION 1.5.2. 0"CODE" 34410; "PROCEDURE" LNGVECVEC(L, U, SHIFT, A, B, C, CC, D, DD); "VALUE" L, U, SHIFT, C, CC; "INTEGER" L,U,SHIFT; "REAL" C, CC, D, DD; "ARRAY" A, B; "BEGIN" "REAL" E, EE; "PROCEDURE" DPMUL(A, B, C, CC); "CODE" 31103; "PROCEDURE" LNGADD(A, AA, B, BB, C, CC); "CODE" 31105; "FOR" L:= L "STEP" 1 "UNTIL" U "DO" "BEGIN" DPMUL(A[L], B[L + SHIFT], E, EE); LNGADD(C, CC, E, EE, C, CC) "END"; D:= C; DD:= CC "END" LNGVECVEC; "EOP" 1SECTION : 1.5.2 (JANUARY 1976) PAGE 18 0"CODE" 34411; "PROCEDURE" LNGMATVEC(L, U, I, A, B, C, CC, D, DD); "VALUE" L, U, I, C, CC; "INTEGER" L, U, I; "REAL" C, CC, D, DD; "ARRAY" A, B; "BEGIN" "REAL" E, EE; "PROCEDURE" DPMUL(A, B, C, CC); "CODE" 31103; "PROCEDURE" LNGADD(A, AA, B, BB, C, CC); "CODE" 31105; "FOR" L:= L "STEP" 1 "UNTIL" U "DO" "BEGIN" DPMUL(A[I,L], B[I], E, EE); LNGADD(C, CC, E, EE, C, CC) "END"; D:= C; DD:= CC "END" LNGMATVEC; "EOP" 0"CODE" 34412; "PROCEDURE" LNGTAMVEC(L, U, I, A, B, C, CC, D, DD); "VALUE" L, U, I, C, CC; "INTEGER" L, U, I; "REAL" C, CC, D, DD; "ARRAY" A, B; "BEGIN" "REAL" E, EE; "PROCEDURE" DPMUL(A, B, C, CC); "CODE" 31103; "PROCEDURE" LNGADD(A, AA, B, BB, C, CC); "CODE" 31105; "FOR" L:= L "STEP" 1 "UNTIL" U "DO" "BEGIN" DPMUL(A[L,I], B[I], E, EE); LNGADD(C, CC, E, EE, C, CC) "END"; D:= C; DD:= CC "END" LNGTAMVEC; "EOP" 0"CODE" 34413; "PROCEDURE" LNGMATMAT(L, U, I, J, A, B, C, CC, D, DD); "VALUE" L, U, I, J, C, CC; "INTEGER" L, U, I, J; "REAL" C, CC, D, DD; "ARRAY" A, B; "BEGIN" "REAL" E, EE; "PROCEDURE" DPMUL(A, B, C, CC); "CODE" 31103; "PROCEDURE" LNGADD(A, AA, B, BB, C, CC); "CODE" 31105; "FOR" L:= L "STEP" 1 "UNTIL" U "DO" "BEGIN" DPMUL(A[I,L], B[L,J], E, EE); LNGADD(C, CC, E, EE, C, CC) "END"; D:= C; DD:= CC "END" LNGMATMAT; "EOP" 1SECTION : 1.5.2 (JANUARY 1976) PAGE 19 0"CODE" 34414; "PROCEDURE" LNGTAMMAT(L, U, I, J, A, B, C, CC, D, DD); "VALUE" L, U, I, J, C, CC; "INTEGER" L, U, I, J; "REAL" C, CC, D, DD; "ARRAY" A, B; "BEGIN" "REAL" E, EE; "PROCEDURE" DPMUL(A, B, C, CC); "CODE" 31103; "PROCEDURE" LNGADD(A, AA, B, BB, C, CC); "CODE" 31105; "FOR" L:= L "STEP" 1 "UNTIL" U "DO" "BEGIN" DPMUL(A[L,I], B[L,J], E, EE); LNGADD(C, CC, E, EE, C, CC) "END"; D:= C; DD:= CC "END" LNGTAMMAT; "EOP" 0"CODE" 34415; "PROCEDURE" LNGMATTAM(L, U, I, J, A, B, C, CC, D, DD); "VALUE" L, U, I, J, C, CC; "INTEGER" L, U, I, J; "REAL" C, CC, D, DD; "ARRAY" A, B; "BEGIN" "REAL" E, EE; "PROCEDURE" DPMUL(A, B, C, CC); "CODE" 31103; "PROCEDURE" LNGADD(A, AA, B, BB, C, CC); "CODE" 31105; "FOR" L:= L "STEP" 1 "UNTIL" U "DO" "BEGIN" DPMUL(A[I,L], B[J,L], E, EE); LNGADD(C, CC, E, EE, C, CC) "END"; D:= C; DD:= CC "END" LNGMATTAM; "EOP" 0"CODE" 34416; "PROCEDURE" LNGSEQVEC(L, U, IL, SHIFT, A, B, C, CC, D, DD); "VALUE" L, U, IL, SHIFT, C, CC; "INTEGER" L, U, IL, SHIFT; "REAL" C, CC, D, DD; "ARRAY" A, B; "BEGIN" "REAL" E, EE; "PROCEDURE" DPMUL(A, B, C, CC); "CODE" 31103; "PROCEDURE" LNGADD(A, AA, B, BB, C, CC); "CODE" 31105; "FOR" L:= L "STEP" 1 "UNTIL" U "DO" "BEGIN" DPMUL(A[IL], B[L + SHIFT], E, EE); IL:= IL + L; LNGADD(C, CC, E, EE, C, CC) "END"; D:= C; DD:= CC "END" LNGSEQVEC; "EOP" 1SECTION : 1.5.2 (JANUARY 1976) PAGE 20 "CODE" 34417; "PROCEDURE" LNGSCAPRD1(LA, SA, LB, SB, N, A, B, C, CC, D, DD); "VALUE" LA, SA, LB, SB, N, C, CC, D, DD; "INTEGER" LA, SA, LB, SB,N; "REAL" C, CC, D, DD; "ARRAY" A, B; "BEGIN" "REAL" E, EE; "INTEGER" K; "PROCEDURE" DPMUL(A, B, C, CC); "CODE" 31103; "PROCEDURE" LNGADD(A, AA, B, BB, C, CC); "CODE" 31105; "FOR" K:= 1 "STEP" 1 "UNTIL" N "DO" "BEGIN" DPMUL(A[LA], B[LB], E, EE); LA:= LA + SA; LB:= LB + SB; LNGADD(C, CC, E, EE, C, CC) "END"; D:= C; DD:= CC "END"; "EOP" "CODE" 34418; "PROCEDURE" LNGSYMMATVEC(L, U, I, A, B, C, CC, D, DD); "VALUE" L, U, I, C, CC; "INTEGER" L, U, I; "REAL" C, CC, D, DD; "ARRAY" A, B; "BEGIN" "INTEGER" K, M; "PROCEDURE" LNGVECVEC(L, U, S, A, B, C, T, D, R); "CODE" 34410; "PROCEDURE" LNGSEQVEC(L, U, IL, S, A, B, C, T, D, R); "CODE" 34416; M:= "IF" L > I "THEN" L "ELSE" I; K:= M * (M - 1) // 2; LNGVECVEC(L, "IF" I <= U "THEN" I - 1 "ELSE" U, K, B, A, C, CC, C, CC); LNGSEQVEC(M, U, K + I, 0, A, B, C, CC, D, DD) "END" LNGSYMMATVEC; "EOP" "CODE" 31505; "PROCEDURE" LNGFULMATVEC(LR, UR, LC, UC, A, B, C); "VALUE" LR, UR, LC, UC, B; "INTEGER" LR, UR, LC, UC; "ARRAY" A, B, C; "BEGIN" "REAL" D, DD; "PROCEDURE" LNGMATVEC(L, U, I, A, B, C, CC, D, DD); "CODE" 34411; "FOR" LR:= LR "STEP" 1 "UNTIL" UR "DO" "BEGIN" LNGMATVEC(LC, UC, LR, A, B, 0, 0, D, DD); C[LR]:= D + DD "END" "END" LNGFULMATVEC; "EOP" "CODE" 31506; "PROCEDURE" LNGFULTAMVEC(LR, UR, LC, UC, A, B, C); "VALUE" LR, UR, LC, UC, B; "INTEGER" LR, UR, LC, UC; "ARRAY" A, B, C; "BEGIN" "REAL" D, DD; "PROCEDURE" LNGTAMVEC(L, U, I, A, B, C, CC, D, DD); "CODE" 34412; "FOR" LC:= LC "STEP" 1 "UNTIL" UC "DO" "BEGIN" LNGTAMVEC(LR, UR, LC, A, B, 0, 0, D, DD); C[LC]:= D + DD "END" "END" LNGFULTAMVEC; "EOP" 1SECTION : 1.5.2 (JANUARY 1976) PAGE 21 0"CODE" 31507; "PROCEDURE" LNGFULSYMMATVEC(LR, UR, LC, UC, A, B, C); "VALUE" LR, UR, LC, UC, B; "INTEGER" LR, UR, LC, UC; "ARRAY" A, B, C; "BEGIN" "REAL" D, DD; "PROCEDURE" LNGSYMMATVEC(L, U, I, A, B, C, CC, D, DD); "CODE" 34418; "FOR" LR:= LR "STEP" 1 "UNTIL" UR "DO" "BEGIN" LNGSYMMATVEC(LC, UC, LR, A, B, 0, 0, D, DD); C[LR]:= D + DD "END" "END" LNGFULSYMMATVEC; "EOP" 0"CODE" 31508; "PROCEDURE" LNGRESVEC(LR, UR, LC, UC, A, B, C, X); "VALUE" LR, UR, LC, UC, X; "INTEGER" LR, UR, LC, UC; "REAL" X; "ARRAY" A, B, C; "BEGIN" "REAL" D, DD, E, EE; "PROCEDURE" DPMUL(X, Y, E, EE); "CODE" 31103; "PROCEDURE" LNGMATVEC(L, U, I, A, B, C, CC, D, DD); "CODE" 34411; "FOR" LR:= LR "STEP" 1 "UNTIL" UR "DO" "BEGIN" DPMUL(C[LR], X, E, EE); LNGMATVEC(LC, UC, LR, A, B, E, EE, D, DD); C[LR]:= D + DD "END" "END" LNGRESVEC; "EOP" 0"CODE" 31509; "PROCEDURE" LNGSYMRESVEC(LR, UR, LC, UC, A, B, C, X); "VALUE" LR, UR, LC, UC, B, X; "INTEGER" LR, UR, LC, UC; "REAL" X; "ARRAY" A, B, C; "BEGIN" "REAL" D, DD, E, EE; "PROCEDURE" DPMUL(X, Y, E, EE); "CODE" 31103; "PROCEDURE" LNGSYMMATVEC(L, U, I, A, B, C, CC, D, DD); "CODE" 34418; "FOR" LR:= LR "STEP" 1 "UNTIL" UR "DO" "BEGIN" DPMUL(C[LR], X, E, EE); LNGSYMMATVEC(LC, UC, LR, A, B, E, EE, D, DD); C[LR]:= D + DD "END" "END" LNGSYMRESVEC; "EOP" 1SECTION : 1.2.7 (DECEMBER 1979) PAGE 1 AUTHORS : J.J.G. ADMIRAAL, C.G. VAN DER LAAN. CONTRIBUTORS : J.J.G. ADMIRAAL, H. FIOLET, C.G. VAN DER LAAN. INSTITUTE: MATHEMATICAL CENTRE, UNIVERSITY OF AMSTERDAM. RECEIVED : 730817. BRIEF DESCRIPTION : THIS SECTION CONTAINS THE PROCEDURES ROTCOMCOL, ROTCOMROW AND CHSH2 ROTCOMCOL REPLACES THE COLUMN VECTOR VR+I*VI GIVEN IN THE ARRAYS AR,AI[L:U,I:I] AND THE COLUMN VECTOR YR+I*YI GIVEN IN THE ARRAYS AR,AI[L:U,J:J] BY THE VECTORS (VR+I*VI)*(CR-I*CI)-(YR+I*YI)*S AND (YR+I*YI)*(CR+I*CI)+(VR+I*VI)*S, RESPECTIVELY. ROTCOMROW REPLACES THE ROW VECTOR VR+I*VI GIVEN IN THE ARRAYS AR,AI[I:I,L:U] AND THE ROW VECTOR YR+I*YI GIVEN IN THE ARRAYS AR,AI[J:J,L:U] BY THE VECTORS (VR+I*VI)*(CR-I*CI)+(YR+I*YI)*S AND (YR+I*YI)*(CR+I*CI)-(VR+I*VI)*S, RESPECTIVELY. CHSH2 COMPUTES THE COMPLEX HOUSEHOLDERMATRIX THAT MAPS THE COMPLEX VECTOR (A1,A2) INTO THE DIRECTION (1,0). WARNING : IN ROTCOMCOL AND ROTCOMROW THE COSINE IS COMPLEX AND THE SINE IS REAL, IN CONTRAST TO THIS, IN CHSH2 THE SINE IS COMPLEX AND THE COSINE IS REAL. KEYWORDS : COMPLEX VECTOR OPERATIONS, ROTATION, HOUSEHOLDER MATRIX. 1SECTION : 1.2.7 (JANUARY 1976) PAGE 2 SUBSECTION : ROTCOMCOL . CALLING SEQUENCE : THE HEADING OF THE PROCEDURE READS : "PROCEDURE" ROTCOMCOL(L, U, I, J, AR, AI, CR, CI, S); "VALUE" L, U, I, J, CR, CI, S; "INTEGER" L, U, I, J; "REAL" CR, CI, S; "ARRAY" AR, AI; "CODE" 34357; THE MEANING OF THE FORMAL PARAMETERS IS : L,U,I,J: ; THE ROTATION IS PERFORMED ON THE COLUMN VECTORS AR,AI[L:U,I:I] AND AR,AI[L:U,J:J]; AR,AI: ; "ARRAY" AR,AI[L:U,I:J]; ENTRY: AR:THE REAL PARTS OF THE COLUMN VECTORS AI:THE IMAGINARY PARTS OF THE COLUMN VECTORS EXIT: THE RESULTING VECTORS (SEE ALSO BRIEF DESCRIPTION); CR,CI,S: ; ENTRY: ROTATION FACTORS; SEE ALSO BRIEF DESCRIPTION. RUNNING TIME : ROUGHLY PROPORTIONAL TO (U-L) . LANGUAGE: ALGOL 60. 1SECTION : 1.2.7 (JANUARY 1976) PAGE 3 SUBSECTION : ROTCOMROW . CALLING SEQUENCE : THE HEADING OF THE PROCEDURE READS : "PROCEDURE" ROTCOMROW(L, U, I, J, AR, AI, CR, CI, S); "VALUE" L, U, I, J, CR, CI, S; "INTEGER" L, U, I, J; "REAL" CR, CI, S; "ARRAY" AR, AI; "CODE" 34358; THE MEANING OF THE FORMAL PARAMETERS IS : L,U,I,J: ; THE ROTATION IS PERFORMED ON THE ROW VECTORS AR,AI[I:I,L:U] AND AR,AI[J:J,L:U]; AR,AI: ; "ARRAY" AR,AI[I:J,L:U]; ENTRY: AR:THE REAL PARTS OF THE ROW VECTORS AI:THE IMAGINARY PARTS OF THE ROW VECTORS EXIT: THE RESULTING VECTORS (SEE ALSO BRIEF DESCRIPTION); CR,CI,S: ; ENTRY: ROTATION FACTORS; SEE ALSO BRIEF DESCRIPTION. PROCEDURES USED : NONE . RUNNING TIME : ROUGHLY PROPORTIONAL TO (U-L) . LANGUAGE: ALGOL 60. 1SECTION : 1.2.7 (JANUARY 1976) PAGE 4 EXAMPLE OF USE : "BEGIN" "COMMENT" EXAMPLE OF USE ROTCOMCOL; "PROCEDURE" ROTCOMCOL(L,U,I,J,AR,AI,CR,CI,S);"CODE" 34357; "REAL" "ARRAY" AR,AI[1:2,1:2]; "INTEGER" I,J; AR[1,1]:=+4;AR[1,2]:=+5;AR[2,1]:=-5;AR[2,2]:=+4; AI[1,1]:=+3;AI[1,2]:= 0;AI[2,1]:= 0;AI[2,2]:=-3; OUTPUT(61,"(""("INPUT MATRIX:")",/")"); OUTPUT(61,"("-D,+D,"("*I")",4B,-D,Z/,BB-D,Z,3B,-D,+D,"("*I")",/")", AR[1,1],AI[1,1],AR[1,2],AI[1,2],AR[2,1],AI[2,1],AR[2,2],AI[2,2]); OUTPUT(61,"("//,"("AFTER POSTMULTIPLICATION WITH:")",/")"); OUTPUT(61,"(""(".08-.06*I -.1")",/, "(" .1 .08+.06*I")",//")"); ROTCOMCOL(1,2,1,2,AR,AI,.08,.06,-.1); OUTPUT(61,"(""("DELIVERS:")",/")"); OUTPUT(61,"("-D,Z,2BD,Z/,BD,Z,B-D,Z")", AR[1,1],AI[1,1],AR[1,2],AI[1,2],AR[2,1],AI[2,1],AR[2,2],AI[2,2]); "END" OUTPUT: INPUT MATRIX: 4+3*I 5 -5 4-3*I AFTER POSTMULTIPLICATION WITH: .08-.06*I -.1 .1 .08+.06*I DELIVERS: 1 0 0 1 . 1SECTION : 1.2.7 (DECEMBER 1979) PAGE 5 SUBSECTION: CHSH2. CALLING SEQUENCE: THE HEADING OF THE PROCEDURE IS: "PROCEDURE" CHSH2(A1R,A1I,A2R,A2I,C,SR,SI); "VALUE" A1R,A1I,A2R,A2I;"REAL" A1R,A1I,A2R,A2I,C,SR,SI; "CODE" 34611; THE MEANING OF THE FORMAL PARAMETERS IS: A1R: ; ENTRY: THE REAL PART OF THE FIRST VECTORCOMPONENT; A1I: ; ENTRY: THE IMAGINARY PART OF THE FIRST VECTORCOMPONENT; A2R: ; ENTRY: THE REAL PART OF THE SECOND VECTORCOMPONENT; A2I: ; ENTRY: THE IMAGINARY PART OF THE SECOND VECTORCOMPONENT; C,SR,SI: ; EXIT: THE FACTORS THAT DETERMINE THE HOUSEHOLDER MATRIX. THE HOUSEHOLDERMATRIX, DEFINED BY: HA = B A = (A1,A2)' B = (-SIGN(A1R)*SQRT(A1*A1+A2*A2),O)', IS DETERMINED BY: ( -C SR+I*SI) (SR+I*SI C ) PROCEDURES USED: NONE; LANGUAGE: ALGOL 60; METHOD AND PERFORMANCE: AFTER A CALL OF CHSH2 YOU ARE ABLE TO ROTATE A COMPLEX VECTOR OF DIMENSION TWO BY MEANS OF THE FACTORS C,SR AND SI. EXAMPLE OF USE: CHSH2 IS USED IN QZI AND QZIVAL,SECTION 3.4 1SECTION : 1.2.7 (JANUARY 1976) PAGE 6 SOURCE TEXT(S) : 0"CODE" 34357; "PROCEDURE" ROTCOMCOL(L, U, I, J, AR, AI, CR, CI, S); "VALUE" L, U, I, J, CR, CI, S; "INTEGER" L, U, I, J; "REAL" CR, CI, S; "ARRAY" AR, AI; "BEGIN" "REAL" ARLI, AILI, ARLJ, AILJ; "FOR" L:= L "STEP" 1 "UNTIL" U "DO" "BEGIN" ARLI:= AR[L,I]; AILI:= AI[L,I]; ARLJ:= AR[L,J]; AILJ:= AI[L,J]; AR[L,I]:= CR * ARLI + CI * AILI - S * ARLJ; AI[L,I]:= CR * AILI - CI * ARLI - S * AILJ; AR[L,J]:= CR * ARLJ - CI * AILJ + S * ARLI; AI[L,J]:= CR * AILJ + CI * ARLJ + S * AILI; "END" "END" ROTCOMCOL; "EOP" 0"CODE" 34358; "PROCEDURE" ROTCOMROW(L, U, I, J, AR, AI, CR, CI, S); "VALUE" L, U, I, J, CR, CI, S; "INTEGER" L, U, I, J; "REAL" CR, CI, S; "ARRAY" AR, AI; "BEGIN" "REAL" ARIL, AIIL, ARJL, AIJL; "FOR" L:= L "STEP" 1 "UNTIL" U "DO" "BEGIN" ARIL:= AR[I,L]; AIIL:= AI[I,L]; ARJL:= AR[J,L]; AIJL:= AI[J,L]; AR[I,L]:= CR * ARIL + CI * AIIL + S * ARJL; AI[I,L]:= CR * AIIL - CI * ARIL + S * AIJL; AR[J,L]:= CR * ARJL - CI * AIJL - S * ARIL; AI[J,L]:= CR * AIJL + CI * ARJL - S * AIIL; "END" "END" ROTCOMROW; "EOP" 0"CODE" 34611; "PROCEDURE"CHSH2(A1R,A1I,A2R,A2I,C,SR,SI); "VALUE"A1R,A1I,A2R,A2I;"REAL"A1R,A1I,A2R,A2I,C,SR,SI; "BEGIN" "REAL" R; "IF" A2R^=0 "OR" A2I^=0 "THEN" "BEGIN" "IF" A1R^=0 "OR" A1I^=0 "THEN" "BEGIN" R:=SQRT(A1R*A1R+A1I*A1I);C:=R; SR:=(A1R*A2R+A1I*A2I)/R;SI:=(A1R*A2I-A1I*A2R)/R; R:=SQRT(C*C+SR*SR+SI*SI);C:=C/R;SR:=SR/R;SI:=SI/R "END" "ELSE" "BEGIN" SI:=C:=0;SR:=1 "END" "END" "ELSE" "BEGIN" C:=1;SR:=SI:=0 "END" "END" CHSH2; "EOP" 1SECTION : 1.5.3 (MARCH 1977) PAGE 1 AUTHORS: T.J.DEKKER & J.KOOPMAN. CONTRIBUTOR: J.KOOPMAN. INSTITUTE: UNIVERSITY OF AMSTERDAM. RECIEVED: 770328 BRIEF DESCRIPTION: LNGREATODECI CONVERTS A DOUBLE-LENGTH NUMBER TO A NUMBER IN DECIMAL FLOATING-POINT REPRESENTATION. THE RESULT CONSISTS OF A MANTISSA MANT OF MAGNITUDE<1 (AND >=.1) AND A DECIMAL EXPONENT EXPO. KEYWORDS: DOUBLE PRECISION ARITHMETIC, CONVERSION, DECIMAL REPRESENTATION. CALLING SEQUENCE: THE DECLARATION OF THE PROCEDURE IN THE CALLING PROGRAM READS: "PROCEDURE" LNGREATODECI(X, XX, S, MANT, EXPO); "VALUE"X,XX,S;"INTEGER"S,EXPO;"REAL"X,XX;"INTEGER""ARRAY"MANT; "CODE" 31100; THE MEANING OF THE FORMAL PARAMETERS IS: X,XX : ; ENTRY: THE HEAD (X) AND TAIL (XX) OF THE NUMBER THAT IS TO BE CONVERTED; S : ; ENTRY: THE DESIRED NUMBER OF SIGNIFICANT DIGITS OF THE CONVERTED VARIABLE. ONE SHOULD NOT CHOOSE S LARGER THAN THE NUMBER OF DIGITS CORRESPONDING TO THE DOUBLE LENGTH MACHINE PRECISION (FOR CDC: S<29 ). OTHERWISE, THE LAST DIGITS ARE USELESS, AS ALL OPERATIONS IN LNGREATODECI ARE PERFORMED IN DOUBLE-LENGTH A- RITHMETIC; IF S IS CHOSEN NONPOSITIVE,ONLY THE SIGN AND THE DECIMAL EXPONENT OF THE CONVERTED NUMBER ARE DELIVERED; 1SECTION : 1.5.3 (MARCH 1977) PAGE 2 MANT : ; "INTEGER""ARRAY"MANT[0:S]; EXIT: MANT[0]: THE SIGN OF THE DECIMAL NUMBER; MANT[I] (I^=0): THE I-TH SIGNIFICANT DIGIT OF THE MANTISSA OF THE CONVERTED NUMBER; I.E. THE VALUE OF THE MANTISSA EQUALS MANT[0]*(MANT[1]/10+MANT[2]/100+...MANT[S]/10**S); EXPO : ; EXIT: THE DECIMAL EXPONENT OF THE CONVERTED NUMBER, I.E. THE DOUBLE-LENGTH NUMBER (X,XX) APPROXIMATELY EQUALS MANTISSA*10**EXPO WITH THE VALUE OF MANTISSA GIVEN IN MANT. PROCEDURES USED: LNG SUB = CP31106. LNG MUL = CP31108. DP POW = CP31109. RUNNING TIME: ROUGHLY PROPORTIONAL TO LN(LN(X))+S. LANGUAGE: ALGOL 60. METHOD AND PERFORMANCE: LNGREATODECI DETERMINES THE DECIMAL EXPONENT EXPO. AFTER THAT, THE LONG REAL NUMBER (X,XX) IS DIVIDED BY 10**EXPO IN DOUBLE PRECISION. BY TRUNCATING THE RESULT, THE FIRST MOST SIGNIFICANT DIGIT OF THE MANTISSA IS OB- TAINED. SUBTRACTING THIS DIGIT FROM (X,XX)/10**EXPO, MULTIPLYING THE RESULT WITH 10, THE NEXT MOST SIGNI- FICANT DIGIT CAN BE OBTAINED BY TRUNCATION. THIS PRO- CESS OF SUBTRACTION, MULTIPLICATION AND TRUNCATION WILL BE REPEATED UNTIL S DIGITS ARE OBTAINED. FINALLY, THE MANTISSA THUS OBTAINED IS PROPERLY ROUNDED. 1SECTION : 1.5.3 (MARCH 1977) PAGE 3 EXAMPLE OF USE: "BEGIN""COMMENT"EXAMPLE OF USE OF LNGREATODECI AND DP POW; "INTEGER"S,EXPO; "REAL"OP,OPL; "PROCEDURE"DP POW(A,EXPON,C,CC);"CODE"31109; "PROCEDURE"LNGREATODECI(X,XX,S,MANT,EXPO);"CODE"31110; "PROCEDURE"PRINT(S,MANT,EXPONENT); "VALUE"S,EXPONENT;"INTEGER"S,EXPONENT;"INTEGER""ARRAY"MANT; "BEGIN""INTEGER"K; OUTCHARACTER(61,"("-++")",MANT[0]+2); "FOR"K:=1"STEP"1"UNTIL"S"DO" "BEGIN""IF"K=1"THEN"OUTPUT(61,"(""(".")"")"); OUTPUT(61,"("D")",MANT[K]) "END";OUTPUT(61,"(""(""")",+3D")",EXPONENT) "END"PRINT; DP POW(2,48,OP,OPL); "FOR"S:=0"STEP"4 UNTIL"28"DO" "BEGIN""INTEGER""ARRAY"MANT[0:S]; LNGREATODECI(OP,OPL,S,MANT,EXPO); PRINT(S,MANT,EXPO); OUTPUT(61,"("/")") "END" "END" DELIVERS: +"+015 +.2815"+015 +.28147498"+015 +.281474976711"+015 +.2814749767106560"+015 +.28147497671065600000"+015 +.281474976710656000000000"+015 +.2814749767106560000000000000"+015 1SECTION : 1.5.3 (MARCH 1977) PAGE 4 SOURCE TEXT(S): "CODE"31100; "PROCEDURE"LNGREATODECI(X,XX,S,MANT,EXPO); "VALUE"X,XX,S;"INTEGER"S,EXPO;"REAL"X,XX;"INTEGER""ARRAY"MANT; "BEGIN""INTEGER"I,K; "REAL"P,PP; "PROCEDURE"LNG SUB(A,AA,B,BB,C,CC);"CODE"31106; "PROCEDURE"LNG MUL(A,AA,B,BB,C,CC);"CODE"31107; "PROCEDURE"DP POW(A,EXPON,C,CC);"CODE"31109; MANT[0]:=SIGN(X);"IF"X<0"THEN""BEGIN"X:=-X;XX:=-XX"END"; "IF"X=0"THEN"EXPO:=0 "ELSE"EXPO:=ENTIER(LN(X)/LN(10))+1; DP POW(10,-EXPO,P,PP); LNG MUL(X,XX,P,PP,X,XX); "FOR"I:=0"WHILE"ENTIER(X)=0 & X^=0 "DO" "BEGIN"LNG MUL(X,XX,10,0,X,XX);EXPO:=EXPO-1"END"; "FOR"I:=1"STEP"1"UNTIL"S"DO" "BEGIN"K:=ENTIER(X);"IF"K>9"THEN"K:=9;MANT[I]:=K; LNG SUB(X,XX,K,0,P,PP);LNG MUL(P,PP,10,0,X,XX) "END"; "IF"ENTIER(X)>=5 "THEN""BEGIN""FOR"I:=S"STEP"-1"UNTIL"1"DO" "BEGIN"K:=MANT[I]+1; "IF"K<10"THEN""BEGIN"MANT[I]:=K;"GOTO"READY "END"; MANT[I]:=0 "END"; EXPO:=EXPO+1; "IF"S>0"THEN"MANT[1]:=1; READY: "END"; EXPO:=EXPO+1 "END" LNGREATODECI; "EOP" 1SECTION : 4.2.3.2 (NOVEMBER 1976) PAGE 1 AUTHORS: M.BAKKER. INSTITUTE: MATHEMATICAL CENTRE, AMSTERDAM. RECEIVED: 760131. BRIEF DESCRIPTION: THIS SECTION CONTAINS THE FOLLOWING PROCEDURES: (1) GSS JAC WGHTS: GIVEN THE TWO PARAMETERS ALFA AND BETA,THIS PROCEDURE CALCULATES THE N ZEROS OF THE N-TH JACOBI POLYNOMIAL AND THE CORRESPONDING GAUSS-CHRISTOFFEL NUMBERS NEEDED FOR THE N-POINT GAUSS-JACOBI QUADRATURE OVER [-1,+1] WITH WEIGHT FUNCTION W(X) = (1-X)**ALFA*(1+X)**BETA; (2) GSS LAG WGHTS: GIVEN THE PARAMETER ALFA,THIS PROCEDURE CALCULATES THE N ZEROS OF THE N-TH LAGUERRE POLYNOMIAL AND THE GAUSS-CHRISTOFFEL NUMBERS NEEDED FOR THE N-POINT GAUSS-LAGUERRE QUADRATURE OF F(X) OVER (0, INFINITY) WITH RESPECT TO THE WEIGHT FUNCTION W(X) = X**ALFA*EXP(-X). THESE PROCEDURES CAN BE USED FOR GAUSSIAN QUADRATURE- RULES OF THE JACOBI AND LAGUERRE TYPE.LET THE WEIGHT FUNCTION W(X) AND THE INTERVAL(A,B) DETERMINE THE SYSTEM OF POLYNOMIALS ORTHOGONAL ON (A,B) WITH RESPECT TO W(X).THEN THE N-POINT GAUSSIAN QUADRATURE RULE APPROXIMATES THE INTEGRAL TO X=B INTEGRAL F(X) W(X) DX FROM X=A BY THE EXPRESSION J=N SUM W[J] F(X[J]) J=1 WHERE THE ABSCISSAS X[J] ARE THE ZEROS OF THE N-TH POLYNOMIAL AND W[J] ARE THE CORESPONDING GAUSS-CHRISTOFFEL NUMBERS. 1SECTION : 4.2.3.2 (NOVEMBER 1976) PAGE 2 KEYWORDS: GAUSSIAN QUADRATURE, ZEROS OF ORTHOGONAL POLYNOMIALS, GAUSS-CHRISTOFFEL NUMBERS, GAUSSIAN WEIGHTS. LANGUAGE: ALGOL 60. REFERENCES: [1] M.ABRAMOWITZ AND I.A. STEGUN. HANDBOOK OF MATHEMATICAL FUNCTIONS, CH.22, [2] J.STOER, EINFUEHRUNG IN DIE NUMERISCHE MATHEMATIK 1, SPRINGER VERLAG, BERLIN, HEIDELBERG, GOETTINGEN. SUBSECTION: GSS JAC WGHTS. CALLING SEQUENCE: THE DECLARATION OF THE PROCEDURE IN THE CALLING PROGRAM READS: "PROCEDURE" GSS JAC WGHTS(N, ALFA, BETA, X, W); "VALUE" N, ALFA, BETA; "INTEGER" N; "REAL" ALFA, BETA; "CODE" 31425; THE MEANING OF THE FORMAL PARAMETERS IS: N: ; THE UPPER BOUND OF THE ARRAYS X AND W; N>=1; ALFA, BETA: ; THE PARAMETERS OF THE WEIGHT FUNCTION FOR THE JACOBI POLYNOMIALS; ALFA, BETA > -1; X; ; "ARRAY" X[1:N]; EXIT: X[I] IS THE I-TH ZERO OF THE N-TH JACOBI POLYNOMIAL; W: ; "ARRAY" W[1:N]; EXIT: W[I] IS THE GAUSS-CHRISTOFFEL NUMBER ASSOCIATED WITH THE I-TH ZERO OF THE N-TH JACOBI POLYNOMIAL. 1SECTION : 4.2.3.2 (NOVEMBER 1976) PAGE 3 PROCEDURES USED: GAMMA = CP 35061; ALL JAC ZER = CP 31370. REQUIRED CENTRAL MEMORY: TWO AUXILIARY ARRAYS OF N REALS ARE USED. RUNNING TIME: ROUGHLY PROPORTIONAL TO N CUBED. METHOD AND PERFORMANCE: AS IS WELL-KNOWN,THE GAUSSIAN QUADRATURE RULES ARE BASED ON THE ZEROS OF ORTHOGONAL POLYNOMIALS. PROCEDURES FOR THE COMPUTATION OF THESE ZEROS CAN BE FOUND IN SECTION 3.6.2. AFTER THE COMPUTATION OF THE ZEROS OF THE JACOBI POLYNOMIAL THE GAUSSIAN WEIGHTS ARE COMPUTED OF THE FORMULA J=N-1 W[I]=1/( SUM P(J, ALFA, BETA, X[I])**2) J=0 WHERE P(J, ALFA, BETA, X[I]) IS THE J-TH ORTHONORMAL JACOBI POLYNOMIAL;SEE FURTHER [2],CH.III. EXAMPLE OF USE: THE PROGRAM "BEGIN""COMMENT" EVALUATION OF THE INTEGRAL TO X=1 INTEGRAL (1+X)**2 * (1-X) * EXP(X) DX FROM X=-1 BY MEANS OF FIVE POINT GAUSS-JACOBI QUADRATURE. THE EXACT VALUE IS 2*EXP(1)-10/EXP(1); "REAL" ALFA, BETA, INT; "INTEGER" N; "ARRAY" X, W[1:5]; "REAL" "PROCEDURE" F(X); "VALUE" X; "REAL" X; F:=EXP(X); "PROCEDURE" GSS JAC WGHTS (N, ALFA, BETA, X, W); "CODE" 31425; ALFA:= 1; BETA:= 2; N:= 5; INT:= 0; GSS JAC WGHTS( N, ALFA, BETA, X, W); "FOR" N:= 1 "STEP" 1 "UNTIL" 5 "DO" INT:= INT + W[N] * F(X[N]); OUTPUT(61, "(" /, 4B+D. 4D"+ZD")", INT - 2 * EXP(1) + 10 / EXP(1)) "END" PRINTS THE FOLOWING RESULT: -1.5932"-10 1SECTION : 4.2.3.2 (NOVEMBER 1976) PAGE 4 SUBSECTION: GSS LAG WGHTS. CALLING SEQUENCE: THE DECLARATION OF THE PROCEDURE IN THE CALLING PROGRAM READS: "PROCEDURE" GSS LAG WGHTS (N, ALFA, X, W); "VALUE" N, ALFA; "INTEGER" N; "REAL" ALFA; "ARRAY" X, W; "CODE" 31427; THE MEANING OF THE FORMAL PARAMETERS IS: N: ; THE UPPER BOUND OF THE ARRAYS X AND W; N>=1; ALFA: ; THE PARAMETER OF THE WEIGHT FUNCTION FOR THE LAGUERRE POLYNOMIALS; ALFA>-1; X: ; "ARRAY" X[1: N]; EXIT: X[I] IS THE I-TH ZERO OF THE N-TH LAGUERRE POLYNOMIAL; W: ; "ARRAY" W[1: N]; EXIT: W[I] IS THE GAUSSIAN WEIGHT CORRESPONDING WITH THE I-TH ZERO OF THE N-TH LAGUERRE POLYNOMIAL. PROCEDURES USED: GAMMA = CP 35061, ALL LAG ZER = CP 31371. REQUIRED CENTRAL MEMORY: TWO AUXILIARY ARRAYS OF N REALS ARE USED. RUNNING TIME: ROUGHLY PROPORTIONAL TO N CUBED. METHOD AND PERFORMANCE: THE ZEROS AND WEIGHTS ARE COMPUTED IN THE SAME WAY AS IN THE PROCEDURE GSS JAC WGHTS. 1SECTION : 4.2.3.2 (NOVEMBER 1976) PAGE 5 EXAMPLE OF USE: THE PROGRAM "BEGIN" "COMMENT" COMPUTATION OF THE INTEGRAL FROM 0 TO INFINITY OF SIN(X)*EXP(-X) BY MEANS OF A TEN POINT GAUSS-LAGUERRE QUADRATURE.THE EXACT VALUE IS 0.5; "REAL" INT;"INTEGER" N;"ARRAY" X, W[1:10]; "REAL""PROCEDURE" F(X); "VALUE"X; "REAL"X; F:=SIN(X); "PROCEDURE" GSS LAG WGHTS(N, ALFA, X, W); "CODE" 31427; GSS LAG WGHTS(10, 0, X, W); INT:=0; "FOR" N:=10 "STEP" -1 "UNTIL" 1 "DO" INT:= INT + W[N] * F(X[N]); OUTPUT(61, "("-D.4D"-ZD")", INT-0.5) "END" PRINTS THE RESULT: 2.0497" -7 1SECTION : 4.2.3.2 (NOVEMBER 1976) PAGE 6 SOURCE TEXTS: "CODE" 31425; "PROCEDURE" GSS JAC WGHTS(N, ALFA, BETA, X, W); "VALUE" N, ALFA, BETA; "INTEGER" N; "REAL" ALFA, BETA; "ARRAY" X, W; "IF" ALFA = BETA "THEN" "BEGIN" "INTEGER" I, J, M; "ARRAY" B[1:N - 1]; "REAL" R0, R1, R2, S, H0, ALFA2, XI; "REAL" "PROCEDURE" GAMMA(X); "CODE" 35061; "PROCEDURE" ALL JAC ZER(N, ALFA, BETA, ZER); "CODE" 31370; ALL JAC ZER(N, ALFA, ALFA, X); ALFA2:= 2*ALFA; H0:= 2**(ALFA2 + 1)*GAMMA(1 + ALFA)**2/GAMMA(ALFA2 + 2); B[1]:= 1/SQRT(3 + ALFA2); M:= N - (N//2); "FOR" I:= 2 "STEP" 1 "UNTIL" N - 1 "DO" B[I]:= SQRT(I*(I + ALFA2)/(4*(I + ALFA)**2 - 1)); "FOR" I:= 1 "STEP" 1 "UNTIL" M "DO" "BEGIN" XI:= ABS(X[I]); R0:= 1; R1:= XI/B[1]; S:= 1 + R1*R1; "FOR" J:= 2 "STEP" 1 "UNTIL" N - 1 "DO" "BEGIN" R2:= (XI*R1 - B[J - 1]*R0)/B[J]; R0:= R1; R1:= R2; S:= S + R2*R2 "END"; W[I]:= W[N + 1 - I]:= H0/S "END" "END" "ELSE" "BEGIN" "INTEGER" I, J; "ARRAY" A, B[0:N]; "REAL" MIN, SUM, H0, R0, R1, R2, XI, ALFABETA; "PROCEDURE" ALL JAC ZER(N, ALFA, BETA, ZER); "CODE" 31370; "REAL" "PROCEDURE" GAMMA(X); "CODE" 35061; ALFABETA:= ALFA + BETA; MIN:= (BETA - ALFA)*ALFABETA; B[0]:= 0; SUM:= ALFABETA + 2; A[0]:= (BETA - ALFA)/SUM; A[1]:= MIN /SUM/(SUM + 2); B[1]:= 2*SQRT((1 + ALFA)*(1 + BETA)/(SUM + 1))/SUM; "FOR" I:= 2 "STEP" 1 "UNTIL" N - 1 "DO" "BEGIN" SUM:= I + I + ALFABETA; A[I]:= MIN/SUM/(SUM + 2); B[I]:= (2/SUM)* SQRT(I*(SUM - I)*(I + ALFA)*(I + BETA)/(SUM*SUM - 1)) "END"; H0:= 2**(ALFABETA + 1)*GAMMA(1 + ALFA)*GAMMA(1 + BETA)/ GAMMA(2 + ALFABETA); ALL JAC ZER(N, ALFA, BETA, X); "FOR" I:= 1 "STEP" 1 "UNTIL" N "DO" "BEGIN" XI:= X[I]; R0:= 1; R1:= (XI - A[0])/B[1]; SUM:= 1 + R1*R1; "FOR" J:= 2 "STEP" 1 "UNTIL" N - 1 "DO" "BEGIN" R2:= ((XI - A[J - 1])*R1 - B[J - 1]*R0)/B[J]; SUM:= SUM + R2*R2; R0:= R1; R1:= R2 "END"; W[I]:= H0/SUM "END" "END" GSS JAC WGHTS; "EOP" 1SECTION : 4.2.3.2 (NOVEMBER 1976) PAGE 7 "CODE" 31427; "PROCEDURE" GSS LAG WGHTS(N, ALFA, X, W); "VALUE" N, ALFA; "INTEGER" N; "REAL" ALFA; "ARRAY" X, W; "BEGIN" "INTEGER" I, J; "REAL" H0, S, R0, R1, R2, XI; "ARRAY" A, B[0:N]; "PROCEDURE" ALL LAG ZER(N, ALFA, X); "CODE" 31371; "REAL" "PROCEDURE" GAMMA(X); "CODE" 35061; A[0]:= 1 + ALFA; A[1]:= 3 + ALFA; B[1]:= SQRT(A[0]); "FOR" I:= 2 "STEP" 1 "UNTIL" N - 1 "DO" "BEGIN" A[I]:= I + I + ALFA + 1; B[I]:= SQRT(I*(I + ALFA)) "END"; ALL LAG ZER(N, ALFA, X); H0:= GAMMA(1 + ALFA); "FOR" I:= 1 "STEP" 1 "UNTIL" N "DO" "BEGIN" XI:= X[I]; R0:= 1; R1:= (XI - A[0])/B[1]; S:= 1 + R1*R1; "FOR" J:= 2 "STEP" 1 "UNTIL" N - 1 "DO" "BEGIN" R2:= ((XI - A[J - 1])*R1 - B[J - 1]*R0)/B[J]; R0:= R1; R1:= R2; S:= S + R2*R2 "END"; W[I]:= H0/S "END" "END" GSS LAG WGHTS; "EOP" 1SECTION : 2.2.2.1 (NOVEMBER 1978) PAGE 1 AUTHOR: C.G. VAN DER LAAN CONTRIBUTORS: C.G. VAN DER LAAN, M. VOORINTHOLT INSTITUTE: REKENCENTRUM RIJKSUNIVERSITEIT GRONINGEN RECEIVED: 780601 BRIEF DESCRIPTION: THIS SECTION CONTAINS SIX PROCEDURES FOR THE EVALUATION OF ORTHOGONAL POLYNOMIALS: ORTPOL,ORTPOLSYM: EVALUATE AN ORTHOGONAL POLYNOMIAL, ALLORTPOL,ALLORTPOLSYM: EVALUATE ALL ORTHOGONAL POLYNOMIALS OF DEGREE LESS THAN A GIVEN POSITIVE INTEGER. SUMORTPOL,SUMORTPOLSYM: EVALUATE A SERIES OF ORTHOGONAL POLYNOMIALS. THE PROCEDURES ENDING WITH SYM ARE EFFICIENT VERSIONS FOR SYMMETRICAL POLYNOMIALS (I.E. ODD INDEXED POLYNOMIALS THAT ARE ODD FUNCTIONS AND EVEN INDEXED POLYNOMIALS THAT ARE EVEN FUNCTIONS); KEYWORDS: ORTHOGONAL POLYNOMIAL, SERIES OF ORTHOGONAL POLYNOMIALS, LINEAR THREE TERM (IN)HOMOGENEOUS RECURRENCE RELATION. DATA AND RESULTS: ORTHOGONAL POLYNOMIALS CAN BE CHARACTERIZED BY A RECURRENCE RELATION OF THE FOLLOWING FORM A1[K] * F[K+1](X) = (A2[K] + X * A3[K]) * F[K](X) - A4[K] * F[K-1](X), WHERE AI[K] ARE REAL NUMBERS. SEE FOR INSTANCE TABLE 22.7 IN ABRAMOWITZ AND STEGUN (1964) FOR THE CLASSICAL ORTHOGONAL POLYNOMIALS. BY AN ELEMENTARY TRANSFORMATION, THE COEFFICIENTS IN THE RECURRENCE RELATION ABOVE CAN BE MADE SUCH THAT THE RECURRENCE RELATION IS GIVEN BY P[K+1](X) = (X - B[K]) * P[K](X) - C[K] * P[K-1](X), P[0](X) = 1, P[1](X) = X - B[0]. 1SECTION : 2.2.2.1 (NOVEMBER 1978) PAGE 2 IN THIS WAY WE OBTAIN A NORMALIZATION OF THE ORTHOGONAL POLYNOMIAL SUCH THAT THE LEADING COEFFICIENT IN THE EXPLICIT REPRESENTATION OF P[K](X) EQUALS 1. AS A CONSEQUENCE THE FOLLOWING RELATION HOLDS BETWEEN THE VALUES OBTAINED BY OUR PROCEDURES (E.G. ORTPOL) AND THE VALUES FROM THE REPRESENTATION IN ABRAMOWITZ AND STEGUN (1964) (I.C. F) N-1 ORTPOL[N](X) = PROD (A1[K]/A3[K]) * F[N](X) , K=0 WHERE A1[K], A3[K], F[N] ARE DETERMINED BY 22.4 AND 22.7 IN ABRAMOWITZ AND STEGUN(1964). WE NOTICE THAT OVERFLOW/UNDERFLOW MAY OCCUR EARLIER AS A CONSEQUENCE OF OUR NORMALIZATION. IN ORDER TO AVOID MISTAKES WHEN OBTAINING THE RECURRENCE COEFFICIENTS THE FOLLOWING TABLE GIVES THE RECURRENCE COEFFICIENTS FOR THE CLASSICAL ORTHOGONAL POLYNOMIALS (NOTE THAT THE FIRST AND SECOND POLYNOMIAL ARE DEFINED BY THE NORMALIZATION AND B[0]): POLYNOMIAL KIND : RECURRENCE COEFFICIENTS ----------------------:-------------------------------------------- : B[K] : C[K] :-----------------------:-------------------- CHEBYSHEV (1-ST KIND) : 0 : 1/2 , K=1 : : 1/4 , K>1 : : CHEBYSHEV (2-ND KIND) : 0 : 1/4 : : LEGENDRE : 0 : K**2/(4*K**2-1) : : JACOBI : -(ALPHA**2-BETA**2)/ : 4*(1+ALPHA)* : ((2*K+ALPHA+BETA)* : (1+BETA)/((ALPHA+ : (2*K+ALPHA+BETA+2)) : BETA+2)**2*(ALPHA+ : : BETA+3)) , K=1 : : : : 4*K*(K+ALPHA)* : : (K+BETA)*(K+ALPHA+ : : BETA)/((2*K+ALPHA+ : : BETA)**2*((2*K+ : : ALPHA+BETA)**2-1)) : : , K>1 : : : : (ALPHA,BETA > -1) : : LAGUERRE : 2*K+ALPHA+1 : K*(K+ALPHA) : : HERMITE : 0 : K/2 1SECTION : 2.2.2.1 (NOVEMBER 1978) PAGE 3 IN GENERAL THE RECURRENCE COEFFICIENTS MAY BE OBTAINED BY USE OF THE PROCEDURE RECCOF. THE ABOVE TABLE IS OBTAINED BY ADAPTION OF TABLE 22.7,ABRAMOWITZ AND STEGUN (1964) P. 782, AS FOLLOWS: (NOTE THAT N>=1; FOR N=0 CONSULT 22.4) B[I]:=-A2[I]/A3[I] , I=0,1,...,N-1 C[I]:=(A4[I]*A1[I-1])/(A3[I]*A3[I-1]) , I=1,2,...,N-1. METHOD AND PERFORMANCE: LET THE ORTHOGONAL POLYNOMIAL BE DEFINED BY P[K+1](X)=(X-B[K])*P[K](X)-C[K]*P[K-1](X) , K=1,2,...N-1 WHERE B[0:N-1], C[1:N-1] CONTAIN THE RECURRENCE COEFFICIENTS AND P[1](X)=X-B[0], P[0](X)=1 (SEE STOER 1972, P. 119). THEN N-1 / X-B[K] 1 \ / 1 \ P[N](X) = (X-B[0],1) * PROD [ ] [ ] , N=1,2,..... K=1 \ -C[K] 0 / \ 0 / AND A[0]+A[1]*P[1](X)+...+A[N]*P[N](X)= N J-1 / X-B[K] 1 \ / A[J] \ A[0] + (X-B[0],1) * SUM PROD [ ] [ ] . J=1 K=1 \ -C[K] 0 / \ 0 / THESE EXPRESSIONS ARE EVALUATED BY A GENERALISED HORNER RULE. (SEE ALSO LUKE, 1969, P. 327). REFERENCES: ABRAMOWITZ, M. & I. STEGUN (1964): HANDBOOK OF MATHEMATICAL FUNCTIONS, NATIONAL BUREAU OF STANDARDS, WASHINGTON D.C. LUKE, Y.L. (1969): THE SPECIAL FUNCTIONS AND THEIR APPROXIMATIONS I, ACADEMIC PRESS, LONDON, NEW YORK. STOER, J. (1972): EINFUEHRUNG IN DIE NUMERISCHE MATHEMATIK I, SPRINGER VERLAG, BERLIN. 1SECTION : 2.2.2.1 (NOVEMBER 1978) PAGE 4 SUBSECTION: ORTPOL. CALLING SEQUENCE: THE DECLARATION OF THE PROCEDURE IN THE CALLING PROGRAM READS: "REAL" "PROCEDURE" ORTPOL(N,X,B,C); "VALUE" N,X; "INTEGER" N; "REAL" X; "ARRAY" B,C; "CODE" 31044; ORTPOL DELIVERS THE VALUE OF THE ORTHOGONAL POLYNOMIAL OF DEGREE N FOR THE ARGUMENT X AS DETERMINED BY THE RECURRENCE COEFFICIENTS. THE MEANING OF THE FORMAL PARAMETERS IS: N: ; ENTRY: THE DEGREE OF THE POLYNOMIAL; X: ; ENTRY: THE ARGUMENT OF THE ORTHOGONAL POLYNOMIAL; B,C: ; "ARRAY" B[0:N-1], C[1:N-1]; ENTRY: THE RECURRENCE COEFFICIENTS (SEE DATA AND RESULTS). SUBSECTION: ORTPOLSYM. CALLING SEQUENCE: THE DECLARATION OF THE PROCEDURE IN THE CALLING PROGRAM READS: "REAL" "PROCEDURE" ORTPOLSYM(N,X,C); "VALUE" N,X; "INTEGER" N; "REAL" X; "ARRAY" C; "CODE" 31048; ORTPOLSYM DELIVERS THE VALUE OF THE ORTHOGONAL POLYNOMIAL OF DEGREE N FOR THE ARGUMENT X AS DETERMINED BY THE RECURRENCE COEFFICIENTS. THE MEANING OF THE FORMAL PARAMETERS IS: N: ; ENTRY: THE DEGREE OF THE POLYNOMIAL; X: ; ENTRY: THE ARGUMENT OF THE ORTHOGONAL POLYNOMIAL; C: ; "ARRAY" C[1:N-1]; ENTRY: THE RECURRENCE COEFFICIENTS (SEE DATA AND RESULTS). 1SECTION : 2.2.2.1 (NOVEMBER 1978) PAGE 5 SUBSECTION: ALLORTPOL. CALLING SEQUENCE: THE DECLARATION OF THE PROCEDURE IN THE CALLING PROGRAM READS: "PROCEDURE" ALLORTPOL(N,X,B,C,P); "VALUE" N,X; "INTEGER" N; "REAL" X; "ARRAY" B,C,P; "CODE" 31045; THE MEANING OF THE FORMAL PARAMETERS IS: N,X,B,C: SEE ORTPOL (THIS SECTION); P: ; "ARRAY" P[0:N]; EXIT: P[K] CONTAINS, FOR THE ARGUMENT, THE VALUE OF THE K-TH ORTHOGONAL POLYNOMIAL AS DEFINED BY THE RECURRENCE COEFFICIENTS. SUBSECTION: ALLORTPOLSYM. CALLING SEQUENCE: THE DECLARATION OF THE PROCEDURE IN THE CALLING PROGRAM READS: "PROCEDURE" ALLORTPOLSYM(N,X,C,P); "VALUE" N,X; "INTEGER" N; "REAL" X; "ARRAY" C,P; "CODE" 31049; THE MEANING OF THE FORMAL PARAMETERS IS: N,X,C: SEE ORTPOLSYM (THIS SECTION); P: ; "ARRAY" P[0:N]; EXIT: P[K] CONTAINS, FOR THE ARGUMENT, THE VALUE OF THE K-TH ORTHOGONAL POLYNOMIAL AS DEFINED BY THE RECURRENCE COEFFICIENTS. SUBSECTION: SUMORTPOL. CALLING SEQUENCE: THE DECLARATION OF THE PROCEDURE IN THE CALLING PROGRAM READS: "REAL" "PROCEDURE" SUMORTPOL(N,X,B,C,A); "VALUE" N,X; "INTEGER" N; "REAL" X; "ARRAY" B,C,A; "CODE" 31047; 1SECTION : 2.2.2.1 (NOVEMBER 1978) PAGE 6 SUMORTPOL : DELIVERS THE VALUE OF THE POLYNOMIAL A[0]+A[1]*P[1](X)+...+A[N]*P[N](X) WHERE P[K](X) IS THE K-TH ORTHOGONAL POLYNOMIAL. THE MEANING OF THE FORMAL PARAMETERS IS: N,X,B,C: SEE ORTPOL (THIS SECTION); A: ; "ARRAY" A[0:N]; ENTRY: THE COEFFICIENTS OF THE SERIES EXPANSION A[0]+A[1]*P[1](X)+...+A[N]*P[N](X) WHERE P[K](X) IS THE K-TH ORTHOGONAL POLYNOMIAL AS DEFINED BY THE RECURRENCE COEFFICIENTS. SUBSECTION: SUMORTPOLSYM. CALLING SEQUENCE: THE DECLARATION OF THE PROCEDURE IN THE CALLING PROGRAM READS: "REAL" "PROCEDURE" SUMORTPOLSYM(N,X,C,A); "VALUE" N,X; "INTEGER" N; "REAL" X; "ARRAY" C,A; "CODE" 31058; SUMORTPOLSYM : DELIVERS THE VALUE OF THE POLYNOMIAL A[0]+A[1]*P[1](X)+...+A[N]*P[N](X) WHERE P[K](X) IS THE K-TH ORTHOGONAL POLYNOMIAL. THE MEANING OF THE FORMAL PARAMETERS IS: N,X,C: SEE ORTPOLSYM (THIS SECTION); A: ; "ARRAY" A[0:N]; ENTRY: THE COEFFICIENTS OF THE SERIES EXPANSION A[0]+A[1]*P[1](X)+...+A[N]*P[N](X) WHERE P[K](X) IS THE K-TH ORTHOGONAL POLYNOMIAL AS DEFINED BY THE RECURRENCE COEFFICIENTS. 1SECTION : 2.2.2.1 (NOVEMBER 1978) PAGE 7 EXAMPLE OF USE: THE FOLLOWING PROGRAM DELIVERS THE VALUES OF THE LAGUERRE POLYNOMIAL OF DEGREES 0,1,2,3,4,5 FOR X=0 BY MEANS OF THE PROCEDURE ALLORTPOL (B[K]=2*K+1, C[K]=K**2): "BEGIN" "ARRAY" B[0:4], C[1:4], P[0:5]; "INTEGER" I; "PROCEDURE" ALLORTPOL(N,X,B,C,P); "CODE" 31045; B[0]:=1; "FOR" I:=1 "STEP" 1 "UNTIL" 4 "DO" "BEGIN" B[I]:=2*I+1; C[I]:=I**2; "END" I; ALLORTPOL(5,0,B,C,P); OUTPUT(61,"("2/,6(2B,+2ZD.D)")",(P[I],I:=0:5)); "END" PROGRAM; RESULTS (NOTE THE DIFFERENCE WITH FIGURE 22.9 IN ABRAMOWITZ AND STEGUN (1964) BECAUSE OF THE NORMALIZATION): +1.0 -1.0 +2.0 -6.0 +24.0 -120.0 1SECTION : 2.2.2.1 (NOVEMBER 1978) PAGE 8 SOURCE TEXTS: "CODE" 31044; "REAL" "PROCEDURE" ORTPOL(N,X,B,C); "VALUE" N,X; "INTEGER" N; "REAL" X; "ARRAY" B,C; "IF" N=0 "THEN" ORTPOL:=1 "ELSE" "BEGIN" "INTEGER" K,L; "REAL" R,S,H; R:=X-B[0]; S:=1; L:=N-1; "FOR" K:=1 "STEP" 1 "UNTIL" L "DO" "BEGIN" H:=R; R:=(X-B[K])*R-C[K]*S; S:=H; "END"; ORTPOL:=R; "END" ORTPOL; "EOP" 0 "CODE" 31048; "REAL" "PROCEDURE" ORTPOLSYM(N,X,C); "VALUE" N,X; "INTEGER" N; "REAL" X; "ARRAY" C; "IF" N=0 "THEN" ORTPOLSYM:=1.0 "ELSE" "BEGIN" "INTEGER" K,L; "REAL" R,S,H; R:=X; S:=1.0; L:=N-1; "FOR" K:=1 "STEP" 1 "UNTIL" L "DO" "BEGIN" H:=R; R:=X*R-C[K]*S; S:=H "END"; ORTPOLSYM:=R "END" ORTPOLSYM; "EOP" 0 "CODE" 31045; "PROCEDURE"ALLORTPOL(N,X,B,C)RESULTS:(P); "VALUE" N,X; "INTEGER" N;"REAL" X;"ARRAY" B,C,P; "IF" N=0 "THEN" P[0]:=1"ELSE" "BEGIN" "INTEGER" K,K1; "REAL" R,S,H; R:=P[1]:=X-B[0]; S:=P[0]:=1;K:=1; "FOR" K1:=2 "STEP" 1 "UNTIL" N "DO" "BEGIN" H:=R; P[K1]:=R:=(X-B[K])*R-C[K]*S; S:=H; K:=K1; "END"; "END" ALLORTPOL "EOP" 1SECTION : 2.2.2.1 (NOVEMBER 1978) PAGE 9 "CODE" 31049; "PROCEDURE" ALLORTPOLSYM(N,X,C)RESULTS:(P); "VALUE" N,X; "INTEGER" N; "REAL" X; "ARRAY" C,P; "IF" N=0 "THEN" P[0]:=1.0 "ELSE" "BEGIN" "INTEGER" K; "REAL" R,S,H; R:=P[1]:=X; S:=P[0]:=1.0; "FOR" K:=2 "STEP" 1 "UNTIL" N "DO" "BEGIN" H:=R; P[K]:=R:=X*R-C[K-1]*S; S:=H "END"; "END" ALLORTPOLSYM; "EOP" 0 "CODE" 31047; "REAL" "PROCEDURE" SUMORTPOL(N,X,B,C,A); "VALUE" N,X; "INTEGER" N; "REAL" X; "ARRAY" B,C,A; "IF" N=0 "THEN" SUMORTPOL:=A[0] "ELSE" "BEGIN" "INTEGER" K; "REAL" H,R,S; R:=A[N]; S:=0; "FOR" K:=N-1 "STEP" -1 "UNTIL" 1 "DO" "BEGIN" H:=R; R:=A[K]+(X-B[K])*R+S; S:= -C[K]*H "END"; SUMORTPOL:=A[0]+(X-B[0])*R+S "END" SUMORTPOL; "EOP" 0 "CODE" 31058; "REAL" "PROCEDURE" SUMORTPOLSYM(N,X,C,A); "VALUE" N,X; "INTEGER" N; "REAL" X; "ARRAY" C,A; "IF" N=0 "THEN" SUMORTPOLSYM:=A[0] "ELSE" "BEGIN" "INTEGER" K; "REAL" H,R,S; R:=A[N]; S:=0; "FOR" K:=N-1 "STEP" -1 "UNTIL" 1 "DO" "BEGIN" H:=R; R:=A[K]+X*R+S; S:= -C[K]*H "END"; SUMORTPOLSYM:=A[0]+X*R+S "END" SUMORTPOLSYM; "EOP" 1SECTION : 5.2.1.1.3 (NOVEMBER 1976) PAGE 1 AUTHORS: P.A. BEENTJES, H.G.J. ROZENHART. INSTITUTE: MATHEMATICAL CENTRE. RECEIVED: 760201 BRIEF DESCRIPTION: ARKMAT SOLVES AN INITIAL VALUE PROBLEM, GIVEN AS A SYSTEM OF FIRST ORDER (NON-LINEAR) DIFFERENTIAL EQUATIONS BY MEANS OF A STABILIZED RUNGE KUTTA METHOD; IN PARTICULAR THIS PROCEDURE IS SUITABLE FOR THE INTEGRATION OF SYSTEMS WHERE THE DEPENDENT VARIABLE AND THE RIGHTHAND SIDE ARE STORED IN A RECTANGULAR ARRAY INSTEAD OF A VECTOR , I.E. DU / DT = F( T, U), WHERE U AND F ARE (N * M) MATRICES ( SEE METHOD AND PERFORMANCE). KEYWORDS: MATRIX DIFFERENTIAL EQUATIONS, INITIAL VALUE PROBLEMS, EXPLICIT ONE-STEP METHODS, STABILIZED RUNGE KUTTA METHODS. CALLING SEQUENCE: THE DECLARATION OF THE PROCEDURE IN THE CALLING PROGRAM READS: "PROCEDURE" ARKMAT(T, TE, M, N, U, DER, TYPE, ORDER, SPR, OUT); "VALUE" M, N, TYPE, ORDER; "INTEGER" M, N, TYPE, ORDER; "REAL" T, TE, SPR; "ARRAY" U; "PROCEDURE" DER, OUT; "CODE" 33066; 1SECTION : 5.2.1.1.3 (FEBRUARY 1979) PAGE 2 THE MEANING OF THE FORMAL PARAMETERS IS T: ; THE INDEPENDENT VARIABLE T; ENTRY: THE INITIAL VALUE T0; EXIT : THE FINAL VALUE TE; TE: ; ENTRY: THE FINAL VALUE OF T; M: ; NUMBER OF COLUMNS OF U; N: ; NUMBER OF ROWS OF U; U: ; "ARRAY" U[1:N,1:M]; ENTRY: THE INITIAL VALUES OF THE SOLUTION OF THE SYSTEM OF DIFFERENTIAL EQUATIONS AT T=T0; EXIT : THE VALUES OF THE SOLUTION AT T=TE; DER: ; THE HEADING OF THIS PROCEDURE READS: "PROCEDURE" DER(T, V, FTV); "VALUE" T; "REAL" T; "ARRAY" V, FTV; THIS PROCEDURE MUST BE GIVEN BY THE USER AND PERFORMS AN EVALUATION OF THE RIGHTHAND SIDE F( T, V) OF THE SYSTEM; UPON COMPLETION OF DER,THE RIGHTHAND SIDE SHOULD BE STORED IN FTV[1:N,1:M]; TYPE: ; ENTRY: THE TYPE OF THE SYSTEM OF DIFFERENTIAL EQUATIONS TO BE SOLVED; THE USER SHOULD SUPPLY ONE OF THE FOLLOWING VALUES; 1: IF NO SPECIFICATION OF THE TYPE CAN BE MADE; 2: IF THE EIGENVALUES OF THE JACOBIAN MATRIX OF THE RIGHTHAND SIDE ARE NEGATIVE REAL; 3: IF THE EIGENVALUES OF THE JACOBIAN MATRIX OF THE RIGHTHAND SIDE ARE PURELY IMAGINARY; ORDER: ; THE ORDER OF THE RUNGE KUTTA METHOD USED; ENTRY: FOR TYPE=2 THE USER MAY CHOOSE ORDER=1 OR ORDER=2; ORDER SHOULD BE 2 FOR THE OTHER TYPES; EXIT : IF ORDER IS SET TO ANOTHER VALUE,IT IS ASSUMED TO BE (IF TYPE=2 "THEN" 1 "ELSE" 2); SPR: ; ENTRY: THE SPECTRAL RADIUS OF THE JACOBIAN MATRIX OF THE RIGHTHAND SIDE, WHEN THE SYSTEM IS WRITTEN IN ONE DIMENSIONAL FORM (I.E. VECTORFORM); THE INTEGRATION STEP WILL EQUAL CONSTANT/SPR (SEE DATA AND RESULTS); IF NECESSARY SPR CAN BE UPDATED (AFTER EACH STEP) BY MEANS OF THE PROCEDURE OUT; OUT: THE HEADING OF THIS PROCEDURE READS: "PROCEDURE" OUT; AFTER EACH INTEGRATION STEP PERFORMED, INFORMATION CAN BE OBTAINED OR UPDATED BY THIS PROCEDURE, E.G. THE VALUES OF T, U[1:N,1:M] AND SPR. 1SECTION : 5.2.1.1.3 (FEBRUARY 1979) PAGE 3 DATA AND RESULTS: IF THE USER WANTS TO PERFORM THE INTEGRATION WITH A PRESCRIBED STEP H, HE HAS TO GIVE SPR THE VALUE CONSTANT/H, WHERE CONSTANT HAS THE FOLLOWING VALUES: CONSTANT= 4.3 IF TYPE=1 AND ORDER=2; CONSTANT= 156 IF TYPE=2 AND ORDER=1; CONSTANT= 64 IF TYPE=2 AND ORDER=2; CONSTANT= 8 IF TYPE=3 AND ORDER=2; PROCEDURES USED: ELMCOL = CP34023, DUPMAT = CP31035. REQUIRED CENTRAL MEMORY: TWO AUXILIARY ARRAYS OF ORDER N*M ARE DECLARED. METHOD AND PERFORMANCE: ARKMAT IS AN IMPLEMENTATION OF LOW ORDER STABILIZED RUNGE KUTTA METHODS (SEE REFERENCE[1]); THE INTEGRATION STEPSIZE USED WILL DEPEND ON: 1. THE TYPE OF SYSTEM TO BE SOLVED (I.E. HYPERBOLIC OR PARABOLIC); 2. THE SPECTRAL RADIUS OF THE JACOBIAN MATRIX OF THE SYSTEM; 3. THE INDICATED ORDER OF THE PARTICULAR RUNGE KUTTA METHOD; THE PROCEDURE ARKMAT IS ESPECIALLY INTENDED FOR SYSTEMS OF DIFFERENTIAL EQUATIONS ARISING FROM INITIAL BOUNDARY VALUE PROBLEMS IN TWO DIMENSIONS, E.G. WHEN THE METHOD OF LINES IS APPLIED TO THIS KIND OF PROBLEMS,THE RIGHTHAND SIDE OF THE RESULTING SYSTEM IS MUCH EASIER TO DESCRIBE IN MATRIX THAN IN VECTOR FORM; BECAUSE OF THIS FACT THE ARRAY OF DEPENDENT VARIABLES U IS A MATRIX, RATHER THAN A VECTOR. REFERENCE: [1]. P.J. VAN DER HOUWEN. STABILIZED RUNGE KUTTA METHOD WITH LIMITED STORAGE REQUIREMENTS. MATH. CENTR. REPORT TW 124/71. 1SECTION : 5.2.1.1.3 (NOVEMBER 1976) PAGE 4 EXAMPLE OF USE: GIVEN THE FOLLOWING SYSTEM OF EQUATIONS: DU / DT = V( T, X, Y), (1) DV / DT = D( DU / DX) / DX + D( DU / DY) / DY, ( ORIGINATING FROM THE INITIAL BOUNDARY VALUE PROBLEM D( DU / DT) / DT = D( DU / DX) / DX + D( DU / DY) / DY, ON THE DOMAIN 0 <= X <= PI , 0 <= Y <= 1 ), WITH THE FOLLOWING BOUNDARY CONDITIONS: U( T, 0, Y) = U( T, PI, Y) = U( T, X, 1) = 0, U( T, X, 0) = SIN( X ) * COS( SQRT( 1 + PI * PI / 4) * T), AND THE INITIAL VALUES: U( 0, X, Y) = SIN( X ) * COS( PI * Y / 2), V( 0, X, Y) = 0; BY APPLYING THE METHOD OF LINES TO PROBLEM (1), USING A TEN BY TEN GRID ON THE INDICATED DOMAIN, THE SYSTEM IS TRANSFORMED TO A MATRIX -DIFFERENTIAL EQUATION; THE SOLUTION OF THE LATTER PROBLEM AT T=1 IS COMPUTED BY THE FOLLOWING PROGRAM, USING A CONSTANT STEPSIZE .1; "BEGIN" "REAL" HPI,H1,H2,H1K,H2K,T,TE; "INTEGER" I,J,N,M,TYP,ORDE,TEL;"ARRAY" U[1:20,1:10]; "PROCEDURE" ARKMAT(T,TE,M,N,U,DER,TYPE,ORDER,SPR,OUT); "CODE" 33066; "PROCEDURE" INIMAT(LR,UR,LC,UC,A,X); "CODE" 31011; "PROCEDURE" DERIV(T,U,DU);"REAL" T;"ARRAY" U,DU; "BEGIN" "FOR" I:=2 "STEP" 1 "UNTIL" N-1 "DO" "FOR" J:=2 "STEP" 1 "UNTIL" M-1 "DO" "BEGIN" DU[I,J]:=U[I+N,J]; DU[I+N,J]:=(U[I,J+1]-2*U[I,J]+U[I,J-1])/H1K+ (U[I+1,J]-2*U[I,J]+U[I-1,J])/H2K "END"; "FOR" J:=1,M "DO" "BEGIN" INIMAT(N+1,N+N,J,J,DU,0); "FOR" I:=1 "STEP" 1 "UNTIL" N "DO" DU[I,J]:=U[N+1,J] "END"; "FOR" I:=1,N "DO" "FOR" J:=2 "STEP" 1 "UNTIL" M-1 "DO" "BEGIN" DU[I,J]:=U[I+N,J]; "IF" I=1 "THEN" DU[N+1,J]:=(U[1,J+1]-2*U[1,J]+U[1,J-1])/H1K+ (2*U[2,J]-2*U[1,J])/H2K "ELSE" DU[2*N,J]:=0 "END" "END" DERIV; 1SECTION : 5.2.1.1.3 (NOVEMBER 1976) PAGE 5 "PROCEDURE" OUT; "BEGIN" TEL:=TEL+1; "IF" T=TE "THEN" "BEGIN" OUTPUT(61,"("//,3B,"("X")",7B,"("Y")",4B, "("U(1,X,Y)")",7B,"("U(1,X,Y)")",/,16B,"("COMPUTED")",7B, "("EXACT")",//")"); OUTPUT(61,"("10(2(-D.3D2B),2(-D.6D6B),/)")", ((I-1)*H1,(I-1)*H2,U[I,I],SIN(H1*(I-1))*COS(HPI*H2*(I-1))* COS(T*SQRT(1+HPI*HPI)),I:=1:10)); OUTPUT(61,"("/,"("NUMBER OF INTEGRATION STEPS: ")" ,ZZZD")",TEL); OUTPUT(61,"("//,"(" TYPE IS:")",ZD,"(" ORDER IS:")", ,ZD")",TYP,ORDE); "END"; "END" OUT; "PROCEDURE" START; "BEGIN" "FOR" J:=1 "STEP" 1 "UNTIL" M "DO" U[N,J]:=SIN(H1*(J-1)); "FOR" I:=1 "STEP" 1 "UNTIL" N "DO" "BEGIN" "REAL" COS1; COS1:=COS(H2*HPI*(I-1)); "FOR" J:=1 "STEP" 1 "UNTIL" M "DO" U[I,J]:=U[N,J]*COS1 "END" INIMAT(N+1,N+N,1,M,U,0) "END" START; HPI:=2*ARCTAN(1);H2:=1/9;H1:=(2*HPI)/9;N:=M:=10; H1K:=H1*H1;H2K:=H2*H2;TEL:=0; T:=0; TE:=1 ; START; TYP:=3; ORDE:=2; ARKMAT(T,TE,M,N+N,U,DERIV,TYP,ORDE,80,OUT) "END" THIS PROGRAM DELIVERS: X Y U(1,X,Y) U(1,X,Y) COMPUTED EXACT 0.000 0.000 0.000000 0.000000 0.349 0.111 -0.095201 -0.096735 0.698 0.222 -0.170723 -0.173474 1.047 0.333 -0.211983 -0.215398 1.396 0.444 -0.213228 -0.216663 1.745 0.556 -0.178920 -0.181802 2.094 0.667 -0.122388 -0.124360 2.443 0.778 -0.062138 -0.063139 2.793 0.889 -0.016787 -0.017057 3.142 1.000 0.000000 -0.000000 NUMBER OF INTEGRATION STEPS: 10 TYPE IS: 3 ORDER IS: 2 1SECTION : 5.2.1.1.3 (NOVEMBER 1976) PAGE 6 SOURCE TEXT(S): "CODE" 33066; "PROCEDURE" ARKMAT( T, TE, M, N, U, DER, TYPE, ORDER, SPR, OUT); "VALUE" M,N,TYPE,ORDER; "INTEGER" M,N,TYPE,ORDER; "REAL" T,TE,SPR; "ARRAY" U; "PROCEDURE" DER,OUT; "BEGIN" "INTEGER" SIG,L; "REAL" TAU; "ARRAY" LAMBDA[1:9],UH,DU[1:N,1:M]; "BOOLEAN" LAST; "PROCEDURE" ELMCOL(L,U,I,J,A,B,X); "CODE" 34023; "PROCEDURE" DUPMAT(L,U,I,J,A,B); "CODE" 31035; "PROCEDURE" ELMMAT(A,B,X); "VALUE" X; "ARRAY" A,B; "REAL" X; "FOR" L:=1 "STEP" 1 "UNTIL" M "DO" ELMCOL(1,N,L,L,A,B,X); "PROCEDURE" INITIALIZE; "BEGIN" "INTEGER" I;"REAL" LBD; "SWITCH" TYPEODE:=NOTSPECIFIED2,PARABOLIC1,PARABOLIC2,HYPERBOLIC2; "IF" TYPE^=2 "AND" TYPE^=3 "THEN" TYPE:=1; "IF" TYPE^=2 "THEN" ORDER:=2 "ELSE" "IF" ORDER^=2 "THEN" ORDER:=1; I:=1; "GOTO" TYPEODE["IF" TYPE=1 "THEN" 1 "ELSE" TYPE+ORDER-1]; NOTSPECIFIED2: "FOR" LBD:=1/9,1/8,1/7,1/6,1/5,1/4,1/3,1/2,4.3 "DO" "BEGIN" LAMBDA[I]:=LBD; I:=I+1 "END"; "GOTO" EXIT; PARABOLIC1: "FOR"LBD:=.1418519249"-2,.3404154076"-2,.0063118569 ,.01082794375,.01842733851,.03278507942, .0653627415,.1691078577,156 "DO" "BEGIN" LAMBDA[I]:=LBD; I:=I+1 "END"; "GOTO" EXIT; PARABOLIC2: "FOR" LBD:=.3534355908"-2,.8532600867"-2,.015956206 ,.02772229155,.04812587964,.08848689452, .1863578961,.5,64 "DO" "BEGIN" LAMBDA[I]:=LBD; I:=I+1 "END"; "GOTO" EXIT; HYPERBOLIC2: "FOR" LBD:=1/8,1/20,5/32,2/17,17/80,5/22,11/32,1/2, 8 "DO" "BEGIN" LAMBDA[I]:=LBD; I:=I+1 "END"; "GOTO" EXIT; "COMMENT" 1SECTION : 5.2.1.1.3 (NOVEMBER 1976) PAGE 7 ; EXIT: SIG:=SIGN(TE-T) "END" INITIALIZE; "PROCEDURE" DIFFERENCE SCHEME; "BEGIN" "INTEGER" I;"REAL" MLT; DER(T,U,DU); "FOR" I:=1 "STEP" 1 "UNTIL" 8 "DO" "BEGIN" MLT:=LAMBDA[I]*TAU; DUPMAT(1,N,1,M,UH,U); ELMMAT(UH,DU,MLT); DER(T+MLT,UH,DU) "END"; ELMMAT(U,DU,TAU); T:="IF" LAST "THEN" TE "ELSE" T+TAU; "END" DIFFERENCE SCHEME; INITIALIZE; LAST:="FALSE"; STEP: TAU:=("IF" SPR=0 "THEN" ABS(TE-T) "ELSE" ABS(LAMBDA[9]/SPR))*SIG; "IF" T+TAU >= TE "EQUIV" TAU>=0 "THEN" "BEGIN" TAU:=TE-T;LAST:="TRUE" "END"; DIFFERENCE SCHEME ; OUT; "IF" "NOT" LAST "THEN" "GOTO" STEP "END" ARKMAT; "EOP" 1SECTION : 3.3.1.1.3.1 (NOVEMBER 1976) PAGE 1 AUTHORS: J.J.G. ADMIRAAL, A.C. YSSELSTEIN. CONTRIBUTOR: J.J.G. ADMIRAAL. INSTITUTE: UNIVERSITY OF AMSTERDAM. RECEIVED: 761101. BRIEF DESCRIPTION: THIS SECTION CONTAINS THREE PROCEDURES FOR SORTING THE ELEMENTS OF A VECTOR AND CORRESPONDINGLY PERMUTING THE ELEMENTS OF A VECTOR OR A MATRIX ROW. A) MERGESORT DELIVERS A PERMUTATION OF INDICES CORRESPONDING TO SORTING THE ELEMENTS OF A GIVEN VECTOR INTO NON-DECREASING ORDER. B) VECPERM PERMUTES THE ELEMENTS OF A GIVEN VECTOR CORRESPONDING TO A GIVEN PERMUTATION OF INDICES. C) ROWPERM PERMUTES THE ELEMENTS OF A GIVEN ROW OF A MATRIX CORRESPONDING TO A GIVEN PERMUTATION OF INDICES. KEYWORDS: SORTING, PERMUTING. SUBSECTION: MERGESORT. CALLING SEQUENCE: THE DECLARATION OF THE PROCEDURE IN THE CALLING PROGRAM READS: "PROCEDURE" MERGESORT(VEC1,VEC2,LOW,UPP); "VALUE" LOW,UPP;"INTEGER" LOW,UPP;"ARRAY" VEC1,VEC2; "CODE" 36405; THE MEANING OF THE FORMAL PARAMETERS IS: VEC1 :; "ARRAY" VEC1[LOW:UPP]; ENTRY: THE VECTOR TO BE SORTED INTO NONDECREASING ORDER; EXIT: THE CONTENTS OF VEC1 ARE LEFT INVARIANT; VEC2 :; "INTEGER""ARRAY" VEC2[LOW:UPP]; EXIT: THE PERMUTATION OF INDICES CORRESPONDING TO SORTING THE ELEMENTS OF VEC1 INTO NON-DECREASING ORDER; LOW : ; THE LOWER INDEX OF THE ARRAYS VEC1 AND VEC2; UPP : ; THE UPPER INDEX OF THE ARRAYS VEC1 AND VEC2; 1SECTION : 3.3.1.1.3.1 (NOVEMBER 1976) PAGE 2 PROCEDURES USED: NONE. REQUIRED CENTRAL MEMORY: ONE LOCAL INTEGER ARRAY OF LENGTH N,WHERE N = UPP - LOW + 1. RUNNING TIME: AVERAGE PROPORTIONAL TO N * LN(N). LANGUAGE: ALGOL 60. METHOD AND PERFORMANCE: SORTING BY MERGING. ([1],[2]) EXAMPLE OF USE: THE PROCEDURE MERGESORT IS USED IN SYMEIGIMP (SECTION 3.3.1.1.3.3). SUBSECTION: VECPERM. CALLING SEQUENCE: THE DECLARATION OF THE PROCEDURE IN THE CALLING PROGRAM READS: "PROCEDURE" VECPERM(PERM,LOW,UPP,VECTOR); "VALUE" LOW,UPP; "INTEGER" LOW,UPP; "INTEGER" "ARRAY" PERM;"ARRAY" VECTOR; "CODE" 36404; THE MEANING OF THE FORMAL PARAMETERS IS: PERM :; "INTEGER" "ARRAY" PERM[LOW:UPP]; ENTRY: A GIVEN PERMUTATION (E.G.AS PRODUCED BY MERGESORT) OF THE NUMBERS IN THE ARRAY VECTOR; LOW :; THE LOWER INDEX OF THE ARRAYS PERM AND VECTOR; UPP :; THE UPPER INDEX OF THE ARRAYS PERM AND VECTOR; VECTOR :; "ARRAY" VECTOR[LOW:UPP]; ENTRY: THE REAL VECTOR TO BE PERMUTED; EXIT: THE PERMUTED VECTOR ELEMENTS. PROCEDURE USED: NONE. REQUIRED CENTRAL MEMORY : ONE LOCAL BOOLEAN ARRAY OF LENGTH N IS DECLARED. RUNNING TIME: PROPORTIONAL TO N. LANGUAGE: ALGOL 60. EXAMPLE OF USE: THE PROCEDURE VECPERM IS USED IN SYMEIGIMP; (SECTION 3.3.1.1.3.3). 1SECTION : 3.3.1.1.3.1 (NOVEMBER 1976) PAGE 3 SUBSECTION: ROWPERM. CALLING SEQUENCE: THE DECLARATION OF THE PROCEDURE IN THE CALLING PROGRAM READS: "PROCEDURE" ROWPERM(PERM,LOW,UPP,I,MATRIX); "VALUE" LOW,UPP,I;"INTEGER" LOW,UPP,I; "INTEGER" "ARRAY" PERM;"ARRAY" MATRIX; "CODE" 36403; THE MEANING OF THE FORMAL PARAMETERS IS: PERM :; "INTEGER" "ARRAY" PERM[LOW:UPP]; ENTRY: A GIVEN PERMUTATION (E.G. AS PRODUCED BY MERGESORT) OF THE NUMBERS IN THE ARRAY VECTOR; LOW :; THE LOWER INDEX OF THE ARRAY PERM; UPP :; THE UPPER INDEX OF THE ARRAY PERM; I :; THE ROW INDEX OF THE MATRIX ELEMENTS; MATRIX :: "ARRAY" MATRIX [I :I, LOW : UP]; ENTRY:MATRIX [I, LOW : UP] SHOULD CONTAIN THE ELEMENTS TO BE PERMUTED; EXIT: MATRIX[I, LOW : UP] CONTAINS THE ROW OF PERMUTED ELEMENTS; PROCEDURES USED: NONE. REQUIRED CENTRAL MEMORY: ONE LOCAL BOOLEAN ARRAY OF LENGTH N IS DECLARED. RUNNING TIME: PROPORTIONAL TO N, WHERE N = UPP - LOW + 1. LANGUAGE: ALGOL 60. EXAMPLE OF USE: THE PROCEDURE ROWPERM IS USED IN SYMEIGIMP. (SECTION 3.3.1.1.3.3). REFERENCES: [1] D.E. KNUTH , THE ART OF COMPUTER PROGRAMMING, VOL. 3/ SORTING AND SEARCHING,ADDISON-WESLEY 1973. (SECTION 5.2.4 P.159-173). [2] A.V. AHO, J.E. HOPCROFT & J.D. ULLMAN, THE DESIGN AND ANALYSIS OF COMPUTER ALGORITHMS, ADDISON-WESLEY 1974. (SECTION I.M, P65-67). 1SECTION : 3.3.1.1.3.1 (NOVEMBER 1976) PAGE 4 SOURCE TEXTS: "CODE" 36405; "PROCEDURE" MERGESORT(A,P,LOW,UP);"VALUE" LOW,UP; "INTEGER" LOW,UP;"ARRAY" A;"INTEGER" "ARRAY" P; "BEGIN" "INTEGER" I,L,R,PL,PR,LO,STEP,STAP,UMLP1,UMSP1,REST,RESTV; "BOOLEAN" ROUT,LOUT; "INTEGER" "ARRAY" HP[LOW:UP]; "PROCEDURE" MERGE(LO,LS,RS);"VALUE" LO,LS,RS;"INTEGER" LO,LS,RS; "BEGIN" L:=LO; R:=LO+LS; LOUT:=ROUT:="FALSE"; "FOR" I:=LO,I+1 "WHILE" ^(LOUT "OR" ROUT) "DO" "BEGIN" PL:=P[L];PR:=P[R];"IF" A[PL]>A[PR] "THEN" "BEGIN" HP[I]:=PR;R:=R+1;ROUT:=R=LO+LS+RS "END" "ELSE" "BEGIN" HP[I]:=PL;L:=L+1;LOUT:=L=LO+LS "END" "END" FOR I; "IF" ROUT "THEN" "BEGIN" "FOR" I:=LO+LS-1 "STEP" -1 "UNTIL" L "DO" P[I+RS]:=P[I];R:=L+RS "END"; "FOR" I:=R-1 "STEP" -1 "UNTIL" LO "DO" P[I]:=HP[I]; "END" MERGE; "FOR" I:=LOW "STEP" 1 "UNTIL" UP "DO" P[I]:=I;RESTV:=0; UMLP1:=UP-LOW+1; "FOR" STEP:=1, STEP*2 "WHILE" STEP < UMLP1 "DO" "BEGIN" STAP:=2*STEP;UMSP1:=UP-STAP+1; "FOR" LO:=LOW "STEP" STAP "UNTIL" UMSP1 "DO" MERGE(LO,STEP,STEP); REST:=UP-LO+1; "IF" REST>RESTV & RESTV>0 "THEN" MERGE(LO,REST-RESTV,RESTV); RESTV:=REST "END" FOR STEP "END" MERGESORT; 1SECTION : 3.3.1.1.3.1 (NOVEMBER 1976) PAGE 5 "CODE" 36404; "PROCEDURE" VECPERM(PERM,LOW,UPP,VECTOR);"VALUE" LOW,UPP; "INTEGER" LOW,UPP;"INTEGER" "ARRAY" PERM;"REAL" "ARRAY" VECTOR; "BEGIN" "INTEGER" T,J,K;"REAL" A;"BOOLEAN" "ARRAY" TODO[LOW:UPP]; "FOR" T:=LOW "STEP" 1 "UNTIL" UPP "DO" TODO[T]:="TRUE"; "FOR" T:=LOW "STEP" 1 "UNTIL" UPP "DO" "BEGIN" "IF" TODO[T] "THEN" "BEGIN" K:=T;A:=VECTOR[K]; "FOR" J:=PERM[K] "WHILE" J^=T "DO" "BEGIN" VECTOR[K]:=VECTOR[J];TODO[K]:="FALSE";K:=J "END";VECTOR[K]:=A;TODO[K]:="FALSE" "END" CYCLE; "END" FOR T; "END" VECPERM; "CODE" 36403; "PROCEDURE" ROWPERM(PERM,LOW,UPP,I,MAT);"VALUE" LOW,UPP,I; "INTEGER" LOW,UPP,I;"INTEGER" "ARRAY" PERM;"REAL" "ARRAY" MAT; "BEGIN" "INTEGER" T,J,K;"REAL" A;"BOOLEAN" "ARRAY" TODO[LOW:UPP]; "FOR" T:=LOW "STEP" 1 "UNTIL" UPP "DO" TODO[T]:="TRUE"; "FOR" T:=LOW "STEP" 1 "UNTIL" UPP "DO" "BEGIN" "IF" TODO[T] "THEN" "BEGIN" K:=T;A:=MAT[I,K]; "FOR" J:=PERM[K] "WHILE" J^=T "DO" "BEGIN" MAT[I,K]:=MAT[I,J];TODO[K]:="FALSE";K:=J "END";MAT[I,K]:=A;TODO[K]:="FALSE" "END" CYCLE; "END" FOR T; "END" ROWPERM; 1SECTION : 3.3.1.1.3.2 (NOVEMBER 1976) PAGE 1 AUTHOR/CONTRIBUTOR: J.J.G. ADMIRAAL. INSTITUTE: UNIVERSITY OF AMSTERDAM. RECEIVED: 761101. BRIEF DESCRIPTION: THE PROCEDURE ORTHOG ORTHOGONALIZES SOME ADJACENT MATRIX COLUMNS ACCORDING TO THE MODIFIED GRAM SCHMIDT METHOD (SEE [1]). KEYWORDS: MATRIX COLUMNS, MODIFIED GRAM SCHMIDT ORTHOGONALIZATION. CALLING SEQUENCE: THE DECLARATION OF THE PROCEDURE IN THE CALLING PROGRAM READS: "PROCEDURE" ORTHOG(N,LC,UC,X); "VALUE" N,LC,UC; "INTEGER" N,LC,UC;"ARRAY" X; "CODE" 36402; THE MEANING OF THE FORMAL PARAMETERS IS: N :; THE ORDER OF THE MATRIX X; LC :; THE LOWER COLUMN INDEX OF THE MATRIX COLUMNS; UC :; THE UPPER COLUMN INDEX OF THE MATRIX COLUMNS; X :; "ARRAY" X[1:N,LC:UC]; ENTRY: THE MATRIX COLUMNS,TO BE ORTHOGONALIZED; EXIT: THE ORTHOGONALIZED MATRIX COLUMNS. PROCEDURES USED: TAMMAT = CP34014, ELMCOL = CP34023. REQUIRED CENTRAL MEMORY: NO LOCAL ARRAYS ARE DECLARED. RUNNING TIME: PROPORTIONAL TO N**3. 1SECTION : 3.3.1.1.3.2 (NOVEMBER 1976) PAGE 2 LANGUAGE: ALGOL 60. METHOD AND PERFORMANCE: THE MODIFIED GRAM SCHMIDT METHOD (SEE [1],CHAPTER 4.54). EXAMPLE OF USE: THE PROCEDURE ORTHOG IS USED IN SYMEIGIMP. (SECTION 3.3.1.1.3.3). REFERENCES: [1] J.H. WILKINSON. THE ALGEBRAIC EIGENVALUE PROBLEM. CLARENDON PRESS,OXFORD,1965. SOURCE TEXT: "CODE" 36402; "PROCEDURE" ORTHOG(N,LC,UC,X);"VALUE" N,LC,UC; "INTEGER" N,LC,UC;"ARRAY" X; "BEGIN" "INTEGER" I,J,K; "REAL" NORMX; "REAL" "PROCEDURE" TAMMAT(L,U,I,J,A,B); "CODE" 34014; "PROCEDURE" ELMCOL(L,U,I,J,A,B,X); "CODE" 34023; "FOR" J:=LC "STEP" 1 "UNTIL" UC "DO" "BEGIN" NORMX:=SQRT(TAMMAT(1,N,J,J,X,X)); "FOR" I:=1 "STEP" 1 "UNTIL" N "DO" X[I,J]:=X[I,J]/NORMX; "FOR" K:=J+1 "STEP" 1 "UNTIL" UC "DO" ELMCOL(1,N,K,J,X,X,-TAMMAT(1,N,K,J,X,X)) "END" "END" ORTHOG;