Implement sparse arrays for every back-end by maleadt · Pull Request #745 · JuliaGPU/GPUArrays.jl

maleadt · 2026-06-23T20:26:11Z

GPUArrays now implements sparse arrays itself. A back-end gets them by naming its storage. Before, GPUArrays only defined abstract types and hooks that every back-end had to implement.

using CUDA, SparseArrays, LinearAlgebra

A = cu(sprand(Float32, 1000, 1000, 0.01))  # CuSparseMatrixCSC{Float32, Int32}
x = CUDA.rand(Float32, 1000)

A * x            # cuSPARSE SpMV
2 .* A .+ A      # broadcasting stays sparse
A * A'           # sparse × sparse
A[1:10, :]       # indexing with ranges, slices and masks
sparse(cu(Int32[1, 2, 2]), cu(Int32[1, 1, 1]), cu(Float32[1, 2, 3]), 2, 2)  # assembly, combining duplicates
SparseMatrixCSC(A)  # back to the host

The same code works on MtlArray, ROCArray, oneArray, CLArray and JLArray, all running on the device.

Design

GPUSparseVector and GPUSparseMatrixCSR/CSC/COO are concrete structs whose fields are the back-end's dense vectors (GPUSparseMatrixCSR{Tv,Ti,Vi,Vv}). The same struct is passed to kernels, so there are no separate device types. Every operation has a portable KernelAbstractions implementation. A back-end names the types for its storage:

const CuSparseMatrixCSR{Tv,Ti} = GPUSparseMatrixCSR{Tv,Ti,<:CuVector{Ti},<:CuVector{Tv}}

To use a vendor library, it overrides LinearAlgebra's storage-level mul! for those aliases and falls back to GPUArrays.generic_spmv!/generic_spmm!/generic_spgemm, as CUDA.jl already does. The interface docs cover the contract, the supported operations, where semantics differ from SparseArrays, and migration.

Supported: products with dense and sparse arrays (including transposes, Symmetric, Hermitian, Diagonal and I), broadcasting and map, reductions, norm/opnorm/dot, ==, indexing, hcat/vcat, kron, exp, triu/tril/diag/tr, dropzeros!/fkeep!, reshape, conversions between all formats and to and from dense arrays, and sparse/sparsevec/spdiagm assembly. Not implemented: hash, blockdiag, and triangular solves without a vendor library.

Breaking

The abstract sparse types can't be subtyped anymore. AbstractGPUSparseMatrixCSR/CSC/COO remain as deprecated aliases of the concrete types, so dispatch on them keeps working.
The back-end hooks (dense_array_type, sparse_array_type, csc_type, csr_type, coo_type, _sptranspose, _spadjoint) and the GPUSparseDevice* types are removed.

None of the back-ends needs a breaking release: they keep their constructors' defaults (e.g. Cint indices on CUDA.jl and AMDGPU.jl) and keep old field names as deprecated properties.

Back-end PRs

Each branch ends with a TEMP commit that points CI at this branch.

Back-end	Branch	Release	Vendor library
CUDA.jl	`tb/sparse-v13`	6.6.0	cuSPARSE (SpMV/SpMM/SpGEMM, triangular solves, BSR)
AMDGPU.jl	`tb/sparse-v13`	2.12.0	rocSPARSE
oneAPI.jl	`tb/sparse-v13`	2.11.0	oneMKL
Metal.jl	JuliaGPU/Metal.jl#982	1.13.0	none (portable kernels)
OpenCL.jl	`tb/sparse-v13`	compat only	none

Validated on an RTX 5080, an AMD iGPU (gfx1036), an Intel Arc A770, an M1 and an M3 Pro, and with PoCL: each back-end's own sparse tests plus GPUArrays' shared sparse testsuite on its arrays.

Fixes #648, #655 and #677. Supersedes #684 and JuliaGPU/Metal.jl#909.

Tested here: JLArray/sparse (14,730 tests) and a check that the sparse code adds no method ambiguities, on Julia 1.10 and 1.13.

kshyatt · 2026-06-23T20:54:11Z

+
+A sparse array can't share the `AbstractGPUArray` supertype — that is a `DenseArray`,
+whereas a sparse array must be an `AbstractSparseArray` — so GPUArrays keeps a parallel
+sparse hierarchy with its own generic functionality. Integrating a back-end has three


Weird phrasing here, let me think a bit about how to make this more natural

kshyatt · 2026-06-23T20:54:45Z

+| `AbstractGPUSparseMatrixCSR{Tv,Ti}` | `rowPtr`, `colVal`, `nzVal`, `dims`, `nnz` |
+| `AbstractGPUSparseMatrixCOO{Tv,Ti}` | `rowInd`, `colInd`, `nzVal`, `dims`, `nnz` |
+
+The pointer/index/value arrays are the back-end's own dense vector type. Provide only the


Should we give an example (e.g. CuVector for CUDA)

kshyatt · 2026-06-23T20:55:21Z

+
+### Interface to implement
+
+  * **Constructors** — from component arrays (`MyCSR(rowPtr, colVal, nzVal, dims)`),


What about to/from dense arrays?

kshyatt · 2026-06-23T20:56:59Z

-JLArray(x::JLSparseVector)    = JLArray(collect(SparseVector(x)))
-JLArray(x::JLSparseMatrixCSC) = JLArray(collect(SparseMatrixCSC(x)))
-JLArray(x::JLSparseMatrixCSR) = JLArray(collect(SparseMatrixCSC(x)))
+JLArray(x::JLSparseVector)    = GPUArrays.to_dense(x)


should we make this an eval block?

kshyatt · 2026-06-23T20:58:50Z

 SparseArrays.SparseMatrixCSC(x::AbstractGPUSparseMatrixCSC) = SparseMatrixCSC(size(x)..., Array(SparseArrays.getcolptr(x)), Array(SparseArrays.rowvals(x)), Array(SparseArrays.nonzeros(x)))

+function check_sparse_target(::Type{ST}, ::Type{Tv}, ::Type{Ti}) where {ST,Tv,Ti}
+    Tv isa Type && Ti isa Type ||


This looks really funky to me. Shouldn't Ti also check it's a <: Integer?

kshyatt · 2026-06-23T20:59:08Z


+function check_sparse_target(::Type{ST}, ::Type{Tv}, ::Type{Ti}) where {ST,Tv,Ti}
+    Tv isa Type && Ti isa Type ||
+        throw(ArgumentError("sparse target type $ST must specify value and index eltypes"))


Should we include Tv and Ti in the error message?

kshyatt · 2026-06-23T21:01:44Z

+@kernel function densify_vector_kernel!(out, iPtr, nzVal)
+    k = @index(Global, Linear)
+    if k <= length(nzVal)
+        @inbounds out[iPtr[k]] = nzVal[k]


it's funny how Claude made this a one liner here but not above

maleadt · 2026-06-23T21:06:16Z

Yeah this hasn't passed my personal sniff test yet, hence draft, but thanks for the review. FWIW, the initial motivation was Metal support + SpMV/SpMM (which happens to be easier to express in COO), with the rest following out of it. The API simplification, replacing *_array_type / dense_array_type / sparse_array_type with something like similar feels like a genuine simplification though. (tried to get rid of the remaining csc_type/csr_type/... too, but that proved too hard)

kshyatt · 2026-06-23T21:12:10Z

Ah I can eff off until it's ready if you prefer, sorry!

gdalle · 2026-06-24T08:27:34Z

I'd love to review this one too once it's ready!

maleadt · 2026-07-07T13:36:44Z

Note to self: add a findnz regression test here to catch issues like JuliaGPU/CUDA.jl#3189

kshyatt · 2026-08-14T15:53:52Z

Can we resurrect this?

gdalle · 2026-09-29T16:22:40Z

This is exciting!!!

gdalle · 2026-11-06T09:22:22Z

Hey @maleadt, need any help to push this one over the finish line? It would be very cool indeed to have sparse linalg on a Metal GPU

maleadt · 2026-11-06T17:34:22Z

I'll spend some tokens on another rebase and update, but I won't have the time to push this over the line right now. It's definitely on the short-term TODO list though.

GPUArrays now implements the sparse formats itself: GPUSparseVector and GPUSparseMatrixCSR/CSC/COO are immutable structs parametrized on the dense vectors that store their indices and values. Generic code allocates results with similar() on the buffers it already has, so the trait hooks (dense_array_type, sparse_array_type, csc_type, csr_type, coo_type, _sptranspose, _spadjoint) are gone, and a back-end needs no sparse-specific code: it can alias the types for its storage and add a transfer policy for its adaptor. The same struct is used on the device, replacing the GPUSparseDevice* types and the per-back-end adapt rules. The formats are canonical, as in SparseArrays (COO in row-major order). Constructors check host metadata only; check_structure verifies the invariants on the device. similar, copy, copyto!, zero and adapt follow SparseArrays and Adapt semantics: copies and conversions never share buffers, and adapting to a storage type keeps the format and types. Conversions between all formats run on the device (generic_regroup for CSR <-> CSC, generic_expand/generic_compress for CSR <-> COO). Broadcasting and reductions are ported to the new types. Broadcasting no longer reuses the input's pointer buffer for its output, and handles matrices without rows or columns. AbstractGPUSparseMatrixCSR/CSC/COO remain as deprecated aliases of the concrete types for downstream dispatch. JLArrays keeps only aliases and constructors on them. The testsuite builds its inputs from dense buffers and covers every format for every back-end. The raw constructors borrow their buffers, like SparseArrays' raw constructors. Operations that change the structure of an array resize its buffers, and check that they can be resized before modifying anything, so arrays backed by views remain usable for reading and for updating values. dataids covers every buffer, so that copyto! and broadcasts unalias their arguments. Scalar indexing reads the buffers, so it requires allowscalar exactly when the storage does.

GPUArrays now implements the sparse formats generically, parametrized on their storage, so Metal only names them for MtlVector storage (MtlSparseVector, MtlSparseMatrixCSR/CSC/COO) and decides what mtl does with host sparse arrays: matrices become CSR, indices Int32 when they fit, and values follow the dense narrowing rules. Structural adapt(MtlArray, S) keeps the format and types without Metal code. This requires GPUArrays 13 (JuliaGPU/GPUArrays.jl#745).

Dense GPU arrays convert to every sparse format on the device, through findall in column-major order (or on the dense transpose for row-major formats), and sparse(A; fmt) selects the format like CUDA.jl and AMDGPU.jl do. copy(transpose(A)), copy(adjoint(A)) and permutedims produce the same format by regrouping, and the format constructors accept transposed operands. findnz returns the coordinates in column-major order like SparseArrays, for every format.

sparse(I, J, V, m, n, combine; fmt) and sparsevec(I, V, n, combine) on GPU vectors sort packed coordinate keys, combine the values of repeated coordinates with a per-run fold (reduce_by_key), and compress. As in SparseArrays, repeated entries are folded in input order, with + by default (| for Bool) and explicit zeros kept, so that a call gives the same result on every back-end. Because the device sort need not be stable, a second sort on keys that are unique per entry restores the input order, unless the order cannot change the result (+ and * on integers, max, min, |, &). spdiagm builds on the same path.

The sparse redesign removes the back-end hooks and changes the sparse types, which is breaking for back-ends. JLArrays requires the new version; its own breaking release (0.5, for KernelAbstractions 0.10) is not registered yet and covers the change of its sparse types.

Products hook into LinearAlgebra's storage-level mul!(C, tA, [tB,] A, B, α, β), which Julia 1.13 calls directly and which older versions reach through the forwarded generic_matvecmul!/generic_matmatmul!, like the dense products do. The generic implementations (generic_spmv!, generic_spmm!) gather over the CSR layout of op(A), or the CSC layout of op(B) for dense × sparse, so they are deterministic and work for every element type; Float16 accumulates in Float32. Transposed operands regroup the matrix first. Symmetric and Hermitian operands, which the old code treated as plain matrices, are computed as two passes over the used triangle and its transpose. LinearAlgebra allocates the destination of dense × sparse like the sparse operand, so * allocates a dense result for those.

The broadcast kernels iterate one compressed format, so arguments of other formats are converted to that of the first sparse matrix, COO matrices being read as CSR (which has the same order of entries) and a COO result converted back. map over sparse arrays is a broadcast, as in SparseArrays. SparseArrays implements +, - and map/broadcast of the basic operators for any AbstractSparseVector through scalar loops; GPUSparseVector gets its own methods.

fkeep!, dropzeros[!], droptol!, triu and tril filter the stored entries on the device and recompress; istriu, istril and isdiag check the nonzero values, and issymmetric/ishermitian compare with the transpose, which also accounts for stored zeros on one side. diag looks up each element with a binary search in its slice, reshape and vec reuse the column-major order of CSC, which reshaping preserves, and kron fills each column of the result from a pair of columns without sorting. The istriu/istril workaround for triangular wrappers of sparse matrices (JuliaLang/julia#55547) is carried over from the old code.

A[I, J] with ranges, colons and integers (giving a sparse vector for a row or column) and x[I] for vectors keep the stored entries in the selection and renumber them on the device. Ranges with a positive step preserve the order of the entries, so the result needs no sorting; other index vectors are rejected. The vector methods replace SparseArrays' own, which index through scalar loops.

dot between sparse vectors multiplies them with a broadcast, dot with a dense vector gathers its elements at the stored indices, and dot(x, A, y) is dot(x, A * y). count, any and all combine the stored values with the implicit zeros. float, complex, real and imag convert the stored values, sparse(A; fmt) and sparsevec(A) convert between sparse formats, and reductions of vectors support dims=1. Several of these replace SparseArrays' methods for any sparse vector, which loop over the elements on the host (Vector(x), imag).

generic_spgemm computes A * B for any pair of formats by expansion, sorting and compression: every product of a stored entry of A with one of B is materialized, each thread finding its entry of A by binary search over the running count of products, and the products are summed per coordinate through the assembly path, which keeps the result deterministic. The memory needed grows with the number of products. The sparse-output mul!(C, A, B, α, β) gives C the union of the structures of α*A*B and β*C, resizing its buffers, rather than assuming that an equal number of stored entries means an equal structure. exp is ported from the old code on top of these, with a fix for matrices whose norm is below one, which it scaled up instead of down.

…enation SparseArrays implements these for any AbstractSparseArray with scalar loops, which a GPU sparse array cannot run. == compares values, so stored zeros do not matter. Broadcasting into a sparse array materializes the result and copies it over, so the destination takes the structure of the result. Arithmetic with transposed, Symmetric and Hermitian wrappers, UniformScaling and Diagonal materializes the wrapper and combines the operands on the device, and hcat/vcat of sparse arrays offset the coordinates and reassemble them. blockdiag throws instead of falling back to scalar indexing.

maleadt · 2026-11-07T07:34:14Z

@gdalle May be worth testing this branch with your favorite back-end (ideally one that doesn't come with sparse array abstractions out of the box, e.g., Metal.jl or so) and test whether the functionality covers what you need.

gdalle · 2026-11-08T18:42:29Z

Thanks! I'll take a look

kshyatt reviewed Jun 23, 2026

View reviewed changes

michel2323 mentioned this pull request Sep 24, 2026

Make oneMKL sparse matrices AbstractGPUSparseArrays JuliaGPU/oneAPI.jl#642

Merged

This was referenced Sep 25, 2026

Compact sparse-vector broadcast output on the device #795

Merged

Fix sparse mapreduce for operators other than + #796

Merged

maleadt force-pushed the tb/sparse branch from 5115577 to 2a50d28 Compare September 25, 2026 18:23

maleadt changed the base branch from main to tb/permutedims-empty September 25, 2026 18:23

maleadt added this pull request to stack #797 September 25, 2026 18:23

This was referenced Sep 25, 2026

Don't launch a kernel to permute empty arrays #798

Merged

Sparse arrays JuliaGPU/Metal.jl#982

Draft

maleadt changed the title ~~Generic sparse arrays: CSC/CSR/COO formats, on-device conversions, SpMM/SpMV~~ Generic GPU sparse arrays: storage-parametric formats owned by GPUArrays Sep 25, 2026

christiangnrd mentioned this pull request Oct 1, 2026

Add GPUArrays sparse array support JuliaGPU/Metal.jl#909

Open

maleadt added the arrays The array abstraction: broadcasting, mapreduce, indexing, views, generic linear algebra. label Oct 19, 2026

This was referenced Oct 19, 2026

fixes the kron implementation for sparse + diagonal matrix JuliaGPU/CUDA.jl#2804

Closed

Fix sparse mul! with integer inputs and in-place CSC products JuliaGPU/CUDA.jl#3416

Merged

maleadt force-pushed the tb/sparse branch from 2a50d28 to 7d2123d Compare October 21, 2026 12:10

maleadt removed this pull request from stack #797 October 21, 2026 13:00

maleadt changed the base branch from tb/permutedims-empty to main October 21, 2026 13:00

maleadt force-pushed the tb/sparse branch from 65ff6ec to dd0ee48 Compare October 21, 2026 15:18

maleadt force-pushed the tb/sparse branch 5 times, most recently from 9a9126f to 76600f3 Compare October 21, 2026 18:27

maleadt force-pushed the tb/sparse branch from 76600f3 to bcd4482 Compare November 6, 2026 19:01

maleadt force-pushed the tb/sparse branch 4 times, most recently from 4cb19d3 to e2eb984 Compare November 6, 2026 20:03

maleadt added 10 commits November 6, 2026 21:15

maleadt force-pushed the tb/sparse branch from e2eb984 to 7f07bb1 Compare November 6, 2026 20:15

maleadt changed the title ~~Generic GPU sparse arrays: storage-parametric formats owned by GPUArrays~~ Implement sparse arrays for every back-end Nov 6, 2026


		### Interface to implement

		* Constructors — from component arrays (`MyCSR(rowPtr, colVal, nzVal, dims)`),

Conversation

maleadt commented Jun 23, 2026 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Design

Breaking

Back-end PRs

Uh oh!

kshyatt Jun 23, 2026

Choose a reason for hiding this comment

Uh oh!

kshyatt Jun 23, 2026

Choose a reason for hiding this comment

Uh oh!

kshyatt Jun 23, 2026

Choose a reason for hiding this comment

Uh oh!

kshyatt Jun 23, 2026

Choose a reason for hiding this comment

Uh oh!

kshyatt Jun 23, 2026

Choose a reason for hiding this comment

Uh oh!

kshyatt Jun 23, 2026

Choose a reason for hiding this comment

Uh oh!

kshyatt Jun 23, 2026

Choose a reason for hiding this comment

Uh oh!

maleadt commented Jun 23, 2026

Uh oh!

kshyatt commented Jun 23, 2026

Uh oh!

gdalle commented Jun 24, 2026

Uh oh!

maleadt commented Jul 7, 2026

Uh oh!

kshyatt commented Aug 14, 2026

Uh oh!

gdalle commented Sep 29, 2026

Uh oh!

gdalle commented Nov 6, 2026

Uh oh!

maleadt commented Nov 6, 2026

Uh oh!

maleadt commented Nov 7, 2026

Uh oh!

gdalle commented Nov 8, 2026

Uh oh!

This branch has not been deployed

Reviewers

Uh oh!

Assignees

Labels

Projects

Milestone

Development

Uh oh!

3 participants

maleadt commented Jun 23, 2026 •

edited

Loading