# CSC kills the prospect of multithreading. Shouldn't Julia use CSR?
**URL:** https://discourse.julialang.org/t/csc-kills-the-prospect-of-multithreading-shouldnt-julia-use-csr/102491
**Category:** Julia at Scale
**Tags:** sparse
**Created:** [August 4, 2023, 3:22pm UTC](https://discourse.julialang.org/t/csc-kills-the-prospect-of-multithreading-shouldnt-julia-use-csr/102491 "2023-08-04T15:22:26Z")
**Posts on this page:** 20
**Page:** 2
### Author:  [@nvenkov1](https://discourse.julialang.org/u/nvenkov1)
#### Post date: [August 5, 2023, 1:32pm UTC](https://discourse.julialang.org/t/csc-kills-the-prospect-of-multithreading-shouldnt-julia-use-csr/102491/21 "2023-08-05T13:32:51Z")
I find it strange that a community who develops a language with a reputation based on performance finds no “significant advantage” of having built-in multithreaded SpMV. Why should dense BLAS be multithreaded and not sparse BLAS? I doubt everyone thinks like you, especially the users whose work is focused on scientific computing.
---
### Author:  [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)
#### Post date: [August 5, 2023, 1:35pm UTC](https://discourse.julialang.org/t/csc-kills-the-prospect-of-multithreading-shouldnt-julia-use-csr/102491/22 "2023-08-05T13:35:21Z")
No, you seem to be missing @mkitti’s point. What is the performance advantage of a built-in library vs a package?
While I think it would be nice for SparseMatricesCSR.jl to be upstreamed into SparseArrays.jl, being in a separate package doesn’t magically change its performance.
---
### Author:  [@nvenkov1](https://discourse.julialang.org/u/nvenkov1)
#### Post date: [August 5, 2023, 1:37pm UTC](https://discourse.julialang.org/t/csc-kills-the-prospect-of-multithreading-shouldnt-julia-use-csr/102491/23 "2023-08-05T13:37:00Z")
Yes, I do agree that making it built-in does not mean performance gain. But it would be nice to have a seamless access to multithreaded SpMV, just like we do with dense BLAS functions.
---
### Author:  [@mkitti](https://discourse.julialang.org/u/mkitti)
#### Post date: [August 5, 2023, 1:38pm UTC](https://discourse.julialang.org/t/csc-kills-the-prospect-of-multithreading-shouldnt-julia-use-csr/102491/24 "2023-08-05T13:38:35Z")
What does built-in mean to you?
To me it means slow development of a package that is difficult to change and update.
---
### Author:  [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)
#### Post date: [August 5, 2023, 1:43pm UTC](https://discourse.julialang.org/t/csc-kills-the-prospect-of-multithreading-shouldnt-julia-use-csr/102491/25 "2023-08-05T13:43:49Z")
Note that SparseArrays is actually not a built-in anymore anyways, so it’s kind of funny that this is the library being discussed. Here is its repo:
> **[GitHub - JuliaSparse/SparseArrays.jl: SparseArrays.jl is a Julia stdlib](https://github.com/JuliaSparse/SparseArrays.jl)**
>
> SparseArrays.jl is a Julia stdlib. Contribute to JuliaSparse/SparseArrays.jl development by creating an account on GitHub.
It was made separate because it pulls in SuiteSparse, so people wanted the option to remove it from the standard build in order to have a GPL-free distribution. SparseArrays couldn’t be fully removed without breaking, so now it’s a standard library that lives in a separate Github repo in a separate organization but is built into the standard image, set to be fully removed (and then just become a normal Julia library) if a breaking change ever happens.
And this is the library that we’re discussing about whether a feature needs to be moved into in order to be “built-in” vs “a package”. I think that’s really splitting hairs. As long as there’s a good easily available library for SparseMatrixCSR I think it’s fine. I think organizationally and documentation wise, having it in SparseArrays would make more sense, and the discoverability of it in GridAp’s org instead of JuliaSparse is a bit odd and IMO should be changed, but there’s nothing about the package itself that is an issue from what I can tell.
---
### Author:  [@nvenkov1](https://discourse.julialang.org/u/nvenkov1)
#### Post date: [August 5, 2023, 1:44pm UTC](https://discourse.julialang.org/t/csc-kills-the-prospect-of-multithreading-shouldnt-julia-use-csr/102491/26 "2023-08-05T13:44:05Z")
If CSR support and multithreaded SpMV were part of SparseArrays.jl it would be easier to not have to use several packages, and it would improve the visibility of those functionalities as SparseArrays.jl is way more popular and used than ThreadedSparseArrays.jl and SparseMatricesCSR.jl. I’m not really conscious of the additional burden it would be, but I would surely appreciate it.
---
### Author:  [@nvenkov1](https://discourse.julialang.org/u/nvenkov1)
#### Post date: [August 5, 2023, 1:45pm UTC](https://discourse.julialang.org/t/csc-kills-the-prospect-of-multithreading-shouldnt-julia-use-csr/102491/27 "2023-08-05T13:45:14Z")
Oh… I didn’t know this.
---
### Author:  [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)
#### Post date: [August 5, 2023, 1:46pm UTC](https://discourse.julialang.org/t/csc-kills-the-prospect-of-multithreading-shouldnt-julia-use-csr/102491/28 "2023-08-05T13:46:33Z")
The minimum change to fix this is probably just to add a reference to SparseMatrixCSR.jl in the SparseArrays documentation, until people get around to the reorganization (though I’m sure that just by putting that reference into the docs it’ll spark the conversation of reorganization and get it done).
---
### Author:  [@stevengj](https://discourse.julialang.org/u/stevengj)
#### Post date: [August 5, 2023, 4:41pm UTC](https://discourse.julialang.org/t/csc-kills-the-prospect-of-multithreading-shouldnt-julia-use-csr/102491/31 "2023-08-05T16:41:10Z")
> [@nvenkov1](#):
>
> It is not possible to reorder an array for free. I cannot believe several people are convinced of such nonsense and try to disseminate this belief
If you have a matrix stored in CSC, I agree that you cannot re-order it for free.
What I said was that you can _construct_ the _transpose_ matrix as CSC initially, and then re-interpret it as CSR via the transpose wrapper. In this way you can use the existing SparseArrays library to construct CSR matrices with no overhead.A thin wrapper could make this API nicer.
---
### Author:  [@nvenkov1](https://discourse.julialang.org/u/nvenkov1)
#### Post date: [August 5, 2023, 5:07pm UTC](https://discourse.julialang.org/t/csc-kills-the-prospect-of-multithreading-shouldnt-julia-use-csr/102491/32 "2023-08-05T17:07:57Z")
Okay. Say all you have is a CSR representation of your matrix. That is, you are equipped with `m::Int`, `n::Int`, `rowptr::Vector{Ti}`, `colval::Vector{Ti}` and `nzval::Vector{Tv}`. If you want to make use of Julia’s built-in CSC functionalities to carry the equivalent of a CSR SpMV with CSC SpMV\_T, you can indeed proceed as follows: `B=SparseMatrixCSC(m, n, rowptr, colval, nzval); transpose(B)*x`. If that’s what you meant, I agree. Then it wouldn’t be a bad idea to parallelize SpMV\_T in SparseArrays.jl.
---
### Author:  [@Oscar\_Smith](https://discourse.julialang.org/u/Oscar_Smith)
#### Post date: [August 5, 2023, 5:15pm UTC](https://discourse.julialang.org/t/csc-kills-the-prospect-of-multithreading-shouldnt-julia-use-csr/102491/33 "2023-08-05T17:15:19Z")
The point is that any sparse matrix is constructed at some point so if you want a CSR rather than a CSC you can always construct the transpose of the matrix you want.
---
### Author:  [@bertschi](https://discourse.julialang.org/u/bertschi)
#### Post date: [August 5, 2023, 5:25pm UTC](https://discourse.julialang.org/t/csc-kills-the-prospect-of-multithreading-shouldnt-julia-use-csr/102491/34 "2023-08-05T17:25:45Z")
The point is that the type `Transpose{SparseMatrixCSC}` has exactly the same memory layout as CSR, i.e., knowing that you can directly construct a matrix accordingly. Here is a simple example:
```julia
julia> I = [1,2,2,3]; J = [4,1,3,1]; V = rand(4);
julia> A = sparse(I, J, V) # construct CSC from COO
3×4 SparseMatrixCSC{Float64, Int64} with 4 stored entries:
⋅ ⋅ ⋅ 0.0265797
0.502663 ⋅ 0.836581 ⋅
0.588429 ⋅ ⋅ ⋅
julia> B = transpose(sparse(J, I, V)) # construct CSR from COO
3×4 LinearAlgebra.Transpose{Float64, SparseMatrixCSC{Float64, Int64}} with 4 stored entries:
⋅ ⋅ ⋅ 0.0265797
0.502663 ⋅ 0.836581 ⋅
0.588429 ⋅ ⋅ ⋅
# Now, get the details of the CSR representation
julia> m, n = size(A);
julia> rowptr = B.parent.colptr; # Rows are stored as columns in transpose
julia> colval = B.parent.rowval; # Same for cols
julia> nzval = B.parent.nzval
# Reconstruct directly from CSR information
julia> C = transpose(SparseMatrixCSC(n, m, rowptr, colval, nzval))
3×4 LinearAlgebra.Transpose{Float64, SparseMatrixCSC{Float64, Int64}} with 4 stored entries:
⋅ ⋅ ⋅ 0.0265797
0.502663 ⋅ 0.836581 ⋅
0.588429 ⋅ ⋅ ⋅
```
Hope that helps.
---
### Author:  [@nvenkov1](https://discourse.julialang.org/u/nvenkov1)
#### Post date: [August 5, 2023, 5:42pm UTC](https://discourse.julialang.org/t/csc-kills-the-prospect-of-multithreading-shouldnt-julia-use-csr/102491/35 "2023-08-05T17:42:49Z")
I’m not sure to follow the point you’re making. The argument was about whether there is a cost associated with the workaround to do CSR-like SpMVs. There seem to be three cases.
First, we have `A::SparseMatrixCSC{Tv,Ti}`, in which case you can do a CSR-like SpMV by doing `B=sparse(transpose(A)); transpose(B)*x`. Then, you are paying a cost when doing `B=sparse(transpose(A))`.
Second, if you know the COO representation of your array with `I::Vector{Ti}`, `J::Vector{Ti}` and `V::Vector{Tv}`, you can do `B=sparse(J, I, V); transpose(B)*b`, in which case you are paying a cost when doing `B=sparse(J, I, V)` as it induces a conversion from COO to CSC.
Finally, if you have a CSR representation of your matrix in the form of `m::Int`, `n::Int`, `rowptr::Vector{Ti}`, `colval::Vector{Ti}`, `nzval::Vector{Tv}`, you can do `B=SparseMatrixCSC(m, n, rowptr, colval, nzval); transpose(B)*x`, in which case there is no other cost than SpMV\_T.
---
### Author:  [@bertschi](https://discourse.julialang.org/u/bertschi)
#### Post date: [August 5, 2023, 5:59pm UTC](https://discourse.julialang.org/t/csc-kills-the-prospect-of-multithreading-shouldnt-julia-use-csr/102491/36 "2023-08-05T17:59:49Z")
The point is that you could just define `SparseMatrixCSR` as a type alias for `Transpose{SparseMatrixCSC}`, together with some convenience constructors. There is no **additional** cost involved when constructing, i.e., CSR from COO has similar cost as CSC from COO and CSR from `rowptr, colval` and `nzval` is just like CSC from `colptr, rowval` and `nzval` as my example shows.
---
### Author:  [@mkitti](https://discourse.julialang.org/u/mkitti)
#### Post date: [August 5, 2023, 6:32pm UTC](https://discourse.julialang.org/t/csc-kills-the-prospect-of-multithreading-shouldnt-julia-use-csr/102491/37 "2023-08-05T18:32:58Z")
> [@nvenkov1](#):
>
> Finally, if you have a CSR representation of your matrix in the form of `m::Int`, `n::Int`, `rowptr::Vector{Ti}`, `colval::Vector{Ti}`, `nzval::Vector{Tv}`, you can do `B=SparseMatrixCSC(m, n, rowptr, colval, nzval); transpose(B)*x`, in which case there is no other cost than SpMV\_T.
I think we’re on the same page. Let me see try to unpack it.
Suppose I have the following CSR representation.
```julia
julia> A
5×5 SparseMatrixCSR{1, Int8, Int64} with 10 stored entries:
[1, 4] = 78
[1, 5] = -87
[2, 1] = 88
[2, 3] = 38
[2, 5] = 77
[3, 5] = 31
[4, 1] = -93
[5, 1] = -25
[5, 2] = -20
[5, 3] = 23
julia> Matrix(A)
5×5 Matrix{Int8}:
0 0 0 78 -87
88 0 38 0 77
0 0 0 0 31
-93 0 0 0 0
-25 -20 23 0 0
```
Construct the CSC matrix, `B`, as follows.
```julia
julia> B = SparseMatrixCSC(A.m, A.n, A.rowptr, A.colval, A.nzval)
5×5 SparseMatrixCSC{Int8, Int64} with 10 stored entries:
⋅ 88 ⋅ -93 -25
⋅ ⋅ ⋅ ⋅ -20
⋅ 38 ⋅ ⋅ 23
78 ⋅ ⋅ ⋅ ⋅
-87 77 31 ⋅ ⋅
julia> B_T = transpose(B)
5×5 Transpose{Int8, SparseMatrixCSC{Int8, Int64}} with 10 stored entries:
⋅ ⋅ ⋅ 78 -87
88 ⋅ 38 ⋅ 77
⋅ ⋅ ⋅ ⋅ 31
-93 ⋅ ⋅ ⋅ ⋅
-25 -20 23 ⋅ ⋅
```
Then we can reconstruct the matrix `A` in CSR form in the following way.
```julia
julia> A2 = SparseMatrixCSR(B_T)
5×5 SparseMatrixCSR{1, Int8, Int64} with 10 stored entries:
[1, 4] = 78
[1, 5] = -87
[2, 1] = 88
[2, 3] = 38
[2, 5] = 77
[3, 5] = 31
[4, 1] = -93
[5, 1] = -25
[5, 2] = -20
[5, 3] = 23
```
That just expands to the following.
```julia
julia> A3 = SparseMatrixCSR{1}(B.m, B.n, B.colptr, B.rowval, B.nzval)
5×5 SparseMatrixCSR{1, Int8, Int64} with 10 stored entries:
[1, 4] = 78
[1, 5] = -87
[2, 1] = 88
[2, 3] = 38
[2, 5] = 77
[3, 5] = 31
[4, 1] = -93
[5, 1] = -25
[5, 2] = -20
[5, 3] = 23
julia> A === A2 && A === A3
true
```
Then we calculate `A*x` via the following function which I think you would refer to as “CRC SpMV”.
>
```julia
function mul!(y::AbstractVector,A::SparseMatrixCSR,v::AbstractVector)
A.n == size(v, 1) || throw(DimensionMismatch())
A.m == size(y, 1) || throw(DimensionMismatch())
fill!(y, zero(eltype(y)))
o = getoffset(A)
for row = 1:size(y, 1)
@inbounds for nz in nzrange(A,row)
col = A.colval[nz]+o
y[row] += A.nzval[nz]*v[col]
end
end
return y
end
```
Are we saying the same thing?
---
### Author:  [@nvenkov1](https://discourse.julialang.org/u/nvenkov1)
#### Post date: [August 6, 2023, 12:18pm UTC](https://discourse.julialang.org/t/csc-kills-the-prospect-of-multithreading-shouldnt-julia-use-csr/102491/38 "2023-08-06T12:18:46Z")
No. CSR from COO, CSR from CSC and CSC from COO all **do have a cost** associated with them which CSR from `rowptr`, `colval`, `nzval` and CSC from `colptr`, `rowval`, `nzval` don’t.
---
### Author:  [@bertschi](https://discourse.julialang.org/u/bertschi)
#### Post date: [August 6, 2023, 12:23pm UTC](https://discourse.julialang.org/t/csc-kills-the-prospect-of-multithreading-shouldnt-julia-use-csr/102491/39 "2023-08-06T12:23:31Z")
I understand, but first, the cost for CSR from COO and CSC from COO are the same, and, secondly – and probably more important to you – CSR from `rowptr`, `colval`, `nzval` and CSC from `colptr`, `rowval`, `nzval` also have the same (much lower) cost.
Now, using the representation that I showed, the last call `transpose(SparseMatrixCSC(n, m, rowptr, colval, nzval))` in my example constructs a CSR representation directly from `rowptr`, `colval`, `nzval` and without any additional cost (compared to constructing a CSC representation from `colptr`, `rowval`, `nzval`).
---
### Author:  [@nvenkov1](https://discourse.julialang.org/u/nvenkov1)
#### Post date: [August 6, 2023, 12:41pm UTC](https://discourse.julialang.org/t/csc-kills-the-prospect-of-multithreading-shouldnt-julia-use-csr/102491/40 "2023-08-06T12:41:13Z")
Yes, because at that stage you had the data structures of the CSR assembled, i.e., `rowptr`, `colval` and `nzval`. But if your starting point is a CSC matrix `A::SparseMatrixCSC{Tv,Ti}` or COO representation given by `I::Vector{Ti}`, `J::Vector{Ti}` and `V::Vector{Tv}`, like at the beginnig of your example, then you have some work to do to get to your CSR.
---
### Author:  [@bertschi](https://discourse.julialang.org/u/bertschi)
#### Post date: [August 6, 2023, 12:54pm UTC](https://discourse.julialang.org/t/csc-kills-the-prospect-of-multithreading-shouldnt-julia-use-csr/102491/41 "2023-08-06T12:54:59Z")
Ok, but as others have already explained as well, you have to construct your sparse matrix somehow in any case. So, why not start with `rowptr`, `colval` and `nzval`, if that’s what you want in the first place?
The example was just an illustration, I have no idea where the sparse matrix comes from in your case.
---
### Author:  [@Barget](https://discourse.julialang.org/u/Barget)
#### Post date: [August 6, 2023, 1:27pm UTC](https://discourse.julialang.org/t/csc-kills-the-prospect-of-multithreading-shouldnt-julia-use-csr/102491/42 "2023-08-06T13:27:59Z")
This thread is interesting, and made me learn a few things about CSR vs CSC.
@nvenkov1 , could you take a step back a talk a bit about your workflow? It’s not clear to me what kind of input do you have (COO? CSR? CSC? other? …), and what are your requirements.
In particular, why do you need your starting point to be CSR? (just to understand)
If mandatory (because … why not), then considering the memory layout is equivalent to `transpose{SparseMatrixCSC}`, why not just “cast” your data to that type?
Granted, as already mention, this would be easier with a thin wrapper that could worth a PR. But all in all, nothing insurmontable or critical, isn’t it?
[Previous page](https://discourse.julialang.org/t/csc-kills-the-prospect-of-multithreading-shouldnt-julia-use-csr/102491.md?page=1)
[Next page](https://discourse.julialang.org/t/csc-kills-the-prospect-of-multithreading-shouldnt-julia-use-csr/102491.md?page=3)