From d992f3dd10602ef52d8c297ebb7c7396a579c963 Mon Sep 17 00:00:00 2001 From: Andreas Noack Date: Wed, 10 Aug 2016 03:18:17 -0400 Subject: [PATCH] Remove redundant uplo argument in chol family. When using Hermitian (#17909) 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. (cherry picked from commit 35df9b806956584043b9ada634124162bccd096e) --- base/deprecated.jl | 3 ++ base/linalg/cholesky.jl | 97 ++++++++++++++++++++++------------------- doc/stdlib/linalg.rst | 12 ++--- test/linalg/cholesky.jl | 10 ++--- 4 files changed, 66 insertions(+), 56 deletions(-) diff --git a/base/deprecated.jl b/base/deprecated.jl index ac4eb4e1c5347..002da08644150 100644 --- a/base/deprecated.jl +++ b/base/deprecated.jl @@ -775,6 +775,9 @@ function transpose(x) return x end +@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) + # During the 0.5 development cycle, do not add any deprecations below this line # To be deprecated in 0.6 diff --git a/base/linalg/cholesky.jl b/base/linalg/cholesky.jl index e961ef09f5a57..ce799f327b782 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,26 +244,25 @@ 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. -The `uplo` argument may be `:L` for using the lower part or `:U` for the upper part of `A`. +and return a `Cholesky` factorization. The matrix `A` can either be a `Symmetric` or `Hermitian` +`StridedMatrix` or a *perfectly* symmetric or Hermitian `StridedMatrix`. In the latter case, +the optional argument `uplo` may be `:L` for using the lower part or `:U` for the upper part of `A`. The default is to use `:U`. The triangular Cholesky factor can be obtained from the factorization `F` with: `F[:L]` and `F[:U]`. The following functions are available for `Cholesky` objects: `size`, `\\`, `inv`, `det`. @@ -262,36 +270,35 @@ 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. -The `uplo` argument may be `:L` for using the lower part or `:U` for the upper part of `A`. +and return a `CholeskyPivoted` factorization. The matrix `A` can either be a `Symmetric` or `Hermitian` +`StridedMatrix` or a *perfectly* symmetric or Hermitian `StridedMatrix`. In the latter case, +the optional argument `uplo` may be `:L` for using the lower part or `:U` for the upper part of `A`. The default is to use `:U`. The triangular Cholesky factor can be obtained from the factorization `F` with: `F[:L]` and `F[:U]`. The following functions are available for `PivotedCholesky` objects: `size`, `\\`, `inv`, `det`, and `rank`. @@ -300,7 +307,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/doc/stdlib/linalg.rst b/doc/stdlib/linalg.rst index 2bd41341c74fb..349869b5a5214 100644 --- a/doc/stdlib/linalg.rst +++ b/doc/stdlib/linalg.rst @@ -308,17 +308,17 @@ Linear algebra functions in Julia are largely implemented by calling functions f Compute the square root of a non-negative number ``x``\ . -.. function:: cholfact(A, uplo::Symbol, Val{false}) -> Cholesky +.. function:: cholfact(A, [uplo::Symbol,] Val{false}) -> Cholesky .. Docstring generated from Julia source - Compute the Cholesky factorization of a dense symmetric positive definite matrix ``A`` and return a ``Cholesky`` factorization. The ``uplo`` argument may be ``:L`` for using the lower part or ``:U`` for the upper part of ``A``\ . The default is to use ``:U``\ . The triangular Cholesky factor can be obtained from the factorization ``F`` with: ``F[:L]`` and ``F[:U]``\ . The following functions are available for ``Cholesky`` objects: ``size``\ , ``\``\ , ``inv``\ , ``det``\ . A ``PosDefException`` exception is thrown in case the matrix is not positive definite. + Compute the Cholesky factorization of a dense symmetric positive definite matrix ``A`` and return a ``Cholesky`` factorization. The matrix ``A`` can either be a ``Symmetric`` or ``Hermitian`` ``StridedMatrix`` or a *perfectly* symmetric or Hermitian ``StridedMatrix``\ . In the latter case, the optional argument ``uplo`` may be ``:L`` for using the lower part or ``:U`` for the upper part of ``A``\ . The default is to use ``:U``\ . The triangular Cholesky factor can be obtained from the factorization ``F`` with: ``F[:L]`` and ``F[:U]``\ . The following functions are available for ``Cholesky`` objects: ``size``\ , ``\``\ , ``inv``\ , ``det``\ . A ``PosDefException`` exception is thrown in case the matrix is not positive definite. -.. function:: cholfact(A, uplo::Symbol, Val{true}; tol = 0.0) -> CholeskyPivoted +.. function:: cholfact(A, [uplo::Symbol,] Val{true}; tol = 0.0) -> CholeskyPivoted .. Docstring generated from Julia source - Compute the pivoted Cholesky factorization of a dense symmetric positive semi-definite matrix ``A`` and return a ``CholeskyPivoted`` factorization. The ``uplo`` argument may be ``:L`` for using the lower part or ``:U`` for the upper part of ``A``\ . The default is to use ``:U``\ . The triangular Cholesky factor can be obtained from the factorization ``F`` with: ``F[:L]`` and ``F[:U]``\ . The following functions are available for ``PivotedCholesky`` objects: ``size``\ , ``\``\ , ``inv``\ , ``det``\ , and ``rank``\ . The argument ``tol`` determines the tolerance for determining the rank. For negative values, the tolerance is the machine precision. + Compute the pivoted Cholesky factorization of a dense symmetric positive semi-definite matrix ``A`` and return a ``CholeskyPivoted`` factorization. The matrix ``A`` can either be a ``Symmetric`` or ``Hermitian`` ``StridedMatrix`` or a *perfectly* symmetric or Hermitian ``StridedMatrix``\ . In the latter case, the optional argument ``uplo`` may be ``:L`` for using the lower part or ``:U`` for the upper part of ``A``\ . The default is to use ``:U``\ . The triangular Cholesky factor can be obtained from the factorization ``F`` with: ``F[:L]`` and ``F[:U]``\ . The following functions are available for ``PivotedCholesky`` objects: ``size``\ , ``\``\ , ``inv``\ , ``det``\ , and ``rank``\ . The argument ``tol`` determines the tolerance for determining the rank. For negative values, the tolerance is the machine precision. .. function:: cholfact(A; shift = 0.0, perm = Int[]) -> CHOLMOD.Factor @@ -344,13 +344,13 @@ Linear algebra functions in Julia are largely implemented by calling functions f This method uses the CHOLMOD library from SuiteSparse, which only supports doubles or complex doubles. Input matrices not of those element types will be converted to ``SparseMatrixCSC{Float64}`` or ``SparseMatrixCSC{Complex128}`` as appropriate. -.. function:: cholfact!(A, uplo::Symbol, Val{false}) -> Cholesky +.. function:: cholfact!(A, [uplo::Symbol,] Val{false}) -> Cholesky .. Docstring generated from Julia source 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 produces a number not representable by the element type of ``A``\ , e.g. for integer types. -.. function:: cholfact!(A, uplo::Symbol, Val{true}; tol = 0.0) -> CholeskyPivoted +.. function:: cholfact!(A, [uplo::Symbol,] Val{true}; tol = 0.0) -> CholeskyPivoted .. Docstring generated from Julia source diff --git a/test/linalg/cholesky.jl b/test/linalg/cholesky.jl index 00e4696b4f12b..bb0657c799838 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