Missing method in `/(::Adjoint{ComplexF64}, Hessenberg{Float64})`

[martinholters]

julia> A = rand(5,5);

julia> b = rand(5);

julia> H = hessenberg(A);

julia> b' / H # adjoint of real vector is ok
1×5 Matrix{Float64}:
 -7.66339  3.09448  7.69176  0.771431  -0.683927

julia> Matrix(complex(b)') / H # complex single-row matrix is ok
1×5 Matrix{ComplexF64}:
 -7.66339+0.0im  3.09448+0.0im  7.69176-0.0im  0.771431+0.0im  -0.683927+0.0im

julia> complex(b)' / H # adjoint of complex vector fails
ERROR: MethodError: no method matching rmul!(::Matrix{ComplexF64}, ::LinearAlgebra.HessenbergQ{Float64, Matrix{Float64}, Vector{Float64}, false})

Is there a good reason why the last call shouldn't give the same result as the one before? (Or maybe even better, an adjoint vector rather than a single-row matrix, also in the real case.)

Noticed while working on #1015 and as a by-product making complex.(b') return an Adjoint (just as broadcast(complex, b') does already) and thus breaking
https://github.com/JuliaLang/julia/blob/958f647d8e2d7d93b9d158f283e0891e8e9cf0b0/stdlib/LinearAlgebra/test/hessenberg.jl#L181

[dkarrasch]

The underlying reason is that all operations with HessenbergQ call a LAPACK routine, which requires the eltypes to be the same. The last but one call works because of https://github.com/JuliaLang/julia/blob/d1be33d4bc9fe88bbcee2055ad3b0fb943673122/stdlib/LinearAlgebra/src/factorization.jl#L154-L158

For this reason, I think we should simply remove the complex part of the test, or add an explicit Matrix call.

[martinholters]

Alright. Looking at that method, the VecOrMat does not make much sense, does it? A column vector for B should probably not work there. And indeed, it does not:

julia> complex(b) / H
ERROR: MethodError: no method matching rdiv!(::Vector{Float64}, ::Hessenberg{Float64, UpperHessenberg{Float64, Matrix{Float64}}, Matrix{Float64}, Vector{Float64}, Bool})

Stacktrace:
 [1] /(B::Vector{ComplexF64}, F::Hessenberg{Float64, UpperHessenberg{Float64, Matrix{Float64}}, Matrix{Float64}, Vector{Float64}, Bool})
   @ LinearAlgebra /usr/share/julia/stdlib/v1.11/LinearAlgebra/src/factorization.jl:156

So I guess it should be just B::Matrix{Complex{T}} or - and that might actually fix this issue - B::Union{Matrix{Complex{T}},Adjoint{Complex{T}, Vector{Complex{T}}}}. Would that make sense?

[dkarrasch]

The method goes back JuliaLang/julia#31853, but that detail is likely only a copy-paste oversight, so we could remove the Vec part and run pkgeval. OTOH, this shows that the test was added on purpose, however, assuming that broadcasting would yield a plain VecOrMat. I'm not sure we can add the Adjoint in the signature, because the reinterpret trick doesn't work:

julia> reinterpret(complex(b)', Float64)
ERROR: MethodError: no method matching reinterpret(::Adjoint{ComplexF64, Vector{ComplexF64}}, ::Type{Float64})

[martinholters]

You typoed that one (reversed the arguments):

julia> reinterpret(Float64, complex(b)')
2×5 reinterpret(Float64, adjoint(::Vector{ComplexF64})):
  0.111025   0.431175   0.848844   0.609683   0.622186
 -0.0       -0.0       -0.0       -0.0       -0.0

So can do:

julia> reinterpret(ComplexF64, reinterpret(Float64, complex(b)') / H)
1×5 reinterpret(ComplexF64, ::Matrix{Float64}):
 -0.108541+0.0im  0.413376+0.0im  0.795175-0.0im  -0.34378-0.0im  0.574402-0.0im

[dkarrasch]

🤦 Thanks for pointing out. Cool, that seems like we should be able to change the signature as you proposed, right?

[martinholters]

Needs to resolve an ambiguity with
https://github.com/JuliaLang/julia/blob/d1be33d4bc9fe88bbcee2055ad3b0fb943673122/stdlib/LinearAlgebra/src/factorization.jl#L177
but seems to work then. I'll also look at Transpose which can probably be treated the same and open a PR once I've resolved that.

[dkarrasch]

Hihi, I'm also working on it currently. I found a few typos in the / methods below the "don't use reinterpret" comment, and realized they are not covered by tests.