Defer HessenbergQ multiplication to LinearAlgebra where available (#185)

* Defer HessenbergQ multiplication to LinearAlgebra where available

LinearAlgebra now ships generic `lmul!`/`rmul!` methods for `HessenbergQ`
with the same signatures as the ones defined here. Redefining them
overwrites the upstream methods, which is an error during precompilation,
so only define the local fallbacks when LinearAlgebra lacks them.

The `HessenbergQ`s built in `BidiagonalFactorization.getproperty` store
their reflectors below the subdiagonal, which is the `'L'` layout that
`hessenberg` produces; mark them as such instead of `'U'`. Also pass the
matching `τ` type parameter for `leftQ`.

Extend the HessenbergQ test to both factors, both bidiagonal
orientations, real and complex BigFloat, all four multiplication
directions, and the reconstruction of the original matrix.

Co-Authored-By: Claude Fable 5.1 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01Vv9VYe4nM23v8ruYwk9AQx

* bump version

---------

Co-authored-by: Claude <noreply@anthropic.com>

Author: dkarrasch

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