fix LAPACK.ormlq! and postmultiplication with 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

@@ -92,17 +92,22 @@ full(Q::LQPackedQ; thin::Bool = true) =
 
 size(A::LQ, dim::Integer) = size(A.factors, dim)
 size(A::LQ) = size(A.factors)
-function size(A::LQPackedQ, dim::Integer)
-    if 0 < dim && dim <= 2
-        return size(A.factors, dim)
-    elseif 0 < dim && dim > 2
-        return 1
-    else
+
+# size(Q::LQPackedQ) yields the shape of Q's square form
+function size(Q::LQPackedQ)
+    n = size(Q.factors, 2)
+    return n, n
+end
+function size(Q::LQPackedQ, dim::Integer)
+    if dim < 1
         throw(BoundsError())
+    elseif dim <= 2 # && 1 <= dim
+        return size(Q.factors, 2)
+    else # 2 < dim
+        return 1
     end
 end
 
-size(A::LQPackedQ) = size(A.factors)
 
 ## Multiplication by LQ
 A_mul_B!(A::LQ{T}, B::StridedVecOrMat{T}) where {T<:BlasFloat} =
@@ -159,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

@@ -105,7 +105,12 @@ bimg  = randn(n,2)/2
 end
 
 @testset "correct form of Q from lq(...) (#23729)" begin
-    # matrices with more rows than columns
+    # where the original matrix (say A) is square or has more rows than columns,
+    # then A's factorization's triangular factor (say L) should have the same shape
+    # as A independent of form (thin, square), and A's factorization's orthogonal
+    # factor (say Q) should be a square matrix of order of A's number of columns
+    # independent of form (thin, square), and L and Q should have
+    # multiplication-compatible shapes.
     m, n = 4, 2
     A = randn(m, n)
     for thin in (true, false)
@@ -114,13 +119,22 @@ end
         @test size(Q) == (n, n)
         @test isapprox(A, L*Q)
     end
-    # matrices with more columns than rows
+    # where the original matrix has strictly fewer rows than columns ...
     m, n = 2, 4
     A = randn(m, n)
+    # ... then, for a thin factorization of A, L should be a square matrix
+    # of order of A's number of rows, Q should have the same shape as A,
+    # and L and Q should have multiplication-compatible shapes
     Lthin, Qthin = lq(A, thin = true)
     @test size(Lthin) == (m, m)
     @test size(Qthin) == (m, n)
     @test isapprox(A, Lthin * Qthin)
+    # ... and, for a non-thin factorization of A, L should have the same shape as A,
+    # Q should be a square matrix of order of A's number of columns, and L and Q
+    # should have multiplication-compatible shape. but instead the L returned has
+    # no zero-padding on the right / is L for the thin factorization, so for
+    # L and Q to have multiplication-compatible shapes, L must be zero-padded
+    # to have the shape of A.
     Lsquare, Qsquare = lq(A, thin = false)
     @test size(Lsquare) == (m, m)
     @test size(Qsquare) == (n, n)
@@ -157,3 +171,50 @@ end
     @test implicitQ[1, n] == explicitQ[1, n]
     @test implicitQ[n, n] == explicitQ[n, n]
 end
+
+@testset "size on LQPackedQ (#23780)" begin
+    # size(Q::LQPackedQ) yields the shape of Q's square form
+    for ((mA, nA), nQ) in (
+        ((3, 3), 3), # A 3-by-3 => square Q 3-by-3
+        ((3, 4), 4), # A 3-by-4 => square Q 4-by-4
+        ((4, 3), 3) )# A 4-by-3 => square Q 3-by-3
+        @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