LU decomposition: Difference between revisions
jq
(added R) |
(jq) |
||
Line 1,349:
3 17 18 1
2 5 7 1</lang>
=={{header|jq}}==
{{Works with|jq|1.4}}
jq currently does not have builtin support for matrices and therefore
some infrastructure is needed to make the following self-contained.
Matrices here are represented as arrays of arrays in the usual way.
'''Infrastructure'''
<lang jq># Create an m x n matrix
def matrix(m; n; init):
if m == 0 then []
elif m == 1 then [range(0;n)] | map(init)
elif m > 0 then
matrix(1;n;init) as $row
| [range(0;m)] | map( $row )
else error("matrix\(m);_;_) invalid")
end ;
def I(n): matrix(n;n;0) as $m
| reduce range(0;n) as $i ($m; . | setpath( [$i,$i]; 1));
def dot_product(a; b):
reduce range(0;a|length) as $i (0; . + (a[$i] * b[$i]) );
# transpose/0 expects its input to be a rectangular matrix
def transpose:
if (.[0] | length) == 0 then []
else [map(.[0])] + (map(.[1:]) | transpose)
end ;
# A and B should both be numeric matrices, A being m by n, and B being n by p.
def multiply(A; B):
(B[0]|length) as $p
| (B|transpose) as $BT
| reduce range(0; A|length) as $i
([];
reduce range(0; $p) as $j
(.;
.[$i][$j] = dot_product( A[$i]; $BT[$j] ) ));
def swap_rows(i;j):
if i == j then .
else .[i] as $i | .[i] = .[j] | .[j] = $i
end ;
# Print a matrix neatly, each cell occupying n spaces, but without truncation
def neatly(n):
def right: tostring | ( " " * (n-length) + .);
. as $in
| length as $length
| reduce range (0;$length) as $i
(""; . + reduce range(0;$length) as $j
(""; "\(.) \($in[$i][$j] | right )" ) + "\n" ) ;</lang>
'''LU decomposition'''
<lang jq># Create the pivot matrix for the input matrix:
def pivotize:
length as $n
| . as $m
| reduce range(0;$n) as $j
(I($n);
# state: [row; max]
(reduce range($j+1; $n) as $i
([$j, $m[$j][$j] ];
if $m[$i][$j] > .[1] then [ $i, $m[$i][$j] ] else . end) | .[0]) as $row
| swap_rows( $j; $row)
) ;
# Decompose the input nxn matrix A by PA=LU and return [L, U, P].
def lup:
def div(i;j):
if j == 0 then if i==0 then 0 else error("\(i)/0") end
else i/j
end;
. as $A
| length as $n
| I($n) as $L # matrix($n; $n; 0.0) as $L
| matrix($n; $n; 0.0) as $U
| ($A|pivotize) as $P
| multiply($P;$A) as $A2
# state: [L, U]
| reduce range(0; $n) as $i ( [$L, $U];
reduce range(0; $n) as $j (.;
.[0] as $L
| .[1] as $U
| if ($j >= $i) then
(reduce range(0;$i) as $k (0; . + ($U[$k][$j] * $L[$i][$k] ))) as $s1
| [$L, ($U| setpath([$i,$j]; ($A2[$i][$j] - $s1))) ]
else
(reduce range(0;$j) as $k (0; . + ($U[$k][$j] * $L[$i][$k]))) as $s2
| [ ($L | setpath([$i,$j]; div(($A2[$i][$j] - $s2) ; $U[$j][$j] ))), $U ]
end ))
| . + [ $P ]
;
</lang>
'''Example 1''':
<lang jq>def a: [[1, 3, 5], [2, 4, 7], [1, 1, 0]];
a | lup[] | neatly(4)
</lang>
{{Out}}
$ /usr/local/bin/jq -M -n -r -f LU.jq
1 0 0
0.5 1 0
0.5 -1 1
2 4 7
0 1 1.5
0 0 -2
0 1 0
1 0 0
0 0 1
'''Example 2''':
<lang jq>def b: [[11,9,24,2],[1,5,2,6],[3,17,18,1],[2,5,7,1]];
b | lup[] | neatly(21)</lang>
{{Out}}
$ /usr/local/bin/jq -M -n -r -f LU.jq
1 0 0 0
0.2727272727272727 1 0 0
0.09090909090909091 0.2875 1 0
0.18181818181818182 0.23124999999999996 0.0035971223021580693 1
11 9 24 2
0 14.545454545454547 11.454545454545455 0.4545454545454546
0 0 -3.4749999999999996 5.6875
0 0 0 0.510791366906476
1 0 0 0
0 0 1 0
0 1 0 0
0 0 0 1
=={{header|Julia}}==
| |||