diff --git a/base/deprecated.jl b/base/deprecated.jl index 555bcc2197013..3145991da957c 100644 --- a/base/deprecated.jl +++ b/base/deprecated.jl @@ -858,3 +858,8 @@ end #13496 @deprecate A_ldiv_B!(A::SparseMatrixCSC, B::StridedVecOrMat) A_ldiv_B!(factorize(A), B) + +@deprecate chol(A::Number, ::Type{Val{:U}}) chol(A) +@deprecate chol(A::AbstractMatrix, ::Type{Val{:U}}) chol(A) +@deprecate chol(A::Number, ::Type{Val{:L}}) ctranspose(chol(A)) +@deprecate chol(A::AbstractMatrix, ::Type{Val{:L}}) ctranspose(chol(A)) diff --git a/base/docs/helpdb.jl b/base/docs/helpdb.jl index 21f90f17bad34..666e44d47b3c2 100644 --- a/base/docs/helpdb.jl +++ b/base/docs/helpdb.jl @@ -5734,13 +5734,6 @@ as ``v0``. In general, this cannot be used with empty collections """ foldr(op, itr) -doc""" - chol(A, [LU]) -> F - -Compute the Cholesky factorization of a symmetric positive definite matrix `A` and return the matrix `F`. If `LU` is `Val{:U}` (Upper), `F` is of type `UpperTriangular` and `A = F'*F`. If `LU` is `Val{:L}` (Lower), `F` is of type `LowerTriangular` and `A = F*F'`. `LU` defaults to `Val{:U}`. -""" -chol - doc""" ParseError(msg) diff --git a/base/linalg/cholesky.jl b/base/linalg/cholesky.jl index fe8579add387a..84926d84c8544 100644 --- a/base/linalg/cholesky.jl +++ b/base/linalg/cholesky.jl @@ -23,24 +23,24 @@ function CholeskyPivoted{T}(A::AbstractMatrix{T}, uplo::Char, piv::Vector{BlasIn CholeskyPivoted{T,typeof(A)}(A, uplo, piv, rank, tol, info) end -function chol!{T<:BlasFloat}(A::StridedMatrix{T}, ::Type{Val{:U}}) +function chol!{T<:BlasFloat}(A::StridedMatrix{T}, ::Type{UpperTriangular}) C, info = LAPACK.potrf!('U', A) return @assertposdef UpperTriangular(C) info end -function chol!{T<:BlasFloat}(A::StridedMatrix{T}, ::Type{Val{:L}}) +function chol!{T<:BlasFloat}(A::StridedMatrix{T}, ::Type{LowerTriangular}) C, info = LAPACK.potrf!('L', A) return @assertposdef LowerTriangular(C) info end -chol!(A::StridedMatrix) = chol!(A, Val{:U}) +chol!(A::StridedMatrix) = chol!(A, UpperTriangular) -function chol!{T}(A::AbstractMatrix{T}, ::Type{Val{:U}}) +function chol!{T}(A::AbstractMatrix{T}, ::Type{UpperTriangular}) n = chksquare(A) @inbounds begin for k = 1:n for i = 1:k - 1 A[k,k] -= A[i,k]'A[i,k] end - Akk = chol!(A[k,k], Val{:U}) + Akk = chol!(A[k,k], UpperTriangular) A[k,k] = Akk AkkInv = inv(Akk') for j = k + 1:n @@ -53,14 +53,14 @@ function chol!{T}(A::AbstractMatrix{T}, ::Type{Val{:U}}) end return UpperTriangular(A) end -function chol!{T}(A::AbstractMatrix{T}, ::Type{Val{:L}}) +function chol!{T}(A::AbstractMatrix{T}, ::Type{LowerTriangular}) n = chksquare(A) @inbounds begin for k = 1:n for i = 1:k - 1 A[k,k] -= A[k,i]*A[k,i]' end - Akk = chol!(A[k,k], Val{:L}) + Akk = chol!(A[k,k], LowerTriangular) A[k,k] = Akk AkkInv = inv(Akk) for j = 1:k @@ -78,14 +78,15 @@ function chol!{T}(A::AbstractMatrix{T}, ::Type{Val{:L}}) return LowerTriangular(A) end +doc""" + chol(A::AbstractMatrix) -> U + +Compute the Cholesky factorization of a positive definite matrix `A` and return the UpperTriangular matrix `U` such that `A = U'U`. +""" function chol{T}(A::AbstractMatrix{T}) S = promote_type(typeof(chol(one(T))), Float32) chol!(copy_oftype(A, S)) end -function chol{T}(A::AbstractMatrix{T}, uplo::Union{Type{Val{:L}}, Type{Val{:U}}}) - S = promote_type(typeof(chol(one(T))), Float32) - chol!(copy_oftype(A, S), uplo) -end function chol!(x::Number, uplo) rx = real(x) if rx != abs(x) @@ -94,50 +95,50 @@ function chol!(x::Number, uplo) rxr = sqrt(rx) convert(promote_type(typeof(x), typeof(rxr)), rxr) end -chol(x::Number, uplo::Symbol=:U) = chol!(x, uplo) - -function cholfact!{T<:BlasFloat}(A::StridedMatrix{T}, uplo::Symbol=:U, - pivot::Union{Type{Val{false}},Type{Val{true}}}=Val{false}; tol=0.0) - _cholfact!(A, pivot, uplo, tol=tol) -end -function _cholfact!{T<:BlasFloat}(A::StridedMatrix{T}, ::Type{Val{false}}, uplo::Symbol=:U; tol=0.0) - return Cholesky(chol!(A, Val{uplo}).data, uplo) +doc""" + chol(x::Number) -> y + +Compute the square root of a non-negative number `x`. +""" +chol(x::Number, args...) = chol!(x, nothing) + +cholfact!{T<:BlasFloat}(A::StridedMatrix{T}, uplo::Symbol, ::Type{Val{false}}) = + _cholfact!(A, Val{false}, uplo) +cholfact!{T<:BlasFloat}(A::StridedMatrix{T}, uplo::Symbol, ::Type{Val{true}}; tol = 0.0) = + _cholfact!(A, Val{true}, uplo; tol = tol) +cholfact!{T<:BlasFloat}(A::StridedMatrix{T}, uplo::Symbol = :U) = + _cholfact!(A, Val{false}, uplo) + +function _cholfact!{T<:BlasFloat}(A::StridedMatrix{T}, ::Type{Val{false}}, uplo::Symbol=:U) + if uplo == :U + return Cholesky(chol!(A, UpperTriangular).data, uplo) + else + return Cholesky(chol!(A, LowerTriangular).data, uplo) + end end function _cholfact!{T<:BlasFloat}(A::StridedMatrix{T}, ::Type{Val{true}}, uplo::Symbol=:U; tol=0.0) uplochar = char_uplo(uplo) A, piv, rank, info = LAPACK.pstrf!(uplochar, A, tol) return CholeskyPivoted{T,StridedMatrix{T}}(A, uplochar, piv, rank, tol, info) end -function cholfact!(A::StridedMatrix, uplo::Symbol=:U, - pivot::Union{Type{Val{false}},Type{Val{true}}}=Val{false}; tol=0.0) - Cholesky(chol!(A, Val{uplo}).data, uplo) +function cholfact!(A::StridedMatrix, uplo::Symbol, ::Type{Val{false}}) + if uplo == :U + Cholesky(chol!(A, UpperTriangular).data, 'U') + else + Cholesky(chol!(A, UpperTriangular).data, 'U') + end end +cholfact!(A::StridedMatrix, uplo::Symbol, ::Type{Val{true}}; tol = 0.0) = throw(ArgumentError("generic pivoted Cholesky fectorization is not implemented yet")) +cholfact!(A::StridedMatrix, uplo::Symbol = :U) = cholfact!(A, uplo, Val{false}) -# function cholfact!{T<:BlasFloat,S}(C::Cholesky{T,S}) -# _, info = LAPACK.potrf!(C.uplo, C.factors) -# info[1]>0 && throw(PosDefException(info[1])) -# C -# end - -# cholfact{T<:BlasFloat}(A::StridedMatrix{T}, uplo::Symbol=:U, pivot::Union{Type{Val{false}}, Type{Val{true}}}=Val{false}; tol=0.0) = -# cholfact!(copy(A), uplo, pivot, tol=tol) - - -cholfact{T}(A::StridedMatrix{T}, uplo::Symbol=:U, pivot::Union{Type{Val{false}}, Type{Val{true}}}=Val{false}; tol=0.0) = - cholfact!(copy_oftype(A, promote_type(typeof(chol(one(T))),Float32)), uplo, pivot; tol = tol) -# _cholfact{T<:BlasFloat}(A::StridedMatrix{T}, pivot::Type{Val{true}}, uplo::Symbol=:U; tol=0.0) = -# cholfact!(A, uplo, pivot, tol = tol) -# _cholfact{T<:BlasFloat}(A::StridedMatrix{T}, pivot::Type{Val{false}}, uplo::Symbol=:U; tol=0.0) = -# cholfact!(A, uplo, pivot, tol = tol) - -# _cholfact{T}(A::StridedMatrix{T}, ::Type{Val{false}}, uplo::Symbol=:U; tol=0.0) = -# cholfact!(A, uplo) -# _cholfact{T}(A::StridedMatrix{T}, ::Type{Val{true}}, uplo::Symbol=:U; tol=0.0) = -# throw(ArgumentError("pivot only supported for Float32, Float64, Complex{Float32} and Complex{Float64}")) - +cholfact{T}(A::StridedMatrix{T}, uplo::Symbol, ::Type{Val{false}}) = + cholfact!(copy_oftype(A, promote_type(typeof(chol(one(T))),Float32)), uplo, Val{false}) +cholfact{T}(A::StridedMatrix{T}, uplo::Symbol, ::Type{Val{true}}; tol = 0.0) = + cholfact!(copy_oftype(A, promote_type(typeof(chol(one(T))),Float32)), uplo, Val{true}; tol = tol) +cholfact{T}(A::StridedMatrix{T}, uplo::Symbol = :U) = cholfact(A, uplo, Val{false}) function cholfact(x::Number, uplo::Symbol=:U) - xf = fill(chol!(x, uplo), 1, 1) + xf = fill(chol(x), 1, 1) Cholesky(xf, uplo) end @@ -155,7 +156,10 @@ convert{T}(::Type{CholeskyPivoted{T}},C::CholeskyPivoted) = convert{T}(::Type{Factorization{T}}, C::CholeskyPivoted) = convert(CholeskyPivoted{T}, C) full{T,S}(C::Cholesky{T,S}) = C.uplo == 'U' ? C[:U]'C[:U] : C[:L]*C[:L]' -full(F::CholeskyPivoted) = (ip=invperm(F[:p]); (F[:L] * F[:U])[ip,ip]) +function full(F::CholeskyPivoted) + ip = invperm(F[:p]) + return (F[:L] * F[:U])[ip,ip] +end copy(C::Cholesky) = Cholesky(copy(C.factors), C.uplo) copy(C::CholeskyPivoted) = CholeskyPivoted(copy(C.factors), C.uplo, C.piv, C.rank, C.tol, C.info) diff --git a/doc/stdlib/linalg.rst b/doc/stdlib/linalg.rst index 76135a1999b08..bc4566b3cae1a 100644 --- a/doc/stdlib/linalg.rst +++ b/doc/stdlib/linalg.rst @@ -111,11 +111,17 @@ Linear algebra functions in Julia are largely implemented by calling functions f ``lufact!`` is the same as :func:`lufact`, but saves space by overwriting the input ``A``, instead of creating a copy. For sparse ``A`` the ``nzval`` field is not overwritten but the index fields, ``colptr`` and ``rowval`` are decremented in place, converting from 1-based indices to 0-based indices. -.. function:: chol(A, [LU]) -> F +.. function:: chol(A::AbstractMatrix) -> U .. Docstring generated from Julia source - Compute the Cholesky factorization of a symmetric positive definite matrix ``A`` and return the matrix ``F``\ . If ``LU`` is ``Val{:U}`` (Upper), ``F`` is of type ``UpperTriangular`` and ``A = F'*F``\ . If ``LU`` is ``Val{:L}`` (Lower), ``F`` is of type ``LowerTriangular`` and ``A = F*F'``\ . ``LU`` defaults to ``Val{:U}``\ . + Compute the Cholesky factorization of a positive definite matrix ``A`` and return the UpperTriangular matrix ``U`` such that ``A = U'U``\ . + +.. function:: chol(x::Number) -> y + + .. Docstring generated from Julia source + + Compute the square root of a non-negative number ``x``\ . .. function:: cholfact(A, [LU=:U[,pivot=Val{false}]][;tol=-1.0]) -> Cholesky diff --git a/test/linalg/cholesky.jl b/test/linalg/cholesky.jl index 2e0e4887a1af8..0b9d7c6418d29 100644 --- a/test/linalg/cholesky.jl +++ b/test/linalg/cholesky.jl @@ -125,8 +125,8 @@ begin # Cholesky factor of Matrix with non-commutative elements, here 2x2-matrices X = Matrix{Float64}[0.1*rand(2,2) for i in 1:3, j = 1:3] - L = full(Base.LinAlg.chol!(X*X', Val{:L})) - U = full(Base.LinAlg.chol!(X*X', Val{:U})) + L = full(Base.LinAlg.chol!(X*X', LowerTriangular)) + U = full(Base.LinAlg.chol!(X*X', UpperTriangular)) XX = full(X*X') @test sum(sum(norm, L*L' - XX)) < eps() @@ -141,7 +141,6 @@ for elty in (Float32, Float64, Complex{Float32}, Complex{Float64}) A = randn(5,5) end A = convert(Matrix{elty}, A'A) - for ul in (:U, :L) - @test_approx_eq full(cholfact(A, ul)[ul]) full(invoke(Base.LinAlg.chol!, Tuple{AbstractMatrix, Type{Val{ul}}},copy(A), Val{ul})) - end + @test_approx_eq full(cholfact(A)[:L]) full(invoke(Base.LinAlg.chol!, Tuple{AbstractMatrix, Type{LowerTriangular}}, copy(A), LowerTriangular)) + @test_approx_eq full(cholfact(A)[:U]) full(invoke(Base.LinAlg.chol!, Tuple{AbstractMatrix, Type{UpperTriangular}}, copy(A), UpperTriangular)) end diff --git a/test/linalg/triangular.jl b/test/linalg/triangular.jl index 51a25a98ce4e3..01f92067d3221 100644 --- a/test/linalg/triangular.jl +++ b/test/linalg/triangular.jl @@ -22,7 +22,7 @@ for elty1 in (Float32, Float64, Complex64, Complex128, BigFloat, Int) (UnitLowerTriangular, :L)) # Construct test matrix - A1 = t1(elty1 == Int ? rand(1:7, n, n) : convert(Matrix{elty1}, (elty1 <: Complex ? complex(randn(n, n), randn(n, n)) : randn(n, n)) |> t -> chol(t't, Val{uplo1}))) + A1 = t1(elty1 == Int ? rand(1:7, n, n) : convert(Matrix{elty1}, (elty1 <: Complex ? complex(randn(n, n), randn(n, n)) : randn(n, n)) |> t -> chol(t't) |> t -> uplo1 == :U ? t : ctranspose(t))) debug && println("elty1: $elty1, A1: $t1") @@ -243,7 +243,7 @@ for elty1 in (Float32, Float64, Complex64, Complex128, BigFloat, Int) debug && println("elty1: $elty1, A1: $t1, elty2: $elty2") - A2 = t2(elty2 == Int ? rand(1:7, n, n) : convert(Matrix{elty2}, (elty2 <: Complex ? complex(randn(n, n), randn(n, n)) : randn(n, n)) |> t-> chol(t't, Val{uplo2}))) + A2 = t2(elty2 == Int ? rand(1:7, n, n) : convert(Matrix{elty2}, (elty2 <: Complex ? complex(randn(n, n), randn(n, n)) : randn(n, n)) |> t -> chol(t't) |> t -> uplo2 == :U ? t : ctranspose(t))) # Convert if elty1 <: Real && !(elty2 <: Integer)