Skip to content

Commit 6844259

Browse files
sostocksostock
authored
Preserve Hessenberg shape when possible (#40039)
* Preserve Hessenberg shape when possible * Remove obsolete comments * Add NEWS * Support for unitful Hessenberg matrices * Mark some tests for unitful matrices as broken
1 parent 5b7f4c5 commit 6844259

File tree

4 files changed

+139
-6
lines changed

4 files changed

+139
-6
lines changed

‎NEWS.md

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -93,6 +93,7 @@ Standard library changes
9393
* On aarch64, OpenBLAS now uses an ILP64 BLAS like all other 64-bit platforms. ([#39436])
9494
* OpenBLAS is updated to 0.3.13. ([#39216])
9595
* SuiteSparse is updated to 5.8.1. ([#39455])
96+
* The shape of an `UpperHessenberg` matrix is preserved under certain arithmetic operations, e.g. when multiplying or dividing by an `UpperTriangular` matrix. ([#40039])
9697
* `cis(A)` now supports matrix arguments ([#40194]).
9798

9899
#### Markdown

‎stdlib/LinearAlgebra/src/hessenberg.jl

Lines changed: 85 additions & 6 deletions
Original file line numberDiff line numberDiff line change
@@ -94,17 +94,96 @@ Base.copy(A::Transpose{<:Any,<:UpperHessenberg}) = tril!(transpose!(similar(A.pa
9494
rmul!(H::UpperHessenberg, x::Number) = (rmul!(H.data, x); H)
9595
lmul!(x::Number, H::UpperHessenberg) = (lmul!(x, H.data); H)
9696

97-
# (future: we could also have specialized routines for UpperHessenberg * UpperTriangular)
98-
9997
fillstored!(H::UpperHessenberg, x) = (fillband!(H.data, x, -1, size(H,2)-1); H)
10098

10199
+(A::UpperHessenberg, B::UpperHessenberg) = UpperHessenberg(A.data+B.data)
102100
-(A::UpperHessenberg, B::UpperHessenberg) = UpperHessenberg(A.data-B.data)
103-
# (future: we could also have specialized routines for UpperHessenberg ± UpperTriangular)
104101

105-
# shift Hessenberg by λI
106-
+(H::UpperHessenberg, J::UniformScaling) = UpperHessenberg(H.data + J)
107-
-(J::UniformScaling, H::UpperHessenberg) = UpperHessenberg(J - H.data)
102+
for T = (:UniformScaling, :Diagonal, :Bidiagonal, :Tridiagonal, :SymTridiagonal,
103+
:UpperTriangular, :UnitUpperTriangular)
104+
for op = (:+, :-)
105+
@eval begin
106+
$op(H::UpperHessenberg, x::$T) = UpperHessenberg($op(H.data, x))
107+
$op(x::$T, H::UpperHessenberg) = UpperHessenberg($op(x, H.data))
108+
end
109+
end
110+
end
111+
112+
for T = (:Number, :UniformScaling, :Diagonal)
113+
@eval begin
114+
*(H::UpperHessenberg, x::$T) = UpperHessenberg(H.data * x)
115+
*(x::$T, H::UpperHessenberg) = UpperHessenberg(x * H.data)
116+
/(H::UpperHessenberg, x::$T) = UpperHessenberg(H.data / x)
117+
\(x::$T, H::UpperHessenberg) = UpperHessenberg(x \ H.data)
118+
end
119+
end
120+
121+
function *(H::UpperHessenberg, U::UpperOrUnitUpperTriangular)
122+
T = typeof(oneunit(eltype(H))*oneunit(eltype(U)))
123+
HH = similar(H.data, T, size(H))
124+
copyto!(HH, H)
125+
rmul!(HH, U)
126+
UpperHessenberg(HH)
127+
end
128+
function *(U::UpperOrUnitUpperTriangular, H::UpperHessenberg)
129+
T = typeof(oneunit(eltype(H))*oneunit(eltype(U)))
130+
HH = similar(H.data, T, size(H))
131+
copyto!(HH, H)
132+
lmul!(U, HH)
133+
UpperHessenberg(HH)
134+
end
135+
136+
function /(H::UpperHessenberg, U::UpperTriangular)
137+
T = typeof(oneunit(eltype(H))/oneunit(eltype(U)))
138+
HH = similar(H.data, T, size(H))
139+
copyto!(HH, H)
140+
rdiv!(HH, U)
141+
UpperHessenberg(HH)
142+
end
143+
function /(H::UpperHessenberg, U::UnitUpperTriangular)
144+
T = typeof(oneunit(eltype(H))/oneunit(eltype(U)))
145+
HH = similar(H.data, T, size(H))
146+
copyto!(HH, H)
147+
rdiv!(HH, U)
148+
UpperHessenberg(HH)
149+
end
150+
151+
function \(U::UpperTriangular, H::UpperHessenberg)
152+
T = typeof(oneunit(eltype(U))\oneunit(eltype(H)))
153+
HH = similar(H.data, T, size(H))
154+
copyto!(HH, H)
155+
ldiv!(U, HH)
156+
UpperHessenberg(HH)
157+
end
158+
function \(U::UnitUpperTriangular, H::UpperHessenberg)
159+
T = typeof(oneunit(eltype(U))\oneunit(eltype(H)))
160+
HH = similar(H.data, T, size(H))
161+
copyto!(HH, H)
162+
ldiv!(U, HH)
163+
UpperHessenberg(HH)
164+
end
165+
166+
function *(H::UpperHessenberg, B::Bidiagonal)
167+
TS = promote_op(matprod, eltype(H), eltype(B))
168+
if B.uplo == 'U'
169+
A_mul_B_td!(UpperHessenberg(zeros(TS, size(H)...)), H, B)
170+
else
171+
A_mul_B_td!(zeros(TS, size(H)...), H, B)
172+
end
173+
end
174+
function *(B::Bidiagonal, H::UpperHessenberg)
175+
TS = promote_op(matprod, eltype(B), eltype(H))
176+
if B.uplo == 'U'
177+
A_mul_B_td!(UpperHessenberg(zeros(TS, size(B)...)), B, H)
178+
else
179+
A_mul_B_td!(zeros(TS, size(B)...), B, H)
180+
end
181+
end
182+
183+
function /(H::UpperHessenberg, B::Bidiagonal)
184+
A = Base.@invoke /(H::AbstractMatrix, B::Bidiagonal)
185+
B.uplo == 'U' ? UpperHessenberg(A) : A
186+
end
108187

109188
# Solving (H+µI)x = b: we can do this in O(m²) time and O(m) memory
110189
# (in-place in x) by the RQ algorithm from:

‎stdlib/LinearAlgebra/test/hessenberg.jl

Lines changed: 52 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -4,6 +4,10 @@ module TestHessenberg
44

55
using Test, LinearAlgebra, Random
66

7+
const BASE_TEST_PATH = joinpath(Sys.BINDIR, "..", "share", "julia", "test")
8+
isdefined(Main, :Furlongs) || @eval Main include(joinpath($(BASE_TEST_PATH), "testhelpers", "Furlongs.jl"))
9+
using .Main.Furlongs
10+
711
# for tuple tests below
812
≅(x,y) = all(p -> p[1] ≈ p[2], zip(x,y))
913

@@ -55,6 +59,54 @@ let n = 10
5559
H = UpperHessenberg(Areal)
5660
@test Array(Hc + H) == Array(Hc) + Array(H)
5761
@test Array(Hc - H) == Array(Hc) - Array(H)
62+
@testset "Preserve UpperHessenberg shape (issue #39388)" begin
63+
for H = (UpperHessenberg(Areal), UpperHessenberg(Furlong.(Areal)))
64+
if eltype(H) <: Furlong
65+
A = Furlong.(rand(n,n))
66+
d = Furlong.(rand(n))
67+
dl = Furlong.(rand(n-1))
68+
du = Furlong.(rand(n-1))
69+
us = Furlong(1)*I
70+
else
71+
A = rand(n,n)
72+
d = rand(n)
73+
dl = rand(n-1)
74+
du = rand(n-1)
75+
us = 1*I
76+
end
77+
@testset "$op" for op = (+,-)
78+
for x = (us, Diagonal(d), Bidiagonal(d,dl,:U), Bidiagonal(d,dl,:L),
79+
Tridiagonal(dl,d,du), SymTridiagonal(d,dl),
80+
UpperTriangular(A), UnitUpperTriangular(A))
81+
@test op(H,x) == op(Array(H),x)
82+
@test op(x,H) == op(x,Array(H))
83+
@test op(H,x) isa UpperHessenberg
84+
@test op(x,H) isa UpperHessenberg
85+
end
86+
end
87+
A = randn(n,n)
88+
d = randn(n)
89+
dl = randn(n-1)
90+
@testset "Multiplication/division" begin
91+
for x = (5, 5I, Diagonal(d), Bidiagonal(d,dl,:U),
92+
UpperTriangular(A), UnitUpperTriangular(A))
93+
@test H*x == Array(H)*x broken = eltype(H) <: Furlong && x isa Bidiagonal
94+
@test x*H == x*Array(H) broken = eltype(H) <: Furlong && x isa Bidiagonal
95+
@test H/x == Array(H)/x broken = eltype(H) <: Furlong && x isa Union{Bidiagonal, Diagonal, UpperTriangular}
96+
@test x\H == x\Array(H) broken = eltype(H) <: Furlong && x isa Union{Bidiagonal, Diagonal, UpperTriangular}
97+
@test H*x isa UpperHessenberg broken = eltype(H) <: Furlong && x isa Bidiagonal
98+
@test x*H isa UpperHessenberg broken = eltype(H) <: Furlong && x isa Bidiagonal
99+
@test H/x isa UpperHessenberg broken = eltype(H) <: Furlong && x isa Union{Bidiagonal, Diagonal}
100+
@test x\H isa UpperHessenberg broken = eltype(H) <: Furlong && x isa Union{Bidiagonal, Diagonal}
101+
end
102+
x = Bidiagonal(d, dl, :L)
103+
@test H*x == Array(H)*x
104+
@test x*H == x*Array(H)
105+
@test H/x == Array(H)/x broken = eltype(H) <: Furlong
106+
@test_broken x\H == x\Array(H) # issue 40037
107+
end
108+
end
109+
end
58110
end
59111

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

‎test/testhelpers/Furlongs.jl

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -36,6 +36,7 @@ Base.floatmin(::Type{Furlong{p,T}}) where {p,T<:AbstractFloat} = Furlong{p}(floa
3636
Base.floatmin(::Furlong{p,T}) where {p,T<:AbstractFloat} = floatmin(Furlong{p,T})
3737
Base.floatmax(::Type{Furlong{p,T}}) where {p,T<:AbstractFloat} = Furlong{p}(floatmax(T))
3838
Base.floatmax(::Furlong{p,T}) where {p,T<:AbstractFloat} = floatmax(Furlong{p,T})
39+
Base.conj(x::Furlong{p,T}) where {p,T} = Furlong{p,T}(conj(x.val))
3940

4041
# convert Furlong exponent p to a canonical form. This
4142
# is not type stable, but it doesn't matter since it is used

0 commit comments

Comments
 (0)