DIM A1(2,2)
      A1() = 1, 3, 5, 2, 4, 7, 1, 1, 0
      PROCLUdecomposition(A1(), L1(), U1(), P1())
      PRINT "L1:" ' FNshowmatrix(L1())
      PRINT "U1:" ' FNshowmatrix(U1())
      PRINT "P1:" ' FNshowmatrix(P1())

      DIM A2(3,3)
      A2() = 11, 9, 24, 2, 1, 5, 2, 6, 3, 17, 18, 1, 2, 5, 7, 1
      PROCLUdecomposition(A2(), L2(), U2(), P2())
      PRINT "L2:" ' FNshowmatrix(L2())
      PRINT "U2:" ' FNshowmatrix(U2())
      PRINT "P2:" ' FNshowmatrix(P2())
      END

      DEF PROCLUdecomposition(a(), RETURN l(), RETURN u(), RETURN p())
      LOCAL i%, j%, k%, n%, s, b() : n% = DIM(a(),2)
      DIM l(n%,n%), u(n%,n%), b(n%,n%)
      PROCpivot(a(), p())
      b() = p() . a()
      FOR j% = 0 TO n%
        l(j%,j%) = 1
        FOR i% = 0 TO j%
          s = 0
          FOR k% = 0 TO i% : s += u(k%,j%) * l(i%,k%) : NEXT
          u(i%,j%) = b(i%,j%) - s
        NEXT
        FOR i% = j% TO n%
          s = 0
          FOR k% = 0 TO j% : s += u(k%,j%) * l(i%,k%) : NEXT
          IF i%<>j% l(i%,j%) = (b(i%,j%) - s) / u(j%,j%)
        NEXT
      NEXT j%
      ENDPROC

      DEF PROCpivot(a(), RETURN p())
      LOCAL i%, j%, m%, n%, r% : n% = DIM(a(),2)
      DIM p(n%,n%) : FOR i% = 0 TO n% : p(i%,i%) = 1 : NEXT
      FOR i% = 0 TO n%
        m% = a(i%,i%)
        r% = i%
        FOR j% = i% TO n%
          IF a(j%,i%) > m% m% = a(j%,i%) : r% = j%
        NEXT
        IF i%<>r% THEN
          FOR j% = 0 TO n% : SWAP p(i%,j%),p(r%,j%) : NEXT
        ENDIF
      NEXT i%
      ENDPROC

      DEF FNshowmatrix(a())
      LOCAL @%, i%, j%, a$
      @% = &102050A
      FOR i% = 0 TO DIM(a(),1)
        FOR j% = 0 TO DIM(a(),2)
          a$ += STR$(a(i%,j%)) + ", "
        NEXT
        a$ = LEFT$(LEFT$(a$)) + CHR$(13) + CHR$(10)
      NEXT i%
      = a$