language-Modula2-0.1: examples/Modula-2_Libraries/andrea-m2/lib/generic/matrixope.mod
IMPLEMENTATION MODULE MatrixOperations;
(* matrix operations using dynamic memory, see the corresponding definition *)
(* V1.2, J. Andrea, Jun.22/93 -add Duplicate *)
(* V1.1, J. Andrea, Jun.11/92 - recursive determinant for any size matrix *)
(* V1.0, J. Andrea, Mar.16/92 *)
(* This code may be freely used and distributed, it may not be sold. *)
FROM SYSTEM IMPORT ADDRESS, TSIZE;
FROM Storage IMPORT ALLOCATE, DEALLOCATE;
TYPE
Matrix = POINTER TO RECORD
rows, cols :CARDINAL; (* dimemsions of the matrix *)
start :ADDRESS; (* memory address of the first item *)
size :CARDINAL; (* actual size in bytes *)
END;
DataType = REAL;
DataPoint = POINTER TO DataType;
VAR
type_size :CARDINAL; (* #bytes in a single item *)
delta :DataType; (* small number for some operations *)
(* -------------------------------------------------------- *)
PROCEDURE InRange( a :Matrix; row, col :CARDINAL ) :BOOLEAN;
(* is the specified row/col item in this matrix ? *)
BEGIN
RETURN ( row >= 1 ) & ( row <= a^.rows ) & ( col >= 1 ) & ( col <= a^.cols );
END InRange;
(* -------------------------------------------------------- *)
PROCEDURE Offset( a :Matrix; row, col :CARDINAL ) :CARDINAL;
(* calculate the memory offset to the row/col item from the first item *)
BEGIN
RETURN type_size * ( ( row - 1 ) * a^.cols + ( col - 1 ) );
END Offset;
(* -------------------------------------------------------- *)
PROCEDURE Build( VAR a :Matrix; n_rows, n_columns :CARDINAL );
BEGIN
NEW( a );
(* get something *)
IF n_rows = 0 THEN n_rows := 1 END;
IF n_columns = 0 THEN n_columns := 1 END;
a^.rows := n_rows;
a^.cols := n_columns;
a^.size := n_rows * n_columns * type_size;
ALLOCATE( a^.start, a^.size );
END Build;
(* -------------------------------------------------------- *)
PROCEDURE Destroy( VAR a :Matrix );
BEGIN
DEALLOCATE( a^.start, a^.size );
DISPOSE( a );
END Destroy;
(* -------------------------------------------------------- *)
PROCEDURE Put( a :Matrix; row, col :CARDINAL; x :DataType );
VAR
adra :DataPoint;
BEGIN
IF InRange( a, row, col ) THEN
adra := a^.start + Offset( a, row, col );
adra^ := x;
END;
END Put;
(* -------------------------------------------------------- *)
PROCEDURE Get( a :Matrix; row, col :CARDINAL ) :DataType;
VAR
adra :DataPoint;
result :DataType;
BEGIN
IF InRange( a, row, col ) THEN
adra := a^.start + Offset( a, row, col );
result := adra^;
ELSE
result := 0.0;
END;
RETURN result;
END Get;
(* -------------------------------------------------------- *)
PROCEDURE Size( a :Matrix; VAR n_rows, n_cols :CARDINAL );
BEGIN
n_rows := a^.rows;
n_cols := a^.cols;
END Size;
(* -------------------------------------------------------- *)
PROCEDURE Min( a :Matrix ) :DataType;
VAR
result :DataType;
adra :DataPoint;
i, k, n :CARDINAL;
BEGIN
n := a^.rows * a^.cols;
adra := a^.start;
result := adra^;
k := type_size;
FOR i := 2 TO n DO
adra := a^.start + k;
k := k + type_size;
IF adra^ < result THEN result := adra^ END;
END;
RETURN result;
END Min;
(* -------------------------------------------------------- *)
PROCEDURE Max( a :Matrix ) :DataType;
VAR
result :DataType;
adra :DataPoint;
i, k, n :CARDINAL;
BEGIN
n := a^.rows * a^.cols;
adra := a^.start;
result := adra^;
k := type_size;
FOR i := 2 TO n DO
adra := a^.start + k;
k := k + type_size;
IF adra^ > result THEN result := adra^ END;
END;
RETURN result;
END Max;
(* -------------------------------------------------------- *)
PROCEDURE Compare( a, b :Matrix ) :BOOLEAN;
VAR
ok :BOOLEAN;
adra, adrb :DataPoint;
i, k, n :CARDINAL;
BEGIN
IF ( a^.rows = b^.rows ) & ( a^.cols = b^.cols ) THEN
ok := TRUE;
n := a^.rows * a^.cols;
i := 1;
k := 0;
WHILE ok & ( i <= n ) DO
adra := a^.start + k;
adrb := b^.start + k;
k := k + type_size;
ok := ABS( adra^ - adrb^ ) <= delta;
i := i + 1;
END;
ELSE
ok := FALSE;
END;
RETURN ok;
END Compare;
(* -------------------------------------------------------- *)
PROCEDURE Assign( a :Matrix; x :DataType );
VAR
adra :DataPoint;
i, k, n :CARDINAL;
BEGIN
n := a^.rows * a^.cols;
k := 0;
FOR i := 1 TO n DO
adra := a^.start + k;
adra^ := x;
k := k + type_size;
END;
END Assign;
(* -------------------------------------------------------- *)
PROCEDURE Ident( a :Matrix );
VAR
adra :DataPoint;
i, j, k, n :CARDINAL;
BEGIN
n := a^.rows * a^.cols;
j := 1; (* item number of first '1' *)
k := 0;
FOR i := 1 TO n DO
adra := a^.start + k;
IF i = j THEN
adra^ := 1.0;
j := j + a^.cols + 1; (* next item number which will be a '1' *)
ELSE
adra^ := 0.0;
END;
k := k + type_size;
END;
END Ident;
(* -------------------------------------------------------- *)
PROCEDURE Copy( a, b :Matrix );
VAR
adra, adrb :DataPoint;
i, k, n :CARDINAL;
BEGIN
IF ( a^.rows = b^.rows ) & ( a^.cols = b^.cols ) THEN
n := a^.rows * a^.cols;
k := 0;
FOR i := 1 TO n DO
adra := a^.start + k;
adrb := b^.start + k;
adrb^ := adra^;
k := k + type_size;
END;
END;
END Copy;
(* -------------------------------------------------------- *)
PROCEDURE Duplicate( a :Matrix; VAR b :Matrix );
BEGIN
Build( b, a^.rows, a^.cols );
Copy( a, b );
END Duplicate;
(* -------------------------------------------------------- *)
PROCEDURE Scale( a :Matrix; x :DataType );
VAR
adra :DataPoint;
i, k, n :CARDINAL;
BEGIN
n := a^.rows * a^.cols;
k := 0;
FOR i := 1 TO n DO
adra := a^.start + k;
adra^ := adra^ * x;
k := k + type_size;
END;
END Scale;
(* -------------------------------------------------------- *)
PROCEDURE Add( a, b :Matrix; VAR c :Matrix; VAR ok :BOOLEAN );
VAR
adra, adrb, adrc :DataPoint;
i, k, n :CARDINAL;
BEGIN
IF ( a^.rows # b^.rows ) OR ( a^.cols # b^.cols ) THEN
ok := FALSE;
ELSE
ok := TRUE;
Build( c, a^.rows, a^.cols );
n := a^.rows * a^.cols;
k := 0;
FOR i := 1 TO n DO
adra := a^.start + k;
adrb := b^.start + k;
adrc := c^.start + k;
adrc^ := adra^ + adrb^;
k := k + type_size;
END;
END;
END Add;
(* -------------------------------------------------------- *)
PROCEDURE Multiply( a, b :Matrix; VAR c :Matrix; VAR ok :BOOLEAN );
VAR
sum :DataType;
adra, adrb, adrc :DataPoint;
i, j, k :CARDINAL;
BEGIN
IF a^.cols # b^.rows THEN
ok := FALSE;
ELSE
ok := TRUE;
Build( c, a^.rows, b^.cols );
FOR i := 1 TO a^.rows DO
FOR j := 1 TO b^.cols DO
sum := 0.0;
FOR k := 1 TO a^.cols DO
adra := a^.start + Offset( a, i, k );
adrb := b^.start + Offset( b, k, j );
sum := sum + adra^ * adrb^;
END;
adrc := c^.start + Offset( c, i, j );
adrc^ := sum;
END;
END;
END;
END Multiply;
(* -------------------------------------------------------- *)
PROCEDURE Transpose( a :Matrix; VAR b :Matrix );
VAR
adra, adrb :DataPoint;
i, j :CARDINAL;
BEGIN
Build( b, a^.cols, a^.rows ); (* flip dimensions *)
FOR i := 1 TO a^.rows DO
FOR j := 1 TO a^.cols DO
adra := a^.start + Offset( a, i, j );
adrb := b^.start + Offset( b, j, i );
adrb^ := adra^;
END;
END;
END Transpose;
(* -------------------------------------------------------- *)
PROCEDURE Invert( b :Matrix; VAR a :Matrix; VAR ok :BOOLEAN ); (* note a <-> b *)
VAR
r, p :DataType;
adra, adrb, adrc :DataPoint;
i, j, k, n :CARDINAL;
BEGIN
ok := TRUE;
IF ( b^.rows # b^.cols ) OR ( b^.rows < 2 ) THEN
ok := FALSE; (* not square, forget it *)
ELSE
Copy( b, a );
n := a^.rows;
k := 1;
WHILE ok & ( k <= n ) DO
adra := a^.start + Offset( a, k, 1 ); (* pivot := mat(k,1) *)
p := adra^;
IF p <= 1.E-10 THEN
ok := FALSE; (* zero pivot point *)
ELSE
FOR j := 1 TO n - 1 DO
adra := a^.start + Offset( a, k, j+1 );
adrb := a^.start + Offset( a, k, j );
adrb^ := adra^ / p; (* mat(k,j) := mat(k,j+1)/p *)
END;
adra := a^.start + Offset( a, k, n );
adra^ := 1.0 / p; (* mat(k,n) := 1/p *)
FOR i := 1 TO n DO
IF i # k THEN
adra := a^.start + Offset( a, i, 1 );
r := adra^; (* r := mat(i,1) *)
FOR j := 1 TO n - 1 DO
adra := a^.start + Offset( a, k, j );
adrb := a^.start + Offset( a, i, j+1 );
adrc := a^.start + Offset( a, i, j );
adrc^ := adrb^ - r * adra^;
(* mat(i,j) := mat(i,j+1) - r * mat(k,j) *)
END;
adra := a^.start + Offset( a, k, n );
adrb := a^.start + Offset( a, i, n );
adrb^ := - r * adra^; (* mat(i,n) := - r * mat(k,n) *)
END;
END;
END;
k := k + 1;
END;
IF NOT ok THEN
Destroy( b );
END;
END;
END Invert;
(* -------------------------------------------------------- *)
PROCEDURE Determinant( a :Matrix; VAR d :DataType; VAR ok :BOOLEAN );
VAR
n, i, j, p, q, k :CARDINAL;
b :Matrix;
f, x, y, z :DataType;
adra, adrb :DataPoint;
BEGIN
ok := TRUE;
n := a^.rows;
IF n # a^.cols THEN
ok := FALSE;
d := 0.0;
ELSE
IF n = 1 THEN
adra := a^.start;
d := adra^;
ELSIF n = 2 THEN
adra := a^.start;
adrb := a^.start + Offset( a, 2, 2 );
d := adra^ * adrb^;
adra := a^.start + Offset( a, 2, 1 );
adrb := a^.start + Offset( a, 1, 2 );
d := d - adra^ * adrb^;
ELSE
z := 0.0;
f := 1.0;
(* create sub matricies, and found their determinants *)
Build( b, n-1, n-1 );
(* run across all the columns in the input matrix *)
FOR k := 1 TO n DO
adra := a^.start + Offset( a, 1, k );
x := adra^;
(* skip any zero element *)
IF x # 0.0 THEN
q := 1;
FOR j := 1 TO n DO
(* skip the column that matches the top selected column *)
IF k # j THEN
(* pick up all the rows in this whole column *)
FOR i := 2 TO n DO
p := i - 1;
adra := a^.start + Offset( a, i, j );
adrb := b^.start + Offset( b, p, q );
adrb^ := adra^;
END;
q := q + 1;
END;
END;
Determinant( b, y, ok );
z := z + f * x * y; (* add in the result from that submatrix *)
END;
f := - f; (* +-+-+-+-... is how the summing goes *)
END;
Destroy( b );
d := z; (* return the final result *)
END;
END;
END Determinant;
BEGIN
delta := 1.E-6;
type_size := TSIZE( DataType );
END MatrixOperations.