LDLFactorizations: Factorization of Symmetric Matrices
A translation of Tim Davis's Concise LDLᵀ Factorization, part of SuiteSparse with several improvements.
Please cite this repository if you use LDLFactorizations.jl in your work: see
This package is appropriate for matrices A that possess a factorization of the form LDLᵀ without pivoting, where L is unit lower triangular and D is diagonal (indefinite in general), including definite and quasi-definite matrices.
LDLFactorizations.jl should not be expected to be as fast, as robust or as accurate as factorization packages such as HSL.jl, MUMPS.jl or Pardiso.jl. Those are multifrontal and/or implement various forms of parallelism, and employ sophisticated pivot strategies.
The main advantages of LDLFactorizations.jl are that
it is very short and has a small footprint;
it is in pure Julia, and so
2.a. it does not require external compiled dependencies;
2.b. it will work with multiple input data types.
Whereas MUMPS.jl, HSL.jl and Pardiso.jl only work with single and double precision reals and complex data types, LDLFactorizations.jl accepts any numerical data type.
julia> ] pkg> add LDLFactorizations
The only exported functions are
ldl() with a dense array converts it to a sparse matrix.
A permutation ordering can be supplied:
ldl(A, p) where
p is an
array representing a permutation of the integers between 1 and the order
If no permutation is supplied, one is automatically computed using AMD.jl.
Only the upper triangle of
A is accessed.
ldl returns a factorization in the form of a
ldiv! methods are implemented for objects of type
LDLFactorization so that
solving a linear system is as easy as
LDLT = ldl(A) # LDLᵀ factorization of A x = LDLT \ b # solves Ax = b ldiv!(LDLT, b) # computes LDLT \ b in-place and overwriting b to store the result y = similar(b) ldiv!(y, LDLT, b) # computes LDLT \ b in-place and store the result in y
The factorization can of course be reused to solve for multiple right-hand sides.
Factors can be accessed as
LDLT.D, and the permutation vector as
Because the L factor is unit lower triangular, its diagonal is not stored.
Thus the factors satisfy: PAPᵀ = (L + I) D (L + I)ᵀ.
Timothy A. Davis. 2005. Algorithm 849: A concise sparse Cholesky factorization package. ACM Trans. Math. Softw. 31, 4 (December 2005), 587-591. DOI 10.1145/1114268.1114277.
Like the original LDL, this package is distributed under the LGPL.