[50G] Question

[Tugdual]

I have reimplemented the LU decomposition using Crout method as described on Wikipedia with one limit, the P matrix is not built (nor used) so it will crash in some cases.

Code:

%%HP: T(3)A(R)F(,);
\<< DUPDUP SIZE LIST\-> DROP2 \-> A L N
  \<< 'L' IDN 1 N 1 -
    FOR n n 1 + N
      FOR i 'A' i n 2 \->LIST GET 'A' n DUP 2 \->LIST GET / NEG
      NEXT N n - 1 2 \->LIST \->ARRY n DUP 1 + SWAP 2 \->LIST A IDN SWAP ROT REPL DUP A * 'A' STO L SWAP INV * 'L' STO
    NEXT L A
  \>>
\>>

Usage:
Have the matrix on 1: and call this function.
Returns 2 matrices, 2:L and 1:U

Comments
Returns the exact same result as the Prime CAS in the few samples I tried.
Much slower than native LU() method but acceptable for normal use.
Also the native LU() has a diagonal of 1 on U while common use is to have it on L which is the case with this implementation