lu.RdComputes the pivoted LU factorization of an \(m \times n\) real matrix \(A\), which has the general form $$P_{1} A P_{2} = L U$$ or (equivalently) $$A = P_{1}' L U P_{2}'$$ where \(P_{1}\) is an \(m \times m\) permutation matrix, \(P_{2}\) is an \(n \times n\) permutation matrix, \(L\) is an \(m \times \min(m,n)\) unit lower trapezoidal matrix, and \(U\) is a \(\min(m,n) \times n\) upper trapezoidal matrix.
Methods for denseMatrix are built on
LAPACK routine dgetrf, which does not permute columns,
so that \(P_{2}\) is an identity matrix.
Methods for sparseMatrix are built on
CSparse routine cs_lu, which requires \(m = n\),
so that \(L\) and \(U\) are triangular matrices.
lu(x, ...)
# S4 method for class 'dgeMatrix'
lu(x, warnSing = TRUE, ...)
# S4 method for class 'dgCMatrix'
lu(x, errSing = TRUE, order = NA_integer_,
tol = 1, ...)
# S4 method for class 'dsyMatrix'
lu(x, cache = TRUE, ...)
# S4 method for class 'dsCMatrix'
lu(x, cache = TRUE, ...)
# S4 method for class 'matrix'
lu(x, ...)a finite matrix or
Matrix to be factorized,
which must be square if sparse.
a logical indicating if a warning should
be signaled for singular x. Used only by methods for
dense matrices.
a logical indicating if an error should
be signaled for singular x. Used only by methods for
sparse matrices.
an integer in 0:3 passed to CSparse routine
cs_sqr, indicating a strategy for choosing the column
permutation \(P_{2}\). 0 means no column permutation.
1, 2, and 3 indicate a fill-reducing ordering of \(A + A'\),
\(\tilde{A}' \tilde{A}\), and \(A' A\),
where \(\tilde{A}\) is \(A\) with “dense” rows
removed.
NA (the default) is equivalent to 2 if tol == 1
and 1 otherwise.
Do not set to 0 unless you know that the column order of \(A\)
is already sensible.
a number. The original pivot element is used
if its absolute value exceeds tol * a,
where a is the maximum in absolute value of the
other possible pivot elements.
Set tol < 1 only if you know what you are doing.
a logical indicating if the result should be
cached in x@factors[["LU"]]. Note that
caching is experimental and that only methods for classes
extending compMatrix will have this
argument.
further arguments passed to or from methods.
An object representing the factorization, inheriting from
virtual class LU. The specific class
is denseLU unless x inherits
from virtual class sparseMatrix,
in which case it is sparseLU.
What happens when x is determined to be near-singular
differs by method. The method for class dgeMatrix
completes the factorization, warning if warnSing = TRUE
and in any case returning a valid denseLU
object. Users of this method can detect singular x with
a suitable warning handler; see tryCatch.
In contrast, the method for class dgCMatrix
abandons further computation, throwing an error if errSing = TRUE
and otherwise returning NA. Users of this method can
detect singular x with an error handler or by setting
errSing = FALSE and testing for a formal result with
is(., "sparseLU").
The LAPACK source code, including documentation; see https://netlib.org/lapack/double/dgetrf.f.
Davis, T. A. (2006). Direct methods for sparse linear systems. Society for Industrial and Applied Mathematics. doi:10.1137/1.9780898718881
Golub, G. H., & Van Loan, C. F. (2013). Matrix computations (4th ed.). Johns Hopkins University Press. doi:10.56021/9781421407944
showMethods("lu", inherited = FALSE)
#> Function: lu (package Matrix)
#> x="denseMatrix"
#> x="dgCMatrix"
#> x="dgRMatrix"
#> x="dgTMatrix"
#> x="dgeMatrix"
#> x="diagonalMatrix"
#> x="dsCMatrix"
#> x="dsRMatrix"
#> x="dsTMatrix"
#> x="dspMatrix"
#> x="dsyMatrix"
#> x="dtCMatrix"
#> x="dtRMatrix"
#> x="dtTMatrix"
#> x="dtpMatrix"
#> x="dtrMatrix"
#> x="matrix"
#> x="sparseMatrix"
#>
set.seed(0)
## ---- Dense ----------------------------------------------------------
(A1 <- Matrix(rnorm(9L), 3L, 3L))
#> 3 x 3 Matrix of class "dgeMatrix"
#> [,1] [,2] [,3]
#> [1,] 1.2629543 1.2724293 -0.928567035
#> [2,] -0.3262334 0.4146414 -0.294720447
#> [3,] 1.3297993 -1.5399500 -0.005767173
(lu.A1 <- lu(A1))
#> LU factorization of Formal class 'denseLU' [package "Matrix"] with 4 slots
#> ..@ x : num [1:9] 1.33 0.95 -0.245 -1.54 2.735 ...
#> ..@ perm : int [1:3] 3 3 3
#> ..@ Dim : int [1:2] 3 3
#> ..@ Dimnames:List of 2
#> .. ..$ : NULL
#> .. ..$ : NULL
(A2 <- round(10 * A1[, -3L]))
#> 3 x 2 Matrix of class "dgeMatrix"
#> [,1] [,2]
#> [1,] 13 13
#> [2,] -3 4
#> [3,] 13 -15
(lu.A2 <- lu(A2))
#> LU factorization of Formal class 'denseLU' [package "Matrix"] with 4 slots
#> ..@ x : num [1:6] 13 1 -0.231 13 -28 ...
#> ..@ perm : int [1:2] 1 3
#> ..@ Dim : int [1:2] 3 2
#> ..@ Dimnames:List of 2
#> .. ..$ : NULL
#> .. ..$ : NULL
## A ~ P1' L U in floating point
str(e.lu.A2 <- expand2(lu.A2), max.level = 2L)
#> List of 3
#> $ P1.:Formal class 'pMatrix' [package "Matrix"] with 5 slots
#> $ L :Formal class 'dgeMatrix' [package "Matrix"] with 4 slots
#> $ U :Formal class 'dtrMatrix' [package "Matrix"] with 5 slots
stopifnot(all.equal(A2, Reduce(`%*%`, e.lu.A2)))
## ---- Sparse ---------------------------------------------------------
A3 <- as(readMM(system.file("external/pores_1.mtx", package = "Matrix")),
"CsparseMatrix")
(lu.A3 <- lu(A3))
#> LU factorization of Formal class 'sparseLU' [package "Matrix"] with 6 slots
#> ..@ L :Formal class 'dtCMatrix' [package "Matrix"] with 7 slots
#> .. .. ..@ i : int [1:141] 0 1 18 22 25 26 1 18 22 25 ...
#> .. .. ..@ p : int [1:31] 0 6 11 13 16 19 22 26 31 35 ...
#> .. .. ..@ Dim : int [1:2] 30 30
#> .. .. ..@ Dimnames:List of 2
#> .. .. .. ..$ : NULL
#> .. .. .. ..$ : NULL
#> .. .. ..@ x : num [1:141] 1 -0.993819 -0.004979 -0.000132 0.000132 ...
#> .. .. ..@ uplo : chr "L"
#> .. .. ..@ diag : chr "N"
#> ..@ U :Formal class 'dtCMatrix' [package "Matrix"] with 7 slots
#> .. .. ..@ i : int [1:206] 0 0 1 2 2 3 2 4 2 3 ...
#> .. .. ..@ p : int [1:31] 0 1 3 4 6 8 12 15 20 25 ...
#> .. .. ..@ Dim : int [1:2] 30 30
#> .. .. ..@ Dimnames:List of 2
#> .. .. .. ..$ : NULL
#> .. .. .. ..$ : NULL
#> .. .. ..@ x : num [1:206] -7178502 -24613411 -18311731 -6399179 715 ...
#> .. .. ..@ uplo : chr "U"
#> .. .. ..@ diag : chr "N"
#> ..@ p : int [1:30] 1 11 29 27 19 28 9 18 17 7 ...
#> ..@ q : int [1:30] 0 1 29 27 19 28 9 18 17 8 ...
#> ..@ Dim : int [1:2] 30 30
#> ..@ Dimnames:List of 2
#> .. ..$ : NULL
#> .. ..$ : NULL
## A ~ P1' L U P2' in floating point
str(e.lu.A3 <- expand2(lu.A3), max.level = 2L)
#> List of 4
#> $ P1.:Formal class 'pMatrix' [package "Matrix"] with 5 slots
#> $ L :Formal class 'dtCMatrix' [package "Matrix"] with 7 slots
#> $ U :Formal class 'dtCMatrix' [package "Matrix"] with 7 slots
#> $ P2.:Formal class 'pMatrix' [package "Matrix"] with 5 slots
stopifnot(all.equal(A3, Reduce(`%*%`, e.lu.A3)))