From 1a3a63349a9b94299c7d84b8082ef32fa745283c Mon Sep 17 00:00:00 2001 From: Sacha Verweij Date: Thu, 21 Sep 2017 00:12:07 -0700 Subject: [PATCH] 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 | 6 ++-- base/linalg/lq.jl | 79 +++++++++++++++++++++++++++---------------- test/linalg/lq.jl | 37 ++++++++++++++++++++ 3 files changed, 90 insertions(+), 32 deletions(-) diff --git a/base/linalg/lapack.jl b/base/linalg/lapack.jl index 4db08f1f57d54..c8a5d54f8fc4f 100644 --- a/base/linalg/lapack.jl +++ b/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")) diff --git a/base/linalg/lq.jl b/base/linalg/lq.jl index 14e90ef113587..d4b6f7fa3e986 100644 --- a/base/linalg/lq.jl +++ b/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) diff --git a/test/linalg/lq.jl b/test/linalg/lq.jl index aaf327efa5b0b..e554319e71e00 100644 --- a/test/linalg/lq.jl +++ b/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