% linalg.s
%
% Purpose: Matrix/vector operations script
% Copyright: R.A. van Engelen, Leiden University, 1997
banner("Linear Algebra Package").
infix ++ is associative(grouping = left).
[X] ++ [Y] := [[X, Y]].
[X | ColX] ++ [Y | ColY] := [[X, Y] | ColX ++ ColY].
[X] ++ [[Y | RowY]] := [[X, Y | RowY]].
[X | ColX] ++ [[Y | RowY] | MatY] := [[X, Y | RowY] | ColX ++ MatY].
[[X | RowX]] ++ [[Y | RowY]] := [[X | RowX] // [Y | RowY]].
([[X | RowX] | MatX] ++ [[Y | RowY] | MatY] :=
[[X | RowX] // [Y | RowY] | MatX ++ MatY])
.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Create zero matrix
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
zero(Rank) := zero(Rank, Rank).
zero(RowRank, ColRank) := fill(RowRank, fill(ColRank, 0)).
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Create identity matrix
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
ident(1) := [[1]].
ident(N) := (iff N <= 0 then error("Cannot create identity matrix of rank " // sprint(N)) else
[[1 | fill(N - 1, 0)] | zero(N - 1, 1) ++ ident(N - 1)]
).
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Matrix transpose
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
prefix transpose is intrinsic.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Matrix-matrix and vector-vector addition
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[[_G2509 | _G2510] | _G2513] + [_G2515 | _G2516] := error("matrix + vector").
[_G2509 | _G2510] + [[_G2512 | _G2513] | _G2516] := error("vector + matrix").
[[_G2509 | _G2510] | _G2513] - [_G2515 | _G2516] := error("matrix - vector").
[_G2509 | _G2510] - [[_G2512 | _G2513] | _G2516] := error("vector - matrix").
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Dot/inner product
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
infix .* has precedence(*) has latex_name("\cdot").
[X] .* [Y] := X * Y.
[X | Xs] .* [Y | Ys] := X * Y + Xs .* Ys.
_G2509 .* _G2510 := error("dot product of scalars").
[ _G2509
| _G2510 ] .* _G2513 := error("dot product of vector with scalar").
_G2512 .* [ _G2509
| _G2510 ] := error("dot product of scalar with vector").
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Cross/outer product
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
infix #* has precedence(*) has latex_name("\times").
[X] #* [Y] := [X * Y].
[X1, X2] #* [Y1, Y2] := [X1 * Y2 - X2 * Y1, X2 * Y1 - X1 * Y2].
([X1, X2, X3] #* [Y1, Y2, Y3] :=
[X2 * Y3 - X3 * Y2, X3 * Y1 - X1 * Y3, X1 * Y2 - X2 * Y1])
.
_G2509 #* _G2510 := error("cross product of scalars").
[ _G2509
| _G2510 ] #* _G2513 := error("cross product of vector with scalar").
_G2512 #* [ _G2509
| _G2510 ] := error("cross product of scalar with vector").
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Matrix-matrix and matrix-vector product
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
infix &* is associative(grouping = left) has latex_name("\,").
[X | Col] &* [Y] := [X | Col] * Y.
[X | Col] &* [[Y | Row]] := matmult_([X | Col], [Y | Row]).
[[X | Row]] &* [Y | Col] := [X | Row] .* [Y | Col].
([[X | Row1] | Rows1] &* [Y | Row2] :=
matmult_([[X | Row1] | Rows1], [[Y | Row2]]))
.
([[X | Row1] | Rows1] &* [[Y | Row2] | Rows2] :=
matmult_([[X | Row1] | Rows1], transpose [[Y | Row2] | Rows2]))
.
[_G2509 | _G2510] &* _G2513 := error("matrix product with scalar").
[ [ _G2509 | _G2510 ]
| _G2513 ] &* _G2516 := error("matrix product with scalar").
_G2512 &* [_G2509 | _G2510] := error("matrix product with scalar").
_G2515 &* [ [ _G2509 | _G2510 ]
| _G2513 ] := error("matrix product with scalar").
% matmult_: matrix product with transpose of matrix
(matmult_([[X | Row] | Rows], [[Y | Col]]) :=
[[X | Row] .* [Y | Col] | matmult_(Rows, [[Y | Col]])])
.
(matmult_([[X | Row] | Rows], [[Y | Col] | Cols]) :=
[ dotprodrows_([X | Row], [[Y | Col] | Cols])
| matmult_(Rows, [[Y | Col] | Cols]) ])
.
matmult_([], [[_G2509 | _G2510] | _G2513]) := [].
(matmult_([[X] | Rows], [Y | Col]) :=
[X * [Y | Col] | matmult_(Rows, [Y | Col])])
.
(matmult_([X | Rows], [Y | Col]) :=
[X * [Y | Col] | matmult_(Rows, [Y | Col])])
.
matmult_([], [_G2509 | _G2510]) := [].
matmult_(_G2509, _G2510) := error("incompatible matrix dimensions").
(dotprodrows_([X | Row1], [[Y | Row2] | Rows2]) :=
[[X | Row1] .* [Y | Row2] | dotprodrows_([X | Row1], Rows2)])
.
dotprodrows_([_G2509 | _G2510], []) := [].
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Gaussian elimination
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
gauss(_G2509) is intrinsic(module = gauss).
% Prolog predicate gauss/2
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Matrix power
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
([[X | Row] | Mat] ^ -1 :=
(map(drop(Rank, gauss([[X | Row] | Mat] ++ ident(Rank)))) where
Rank = length([[X | Row] | Mat])
))
.
[[X | Row] | _G2513] ^ 0 := ident(length([X | Row])).
[[X | Row] | Mat] ^ 1 := [[X | Row] | Mat].
[[X | Row] | Mat] ^ 2 := [[X | Row] | Mat] &* [[X | Row] | Mat].
([[X | Row] | Mat] ^ N :=
(iff op(0, N) == "T" then transpose [[X | Row] | Mat] else
iff N mod 2 == 0 then ( [[X | Row] | Mat]
^ (N div 2)
) ^ 2 else
[[X | Row] | Mat] &* [[X | Row] | Mat] ^ (N - 1)
))
.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Solve (MATLAB left division)
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[[X | Row] | Mat] \ [Y | Vec] := [[X | Row] | Mat] ^ -1 &* [Y | Vec].