Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

Adjust signatures for bkfact to use Symmetric and Hermitian #22605

Merged
merged 2 commits into from
Jun 30, 2017
Merged

Adjust signatures for bkfact to use Symmetric and Hermitian #22605

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
31cdb97
Adjust signatures for bkfact to use Symmetric and Hermitian
andreasnoack Jun 29, 2017
d605eea
Add NEWS.md entry
andreasnoack Jun 29, 2017
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
3 changes: 3 additions & 0 deletions NEWS.md
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -105,6 +105,9 @@ Deprecated or removed
have been deprecated in favor of `isposdef(Hermitian(A, UL))` and `isposdef!(Hermitian(A, UL))`
respectively ([#22245]).

* The `bkfact`/`bkfact!` methods that accepted `uplo` and `issymmetric` symbols have been deprecated
in favor of using `Hermitian` (or `Symmetric`) views ([#22605]).

* The function `current_module` is deprecated and replaced with `@__MODULE__` ([#22064]).
This caused the deprecation of some reflection methods (such as `macroexpand` and `isconst`),
which now require a module argument.
Expand Down
12 changes: 12 additions & 0 deletions base/deprecated.jl
View file
Open in desktop
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
55 changes: 18 additions & 37 deletions base/linalg/bunchkaufman.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -22,36 +22,22 @@ 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"))
end
if rook
LD, ipiv, info = LAPACK.sytrf_rook!(char_uplo(uplo), A)
else
LD, ipiv, info = LAPACK.sytrf!(char_uplo(uplo), A)
end
BunchKaufman(LD, ipiv, char_uplo(uplo), symmetric, rook, info)
function bkfact!(A::RealHermSymComplexSym{T,S} where {T<:BlasReal,S<:StridedMatrix}, rook::Bool = false)
LD, ipiv, info = rook ? LAPACK.sytrf_rook!(A.uplo, A.data) : LAPACK.sytrf!(A.uplo, A.data)
BunchKaufman(LD, ipiv, A.uplo, true, rook, info)
end
function bkfact!(A::StridedMatrix{<:BlasComplex}, uplo::Symbol=:U,
symmetric::Bool=issymmetric(A), 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
function bkfact!(A::Hermitian{T,S} where {T<:BlasComplex,S<:StridedMatrix{T}}, rook::Bool = false)
LD, ipiv, info = rook ? LAPACK.hetrf_rook!(A.uplo, A.data) : LAPACK.hetrf!(A.uplo, A.data)
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
if symmetric
LD, ipiv, info = LAPACK.sytrf!(char_uplo(uplo), A)
else
LD, ipiv, info = LAPACK.hetrf!(char_uplo(uplo), A)
end
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 +55,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)
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Removing the former promotion with Float32 to support e.g. Float16? :)

Copy link
Contributor

@tkelman tkelman Jun 29, 2017 •
edited
Loading

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

good eye, missed that change - might want to throw in a test for this bit if there's any behavior that this enables right now

Copy link
Member Author

@andreasnoack andreasnoack Jun 29, 2017 •
edited
Loading

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Hm. I had integer and rational input in mind and then sqrt(one(T)) should be sufficient to get a BLAS type but you are right that a matrix with Float16 would stay Float16. Probably not a big problem, though. If you do arithmetic with Float16 you probably would like to avoid that your elements get promoted silently.

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

so will integer or rational input now behave differently with this?

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

No. They are unaffected since sqrt(one(T)) isa Float64 for T in (Rational{Int},Int). A Matrix{Float16} will now fail instead of being promoted. I think that is reasonable because people using Float16 most likely prefer an error over silent promotion and I really don't think anybody was ever called bkfact with a Matrix{Float16} in serious code.

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Float16 is used somewhat frequently in the gpu community. What will the error be in that case? We silently promote for intermediate calculations pretty often and convert back to the destination type when an implementation isn't available in a specific element type.

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

We don't do this for factorizations and there hasn't been expressed any demand for this so I think we good for now. We can revisit if Float16 users start asking for this behavior.

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

can you answer the question at least?

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

don't all the other factorizations have generic implementations so they wouldn't error?

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

You'll get a MethodError in bkfact!. Not all factorizations have generic fallbacks, e.g. eig and svd.


convert(::Type{BunchKaufman{T}}, B::BunchKaufman{T}) where {T} = B
convert(::Type{BunchKaufman{T}}, B::BunchKaufman) where {T} =
Expand All @@ -90,9 +71,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
Open in desktop
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
Open in desktop
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