Skip to content

Commit

Permalink
Make LQPackedQ's postmultiplication behavior context-dependent only w…
Browse files Browse the repository at this point in the history
…here necessary.
  • Loading branch information
Sacha0 committed Sep 22, 2017
1 parent ff9fb48 commit 557fe05
Show file tree
Hide file tree
Showing 2 changed files with 39 additions and 19 deletions.
54 changes: 38 additions & 16 deletions base/linalg/lq.jl
Original file line number Diff line number Diff line change
Expand Up @@ -176,17 +176,41 @@ A_mul_Bc!(A::StridedMatrix{T}, B::LQPackedQ{T}) where {T<:BlasReal} =
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
# out-of-place right application of LQPackedQs
#
# LQPackedQ's out-of-place multiplication behavior is context dependent. specifically,
# if the inner dimension in the multiplication is the LQPackedQ's second dimension,
# the LQPackedQ behaves like its full/square form. if the inner dimension in the
# multiplication is the LQPackedQ's first dimension, the LQPackedQ behaves like either
# its full/square form or its thin form depending on the shape of the other object
# involved in the multiplication. we treat these cases separately.
#
# (1) the inner dimension in the multiplication is the LQPackedQ's second dimension.
# in this case, the LQPackedQ behaves like its square/full form.
#
function A_mul_Bc(A::StridedVecOrMat, Q::LQPackedQ)
TR = promote_type(eltype(A), eltype(Q))
return A_mul_Bc!(copy_oftype(A, TR), convert(AbstractMatrix{TR}, Q))
end
function Ac_mul_Bc(A::StridedMatrix, Q::LQPackedQ)
TR = promote_type(eltype(A), eltype(Q))
C = adjoint!(similar(A, TR, reverse(size(A))), A)
return A_mul_Bc!(C, convert(AbstractMatrix{TR}, Q))
end
#
# (2) the inner dimension in the multiplication is the LQPackedQ's first dimension.
# in this case, the LQPackedQ behaves like either its square/full form or its
# thin form depending on the shape of the other object in the multiplication.
#
# 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),
# and if so effectively apply Q's square form to A without additional shenanigans; and
# (2) if the preceding dimensions do not match, 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 implicitly apply Q's thin form to A by zero extending
# A as necessary for check (1) to pass (if possible) and then applying Q's square form
#
function *(A::StridedVecOrMat, Q::LQPackedQ)
TR = promote_type(eltype(A), eltype(Q))
if size(A, 2) == size(Q.factors, 2)
C = copy_oftype(A, TR)
Expand All @@ -196,11 +220,9 @@ function _A_mul_Bq(A_mul_Bop!::FT, A::StridedVecOrMat, Q::LQPackedQ) where FT<:F
else
_rightappdimmismatch("columns")
end
return A_mul_Bop!(C, convert(AbstractMatrix{TR}, Q))
return A_mul_B!(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
function Ac_mul_B(A::StridedMatrix, Q::LQPackedQ)
TR = promote_type(eltype(A), eltype(Q))
if size(A, 1) == size(Q.factors, 2)
C = adjoint!(similar(A, TR, reverse(size(A))), A)
Expand All @@ -210,7 +232,7 @@ function _Ac_mul_Bq(A_mul_Bop!::FT, A::StridedMatrix, Q::LQPackedQ) where FT<:Fu
else
_rightappdimmismatch("rows")
end
return A_mul_Bop!(C, convert(AbstractMatrix{TR}, Q))
return A_mul_B!(C, convert(AbstractMatrix{TR}, Q))
end
_rightappdimmismatch(rowsorcols) =
throw(DimensionMismatch(string("the number of $(rowsorcols) of the matrix on the left ",
Expand Down
4 changes: 1 addition & 3 deletions test/linalg/lq.jl
Original file line number Diff line number Diff line change
Expand Up @@ -207,14 +207,12 @@ end
# 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.
# A*_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

0 comments on commit 557fe05

Please sign in to comment.