Opt-in sparse results for reductions along a dimension; reduce column-range views through the sparse kernels - #798

Fixes #43. Fixes #377.

sum(A; dims) and the other reductions along a dimension of a SparseMatrixCSC keep returning a dense Matrix, as they do today: the result has one dimension fewer and is usually dense, and downstream code expects it that way (see the discussion below). The sparse result asked for in #43 is opt-in by reducing into a sparse destination:

sum!(spzeros(size(A, 1), 1), A)                    # row sums, stored only for rows of A that store something
Base.mapreducedim!(abs2, +, spzeros(1, size(A, 2)), A)

This works for sum!, prod!, maximum!, minimum!, count!, any!, all! and Base.mapreducedim!. The result stores an entry only for the rows or columns that store one themselves, unless the reduction of a structurally empty slice is nonzero (as for mapreduce(x -> x + 1, +, A; dims)), in which case it is fully stored. A destination that stores only zeros, as a reused sum! destination does after its fill!, folds like an empty one, so preallocating once and calling sum! repeatedly stays on the fast path. The manual gets a short paragraph on this.

How. A _mapreducedim! for a sparse destination reduces a fully stored one as the dense array its values form, fills an empty one without touching the slices that store nothing, and sends anything in between through the existing element-wise kernel. Row reductions build the 1 x n result column by column. Column reductions use the result's value vector as a dense workspace compressed in place when there are enough stored entries, and otherwise sort the stored entries by row so that only rows storing something are visited.

Measured on the 1.14-DEV build, min over BenchmarkTools samples:

Views of a column range (fixes #377). sum(view(A, :, j:k); dims) and the other reductions of such a view went through Base's element-wise fallback, indexing the parent once per element, and returned a sparse result because Base's similar on a SubArray defers to the sparse parent. The kernels above only need the column pointers, row indices and values, which a SparseMatrixCSCView already exposes off the parent's storage, so they now accept SparseMatrixCSCUnion, and reducedim_initarray covers the view so it returns the same dense Matrix as its copy. nnz of such a view is O(1) instead of a pass over its columns. The issue's example, 101 x 99 at density 0.01, columns 5:11:

Tests compare every reduction against the dense result over shapes, densities and dims = 1, 2, (1, 2), 3, for the matrix, its column-range view (same type and values as the copy) and into sparse destinations; check the stored pattern, empty dimensions, the hypersparse path (an allocation bound on a 10^6-row matrix with a reused destination), and sum!/mapreducedim! into dense, empty, partially and fully stored destinations. The #27836 test asserts a Matrix again. The full suite passes locally on 1.14-DEV.

A local sweep over the exported in-place reductions, shapes, densities, index types, stored and negative zeros, column-range views and reused destinations found the opt-in results agreeing with dense up to summation order, and the default results matching dense in type and value. Two pre-existing corners it hit are left alone: fill!(spzeros(0, n), x) with x != 0 throws "step cannot be zero" in _fillnonzero!, which prod! or all! into a 0-row sparse destination reach, and minimum(spzeros(3, 0); dims = 1) returns a sparse 1 x 0 result through Base's map over an empty view.

🤖 Generated with Claude Code

https://claude.ai/code/session_01DGaTC2P39YrE2aAU5es55x