Skip to content

Commit

Permalink
Fix Matrix(QRSparseQ) constructor (JuliaLang#26392)
Browse files Browse the repository at this point in the history
* Fix Matrix(QRSparseQ) constructor

Fixes JuliaLang#26367

* Make sure that size(Q, i) definition covers HessenbergQ
  • Loading branch information
andreasnoack committed Nov 28, 2018
1 parent 9c8560f commit 076ea45
Show file tree
Hide file tree
Showing 4 changed files with 25 additions and 18 deletions.
29 changes: 14 additions & 15 deletions stdlib/LinearAlgebra/src/qr.jl
Original file line number Diff line number Diff line change
Expand Up @@ -502,23 +502,22 @@ AbstractMatrix{T}(Q::QRPackedQ) where {T} = QRPackedQ{T}(Q)
QRCompactWYQ{S}(Q::QRCompactWYQ) where {S} = QRCompactWYQ(convert(AbstractMatrix{S}, Q.factors), convert(AbstractMatrix{S}, Q.T))
AbstractMatrix{S}(Q::QRCompactWYQ{S}) where {S} = Q
AbstractMatrix{S}(Q::QRCompactWYQ) where {S} = QRCompactWYQ{S}(Q)
Matrix{T}(A::AbstractQ) where {T} = lmul!(A, Matrix{T}(I, size(A.factors, 1), min(size(A.factors)...)))
Matrix(A::AbstractQ{T}) where {T} = Matrix{T}(A)
Array{T}(A::AbstractQ) where {T} = Matrix{T}(A)
Array(A::AbstractQ) = Matrix(A)

size(A::Union{QR,QRCompactWY,QRPivoted}, dim::Integer) = size(getfield(A, :factors), dim)
size(A::Union{QR,QRCompactWY,QRPivoted}) = size(getfield(A, :factors))
size(A::AbstractQ, dim::Integer) = 0 < dim ? (dim <= 2 ? size(getfield(A, :factors), 1) : 1) : throw(BoundsError())
size(A::AbstractQ) = size(A, 1), size(A, 2)


function getindex(A::AbstractQ, i::Integer, j::Integer)
x = zeros(eltype(A), size(A, 1))
Matrix{T}(Q::AbstractQ) where {T} = lmul!(Q, Matrix{T}(I, size(Q, 1), min(size(Q.factors)...)))
Matrix(Q::AbstractQ{T}) where {T} = Matrix{T}(Q)
Array{T}(Q::AbstractQ) where {T} = Matrix{T}(Q)
Array(Q::AbstractQ) = Matrix(Q)

size(F::Union{QR,QRCompactWY,QRPivoted}, dim::Integer) = size(getfield(F, :factors), dim)
size(F::Union{QR,QRCompactWY,QRPivoted}) = size(getfield(F, :factors))
size(Q::AbstractQ, dim::Integer) = size(getfield(Q, :factors), dim == 2 ? 1 : dim)
size(Q::AbstractQ) = size(Q, 1), size(Q, 2)

function getindex(Q::AbstractQ, i::Integer, j::Integer)
x = zeros(eltype(Q), size(Q, 1))
x[i] = 1
y = zeros(eltype(A), size(A, 2))
y = zeros(eltype(Q), size(Q, 2))
y[j] = 1
return dot(x, lmul!(A, y))
return dot(x, lmul!(Q, y))
end

## Multiplication by Q
Expand Down
4 changes: 2 additions & 2 deletions stdlib/LinearAlgebra/test/qr.jl
Original file line number Diff line number Diff line change
Expand Up @@ -145,7 +145,7 @@ rectangularQ(Q::LinearAlgebra.AbstractQ) = convert(Array, Q)
@test_throws DimensionMismatch rmul!(Matrix{eltya}(I, n+1, n+1),q)
@test rmul!(squareQ(q), adjoint(q)) Matrix(I, n, n)
@test_throws DimensionMismatch rmul!(Matrix{eltya}(I, n+1, n+1), adjoint(q))
@test_throws BoundsError size(q,-1)
@test_throws ErrorException size(q,-1)
@test_throws DimensionMismatch LinearAlgebra.lmul!(q,zeros(eltya,n1+1))
@test_throws DimensionMismatch LinearAlgebra.lmul!(adjoint(q), zeros(eltya,n1+1))

Expand All @@ -155,7 +155,7 @@ rectangularQ(Q::LinearAlgebra.AbstractQ) = convert(Array, Q)
@test_throws DimensionMismatch rmul!(Matrix{eltya}(I, n+1, n+1),q)
@test rmul!(squareQ(q), adjoint(q)) Matrix(I, n, n)
@test_throws DimensionMismatch rmul!(Matrix{eltya}(I, n+1, n+1),adjoint(q))
@test_throws BoundsError size(q,-1)
@test_throws ErrorException size(q,-1)
@test_throws DimensionMismatch q * Matrix{Int8}(I, n+4, n+4)
end
end
Expand Down
5 changes: 4 additions & 1 deletion stdlib/SuiteSparse/src/spqr.jl
Original file line number Diff line number Diff line change
Expand Up @@ -130,11 +130,14 @@ Base.axes(F::QRSparse) = map(Base.OneTo, size(F))
struct QRSparseQ{Tv<:CHOLMOD.VTypes,Ti<:Integer} <: LinearAlgebra.AbstractQ{Tv}
factors::SparseMatrixCSC{Tv,Ti}
τ::Vector{Tv}
n::Int # Number of columns in original matrix
end

Base.size(Q::QRSparseQ) = (size(Q.factors, 1), size(Q.factors, 1))
Base.axes(Q::QRSparseQ) = map(Base.OneTo, size(Q))

Matrix{T}(Q::QRSparseQ) where {T} = lmul!(Q, Matrix{T}(I, size(Q, 1), min(size(Q, 1), Q.n)))

# From SPQR manual p. 6
_default_tol(A::SparseMatrixCSC) =
20*sum(size(A))*eps(real(eltype(A)))*maximum(norm(view(A, :, i)) for i in 1:size(A, 2))
Expand Down Expand Up @@ -305,7 +308,7 @@ julia> F.pcol
"""
@inline function Base.getproperty(F::QRSparse, d::Symbol)
if d == :Q
return QRSparseQ(F.factors, F.τ)
return QRSparseQ(F.factors, F.τ, size(F, 2))
elseif d == :prow
return invperm(F.rpivinv)
elseif d == :pcol
Expand Down
5 changes: 5 additions & 0 deletions stdlib/SuiteSparse/test/spqr.jl
Original file line number Diff line number Diff line change
Expand Up @@ -81,6 +81,11 @@ end
@test A*xs A*xd
end

@testset "Issue 26367" begin
A = sparse([0.0 1 0 0; 0 0 0 0])
@test Matrix(qr(A).Q) == Matrix(qr(Matrix(A)).Q) == Matrix(I, 2, 2)
end

@testset "Issue 26368" begin
A = sparse([0.0 1 0 0; 0 0 0 0])
F = qr(A)
Expand Down

0 comments on commit 076ea45

Please sign in to comment.