summaryrefslogtreecommitdiff
path: root/system/std.zusatz/1.7.3/src/matrix
diff options
context:
space:
mode:
Diffstat (limited to 'system/std.zusatz/1.7.3/src/matrix')
-rw-r--r--system/std.zusatz/1.7.3/src/matrix470
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;