lu: Triangular Decomposition of a Square Matrix

Description

Computes triangular decompositions of square matrices.

Usage

lu(x, ...)
## S3 method for class 'dgeMatrix':
lu(x, warnSing = TRUE, \dots)

Arguments

a matrix of square dimension. No missing values or IEEE special values are allowed.

warnSing

(when x is a "denseMatrix") logical specifying if a warning should be signalled when x is singular.

...

further arguments passed to or from other methods.

Value

An object of class "LU", i.e., "denseLU" or "sparseLU", see sparseLU; this is a representation of a triangular decomposition of x.

Details

lu() is a generic function with special methods for different types of matrices. Use showMethods("lu") to list all the methods for the lu generic.

The method for class dgeMatrix (and all dense matrices) is based on LAPACK's "dgetrf" subroutine. It returns a decomposition also for singular matrices.

The method for class dgCMatrix (and all sparse matrices) is based on functions from the CSparse library. It signals an error when the decomposition algorithm fails, as when x is (too close to) singular.

References

Golub, G., and Van Loan, C. F. (1989). Matrix Computations, 2nd edition, Johns Hopkins, Baltimore.

Tim Davis (2005) http://www.cise.ufl.edu/research/sparse/CSparse/

Timothy A. Davis (2006) Direct Methods for Sparse Linear Systems, SIAM Series Fundamentals of Algorithms.

Examples

Run this code

##--- 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

Run the code above in your browser using DataLab

Last chance! 50% off unlimited learning