GENERIC MODULE Matrix(R, V);

Arithmetic for Modula-3, see doc for details

<* UNUSED *>
CONST
  Module = "Matrix.";

PROCEDURE New (m, n: CARDINAL; ): T =
  BEGIN
    RETURN NEW(T, m, n);
  END New;

PROCEDURE FromArray (READONLY x: TBody; ): T =
  VAR
    m := NUMBER(x);
    n := NUMBER(x[0]);
    z := NEW(T, m, n);
  BEGIN
    z^ := x;
    RETURN z;
  END FromArray;

PROCEDURE FromMatrixArray (READONLY x: TMBody; ): T =
  BEGIN
    IF NUMBER(x) = 0 OR NUMBER(x[0]) = 0 THEN
      RETURN New(0, 0);
    ELSE
      VAR m, n: CARDINAL := 0;
      BEGIN
        (* check matching row numbers and sum them up *)
        FOR i := 0 TO LAST(x) DO
          WITH size = NUMBER(x[i, 0]^) DO
            FOR j := 1 TO LAST(x[0]) DO
              <* ASSERT size = NUMBER(x[i, j]^), "Row numbers don't match." *>
            END;
            INC(m, size);
          END;
        END;
        (* check matching column numbers and sum them up *)
        FOR j := 0 TO LAST(x[0]) DO
          WITH size = NUMBER(x[0, j][0]) DO
            FOR i := 1 TO LAST(x) DO
              <* ASSERT size = NUMBER(x[i, j][0]),
                          "Column numbers don't match." *>
            END;
            INC(n, size);
          END;
        END;

        VAR
          z                := New(m, n);
          iz, jz: CARDINAL := 0;
        BEGIN
          FOR i := 0 TO LAST(x) DO
            jz := 0;
            FOR j := 0 TO LAST(x[0]) DO
              WITH y = x[i, j] DO
                FOR k := 0 TO LAST(y^) DO
                  SUBARRAY(z[iz + k], jz, NUMBER(y[0])) := y[k];
                END;
                INC(jz, NUMBER(y[0]));
              END;
            END;
            INC(iz, NUMBER(x[i, 0]^));
          END;
          RETURN z;
        END;
      END;
    END;
  END FromMatrixArray;

PROCEDURE RowFromArray (READONLY x: V.TBody; ): T =
  VAR z := NEW(T, 1, NUMBER(x));
  BEGIN
    z[0] := x;
    RETURN z;
  END RowFromArray;

PROCEDURE ColumnFromArray (READONLY x: V.TBody; ): T =
  VAR z := NEW(T, NUMBER(x), 1);
  BEGIN
    FOR i := 0 TO LAST(x) DO z[i, 0] := x[i]; END;
    RETURN z;
  END ColumnFromArray;

PROCEDURE DiagonalFromArray (READONLY x: V.TBody; ): T =
  VAR z := NEW(T, NUMBER(x), NUMBER(x));
  BEGIN
    FOR i := FIRST(x) TO LAST(x) DO
      z[i, i] := x[i];
      FOR j := FIRST(x) TO i - 1 DO
        z[i, j] := R.Zero;
        z[j, i] := R.Zero;
      END;
    END;
    RETURN z;
  END DiagonalFromArray;

PROCEDURE RowFromVector (x: V.T; ): T =
  BEGIN
    RETURN RowFromArray(x^);
  END RowFromVector;

PROCEDURE ColumnFromVector (x: V.T; ): T =
  BEGIN
    RETURN ColumnFromArray(x^);
  END ColumnFromVector;

PROCEDURE DiagonalFromVector (x: V.T; ): T =
  BEGIN
    RETURN DiagonalFromArray(x^);
  END DiagonalFromVector;

PROCEDURE FromScalar (x: R.T; ): T =
  VAR z := NEW(T, 1, 1);
  BEGIN
    z[0, 0] := x;
    RETURN z;
  END FromScalar;

PROCEDURE Copy (x: T; ): T =
  VAR
    m := NUMBER(x^);
    n := NUMBER(x[0]);
    z := NEW(T, m, n);
  BEGIN
    z^ := x^;
    RETURN z;
  END Copy;

PROCEDURE Cyclic (x: V.T; size: CARDINAL; shift: INTEGER; ): T =
  VAR
    z             := New(size, NUMBER(x^));
    rem: CARDINAL;
  BEGIN
    shift := shift MOD NUMBER(x^);
    rem := NUMBER(x^) - shift;
    IF size > 0 THEN
      z[0] := x^;
      FOR i := 1 TO LAST(z^) DO
        SUBARRAY(z[i], 0, shift) := SUBARRAY(z[i - 1], rem, shift);
        SUBARRAY(z[i], shift, rem) := SUBARRAY(z[i - 1], 0, rem);
      END;
    END;
    RETURN z;
  END Cyclic;

PROCEDURE GetRow (x: T; k: CARDINAL; ): V.T =
  VAR y := V.New(NUMBER(x[0]));
  BEGIN
    y^ := x[k];
    RETURN y;
  END GetRow;

PROCEDURE GetColumn (x: T; k: CARDINAL; ): V.T =
  VAR y := V.New(NUMBER(x^));
  BEGIN
    FOR j := 0 TO LAST(y^) DO y[j] := x[j, k]; END;
    RETURN y;
  END GetColumn;

PROCEDURE Apply (x: T; f: ApplyFtn; ) =
  BEGIN
    FOR i := 0 TO LAST(x^) DO
      FOR j := 0 TO LAST(x[0]) DO f(x[i, j]); END;
    END;
  END Apply;

PROCEDURE Map (x: T; f: MapFtn; ): T =
  VAR y := NEW(T, NUMBER(x^), NUMBER(x[0]));
  BEGIN
    FOR i := 0 TO LAST(x^) DO
      FOR j := 0 TO LAST(x[0]) DO y[i, j] := f(x[i, j]); END;
    END;
    RETURN y;
  END Map;

PROCEDURE ReduceRows (x: T; f: ReduceFtn; READONLY init: V.TBody; ): V.T =
  VAR y := NEW(V.T, NUMBER(init));
  BEGIN
    y^ := init;
    FOR i := 0 TO LAST(x^) DO
      FOR j := 0 TO LAST(x[0]) DO x[i, j] := f(y[i], x[i, j]); END;
    END;
    RETURN y;
  END ReduceRows;

PROCEDURE ReduceColumns (x: T; f: ReduceFtn; READONLY init: V.TBody; ):
  V.T =
  VAR y := NEW(V.T, NUMBER(init));
  BEGIN
    y^ := init;
    FOR i := 0 TO LAST(x^) DO
      FOR j := 0 TO LAST(x[0]) DO x[i, j] := f(y[j], x[i, j]); END;
    END;
    RETURN y;
  END ReduceColumns;

BEGIN
END Matrix.