From a9dc0d5b12bee6b9d644acc6f80625fc23f286b3 Mon Sep 17 00:00:00 2001 From: Mikael Slevinsky Date: Mon, 11 Feb 2019 15:16:52 -0600 Subject: [PATCH 1/4] add ldlt --- src/BandedMatrices.jl | 5 ++-- src/symbanded/ldlt.jl | 57 ++++++++++++++++++++++++++++++++++++++++++ test/test_symbanded.jl | 25 ++++++++++++++++++ 3 files changed, 85 insertions(+), 2 deletions(-) create mode 100644 src/symbanded/ldlt.jl diff --git a/src/BandedMatrices.jl b/src/BandedMatrices.jl index 0def4240..5f19892f 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, 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..ce9ea84f --- /dev/null +++ b/src/symbanded/ldlt.jl @@ -0,0 +1,57 @@ +## 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(A::Symmetric{T,M}) where {T<:Real,M<:BandedMatrix{T}} + F = Symmetric(similar(A.data), :L) + B = F.data + n = size(A, 1) + b = bandwidth(A, 1) + for j = 1:n + B[j,j] = A[j,j] + for k = max(1,j-b-1):j-1 + B[j,j] -= B[j,k]^2*B[k,k] + end + for i = j+1:min(n,j+b) + B[i,j] = A[i,j] + for k = max(1,j-b-1):j-1 + B[i,j] -= B[i,k]*B[j,k]*B[k,k] + end + B[i,j] /= B[j,j] + end + end + return LDLt{T,Symmetric{T,M}}(F) +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..2e07fabe 100644 --- a/test/test_symbanded.jl +++ b/test/test_symbanded.jl @@ -119,6 +119,31 @@ 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 subtypes(AbstractFloat) + A = BandedMatrix{T}(undef,(10,10),(2,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[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 +end + + @testset "Cholesky" begin A = Symmetric(BandedMatrix(0 => 1 ./ [12, 6, 6, 6, 12], 1 => ones(4) ./ 24)) From d3b592a086c0971a1b9b409fbd858ecf87e4dd1b Mon Sep 17 00:00:00 2001 From: Mikael Slevinsky Date: Mon, 11 Feb 2019 15:52:18 -0600 Subject: [PATCH 2/4] make it in-place --- src/BandedMatrices.jl | 2 +- src/symbanded/ldlt.jl | 41 +++++++++++++++++++++++++++++------------ 2 files changed, 30 insertions(+), 13 deletions(-) diff --git a/src/BandedMatrices.jl b/src/BandedMatrices.jl index 5f19892f..3e9f927e 100644 --- a/src/BandedMatrices.jl +++ b/src/BandedMatrices.jl @@ -7,7 +7,7 @@ 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!, ldlt, AbstractTriangular, has_offset_axes, + qr, axpy!, ldiv!, mul!, lu, lu!, ldlt, ldlt!, AbstractTriangular, has_offset_axes, chkstride1, kron, lmul!, rmul!, factorize, StructuredMatrixStyle, logabsdet import SparseArrays: sparse diff --git a/src/symbanded/ldlt.jl b/src/symbanded/ldlt.jl index ce9ea84f..7e4618da 100644 --- a/src/symbanded/ldlt.jl +++ b/src/symbanded/ldlt.jl @@ -4,27 +4,44 @@ factorize(S::Symmetric{<:Any,<:BandedMatrix}) = ldlt(S) -function ldlt(A::Symmetric{T,M}) where {T<:Real,M<:BandedMatrix{T}} - F = Symmetric(similar(A.data), :L) - B = F.data +function ldlt!(F::Symmetric{T,M}) where {T<:Real,M<:BandedMatrix{T}} + A = F.data n = size(A, 1) b = bandwidth(A, 1) - for j = 1:n - B[j,j] = A[j,j] - for k = max(1,j-b-1):j-1 - B[j,j] -= B[j,k]^2*B[k,k] + 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 - for i = j+1:min(n,j+b) - B[i,j] = A[i,j] - for k = max(1,j-b-1):j-1 - B[i,j] -= B[i,k]*B[j,k]*B[k,k] + 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 - B[i,j] /= B[j,j] 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)) + return S == T ? ldlt!(copy(A)) : ldlt!(Symmetric(BandedMatrix{S}(M), 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) From 3b786e45ea29ca43ed43fb843f827374bd83f5fb Mon Sep 17 00:00:00 2001 From: Mikael Slevinsky Date: Mon, 11 Feb 2019 15:59:43 -0600 Subject: [PATCH 3/4] `subtypes` has been moved out of base --- test/test_symbanded.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/test_symbanded.jl b/test/test_symbanded.jl index 2e07fabe..54b562e5 100644 --- a/test/test_symbanded.jl +++ b/test/test_symbanded.jl @@ -120,7 +120,7 @@ end end @testset "LDLᵀ" begin - for T in subtypes(AbstractFloat) + for T in (Float16, Float32, Float64, BigFloat) A = BandedMatrix{T}(undef,(10,10),(2,2)) A[band(0)] .= 4 A[band(1)] .= -one(T)/4 From 0b71eeb3e6016b4b7227966bca8da1349ff260b9 Mon Sep 17 00:00:00 2001 From: Mikael Slevinsky Date: Tue, 12 Feb 2019 11:04:53 -0600 Subject: [PATCH 4/4] fix two errors add bandwidth tests add type promotion tests --- src/symbanded/ldlt.jl | 5 +++-- test/test_symbanded.jl | 15 ++++++++++++--- 2 files changed, 15 insertions(+), 5 deletions(-) diff --git a/src/symbanded/ldlt.jl b/src/symbanded/ldlt.jl index 7e4618da..092cae63 100644 --- a/src/symbanded/ldlt.jl +++ b/src/symbanded/ldlt.jl @@ -7,7 +7,7 @@ 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(A, 1) + b = bandwidth(F) if F.uplo == 'L' @inbounds for j = 1:n @simd for k = max(1,j-b-1):j-1 @@ -38,7 +38,8 @@ end function ldlt(A::Symmetric{T,M}) where {T<:Real,M<:BandedMatrix{T}} S = typeof(zero(T)/one(T)) - return S == T ? ldlt!(copy(A)) : ldlt!(Symmetric(BandedMatrix{S}(M), A.uplo)) + uplo = A.uplo == 'U' ? :U : :L + return S == T ? ldlt!(copy(A)) : ldlt!(Symmetric(BandedMatrix{S}(A), uplo)) end diff --git a/test/test_symbanded.jl b/test/test_symbanded.jl index 54b562e5..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 @@ -120,8 +120,8 @@ end end @testset "LDLᵀ" begin - for T in (Float16, Float32, Float64, BigFloat) - A = BandedMatrix{T}(undef,(10,10),(2,2)) + 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 @@ -131,6 +131,8 @@ end 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) @@ -141,6 +143,13 @@ end @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