Skip to content

Commit

Permalink
Adjust signatures for bkfact to use Symmetric and Hermitian
Browse files Browse the repository at this point in the history
  • Loading branch information
@andreasnoack
andreasnoack committed Jun 29, 2017
1 parent 2c85595 commit 3563f5e
Show file tree
Hide file tree
Showing 4 changed files with 71 additions and 56 deletions.
12 changes: 12 additions & 0 deletions base/deprecated.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1491,6 +1491,18 @@ export conv, conv2, deconv, filt, filt!, xcorr
@deprecate cov(X::AbstractVector, Y::AbstractVector, corrected::Bool) cov(X, Y, corrected=corrected)
@deprecate cov(X::AbstractVecOrMat, Y::AbstractVecOrMat, vardim::Int, corrected::Bool) cov(X, Y, vardim, corrected=corrected)

# bkfact
function bkfact(A::StridedMatrix, uplo::Symbol, symmetric::Bool = issymmetric(A), rook::Bool = false)
depwarn("bkfact with uplo and symmetric arguments deprecated. Please use bkfact($(symmetric ? "Symmetric(" : "Hermitian(")A, :$uplo))",
:bkfact)
return bkfact(symmetric ? Symmetric(A, uplo) : Hermitian(A, uplo), rook)
end
function bkfact!(A::StridedMatrix, uplo::Symbol, symmetric::Bool = issymmetric(A), rook::Bool = false)
depwarn("bkfact! with uplo and symmetric arguments deprecated. Please use bkfact!($(symmetric ? "Symmetric(" : "Hermitian(")A, :$uplo))",
:bkfact!)
return bkfact!(symmetric ? Symmetric(A, uplo) : Hermitian(A, uplo), rook)
end

# END 0.7 deprecations

# BEGIN 1.0 deprecations
Expand Down
61 changes: 29 additions & 32 deletions base/linalg/bunchkaufman.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -22,36 +22,38 @@ BunchKaufman(A::AbstractMatrix{T}, ipiv::Vector{BlasInt}, uplo::Char, symmetric:
`bkfact!` is the same as [`bkfact`](@ref), but saves space by overwriting the
input `A`, instead of creating a copy.
"""
function bkfact!(A::StridedMatrix{<:BlasReal}, uplo::Symbol = :U,
symmetric::Bool = issymmetric(A), rook::Bool = false)

if !symmetric
throw(ArgumentError("Bunch-Kaufman decomposition is only valid for symmetric matrices"))
function bkfact!(A::HermOrSym{T,S} where {T<:BlasReal,S<:StridedMatrix{T}}, rook::Bool = false)
if rook
LD, ipiv, info = LAPACK.sytrf_rook!(A.uplo, A.data)
else
LD, ipiv, info = LAPACK.sytrf!(A.uplo, A.data)
end
BunchKaufman(LD, ipiv, A.uplo, true, rook, info)
end
function bkfact!(A::Symmetric{T,S} where {T<:BlasComplex,S<:StridedMatrix{T}}, rook::Bool = false)
if rook
LD, ipiv, info = LAPACK.sytrf_rook!(char_uplo(uplo), A)
LD, ipiv, info = LAPACK.sytrf_rook!(A.uplo, A.data)
else
LD, ipiv, info = LAPACK.sytrf!(char_uplo(uplo), A)
LD, ipiv, info = LAPACK.sytrf!(A.uplo, A.data)
end
BunchKaufman(LD, ipiv, char_uplo(uplo), symmetric, rook, info)
BunchKaufman(LD, ipiv, A.uplo, true, rook, info)
end
function bkfact!(A::StridedMatrix{<:BlasComplex}, uplo::Symbol=:U,
symmetric::Bool=issymmetric(A), rook::Bool=false)

function bkfact!(A::Hermitian{T,S} where {T<:BlasComplex,S<:StridedMatrix{T}}, rook::Bool = false)
if rook
if symmetric
LD, ipiv, info = LAPACK.sytrf_rook!(char_uplo(uplo), A)
else
LD, ipiv, info = LAPACK.hetrf_rook!(char_uplo(uplo), A)
end
LD, ipiv, info = LAPACK.hetrf_rook!(A.uplo, A.data)
else
if symmetric
LD, ipiv, info = LAPACK.sytrf!(char_uplo(uplo), A)
else
LD, ipiv, info = LAPACK.hetrf!(char_uplo(uplo), A)
end
LD, ipiv, info = LAPACK.hetrf!(A.uplo, A.data)
end
BunchKaufman(LD, ipiv, A.uplo, false, rook, info)
end
function bkfact!(A::StridedMatrix{<:BlasFloat}, rook::Bool = false)
if ishermitian(A)
return bkfact!(Hermitian(A), rook)
elseif issymmetric(A)
return bkfact!(Symmetric(A), rook)
else
throw(ArgumentError("Bunch-Kaufman decomposition is only valid for symmetric or Hermitian matrices"))
end
BunchKaufman(LD, ipiv, char_uplo(uplo), symmetric, rook, info)
end

"""
Expand All @@ -69,13 +71,8 @@ The following functions are available for
[^Bunch1977]: J R Bunch and L Kaufman, Some stable methods for calculating inertia and solving symmetric linear systems, Mathematics of Computation 31:137 (1977), 163-179. [url](http://www.ams.org/journals/mcom/1977-31-137/S0025-5718-1977-0428694-0/).
"""
bkfact(A::StridedMatrix{<:BlasFloat}, uplo::Symbol=:U, symmetric::Bool=issymmetric(A),
rook::Bool=false) =
bkfact!(copy(A), uplo, symmetric, rook)
bkfact(A::StridedMatrix{T}, uplo::Symbol=:U, symmetric::Bool=issymmetric(A),
rook::Bool=false) where {T} =
bkfact!(convert(Matrix{promote_type(Float32, typeof(sqrt(one(T))))}, A),
uplo, symmetric, rook)
bkfact(A::AbstractMatrix{T}, rook::Bool=false) where {T} =
bkfact!(copy_oftype(A, typeof(sqrt(one(T)))), rook)

convert(::Type{BunchKaufman{T}}, B::BunchKaufman{T}) where {T} = B
convert(::Type{BunchKaufman{T}}, B::BunchKaufman) where {T} =
Expand All @@ -90,9 +87,9 @@ ishermitian(B::BunchKaufman) = !B.symmetric

issuccess(B::BunchKaufman) = B.info == 0

function Base.show(io::IO, F::BunchKaufman)
println(io, summary(F))
print(io, "successful: $(issuccess(F))")
function Base.show(io::IO, mime::MIME{Symbol("text/plain")}, B::BunchKaufman)
println(io, summary(B))
print(io, "successful: $(issuccess(B))")
end

function inv(B::BunchKaufman{<:BlasReal})
Expand Down
1 change: 0 additions & 1 deletion base/linalg/symmetric.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -330,7 +330,6 @@ for T in (:Symmetric, :Hermitian), op in (:+, :-, :*, :/)
@eval ($op)(A::$T, x::$S) = ($T)(($op)(A.data, x), Symbol(A.uplo))
end

bkfact(A::HermOrSym) = bkfact(A.data, Symbol(A.uplo), issymmetric(A))
factorize(A::HermOrSym) = bkfact(A)

det(A::RealHermSymComplexHerm) = real(det(bkfact(A)))
Expand Down
53 changes: 30 additions & 23 deletions test/linalg/bunchkaufman.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -23,16 +23,17 @@ bimg = randn(n,2)/2
a = eltya == Int ? rand(1:7, n, n) : convert(Matrix{eltya}, eltya <: Complex ? complex.(areal, aimg) : areal)
a2 = eltya == Int ? rand(1:7, n, n) : convert(Matrix{eltya}, eltya <: Complex ? complex.(a2real, a2img) : a2real)
@testset for atype in ("Array", "SubArray")
asym = a'+a # symmetric indefinite
apd = a'*a # symmetric positive-definite
asym = a.'+ a # symmetric indefinite
aher = a' + a # Hermitian indefinite
apd = a' * a # Positive-definite
if atype == "Array"
a = a
a = a
a2 = a2
else
a = view(a, 1:n, 1:n)
a2 = view(a2, 1:n, 1:n)
asym = view(asym, 1:n, 1:n)
apd = view(apd, 1:n, 1:n)
a = view(a , 1:n, 1:n)
a2 = view(a2 , 1:n, 1:n)
aher = view(aher, 1:n, 1:n)
apd = view(apd , 1:n, 1:n)
end
ε = εa = eps(abs(float(one(eltya))))

Expand All @@ -48,36 +49,40 @@ bimg = randn(n,2)/2
εb = eps(abs(float(one(eltyb))))
ε = max(εa,εb)

@testset "(Automatic) Bunch-Kaufman factor of indefinite matrix" begin
bc1 = factorize(asym)
# check that factorize gives a Bunch-Kaufman
@test isa(factorize(asym), LinAlg.BunchKaufman)
@test isa(factorize(aher), LinAlg.BunchKaufman)

@testset "$uplo Bunch-Kaufman factor of indefinite matrix" for uplo in (:L, :U)
bc1 = bkfact(Hermitian(aher, uplo))
@test LinAlg.issuccess(bc1)
@test logabsdet(bc1)[1] ≈ log(abs(det(bc1)))
if eltya <: Real
@test logabsdet(bc1)[2] == sign(det(bc1))
else
@test logabsdet(bc1)[2] ≈ sign(det(bc1))
end
@test inv(bc1)*asym ≈ eye(n)
@test asym*(bc1\b) ≈ b atol=1000ε
@test inv(bc1)*aher ≈ eye(n)
@test aher*(bc1\b) ≈ b atol=1000ε
@testset for rook in (false, true)
@test inv(bkfact(a.'+a, :U, true, rook))*(a.'+a) ≈ eye(n)
@test inv(bkfact(Symmetric(a.' + a, uplo), rook))*(a.' + a) ≈ eye(n)
@test size(bc1) == size(bc1.LD)
@test size(bc1,1) == size(bc1.LD,1)
@test size(bc1,2) == size(bc1.LD,2)
@test size(bc1, 1) == size(bc1.LD, 1)
@test size(bc1, 2) == size(bc1.LD, 2)
if eltya <: BlasReal
@test_throws ArgumentError bkfact(a)
end
end
end

@testset "Bunch-Kaufman factors of a pos-def matrix" begin
@testset for rook in (false, true)
bc2 = bkfact(apd, :U, issymmetric(apd), rook)
@testset "$uplo Bunch-Kaufman factors of a pos-def matrix" for uplo in (:U, :L)
@testset "rook pivoting: $rook" for rook in (false, true)
bc2 = bkfact(Hermitian(apd, uplo), rook)
@test LinAlg.issuccess(bc2)
@test logdet(bc2) ≈ log(det(bc2))
@test logabsdet(bc2)[1] ≈ log(abs(det(bc2)))
@test logabsdet(bc2)[2] == sign(det(bc2))
@test inv(bc2)*apd ≈ eye(n)
@test inv(bc2)*apd ≈ eye(eltyb, n)
@test apd*(bc2\b) ≈ b atol=150000ε
@test ishermitian(bc2) == !issymmetric(bc2)
end
Expand All @@ -104,11 +109,13 @@ end
end

@testset for rook in (false, true)
F = bkfact(As, :U, issymmetric(As), rook)
@test !LinAlg.issuccess(F)
@test det(F) == 0
@test_throws LinAlg.SingularException inv(F)
@test_throws LinAlg.SingularException F \ ones(size(As, 1))
@testset for uplo in (:L, :U)
F = bkfact(issymmetric(As) ? Symmetric(As, uplo) : Hermitian(As, uplo), rook)
@test !LinAlg.issuccess(F)
@test det(F) == 0
@test_throws LinAlg.SingularException inv(F)
@test_throws LinAlg.SingularException F \ ones(size(As, 1))
end
end
end
end
Expand Down

0 comments on commit 3563f5e

Please sign in to comment.