Add bidiagonal cholesky version

Author: dkarrasch

src/cholesky.jl

@@ -220,6 +220,32 @@ function _chol!(A::AbstractMatrix, ::Type{UpperTriangular})
     end
     return UpperTriangular(A), 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::AbstractMatrix, ::Type{LowerTriangular})
     require_one_based_indexing(A)
     n = checksquare(A)
@@ -250,6 +276,32 @@ function _chol!(A::AbstractMatrix, ::Type{LowerTriangular})
      end
     return LowerTriangular(A), 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
 
 ## Numbers
 function _chol!(x::Number, _)

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,16 @@ 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 M in (SymTridiagonal(fill(2.0, 4), ones(3)),
+        Tridiagonal(ones(3), fill(2.0, 4), ones(3)),
+        )
+        C = cholesky(M)
+        @test C.L * C.U ≈ M
+    end
+    # test LowerTriangular version
+    M = Hermitian(Bidiagonal(fill(2.0, 4), ones(3), :L), :L)
+    C = cholesky!(copy(M))
+    @test C.L * C.U ≈ M
 end
 
 @testset "constructor with non-BlasInt arguments" begin