Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

eliminate uses of full from base/linalg/lq.jl #23729

Merged
merged 1 commit into from
Sep 19, 2017
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
32 changes: 13 additions & 19 deletions base/linalg/lq.jl
Original file line number Diff line number Diff line change
Expand Up @@ -38,11 +38,12 @@ lqfact(x::Number) = lqfact(fill(x,1,1))
Perform an LQ factorization of `A` such that `A = L*Q`. The
default is to compute a thin factorization. The LQ factorization
is the QR factorization of `A.'`. `L` is not extended with
zeros if the full `Q` is requested.
zeros if the explicit, square form of `Q` is requested via `thin = false`.
"""
function lq(A::Union{Number, AbstractMatrix}; thin::Bool=true)
function lq(A::Union{Number,AbstractMatrix}; thin::Bool = true)
F = lqfact(A)
F[:L], full(F[:Q], thin=thin)
L, Q = F[:L], F[:Q]
return L, thin ? Array(Q) : A_mul_B!(Q, eye(eltype(Q), size(Q.factors, 2)))
end

copy(A::LQ) = LQ(copy(A.factors), copy(A.τ))
Expand Down Expand Up @@ -90,19 +91,9 @@ convert(::Type{LQPackedQ{T}}, Q::LQPackedQ) where {T} = LQPackedQ(convert(Abstra
convert(::Type{AbstractMatrix{T}}, Q::LQPackedQ) where {T} = convert(LQPackedQ{T}, Q)
convert(::Type{Matrix}, A::LQPackedQ) = LAPACK.orglq!(copy(A.factors),A.τ)
convert(::Type{Array}, A::LQPackedQ) = convert(Matrix, A)
function full(A::LQPackedQ{T}; thin::Bool = true) where T
#= We construct the full eye here, even though it seems inefficient, because
every element in the output matrix is a function of all the elements of
the input matrix. The eye is modified by the elementary reflectors held
in A, so this is not just an indexing operation. Note that in general
explicitly constructing Q, rather than using the ldiv or mult methods,
may be a wasteful allocation. =#
if thin
convert(Array, A)
else
A_mul_B!(A, eye(T, size(A.factors,2), size(A.factors,1)))
end
end

full(Q::LQPackedQ; thin::Bool = true) =
thin ? Array(Q) : A_mul_B!(Q, eye(eltype(Q), size(Q.factors, 2)))

size(A::LQ, dim::Integer) = size(A.factors, dim)
size(A::LQ) = size(A.factors)
Expand All @@ -119,9 +110,12 @@ end
size(A::LQPackedQ) = size(A.factors)

## Multiplication by LQ
A_mul_B!(A::LQ{T}, B::StridedVecOrMat{T}) where {T<:BlasFloat} = A[:L]*LAPACK.ormlq!('L','N',A.factors,A.τ,B)
A_mul_B!(A::LQ{T}, B::QR{T}) where {T<:BlasFloat} = A[:L]*LAPACK.ormlq!('L','N',A.factors,A.τ,full(B))
A_mul_B!(A::QR{T}, B::LQ{T}) where {T<:BlasFloat} = A_mul_B!(zeros(full(A)), full(A), full(B))
A_mul_B!(A::LQ{T}, B::StridedVecOrMat{T}) where {T<:BlasFloat} =
A[:L] * LAPACK.ormlq!('L', 'N', A.factors, A.τ, B)
A_mul_B!(A::LQ{T}, B::QR{T}) where {T<:BlasFloat} =
A[:L] * LAPACK.ormlq!('L', 'N', A.factors, A.τ, Matrix(B))
A_mul_B!(A::QR{T}, B::LQ{T}) where {T<:BlasFloat} =
A_mul_B!(zeros(eltype(A), size(A)), Matrix(A), Matrix(B))
function *(A::LQ{TA}, B::StridedVecOrMat{TB}) where {TA,TB}
TAB = promote_type(TA, TB)
A_mul_B!(convert(Factorization{TAB},A), copy_oftype(B, TAB))
Expand Down
25 changes: 24 additions & 1 deletion test/linalg/lq.jl
Original file line number Diff line number Diff line change
Expand Up @@ -95,11 +95,34 @@ bimg = randn(n,2)/2
lqa = lqfact(a[:,1:n1])
l,q = lqa[:L], lqa[:Q]
@test full(q)*full(q)' ≈ eye(eltya,n1)
@test (full(q,thin=false)'*full(q,thin=false))[1:n1,:] ≈ eye(eltya,n1,n)
@test full(q,thin=false)'*full(q,thin=false) ≈ eye(eltya, n1)
@test_throws DimensionMismatch A_mul_B!(eye(eltya,n+1),q)
@test Ac_mul_B!(q,full(q)) ≈ eye(eltya,n1)
@test_throws DimensionMismatch A_mul_Bc!(eye(eltya,n+1),q)
@test_throws BoundsError size(q,-1)
end
end
end

@testset "correct form of Q from lq(...) (#23729)" begin
# matrices with more rows than columns
m, n = 4, 2
A = randn(m, n)
for thin in (true, false)
L, Q = lq(A, thin = thin)
@test size(L) == (m, n)
@test size(Q) == (n, n)
@test isapprox(A, L*Q)
end
# matrices with more columns than rows
m, n = 2, 4
A = randn(m, n)
Lthin, Qthin = lq(A, thin = true)
@test size(Lthin) == (m, m)
@test size(Qthin) == (m, n)
@test isapprox(A, Lthin * Qthin)
Lsquare, Qsquare = lq(A, thin = false)
@test size(Lsquare) == (m, m)
@test size(Qsquare) == (n, n)
@test isapprox(A, [Lsquare zeros(m, n - m)] * Qsquare)
end