Skip to content

Commit

Permalink
Eliminate uses of full from base/linalg/lq.jl and fix/test lq(A, thin…
Browse files Browse the repository at this point in the history
…=false) shape. (#23729)
  • Loading branch information
Sacha0 authored and andreasnoack committed Sep 19, 2017
1 parent aa79148 commit 417331c
Show file tree
Hide file tree
Showing 2 changed files with 37 additions and 20 deletions.
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

0 comments on commit 417331c

Please sign in to comment.