1 files changed, 470 insertions, 0 deletions
diff --git a/system/std.zusatz/1.7.3/src/matrix b/system/std.zusatz/1.7.3/src/matrix
new file mode 100644
index 0000000..fbc5ffc
--- /dev/null
+++ b/system/std.zusatz/1.7.3/src/matrix
@@ -0,0 +1,470 @@
+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;