Skip to content

Commit

Permalink
Deprecate chol(A, Val{:U/:L}) and move documentation to method defini…
Browse files Browse the repository at this point in the history 
…tion
  • Loading branch information
@andreasnoack
andreasnoack committed Oct 21, 2015
1 parent 5d7f2de commit 6457fec
Show file tree
Hide file tree
Showing 6 changed files with 69 additions and 62 deletions.
5 changes: 5 additions & 0 deletions base/deprecated.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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))
7 changes: 0 additions & 7 deletions base/docs/helpdb.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand Down
96 changes: 50 additions & 46 deletions base/linalg/cholesky.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand All @@ -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
Expand All @@ -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)
Expand All @@ -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

Expand All @@ -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)
Expand Down
10 changes: 8 additions & 2 deletions doc/stdlib/linalg.rst
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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

Expand Down
9 changes: 4 additions & 5 deletions test/linalg/cholesky.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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()
Expand All @@ -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
4 changes: 2 additions & 2 deletions test/linalg/triangular.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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")

Expand Down Expand Up @@ -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)
Expand Down

0 comments on commit 6457fec

Please sign in to comment.