Add cholesky(::Diagonal), with tests. (JuliaLang#27823)

* Add cholesky(::Diagonal), with tests.

* Wrap value in Cholesky, fix tests, include order.

* Move diagonal.jl include near the other types (cosmetic).

* Handle check=false, add complex w/tests, move include back.

1. Code is moved to diagonal.jl, this looked like the simplest
solution to have all types available for functions.

2. Added case for `check = false`, which does not fail if the diagonal
value is ≤ 0, but records it in `info`.

3. Allow for complex values, only check for `== 0`.

4. Unit tests for `check = false`, also for complex case.

* Fix deprecated syntax in tests.

* Add cholesky!(::Diagonal), generalize tests for positive definiteness.

1. Control flow is simplified by using a loop.

2. Allow for the possibility of types other than Real or Complex (eg
Quaternion).

3. Add tests `InexactError` for some `cholesky!(::Diagonal{Int})`.

* Fix missing default.

stdlib/LinearAlgebra/src/cholesky.jl

@@ -174,7 +174,6 @@ function cholesky!(A::StridedMatrix, ::Val{false}=Val(false); check::Bool = true
  end
 end
 
-
 ## With pivoting
 ### BLAS/LAPACK element types
 function cholesky!(A::RealHermSymComplexHerm{<:BlasReal,<:StridedMatrix},

stdlib/LinearAlgebra/src/diagonal.jl

@@ -487,3 +487,23 @@ end
 *(x::Transpose{<:Any,<:AbstractVector}, D::Diagonal, y::AbstractVector) =
  mapreduce(t -> t[1]*t[2]*t[3], +, zip(x, D.diag, y))
 # TODO: these methods will yield row matrices, rather than adjoint/transpose vectors
+
+function cholesky!(A::Diagonal, ::Val{false} = Val(false); check::Bool = true)
+ info = 0
+ diagonal = A.diag
+ for i in axes(diagonal, 1)
+ d = diagonal[i]
+ if !(d == 0 || (isreal(d) && d < 0))
+ diagonal[i] = √d
+ elseif check
+ throw(PosDefException(i))
+ else
+ info = i
+ break
+ end
+ end
+ Cholesky(A, 'U', convert(BlasInt, info))
+end
+
+cholesky(A::Diagonal, ::Val{false} = Val(false); check::Bool = true) =
+ cholesky!(cholcopy(A), Val(false); check = check)

stdlib/LinearAlgebra/test/cholesky.jl

@@ -271,4 +271,39 @@ end
  @test_throws ArgumentError cholesky!(Hermitian(rand(Float16, 5,5)), Val(true))
 end
 
+@testset "cholesky Diagonal" begin
+ # real
+ d = abs.(randn(3)) .+ 0.1
+ D = Diagonal(d)
+ CD = cholesky(D)
+ @test CD isa Cholesky{Float64}
+ @test CD.U isa UpperTriangular{Float64}
+ @test CD.U == Diagonal(.√d)
+ @test CD.info == 0
+
+ # real, failing
+ @test_throws PosDefException cholesky(Diagonal([1.0, -2.0]))
+ Dnpd = cholesky(Diagonal([1.0, -2.0]); check = false)
+ @test Dnpd.info == 2
+
+ # complex
+ d = cis.(rand(3) .* 2*π)
+ d .*= abs.(randn(3) .+ 0.1)
+ D = Diagonal(d)
+ CD = cholesky(D)
+ @test CD isa Cholesky{Complex{Float64}}
+ @test CD.U isa UpperTriangular{Complex{Float64}}
+ @test CD.U == Diagonal(.√d)
+ @test CD.info == 0
+
+ # complex, failing
+ D[2, 2] = 0.0 + 0im
+ @test_throws PosDefException cholesky(D)
+ Dnpd = cholesky(D; check = false)
+ @test Dnpd.info == 2
+
+ # InexactError for Int
+ @test_throws InexactError cholesky!(Diagonal([2, 1]))
+end
+
 end # module TestCholesky