added sparse matrix inner product #27470
Conversation
ararslan
added
needs tests
|
Tests added. If anyone can, please cancel the two earlier jobs running on AppVeyor to let the last one go through. Travis seems to have done that automatically. edit: jobs can be removed on CircleCI as well |
stdlib/SparseArrays/src/linalg.jl
Outdated
| @@ -203,6 +203,28 @@ function spmatmul(A::SparseMatrixCSC{Tv,Ti}, B::SparseMatrixCSC{Tv,Ti}; | |||
| return C | |||
| end | |||
|
|
|||
| # Frobenius inner product: trace(A†B) | |||
There was a problem hiding this comment.
Why not just trace(A'*B)?
There was a problem hiding this comment.
A dagger, for the conjugate transpose. That's how I had it in my package, but I can change it to whatever notation is preferred (maybe ')
There was a problem hiding this comment.
It is ' in Julia, so that would make sense.
stdlib/SparseArrays/src/linalg.jl
Outdated
| for i1 = A.colptr[j]:A.colptr[j+1]-1 | ||
| ra = A.rowval[i1] | ||
| for i2 = B.colptr[j]:B.colptr[j+1]-1 | ||
| rb = B.rowval[i2] |
There was a problem hiding this comment.
I realise it's a bit pedantic, but maybe keep the indices consistent? i.e. ia instead of i1 and ib instead of i2.
stdlib/SparseArrays/src/linalg.jl
Outdated
| @inbounds for j = 1:n | ||
| for i1 = A.colptr[j]:A.colptr[j+1]-1 | ||
| ra = A.rowval[i1] | ||
| for i2 = B.colptr[j]:B.colptr[j+1]-1 |
There was a problem hiding this comment.
you shouldn't need to restart this iterator
There was a problem hiding this comment.
ah you are right, I'll improve the algorithm
There was a problem hiding this comment.
it's not the prettiest, but you could do something like:
ia = ia_nxt; ib = ib_nxt
ia_nxt = A.colptr[j+1]; ib_nxt = B.colptr[j+1]
if ia < ia_nxt && ib < ib_nxt
ra = A.rowval[ia]; rb = B.rowval[ib]
while true
if ra < rb
ia += 1
ia < ia_nxt || break
ra = A.rowval[ia];
elseif ra > rb
ib += 1
ib < ib_nxt || break
rb = B.rowval[ib]
else # ra == rb
r += A.nzval[ia]' * B.nzval[ib]
ia += 1; ib += 1
ia < ia_nxt && ib < ib_nxt || break
ra = A.rowval[ia]; rb = B.rowval[ib]
end
end
endThere was a problem hiding this comment.
I compared yours (which has a small typo on the first line) with one where I just updated the lower bound of the loop and yours was faster
stdlib/SparseArrays/src/linalg.jl
Outdated
| ra = A.rowval[ia]; rb = B.rowval[ib] | ||
| while true | ||
| if ra < rb | ||
| ia += 1 |
There was a problem hiding this comment.
This increment (and similar increments below) may yield type instability when the relevant SparseMatrixCSC's colptr type is not Int?
stdlib/SparseArrays/src/linalg.jl
Outdated
| function vecdot(A::SparseMatrixCSC{T1,S1},B::SparseMatrixCSC{T2,S2}) where {T1,T2,S1,S2} | ||
| m, n = size(A) | ||
| size(B) == (m,n) || throw(DimensionMismatch("matrices must have the same dimensions")) | ||
| r = zero(promote_type(T1,T2)) |
There was a problem hiding this comment.
Rather than promote_type(T1, T2), initialization of r should probably reflect the actual operation? (Edit: See, e.g., the generic implementation in stdlib/LinearAlgebra/src/generic.jl:651 .)
stdlib/SparseArrays/src/linalg.jl
Outdated
| ib < ib_nxt || break | ||
| rb = B.rowval[ib] | ||
| else # ra == rb | ||
| r += A.nzval[ia]' * B.nzval[ib] |
There was a problem hiding this comment.
Update of r should probably call (vec)dot? See, e.g., the generic implementation in stdlib/LinearAlgebra/src/generic.jl:651.
|
Comments addressed. I'd like to change |
fredrikekre
removed
needs news
stdlib/SparseArrays/src/linalg.jl
Outdated
| ra = A.rowval[ia]; rb = B.rowval[ib] | ||
| while true | ||
| if ra < rb | ||
| ia += one(S1) |
There was a problem hiding this comment.
Here you need an additive identity, so oneunit rather than one :).
(Edit: And likewise below.)
|
Barring style changes to the |
|
thanks! |
Provides a nice speedup relative to the generic method currently dispatched by
vecdot.Please review the style and tell me if anything needs to be changed.
Tests to be added andfinal name TBD pending resolution of #25565.