Skip to content

Commit

Permalink
This commit does not belong to any branch on this repository, and may belong to a fork outside of the repository.
Remove redundant uplo argument in chol family. When using Hermitian
Browse files Browse the repository at this point in the history
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.
andreasnoack committed Aug 9, 2016

Verified

This commit was created on GitHub.com and signed with GitHub’s verified signature.
1 parent 2037464 commit b0cab7e
Showing 3 changed files with 54 additions and 46 deletions.
3 changes: 3 additions & 0 deletions base/deprecated.jl
Original file line number Diff line number Diff line change
@@ -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
87 changes: 46 additions & 41 deletions base/linalg/cholesky.jl
Original file line number Diff line number Diff line change
@@ -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
10 changes: 5 additions & 5 deletions test/linalg/cholesky.jl
Original file line number Diff line number Diff line change
@@ -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

0 comments on commit b0cab7e

Please sign in to comment.