(*--------------------------------------------------------------------

:Program.       MATRIX.mod
:Author.        Michael Meyer
:Address.       Würmseestrasse 28, 81476 München
:EMail.         ZNet: M.MEYER@AMC.ZER.SUB.ORG
:EMail.         Fido: Michael Meyer 2:2480/21.100
:Version.       1.0
:Date.          6 Feb 1994
:Copyright.     1994 Michael Meyer
:Contents.      Dieses Modul unterstüzt das Rechnen mit Matrizen
:Language.      Oberon-2
:Translator.    Amiga Oberon 3.00d

--------------------------------------------------------------------*)

(* Matrix = elements( [0,  0] [0,  1] [0,  2] ... [0,  n-1]
 *                    [1,  0] [1,  1] [1,  2] ... [1,  n-1]
 *                    [2,  0] [2,  1] [2,  2] ... [2,  n-1]
 *                       .       .       .           .
 *                       .       .       .           .
 *                       .       .       .           .
 *                    [m-1,0] [m-1,1] [m-1,2] ... [m-1,n-1] )
 *)

MODULE MATRIX;

IMPORT BT * := BasicTypes;

TYPE
  MATRIX * = POINTER TO MATRIXDesc;
  MATRIXDesc * = RECORD (BT.GROUPDesc)
    m, n: LONGINT;
    elements: POINTER TO ARRAY OF ARRAY OF LONGREAL;
  END;

PROCEDURE Create * (m, n: LONGINT): MATRIX;
(*
 * require
 *   m>0 & n>0
 *)
VAR
  new: MATRIX;
  i, j: LONGINT;
BEGIN
  NEW(new);
  NEW(new.elements, m, n);
  new.m := m;
  new.n := n;
  FOR i := 0 TO m-1 DO
    FOR j := 0 TO n-1 DO
      new.elements[j, i] := 0.0;
    END;
  END;
  RETURN new;
END Create;

PROCEDURE NeutralMatrix * (m: LONGINT): MATRIX;
(*
 * require
 *   m>0 & n>0
 *   a.DimM() = a.DimN()
 *)
VAR
  new: MATRIX;
  i, j: LONGINT;
BEGIN
  NEW(new);
  NEW(new.elements, m, m);
  new.m := m;
  new.n := m;
  FOR i := 0 TO m-1 DO
    FOR j := 0 TO m-1 DO
      new.elements[j, i] := 0.0;
    END;
  END;
  FOR i := 0 TO m-1 DO
    new.elements[i, i] := 1.0;
  END;
  RETURN new;
END NeutralMatrix;

PROCEDURE (a: MATRIX) Get * (m, n: LONGINT): LONGREAL;
(*
 * require
 *   a.DimM()>m &
 *   a.DimN()>n;
 *   m>0 &
 *   n>0
 *)
BEGIN
  RETURN a.elements[m, n];
END Get;

PROCEDURE (a: MATRIX) Put * (m, n: LONGINT; x: LONGREAL);
(*
 * require
 *   a.DimM()>m &
 *   a.DimN()>n;
 *   m>0 &
 *   n>0
 *)
BEGIN
  a.elements[m, n] := x;
END Put;

PROCEDURE (a: MATRIX) DimM * (): LONGINT;
BEGIN
  RETURN a.m;
END DimM;

PROCEDURE (a: MATRIX) DimN * (): LONGINT;
BEGIN
  RETURN a.n;
END DimN;

PROCEDURE (a: MATRIX) Add * (b: BT.GROUP): BT.GROUP;
(*
 * require
 *   a.DimM() = b(MATRIX).DimM()
 *   a.DimN() = b(MATRIX).DimN()
 *)
VAR
  new: MATRIX;
  i, j: LONGINT;
BEGIN
  WITH b: MATRIX DO
    new := Create(a.m, a.n);
    FOR i := 0 TO a.m-1 DO
      FOR j := 0 TO a.n-1 DO
        new.elements[i, j] := a.elements[i, j] + b.elements[i, j]
      END;
    END;
  END;
  RETURN new;
END Add;

PROCEDURE (a: MATRIX) Sub * (b: BT.GROUP): BT.GROUP;
(*
 * require
 *   a.DimM() = b(MATRIX).DimM()
 *   a.DimN() = b(MATRIX).DimN()
 *)
VAR
  new: MATRIX;
  i, j: LONGINT;
BEGIN
  WITH b: MATRIX DO
    new := Create(a.m, a.n);
    FOR i := 0 TO a.m-1 DO
      FOR j := 0 TO a.n-1 DO
        new.elements[i, j] := a.elements[i, j] - b.elements[i, j]
      END;
    END;
  END;
  RETURN new;
END Sub;

PROCEDURE (a: MATRIX) Mul * (b: BT.GROUP): BT.GROUP;
(*
 * c = a * b
 *
 * require
 * a.DimN() = b(MATRIX).DimM()
 *
 * => c.DimM() = a.DimM()
 *    c.DimN() = b(MATRIX).DimN()
 *)
VAR
  new: MATRIX;
  i, j, o: LONGINT;
  c: LONGREAL;
BEGIN
  WITH b: MATRIX DO
    new := Create(a.m, b.n);
    FOR i := 0 TO a.m-1 DO
      FOR j := 0 TO b.n-1 DO
        c := 0.0;
        FOR o := 0 TO a.n-1 DO
          c := c + a.elements[i, o] * b.elements[o, j];
        END;
        new.elements[i, j] := c;
      END;
    END;
  END;
  RETURN new;
END Mul;

PROCEDURE (a: MATRIX) SkalarMul * (r: LONGREAL): BT.GROUP;
(*
 * multiply each element of a by r.
 *)
VAR
  new: MATRIX;
  i, j: LONGINT;
BEGIN
  WITH a: MATRIX DO
    new := Create(a.m, a.n);
    FOR i := 0 TO a.m-1 DO
      FOR j := 0 TO a.n-1 DO
        new.elements[i, j] := r * a.elements[i, j];
      END;
    END;
  END;
  RETURN new;
END SkalarMul;

PROCEDURE (a: MATRIX) isEqual * (b: BT.COMPAREABLE): BOOLEAN;
VAR
  i, j: LONGINT;
BEGIN
  WITH b: MATRIX DO
    IF ~((a.m = b.m) & (a.n = b.n)) THEN RETURN FALSE END;
    FOR i := 0 TO a.m-1 DO
      FOR j := 0 TO a.n-1 DO
        IF a.elements[i, j] # b.elements[i, j] THEN RETURN FALSE END;
      END;
    END;
  END;
  RETURN TRUE;
END isEqual;

PROCEDURE (a: MATRIX) isDimEqual * (b: BT.COMPAREABLE): BOOLEAN;
BEGIN
  WITH b: MATRIX DO
    IF ~((a.m = b.m) & (a.n = b.n)) THEN
      RETURN FALSE;
    ELSE
      RETURN TRUE;
    END;
  END;
END isDimEqual;

PROCEDURE (a: MATRIX) TransposeMatrix * (): BT.GROUP;
VAR
  i, j: LONGINT;
  new: MATRIX;
BEGIN
  new := Create(a.n, a.m);
  FOR i := 0 TO a.m-1 DO
    FOR j := 0 TO a.n-1 DO
      new.elements[j, i] := a.elements[i, j];
    END;
  END;
  RETURN new;
END TransposeMatrix;

PROCEDURE FindPivot (a: MATRIX; c, l : LONGINT): LONGINT;
(*
 * find Pivot-element
 * in column c
 * from line l on
 *
 * returns line of Pivot-element
 *)
VAR
  i : LONGINT;
  pivot, next : LONGREAL;
  line : LONGINT;
BEGIN
  pivot := ABS(a.elements[l, c]);
  line := l;
  FOR i := l+1 TO a.m-1 DO
    next := ABS(a.elements[i, c]);
    IF next > pivot THEN
      pivot := next;
      line := i;
    END;
  END;
RETURN line;
END FindPivot;

PROCEDURE ChangeTwoLines (a: MATRIX; l1, l2, c: LONGINT);
(*
 * change two lines
 * from column c to column a.DimN()-1
 *)
VAR
  i : LONGINT;
  x : LONGREAL;
BEGIN
  FOR i := c TO a.n-1 DO
    x := a.elements[l1, i];
    a.elements[l1, i] := a.elements[l2, i];
    a.elements[l2, i] := x;
  END;
END ChangeTwoLines;

PROCEDURE MulLineAddLine(a: MATRIX; l1, l2, c: LONGINT; s: LONGREAL);
(*
 * multiplicate each element of line l1 from column c with
 * skalar s and add the result to line l2
 *)
VAR
  i : LONGINT;
  x : LONGREAL;
BEGIN
  FOR i := c TO a.n-1 DO
    a.elements[l2, i] := a.elements[l2, i] + a.elements[l1, i] * s;
  END;
END MulLineAddLine;


PROCEDURE (a: MATRIX) Neg * (): BT.GROUP;
(*
 * require
 *   a.DimM() = a.DimN()
 *   a.Det() # 0
 *)
VAR
  new, newN, return: MATRIX;
  i, j, k: LONGINT;
  ln: LONGINT;
  t1, t2, skalar, pivot, det: LONGREAL;
  sign: BOOLEAN;
BEGIN
  new := a;
  newN := NeutralMatrix(a.m);
  return := Create(a.m, a.n);
(*------------------ Gauss ------------------*)
  sign := TRUE;
  FOR i := 0 TO new.m-2 DO
    ln := FindPivot(new, i, i);
    pivot := -1.0 * new.elements[ln, i];
    IF ~(ln = i) THEN
      sign := ~sign;
      ChangeTwoLines(new, i, ln, i);
      ChangeTwoLines(newN, i, ln, 0);
    END;
    FOR j := i+1 TO new.m-1 DO
      skalar := new.elements[j, i]/pivot;
      MulLineAddLine(new, i, j, i, skalar);
      MulLineAddLine(newN, i, j, 0, skalar);
    END;
  END;
  IF ~sign THEN
    FOR i := 0 TO new.n-1 DO
      new.elements[new.m-1, i] := -1.0 * new.elements[new.m-1, i];
      newN.elements[new.m-1, i] := -1.0 * newN.elements[new.m-1, i];
    END;

  END;
(*------------------ Jordan ------------------*)
  FOR i := new.m-1 TO 1 BY -1 DO
    t1 := -1.0 * new.elements[i, i];
    FOR j := i-1 TO 0 BY -1 DO
      t2 := new.elements[j, i] / t1;
      MulLineAddLine(new, i, j, j, t2);
      MulLineAddLine(newN, i, j, 0, t2);
    END;
  END;
  FOR i := 0 TO new.m-1 DO
    FOR j := 0 TO new.n-1 DO
      return.elements[i, j] := newN.elements[i, j]/new.elements[i, i];
    END;
  END;
  RETURN return;
END Neg;

PROCEDURE (a: MATRIX) Gauss * (): BT.GROUP;
(*
 * require
 *   a.DimM() = a.DimN()
 *)
VAR
  new: MATRIX;
  sign: BOOLEAN;  (* sign of the matrix: TRUE=plus and FALSE=minus *)
  i, j: LONGINT;
  ln: LONGINT;  (* line-number in which Pivot-element is *)
  skalar, pivot: LONGREAL;
BEGIN
  new := Create(a.m, a.n);
  new := a;
  sign := TRUE;
  FOR i := 0 TO new.m-2 DO
    ln := FindPivot(new, i, i);
    pivot := -1.0 * new.elements[ln, i];
    IF ~(ln = i) THEN
      sign := ~sign;
      ChangeTwoLines(new, i, ln, i);
    END;
    FOR j := i+1 TO new.m-1 DO
      skalar := new.elements[j, i]/pivot;
      MulLineAddLine(new, i, j, i, skalar);
    END;
  END;
  IF ~sign THEN
    FOR i := 0 TO new.n-1 DO
      new.elements[new.m-1, i] := -1.0 * new.elements[new.m-1, i];
    END;
  END;
  RETURN new;
END Gauss;

PROCEDURE (a: MATRIX) Det * (): LONGREAL;
(*
 * require
 *   a.DimM() = a.DimN()
 *)
VAR
  new: MATRIX;
  i: LONGINT;
  result: LONGREAL;
BEGIN
  new := a.Gauss()(MATRIX);
  result := new.elements[0, 0];
  FOR i := 1 TO new.m-1 DO
    result := result * new.elements[i, i];
  END;
  RETURN result;
END Det;

PROCEDURE (a: MATRIX) Norm * (): LONGREAL;
BEGIN
  RETURN a.Det();
END Norm;

PROCEDURE (a: MATRIX) Compare * (b: BT.COMPAREABLE): LONGINT;
(*
 * require
 *   a.DimM() = a.DimN()
 *   b(MATRIX).DimM() = b(MATRIX).DimN()
 *)
VAR
  i, j: LONGINT;
  an, bn: LONGREAL;
BEGIN
  WITH b: MATRIX DO
    an := a.Norm(); bn := b.Norm();
    IF    an = bn THEN RETURN  0
    ELSIF an > bn THEN RETURN  1
    ELSE               RETURN -1 END;
  END;
END Compare;

END MATRIX.