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

Preserve Hessenberg shape when possible #40039

Merged
merged 7 commits into from
Apr 11, 2021
Merged

Preserve Hessenberg shape when possible #40039

Show file tree
Hide file tree
Changes from 1 commit
Commits
Show all changes
7 commits
Select commit Hold shift + click to select a range
7c7b829
Preserve Hessenberg shape when possible
sostock Mar 15, 2021
8836f29
Remove obsolete comments
sostock Apr 7, 2021
6ccf862
Add NEWS
sostock Apr 7, 2021
fc6da6f
Merge branch 'master' into preserve_hessenberg
sostock Apr 7, 2021
f57d816
Support for unitful Hessenberg matrices
sostock Apr 10, 2021
4bdfcdf
Merge branch 'master' into preserve_hessenberg
sostock Apr 10, 2021
9d02340
Mark some tests for unitful matrices as broken
sostock Apr 10, 2021
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
88 changes: 85 additions & 3 deletions stdlib/LinearAlgebra/src/hessenberg.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -102,9 +102,91 @@ fillstored!(H::UpperHessenberg, x) = (fillband!(H.data, x, -1, size(H,2)-1); H)
-(A::UpperHessenberg, B::UpperHessenberg) = UpperHessenberg(A.data-B.data)
# (future: we could also have specialized routines for UpperHessenberg ± UpperTriangular)

# shift Hessenberg by λI
+(H::UpperHessenberg, J::UniformScaling) = UpperHessenberg(H.data + J)
-(J::UniformScaling, H::UpperHessenberg) = UpperHessenberg(J - H.data)
for T = (:UniformScaling, :Diagonal, :Bidiagonal, :Tridiagonal, :SymTridiagonal,
:UpperTriangular, :UnitUpperTriangular)
for op = (:+, :-)
@eval begin
$op(H::UpperHessenberg, x::$T) = UpperHessenberg($op(H.data, x))
$op(x::$T, H::UpperHessenberg) = UpperHessenberg($op(x, H.data))
end
end
end

for T = (:Number, :UniformScaling, :Diagonal)
@eval begin
*(H::UpperHessenberg, x::$T) = UpperHessenberg(H.data * x)
*(x::$T, H::UpperHessenberg) = UpperHessenberg(x * H.data)
/(H::UpperHessenberg, x::$T) = UpperHessenberg(H.data / x)
\(x::$T, H::UpperHessenberg) = UpperHessenberg(x \ H.data)
end
end

function *(H::UpperHessenberg, U::UpperOrUnitUpperTriangular)
T = _inner_type_promotion(H,U)
HH = similar(H.data, T, size(H))
copyto!(HH, H)
rmul!(HH, convert(AbstractArray{T}, U))
UpperHessenberg(HH)
end
function *(U::UpperOrUnitUpperTriangular, H::UpperHessenberg)
T = _inner_type_promotion(U,H)
HH = similar(H.data, T, size(H))
copyto!(HH, H)
lmul!(convert(AbstractArray{T}, U), HH)
UpperHessenberg(HH)
end

function /(H::UpperHessenberg, U::UpperTriangular)
T = typeof((zero(eltype(H))*zero(eltype(U)) + zero(eltype(H))*zero(eltype(U)))/one(eltype(H)))
HH = similar(H.data, T, size(H))
copyto!(HH, H)
rdiv!(HH, convert(AbstractArray{T}, U))
UpperHessenberg(HH)
end
function /(H::UpperHessenberg, U::UnitUpperTriangular)
T = _inner_type_promotion(H,U)
HH = similar(H.data, T, size(H))
copyto!(HH, H)
rdiv!(HH, convert(AbstractArray{T}, U))
UpperHessenberg(HH)
end

function \(U::UpperTriangular, H::UpperHessenberg)
T = typeof((zero(eltype(U))*zero(eltype(H)) + zero(eltype(U))*zero(eltype(H)))/one(eltype(U)))
HH = similar(H.data, T, size(H))
copyto!(HH, H)
ldiv!(convert(AbstractArray{T}, U), HH)
UpperHessenberg(HH)
end
function \(U::UnitUpperTriangular, H::UpperHessenberg)
T = _inner_type_promotion(U,H)
HH = similar(H.data, T, size(H))
copyto!(HH, H)
ldiv!(convert(AbstractArray{T}, U), HH)
UpperHessenberg(HH)
end

function *(H::UpperHessenberg, B::Bidiagonal)
TS = promote_op(matprod, eltype(H), eltype(B))
if B.uplo == 'U'
A_mul_B_td!(UpperHessenberg(zeros(TS, size(H)...)), H, B)
else
A_mul_B_td!(zeros(TS, size(H)...), H, B)
end
end
function *(B::Bidiagonal, H::UpperHessenberg)
TS = promote_op(matprod, eltype(B), eltype(H))
if B.uplo == 'U'
A_mul_B_td!(UpperHessenberg(zeros(TS, size(B)...)), B, H)
else
A_mul_B_td!(zeros(TS, size(B)...), B, H)
end
end

function /(H::UpperHessenberg, B::Bidiagonal)
A = Base.@invoke /(H::AbstractMatrix, B::Bidiagonal)
B.uplo == 'U' ? UpperHessenberg(A) : A
end

# Solving (H+µI)x = b: we can do this in O(m²) time and O(m) memory
# (in-place in x) by the RQ algorithm from:
Expand Down
34 changes: 34 additions & 0 deletions stdlib/LinearAlgebra/test/hessenberg.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -55,6 +55,40 @@ let n = 10
H = UpperHessenberg(Areal)
@test Array(Hc + H) == Array(Hc) + Array(H)
@test Array(Hc - H) == Array(Hc) - Array(H)
@testset "Preserve UpperHessenberg shape (issue #39388)" begin
A = rand(n,n)
d = rand(n)
dl = rand(n-1)
du = rand(n-1)
@testset "$op" for op = (+,-)
for x = (I, Diagonal(d), Bidiagonal(d,dl,:U), Bidiagonal(d,dl,:L),
Tridiagonal(dl,d,du), SymTridiagonal(d,dl),
UpperTriangular(A), UnitUpperTriangular(A))
@test op(H,x) == op(Array(H),x)
@test op(x,H) == op(x,Array(H))
@test op(H,x) isa UpperHessenberg
@test op(x,H) isa UpperHessenberg
end
end
@testset "Multiplication/division" begin
for x = (5, 5I, Diagonal(d), Bidiagonal(d,dl,:U),
UpperTriangular(A), UnitUpperTriangular(A))
@test H*x == Array(H)*x
@test x*H == x*Array(H)
@test H/x == Array(H)/x
@test x\H == x\Array(H)
@test H*x isa UpperHessenberg
@test x*H isa UpperHessenberg
@test H/x isa UpperHessenberg
@test x\H isa UpperHessenberg
end
x = Bidiagonal(d, dl, :L)
@test H*x == Array(H)*x
@test x*H == x*Array(H)
@test H/x == Array(H)/x
@test_broken x\H == x\Array(H) # issue 40037
end
end
end

@testset for eltya in (Float32, Float64, ComplexF32, ComplexF64, Int), herm in (false, true)
Expand Down