fix LAPACK.ormlq! and postmultiplication with LQPackedQ (#23803)

* Fix size on LQPackedQ and expand related tests' comments.

* Fix methods for postmultiplication with / right-application of an LQPackedQ.

base/linalg/lapack.jl

@@ -2600,13 +2600,13 @@ for (orglq, orgqr, orgql, orgrq, ormlq, ormqr, ormql, ormrq, gemqrt, elty) in
             chkside(side)
             chkstride1(A, C, tau)
             m,n = ndims(C) == 2 ? size(C) : (size(C, 1), 1)
-            mA, nA  = size(A)
+            nA = size(A, 2)
             k   = length(tau)
             if side == 'L' && m != nA
                 throw(DimensionMismatch("for a left-sided multiplication, the first dimension of C, $m, must equal the second dimension of A, $nA"))
             end
-            if side == 'R' && n != mA
-                throw(DimensionMismatch("for a right-sided multiplication, the second dimension of C, $n, must equal the first dimension of A, $mA"))
+            if side == 'R' && n != nA
+                throw(DimensionMismatch("for a right-sided multiplication, the second dimension of C, $n, must equal the second dimension of A, $nA"))
             end
             if side == 'L' && k > m
                 throw(DimensionMismatch("invalid number of reflectors: k = $k should be <= m = $m"))

base/linalg/lq.jl

@@ -164,39 +164,60 @@ for (f1, f2) in ((:A_mul_Bc, :A_mul_B!),
     end
 end
 
-### AQ
-A_mul_B!(A::StridedMatrix{T}, B::LQPackedQ{T}) where {T<:BlasFloat} = LAPACK.ormlq!('R', 'N', B.factors, B.τ, A)
-function *(A::StridedMatrix{TA}, B::LQPackedQ{TB}) where {TA,TB}
-    TAB = promote_type(TA,TB)
-    if size(B.factors,2) == size(A,2)
-        A_mul_B!(copy_oftype(A, TAB),convert(AbstractMatrix{TAB},B))
-    elseif size(B.factors,1) == size(A,2)
-        A_mul_B!( [A zeros(TAB, size(A,1), size(B.factors,2)-size(B.factors,1))], convert(AbstractMatrix{TAB},B))
+# in-place right-application of LQPackedQs
+# these methods require that the applied-to matrix's (A's) number of columns
+# match the number of columns (nQ) of the LQPackedQ (Q) (necessary for in-place
+# operation, and the underlying LAPACK routine (ormlq) treats the implicit Q
+# as its (nQ-by-nQ) square form)
+A_mul_B!(A::StridedMatrix{T}, B::LQPackedQ{T}) where {T<:BlasFloat} =
+    LAPACK.ormlq!('R', 'N', B.factors, B.τ, A)
+A_mul_Bc!(A::StridedMatrix{T}, B::LQPackedQ{T}) where {T<:BlasReal} =
+    LAPACK.ormlq!('R', 'T', B.factors, B.τ, A)
+A_mul_Bc!(A::StridedMatrix{T}, B::LQPackedQ{T}) where {T<:BlasComplex} =
+    LAPACK.ormlq!('R', 'C', B.factors, B.τ, A)
+
+# out-of-place right-application of LQPackedQs
+# unlike their in-place equivalents, these methods: (1) check whether the applied-to
+# matrix's (A's) appropriate dimension (columns for A_*, rows for Ac_*) matches the
+# number of columns (nQ) of the LQPackedQ (Q), as the underlying LAPACK routine (ormlq)
+# treats the implicit Q as its (nQ-by-nQ) square form; and (2) if the preceding dimensions
+# do not match, these methods check whether the appropriate dimension of A instead matches
+# the number of rows of the matrix of which Q is a factor (i.e. size(Q.factors, 1)),
+# and if so zero-extends A as necessary for check (1) to pass (if possible).
+*(A::StridedVecOrMat, Q::LQPackedQ) = _A_mul_Bq(A_mul_B!, A, Q)
+A_mul_Bc(A::StridedVecOrMat, Q::LQPackedQ) = _A_mul_Bq(A_mul_Bc!, A, Q)
+function _A_mul_Bq(A_mul_Bop!::FT, A::StridedVecOrMat, Q::LQPackedQ) where FT<:Function
+    TR = promote_type(eltype(A), eltype(Q))
+    if size(A, 2) == size(Q.factors, 2)
+        C = copy_oftype(A, TR)
+    elseif size(A, 2) == size(Q.factors, 1)
+        C = zeros(TR, size(A, 1), size(Q.factors, 2))
+        copy!(C, 1, A, 1, length(A))
     else
-        throw(DimensionMismatch("second dimension of A, $(size(A,2)), must equal one of the dimensions of B, $(size(B))"))
+        _rightappdimmismatch("columns")
     end
-end
-
-### AQc
-A_mul_Bc!(A::StridedMatrix{T}, B::LQPackedQ{T}) where {T<:BlasReal}    = LAPACK.ormlq!('R','T',B.factors,B.τ,A)
-A_mul_Bc!(A::StridedMatrix{T}, B::LQPackedQ{T}) where {T<:BlasComplex} = LAPACK.ormlq!('R','C',B.factors,B.τ,A)
-function A_mul_Bc(A::StridedVecOrMat{TA}, B::LQPackedQ{TB}) where {TA<:Number,TB<:Number}
-    TAB = promote_type(TA,TB)
-    A_mul_Bc!(copy_oftype(A, TAB), convert(AbstractMatrix{TAB},(B)))
-end
-
-### AcQ/AcQc
-for (f1, f2) in ((:Ac_mul_B, :A_mul_B!),
-                 (:Ac_mul_Bc, :A_mul_Bc!))
-    @eval begin
-        function ($f1)(A::StridedMatrix, B::LQPackedQ)
-            TAB = promote_type(eltype(A), eltype(B))
-            AA = similar(A, TAB, (size(A, 2), size(A, 1)))
-            adjoint!(AA, A)
-            return ($f2)(AA, B)
-        end
+    return A_mul_Bop!(C, convert(AbstractMatrix{TR}, Q))
+end
+Ac_mul_B(A::StridedMatrix, Q::LQPackedQ) = _Ac_mul_Bq(A_mul_B!, A, Q)
+Ac_mul_Bc(A::StridedMatrix, Q::LQPackedQ) = _Ac_mul_Bq(A_mul_Bc!, A, Q)
+function _Ac_mul_Bq(A_mul_Bop!::FT, A::StridedMatrix, Q::LQPackedQ) where FT<:Function
+    TR = promote_type(eltype(A), eltype(Q))
+    if size(A, 1) == size(Q.factors, 2)
+        C = adjoint!(similar(A, TR, reverse(size(A))), A)
+    elseif size(A, 1) == size(Q.factors, 1)
+        C = zeros(TR, size(A, 2), size(Q.factors, 2))
+        adjoint!(view(C, :, 1:size(A, 1)), A)
+    else
+        _rightappdimmismatch("rows")
     end
+    return A_mul_Bop!(C, convert(AbstractMatrix{TR}, Q))
 end
+_rightappdimmismatch(rowsorcols) =
+    throw(DimensionMismatch(string("the number of $(rowsorcols) of the matrix on the left ",
+        "must match either (1) the number of columns of the (LQPackedQ) matrix on the right ",
+        "or (2) the number of rows of that (LQPackedQ) matrix's internal representation ",
+        "(the factorization's originating matrix's number of rows)")))
+
 
 function (\)(A::LQ{TA}, b::StridedVector{Tb}) where {TA,Tb}
     S = promote_type(TA,Tb)

test/linalg/lq.jl

@@ -181,3 +181,40 @@ end
         @test size(lqfact(randn(mA, nA))[:Q]) == (nQ, nQ)
     end
 end
+
+@testset "postmultiplication with / right-application of LQPackedQ (#23779)" begin
+    function getqs(F::Base.LinAlg.LQ)
+        implicitQ = F[:Q]
+        explicitQ = A_mul_B!(implicitQ, eye(eltype(implicitQ), size(implicitQ)...))
+        return implicitQ, explicitQ
+    end
+    # for any shape m-by-n of LQ-factored matrix, where Q is an LQPackedQ
+    # A_mul_B*(C, Q) (Ac_mul_B*(C, Q)) operations should work for
+    # *-by-n (n-by-*) C, which we test below via n-by-n C
+    for (mA, nA) in ((3, 3), (3, 4), (4, 3))
+        implicitQ, explicitQ = getqs(lqfact(randn(mA, nA)))
+        C = randn(nA, nA)
+        @test *(C, implicitQ) ≈ *(C, explicitQ)
+        @test A_mul_Bc(C, implicitQ) ≈ A_mul_Bc(C, explicitQ)
+        @test Ac_mul_B(C, implicitQ) ≈ Ac_mul_B(C, explicitQ)
+        @test Ac_mul_Bc(C, implicitQ) ≈ Ac_mul_Bc(C, explicitQ)
+    end
+    # where the LQ-factored matrix has at least as many rows m as columns n,
+    # Q's square and thin forms have the same shape (n-by-n). hence we expect
+    # _only_ *-by-n (n-by-*) C to work in A_mul_B*(C, Q) (Ac_mul_B*(C, Q)) ops.
+    # and hence the n-by-n C tests above suffice.
+    #
+    # where the LQ-factored matrix has more columns n than rows m,
+    # Q's square form is n-by-n whereas its thin form is m-by-n.
+    # hence we need also test *-by-m (m-by-*) C with
+    # A_mul_B*(C, Q) (Ac_mul_B*(C, Q)) ops, as below via m-by-m C.
+    mA, nA = 3, 4
+    implicitQ, explicitQ = getqs(lqfact(randn(mA, nA)))
+    C = randn(mA, mA)
+    zeroextCright = hcat(C, zeros(eltype(C), mA))
+    zeroextCdown = vcat(C, zeros(eltype(C), (1, mA)))
+    @test *(C, implicitQ) ≈ *(zeroextCright, explicitQ)
+    @test A_mul_Bc(C, implicitQ) ≈ A_mul_Bc(zeroextCright, explicitQ)
+    @test Ac_mul_B(C, implicitQ) ≈ Ac_mul_B(zeroextCdown, explicitQ)
+    @test Ac_mul_Bc(C, implicitQ) ≈ Ac_mul_Bc(zeroextCdown, explicitQ)
+end