75 lines
1.5 KiB
Plaintext
75 lines
1.5 KiB
Plaintext
func is_square(m) { m.all { .len == m.len } }
|
||
func matrix_zero(n, m=n) { m.of { n.of(0) } }
|
||
func matrix_ident(n) { n.of {|i| n.of {|j| i==j ? 1 : 0 } } }
|
||
|
||
func pivotize(m) {
|
||
var size = m.len
|
||
var id = matrix_ident(size)
|
||
for i (^size) {
|
||
var max = m[i][i]
|
||
var row = i
|
||
for j (i .. size-1) {
|
||
if (m[j][i] > max) {
|
||
max = m[j][i]
|
||
row = j
|
||
}
|
||
}
|
||
if (row != i) {
|
||
id.swap(row, i)
|
||
}
|
||
}
|
||
return id
|
||
}
|
||
|
||
func mmult(a, b) {
|
||
var p = []
|
||
for r,c (^a ~X ^b[0]) {
|
||
for i (^b) {
|
||
p[r][c] := 0 += (a[r][i] * b[i][c])
|
||
}
|
||
}
|
||
return p
|
||
}
|
||
|
||
func lu(a) {
|
||
is_square(a) || die "Defined only for square matrices!";
|
||
var n = a.len
|
||
var P = pivotize(a)
|
||
var Aʼ = mmult(P, a)
|
||
var L = matrix_ident(n)
|
||
var U = matrix_zero(n)
|
||
for i,j (^n ~X ^n) {
|
||
if (j >= i) {
|
||
U[i][j] = (Aʼ[i][j] - ({ U[_][j] * L[i][_] }.map(^i).sum))
|
||
} else {
|
||
L[i][j] = (Aʼ[i][j] - ({ U[_][j] * L[i][_] }.map(^j).sum))/U[j][j]
|
||
}
|
||
}
|
||
return [P, Aʼ, L, U]
|
||
}
|
||
|
||
func say_it(message, array) {
|
||
say "\n#{message}"
|
||
array.each { |row|
|
||
say row.map{"%7s" % .as_rat}.join(' ')
|
||
}
|
||
}
|
||
|
||
var t = [[
|
||
%n(1 3 5),
|
||
%n(2 4 7),
|
||
%n(1 1 0),
|
||
],[
|
||
%n(11 9 24 2),
|
||
%n( 1 5 2 6),
|
||
%n( 3 17 18 1),
|
||
%n( 2 5 7 1),
|
||
]]
|
||
|
||
for test (t) {
|
||
say_it('A Matrix', test);
|
||
for a,b (['P Matrix', 'Aʼ Matrix', 'L Matrix', 'U Matrix'] ~Z lu(test)) {
|
||
say_it(a, b)
|
||
}
|
||
}
|