I think I prefer LU to be explicit.
But I think I'd write rop tacitly. Perhaps:
rop=: (({~ {.) - ({~ {:) * 1{])~`({.@[)`]}
Or, just for fun:
rop=: ] -"1 ] +/ .*~ 1&{@[ * 1 0 1&#@[ */ .="1(#: i.)@$@]
Thanks,
--
Raul
On Sun, Jan 12, 2014 at 9:00 PM, km <[email protected]> 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
