Preserve Hessenberg shape when possible

NEWS.md

@@ -93,6 +93,7 @@ Standard library changes
 * On aarch64, OpenBLAS now uses an ILP64 BLAS like all other 64-bit platforms. ([#39436])
 * OpenBLAS is updated to 0.3.13. ([#39216])
 * SuiteSparse is updated to 5.8.1. ([#39455])
+* The shape of an `UpperHessenberg` matrix is preserved under certain arithmetic operations, e.g. when multiplying or dividing by an `UpperTriangular` matrix. ([#40039])
 * `cis(A)` now supports matrix arguments ([#40194]).
 
 #### Markdown

stdlib/LinearAlgebra/src/hessenberg.jl

@@ -94,17 +94,96 @@ Base.copy(A::Transpose{<:Any,<:UpperHessenberg}) = tril!(transpose!(similar(A.pa
 rmul!(H::UpperHessenberg, x::Number) = (rmul!(H.data, x); H)
 lmul!(x::Number, H::UpperHessenberg) = (lmul!(x, H.data); H)
 
-# (future: we could also have specialized routines for UpperHessenberg * UpperTriangular)
-
 fillstored!(H::UpperHessenberg, x) = (fillband!(H.data, x, -1, size(H,2)-1); H)
 
 +(A::UpperHessenberg, B::UpperHessenberg) = UpperHessenberg(A.data+B.data)
 -(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 = typeof(oneunit(eltype(H))*oneunit(eltype(U)))
+    HH = similar(H.data, T, size(H))
+    copyto!(HH, H)
+    rmul!(HH, U)
+    UpperHessenberg(HH)
+end
+function *(U::UpperOrUnitUpperTriangular, H::UpperHessenberg)
+    T = typeof(oneunit(eltype(H))*oneunit(eltype(U)))
+    HH = similar(H.data, T, size(H))
+    copyto!(HH, H)
+    lmul!(U, HH)
+    UpperHessenberg(HH)
+end
+
+function /(H::UpperHessenberg, U::UpperTriangular)
+    T = typeof(oneunit(eltype(H))/oneunit(eltype(U)))
+    HH = similar(H.data, T, size(H))
+    copyto!(HH, H)
+    rdiv!(HH, U)
+    UpperHessenberg(HH)
+end
+function /(H::UpperHessenberg, U::UnitUpperTriangular)
+    T = typeof(oneunit(eltype(H))/oneunit(eltype(U)))
+    HH = similar(H.data, T, size(H))
+    copyto!(HH, H)
+    rdiv!(HH, U)
+    UpperHessenberg(HH)
+end
+
+function \(U::UpperTriangular, H::UpperHessenberg)
+    T = typeof(oneunit(eltype(U))\oneunit(eltype(H)))
+    HH = similar(H.data, T, size(H))
+    copyto!(HH, H)
+    ldiv!(U, HH)
+    UpperHessenberg(HH)
+end
+function \(U::UnitUpperTriangular, H::UpperHessenberg)
+    T = typeof(oneunit(eltype(U))\oneunit(eltype(H)))
+    HH = similar(H.data, T, size(H))
+    copyto!(HH, H)
+    ldiv!(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:

stdlib/LinearAlgebra/test/hessenberg.jl

@@ -4,6 +4,10 @@ module TestHessenberg
 
 using Test, LinearAlgebra, Random
 
+const BASE_TEST_PATH = joinpath(Sys.BINDIR, "..", "share", "julia", "test")
+isdefined(Main, :Furlongs) || @eval Main include(joinpath($(BASE_TEST_PATH), "testhelpers", "Furlongs.jl"))
+using .Main.Furlongs
+
 # for tuple tests below
 ≅(x,y) = all(p -> p[1] ≈ p[2], zip(x,y))
 
@@ -55,6 +59,54 @@ 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
+            for H = (UpperHessenberg(Areal), UpperHessenberg(Furlong.(Areal)))
+                if eltype(H) <: Furlong
+                    A = Furlong.(rand(n,n))
+                    d = Furlong.(rand(n))
+                    dl = Furlong.(rand(n-1))
+                    du = Furlong.(rand(n-1))
+                    us = Furlong(1)*I
+                else
+                    A = rand(n,n)
+                    d = rand(n)
+                    dl = rand(n-1)
+                    du = rand(n-1)
+                    us = 1*I
+                end
+                @testset "$op" for op = (+,-)
+                    for x = (us, 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
+                A = randn(n,n)
+                d = randn(n)
+                dl = randn(n-1)
+                @testset "Multiplication/division" begin
+                    for x = (5, 5I, Diagonal(d), Bidiagonal(d,dl,:U),
+                             UpperTriangular(A), UnitUpperTriangular(A))
+                        @test H*x == Array(H)*x broken = eltype(H) <: Furlong && x isa Bidiagonal
+                        @test x*H == x*Array(H) broken = eltype(H) <: Furlong && x isa Bidiagonal
+                        @test H/x == Array(H)/x broken = eltype(H) <: Furlong && x isa Union{Bidiagonal, Diagonal, UpperTriangular}
+                        @test x\H == x\Array(H) broken = eltype(H) <: Furlong && x isa Union{Bidiagonal, Diagonal, UpperTriangular}
+                        @test H*x isa UpperHessenberg broken = eltype(H) <: Furlong && x isa Bidiagonal
+                        @test x*H isa UpperHessenberg broken = eltype(H) <: Furlong && x isa Bidiagonal
+                        @test H/x isa UpperHessenberg broken = eltype(H) <: Furlong && x isa Union{Bidiagonal, Diagonal}
+                        @test x\H isa UpperHessenberg broken = eltype(H) <: Furlong && x isa Union{Bidiagonal, Diagonal}
+                    end
+                    x = Bidiagonal(d, dl, :L)
+                    @test H*x == Array(H)*x
+                    @test x*H == x*Array(H)
+                    @test H/x == Array(H)/x broken = eltype(H) <: Furlong
+                    @test_broken x\H == x\Array(H) # issue 40037
+                end
+            end
+        end
     end
 
     @testset for eltya in (Float32, Float64, ComplexF32, ComplexF64, Int), herm in (false, true)

test/testhelpers/Furlongs.jl

@@ -36,6 +36,7 @@ Base.floatmin(::Type{Furlong{p,T}}) where {p,T<:AbstractFloat} = Furlong{p}(floa
 Base.floatmin(::Furlong{p,T}) where {p,T<:AbstractFloat} = floatmin(Furlong{p,T})
 Base.floatmax(::Type{Furlong{p,T}}) where {p,T<:AbstractFloat} = Furlong{p}(floatmax(T))
 Base.floatmax(::Furlong{p,T}) where {p,T<:AbstractFloat} = floatmax(Furlong{p,T})
+Base.conj(x::Furlong{p,T}) where {p,T} = Furlong{p,T}(conj(x.val))
 
 # convert Furlong exponent p to a canonical form.  This
 # is not type stable, but it doesn't matter since it is used