PACKET matrix DEFINES MATRIX, matrix, idn,  (* Stand : 21.10.83      *)
                      :=, sub,              (* Autor : H.Indenbirken *)
                      row, column,
                      COLUMNS,
                      ROWS,
                      DET, 
                      INV, 
                      TRANSP, 
                      transp,
                      replace row, replace column,
                      replace element,
                      get, put,
                      =, <>,
                      +, -, * :

TYPE MATRIX = STRUCT (INT rows, columns, VECTOR elems);
TYPE INITMATRIX = STRUCT (INT rows, columns, REAL value, BOOL idn);

MATRIX VAR a :: idn (1);
INT VAR i;

(****************************************************************************
PROC dump (MATRIX CONST m) :
  put line (text (m.rows) + " Reihen, " + text (m.columns) + " Spalten.");
  dump (m.elems)  .

END PROC dump;
****************************************************************************)

OP := (MATRIX VAR l, MATRIX CONST r) :
  CONCR (l) := CONCR (r);
END OP :=;

OP := (MATRIX VAR l, INITMATRIX CONST r) :
  l.rows    := r.rows;
  l.columns := r.columns;
  l.elems   := vector (r.rows*r.columns, r.value);
  IF r.idn
  THEN idn FI  .

idn :
  INT VAR i;
  FOR i FROM 1 UPTO r.rows
  REP replace (l.elems, calc pos (l.columns, i, i), 1.0) PER

END OP :=;
 
INITMATRIX PROC matrix (INT CONST rows, columns, REAL CONST value) :
  IF rows <= 0
  THEN errorstop ("PROC matrix : rows <= 0")
  ELIF columns <= 0
  THEN errorstop ("PROC matrix : columns <= 0") FI;

  INITMATRIX : (rows, columns, value, FALSE)

END PROC matrix;

INITMATRIX PROC matrix (INT CONST rows, columns) :
  matrix (rows, columns, 0.0)

END PROC matrix;

INITMATRIX PROC idn (INT CONST size) :
  IF size <= 0
  THEN errorstop ("MATRIX PROC idn : size <= 0") FI;

  INITMATRIX : (size, size, 0.0, TRUE)

END PROC idn;

VECTOR PROC row (MATRIX CONST m, INT CONST i) :
  VECTOR VAR v :: vector (m.columns);
  INT VAR j, k :: 1, pos :: (i-1) * m.columns;
  FOR j FROM pos+1 UPTO pos + m.columns
  REP replace (v, k, m.elems SUB j);
      k INCR 1
  PER;
  v

END PROC row;

VECTOR PROC column (MATRIX CONST m, INT CONST j) :
  VECTOR VAR v :: vector (m.rows);
  INT VAR i, k :: j;
  FOR i FROM 1 UPTO m.rows
  REP replace (v, i, m.elems SUB k);
      k INCR m.columns
  PER;
  v

END PROC column;

INT OP COLUMNS (MATRIX CONST m) :
  m.columns

END OP COLUMNS;

INT OP ROWS (MATRIX CONST m) :
  m.rows

END OP ROWS;

REAL PROC sub (MATRIX CONST a, INT CONST row, column) :
  a.elems SUB calc pos (a.columns, row, column)

END PROC sub;

PROC replace row (MATRIX VAR m, INT CONST rowindex, VECTOR CONST rowvalue) :
  test ("PROC replace row : ", "LENGTH rowvalue", "COLUMNS m",
        LENGTH rowvalue, m.columns);
  test ("PROC replace row : row ", rowindex, m.rows);

  INT VAR i, pos :: (rowindex-1) * m.columns;
  FOR i FROM 1 UPTO m.columns
  REP replace (m.elems, pos+i, rowvalue SUB i) PER

END PROC replace row;

PROC replace column (MATRIX VAR m, INT CONST columnindex,
                     VECTOR CONST columnvalue) :
  test ("PROC replace column : ", "LENGTH columnvalue", "ROWS m",
        LENGTH columnvalue, m.rows);
  test ("PROC replace column : column ", columnindex, m.columns);

  INT VAR i;
  FOR i FROM 1 UPTO m.rows
  REP replace (m.elems, calc pos (m.columns, i, columnindex),
               columnvalue SUB i) PER

END PROC replace column;

PROC replace element (MATRIX VAR a, INT CONST row, column, REAL CONST x) :
  test ("PROC replace element : row ", row, a.rows);
  test ("PROC replace element : column ", column, a.columns);
  replace (a.elems, calc pos (a.columns, row, column), x)

END PROC replace element;

BOOL OP = (MATRIX CONST l, r) :
  IF l.rows <> r.rows
  THEN FALSE
  ELIF l.columns <> r.columns
  THEN FALSE
  ELSE l.elems = r.elems FI

END OP =;

BOOL OP <> (MATRIX CONST l, r) :
  IF l.rows <> r.rows
  THEN TRUE
  ELIF l.columns <> r.columns
  THEN TRUE
  ELSE l.elems <> r.elems FI

END OP <>;

INT PROC calc pos (INT CONST columns, z, s) :
 (z-1) * columns + s
END PROC calc pos;

MATRIX OP + (MATRIX CONST m) :
  m

END OP +;

MATRIX OP + (MATRIX CONST l, r) :
  test ("MATRIX OP + : ", "ROWS l", "ROWS r", l.rows, r.rows);
  test ("MATRIX OP + : ", "COLUMNS l", "COLUMNS r", l.columns, r.columns);

  a := l;
  INT VAR i;
  FOR i FROM 1 UPTO l.rows * l.columns
  REP replace (a.elems, i, (l.elems SUB i) + (r.elems SUB i))
  PER;
  a

END OP +;

MATRIX OP - (MATRIX CONST m) :
  a := m;
  INT VAR i;
  FOR i FROM 1 UPTO m.rows * m.columns
  REP replace (a.elems, i, -a.elems SUB i)
  PER;
  a

END OP -;

MATRIX OP - (MATRIX CONST l, r) :
  test ("MATRIX OP - : ", "ROWS l", "ROWS r", l.rows, r.rows);
  test ("MATRIX OP - : ", "COLUMNS l", "COLUMNS r", l.columns, r.columns);

  a := l;
  INT VAR i;
  FOR i FROM 1 UPTO l.rows * l.columns
  REP replace (a.elems, i, (l.elems SUB i) - (r.elems SUB i))
  PER;
  a

END OP -;

MATRIX OP * (REAL CONST x, MATRIX CONST m) :
 m*x

END OP *;

MATRIX OP * (MATRIX CONST m, REAL CONST x) :
  a := m;
  INT VAR i;
  FOR i FROM 1 UPTO m.rows * m.columns
  REP replace (a.elems, i, x*m.elems SUB i) PER;
  a

END OP *;

VECTOR OP * (VECTOR CONST v, MATRIX CONST m) :
  test ("VECTOR OP * : ", "LENGTH v", "ROWS m", LENGTH v, m.rows);
  VECTOR VAR result :: vector (m.rows);
  INT VAR i;
  FOR i FROM 1 UPTO m.rows
  REP replace (result, i, v * column (m, i)) PER;
  result  .

END OP *;

VECTOR OP * (MATRIX CONST m, VECTOR CONST v) :
  test ("VECTOR OP * : ", "COLUMNS m", "LENGTH v", COLUMNS m, LENGTH v);
  VECTOR VAR result :: vector (m.columns);
  INT VAR i;
  FOR i FROM 1 UPTO m.columns
  REP replace (result, i, row (m, i) * v) PER;
  result  .

END OP *;

MATRIX OP * (MATRIX CONST l, r) :
  test ("MATRIX OP * : ","COLUMNS l","ROWS r", l.columns, r.rows);

  a.rows := l.rows;
  a.columns := r.columns;
  a.elems := vector (a.rows*a.columns)
  INT VAR i, j;
  FOR i FROM 1 UPTO a.rows
  REP FOR j FROM 1 UPTO a.columns
      REP VECTOR VAR rl :: row (l, i), cr :: column (r, j);
          replace (a.elems, calc pos (a.columns, i, j), rl * cr)
      PER
  PER;
  a  .

END OP *;

PROC get (MATRIX VAR a, INT CONST rows, columns) :
 
  a := matrix (rows,columns);
  INT VAR i, j;
  VECTOR VAR v;
  FOR i FROM 1 UPTO rows
  REP get (v, columns);
      store row
  PER  .

store row :
  FOR j FROM 1 UPTO a.columns
  REP replace (a.elems, calc pos (a.columns, i, j), v SUB j)
  PER  .

END PROC get;

PROC put (MATRIX CONST a, INT CONST length, fracs) :
  INT VAR i, j;
  FOR i FROM 1 UPTO a.rows
  REP FOR j FROM 1 UPTO a.columns
      REP put (text (sub (a, i, j), length, fracs)) PER;
      line (2);
  PER

END PROC put;

PROC put (MATRIX CONST a) :
  INT VAR i, j;
  FOR i FROM 1 UPTO a.rows
  REP FOR j FROM 1 UPTO a.columns
      REP TEXT CONST number :: "                 " + text (sub (a, i, j));
          put (subtext (number, LENGTH number - 15))
      PER;
      line (2);
  PER

END PROC put;

TEXT VAR error :: "";
PROC test (TEXT CONST proc, l text, r text, INT CONST left, right) :
  IF left <> right
  THEN error := proc;
       error CAT l text;
       error CAT " (";
       error CAT text (left);
       error CAT ") <> ";
       error CAT r text;
       error CAT " (";
       error CAT text (right);
       error CAT ")";
       errorstop (error)
  FI  .

END PROC test;

PROC test (TEXT CONST proc, INT CONST i, n) :
  IF i < 1
  THEN error := proc;
       error CAT "subscript underflow (";
       error CAT text (i);
       error CAT ")";
       errorstop (error)
  ELIF i > n
  THEN error := proc;
       error CAT "subscript overflow (i=";
       error CAT text (i);
       error CAT ", max=";
       IF n <= 0
       THEN error CAT "undefined"
       ELSE error CAT text (n) FI;
       error CAT ")";
       errorstop (error)
  FI

END PROC test;


MATRIX OP TRANSP (MATRIX CONST m) :
  MATRIX VAR a :: m;
  transp (a);
  a

END OP TRANSP;

PROC transp (MATRIX VAR m) :
  INT VAR k :: 1, n :: m.rows*m.columns;
  a := m;
  FOR i FROM 2 UPTO n
  REP replace (m.elems, i, a.elems SUB position) PER;
  a := idn (1);
  i := m.rows;
  m.rows := m.columns;
  m.columns := i  .

position :
  k INCR m.columns;
  IF k > n
  THEN k DECR (n-1) FI;
  k  .
END PROC transp;

MATRIX OP INV (MATRIX CONST m) :
  a := m;
  ROW 32 INT VAR pivots;
  INT VAR i, j, k :: ROWS a, n :: COLUMNS a, pos;

  IF n <> k
  THEN errorstop ("MATRIX OP INV : no square matrix") FI;

  initialisiere die pivotpositionen;

  FOR j FROM 1 UPTO n
  REP pivotsuche (a, j, pos);
      IF sub (a, pos, pos) = 0.0
      THEN errorstop ("MATRIX OP INV : singular matrix") FI;
      zeilentausch (a, j, pos);
      merke dir die vertauschung;
      transformiere die matrix
  PER;

  spaltentausch;
  a  .

initialisiere die pivotpositionen :
  FOR i FROM 1 UPTO n
  REP pivots [i] := i PER  .

merke dir die vertauschung :
  IF pos > j
  THEN INT VAR hi :: pivots [j];
       pivots [j] := pivots [pos];
       pivots [pos] := hi
  FI  .

transformiere die matrix :
  REAL VAR h := 1.0/sub (a, j, j);

  FOR k FROM 1 UPTO n
  REP IF k <> j
      THEN FOR i FROM 1 UPTO n
           REP IF i <> j
               THEN replace element (a, i, k, sub (a, i, k) -
                                     sub (a, i, j)*sub (a, j, k)*h);
               FI
           PER;
       FI
  PER;

  FOR k FROM 1 UPTO n
  REP replace element (a, j, k, -h*sub (a, j, k));
      replace element (a, k, j,  h*sub (a, k, j))
  PER;
  replace element (a, j, j, h)  .

spaltentausch :
  VECTOR VAR v :: vector (n);
  FOR i FROM 1 UPTO n
  REP FOR k FROM 1 UPTO n
      REP replace (v, pivots [k], sub (a, i, k)) PER;
      replace row (a, i, v)
  PER  .

END OP INV;

REAL OP DET (MATRIX CONST m) :
  IF COLUMNS m <> ROWS m
  THEN errorstop ("REAL OP DET : no square matrix") FI;

  a := m;
  INT VAR i, j, k, n :: COLUMNS m, pos;
  FOR j FROM 1 UPTO n
  REP pivotsuche (a, j, pos);
      zeilentausch (a, j, pos);
      transformiere die matrix
  PER;
  produkt der pivotelemente  .

transformiere die matrix :
  REAL VAR h := 1.0/sub (a, j, j);
  FOR i FROM j+1 UPTO n
  REP FOR k FROM j+1 UPTO n
      REP replace element (a, i, k, sub (a, i, k) -
                           sub (a, i, j)*h*sub (a, j, k))
      PER
  PER  .

produkt der pivotelemente :
  REAL VAR produkt :: sub (a, 1, 1);
  FOR j FROM 2 UPTO n
  REP produkt := produkt * sub (a, j, j) PER;
  a := idn (1);
  produkt  .

END OP DET;

PROC pivotsuche (MATRIX CONST a, INT CONST start pos, INT VAR pos) :
  REAL VAR max :: abs (sub (a, start pos, start pos));
  INT VAR i;
  pos := start pos;

  FOR i FROM start pos+1 UPTO COLUMNS a
  REP IF abs (sub (a, i, start pos)) > max
      THEN max := abs (sub (a, i, start pos));
           pos := i
      FI
  PER  .

END PROC pivotsuche;

PROC zeilentausch (MATRIX VAR a, INT CONST old pos, pos) :
  VECTOR VAR v := row (a, pos);
  replace row (a, pos, row (a, old pos));
  replace row (a, old pos, v)  .

END PROC zeilentausch;

END PACKET matrix;