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

Speed up SparseMatrixCSC(::AbstractQ) constructor and multiplication #196

Merged
merged 8 commits into from
Aug 2, 2022
Merged
Show file tree
Hide file tree
Changes from 5 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
16 changes: 9 additions & 7 deletions src/solvers/spqr.jl
Original file line number Diff line number Diff line change
Expand Up @@ -2,7 +2,7 @@

module SPQR

import Base: \
import Base: \, *
using Base: require_one_based_indexing
using LinearAlgebra
using ..LibSuiteSparse: SuiteSparseQR_C
Expand Down Expand Up @@ -132,7 +132,7 @@ struct QRSparseQ{Tv<:CHOLMOD.VTypes,Ti<:Integer} <: LinearAlgebra.AbstractQ{Tv}
τ::Vector{Tv}
n::Int # Number of columns in original matrix
end

Base.print_array(io::IO, Q::QRSparseQ) = Base.print_array(io, sparse(Q))
SobhanMP marked this conversation as resolved.
Show resolved Hide resolved
Base.size(Q::QRSparseQ) = (size(Q.factors, 1), size(Q.factors, 1))
Base.axes(Q::QRSparseQ) = map(Base.OneTo, size(Q))

Expand Down Expand Up @@ -167,11 +167,7 @@ julia> A = sparse([1,2,3,4], [1,1,2,2], [1.0,1.0,1.0,1.0])
julia> qr(A)
SparseArrays.SPQR.QRSparse{Float64, Int64}
Q factor:
4×4 SparseArrays.SPQR.QRSparseQ{Float64, Int64}:
-0.707107 0.0 0.0 -0.707107
0.0 -0.707107 -0.707107 0.0
0.0 -0.707107 0.707107 0.0
-0.707107 0.0 0.0 0.707107
4×4 SparseArrays.SPQR.QRSparseQ{Float64, Int64}
R factor:
2×2 SparseMatrixCSC{Float64, Int64} with 2 stored entries:
-1.41421 ⋅
Expand Down Expand Up @@ -284,6 +280,9 @@ function LinearAlgebra.rmul!(A::StridedMatrix, adjQ::Adjoint{<:Any,<:QRSparseQ})
return A
end

(*)(Q::QRSparseQ, B::SparseMatrixCSC) = sparse(Q) * B
(*)(A::SparseMatrixCSC, Q::QRSparseQ) = A * sparse(Q)

@inline function Base.getproperty(F::QRSparse, d::Symbol)
if d === :Q
return QRSparseQ(F.factors, F.τ, size(F, 2))
Expand Down Expand Up @@ -312,6 +311,9 @@ function Base.show(io::IO, mime::MIME{Symbol("text/plain")}, F::QRSparse)
println(io, "\nColumn permutation:")
show(io, mime, F.pcol)
end
function Base.show(io::IO, ::MIME{Symbol("text/plain")}, Q::QRSparseQ)
summary(io, Q)
end

# With a real lhs and complex rhs with the same precision, we can reinterpret
# the complex rhs as a real rhs with twice the number of columns
Expand Down
45 changes: 35 additions & 10 deletions src/sparsematrix.jl
Original file line number Diff line number Diff line change
Expand Up @@ -748,6 +748,31 @@ SparseMatrixCSC{Tv}(M::Transpose{<:Any,<:AbstractSparseMatrixCSC}) where {Tv} =
SparseMatrixCSC{Tv,Ti}(M::Adjoint{<:Any,<:AbstractSparseMatrixCSC}) where {Tv,Ti} = SparseMatrixCSC{Tv,Ti}(copy(M))
SparseMatrixCSC{Tv,Ti}(M::Transpose{<:Any,<:AbstractSparseMatrixCSC}) where {Tv,Ti} = SparseMatrixCSC{Tv,Ti}(copy(M))

# we can only view AbstractQs as columns
function SparseMatrixCSC{Tv,Ti}(Q::LinearAlgebra.AbstractQ) where {Tv,Ti}
colptr = zeros(Ti, size(Q, 2) + 1)
nzval = Tv[]
rowval = Ti[]
col = zeros(eltype(Q), size(Q, 1))

colptr[1] = 1
ind = 1
for j in axes(Q, 2)
fill!(col, false)
col[j] = one(Tv)
lmul!(Q, col)
for (i, v) in enumerate(col)
if iszero(v) == false # should be _isnotzero(v)
push!(nzval, v)
push!(rowval, i)
ind += 1
end
end
colptr[j + 1] = ind
end
return SparseMatrixCSC{Tv,Ti}(size(Q)..., colptr, rowval, nzval)
end

# converting from SparseMatrixCSC to other matrix types
function Matrix(S::AbstractSparseMatrixCSC{Tv}) where Tv
_checkbuffers(S)
Expand Down Expand Up @@ -1137,9 +1162,9 @@ end
"""
ftranspose!(X::AbstractSparseMatrixCSC{Tv,Ti}, A::AbstractSparseMatrixCSC{Tv,Ti}, f::Function) where {Tv,Ti}

Transpose `A` and store it in `X` while applying the function `f` to the non-zero elements.
Transpose `A` and store it in `X` while applying the function `f` to the non-zero elements.
Does not remove the zeros created by `f`. `size(X)` must be equal to `size(transpose(A))`.
No additonal memory is allocated other than resizing the rowval and nzval of `X`, if needed.
No additonal memory is allocated other than resizing the rowval and nzval of `X`, if needed.

See `halfperm!`
"""
Expand All @@ -1158,9 +1183,9 @@ end
"""
transpose!(X::AbstractSparseMatrixCSC{Tv,Ti}, A::AbstractSparseMatrixCSC{Tv,Ti}) where {Tv,Ti}

Transpose the matrix `A` and stores it in the matrix `X`.
Transpose the matrix `A` and stores it in the matrix `X`.
`size(X)` must be equal to `size(transpose(A))`.
No additonal memory is allocated other than resizing the rowval and nzval of `X`, if needed.
No additonal memory is allocated other than resizing the rowval and nzval of `X`, if needed.

See `halfperm!`
"""
Expand All @@ -1170,8 +1195,8 @@ transpose!(X::AbstractSparseMatrixCSC{Tv,Ti}, A::AbstractSparseMatrixCSC{Tv,Ti})
adjoint!(X::AbstractSparseMatrixCSC{Tv,Ti}, A::AbstractSparseMatrixCSC{Tv,Ti}) where {Tv,Ti}

Transpose the matrix `A` and stores the adjoint of the elements in the matrix `X`.
`size(X)` must be equal to `size(transpose(A))`.
No additonal memory is allocated other than resizing the rowval and nzval of `X`, if needed.
`size(X)` must be equal to `size(transpose(A))`.
No additonal memory is allocated other than resizing the rowval and nzval of `X`, if needed.

See `halfperm!`
"""
Expand Down Expand Up @@ -3863,13 +3888,13 @@ end

## Uniform matrix arithmetic

(+)(A::AbstractSparseMatrixCSC{Tv, Ti}, J::UniformScaling{T}) where {T<:Number, Tv, Ti} =
(+)(A::AbstractSparseMatrixCSC{Tv, Ti}, J::UniformScaling{T}) where {T<:Number, Tv, Ti} =
A + sparse(T, Ti, J, size(A)...)
(+)(J::UniformScaling{T}, A::AbstractSparseMatrixCSC{Tv, Ti}) where {T<:Number, Tv, Ti} =
(+)(J::UniformScaling{T}, A::AbstractSparseMatrixCSC{Tv, Ti}) where {T<:Number, Tv, Ti} =
sparse(T, Ti, J, size(A)...) + A
(-)(A::AbstractSparseMatrixCSC{Tv, Ti}, J::UniformScaling{T}) where {T<:Number, Tv, Ti} =
(-)(A::AbstractSparseMatrixCSC{Tv, Ti}, J::UniformScaling{T}) where {T<:Number, Tv, Ti} =
A - sparse(T, Ti, J, size(A)...)
(-)(J::UniformScaling{T}, A::AbstractSparseMatrixCSC{Tv, Ti}) where {T<:Number, Tv, Ti} =
(-)(J::UniformScaling{T}, A::AbstractSparseMatrixCSC{Tv, Ti}) where {T<:Number, Tv, Ti} =
sparse(T, Ti, J, size(A)...) - A


Expand Down
13 changes: 12 additions & 1 deletion test/spqr.jl
Original file line number Diff line number Diff line change
Expand Up @@ -90,12 +90,14 @@ 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)
@test sparse(qr(A).Q) == sparse(qr(Matrix(A)).Q) == Matrix(I, 2, 2)
@test (sparse(I, 2, 2) * F.Q)::SparseMatrixCSC == sparse(qr(A).Q) == sparse(I, 2, 2)
end

@testset "Issue 26368" begin
A = sparse([0.0 1 0 0; 0 0 0 0])
F = qr(A)
@test F.Q*F.R == A[F.prow,F.pcol]
@test (F.Q*F.R)::SparseMatrixCSC == A[F.prow,F.pcol]
end

@testset "select ordering overdetermined" begin
Expand Down Expand Up @@ -137,6 +139,15 @@ end
@test all(iszero, (rank(spzeros(10, i)) for i in 1:10))
end


@testset "sparse" begin
A = I + sprandn(100, 100, 0.01)
Q = qr(A).Q
sQ = sparse(Q)
@test sQ == sparse(Matrix(Q))
@test sQ ≈ sparse(qr(Matrix(A)).Q)
dkarrasch marked this conversation as resolved.
Show resolved Hide resolved
end

end

end # Base.USE_GPL_LIBS
Expand Down