diff --git a/src/BandedMatrices.jl b/src/BandedMatrices.jl index 0def4240..3e9f927e 100644 --- a/src/BandedMatrices.jl +++ b/src/BandedMatrices.jl @@ -7,8 +7,8 @@ import LinearAlgebra: BlasInt, BlasReal, BlasFloat, BlasComplex, axpy!, import LinearAlgebra.BLAS: libblas import LinearAlgebra.LAPACK: liblapack, chkuplo, chktrans import LinearAlgebra: cholesky, cholesky!, norm, diag, eigvals!, eigvals, eigen!, eigen, - qr, axpy!, ldiv!, mul!, lu, lu!, AbstractTriangular, has_offset_axes, - chkstride1, kron, lmul!, rmul!, factorize, StructuredMatrixStyle + qr, axpy!, ldiv!, mul!, lu, lu!, ldlt, ldlt!, AbstractTriangular, has_offset_axes, + chkstride1, kron, lmul!, rmul!, factorize, StructuredMatrixStyle, logabsdet import SparseArrays: sparse import Base: getindex, setindex!, *, +, -, ==, <, <=, >, @@ -74,6 +74,7 @@ include("banded/gbmm.jl") include("banded/linalg.jl") include("symbanded/symbanded.jl") +include("symbanded/ldlt.jl") include("symbanded/BandedCholesky.jl") include("tribanded.jl") diff --git a/src/symbanded/ldlt.jl b/src/symbanded/ldlt.jl new file mode 100644 index 00000000..092cae63 --- /dev/null +++ b/src/symbanded/ldlt.jl @@ -0,0 +1,75 @@ +## Banded LDLᵀ decomposition + +@inline bandwidth(A::LDLt{T,Symmetric{T,M}}, k) where {T,M<:BandedMatrix{T}} = bandwidth(A.data, k) + +factorize(S::Symmetric{<:Any,<:BandedMatrix}) = ldlt(S) + +function ldlt!(F::Symmetric{T,M}) where {T<:Real,M<:BandedMatrix{T}} + A = F.data + n = size(A, 1) + b = bandwidth(F) + if F.uplo == 'L' + @inbounds for j = 1:n + @simd for k = max(1,j-b-1):j-1 + A[j,j] -= abs2(A[j,k])*A[k,k] + end + for i = j+1:min(n,j+b) + @simd for k = max(1,j-b-1):j-1 + A[i,j] -= A[i,k]*A[j,k]*A[k,k] + end + A[i,j] /= A[j,j] + end + end + elseif F.uplo == 'U' + @inbounds for j = 1:n + @simd for k = max(1,j-b-1):j-1 + A[j,j] -= abs2(A[k,j])*A[k,k] + end + for i = j+1:min(n,j+b) + @simd for k = max(1,j-b-1):j-1 + A[j,i] -= A[k,i]*A[k,j]*A[k,k] + end + A[j,i] /= A[j,j] + end + end + end + return LDLt{T,Symmetric{T,M}}(F) +end + +function ldlt(A::Symmetric{T,M}) where {T<:Real,M<:BandedMatrix{T}} + S = typeof(zero(T)/one(T)) + uplo = A.uplo == 'U' ? :U : :L + return S == T ? ldlt!(copy(A)) : ldlt!(Symmetric(BandedMatrix{S}(A), uplo)) +end + + +function ldiv!(S::LDLt{T,Symmetric{T,M}}, B::AbstractVecOrMat{T}) where {T,M<:BandedMatrix{T}} + @assert !has_offset_axes(B) + n, nrhs = size(B, 1), size(B, 2) + if size(S, 1) != n + throw(DimensionMismatch("Matrix has dimensions $(size(S)) but right hand side has first dimension $n")) + end + A = S.data + b = bandwidth(A, 1) + @inbounds for l = 1:nrhs + for j = 1:n + @simd for k = max(1,j-b-1):j-1 + B[j,l] -= A[j,k]*B[k,l] + end + end + @simd for j = 1:n + B[j,l] /= A[j,j] + end + for j = n:-1:1 + @simd for k = j+1:min(n,j+b) + B[j,l] -= A[k,j]*B[k,l] + end + end + end + B +end + +function logabsdet(F::LDLt{T,Symmetric{T,M}}) where {T,M<:BandedMatrix{T}} + it = (F.data[i,i] for i in 1:size(F, 1)) + return sum(log∘abs, it), prod(sign, it) +end diff --git a/test/test_symbanded.jl b/test/test_symbanded.jl index 02e00d31..875507f8 100644 --- a/test/test_symbanded.jl +++ b/test/test_symbanded.jl @@ -51,7 +51,7 @@ using BandedMatrices, LinearAlgebra, LazyArrays, Random, Test end function Bn(::Type{T}, N::Int) where {T} - B = Symmetric(BandedMatrix(Zeros{T}(N,N), (0,2))) + B = Symmetric(BandedMatrix(Zeros{T}(N,N), (0, 2))) for n = 0:N-1 parent(B).data[3,n+1] = T(2*(n+1)*(n+2))/T((2n+1)*(2n+5)) end @@ -119,6 +119,40 @@ end BLAS.hbmv!('L', 1, one(T), view(parent(A).data,3:4,:), x, zero(T), similar(x))) end +@testset "LDLᵀ" begin + for T in (Float16, Float32, Float64, BigFloat, Rational{BigInt}) + A = BandedMatrix{T}(undef,(10,10),(0,2)) + A[band(0)] .= 4 + A[band(1)] .= -one(T)/4 + A[band(2)] .= -one(T)/16 + SAU = Symmetric(A, :U) + F = ldlt(SAU) + b = collect(one(T):size(F, 1)) + x = Matrix(SAU)\b + y = F\b + @test x ≈ y + A = BandedMatrix{T}(undef,(10,10),(2,0)) + A[band(0)] .= 4 + A[band(-1)] .= -one(T)/3 + A[band(-2)] .= -one(T)/9 + SAL = Symmetric(A, :L) + x = Matrix(SAL)\b + F = ldlt(SAL) + y = F\b + @test x ≈ y + @test_throws DimensionMismatch F\[b;b] + @test det(F) ≈ det(SAL) + end + for T in (Int16, Int32, Int64, BigInt) + A = BandedMatrix{T}(undef, (4,4), (1,1)) + A[band(0)] .= 3 + A[band(1)] .= 1 + F = ldlt(Symmetric(A)) + @test eltype(F) == float(T) + end +end + + @testset "Cholesky" begin A = Symmetric(BandedMatrix(0 => 1 ./ [12, 6, 6, 6, 12], 1 => ones(4) ./ 24))