Skip to content

Commit

Permalink
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.
  • Loading branch information
@andreasnoack
andreasnoack committed Aug 9, 2016
1 parent 2037464 commit 79df190
Show file tree
Hide file tree
Showing 4 changed files with 66 additions and 56 deletions.
3 changes: 3 additions & 0 deletions base/deprecated.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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

Expand Down
97 changes: 52 additions & 45 deletions base/linalg/cholesky.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand All @@ -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
Expand Down Expand Up @@ -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')
Expand All @@ -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
Expand All @@ -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
Expand All @@ -235,63 +244,61 @@ 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`.
A `PosDefException` exception is thrown in case the matrix is not positive definite.
"""
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`.
Expand All @@ -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
Expand Down
12 changes: 6 additions & 6 deletions doc/stdlib/linalg.rst
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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

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

0 comments on commit 79df190

Please sign in to comment.