Merge pull request #90 from JuliaMatrices/feat-ldlt

Feat ldlt

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")

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

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