Add bidiagonal cholesky version

src/bidiag.jl

@@ -1467,6 +1467,92 @@ function inv(B::Bidiagonal{T}) where T
     return B.uplo == 'U' ? UpperTriangular(dest) : LowerTriangular(dest)
 end
 
+# cholesky-version for (sym)tridiagonal matrices
+function _chol!(A::Bidiagonal{<:BlasFloat,<:StridedVector}, ::Type{UpperTriangular})
+    d = real(A.dv)
+    e = A.ev
+    dv, ev = LAPACK.pttrf!(d, e)
+    for k in eachindex(dv)
+        Akk = dv[k]
+        Akk, info = _chol!(Akk, UpperTriangular)
+        if info != 0
+            return UpperTriangular(A), convert(BlasInt, k)
+        end
+        dv[k] = Akk
+    end
+    @views ev .*= dv[1:end-1]
+    U = Bidiagonal(dv, ev, :U)
+    return UpperTriangular(U), convert(BlasInt, 0)
+end
+function _chol!(A::Bidiagonal, ::Type{UpperTriangular})
+    require_one_based_indexing(A)
+    n = checksquare(A)
+    realdiag = eltype(A) <: Complex
+    dv = A.dv
+    ev = A.ev
+    @inbounds begin
+        for k = 1:n
+            Akk = realdiag ? real(dv[k]) : dv[k]
+            if k > 1
+                Akk -= realdiag ? abs2(ev[k-1]) : ev[k-1]'ev[k-1]
+            end
+            dv[k] = Akk
+            Akk, info = _chol!(Akk, UpperTriangular)
+            if info != 0
+                return UpperTriangular(A), convert(BlasInt, k)
+            end
+            dv[k] = Akk
+            AkkInv = inv(copy(Akk'))
+            if k < n
+                ev[k] = AkkInv*ev[k]
+            end
+        end
+    end
+    return UpperTriangular(A), convert(BlasInt, 0)
+end
+function _chol!(A::Bidiagonal{<:BlasFloat,<:StridedVector}, ::Type{LowerTriangular})
+    d = real(A.dv)
+    e = A.ev
+    dv, ev = LAPACK.pttrf!(d, e)
+    for k in eachindex(dv)
+        Akk = dv[k]
+        Akk, info = _chol!(Akk, LowerTriangular)
+        if info != 0
+            return LowerTriangular(A), convert(BlasInt, k)
+        end
+        dv[k] = Akk
+    end
+    @views ev .*= dv[1:end-1]
+    L = Bidiagonal(dv, ev, :L)
+    return LowerTriangular(L), convert(BlasInt, 0)
+end
+function _chol!(A::Bidiagonal, ::Type{LowerTriangular})
+    require_one_based_indexing(A)
+    n = checksquare(A)
+    realdiag = eltype(A) <: Complex
+    dv = A.dv
+    ev = A.ev
+    @inbounds begin
+        for k = 1:n
+            Akk = realdiag ? real(dv[k]) : dv[k]
+            if k > 1
+                Akk -= realdiag ? abs2(ev[k-1]) : ev[k-1]*ev[k-1]'
+            end
+            dv[k] = Akk
+            Akk, info = _chol!(Akk, LowerTriangular)
+            if info != 0
+                return LowerTriangular(A), convert(BlasInt, k)
+            end
+            dv[k] = Akk
+            AkkInv = inv(copy(Akk'))
+            if k < n
+                ev[k] *= AkkInv
+            end
+        end
+     end
+    return LowerTriangular(A), convert(BlasInt, 0)
+end
+
 # Eigensystems
 eigvals(M::Bidiagonal) = copy(M.dv)
 function eigvecs(M::Bidiagonal{T}) where T

src/special.jl

@@ -605,8 +605,8 @@ function _hvcat(rows::Tuple{Vararg{Int}}, A::Union{AbstractArray,AbstractQ,Unifo
     end
 end
 
-# factorizations
-function cholesky(S::RealHermSymComplexHerm{<:Real,<:SymTridiagonal}, ::NoPivot = NoPivot(); check::Bool = true)
+# tridiagonal cholesky factorization
+function cholesky(S::RealSymHermitian{<:BiTriSym}, ::NoPivot = NoPivot(); check::Bool = true)
     T = choltype(S)
     B = Bidiagonal{T}(diag(S, 0), diag(S, S.uplo == 'U' ? 1 : -1), sym_uplo(S.uplo))
     cholesky!(Hermitian(B, sym_uplo(S.uplo)), NoPivot(); check = check)

src/tridiag.jl

@@ -1039,7 +1039,7 @@ function dot(x::AbstractVector, A::Tridiagonal, y::AbstractVector)
     return r
 end
 
-function cholesky(S::SymTridiagonal, ::NoPivot = NoPivot(); check::Bool = true)
+function cholesky(S::Union{SymTridiagonal,Tridiagonal}, ::NoPivot = NoPivot(); check::Bool = true)
     if !ishermitian(S)
         check && checkpositivedefinite(-1)
         return Cholesky(S, 'U', convert(BlasInt, -1))

test/cholesky.jl

@@ -492,9 +492,37 @@ end
 end
 
 @testset "Cholesky for AbstractMatrix" begin
-    S = SymTridiagonal(fill(2.0, 4), ones(3))
-    C = cholesky(S)
-    @test C.L * C.U ≈ S
+    for T in (identity, big), M in (SymTridiagonal(fill(T(2.0), 4), ones(3)),
+        Symmetric(SymTridiagonal(fill(T(2.0), 4), ones(3)), :U),
+        Symmetric(SymTridiagonal(fill(T(2.0), 4), ones(3)), :L),
+        Tridiagonal(ones(3), fill(T(2.0), 4), ones(3)),
+        Hermitian(Tridiagonal(ones(3), fill(T(2.0), 4), ones(3)), :U),
+        Hermitian(Bidiagonal(fill(T(2.0), 4), ones(3), :U), :U),
+        Hermitian(Bidiagonal(fill(T(2.0), 4), ones(3), :U), :L),
+        Hermitian(Bidiagonal(fill(T(2.0), 4), ones(3), :L), :L),
+        )
+        C = cholesky(M)
+        @test C.L * C.U ≈ M
+        @test parent(C.U) isa Bidiagonal
+    end
+    for T in (identity, big), M in (SymTridiagonal(fill(T(-2.0), 4), ones(3)),
+        Symmetric(SymTridiagonal(fill(T(-2.0), 4), ones(3)), :U),
+        Symmetric(SymTridiagonal(fill(T(-2.0), 4), ones(3)), :L),
+        Tridiagonal(ones(3), fill(T(-2.0), 4), ones(3)),
+        Hermitian(Tridiagonal(ones(3), fill(T(-2.0), 4), ones(3)), :U),
+        Hermitian(Bidiagonal(fill(T(-2.0), 4), ones(3), :U), :U),
+        )
+        @test_throws (T == identity ? LAPACKException : PosDefException) cholesky(M)
+    end
+    # test LowerTriangular version
+    M = Hermitian(Bidiagonal(fill(2.0, 4), im * ones(3), :L), :L)
+    C = cholesky!(copy(M))
+    @test C.L * C.U ≈ M
+    # non-(RealOrComplex) eltype
+    A = Tridiagonal(randn(Quaternion{Float64}, 4, 4) |> t -> t't)
+    C = cholesky(A)
+    @test C.L * C.U ≈ A
+    @test parent(C.U) isa Bidiagonal
 end
 
 @testset "constructor with non-BlasInt arguments" begin