From 4cdb07946f755ad9c157300d66709d39b43f0b09 Mon Sep 17 00:00:00 2001 From: Jishnu Bhattacharya Date: Mon, 12 May 2025 12:22:25 +0530 Subject: [PATCH] Resolve Hessenberg ambiguity --- src/hessenberg.jl | 23 +++++++++++++++-------- 1 file changed, 15 insertions(+), 8 deletions(-) diff --git a/src/hessenberg.jl b/src/hessenberg.jl index 297ea737..d705c00e 100644 --- a/src/hessenberg.jl +++ b/src/hessenberg.jl @@ -177,28 +177,35 @@ function \(U::UnitUpperTriangular, H::UpperHessenberg) UpperHessenberg(HH) end -function (\)(H::Union{UpperHessenberg,AdjOrTrans{<:Any,<:UpperHessenberg}}, B::AbstractVecOrMat) +AdjUpperHessenberg{T,S<:UpperHessenberg{T}} = Adjoint{T, S} +TransUpperHessenberg{T,S<:UpperHessenberg{T}} = Transpose{T, S} +AdjOrTransUpperHessenberg{T,S<:UpperHessenberg{T}} = AdjOrTrans{T, S} + +function (\)(H::Union{UpperHessenberg,AdjOrTransUpperHessenberg}, 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}}) +(/)(B::AbstractMatrix, H::UpperHessenberg) = _rdiv(B, H) +(/)(B::AbstractMatrix, H::AdjUpperHessenberg) = _rdiv(B, H) +(/)(B::AbstractMatrix, H::TransUpperHessenberg) = _rdiv(B, H) +function _rdiv(B, H) TFB = typeof(oneunit(eltype(B)) / oneunit(eltype(H))) return rdiv!(copy_similar(B, TFB), H) end -ldiv!(H::AdjOrTrans{<:Any,<:UpperHessenberg}, B::AbstractVecOrMat) = +ldiv!(H::AdjOrTransUpperHessenberg, B::AbstractVecOrMat) = (rdiv!(wrapperop(H)(B), parent(H)); B) -rdiv!(B::AbstractVecOrMat, H::AdjOrTrans{<:Any,<:UpperHessenberg}) = +rdiv!(B::AbstractVecOrMat, H::AdjOrTransUpperHessenberg) = (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) +/(u::AdjointAbsVec, A::AdjUpperHessenberg) = adjoint(adjoint(A) \ u.parent) +/(u::TransposeAbsVec, A::TransUpperHessenberg) = transpose(transpose(A) \ u.parent) +/(u::AdjointAbsVec, A::TransUpperHessenberg) = adjoint(conj(A.parent) \ u.parent) # technically should be adjoint(copy(adjoint(copy(A))) \ u.parent) +/(u::TransposeAbsVec, A::AdjUpperHessenberg) = 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: