`\(::Bidiagonal, ::T)` errors for some special matrix types `T`

[sostock]

julia> H = UpperHessenberg(rand(5,5));

julia> U = Bidiagonal(rand(5), rand(4), :U);

julia> L = Bidiagonal(rand(5), rand(4), :L);

julia> U\H
5×5 UpperHessenberg{Float64, Matrix{Float64}}:
 -2.57965  0.696361   1.53155    0.761887    2.2419
  1.1582   0.203637   0.560071   0.0709049  -0.669847
   ⋅       3.97226   -0.0898555  0.164623    5.72021
   ⋅        ⋅         1.02833    0.446765   -0.473328
   ⋅        ⋅          ⋅         0.798988    1.11819

julia> L\H
ERROR: ArgumentError: cannot set index in the lower triangular part (3, 1) of an UpperHessenberg matrix to a nonzero value (-2.050434125517565)
[...]

While U\H is actually of UpperHessenberg shape, L\H is not and must return a full Matrix.

Edit: It also happens for UpperTriangular instead of UpperHessenberg. For UnitUpperTriangular, it errors for both kinds of Bidiagonals.

[inkydragon]

Test code:

using LinearAlgebra, Random
Random.seed!(1234321)

H = UpperHessenberg(rand(5,5))
U = Bidiagonal(rand(5), rand(4), :U)
L = Bidiagonal(rand(5), rand(4), :L)

U\H
L\H  # Error

( Array(U) \ Array(H) ) ≈ Array(U\H)
Array(L) \ Array(H)

Julia REPL output

``` julia> using LinearAlgebra, Random julia> Random.seed!(1234321) MersenneTwister(1234321) julia> julia> H = UpperHessenberg(rand(5,5)) 5×5 UpperHessenberg{Float64, Matrix{Float64}}: 0.0944218 0.87151 0.735234 0.843184 0.969904 0.936611 0.041553 0.212451 0.143325 0.185209 . 0.968779 0.845895 0.914686 0.460513 . . 0.0171136 0.715679 0.957225 . . . 0.757121 0.958056 julia> U = Bidiagonal(rand(5), rand(4), :U) 5×5 Bidiagonal{Float64, Vector{Float64}}: 0.484397 0.963961 . . . . 0.805794 0.387011 . . . . 0.372153 0.152221 . . . . 0.601021 0.538903 . . . . 0.804712 julia> L = Bidiagonal(rand(5), rand(4), :L) 5×5 Bidiagonal{Float64, Vector{Float64}}: 0.558128 . . . . 0.427103 0.221613 . . . . 0.465516 0.404656 . . . . 0.45674 0.0988376 . . . . 0.0239988 0.0703304 julia> julia> U\H 5×5 UpperHessenberg{Float64, Matrix{Float64}}: -2.11817 4.1846 3.15448 3.60014 2.52229 1.16235 -1.1987 -0.82243 -0.934388 -0.261304 . 2.60317 2.26133 2.31582 1.02262 . . 0.0284743 0.347155 0.525158 . . . 0.940859 1.19056 julia> L\H # Error ERROR: ArgumentError: cannot set index in the lower triangular part (3, 1) of an UpperHessenberg matrix to a nonzero value (-4.4869041057918695) Stacktrace: [1] setindex! @ C:\buildbot\worker\package_win64\build\usr\share\julia\stdlib\v1.6\LinearAlgebra\src\hessenberg.jl:78 [inlined] [2] _setindex! @ .\abstractarray.jl:1302 [inlined] [3] setindex! @ .\abstractarray.jl:1267 [inlined] [4] copyto!(dest::UpperHessenberg{Float64, Matrix{Float64}}, dstart::Int64, src::Vector{Float64}, sstart::Int64, n::Int64) @ Base .\abstractarray.jl:1018 [5] ldiv!(A::Bidiagonal{Float64, Vector{Float64}}, B::UpperHessenberg{Float64, Matrix{Float64}}) @ LinearAlgebra C:\buildbot\worker\package_win64\build\usr\share\julia\stdlib\v1.6\LinearAlgebra\src\bidiag.jl:782 [6] \(A::Bidiagonal{Float64, Vector{Float64}}, B::UpperHessenberg{Float64, Matrix{Float64}}) @ LinearAlgebra C:\buildbot\worker\package_win64\build\usr\share\julia\stdlib\v1.6\LinearAlgebra\src\bidiag.jl:864 [7] top-level scope @ REPL[7]:1 julia> julia> ( Array(U) \ Array(H) ) ≈ Array(U\H) true julia> Array(L) \ Array(H) 5×5 Matrix{Float64}: 0.169176 1.56149 1.31732 1.51073 1.73778 3.9003 -2.82187 -1.58015 -2.26482 -2.5134 -4.4869 5.64036 3.90821 4.86586 4.02945 20.7345 -26.0647 -17.8871 -15.2447 -8.9357 -7.07521 8.89405 6.10361 15.9672 16.6713 ```

---

julia> methods(\, (Bidiagonal, UpperHessenberg))
# 2 methods for generic function "\":
[1] \(A::Bidiagonal{var"#s814", V} where {var"#s814"<:Number, V<:AbstractVector{var"#s814"}}, B::AbstractVecOrMat{var"#s813"} where var"#s813"<:Number) in LinearAlgebra at C:\Users\woclass\AppData\Local\Programs\Julia-1.6.3\share\julia\stdlib\v1.6\LinearAlgebra\src\bidiag.jl:861
[2] \(A::Bidiagonal, B::AbstractVecOrMat{T} where T) in LinearAlgebra at C:\Users\woclass\AppData\Local\Programs\Julia-1.6.3\share\julia\stdlib\v1.6\LinearAlgebra\src\bidiag.jl:866

https://github.com/JuliaLang/julia/blob/ae8452a9e0b973991c30f27beb2201db1b0ea0d3/stdlib/LinearAlgebra/src/bidiag.jl#L861-L866

Not all return values are of the same type as B.

https://github.com/JuliaLang/julia/blob/ae8452a9e0b973991c30f27beb2201db1b0ea0d3/stdlib/LinearAlgebra/src/bidiag.jl#L4-L7

Maybe we should move uplo (upper/lower) to type parameter, or add 2 subtypes: BidiagonalU, BidiagonalL

[andreasnoack]

Did this come up in a real use case or are you systematically going through the possible combinations of special matrices? My gut feeling is that splitting Bidiagonal in two isn't worth the overhead associated with introducing new types.

[sostock]

Did this come up in a real use case

I don’t think so, but I think it’s worth fixing regardless. We don’t need to introduce new types to fix this, but it would lead to more type stability since there are quite a few methods where different types are returned for different uplos, for example when multiplying a Bidiagonal and an AbstractTriangular.

[dkarrasch]

For the Hessenberg part of the OP, this should be fixed by #43779.