From b0cab7e733005e97f911a334b47ee2b6698de12f Mon Sep 17 00:00:00 2001 From: Andreas Noack Date: Mon, 8 Aug 2016 23:10:51 -0400 Subject: [PATCH] Remove redundant uplo argument in chol family. When using Hermitian and Symmetric, it is redundant to also have an uplo argument and occasionally, it also gave the wrong result. This can therefore be considered a bugfix. --- base/deprecated.jl | 3 ++ base/linalg/cholesky.jl | 87 ++++++++++++++++++++++------------------- test/linalg/cholesky.jl | 10 ++--- 3 files changed, 54 insertions(+), 46 deletions(-) diff --git a/base/deprecated.jl b/base/deprecated.jl index ac4eb4e1c53471..939dd26940a624 100644 --- a/base/deprecated.jl +++ b/base/deprecated.jl @@ -780,4 +780,7 @@ end const _oldstyle_array_vcat_ = false +@deprecate cholfact!(A::Base.LinAlg.HermOrSym, uplo::Symbol, ::Type{Val{false}}) cholfact!(A, Val{false}) +@deprecate cholfact!(A::Base.LinAlg.HermOrSym, uplo::Symbol = :U) cholfact!(A) + # End deprecations scheduled for 0.6 diff --git a/base/linalg/cholesky.jl b/base/linalg/cholesky.jl index e961ef09f5a57d..8a6caeba9f067b 100644 --- a/base/linalg/cholesky.jl +++ b/base/linalg/cholesky.jl @@ -121,8 +121,10 @@ non_hermitian_error(f) = throw(ArgumentError("matrix is not symmetric/" * # chol!. Destructive methods for computing Cholesky factor of real symmetric or Hermitian # matrix -chol!(A::Hermitian) = _chol!(A.data, UpperTriangular) -chol!{T<:Real,S<:StridedMatrix}(A::Symmetric{T,S}) = _chol!(A.data, UpperTriangular) +chol!(A::Hermitian) = + _chol!(A.uplo == 'U' ? A.data : LinAlg.copytri!(A.data, 'L', true), UpperTriangular) +chol!{T<:Real,S<:StridedMatrix}(A::Symmetric{T,S}) = + _chol!(A.uplo == 'U' ? A.data : LinAlg.copytri!(A.data, 'L', true), UpperTriangular) function chol!(A::StridedMatrix) ishermitian(A) || non_hermitian_error("chol!") return _chol!(A, UpperTriangular) @@ -135,14 +137,22 @@ end function chol(A::Hermitian) T = promote_type(typeof(chol(one(eltype(A)))), Float32) AA = similar(A, T, size(A)) - copy!(AA, A.data) - chol!(Hermitian(AA)) + if A.uplo == 'U' + copy!(AA, A.data) + else + Base.ccopy!(AA, A.data) + end + chol!(Hermitian(AA, :U)) end function chol{T<:Real,S<:AbstractMatrix}(A::Symmetric{T,S}) TT = promote_type(typeof(chol(one(T))), Float32) AA = similar(A, TT, size(A)) - copy!(AA, A.data) - chol!(Hermitian(AA)) + if A.uplo == 'U' + copy!(AA, A.data) + else + Base.ccopy!(AA, A.data) + end + chol!(Hermitian(AA, :U)) end ## for StridedMatrices, check that matrix is symmetric/Hermitian @@ -170,15 +180,15 @@ chol(x::Number, args...) = _chol!(x, nothing) # cholfact!. Destructive methods for computing Cholesky factorization of real symmetric # or Hermitian matrix ## No pivoting -function cholfact!(A::Hermitian, uplo::Symbol, ::Type{Val{false}}) - if uplo == :U +function cholfact!(A::Hermitian, ::Type{Val{false}}) + if A.uplo == :U Cholesky(_chol!(A.data, UpperTriangular).data, 'U') else Cholesky(_chol!(A.data, LowerTriangular).data, 'L') end end -function cholfact!{T<:Real,S}(A::Symmetric{T,S}, uplo::Symbol, ::Type{Val{false}}) - if uplo == :U +function cholfact!{T<:Real,S}(A::Symmetric{T,S}, ::Type{Val{false}}) + if A.uplo == :U Cholesky(_chol!(A.data, UpperTriangular).data, 'U') else Cholesky(_chol!(A.data, LowerTriangular).data, 'L') @@ -187,7 +197,7 @@ end ### for StridedMatrices, check that matrix is symmetric/Hermitian """ - cholfact!(A, uplo::Symbol, Val{false}) -> Cholesky + cholfact!(A, [uplo::Symbol,] Val{false}) -> Cholesky The same as `cholfact`, but saves space by overwriting the input `A`, instead of creating a copy. An `InexactError` exception is thrown if the factorisation @@ -196,37 +206,36 @@ integer types. """ function cholfact!(A::StridedMatrix, uplo::Symbol, ::Type{Val{false}}) ishermitian(A) || non_hermitian_error("cholfact!") - return cholfact!(Hermitian(A), uplo, Val{false}) + return cholfact!(Hermitian(A, uplo), Val{false}) end ### Default to no pivoting (and storing of upper factor) when not explicit -cholfact!(A::Hermitian, uplo::Symbol = :U) = cholfact!(A, uplo, Val{false}) -cholfact!{T<:Real,S}(A::Symmetric{T,S}, uplo::Symbol = :U) = cholfact!(A, uplo, Val{false}) +cholfact!(A::Hermitian) = cholfact!(A, Val{false}) +cholfact!{T<:Real,S}(A::Symmetric{T,S}) = cholfact!(A, Val{false}) #### for StridedMatrices, check that matrix is symmetric/Hermitian function cholfact!(A::StridedMatrix, uplo::Symbol = :U) ishermitian(A) || non_hermitian_error("cholfact!") - return cholfact!(Hermitian(A), uplo) + return cholfact!(Hermitian(A, uplo)) end ## With pivoting ### BLAS/LAPACK element types function cholfact!{T<:BlasReal,S<:StridedMatrix}(A::RealHermSymComplexHerm{T,S}, - uplo::Symbol, ::Type{Val{true}}; tol = 0.0) - uplochar = char_uplo(uplo) - AA, piv, rank, info = LAPACK.pstrf!(uplochar, A.data, tol) - return CholeskyPivoted{eltype(AA),typeof(AA)}(AA, uplochar, piv, rank, tol, info) + ::Type{Val{true}}; tol = 0.0) + AA, piv, rank, info = LAPACK.pstrf!(A.uplo, A.data, tol) + return CholeskyPivoted{eltype(AA),typeof(AA)}(AA, A.uplo, piv, rank, tol, info) end ### Non BLAS/LAPACK element types (generic). Since generic fallback for pivoted Cholesky ### is not implemented yet we throw an error -cholfact!{T<:Real,S}(A::RealHermSymComplexHerm{T,S}, uplo::Symbol, ::Type{Val{true}}; +cholfact!{T<:Real,S}(A::RealHermSymComplexHerm{T,S}, ::Type{Val{true}}; tol = 0.0) = throw(ArgumentError("generic pivoted Cholesky factorization is not implemented yet")) ### for StridedMatrices, check that matrix is symmetric/Hermitian """ - cholfact!(A, uplo::Symbol, Val{true}; tol = 0.0) -> CholeskyPivoted + cholfact!(A, [uplo::Symbol,] Val{true}; tol = 0.0) -> CholeskyPivoted The same as `cholfact`, but saves space by overwriting the input `A`, instead of creating a copy. An `InexactError` exception is thrown if the @@ -235,22 +244,20 @@ e.g. for integer types. """ function cholfact!(A::StridedMatrix, uplo::Symbol, ::Type{Val{true}}; tol = 0.0) ishermitian(A) || non_hermitian_error("cholfact!") - return cholfact!(Hermitian(A), uplo, Val{true}; tol = tol) + return cholfact!(Hermitian(A, uplo), Val{true}; tol = tol) end - - # cholfact. Non-destructive methods for computing Cholesky factorization of real symmetric # or Hermitian matrix ## No pivoting -cholfact(A::Hermitian, uplo::Symbol, ::Type{Val{false}}) = - cholfact!(copy_oftype(A, promote_type(typeof(chol(one(eltype(A)))),Float32)), uplo, Val{false}) -cholfact{T<:Real,S<:StridedMatrix}(A::Symmetric{T,S}, uplo::Symbol, ::Type{Val{false}}) = - cholfact!(copy_oftype(A, promote_type(typeof(chol(one(eltype(A)))),Float32)), uplo, Val{false}) +cholfact(A::Hermitian, ::Type{Val{false}}) = + cholfact!(copy_oftype(A, promote_type(typeof(chol(one(eltype(A)))),Float32)), Val{false}) +cholfact{T<:Real,S<:StridedMatrix}(A::Symmetric{T,S}, ::Type{Val{false}}) = + cholfact!(copy_oftype(A, promote_type(typeof(chol(one(eltype(A)))),Float32)), Val{false}) ### for StridedMatrices, check that matrix is symmetric/Hermitian """ - cholfact(A, uplo::Symbol, Val{false}) -> Cholesky + cholfact(A, [uplo::Symbol,] Val{false}) -> Cholesky Compute the Cholesky factorization of a dense symmetric positive definite matrix `A` and return a `Cholesky` factorization. @@ -262,32 +269,30 @@ A `PosDefException` exception is thrown in case the matrix is not positive defin """ function cholfact(A::StridedMatrix, uplo::Symbol, ::Type{Val{false}}) ishermitian(A) || non_hermitian_error("cholfact") - return cholfact(Hermitian(A), uplo, Val{false}) + return cholfact(Hermitian(A, uplo), Val{false}) end ### Default to no pivoting (and storing of upper factor) when not explicit -cholfact(A::Hermitian, uplo::Symbol = :U) = cholfact(A, uplo, Val{false}) -cholfact{T<:Real,S<:StridedMatrix}(A::Symmetric{T,S}, uplo::Symbol = :U) = - cholfact(A, uplo, Val{false}) +cholfact(A::Hermitian) = cholfact(A, Val{false}) +cholfact{T<:Real,S<:StridedMatrix}(A::Symmetric{T,S}) = cholfact(A, Val{false}) #### for StridedMatrices, check that matrix is symmetric/Hermitian function cholfact(A::StridedMatrix, uplo::Symbol = :U) ishermitian(A) || non_hermitian_error("cholfact") - return cholfact(Hermitian(A), uplo) + return cholfact(Hermitian(A, uplo)) end ## With pivoting -cholfact(A::Hermitian, uplo::Symbol, ::Type{Val{true}}; tol = 0.0) = +cholfact(A::Hermitian, ::Type{Val{true}}; tol = 0.0) = cholfact!(copy_oftype(A, promote_type(typeof(chol(one(eltype(A)))),Float32)), - uplo, Val{true}; tol = tol) -cholfact{T<:Real,S<:StridedMatrix}(A::RealHermSymComplexHerm{T,S}, uplo::Symbol, - ::Type{Val{true}}; tol = 0.0) = + Val{true}; tol = tol) +cholfact{T<:Real,S<:StridedMatrix}(A::RealHermSymComplexHerm{T,S}, ::Type{Val{true}}; tol = 0.0) = cholfact!(copy_oftype(A, promote_type(typeof(chol(one(eltype(A)))),Float32)), - uplo, Val{true}; tol = tol) + Val{true}; tol = tol) ### for StridedMatrices, check that matrix is symmetric/Hermitian """ - cholfact(A, uplo::Symbol, Val{true}; tol = 0.0) -> CholeskyPivoted + cholfact(A, [uplo::Symbol,] Val{true}; tol = 0.0) -> CholeskyPivoted Compute the pivoted Cholesky factorization of a dense symmetric positive semi-definite matrix `A` and return a `CholeskyPivoted` factorization. @@ -300,7 +305,7 @@ For negative values, the tolerance is the machine precision. """ function cholfact(A::StridedMatrix, uplo::Symbol, ::Type{Val{true}}; tol = 0.0) ishermitian(A) || non_hermitian_error("cholfact") - return cholfact(Hermitian(A), uplo, Val{true}; tol = tol) + return cholfact(Hermitian(A, uplo), Val{true}; tol = tol) end ## Number diff --git a/test/linalg/cholesky.jl b/test/linalg/cholesky.jl index 00e4696b4f12ba..bb0657c7998389 100644 --- a/test/linalg/cholesky.jl +++ b/test/linalg/cholesky.jl @@ -74,15 +74,15 @@ for eltya in (Float32, Float64, Complex64, Complex128, BigFloat, Int) @test full(lapd) ≈ apd l = lapd[:L] @test l*l' ≈ apd - @test triu(capd.factors) ≈ lapd[:U] - @test tril(lapd.factors) ≈ capd[:L] + @test capd[:U] ≈ lapd[:U] + @test lapd[:L] ≈ capd[:L] if eltya <: Real capds = cholfact(apds) - lapds = cholfact(apds, :L) + lapds = cholfact(full(apds), :L) ls = lapds[:L] @test ls*ls' ≈ apd - @test triu(capds.factors) ≈ lapds[:U] - @test tril(lapds.factors) ≈ capds[:L] + @test capds[:U] ≈ lapds[:U] + @test lapds[:L] ≈ capds[:L] end #pivoted upper Cholesky