Skip to content

Commit

Permalink
RFC: Extract all parts of QR factorization from SPQR (#22626)
Browse files Browse the repository at this point in the history 
* Wrap full SuiteSparseQR function

* Introduce new QRSparse type and write Q functions in Julia

* Remove unused code and adjust tests

* Update docstrings and add some comments

* Define some aliases, fix nonzeros for views, and some minor fixes.

* Fix ambiguities and finish docstring for getindex

* Fix segfault by avoiding that Sparse object is getting freed while calling
SPQR. This is accomplished by passing Sparse and define unsafe_convert instead
of just passing the pointer.

Also fix memory leak.

* New ntuple use

* Add show method for QRSparse

Minor adjustments following review recommendations
  • Loading branch information
@andreasnoack
andreasnoack authored Jul 10, 2017
1 parent f2f0a75 commit 3334724
Show file tree
Hide file tree
Showing 7 changed files with 504 additions and 182 deletions.
62 changes: 36 additions & 26 deletions base/linalg/qr.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -422,6 +422,8 @@ function getindex(A::QRPivoted{T}, d::Symbol) where T
end
end

abstract type AbstractQ{T} <: AbstractMatrix{T} end

# Type-stable interface to get Q
getq(A::QRCompactWY) = QRCompactWYQ(A.factors,A.T)
getq(A::Union{QR, QRPivoted}) = QRPackedQ(A.factors,A.τ)
Expand All @@ -432,7 +434,7 @@ getq(A::Union{QR, QRPivoted}) = QRPackedQ(A.factors,A.τ)
The orthogonal/unitary ``Q`` matrix of a QR factorization stored in [`QR`](@ref) or
[`QRPivoted`](@ref) format.
"""
struct QRPackedQ{T,S<:AbstractMatrix} <: AbstractMatrix{T}
struct QRPackedQ{T,S<:AbstractMatrix} <: AbstractQ{T}
factors::S
τ::Vector{T}
QRPackedQ{T,S}(factors::AbstractMatrix{T}, τ::Vector{T}) where {T,S<:AbstractMatrix} = new(factors, τ)
Expand All @@ -445,7 +447,7 @@ QRPackedQ(factors::AbstractMatrix{T}, τ::Vector{T}) where {T} = QRPackedQ{T,typ
The orthogonal/unitary ``Q`` matrix of a QR factorization stored in [`QRCompactWY`](@ref)
format.
"""
struct QRCompactWYQ{S, M<:AbstractMatrix} <: AbstractMatrix{S}
struct QRCompactWYQ{S, M<:AbstractMatrix} <: AbstractQ{S}
factors::M
T::Matrix{S}
QRCompactWYQ{S,M}(factors::AbstractMatrix{S}, T::Matrix{S}) where {S,M<:AbstractMatrix} = new(factors, T)
Expand All @@ -458,11 +460,11 @@ convert(::Type{AbstractMatrix{T}}, Q::QRPackedQ) where {T} = convert(QRPackedQ{T
convert(::Type{QRCompactWYQ{S}}, Q::QRCompactWYQ) where {S} = QRCompactWYQ(convert(AbstractMatrix{S}, Q.factors), convert(AbstractMatrix{S}, Q.T))
convert(::Type{AbstractMatrix{S}}, Q::QRCompactWYQ{S}) where {S} = Q
convert(::Type{AbstractMatrix{S}}, Q::QRCompactWYQ) where {S} = convert(QRCompactWYQ{S}, Q)
convert(::Type{Matrix}, A::Union{QRPackedQ{T},QRCompactWYQ{T}}) where {T} = A_mul_B!(A, eye(T, size(A.factors, 1), minimum(size(A.factors))))
convert(::Type{Array}, A::Union{QRPackedQ,QRCompactWYQ}) = convert(Matrix, A)
convert(::Type{Matrix}, A::AbstractQ{T}) where {T} = A_mul_B!(A, eye(T, size(A.factors, 1), minimum(size(A.factors))))
convert(::Type{Array}, A::AbstractQ) = convert(Matrix, A)

"""
full(A::Union{QRPackedQ,QRCompactWYQ}; thin::Bool=true) -> Matrix
full(A::AbstractQ; thin::Bool=true) -> Matrix
Converts an orthogonal or unitary matrix stored as a `QRCompactWYQ` object, i.e. in the
compact WY format [^Bischof1987], or in the `QRPackedQ` format, to a dense matrix.
Expand All @@ -472,7 +474,7 @@ rows of `R` in the QR factorization that are zero. The resulting matrix is the `
QR factorization (sometimes called the reduced QR factorization). If `false`, returns a `Q`
that spans all rows of `R` in its corresponding QR factorization.
"""
function full(A::Union{QRPackedQ{T},QRCompactWYQ{T}}; thin::Bool = true) where T
function full(A::AbstractQ{T}; thin::Bool = true) where T
if thin
convert(Array, A)
else
Expand All @@ -482,11 +484,11 @@ end

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


function getindex(A::Union{QRPackedQ,QRCompactWYQ}, i::Integer, j::Integer)
function getindex(A::AbstractQ, i::Integer, j::Integer)
x = zeros(eltype(A), size(A, 1))
x[i] = 1
y = zeros(eltype(A), size(A, 2))
Expand All @@ -496,8 +498,10 @@ end

## Multiplication by Q
### QB
A_mul_B!(A::QRCompactWYQ{T}, B::StridedVecOrMat{T}) where {T<:BlasFloat} = LAPACK.gemqrt!('L','N',A.factors,A.T,B)
A_mul_B!(A::QRPackedQ{T}, B::StridedVecOrMat{T}) where {T<:BlasFloat} = LAPACK.ormqr!('L','N',A.factors,A.τ,B)
A_mul_B!(A::QRCompactWYQ{T,S}, B::StridedVecOrMat{T}) where {T<:BlasFloat, S<:StridedMatrix} =
LAPACK.gemqrt!('L','N',A.factors,A.T,B)
A_mul_B!(A::QRPackedQ{T,S}, B::StridedVecOrMat{T}) where {T<:BlasFloat, S<:StridedMatrix} =
LAPACK.ormqr!('L','N',A.factors,A.τ,B)
function A_mul_B!(A::QRPackedQ, B::AbstractVecOrMat)
mA, nA = size(A.factors)
mB, nB = size(B,1), size(B,2)
Expand All @@ -523,7 +527,7 @@ function A_mul_B!(A::QRPackedQ, B::AbstractVecOrMat)
B
end

function (*)(A::Union{QRPackedQ,QRCompactWYQ}, b::StridedVector)
function (*)(A::AbstractQ, b::StridedVector)
TAb = promote_type(eltype(A), eltype(b))
Anew = convert(AbstractMatrix{TAb}, A)
if size(A.factors, 1) == length(b)
Expand All @@ -535,7 +539,7 @@ function (*)(A::Union{QRPackedQ,QRCompactWYQ}, b::StridedVector)
end
A_mul_B!(Anew, bnew)
end
function (*)(A::Union{QRPackedQ,QRCompactWYQ}, B::StridedMatrix)
function (*)(A::AbstractQ, B::StridedMatrix)
TAB = promote_type(eltype(A), eltype(B))
Anew = convert(AbstractMatrix{TAB}, A)
if size(A.factors, 1) == size(B, 1)
Expand All @@ -549,10 +553,14 @@ function (*)(A::Union{QRPackedQ,QRCompactWYQ}, B::StridedMatrix)
end

### QcB
Ac_mul_B!(A::QRCompactWYQ{T}, B::StridedVecOrMat{T}) where {T<:BlasReal} = LAPACK.gemqrt!('L','T',A.factors,A.T,B)
Ac_mul_B!(A::QRCompactWYQ{T}, B::StridedVecOrMat{T}) where {T<:BlasComplex} = LAPACK.gemqrt!('L','C',A.factors,A.T,B)
Ac_mul_B!(A::QRPackedQ{T}, B::StridedVecOrMat{T}) where {T<:BlasReal} = LAPACK.ormqr!('L','T',A.factors,A.τ,B)
Ac_mul_B!(A::QRPackedQ{T}, B::StridedVecOrMat{T}) where {T<:BlasComplex} = LAPACK.ormqr!('L','C',A.factors,A.τ,B)
Ac_mul_B!(A::QRCompactWYQ{T,S}, B::StridedVecOrMat{T}) where {T<:BlasReal,S<:StridedMatrix} =
LAPACK.gemqrt!('L','T',A.factors,A.T,B)
Ac_mul_B!(A::QRCompactWYQ{T,S}, B::StridedVecOrMat{T}) where {T<:BlasComplex,S<:StridedMatrix} =
LAPACK.gemqrt!('L','C',A.factors,A.T,B)
Ac_mul_B!(A::QRPackedQ{T,S}, B::StridedVecOrMat{T}) where {T<:BlasReal,S<:StridedMatrix} =
LAPACK.ormqr!('L','T',A.factors,A.τ,B)
Ac_mul_B!(A::QRPackedQ{T,S}, B::StridedVecOrMat{T}) where {T<:BlasComplex,S<:StridedMatrix} =
LAPACK.ormqr!('L','C',A.factors,A.τ,B)
function Ac_mul_B!(A::QRPackedQ, B::AbstractVecOrMat)
mA, nA = size(A.factors)
mB, nB = size(B,1), size(B,2)
Expand All @@ -577,7 +585,7 @@ function Ac_mul_B!(A::QRPackedQ, B::AbstractVecOrMat)
end
B
end
function Ac_mul_B(Q::Union{QRPackedQ,QRCompactWYQ}, B::StridedVecOrMat)
function Ac_mul_B(Q::AbstractQ, B::StridedVecOrMat)
TQB = promote_type(eltype(Q), eltype(B))
return Ac_mul_B!(convert(AbstractMatrix{TQB}, Q), copy_oftype(B, TQB))
end
Expand All @@ -586,7 +594,7 @@ end
for (f1, f2) in ((:A_mul_Bc, :A_mul_B!),
(:Ac_mul_Bc, :Ac_mul_B!))
@eval begin
function ($f1)(Q::Union{QRPackedQ,QRCompactWYQ}, B::StridedVecOrMat)
function ($f1)(Q::AbstractQ, B::StridedVecOrMat)
TQB = promote_type(eltype(Q), eltype(B))
Bc = similar(B, TQB, (size(B, 2), size(B, 1)))
ctranspose!(Bc, B)
Expand All @@ -596,8 +604,10 @@ for (f1, f2) in ((:A_mul_Bc, :A_mul_B!),
end

### AQ
A_mul_B!(A::StridedVecOrMat{T}, B::QRCompactWYQ{T}) where {T<:BlasFloat} = LAPACK.gemqrt!('R','N', B.factors, B.T, A)
A_mul_B!(A::StridedVecOrMat{T}, B::QRPackedQ{T}) where {T<:BlasFloat} = LAPACK.ormqr!('R', 'N', B.factors, B.τ, A)
A_mul_B!(A::StridedVecOrMat{T}, B::QRCompactWYQ{T,S}) where {T<:BlasFloat,S<:StridedMatrix} =
LAPACK.gemqrt!('R','N', B.factors, B.T, A)
A_mul_B!(A::StridedVecOrMat{T}, B::QRPackedQ{T,S}) where {T<:BlasFloat,S<:StridedMatrix} =
LAPACK.ormqr!('R', 'N', B.factors, B.τ, A)
function A_mul_B!(A::StridedMatrix,Q::QRPackedQ)
mQ, nQ = size(Q.factors)
mA, nA = size(A,1), size(A,2)
Expand All @@ -623,7 +633,7 @@ function A_mul_B!(A::StridedMatrix,Q::QRPackedQ)
A
end

function (*)(A::StridedMatrix, Q::Union{QRPackedQ,QRCompactWYQ})
function (*)(A::StridedMatrix, Q::AbstractQ)
TAQ = promote_type(eltype(A), eltype(Q))
return A_mul_B!(copy_oftype(A, TAQ), convert(AbstractMatrix{TAQ}, Q))
end
Expand All @@ -633,7 +643,7 @@ A_mul_Bc!(A::StridedVecOrMat{T}, B::QRCompactWYQ{T}) where {T<:BlasReal} = LAPAC
A_mul_Bc!(A::StridedVecOrMat{T}, B::QRCompactWYQ{T}) where {T<:BlasComplex} = LAPACK.gemqrt!('R','C',B.factors,B.T,A)
A_mul_Bc!(A::StridedVecOrMat{T}, B::QRPackedQ{T}) where {T<:BlasReal} = LAPACK.ormqr!('R','T',B.factors,B.τ,A)
A_mul_Bc!(A::StridedVecOrMat{T}, B::QRPackedQ{T}) where {T<:BlasComplex} = LAPACK.ormqr!('R','C',B.factors,B.τ,A)
function A_mul_Bc!(A::AbstractMatrix,Q::QRPackedQ)
function A_mul_Bc!(A::StridedMatrix,Q::QRPackedQ)
mQ, nQ = size(Q.factors)
mA, nA = size(A,1), size(A,2)
if nA != mQ
Expand All @@ -657,7 +667,7 @@ function A_mul_Bc!(A::AbstractMatrix,Q::QRPackedQ)
end
A
end
function A_mul_Bc(A::AbstractMatrix, B::Union{QRCompactWYQ,QRPackedQ})
function A_mul_Bc(A::StridedMatrix, B::AbstractQ)
TAB = promote_type(eltype(A),eltype(B))
BB = convert(AbstractMatrix{TAB}, B)
if size(A,2) == size(B.factors, 1)
Expand All @@ -670,14 +680,14 @@ function A_mul_Bc(A::AbstractMatrix, B::Union{QRCompactWYQ,QRPackedQ})
throw(DimensionMismatch("matrix A has dimensions $(size(A)) but matrix B has dimensions $(size(B))"))
end
end
@inline A_mul_Bc(rowvec::RowVector, B::Union{LinAlg.QRCompactWYQ,LinAlg.QRPackedQ}) = ctranspose(B*ctranspose(rowvec))
@inline A_mul_Bc(rowvec::RowVector, B::AbstractQ) = ctranspose(B*ctranspose(rowvec))


### AcQ/AcQc
for (f1, f2) in ((:Ac_mul_B, :A_mul_B!),
(:Ac_mul_Bc, :A_mul_Bc!))
@eval begin
function ($f1)(A::StridedVecOrMat, Q::Union{QRPackedQ,QRCompactWYQ})
function ($f1)(A::StridedVecOrMat, Q::AbstractQ)
TAQ = promote_type(eltype(A), eltype(Q))
Ac = similar(A, TAQ, (size(A, 2), size(A, 1)))
ctranspose!(Ac, A)
Expand Down
2 changes: 2 additions & 0 deletions base/sparse/cholmod.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -248,6 +248,8 @@ mutable struct Sparse{Tv<:VTypes} <: AbstractSparseMatrix{Tv,SuiteSparse_long}
end
Sparse(p::Ptr{C_Sparse{Tv}}) where {Tv<:VTypes} = Sparse{Tv}(p)

Base.unsafe_convert(::Type{Ptr{Tv}}, A::Sparse{Tv}) where {Tv} = A.p

# Factor

if build_version >= v"2.1.0" # CHOLMOD version 2.1.0 or later
Expand Down
24 changes: 12 additions & 12 deletions base/sparse/linalg.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -194,17 +194,17 @@ function spmatmul(A::SparseMatrixCSC{Tv,Ti}, B::SparseMatrixCSC{Tv,Ti};
end

## solvers
function fwdTriSolve!(A::SparseMatrixCSC, B::AbstractVecOrMat)
function fwdTriSolve!(A::SparseMatrixCSCUnion, B::AbstractVecOrMat)
# forward substitution for CSC matrices
nrowB, ncolB = size(B, 1), size(B, 2)
ncol = LinAlg.checksquare(A)
if nrowB != ncol
throw(DimensionMismatch("A is $(ncol) columns and B has $(nrowB) rows"))
end

aa = A.nzval
ja = A.rowval
ia = A.colptr
aa = getnzval(A)
ja = getrowval(A)
ia = getcolptr(A)

joff = 0
for k = 1:ncolB
Expand Down Expand Up @@ -239,17 +239,17 @@ function fwdTriSolve!(A::SparseMatrixCSC, B::AbstractVecOrMat)
B
end

function bwdTriSolve!(A::SparseMatrixCSC, B::AbstractVecOrMat)
function bwdTriSolve!(A::SparseMatrixCSCUnion, B::AbstractVecOrMat)
# backward substitution for CSC matrices
nrowB, ncolB = size(B, 1), size(B, 2)
ncol = LinAlg.checksquare(A)
if nrowB != ncol
throw(DimensionMismatch("A is $(ncol) columns and B has $(nrowB) rows"))
end

aa = A.nzval
ja = A.rowval
ia = A.colptr
aa = getnzval(A)
ja = getrowval(A)
ia = getcolptr(A)

joff = 0
for k = 1:ncolB
Expand Down Expand Up @@ -284,11 +284,11 @@ function bwdTriSolve!(A::SparseMatrixCSC, B::AbstractVecOrMat)
B
end

A_ldiv_B!(L::LowerTriangular{T,<:SparseMatrixCSC{T}}, B::StridedVecOrMat) where {T} = fwdTriSolve!(L.data, B)
A_ldiv_B!(U::UpperTriangular{T,<:SparseMatrixCSC{T}}, B::StridedVecOrMat) where {T} = bwdTriSolve!(U.data, B)
A_ldiv_B!(L::LowerTriangular{T,<:SparseMatrixCSCUnion{T}}, B::StridedVecOrMat) where {T} = fwdTriSolve!(L.data, B)
A_ldiv_B!(U::UpperTriangular{T,<:SparseMatrixCSCUnion{T}}, B::StridedVecOrMat) where {T} = bwdTriSolve!(U.data, B)

(\)(L::LowerTriangular{T,<:SparseMatrixCSC{T}}, B::SparseMatrixCSC) where {T} = A_ldiv_B!(L, Array(B))
(\)(U::UpperTriangular{T,<:SparseMatrixCSC{T}}, B::SparseMatrixCSC) where {T} = A_ldiv_B!(U, Array(B))
(\)(L::LowerTriangular{T,<:SparseMatrixCSCUnion{T}}, B::SparseMatrixCSC) where {T} = A_ldiv_B!(L, Array(B))
(\)(U::UpperTriangular{T,<:SparseMatrixCSCUnion{T}}, B::SparseMatrixCSC) where {T} = A_ldiv_B!(U, Array(B))

function A_rdiv_B!(A::SparseMatrixCSC{T}, D::Diagonal{T}) where T
dd = D.diag
Expand Down
18 changes: 18 additions & 0 deletions base/sparse/sparsematrix.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -35,6 +35,24 @@ end

size(S::SparseMatrixCSC) = (S.m, S.n)

# Define an alias for views of a SparseMatrixCSC which include all rows and a unit range of the columns.
# Also define a union of SparseMatrixCSC and this view since many methods can be defined efficiently for
# this union by extracting the fields via the get function: getcolptr, getrowval, and getnzval. The key
# insight is that getcolptr on a SparseMatrixCSCView returns an offset view of the colptr of the
# underlying SparseMatrixCSC
const SparseMatrixCSCView{Tv,Ti} =
SubArray{Tv,2,SparseMatrixCSC{Tv,Ti},
Tuple{Base.Slice{Base.OneTo{Int}},I}} where {I<:AbstractUnitRange}
const SparseMatrixCSCUnion{Tv,Ti} = Union{SparseMatrixCSC{Tv,Ti}, SparseMatrixCSCView{Tv,Ti}}

getcolptr(S::SparseMatrixCSC) = S.colptr
getcolptr(S::SparseMatrixCSCView) = view(S.parent.colptr, first(indices(S, 2)):(last(indices(S, 2)) + 1))
getrowval(S::SparseMatrixCSC) = S.rowval
getrowval(S::SparseMatrixCSCView) = S.parent.rowval
getnzval( S::SparseMatrixCSC) = S.nzval
getnzval( S::SparseMatrixCSCView) = S.parent.nzval


"""
nnz(A)
Expand Down
45 changes: 38 additions & 7 deletions base/sparse/sparsevector.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -30,15 +30,34 @@ end
SparseVector(n::Integer, nzind::Vector{Ti}, nzval::Vector{Tv}) where {Tv,Ti} =
SparseVector{Tv,Ti}(n, nzind, nzval)

# Define an alias for a view of a whole column of a SparseMatrixCSC. Many methods can be written for the
# union of such a view and a SparseVector so we define an alias for such a union as well
const SparseColumnView{T} = SubArray{T,1,<:SparseMatrixCSC,Tuple{Base.Slice{Base.OneTo{Int}},Int},false}
const SparseVectorUnion{T} = Union{SparseVector{T}, SparseColumnView{T}}

### Basic properties

length(x::SparseVector) = x.n
size(x::SparseVector) = (x.n,)
nnz(x::SparseVector) = length(x.nzval)
countnz(x::SparseVector) = countnz(x.nzval)
count(x::SparseVector) = count(x.nzval)

nonzeros(x::SparseVector) = x.nzval
function nonzeros(x::SparseColumnView)
rowidx, colidx = parentindexes(x)
A = parent(x)
@inbounds y = view(A.nzval, nzrange(A, colidx))
return y
end

nonzeroinds(x::SparseVector) = x.nzind
function nonzeroinds(x::SparseColumnView)
rowidx, colidx = parentindexes(x)
A = parent(x)
@inbounds y = view(A.rowval, nzrange(A, colidx))
return y
end

similar(x::SparseVector, Tv::Type=eltype(x)) = SparseVector(x.n, copy(x.nzind), Vector{Tv}(length(x.nzval)))
function similar(x::SparseVector, ::Type{Tv}, ::Type{Ti}) where {Tv,Ti}
Expand Down Expand Up @@ -1396,7 +1415,7 @@ for f in [:sum, :maximum, :minimum], op in [:abs, :abs2]
end
end

vecnorm(x::AbstractSparseVector, p::Real=2) = vecnorm(nonzeros(x), p)
vecnorm(x::SparseVectorUnion, p::Real=2) = vecnorm(nonzeros(x), p)

### linalg.jl

Expand All @@ -1409,7 +1428,7 @@ vecnorm(x::AbstractSparseVector, p::Real=2) = vecnorm(nonzeros(x), p)

# axpy

function LinAlg.axpy!(a::Number, x::AbstractSparseVector, y::StridedVector)
function LinAlg.axpy!(a::Number, x::SparseVectorUnion, y::StridedVector)
length(x) == length(y) || throw(DimensionMismatch())
nzind = nonzeroinds(x)
nzval = nonzeros(x)
Expand Down Expand Up @@ -1454,8 +1473,7 @@ broadcast(::typeof(*), a::Number, x::AbstractSparseVector) = a * x
broadcast(::typeof(/), x::AbstractSparseVector, a::Number) = x / a

# dot

function dot(x::StridedVector{Tx}, y::AbstractSparseVector{Ty}) where {Tx<:Number,Ty<:Number}
function dot(x::StridedVector{Tx}, y::SparseVectorUnion{Ty}) where {Tx<:Number,Ty<:Number}
n = length(x)
length(y) == n || throw(DimensionMismatch())
nzind = nonzeroinds(y)
Expand All @@ -1467,13 +1485,13 @@ function dot(x::StridedVector{Tx}, y::AbstractSparseVector{Ty}) where {Tx<:Numbe
return s
end

function dot(x::AbstractSparseVector{Tx}, y::AbstractVector{Ty}) where {Tx<:Number,Ty<:Number}
function dot(x::SparseVectorUnion{Tx}, y::AbstractVector{Ty}) where {Tx<:Number,Ty<:Number}
n = length(y)
length(x) == n || throw(DimensionMismatch())
nzind = nonzeroinds(x)
nzval = nonzeros(x)
s = zero(Tx) * zero(Ty)
for i = 1:length(nzind)
@inbounds for i = 1:length(nzind)
s += conj(nzval[i]) * y[nzind[i]]
end
return s
Expand All @@ -1500,7 +1518,7 @@ function _spdot(f::Function,
s
end

function dot(x::AbstractSparseVector{<:Number}, y::AbstractSparseVector{<:Number})
function dot(x::SparseVectorUnion{<:Number}, y::SparseVectorUnion{<:Number})
x === y && return sum(abs2, x)
n = length(x)
length(y) == n || throw(DimensionMismatch())
Expand All @@ -1518,6 +1536,19 @@ end

### BLAS-2 / dense A * sparse x -> dense y

# lowrankupdate (BLAS.ger! like)
function LinAlg.lowrankupdate!(A::StridedMatrix, x::StridedVector, y::SparseVectorUnion, α::Number = 1)
nzi = nonzeroinds(y)
nzv = nonzeros(y)
@inbounds for (j,v) in zip(nzi,nzv)
αv = α*v'
for i in indices(x, 1)
A[i,j] += x[i]*αv
end
end
return A
end

# A_mul_B

function (*)(A::StridedMatrix{Ta}, x::AbstractSparseVector{Tx}) where {Ta,Tx}
Expand Down
Loading

2 comments on commit 3334724

@nanosoldier
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Executing the daily benchmark build, I will reply here when finished:

@nanosoldier runbenchmarks(ALL, isdaily = true)

@nanosoldier
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Your benchmark job has completed - possible performance regressions were detected. A full report can be found here. cc @jrevels

Please sign in to comment.