Defer HessenbergQ multiplication to LinearAlgebra where available

Project.toml

@@ -1,6 +1,6 @@
 name = "GenericLinearAlgebra"
 uuid = "14197337-ba66-59df-a3e3-ca00e7dcff7a"
-version = "0.4.1"
+version = "0.4.2"
 
 [deps]
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"

src/svd.jl

@@ -2,8 +2,6 @@ using LinearAlgebra
 
 import LinearAlgebra: mul!, rmul!
 
-AdjointQtype = isdefined(LinearAlgebra, :AdjointQ) ? LinearAlgebra.AdjointQ : Adjoint
-
 lmul!(G::LinearAlgebra.Givens, ::Nothing) = nothing
 rmul!(::Nothing, G::LinearAlgebra.Givens) = nothing
 
@@ -399,7 +397,7 @@ function Base.getproperty(F::BidiagonalFactorization, s::Symbol)
                 typeof(τr),
                 false,
             }(
-                'U',
+                'L',
                 factors,
                 τr,
             )
@@ -412,10 +410,10 @@ function Base.getproperty(F::BidiagonalFactorization, s::Symbol)
             return LinearAlgebra.HessenbergQ{
                 eltype(factors),
                 typeof(factors),
-                typeof(τr),
+                typeof(τl),
                 false,
             }(
-                'U',
+                'L',
                 factors,
                 τl,
             )
@@ -429,7 +427,13 @@ function Base.getproperty(F::BidiagonalFactorization, s::Symbol)
     end
 end
 
-# For now, we define a generic lmul! and rmul! for HessenbergQ here
+# LinearAlgebra's generic `HessenbergQ` multiplication (available in versions that ship it)
+# handles non-BLAS element types and the `uplo` field, so these fallbacks are only needed for
+# older versions. Defining them unconditionally would overwrite the upstream methods, which is
+# an error during precompilation.
+@static if !hasmethod(lmul!, Tuple{LinearAlgebra.HessenbergQ,AbstractVecOrMat})
+AdjointQtype = isdefined(LinearAlgebra, :AdjointQ) ? LinearAlgebra.AdjointQ : Adjoint
+
 function lmul!(Q::LinearAlgebra.HessenbergQ, B::AbstractVecOrMat)
 
     m, n = size(B, 1), size(B, 2)
@@ -530,6 +534,7 @@ function rmul!(A::AbstractMatrix, adjQ::AdjointQtype{<:Any,<:LinearAlgebra.Hesse
 
     return A
 end
+end
 
 # FIXME! Commented out for now because we currently don't produce
 # any LQPackedQ matrices.

test/svd.jl

@@ -249,10 +249,27 @@ using Test, GenericLinearAlgebra, LinearAlgebra, Quaternions, DoubleFloats
         @test diag(F.U' * B * F.Vt') ≈ F.S rtol = 5e-15
     end
 
-    @testset "Generic HessenbergQ multiplication" begin
-        A = big.(randn(10, 10))
+    # the tall case stores the right reflectors in a HessenbergQ, the wide case the left ones
+    @testset "Generic HessenbergQ multiplication: big($T), ($m, $n)" for T in (Float64, ComplexF64),
+        (m, n) in ((10, 10), (10, 6), (6, 10))
+
+        A = big.(randn(T, m, n))
+        TB = eltype(A)
         BF = GenericLinearAlgebra.bidiagonalize!(copy(A))
-        @test (BF.rightQ' * Matrix(I, size(A)...)) * BF.rightQ ≈ I
+        for Q in (BF.leftQ, BF.rightQ)
+            k = size(Q, 1)
+            QI = Q * Matrix{TB}(I, k, k)
+            @test Q' * QI ≈ I
+            B = big.(randn(T, k, 3))
+            C = big.(randn(T, 3, k))
+            @test lmul!(Q, copy(B)) ≈ QI * B
+            @test lmul!(Q', copy(B)) ≈ QI' * B
+            @test rmul!(copy(C), Q) ≈ C * QI
+            @test rmul!(copy(C), Q') ≈ C * QI'
+        end
+        Bd = zeros(TB, m, n)
+        Bd[1:min(m, n), 1:min(m, n)] = BF.bidiagonal
+        @test BF.leftQ * Bd * BF.rightQ' ≈ A
     end
 
     @testset "Issue 119" begin