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.

Sign up for GitHub

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

stdlib/SparseArrays: add rmul! and lmul! of sparse matrix with Diagonal #29296

Merged
merged 5 commits into from
Sep 29, 2018
Merged

stdlib/SparseArrays: add rmul! and lmul! of sparse matrix with Diagonal #29296

Changes from 2 commits
Commits
Show all changes
5 commits
Select commit Hold shift + click to select a range
0537f6b
Add rmul! and lmul! of sparse matrix with Diagonal
dkarrasch Sep 21, 2018
b93376c
Merge pull request #1 from dkarrasch/dk_rmul_lmul_sp_diag
dkarrasch Sep 21, 2018
43af178
Some optimization
dkarrasch Sep 21, 2018
2efc646
merge loops again for consistency
dkarrasch Sep 21, 2018
2edf05f
remove unused parametrized type
dkarrasch Sep 28, 2018
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
22 changes: 22 additions & 0 deletions stdlib/SparseArrays/src/linalg.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -1017,11 +1017,33 @@ function rmul!(A::SparseMatrixCSC, b::Number)
rmul!(A.nzval, b)
return A
end

function lmul!(b::Number, A::SparseMatrixCSC)
lmul!(b, A.nzval)
return A
end

function rmul!(A::SparseMatrixCSC, D::Diagonal{T}) where T
Copy link
Contributor

Choose a reason for hiding this comment

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

T is never used in the function body so I think you can make this signature

function rmul!(A::SparseMatrixCSC, D::Diagonal). (same for lmul)

m, n = size(A)
(n == size(D, 1)) || throw(DimensionMismatch())
Anzval = A.nzval
for col = 1:n, p = A.colptr[col]:(A.colptr[col+1]-1)
@inbounds Anzval[p] *= D.diag[col]
Copy link
Member

Choose a reason for hiding this comment

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

It might be possible to gain a bit of speed here by manually hoisting the D.diag[col] lookup outside of the inner loop. Of course, the compiler may well be able to do that automatically. Depends on how much it knows about whether D.diag can change or move.

Copy link
Member

Choose a reason for hiding this comment

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

From my experience with the SparseArray code, this is only needed when the object one getfield into is mutable. There is a lot of hoisting out fields for SparseMatrixCSC that was put there before #16371 (see discussion in #15668) but it shouldn't be needed now. Always worth benchmarking though...

Copy link
Member Author

Choose a reason for hiding this comment

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

I'm sorry I can't follow. I did the change @StefanKarpinski suggested. Was that just (potentially) unnecessary or even "wrong"?

Copy link
Member Author

Choose a reason for hiding this comment

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

I benchmarked the previous ("my") against the latest version ("Stefan's") on these:

P = sprand(12000, 12000, 0.01)
q = Diagonal(dropdims(sum(P, dims=2), dims=2))

and Stefan's won with 607.996 μs vs. 631.496 μs. So I guess it was worth it.

Copy link
Member Author

Choose a reason for hiding this comment

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

Hm, might be neglible. For

P = sprand(120000, 120000, 0.01)

(an order of magnitude larger) and q as before I get 104.832 ms vs. 103.934 ms.

Copy link
Member

Choose a reason for hiding this comment

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

The question is if the compiler tries to reload D.diag` every time in the loop or if it realizes that this will never change and hoist the load outside the loop.

Copy link
Member Author

Choose a reason for hiding this comment

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

Seems like it doesn't reload, at least it doesn't show up in performance loss. Should I, for better readability, then merge the two loops again as before, but leave the @inbounds upfront like in the last commit?

Copy link
Member

Choose a reason for hiding this comment

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

Doesn't matter that much but perhaps it makes sense to at least have the two methods in this PR be consistent with eachother.

end
return A
end

function lmul!(D::Diagonal{T}, A::SparseMatrixCSC) where T
m, n = size(A)
(m == size(D, 2)) || throw(DimensionMismatch())
Anzval = A.nzval
Arowval = A.rowval
for col = 1:n, p = A.colptr[col]:(A.colptr[col+1]-1)
@inbounds Anzval[p] *= D.diag[Arowval[p]]
end
return A
end

function \(A::SparseMatrixCSC, B::AbstractVecOrMat)
@assert !has_offset_axes(A, B)
m, n = size(A)
Expand Down