## arithmetic/src/linearalgebra/matrix/MatrixFast.mg

```GENERIC MODULE MatrixFast(R, V, MS, MD);
```
Arithmetic for Modula-3, see doc for details
```
IMPORT Arithmetic;

<* INLINE *>
PROCEDURE AssertEqualSize (x, y: T; ) =
BEGIN
<* ASSERT NUMBER(x^) = NUMBER(y^) AND NUMBER(x[0]) = NUMBER(y[0]),
"Sizes of matrices must match." *>
END AssertEqualSize;

PROCEDURE IsZero (x: T; ): BOOLEAN =
BEGIN
FOR i := FIRST(x^) TO LAST(x^) DO
FOR j := FIRST(x[0]) TO LAST(x[0]) DO
IF x[i, j] # R.Zero THEN RETURN FALSE; END;
END;
END;
RETURN TRUE;
END IsZero;

PROCEDURE Equal (x, y: T; ): BOOLEAN =
BEGIN
AssertEqualSize(x, y);

FOR i := FIRST(x^) TO LAST(x^) DO
FOR j := FIRST(x[0]) TO LAST(x[0]) DO
IF x[i, j] # y[i, j] THEN RETURN FALSE; END;
END;
END;
RETURN TRUE;
END Equal;

PROCEDURE Add (x, y: T; ): T =
BEGIN
AssertEqualSize(x, y);

WITH z = NEW(T, NUMBER(x^), NUMBER(x[0])) DO
FOR i := FIRST(x^) TO LAST(x^) DO
FOR j := FIRST(x[0]) TO LAST(x[0]) DO
z[i, j] := x[i, j] + y[i, j];
END;
END;
RETURN z;
END;

PROCEDURE Sub (x, y: T; ): T =
BEGIN
AssertEqualSize(x, y);

WITH z = NEW(T, NUMBER(x^), NUMBER(x[0])) DO
FOR i := FIRST(x^) TO LAST(x^) DO
FOR j := FIRST(x[0]) TO LAST(x[0]) DO
z[i, j] := x[i, j] - y[i, j];
END;
END;
RETURN z;
END;
END Sub;

PROCEDURE Scale (x: T; y: R.T; ): T =
VAR z := NEW(T, NUMBER(x^), NUMBER(x[0]));
BEGIN
FOR i := FIRST(x^) TO LAST(x^) DO
FOR j := FIRST(x[0]) TO LAST(x[0]) DO z[i, j] := x[i, j] * y; END;
END;
RETURN z;
END Scale;

PROCEDURE Mul (x, y: T; ): T =
VAR
m    := NUMBER(x^);
n    := NUMBER(x[0]);
p    := NUMBER(y[0]);
z: T;

BEGIN
MS.AssertEqualWidth(NUMBER(y^), n);
z := NEW(T, m, p);
FOR i := FIRST(x^) TO LAST(x^) DO
FOR j := FIRST(y[0]) TO LAST(y[0]) DO
VAR sum := R.Zero;
BEGIN
FOR k := FIRST(x[0]) TO LAST(x[0]) DO
sum := sum + x[i, k] * y[k, j];
END;
z[i, j] := sum;
END;
END;
END;
RETURN z;
END Mul;

PROCEDURE MulV (x: T; y: V.T; ): V.T =
VAR z := NEW(V.T, NUMBER(x^));
BEGIN
MS.AssertEqualWidth(NUMBER(x[0]), NUMBER(y^));

FOR i := FIRST(x^) TO LAST(x^) DO
VAR sum := R.Zero;
BEGIN
FOR j := FIRST(x[0]) TO LAST(x[0]) DO
sum := sum + y[j] * x[i, j];
END;
z[i] := sum;
END;
END;
RETURN z;
END MulV;

PROCEDURE MulTV (x: T; y: V.T; ): V.T =
VAR z := NEW(V.T, NUMBER(x[0]));
BEGIN
MS.AssertEqualWidth(NUMBER(x^), NUMBER(y^));

FOR i := FIRST(x[0]) TO LAST(x[0]) DO
VAR sum := R.Zero;
BEGIN
FOR j := FIRST(x^) TO LAST(x^) DO sum := sum + x[j, i] * y[j]; END;
z[i] := sum;
END;
END;
RETURN z;
END MulTV;

PROCEDURE MulMAM (x: T; ): T =
VAR z := NEW(T, NUMBER(x[0]), NUMBER(x[0]));
BEGIN
FOR i := 0 TO LAST(x[0]) DO
FOR j := i TO LAST(x[0]) DO
VAR sum := R.Zero;
BEGIN
FOR k := 0 TO LAST(x^) DO sum := sum + x[k, i] * x[k, j]; END;
z[i, j] := sum;
z[j, i] := sum;
END;
END;
END;
RETURN z;
END MulMAM;

PROCEDURE MulMMA (x: T; ): T =
VAR z := NEW(T, NUMBER(x^), NUMBER(x^));
BEGIN
FOR i := 0 TO LAST(x^) DO
FOR j := i TO LAST(x^) DO
VAR sum := R.Zero;
BEGIN
FOR k := 0 TO LAST(x[0]) DO sum := sum + x[i, k] * x[j, k]; END;
z[i, j] := sum;
z[j, i] := sum;
END;
END;
END;
RETURN z;
END MulMMA;

PROCEDURE Trace (x: T; ): R.T =
VAR y := R.Zero;
BEGIN
FOR j := 0 TO MIN(LAST(x^), LAST(x[0])) DO y := y + x[j, j]; END;
RETURN y;
END Trace;
```

```PROCEDURE Determinant (x: T; ): R.T =
BEGIN
TRY
RETURN MD.LUDet(MD.LUFactor(x));
EXCEPT
| Arithmetic.Error =>
(* Although the determinant is always defined, the fast Gauß
algorithm fails sometimes.  In this cases we fall back to the
division-free algorithm. *)
RETURN MS.Determinant(x);
END;
END Determinant;

BEGIN
END MatrixFast.
```
```

```