From 4578357b4f7e9e57a2b29811d87feed08e49b7fb Mon Sep 17 00:00:00 2001
From: "Steven G. Johnson" <stevenj@alum.mit.edu>
Date: Tue, 29 Apr 2025 14:37:58 -0400
Subject: [PATCH 1/2] add missing methods for division of Hessenberg matrices

---
 src/hessenberg.jl  | 23 +++++++++++++++++++++++
 test/hessenberg.jl |  6 ++++++
 2 files changed, 29 insertions(+)

diff --git a/src/hessenberg.jl b/src/hessenberg.jl
index 056b97e3..442e8aa2 100644
--- a/src/hessenberg.jl
+++ b/src/hessenberg.jl
@@ -177,6 +177,29 @@ function \(U::UnitUpperTriangular, H::UpperHessenberg)
     UpperHessenberg(HH)
 end
 
+function (\)(H::Union{UpperHessenberg,AdjOrTrans{<:Any,<:UpperHessenberg}}, B::AbstractVecOrMat)
+    TFB = typeof(oneunit(eltype(H)) \ oneunit(eltype(B)))
+    return ldiv!(H, copy_similar(B, TFB))
+end
+
+function (/)(B::AbstractMatrix, H::Union{UpperHessenberg,AdjOrTrans{<:Any,<:UpperHessenberg}})
+    TFB = typeof(oneunit(eltype(B)) / oneunit(eltype(H)))
+    return rdiv!(copy_similar(B, TFB), H)
+end
+
+ldiv!(H::AdjOrTrans{<:Any,<:UpperHessenberg}, B::AbstractVecOrMat) =
+    (rdiv!(wrapperop(H)(B), parent(H)); B)
+rdiv!(B::AbstractVecOrMat, H::AdjOrTrans{<:Any,<:UpperHessenberg}) =
+    (ldiv!(parent(H), wrapperop(H)(B)); B)
+
+# fix method ambiguities for right division, from adjtrans.jl:
+/(u::AdjointAbsVec, A::UpperHessenberg) = adjoint(adjoint(A) \ u.parent)
+/(u::TransposeAbsVec, A::UpperHessenberg) = transpose(transpose(A) \ u.parent)
+/(u::AdjointAbsVec, A::Adjoint{<:Any,<:UpperHessenberg}) = adjoint(adjoint(A) \ u.parent)
+/(u::TransposeAbsVec, A::Transpose{<:Any,<:UpperHessenberg}) = transpose(transpose(A) \ u.parent)
+/(u::AdjointAbsVec, A::Transpose{<:Any,<:UpperHessenberg}) = adjoint(conj(A.parent) \ u.parent) # technically should be adjoint(copy(adjoint(copy(A))) \ u.parent)
+/(u::TransposeAbsVec, A::Adjoint{<:Any,<:UpperHessenberg}) = transpose(conj(A.parent) \ u.parent)
+
 # 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:
 #
diff --git a/test/hessenberg.jl b/test/hessenberg.jl
index 91f6d155..bf5d8785 100644
--- a/test/hessenberg.jl
+++ b/test/hessenberg.jl
@@ -66,6 +66,12 @@ let n = 10
         H = UpperHessenberg(Areal)
         @test Array(Hc + H) == Array(Hc) + Array(H)
         @test Array(Hc - H) == Array(Hc) - Array(H)
+        @testset "ldiv and rdiv" begin
+            for b in (b_, B_), H in (H, Hc, H', Hc', transpose(Hc))
+                @test H * (H \ b) ≈ b
+                @test (b' / H) * H ≈ (Matrix(b') / H) * H ≈ b'
+            end
+        end
         @testset "Preserve UpperHessenberg shape (issue #39388)" begin
             H = UpperHessenberg(Areal)
             A = rand(n,n)

From 484467044b26c992078e9d4077a381fa45a4db40 Mon Sep 17 00:00:00 2001
From: "Steven G. Johnson" <stevenj@alum.mit.edu>
Date: Wed, 30 Apr 2025 12:32:18 -0400
Subject: [PATCH 2/2] Update test/hessenberg.jl

Co-authored-by: Jishnu Bhattacharya <jishnub.github@gmail.com>
---
 test/hessenberg.jl | 1 +
 1 file changed, 1 insertion(+)

diff --git a/test/hessenberg.jl b/test/hessenberg.jl
index bf5d8785..69a91020 100644
--- a/test/hessenberg.jl
+++ b/test/hessenberg.jl
@@ -70,6 +70,7 @@ let n = 10
             for b in (b_, B_), H in (H, Hc, H', Hc', transpose(Hc))
                 @test H * (H \ b) ≈ b
                 @test (b' / H) * H ≈ (Matrix(b') / H) * H ≈ b'
+                @test (transpose(b) / H) * H ≈ (Matrix(transpose(b)) / H) * H ≈ transpose(b)
             end
         end
         @testset "Preserve UpperHessenberg shape (issue #39388)" begin