### Name: lu
### Title: Triangular Decomposition of a Square Matrix
### Aliases: lu lu,dgeMatrix-method
### Keywords: array algebra

### ** Examples


##--- Dense  -------------------------
x <- Matrix(rnorm(9), 3, 3)
lu(x)

##--- Sparse ------------------------

pm <- as(readMM(system.file("external/pores_1.mtx",
                            package = "Matrix")),
         "CsparseMatrix")
str(pmLU <- lu(pm))             # p is a 0-based permutation of the rows
                                # q is a 0-based permutation of the columns
## permute rows and columns of original matrix
ppm <- pm[pmLU@p + 1L, pmLU@q + 1L]
pLU <- pmLU@L %*% pmLU@U
## equal up to "rounding"
ppm[1:14, 1:5]
pLU[1:14, 1:5]  # product can have extra zeros
## "prove" consistency (up to rounding):
i0 <- ppm != pLU & ppm == 0
iN <- ppm != pLU & ppm != 0
stopifnot(all(abs((ppm - pLU)[i0]) < 1e-7), # absolute error for true 0
          all(abs((ppm - pLU)[iN]/ppm[iN]) < 1e-9)) # relative error