Matrix.xtl

1⍝# Matrix: matrices -- identity, diagonal, trace, product,
2⍝# determinant, inverse, solving linear systems.
3⍝# Import it with an alias of your choice: "mx:" u_se< "Matrix".
4⍝# Put libs/Matrix/src on XETAL_PATH ("just path"); the reference is libs/Matrix/docs.
5⍝# Names with l: are exported; those under h: are private to this file.
6⍝#
7⍝# A matrix is a rank-2 array, rows first. The numerical functions
8⍝# (d_et, i_nverse, s_olve) work in Floats, by Gaussian elimination with
9⍝# partial pivoting; compare their results with Check's k:n_ear. The
10⍝# transpose is the built-in o_\.
11
12⍝## Basics
13
14⍝# i_dentity n: the n by n identity matrix (1 then n zeros, cycled).
15ˡi̲dentity ← { n → (n c̲at n) r̲eshape 1 c̲at n r̲eshape 0 }
16
17⍝# d_iag m: the items on the main diagonal.
18ˡd̲iag ← { m →
19  k ← 'm̲in r̲/ s̲hape m
20  (1 + (o̲ffsets k) × 1 + f̲irst r̲ev s̲hape m) s̲elect r̲avel m
21}
22
23⍝# t_race m: the sum of the diagonal.
24ˡt̲race ← { m → '+ r̲/ ˡd̲iag m }
25
26⍝# a m_ul b: the matrix product (a matrix and a vector work too).
27ˡm̲ul ← { a b → a '+ '× i̲nner b }
28
29⍝## Solving
30
31⍝# The permutation of 1..n that swaps i and j.
32ʰs̲wapping ← { n i j →
33  r ← r̲ange n
34  r + ((r = i) × j − i) + (r = j) × i − j
35}
36
37⍝# The row (k..n) whose item in column k is largest in size.
38ʰp̲ivot ← { k m →
39  c ← a̲bs (k − 1) d̲rop k s̲elect₂ m
40  k + (f̲irst g̲rade n̲eg c) − 1
41}
42
43⍝# Gauss-Jordan on m (Float) for columns k..n of its first n: each
44⍝# column's pivot row moved up, scaled to 1, and the column cleared
45⍝# from every other row.
46ʰg̲j ← { n k m →
47  k > n ? m
48  m ← ((n ʰs̲wapping k)_ k ʰp̲ivot m) s̲elect m
49  r ← r̲ange t̲ally m
50  v ← f̲irst k s̲elect k s̲elect₂ m              ⍝ the pivot
51  0.000000000001 > a̲bs v ? @ p̲anic< "the matrix is singular: it has no inverse, and its system no unique solution"
p̲anic< expands to(⎕P̲ANIC ("the matrix is singular: it has no inverse, and its system no unique solution"))
52  row ← (k s̲elect m) ÷ v
53  c ← (k s̲elect₂ m) − f̲loat r = k
54  ((n ʰg̲j k + 1)_ m − c '× t̲able row)
55}
56
57⍝# b s_olve a: x with a m_ul x = b, for a square a; b is a vector or a
58⍝# matrix of right-hand sides (one per column). A singular a divides by
59⍝# zero (error[division-by-zero]).
60ˡs̲olve ← { b a →
61  n ← t̲ally a
62  x ← n d̲rop₂ ((n ʰg̲j 1)_ f̲loat a c̲at₂ b)
63  (s̲hape b) r̲eshape x
64}
65
66⍝# i_nverse a: the inverse of a square matrix.
67ˡi̲nverse ← { a → (ˡi̲dentity t̲ally a) ˡs̲olve a }
68
69⍝## Determinant
70
71⍝# The determinant of a Float matrix m, by elimination: the first pivot
72⍝# times the determinant of what is left, negated for a row swap.
73ʰd̲etF ← { m →
74  n ← t̲ally m
75  n = 1 ? f̲irst r̲avel m
76  p ← 1 ʰp̲ivot m
77  m ← ((n ʰs̲wapping 1)_ p) s̲elect m
78  v ← f̲irst r̲avel m
79  v = 0 ? 0.0
80  s ← f̲loat (p = 1) − p ≠ 1
81  r ← (1 d̲rop₂ 1 d̲rop m) − ((1 d̲rop f̲irst₂ m) ÷ v) '× t̲able 1 d̲rop f̲irst m
82  s × v × ʰd̲etF r
83}
84
85⍝# d_et a: the determinant of a square matrix.
86ˡd̲et ← { a → ʰd̲etF f̲loat a }