diff --git a/stdlib/LinearAlgebra/src/cholesky.jl b/stdlib/LinearAlgebra/src/cholesky.jl index 794d6f5959b94..81cc017e2c2f7 100644 --- a/stdlib/LinearAlgebra/src/cholesky.jl +++ b/stdlib/LinearAlgebra/src/cholesky.jl @@ -514,6 +514,32 @@ function ldiv!(C::CholeskyPivoted, B::StridedMatrix) end end +function rdiv!(B::StridedMatrix, C::Cholesky{<:Any,<:AbstractMatrix}) + if C.uplo == 'L' + return rdiv!(rdiv!(B, adjoint(LowerTriangular(C.factors))), LowerTriangular(C.factors)) + else + return rdiv!(rdiv!(B, UpperTriangular(C.factors)), adjoint(UpperTriangular(C.factors))) + end +end + +function LinearAlgebra.rdiv!(B::StridedMatrix, C::CholeskyPivoted) + n = size(C, 2) + for i in 1:size(B, 1) + permute!(view(B, i, 1:n), C.piv) + end + if C.uplo == 'L' + rdiv!(rdiv!(B, adjoint(LowerTriangular(C.factors))), + LowerTriangular(C.factors)) + else + rdiv!(rdiv!(B, UpperTriangular(C.factors)), + adjoint(UpperTriangular(C.factors))) + end + for i in 1:size(B, 1) + invpermute!(view(B, i, 1:n), C.piv) + end + B +end + isposdef(C::Union{Cholesky,CholeskyPivoted}) = C.info == 0 function det(C::Cholesky) diff --git a/stdlib/LinearAlgebra/test/cholesky.jl b/stdlib/LinearAlgebra/test/cholesky.jl index 558b1439dbb25..dccfb3ec550d1 100644 --- a/stdlib/LinearAlgebra/test/cholesky.jl +++ b/stdlib/LinearAlgebra/test/cholesky.jl @@ -168,6 +168,42 @@ end end end end + + @testset "solve with generic Cholesky" begin + Breal = convert(Matrix{BigFloat}, randn(n,n)/2) + Bimg = convert(Matrix{BigFloat}, randn(n,n)/2) + B = eltya <: Complex ? complex.(Breal, Bimg) : Breal + εb = eps(abs(float(one(eltype(B))))) + ε = max(εa,εb) + + for B in (B, view(B, 1:n, 1:n)) # Array and SubArray + + # Test error bound on linear solver: LAWNS 14, Theorem 2.1 + # This is a surprisingly loose bound + cpapd = cholesky(eltya <: Complex ? apdh : apds) + BB = copy(B) + rdiv!(BB, cpapd) + @test norm(B / apd - BB, 1) / norm(BB, 1) <= (3n^2 + n + n^3*ε)*ε/(1-(n+1)*ε)*κ + @test norm(BB * apd - B, 1) / norm(B, 1) <= (3n^2 + n + n^3*ε)*ε/(1-(n+1)*ε)*κ + cpapd = cholesky(eltya <: Complex ? apdhL : apdsL) + BB = copy(B) + rdiv!(BB, cpapd) + @test norm(B / apd - BB, 1) / norm(BB, 1) <= (3n^2 + n + n^3*ε)*ε/(1-(n+1)*ε)*κ + @test norm(BB * apd - B, 1) / norm(B, 1) <= (3n^2 + n + n^3*ε)*ε/(1-(n+1)*ε)*κ + if eltya != BigFloat + cpapd = cholesky(eltya <: Complex ? apdh : apds, Val(true)) + BB = copy(B) + rdiv!(BB, cpapd) + @test norm(B / apd - BB, 1) / norm(BB, 1) <= (3n^2 + n + n^3*ε)*ε/(1-(n+1)*ε)*κ + @test norm(BB * apd - B, 1) / norm(B, 1) <= (3n^2 + n + n^3*ε)*ε/(1-(n+1)*ε)*κ + cpapd = cholesky(eltya <: Complex ? apdhL : apdsL, Val(true)) + BB = copy(B) + rdiv!(BB, cpapd) + @test norm(B / apd - BB, 1) / norm(BB, 1) <= (3n^2 + n + n^3*ε)*ε/(1-(n+1)*ε)*κ + @test norm(BB * apd - B, 1) / norm(B, 1) <= (3n^2 + n + n^3*ε)*ε/(1-(n+1)*ε)*κ + end + end + end if eltya <: BlasFloat @testset "generic cholesky!" begin if eltya <: Complex