calcL =: (% {.)@:({."1)
calcU =: (}.@[ - {.@[ *"1 0 ])&:(}."1)
lumain =: ((({.@[ ,: ]) ,&.:(|."1) $:@calcU) calcL)^:(*@#)
lu =: [: (|:@] ,: [)/ 1 0 2 |: _2 ]\ lumain
NB. Half this code is handling joining ragged lists.
NB. Is there a better way?
saveAA =: 3 3 $ 2 1 4 _4 _1 _11 2 4 _2
lu saveAA
1 0 0
_2 1 0
1 3 1
2 1 4
0 1 _3
0 0 3
I suspect that a vectorized explicit version is a better way to go.
This version has memory requirements of O(n^3).
Henry Rich
On 1/12/2014 9:00 PM, km wrote:
Verb LU below produces the matrices L and U of the LU decomposition of a square
matrix A. L is lower triangular, U is upper triangular, and A is L +/ . * U .
Should one attempt a tacit version?
eye =: =@i.@] NB. eye 3 is a 3 by 3 identity matrix
rop =: 3 : 0 NB. row op: subtract c times row i0 from row i1
:
'i1 c i0' =. x
( (i1 { y) - c * i0 { y ) i1 } y
)
LU =: 3 : 0 NB. square matrices L and U for y -: L +/ . * U
m =. # y
L =. eye(m)
U =. y
for_j. i. <: m do.
p =. (< j , j) { U
for_i. j + >: i. <: m - j do.
c =. p %~ (< i , j) { U
L =. c (< i , j) } L
U =. (i, c, j) rop U
end.
end.
L ,: U
)
saveAA
2 1 4
_4 _1 _11
2 4 _2
LU saveAA
1 0 0
_2 1 0
1 3 1
2 1 4
0 1 _3
0 0 3
saveAA -: +/ . */ LU saveAA
1
--Kip Murray
Sent from my iPad
----------------------------------------------------------------------
For information about J forums see http://www.jsoftware.com/forums.htm
----------------------------------------------------------------------
For information about J forums see http://www.jsoftware.com/forums.htm
