add a LAPACK call path

Author: dkarrasch

src/bidiag.jl

@@ -1468,6 +1468,15 @@ function inv(B::Bidiagonal{T}) where T
 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)
+    map!(sqrt, dv)
+    @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)
@@ -1494,6 +1503,15 @@ function _chol!(A::Bidiagonal, ::Type{UpperTriangular})
     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)
+    map!(sqrt, dv)
+    @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)

test/cholesky.jl

@@ -492,14 +492,14 @@ end
 end
 
 @testset "Cholesky for AbstractMatrix" begin
-    for M in (SymTridiagonal(fill(2.0, 4), ones(3)),
-        Symmetric(SymTridiagonal(fill(2.0, 4), ones(3)), :U),
-        Symmetric(SymTridiagonal(fill(2.0, 4), ones(3)), :L),
-        Tridiagonal(ones(3), fill(2.0, 4), ones(3)),
-        Hermitian(Tridiagonal(ones(3), fill(2.0, 4), ones(3)), :U),
-        Hermitian(Bidiagonal(fill(2.0, 4), ones(3), :U), :U),
-        Hermitian(Bidiagonal(fill(2.0, 4), ones(3), :U), :L),
-        Hermitian(Bidiagonal(fill(2.0, 4), ones(3), :L), :L),
+    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