diff --git a/src/bidiag.jl b/src/bidiag.jl index b1fe15e9..cc5e6de7 100644 --- a/src/bidiag.jl +++ b/src/bidiag.jl @@ -1467,6 +1467,35 @@ function inv(B::Bidiagonal{T}) where T return B.uplo == 'U' ? UpperTriangular(dest) : LowerTriangular(dest) end +# cholesky-version for (sym)tridiagonal matrices +for (T, uplo) in ((:UpperTriangular, :(:U)), (:LowerTriangular, :(:L))) + @eval function _chol!(A::Bidiagonal, ::Type{$T}) + dv = real(A.dv) + ev = A.ev + n = length(dv) + @inbounds for i in 1:n-1 + iszero(dv[i]) && throw(ZeroPivotException(i)) + ev[i] /= dv[i] + dv[i+1] -= abs2(ev[i])*dv[i] + Akk = dv[i] + Akk, info = _chol!(Akk, UpperTriangular) + if info != 0 + return $T(A), convert(BlasInt, i) + end + dv[i] = Akk + ev[i] *= dv[i] + end + Akk = dv[n] + Akk, info = _chol!(Akk, $T) + if info != 0 + return $T(A), convert(BlasInt, n) + end + dv[n] = Akk + B = Bidiagonal(dv, ev, $uplo) + return $T(B), convert(BlasInt, 0) + end +end + # Eigensystems eigvals(M::Bidiagonal) = copy(M.dv) function eigvecs(M::Bidiagonal{T}) where T diff --git a/src/special.jl b/src/special.jl index 1327af67..83880ca5 100644 --- a/src/special.jl +++ b/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) diff --git a/src/tridiag.jl b/src/tridiag.jl index e31709fa..519be750 100644 --- a/src/tridiag.jl +++ b/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)) diff --git a/test/cholesky.jl b/test/cholesky.jl index 3eacfc0a..8c8919c7 100644 --- a/test/cholesky.jl +++ b/test/cholesky.jl @@ -492,9 +492,41 @@ 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 + M[1,1] *= -1 + @test_throws PosDefException cholesky(M) + C = cholesky(M, check=false) + @test C.info > 0 + M[1,1] *= -1 + M[end,end] *= -1 + @test_throws PosDefException cholesky(M) + C = cholesky(M, check=false) + @test C.info > 0 + 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 + A = Hermitian(Tridiagonal(randn(Quaternion{Float64}, 4, 4) |> t -> t't), :L) + C = cholesky(A) + @test C.L * C.U ≈ A + @test parent(C.U) isa Bidiagonal end @testset "constructor with non-BlasInt arguments" begin