12345678910111213141516171819202122232425262728293031323334353637383940414243444546474849505152535455565758596061626364656667686970717273747576777879 |
- #!/usr/bin/ruby
- #
- ## https://rosettacode.org/wiki/LU_decomposition
- #
- 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| [i.of(0)..., 1, (n - i - 1).of(0)...] } }
- func pivotize(m) {
- var size = m.len
- var id = matrix_ident(size)
- for i in (^size) {
- var max = m[i][i]
- var row = i
- for j in (i ..^ size) {
- 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 in ^a, c in ^b[0], i in ^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 in ^n, j in ^n {
- if (j >= i) {
- U[i][j] = (Aʼ[i][j] - sum(^i, { U[_][j] * L[i][_] }))
- } else {
- L[i][j] = ((Aʼ[i][j] - sum(^j, { U[_][j] * L[i][_] })) / 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),
- ]]
- t.each { |test|
- say_it('A Matrix', test)
- for a,b in (['P Matrix', 'Aʼ Matrix', 'L Matrix', 'U Matrix'] ~Z lu(test)) {
- say_it(a, b)
- }
- }
|