The following set of procedures can be used for general Algol arrays. The notation used is as follows:-
A,B,C are real arrays, B and C (where used) must have the same number of dimensions as A and corresponding subscripts must have the same number of possible values. For example:-
A [1:3, 10:13, 19:25], B [2 :4, -7:-4, 6:12], C [0:2, 3:6, -7:-1]
is correct, but
A[1:2, 1:2], B[1:2, 1:2], C[1:2, 1:3] and A [1:2, 1:1], B[1:2, 1:1], C[1:2 ]
are incorrect.
If the arrays substituted for A, B, C do not satisfy the above restrictions, there will be a run-time fault. Any or all of the arrays substituted for A, B, C may be the same. Declarations for A, B, C are omitted throughout, in this document, but are of course included in the actual procedures.
The procedures have code bodies.
integer procedure dim(A);
Sets dim = dimension of the array A.
integer procedure bound (A,n); value n;integer n;
Sets bound = nth, bound of the array A. This bound will be the upper or lower bound of the mth subscript according as n = 2m -1,2m . If the nth bound does not exist, there will be a run-time fault.
procedure readarray(A);
Reads data and stores it by rows in A.
procedure printarray (A,m,n); value m,n; integer m,n;
Outputs the a-dimensional array A by rows in style m,n, inserting (2i+1) blank lines when the (a-i)th subscript is incremented.
procedure xarray(A,x); value x; real x;
Sets each element of A equal to x
procedure addxarray (A, B, x, C) ;value x; real x;
Sets A= B + x*C. Efficient code is obeyed for the special cases:
x = 0 giving A = B x = 1 giving A = B + C C ≠ B, x = -1 giving A = B - C C = B, x = -1 giving A = null array C = B, x = -2 giving A = -B C = B giving A = (1 + x) B
A set of additional procedures which call procedures in the ICT5 set:-
procedure nullarray (A);
Sets A = null array = xarray (A,0);
procedure equarray (A,B);
Sets A = B = addxarray (A,B,0,B);
procedure negarray (A,B);
Sets A = -B = addxarray (A,B,-2,B);
procedure addarray (A,B,C);
Sets A=B+C = addxarray (A,B,1,C);
procedure subarray (A,B,C) ;
Sets A = B-C = addxarray (A,B,-1,C);
procedure outarray (A) ;
Outputs the array A = printarray (A,0,10);
ICT7 to ICT11 define a set of procedures which operate on matrices. The notation used is as follows:-
A, B, C are matrices of sizes (a1 × a2), (b1 × b2), (c1 × c2) respectively. D is a diagonal matrix of size (d × d) stored as a vector of length d, I is the unit matrix (with size determined by context).
If A, B, C are not matrices, D is not a vector, or any of the restrictions specified for each procedure are not satisfied, then there will be a run-time fault. Any or all of the matrices substituted for A, B, C may be the same unless stated otherwise. Declarations for A, B, C, D are omitted throughout.
procedure transmx(A);
Restrictions: a1 = a2
Sets A = transpose (A).
procedure equtransmx(A,B);
Restrictions: a1 = b2, a2 = b1
Sets B = transpose (A)
procedure xunitmx(A,x); value x; real x;
Restrictions: a1 = a2
Sets A = x*I.
procedure addxunitmx(A,B,x); value x; real x;
Restrictions: a1 = a2 = b1 = b2
Sets A = B + x * I.
Efficient code is obeyed in the special case A = B giving A = A + x * I
procedure unitmx(A);
Restrictions: a1 = a2
Sets A = I = xunitmx(A,1);
procedure addunitmx(A,B);
Restrictions: a1 = a2 = b1 = b2
Sets A = B + I = addxunitmx(A,B,1);
procedure subunitmx(A,B);
Restrictions: a1 = a2 = b1 = b2
Sets A= B - I = addxunitmx(A,B,-1);
procedure gendiagmx(A,B,D,x, type);
value x, type; real x; integer type;
Restrictions: a1 = a2 = b1 = b2 = d
If type = 0, sets D = x * (leading diagonal of A)
B is not used and should be set as A.
If type ≠ 0, sets A = B + x * D.
Efficient code is obeyed for the special cases:-
type ≠ 0, A = B, x = 1 giving A = A + D type ≠ 0, A = B, x = -1 giving A = A - D type ≠ 0, A = B giving A = A + x*D type = 0, x = 1 giving D = (l. d. of A) type = 0, x = -1 giving D = -(l. d. of A)
procedure addiagmx(A,B,D);
Restrictions: a1 = a2 = b1 = b2 = d
Sets A = B + D = gendiagmx(A,B,D,1,1);
procedure subdiagmx(A,B,D);
Restrictions: a1 = a2 = b1 = b2 = d
Sets A = B - D = gendiamx(A,B,D,-1,1);
procedure equidiagmx(D,A);
Restrictions: a1 = a2 = d
Sets D = (leading diagonal of A) = gendiagmx(A,A,D,1,0);
procedure negdiagmx(D,A);
Restrictions: a1 = a2 = d
Sets D = (leading diagonal of A) = gendiagmx(A,A,D,-1,0);
procedure addxdiagmx(A,B,x,D); real x;
Restrictions: a1 = a2 = b1 = b2 =d
Sets A = B + x D = gendiagmx(A,B,D,x,1);
procedure xdiagmx(D,x,A); real x;
Restrictions: a1 = a2 = b1 = b2 =d
Sets D=x * (leading diagonal of A) = gendlagmx(A,A,D,x,0);
procedure multmx(A,B,C);
Restrictions: A, B, C different a1=b1 a2=c2 b2=c1
Sets A = B*C = genmultmx(A,B,C,"ex");
procedure multtransmx(A,B,C);
Restrictions: A, B, C different a1=b1 a2=c1 b2=x2
Sets A = B* transpose(C) = genmultmx(A,B,C,"ext");
procedure transmultmx(A,B,C);
Restrictions: A, B, C different a1=b2 a2=c2 1=c1
Sets A = transpose (B)* C = genmultmx(A,B,C, "etx");
procedure genmultmx(A,B,C,s); string s;
Restrictions: A must be different from both B and C, but B and C may be the same
The operation is determined by the string s which must be one of the following:
ex giving A = B*C Restriction: a1=b1 a2=c2 b2=c1 ex giving A = -B*C Restriction: a1=b1 a2=c2 b2=c1 px giving A = A + B*C Restriction: a1=b1 a2=c2 b2=c1 mx glvlng A = A - B*C Restriction: a1=b1 a2=c2 b2=c1 ext giving A = B*transpose(C) Restriction: a1=b1 a2=c1 b2=c2 ext giving A = -B*transpose(C) Restriction: a1=b1 a2=c1 b2=c2 pxt giving A = A + B*transpose(C) Restriction: a1=b1 a2=c1 b2=c2 mxt giving A = A- B*transpose( C) Restriction: a1=b1 a2=c1 b2=c2 etx giving A = transpose(B)*C Restriction: a1=b2 a2=c2 b1=c1 etx giving A = -transpose (B)*C Restriction: a1=b2 a2=c2 b1=c1 ptx giving A = A + transpose(B)*C Restriction: a1=b2 a2=c2 b1=c1 mtx giving A = A - transpose(B)*C Restriction: a1=b2 a2=c2 b1=c1
if s is not any of the above, then the action is undefined.
A set of procedures for doing matrix inversion, division and finding determinant.
real procedure invert (A,eps);
value eps; real eps;
Restrictions: a1 = a2
Inverts the matrix A on itself by Gauss-Jordan elimination. If, during elimination, the absolute value of any pivot is eps, there will be a run-time fault. If the inversion is successful, invert will contain the determinant of the original matrix.
real procedure invert part(A,p,q,n,eps);
value p,q,n,eps;
integer p,q,n; real eps;
Restrictions: a1≤p a1≤q a2≥p+n-1 a2≥q+n-1
This procedure acts in the same way as invert except that it inverts the submatrix of A given by A[i,j], p≤i≤p+m-1, q≤j≤q+n-1.
procedure invmx;
This procedure is required for use with invert and/or invertpart.
real procedure matdiv(A,B,eps);
value eps; real eps;
Restrictions: a1 = a2 = b1
Uses a Gauss elimination method to set B = inverse(A)*B.
If, during elimination, the absolute value of any pivot is eps, there will be a run-time fault. Otherwise matdiv will be set to the determinant of A. B may be a matrix or vector. A is destroyed.
real procedure det(A,eps);
value eps; real eps;
Uses a Gauss elimination method to set det to the determinant of A. If, during elimination, the absolute value of any pivot is ≤ eps, then det = 0. A is destroyed.
procedure divmx;
This procedure is required for use with matdiv and/or det.
procedure unfault invmx;
Alters procedure invmx so that if invertor invertpart attempts to invert a singular matrix (i.e. the absolute value of some pivot is ≤ eps), then instead of giving a run-time fault, the procedure exits normally and invert or invertpart is set = 0.
procedure unfault divmx;
Alters procedure divmx so that if matdiv attempts to divide by a singular matrix (i.e. the absolute value of some pivot is ≤ eps), then instead of giving a run-time fault, the procedure exits normally and 'matdiv' is set = 0.
If invmx or divmx is altered by the appropriate above procedure, it cannot be restored to its original condition.