From a1924727776c5e8176120e8f7a488f86ab64971f Mon Sep 17 00:00:00 2001 From: "Viral B. Shah" Date: Tue, 20 Oct 2026 20:51:43 +0200 Subject: [PATCH 01/30] Skip the dense copy of the sparse arguments in a sparse-plus-scalar broadcast such as `A .+ 1` (#1105) MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Since #1004, a broadcast sum of a sparse array and a scalar, such as `A .+ 1`, is dense. It gets there by copying every sparse argument into a dense array, then broadcasting over the copies into a second dense array. For a sparse argument with few stored entries, that copy is half the memory and more than half the time. **Before** (`main`; random `nร—n` sparse `A`, `B` of density `d`; minimum time over 20 runs) | n, d | `A .+ 1` | `A .+ B .- 2` | |---|---|---| | 1000, 0.001 | 2.11 ms, 15.3 MiB | 3.67 ms, 22.9 MiB | | 1000, 0.01 | 2.20 ms, 15.3 MiB | 3.65 ms, 22.9 MiB | | 5000, 0.001 | 34.3 ms, 381.5 MiB | 60.5 ms, 572.2 MiB | | 5000, 0.01 | 34.8 ms, 381.5 MiB | 62.4 ms, 572.2 MiB | | 2000, 0.03 | 6.04 ms, 61.0 MiB | 9.62 ms, 91.6 MiB | | 2000, 0.5 | 7.42 ms, 61.0 MiB | 12.4 ms, 91.6 MiB | **After** | n, d | `A .+ 1` | `A .+ B .- 2` | |---|---|---| | 1000, 0.001 | 0.87 ms, 7.6 MiB | 0.96 ms, 7.6 MiB | | 1000, 0.01 | 1.05 ms, 7.6 MiB | 3.67 ms, 22.9 MiB | | 5000, 0.001 | 14.6 ms, 190.7 MiB | 16.2 ms, 190.7 MiB | | 5000, 0.01 | 22.2 ms, 190.7 MiB | 61.9 ms, 572.2 MiB | | 2000, 0.03 | 3.85 ms, 30.5 MiB | 9.51 ms, 91.6 MiB | | 2000, 0.5 | 7.45 ms, 61.0 MiB | 12.3 ms, 91.6 MiB | Sums with a dense array, such as `A .+ M` and `v .+ M`, are unchanged. The new path evaluates the broadcast once with each sparse argument replaced by its implicit zero. That fills every entry of the result that no stored entry reaches. It then evaluates the broadcast again only at the entries a stored entry is broadcast to. The function is a tree of `+` and `-`, so evaluating it twice at an entry is unobservable. Results match dense broadcast in type and `isequal`, including NaN and `-0.0`. I checked that against dense on 36 expressions over Float64, ComplexF64, Bool and Int entries, with sparse vectors and 1-row and 1-column matrices broadcast across the result. Each re-evaluated entry looks up every sparse argument. So the new path runs only while the entries reached, times the number of sparse arguments, stay below 1/32 of the result; otherwise the dense copy is cheaper. Above that density, as in the 0.5 row and for `A .+ B .- 2` at 0.01, the old path runs. Left out: - Sums with a dense array argument. Measured near the cutoff, the copy was not their main cost, and the new path was slower. - Adjoints, transposes and views of sparse arrays, which are still copied. - Element types that are not `Number`, and abstract ones such as `Real`. Their implicit zero can have another type than the result that dense broadcast infers from the values, so the old path runs whenever the two types differ. ๐Ÿค– Generated with [Claude Code](https://claude.com/claude-code) --------- Co-authored-by: Claude Opus 5.5 (cherry picked from commit ed2edcce78d4523772a3ba74d5b0f71fb06e6353) --- src/higherorderfns.jl | 76 ++++++++++++++++++++++++++++++++++++++++-- test/higherorderfns.jl | 8 +++++ 2 files changed, 82 insertions(+), 2 deletions(-) diff --git a/src/higherorderfns.jl b/src/higherorderfns.jl index 3f47cdaa..be693c63 100644 --- a/src/higherorderfns.jl +++ b/src/higherorderfns.jl @@ -1216,7 +1216,7 @@ end # broadcast entry points for combinations of sparse arrays and other (scalar) types @inline function copy(bc::Broadcasted{<:SparseVecOrMatStyle}) - _isdensesum(bc) && return _densebroadcast(bc) + _isdensesum(bc) && return _densesum(bc) bcf = flatten(bc) _asksthezeros(bc) && _fillszeros(bcf.f, bcf.args...) && return _densenonzero(bc, bcf.f, bcf.args...) return _copy(bcf.f, bcf.args...) @@ -1360,7 +1360,7 @@ broadcast(f::Tf, A::AbstractSparseMatrixCSC, ::Type{T}) where {Tf,T} = broadcast # which is sparse, and inferable, for every element type. function copy(bc::Broadcasted{PromoteToSparse}) - _isdensesum(bc) && return _densebroadcast(bc) + _isdensesum(bc) && return _densesum(bc) bcf = flatten(bc) _asksthezeros(bc) && _fillszeros(bcf.f, bcf.args...) && return _densenonzero(bc, bcf.f, bcf.args...) if can_skip_sparsification(bcf.f, bcf.args...) @@ -1382,6 +1382,78 @@ end end _isdensesum(bc::Broadcasted) = _issumtree(bc) && _anydensifying(bc.args...) + +# A sum of sparse arrays and scalars with a dense result, such as `A .+ 1`, skips the dense +# copy of each sparse argument. It broadcasts once with every sparse argument replaced by its +# implicit zero, which is already right wherever no stored entry reaches, then recomputes the +# entries that a stored entry reaches. The function is made of `+` and `-`, so computing an +# entry twice is harmless. The rest goes to `_densebroadcast`: sums with a dense array, where +# the copy costs little, non-`Number` elements, whose zero may be of another type, and results +# whose element type depends on the values. +function _densesum(bc::Broadcasted) + ax = axes(bc) + (any(isempty, ax) || !_sparseandscalars(bc)) && return _densebroadcast(bc) + # recomputing an entry looks up every sparse argument, which costs about as much as copying + # 32 entries, so beyond that many reached entries the dense copy is cheaper + len = prod(map(length, ax)) + nsparse, reached = _sparsereach(len, bc.args...) + 32 * nsparse * reached <= len || return _densebroadcast(bc) + zbc = Broadcasted{Broadcast.DefaultArrayStyle{length(ax)}}(bc.f, map(_zerosparse, bc.args), ax) + # the zeros alone must give the element type that dense broadcast infers from all the values + T = Broadcast.combine_eltypes(zbc.f, zbc.args) + isconcretetype(T) && T === Broadcast.combine_eltypes(bc.f, bc.args) || return _densebroadcast(bc) + C = copy(Broadcast.instantiate(zbc)) + return _storedentries!(C, bc, bc.args...) +end +_zerosparse(bc::Broadcasted) = Broadcasted(bc.f, map(_zerosparse, bc.args)) +_zerosparse(S::SparseVecOrMat) = _densezero(S) +_zerosparse(x) = x +# the number of sparse arguments, and of the result entries their stored entries reach +_sparsereach(len) = (0, 0) +function _sparsereach(len, x, rest::Vararg{Any,N}) where {N} + n, r = _sparsereach(len, x) + n2, r2 = _sparsereach(len, rest...) + return (n + n2, r + r2) +end +_sparsereach(len, x::Broadcasted) = _sparsereach(len, x.args...) +_sparsereach(len, S::SparseVecOrMat) = (1, nnz(S) * (len รท length(S))) +_sparsereach(len, x) = (0, 0) +_sparseandscalars(bc::Broadcasted) = all(_sparseandscalars, bc.args) +_sparseandscalars(S::SparseVecOrMat) = eltype(S) <: Number +_sparseandscalars(x::Union{AbstractArray,Tuple}) = false +_sparseandscalars(x) = true + +# recompute the entries of `C` that the stored entries of the sparse arguments reach +_storedentries!(C, bc) = C +function _storedentries!(C, bc, x, rest::Vararg{Any,N}) where {N} + _storedentries!(C, bc, x) + return _storedentries!(C, bc, rest...) +end +_storedentries!(C, bc, x::Broadcasted) = _storedentries!(C, bc, x.args...) +_storedentries!(C, bc, x) = C +function _storedentries!(C, bc, S::AbstractCompressedVector) + for i in nonzeroinds(S) + _reevaluate!(C, bc, S, (Int(i),)) + end + return C +end +function _storedentries!(C, bc, S::AbstractSparseMatrixCSC) + rv = getrowval(S) + for j in axes(S, 2), k in nzrange(S, j) + _reevaluate!(C, bc, S, (Int(rv[k]), j)) + end + return C +end +# the stored entry at `I` of `S` reaches all of `C` along the dimensions where `S` has size 1, +# and index `I[d]` along each other dimension `d` +function _reevaluate!(C, bc, S, I::Tuple{Vararg{Int}}) + ranges = ntuple(d -> size(S, d) == 1 ? UnitRange(axes(C, d)) : (I[d]:I[d]), Val(ndims(C))) + for J in CartesianIndices(ranges) + C[J] = bc[J] + end + return C +end + function _densebroadcast(bc::Broadcasted) # an empty result reads no entry, so the sparse arguments are left as they are dbc = any(isempty, axes(bc)) ? bc : Broadcast.instantiate(_densifysparse(bc)) diff --git a/test/higherorderfns.jl b/test/higherorderfns.jl index 7ef4578d..c8ad4bdd 100644 --- a/test/higherorderfns.jl +++ b/test/higherorderfns.jl @@ -808,6 +808,10 @@ densesum(x, y) = x .+ y @test @inferred(broadcast(+, s, d))::Vector{Float64} == fs + d @test @inferred(broadcast(-, D, S))::Matrix{Float64} == D - fS @test @inferred(broadcast(+, s, 1))::Vector{Float64} == fs .+ 1 + # a sum of a sparse array with few stored entries and a scalar makes no dense copy of it + L = sparse([2, 400], [1, 3], [-0.0, NaN], 400, 3) + @test isequal(@inferred(densesum(L, 1))::Matrix{Float64}, Array(L) .+ 1) + @test @allocated(densesum(L, 1)) < 2 * sizeof(Array(L)) # the other broadcasts with a dense array or a scalar, and sums without one, stay sparse @test (S .* D)::SparseMatrixCSC == fS .* D @test (S .- S)::SparseMatrixCSC == fS .- fS @@ -830,6 +834,10 @@ densesum(x, y) = x .+ y @test @inferred(broadcast(-, 1, S))::Matrix{Float64} == 1 .- fS @test (s .+ fill(1.0))::Vector{Float64} == fs .+ 1 @test (S .+ s .- Ref(2))::Matrix{Float64} == fS .+ fs .- 2 + l = sparsevec([5], [2.0], 400) + @test isequal((L .+ l .- 2)::Matrix{Float64}, Array(L) .+ Array(l) .- 2) + R = SparseMatrixCSC{Real,Int}(400, 3, [1, 2, 2, 2], [7], Real[1.5]) + @test (R .+ 1)::Matrix{Real} == Array(R) .+ 1 @test (.-s)::SparseVector == .-fs # a sparse destination keeps the result sparse x = copy(s); x .+= d; x .-= 1 From 2759021cd2cca9006f9c1f2ac5b968daa08ed72a Mon Sep 17 00:00:00 2001 From: "Viral B. Shah" Date: Tue, 20 Oct 2026 21:03:20 +0200 Subject: [PATCH 02/30] Deprecate implementing `rowvals` and `nonzeros` on a subtype of `AbstractSparseMatrixCSC` (#1011) MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Stacked on #1008, which makes `getrowval`/`getnzval` the interface of `AbstractSparseMatrixCSC` and keeps a fallback to `rowvals`/`nonzeros` for subtypes that implement only those (the legacy interface). This PR deprecates the legacy interface. Its base is the #1008 branch, so it should merge after #1008. **Before** (#1008), under `--depwarn=yes`, a subtype `L` implementing only `rowvals`/`nonzeros` is used silently. **After**: ``` implementing `rowvals(::L)` as a storage accessor of a subtype of `AbstractSparseMatrixCSC` is deprecated, implement `SparseArrays.getrowval(::L)` instead. ``` The fallback methods of `getrowval`/`getnzval` call `Base.depwarn` with a message naming the type and the method to define, and the `AbstractSparseMatrixCSC` docstring says the fallback will be removed in SparseArrays 2.0. `Base.depwarn` walks the backtrace on every call and kernels reach the fallback once per entry, so it is skipped when warnings are off and made once per type and accessor under `--depwarn=yes`; under `--depwarn=error` every use throws. 80,000 scalar `getindex` calls on a 2000ร—2000 matrix: | `--depwarn` | legacy subtype | `SparseMatrixCSC` | |---|---|---| | `no` (default) | 0.77 ms, 0 allocations | 0.74 ms | | `yes` | 3.6 ms, 0 allocations | 0.74 ms | This should go in the release notes: packages with a legacy subtype see the warning when they test with depwarns on. Julia's own CI runs with `--depwarn=error`, where the test checks only that the accessors throw. Left out: `rowvals` and `nonzeros` as user-facing functions are not deprecated. ๐Ÿค– Generated with [Claude Code](https://claude.com/claude-code) Co-authored-by: Viral B. Shah Co-authored-by: Claude Opus 5.5 (cherry picked from commit 48204b263a4a1da16b5c6afed7efee1935b297e3) --- src/SparseArrays.jl | 10 +++++++--- src/sparsematrix.jl | 35 +++++++++++++++++++++++++++++++---- test/sparsematrix.jl | 26 ++++++++++++++++---------- 3 files changed, 54 insertions(+), 17 deletions(-) diff --git a/src/SparseArrays.jl b/src/SparseArrays.jl index 5c8ce1dd..04c8cf3b 100644 --- a/src/SparseArrays.jl +++ b/src/SparseArrays.jl @@ -122,9 +122,13 @@ are sorted, unique and in `1:size(S, 1)`. Everything else, such as [`nnz`](@ref) [`nzrange`](@ref), indexing, `copy`, `similar` and the sparse linear algebra, is derived from these. [`rowvals`](@ref) and [`nonzeros`](@ref) are not defined for such a subtype. -Subtypes written while `rowvals` and `nonzeros` were the only exported accessors implement -those in place of `getrowval` and `getnzval`, which fall back to them. This remains -supported, but is planned to stop being supported in SparseArrays 2.0. +!!! note "Legacy interface" + A subtype that implements `rowvals` and `nonzeros` in place of `getrowval` and + `getnzval` uses the legacy interface, which is deprecated. It keeps working, because + `getrowval` and `getnzval` fall back to `rowvals` and `nonzeros`, with a deprecation + warning when Julia runs with `--depwarn=yes` and an error with `--depwarn=error`. The + fallback will be removed in SparseArrays 2.0. To migrate, rename the two methods to + `SparseArrays.getrowval` and `SparseArrays.getnzval`. """ abstract type AbstractSparseMatrixCSC{Tv,Ti<:Integer} <: AbstractSparseMatrix{Tv,Ti} end diff --git a/src/sparsematrix.jl b/src/sparsematrix.jl index 584fda58..9fc496df 100644 --- a/src/sparsematrix.jl +++ b/src/sparsematrix.jl @@ -187,6 +187,29 @@ const _SparseTriOrSymHerm{Tv,Ti} = Union{UpperTriangular{Tv,<:SparseMatrixCSCOrV SparseMatrixCSCSymmHerm{Tv,Ti}} _storageparent(S) = S _storageparent(S::_SparseTriOrSymHerm) = parent(S) + +# The legacy interface of `AbstractSparseMatrixCSC`: a subtype that implements `rowvals` and +# `nonzeros` in place of `getrowval` and `getnzval`. Kernels reach these fallbacks once per +# entry and `Base.depwarn` walks the backtrace on every call, so it is skipped when warnings +# are off and made once per type and accessor when they are on. +const _legacy_accessor_warned = Set{Tuple{Core.TypeName,Symbol}}() +const _legacy_accessor_lock = ReentrantLock() +@inline _depwarn_legacy_accessor(S, get::Symbol, legacy::Symbol) = + Base.JLOptions().depwarn == 0 ? nothing : _depwarn_legacy_accessor_on(S, get, legacy) +@noinline function _depwarn_legacy_accessor_on(S, get::Symbol, legacy::Symbol) + if Base.JLOptions().depwarn == 1 + key = (Base.typename(typeof(S)), get) + seen = @lock _legacy_accessor_lock (key in _legacy_accessor_warned || + (push!(_legacy_accessor_warned, key); false)) + seen && return nothing + end + T = nameof(typeof(S)) + Base.depwarn("implementing `$legacy(::$T)` as a storage accessor of a subtype of " * + "`AbstractSparseMatrixCSC` is deprecated, implement `SparseArrays.$get(::$T)` instead.", + get) + return nothing +end + """ getrowval(A) @@ -203,6 +226,9 @@ the parent's vector, of which `rowvals(A)` is the part belonging to `A`. It is n defined for a triangular, `Symmetric` or `Hermitian` wrapper of a sparse matrix; use `getrowval(parent(A))`, which also holds the entries outside the wrapper's triangle. +A subtype of [`AbstractSparseMatrixCSC`](@ref) implements `getrowval`. Calling it on a subtype +that implements only `rowvals` is deprecated. + # Examples ```jldoctest julia> A = sparse(2I, 3, 3) @@ -225,8 +251,7 @@ julia> getrowval(sparsevec([2, 5], [3.0, 4.0])) """ getrowval(S::SparseMatrixCSC) = getfield(S, :rowval) getrowval(S::FixedSparseCSC) = getfield(S, :rowval) -# a subtype may implement the longer-exported `rowvals` instead -getrowval(S::AbstractSparseMatrixCSC) = rowvals(S) +getrowval(S::AbstractSparseMatrixCSC) = (_depwarn_legacy_accessor(S, :getrowval, :rowvals); rowvals(S)) getrowval(S::SparseMatrixCSCColumnSubset) = getrowval(parent(S)) """ @@ -244,6 +269,9 @@ the parent's vector, of which `nonzeros(A)` is the part belonging to `A`. It is defined for a triangular, `Symmetric` or `Hermitian` wrapper of a sparse matrix; use `getnzval(parent(A))`, which also holds the entries outside the wrapper's triangle. +A subtype of [`AbstractSparseMatrixCSC`](@ref) implements `getnzval`. Calling it on a subtype +that implements only `nonzeros` is deprecated. + # Examples ```jldoctest julia> A = sparse(2I, 3, 3) @@ -266,8 +294,7 @@ julia> getnzval(sparsevec([2, 5], [3.0, 4.0])) """ getnzval(S::SparseMatrixCSC) = getfield(S, :nzval) getnzval(S::FixedSparseCSC) = getfield(S, :nzval) -# a subtype may implement the longer-exported `nonzeros` instead -getnzval(S::AbstractSparseMatrixCSC) = nonzeros(S) +getnzval(S::AbstractSparseMatrixCSC) = (_depwarn_legacy_accessor(S, :getnzval, :nonzeros); nonzeros(S)) getnzval(S::SparseMatrixCSCColumnSubset) = getnzval(parent(S)) nzvalview(S::AbstractSparseMatrixCSC) = view(getnzval(S), 1:nnz(S)) nzvalview(S::SparseMatrixCSCColumnSubset) = nonzeros(S) diff --git a/test/sparsematrix.jl b/test/sparsematrix.jl index 3e6b169a..51bad629 100644 --- a/test/sparsematrix.jl +++ b/test/sparsematrix.jl @@ -208,18 +208,24 @@ end @test isone(SparseMatrixCSC(2, 2, [1, 3, 4], [1, 2, 2], [1, 0, 1])) # stored zero off-diagonal is fine end -@testset "a subtype implementing rowvals and nonzeros but not the get* accessors" begin +@testset "the legacy interface: a subtype implementing rowvals and nonzeros but not the get* accessors" begin S = sparse([1.0 0 2; 0 3 0; 4 0 5]) A = LegacyCSC(S) - @test SparseArrays.getrowval(A) === rowvals(S) && SparseArrays.getnzval(A) === nonzeros(S) - @test nonzeros(view(A, :, 2:3)) == nonzeros(view(S, :, 2:3)) - @test A[2, 2] == 3 - @test copy(A) == S - x = [1.0, 2, 3] - @test A * x == S * x - @test A' * x == S' * x - @test A * S == S * A == A * A == S * S - @test A * Matrix(S) == S * Matrix(S) + # the warning is given once per type and accessor, so these come before any other use + @test_deprecated r"getrowval" SparseArrays.getrowval(A) + @test_deprecated r"getnzval" SparseArrays.getnzval(A) + # under --depwarn=error, which Julia's own CI uses, every use of the type throws + if Base.JLOptions().depwarn != 2 + @test SparseArrays.getrowval(A) === rowvals(S) && SparseArrays.getnzval(A) === nonzeros(S) + @test nonzeros(view(A, :, 2:3)) == nonzeros(view(S, :, 2:3)) + @test A[2, 2] == 3 + @test copy(A) == S + x = [1.0, 2, 3] + @test A * x == S * x + @test A' * x == S' * x + @test A * S == S * A == A * A == S * S + @test A * Matrix(S) == S * Matrix(S) + end end @testset "indtype" begin From 29be2059ad8deef61b10d8919862e5cd3a97f443 Mon Sep 17 00:00:00 2001 From: "Viral B. Shah" Date: Mon, 26 Oct 2026 13:26:16 +0000 Subject: [PATCH 03/30] Convert the values of a sparse matrix before `copyto!` writes them into a sparse vector (#1107) MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit `copyto!(::AbstractCompressedVector, ::AbstractSparseMatrixCSC)` resized the destination and wrote the matrix's row indices into it before converting the values, so a value that does not convert to the vector's eltype left the vector with the matrix's pattern and uninitialised values. **Before** ```julia julia> a = sparsevec([1], [1], 4); julia> copyto!(a, sparse([1, 2], [1, 2], [0.5, 1.0])) ERROR: InexactError: Int64(0.5) julia> (SparseArrays.nonzeroinds(a), nonzeros(a)) ([1, 4], [1, 123754957609376]) ``` **After** ```julia julia> copyto!(a, sparse([1, 2], [1, 2], [0.5, 1.0])) ERROR: InexactError: Int64(0.5) julia> (SparseArrays.nonzeroinds(a), nonzeros(a)) ([1], [1]) ``` The method now runs `_checkconvertible` on the matrix's values before `prep_sparsevec_copy_dest!` touches the destination, as the matrix-to-matrix `copyto!` does, and the `BoundsError` check on the source's column pointers moves ahead of it too. The other preconditions already held before the first write: the length check throws in `prep_sparsevec_copy_dest!` before it resizes, every index fits the vector's index type because its length does, and a fixed-pattern destination goes through `_copyto_fixed!` on a converted copy. ๐Ÿค– Generated with [Claude Code](https://claude.com/claude-code) Co-authored-by: Claude Opus 5.5 (cherry picked from commit b6282ef9cf55eb2f05f743fbcf534f3d58acaaa3) --- src/sparsevector.jl | 3 ++- test/sparsevector.jl | 2 ++ 2 files changed, 4 insertions(+), 1 deletion(-) diff --git a/src/sparsevector.jl b/src/sparsevector.jl index 9c10b931..5d8f997a 100644 --- a/src/sparsevector.jl +++ b/src/sparsevector.jl @@ -669,11 +669,12 @@ end function copyto!(A::AbstractCompressedVector, B::AbstractSparseMatrixCSC) _is_fixed(A) && return _copyto_fixed!(A, copyto!(spzeros(eltype(A), indtype(A), length(B)), B)) + maximum(getcolptr(B))-1 <= length(getrowval(B)) || throw(BoundsError()) + _checkconvertible(eltype(A), getnzval(B)) prep_sparsevec_copy_dest!(A, length(B), nnz(B)) ptr = 1 @assert length(nonzeroinds(A)) >= length(getrowval(B)) - maximum(getcolptr(B))-1 <= length(getrowval(B)) || throw(BoundsError()) @inbounds for col=1:length(getcolptr(B))-1 offsetA = (col - 1) * size(B, 1) while ptr <= getcolptr(B)[col+1]-1 diff --git a/test/sparsevector.jl b/test/sparsevector.jl index 07312934..303195a4 100644 --- a/test/sparsevector.jl +++ b/test/sparsevector.jl @@ -679,6 +679,8 @@ end @test mismatch(dst, [1, 0, 0, 0, 0]) === nothing && nnz(dst) == 1 @test_throws InexactError copyto!(dst, src) @test mismatch(dst, [1, 0, 0, 0, 0]) === nothing && nnz(dst) == 1 + @test_throws InexactError copyto!(dst, sparse([1, 2], [1, 2], [0.5, 1.0])) + @test mismatch(dst, [1, 0, 0, 0, 0]) === nothing && nnz(dst) == 1 @test_throws InexactError copy!(dst, [1.0, 2.5, 0.0, 0.0, 0.0]) @test mismatch(dst, [1, 0, 0, 0, 0]) === nothing && nnz(dst) == 1 @test mismatch(copy!(dst, sparsevec([2, 4], [3.0, 4.0], 5)), [0, 3, 0, 4, 0]) === nothing From 74cbf158c0c9af776290304c84cbf0ed9ea816a6 Mon Sep 17 00:00:00 2001 From: "Viral B. Shah" Date: Mon, 26 Oct 2026 13:26:50 +0000 Subject: [PATCH 04/30] Check that a `CHOLMOD.Sparse` has one column in `sparsevec` (#1108) MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit `sparsevec(::CHOLMOD.Sparse)` never checks the column count, so a matrix with more than one column comes back as a `SparseVector` whose indices repeat, and a symmetric `stype` fails a user-reachable `@assert`. **Before** ```julia julia> using SparseArrays, LinearAlgebra; using SparseArrays: CHOLMOD, nonzeroinds julia> v = sparsevec(CHOLMOD.Sparse(sparse([1.0 2; 0 3]))); (length(v), nonzeroinds(v)) (2, [1, 1, 2]) julia> sparsevec(CHOLMOD.Sparse(Symmetric(sparse([2.0 1; 1 2])))) ERROR: AssertionError: s.stype == 0 ``` **After** ```julia julia> sparsevec(CHOLMOD.Sparse(sparse([1.0 2; 0 3]))) ERROR: DimensionMismatch: cannot convert a CHOLMOD.Sparse with 2 columns to SparseVector: it must have exactly one column julia> sparsevec(CHOLMOD.Sparse(Symmetric(sparse([2.0 1; 1 2])))) ERROR: ArgumentError: matrix has stype != 0. Convert to matrix with stype == 0 before converting to SparseVector ``` `SparseVector{Tv,Ti}(::CHOLMOD.Sparse)` checked the `stype` but then copied the matrix's row indices and values into a vector of length `nrow`, whatever `ncol` was. It now throws `DimensionMismatch` when `ncol != 1`, before copying anything. `sparsevec` had its own `@assert` on the `stype` ahead of that constructor; it now just calls the constructor, which already throws an `ArgumentError` with the reason. There is no other entry point: `SparseVector(::CHOLMOD.Sparse)` and `convert(SparseVector, ...)` have no methods. Left out: `sparsevec(::SparseMatrixCSC)` vectorizes a matrix with several columns (`A[:]`). Doing the same for a `CHOLMOD.Sparse` would add behaviour, and throwing leaves that option open. ๐Ÿค– Generated with [Claude Code](https://claude.com/claude-code) Co-authored-by: Claude Opus 5.5 (cherry picked from commit 71e5eed08cac7c976832d628e6f7bf4158f4e5de) --- src/solvers/cholmod.jl | 10 +++++----- test/solvers/issues.jl | 5 +++++ 2 files changed, 10 insertions(+), 5 deletions(-) diff --git a/src/solvers/cholmod.jl b/src/solvers/cholmod.jl index e0d03046..977da253 100644 --- a/src/solvers/cholmod.jl +++ b/src/solvers/cholmod.jl @@ -1333,6 +1333,10 @@ function SparseVector{Tv, Ti}(A::Sparse{Tv, Ti}) where {Tv, Ti<:ITypes} throw(ArgumentError("matrix has stype != 0. Convert to matrix " * "with stype == 0 before converting to SparseVector")) end + if s.ncol != 1 + throw(DimensionMismatch("cannot convert a CHOLMOD.Sparse with $(Int(s.ncol)) columns " * + "to SparseVector: it must have exactly one column")) + end args = _extract_args(A, Tv) s.sorted == 0 && _sort_buffers!(args...); _trim_nz_builder!(args...) @@ -1384,11 +1388,7 @@ function Hermitian{Tv,SparseMatrixCSC{Tv, Ti}}(A::Sparse{Tv, Ti}) where {Tv<:VTy end convert(T::Type{Hermitian{Tv,SparseMatrixCSC{Tv,Ti}}}, A::Sparse{Tv, Ti}) where {Tv<:VTypes, Ti<:ITypes} = T(A) -function sparsevec(A::Sparse{Tv, Ti}) where {Tv, Ti} - s = unsafe_load(pointer(A)) - @assert s.stype == 0 - return SparseVector{Tv, Ti}(A) -end +sparsevec(A::Sparse{Tv, Ti}) where {Tv, Ti} = SparseVector{Tv, Ti}(A) function sparse(A::Sparse{Tv, Ti}) where {Tv<:VRealTypes, Ti} # Notice! Cannot be type stable because of stype s = unsafe_load(pointer(A)) diff --git a/test/solvers/issues.jl b/test/solvers/issues.jl index 13bb94e5..fd6b6879 100644 --- a/test/solvers/issues.jl +++ b/test/solvers/issues.jl @@ -212,4 +212,9 @@ end @test eltype(F.Q) == eltype(F.R) == eltyA end +@testset "sparsevec of a CHOLMOD.Sparse checks its shape and stype" begin + @test_throws DimensionMismatch sparsevec(CHOLMOD.Sparse(sparse([1.0 2; 0 3]))) + @test_throws ArgumentError sparsevec(CHOLMOD.Sparse(Symmetric(sparse([2.0 1; 1 2])))) +end + end # module From f381159226a5bd32a0c07ba39f67bbeb1739d271 Mon Sep 17 00:00:00 2001 From: "Viral B. Shah" Date: Mon, 26 Oct 2026 13:27:45 +0000 Subject: [PATCH 05/30] Check that the result of sparse `vcat`, `hcat` and `repeat` fits the index type (#1110) MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit `vcat` and `repeat` of sparse matrices write out of bounds when the result has more stored entries than the index type can count. With bounds checks off, as by default, the writes land at negative indices under `@inbounds` and corrupt memory, which shows up as a segfault in a later GC. **Before** ```julia julia> A = SparseMatrixCSC{Bool,Int16}(sparse(trues(130, 130))); # 16900 stored entries julia> vcat(A, A) # with --check-bounds=yes ERROR: BoundsError: attempt to access 33800-element Vector{Int16} at index [-32645] julia> repeat(A, 2) # with --check-bounds=yes ERROR: BoundsError: attempt to access 33800-element Vector{Int16} at index [-32645] julia> vcat(A, A) # default bounds checking, after writing out of bounds ERROR: ArgumentError: 32761 == colptr[127] > colptr[128] == -32515 julia> repeat(A, 2) ERROR: AssertionError: colptr[var"end"] == nnz_new + 1 julia> hcat(A, A) # and likewise repeat(A, 1, 2) ERROR: InexactError: trunc(Int16, 32891) julia> v = SparseVector{Bool,Int16}(sparsevec(trues(20000))); hcat(v, v) ERROR: InexactError: trunc(Int16, 40001) ``` The same happens with `Int32` indices once a result exceeds 2^31 stored entries: `vcat` of two `SparseMatrixCSC{Bool,Int32}` with 1.08e9 entries each writes at index -2147483647. **After** ```julia julia> vcat(A, A) # and likewise repeat(A, 2), repeat(A, 1, 2), hcat(A, A), hcat(v, v) ERROR: ArgumentError: the index type Int16 cannot hold 33800 stored entries; use a larger index type ``` `_vcat_csc` and `repeat(A, m)` carry the running column pointer in the index type, from `colptr[c]`, so it wraps past `typemax(Ti)`, and `stuffcol!` then writes `rowval[k]` and `nzval[k]` at the wrapped index. Each now checks the total with `_checknnzfits` before allocating, as `blockdiag` and `hvcat` already do. `hcat` of sparse matrices, `repeat(A, m, n)` and `hcat`/`stack` of sparse vectors could not corrupt memory, since they compute in `Int` and fail on the conversion, but they get the same check so that every concatenation gives the same error. The `@assert` in `repeat`, which this made unreachable, is gone. Left out: `vcat` of sparse vectors still throws an `InexactError` when the result's length does not fit the index type, which happens before its stored entries can overflow; that is a length check, not this one. ๐Ÿค– Generated with [Claude Code](https://claude.com/claude-code) Co-authored-by: Claude Opus 5.5 (cherry picked from commit 8f7bd49dd9cb2164cb1396a1cbb5bc746445fc30) --- src/concatenation.jl | 6 +++++- test/concatenation.jl | 4 ++++ 2 files changed, 9 insertions(+), 1 deletion(-) diff --git a/src/concatenation.jl b/src/concatenation.jl index 54c0b409..8646554b 100644 --- a/src/concatenation.jl +++ b/src/concatenation.jl @@ -24,6 +24,7 @@ function _vcat_csc(::Type{Tv}, ::Type{Ti}, X::AbstractSparseMatrixCSC...) where nnzX = Int[ nnz(x) for x in X ] nnz_res = sum(nnzX) + _checknnzfits(Ti, nnz_res) colptr = Vector{Ti}(undef, n+1) rowval = Vector{Ti}(undef, nnz_res) nzval = Vector{Tv}(undef, nnz_res) @@ -72,6 +73,7 @@ function hcat(X::AbstractSparseMatrixCSC...) colptr = Vector{Ti}(undef, n+1) nnzX = Int[ nnz(x) for x in X ] nnz_res = sum(nnzX) + _checknnzfits(Ti, nnz_res) rowval = Vector{Ti}(undef, nnz_res) nzval = Vector{Tv}(undef, nnz_res) @@ -99,6 +101,7 @@ end function Base.repeat(A::AbstractSparseMatrixCSC, m) nnz_new = nnz(A) * m + _checknnzfits(indtype(A), nnz_new) colptr = similar(getcolptr(A), length(getcolptr(A))) rowval = similar(getrowval(A), nnz_new) nzval = similar(getnzval(A), nnz_new) @@ -115,12 +118,12 @@ function Base.repeat(A::AbstractSparseMatrixCSC, m) end colptr[c + 1] = ptr_res end - @assert colptr[end] == nnz_new + 1 SparseMatrixCSC(size(A, 1) * m, size(A, 2), colptr, rowval, nzval) end function Base.repeat(A::AbstractSparseMatrixCSC, m, n) + _checknnzfits(indtype(A), nnz(A) * m * n) B = repeat(A, m) nB = size(B, 2) nnzB = nnz(B) @@ -221,6 +224,7 @@ function Base._typed_stack(::Colon, ::Type{Tv}, ::Type{S}, A, Aax::Tuple{Any}) w throw(DimensionMismatch("Inconsistent column lengths.")) tnnz += nnz(x) end + _checknnzfits(Ti, tnnz) colptr = Vector{Ti}(undef, n+1) nzrow = Vector{Ti}(undef, tnnz) diff --git a/test/concatenation.jl b/test/concatenation.jl index 726cd3fb..ef1cf303 100644 --- a/test/concatenation.jl +++ b/test/concatenation.jl @@ -36,6 +36,10 @@ include("testhelpers.jl") se33_32bit = convert(SparseMatrixCSC{Float32,Int32}, se33) @test [se33; se33_32bit] == [Array(se33); Array(se33_32bit)] @test length(nonzeros([sp33; 0I])) == 3 + # the 128 stored entries of the result do not fit `Int8` + A8 = SparseMatrixCSC{Float64,Int8}(sparse(ones(8, 8))) + @test_throws ArgumentError vcat(A8, A8) + @test_throws ArgumentError repeat(A8, 2) end end From 50198881315d43cf8cf2b06610c1a1159e2c277e Mon Sep 17 00:00:00 2001 From: "Viral B. Shah" Date: Mon, 26 Oct 2026 13:43:54 +0000 Subject: [PATCH 06/30] Check the workspace lengths in `sparse!` before writing (#1111) MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit The expert driver `SparseArrays.sparse!(I, J, V, m, n, combine, klasttouch, csrrowptr, csrcolval, csrnzval, ...)` documents that its workspace needs `length(csrrowptr) >= m + 1`, `length(csrcolval) >= length(I)`, `length(csrnzval) >= length(I)` and `length(klasttouch) >= n`, but never checked, and writes all four under `@inbounds`. A short workspace corrupts the heap or segfaults. **Before** ```julia julia> SparseArrays.sparse!([1, 2], [1, 2], [1.0, 2.0], 2, 2, +, Int[], [0], Int[], Float64[], Int[], Int[], Float64[]) [2699137] signal 11 (128): Segmentation fault sparse! at src/constructors.jl:506 ``` **After** ```julia julia> SparseArrays.sparse!([1, 2], [1, 2], [1.0, 2.0], 2, 2, +, Int[], [0], Int[], Float64[], Int[], Int[], Float64[]) ERROR: ArgumentError: 0 == length(klasttouch) < 2 ``` The four lengths are now checked in O(1) with the existing `sparse_check_length`, next to the checks on `J` and `V` and before `csrrowptr` is first written, so a throw leaves the workspace untouched. `csrnzval` is not checked when `combine === nothing`, the `spzeros!` path, where it aliases `cscnzval` and is never read. The `spzeros!` methods and the shorter `sparse!` methods forward here, so they are covered too. The docstring now says the lengths are checked, and fixes the misplaced parenthesis in `length(csrnzval) >= length(I)`. ๐Ÿค– Generated with [Claude Code](https://claude.com/claude-code) --------- Co-authored-by: Claude Opus 5.5 (cherry picked from commit 353601e5467a9d17b5cb624a1e60f9f8207b6f05) --- src/constructors.jl | 15 ++++++++++----- test/constructors.jl | 8 ++++++++ 2 files changed, 18 insertions(+), 5 deletions(-) diff --git a/src/constructors.jl b/src/constructors.jl index ca8db3dd..fbfc69e5 100644 --- a/src/constructors.jl +++ b/src/constructors.jl @@ -429,7 +429,7 @@ preceding CSR form into a fully-sorted CSC form with no repeated entries. Input arrays `csrrowptr`, `csrcolval`, and `csrnzval` constitute storage for the intermediate CSR forms and require `length(csrrowptr) >= m + 1`, -`length(csrcolval) >= length(I)`, and `length(csrnzval >= length(I))`. Input +`length(csrcolval) >= length(I)`, and `length(csrnzval) >= length(I)`. Input array `klasttouch`, workspace for the second stage, requires `length(klasttouch) >= n`. Optional input arrays `csccolptr`, `cscrowval`, and `cscnzval` constitute storage for the returned CSC form `S`. If necessary, these are resized automatically to satisfy @@ -447,8 +447,8 @@ You may reuse the input arrays' storage (`I`, `J`, `V`) for the output arrays Note that they will be resized to satisfy the conditions above. For the sake of efficiency, this method performs no argument checking beyond -`1 <= I[k] <= m` and `1 <= J[k] <= n`. Use with care. Testing with `--check-bounds=yes` -is wise. +`1 <= I[k] <= m`, `1 <= J[k] <= n` and the lengths of the input arrays. Use with care. +Testing with `--check-bounds=yes` is wise. This method runs in `O(m, n, length(I))` time. The HALFPERM algorithm described in F. Gustavson, "Two fast algorithms for sparse matrices: multiplication and permuted @@ -470,11 +470,16 @@ function sparse!(I::AbstractVector{Ti}, J::AbstractVector{Ti}, V::AbstractVector # alias cscnzval, which will be resized and filled with zero(Tv). only_sparsity_pattern = combine === nothing - # Compute the CSR form's row counts and store them shifted forward by one in csrrowptr - fill!(csrrowptr, Tj(0)) coolen = length(I) length(J) >= coolen || throw(ArgumentError("J need length >= length(I) = $coolen")) only_sparsity_pattern || length(V) >= coolen || throw(ArgumentError("V need length >= length(I) = $coolen")) + sparse_check_length("klasttouch", klasttouch, n, String) + sparse_check_length("csrrowptr", csrrowptr, m + 1, String) + sparse_check_length("csrcolval", csrcolval, coolen, String) + only_sparsity_pattern || sparse_check_length("csrnzval", csrnzval, coolen, String) + + # Compute the CSR form's row counts and store them shifted forward by one in csrrowptr + fill!(csrrowptr, Tj(0)) @inbounds for k in 1:coolen Ik = I[k] diff --git a/test/constructors.jl b/test/constructors.jl index 91be8dad..1cc66d3d 100644 --- a/test/constructors.jl +++ b/test/constructors.jl @@ -494,6 +494,14 @@ end end end end + + @static if COMPREHENSIVE + # a short workspace throws before anything is written + csrrowptr = [7, 7] + @test_throws ArgumentError sparse!([1, 2], [1, 2], [1.0, 2.0], 2, 2, +, + [0, 0], csrrowptr, [0, 0], [0.0, 0.0]) + @test csrrowptr == [7, 7] + end end end # module From 7839ead9475d71113600af615e7b22ecfebc62f5 Mon Sep 17 00:00:00 2001 From: "Viral B. Shah" Date: Mon, 26 Oct 2026 13:57:18 +0000 Subject: [PATCH 07/30] Check that the values convert before `lu!(F::UmfpackLU, S)` writes to `F` (#1114) MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit If converting the values of `S` throws, `lu!(F::UmfpackLU, S)` has already written the new size and pattern to `F`. The old numeric factor stays in place with a successful status, so `F` reports the new size, and a solve with it reads past the end of what UMFPACK factored. **Before** ```julia julia> A = sparse(ComplexF64[4 1 0 0; 1 4 1 0; 0 1 4 1; 0 0 1 4]); F = lu(A); julia> S = SparseMatrixCSC{Any,Int}(2, 2, [1, 2, 4], [1, 1, 2], Any[1.0, "x", 1.0]); julia> lu!(F, S) ERROR: MethodError: Cannot `convert` an object of type String to an object of type ComplexF64 julia> size(F) (2, 2) julia> F \ ones(ComplexF64, 2) 2-element Vector{ComplexF64}: 0.19617224880382775 + 9.683892736703e-312im 0.21531100478468898 - 3.87355709468e-311im julia> F \ ones(ComplexF64, 4) ERROR: DimensionMismatch: UmfpackLU is 2ร—2 but the right-hand side has 4 rows ``` **After** ```julia julia> lu!(F, S) ERROR: MethodError: Cannot `convert` an object of type String to an object of type ComplexF64 julia> size(F) (4, 4) julia> F \ ones(ComplexF64, 2) ERROR: DimensionMismatch: UmfpackLU is 4ร—4 but the right-hand side has 2 rows julia> F \ ones(ComplexF64, 4) โ‰ˆ Matrix(A) \ ones(ComplexF64, 4) true ``` `lu!` already checks the index range and the real/complex match before it takes the lock. The value conversion was the only step that could still throw after the first write, in `F.nzval .= getnzval(S)`. It now runs `_checkconvertible(Tv, getnzval(S))` along with the other checks, so a matrix whose values do not convert throws the same error with `F` unchanged. For an eltype that is already a subtype of `Tv`, the check compiles to nothing. ๐Ÿค– Generated with [Claude Code](https://claude.com/claude-code) Co-authored-by: Claude Opus 5.5 (cherry picked from commit a5211ce67b269d1b4952dc9825219fe6afd61835) --- src/solvers/umfpack.jl | 3 ++- test/solvers/issues.jl | 12 ++++++++++++ 2 files changed, 14 insertions(+), 1 deletion(-) diff --git a/src/solvers/umfpack.jl b/src/solvers/umfpack.jl index 68ffb597..f9e8e6db 100644 --- a/src/solvers/umfpack.jl +++ b/src/solvers/umfpack.jl @@ -13,7 +13,7 @@ import LinearAlgebra: Factorization, AdjointFactorization, TransposeFactorizatio checksquare, det, logabsdet, lu, lu!, ldiv! using SparseArrays -using SparseArrays: getcolptr, getrowval, getnzval, AbstractSparseMatrixCSC +using SparseArrays: getcolptr, getrowval, getnzval, AbstractSparseMatrixCSC, _checkconvertible import SparseArrays: nnz import Serialization: AbstractSerializer, deserialize, serialize @@ -531,6 +531,7 @@ function lu!(F::UmfpackLU{Tv, Ti}, S::AbstractSparseMatrixCSC; if Tv <: Real && !(eltype(S) <: Real) throw(ArgumentError("cannot refactorize the real $(typeof(F)) with a matrix of eltype $(eltype(S)); use lu(S) instead")) end + _checkconvertible(Tv, getnzval(S)) @lock F.lock begin # the symbolic factors are those of a matrix of another size diff --git a/test/solvers/issues.jl b/test/solvers/issues.jl index fd6b6879..e48e603f 100644 --- a/test/solvers/issues.jl +++ b/test/solvers/issues.jl @@ -195,6 +195,18 @@ end @test all(logabsdet(F) .โ‰ˆ logabsdet(Matrix(A))) end +@testset "lu!(F, S) leaves F unchanged when the values of S do not convert" begin + A = sparse(ComplexF64[4 1 0 0; 1 4 1 0; 0 1 4 1; 0 0 1 4]) + F = lu(A) + colptr, rowval, nzval = copy(F.colptr), copy(F.rowval), copy(F.nzval) + S = SparseMatrixCSC{Any,Int}(2, 2, [1, 2, 4], [1, 1, 2], Any[1.0, "x", 1.0]) + @test_throws MethodError lu!(F, S) + @test size(F) == (4, 4) + @test (F.colptr, F.rowval, F.nzval) == (colptr, rowval, nzval) + b = ComplexF64[1, 2, 3, 4] + @test F \ b โ‰ˆ Matrix(A) \ b +end + @testset "Issues #26367 and #26368" begin A = sparse([0.0 1 0 0; 0 0 0 0]) F = qr(A) From 2cc94fa80c4353d0d5a6908008276d2e9fddcb89 Mon Sep 17 00:00:00 2001 From: "Viral B. Shah" Date: Mon, 26 Oct 2026 13:57:48 +0000 Subject: [PATCH 08/30] Solve a real `UmfpackLU` against a `ComplexF32` vector in `ldiv!` (#1112) MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit `ldiv!` of a real `UmfpackLU` with a `ComplexF32` vector throws a `MethodError`; the matrix case and `\` work. **Before** ```julia julia> F = lu(sparse([2.0 1; 1 3])); ldiv!(F, ComplexF32[1, 2]) ERROR: MethodError: no method matching solve!(::Vector{Float64}, ::SparseArrays.UMFPACK.UmfpackLU{Float64, Int64}, ::Vector{Float32}, ::Int64; workspace::SparseArrays.UMFPACK.UmfpackWS{Int64}) ``` **After** ```julia julia> F = lu(sparse([2.0 1; 1 3])); ldiv!(F, ComplexF32[1, 2]) 2-element Vector{ComplexF32}: 0.2f0 + 0.0f0im 0.6f0 + 0.0f0im ``` The vector kernel for a complex right-hand side and a real factor solves the real and imaginary parts separately, through a scratch buffer it built with `c = real.(b)`. That buffer has the real eltype of `b`, `Float32` here, and UMFPACK only solves in `Float64`. It is now allocated as `Float64` and filled in place, as the matrix kernel already does. The 3-argument form and `transpose(F)` and `F'` go through the same kernel. ๐Ÿค– Generated with [Claude Code](https://claude.com/claude-code) Co-authored-by: Claude Opus 5.5 (cherry picked from commit 8036bf350886574b7b1131d12c31f7a99d80f1c8) --- src/solvers/umfpack.jl | 3 ++- test/solvers/issues.jl | 5 +++++ 2 files changed, 7 insertions(+), 1 deletion(-) diff --git a/src/solvers/umfpack.jl b/src/solvers/umfpack.jl index f9e8e6db..ed9fce7f 100644 --- a/src/solvers/umfpack.jl +++ b/src/solvers/umfpack.jl @@ -1164,7 +1164,8 @@ function _AqldivB_kernel!(x::StridedVector{Tb}, lu::UmfpackLU{Float64}, b::StridedVector{Tb}, transposeoptype, workspace) where Tb<:Complex r = similar(b, Float64) i = similar(b, Float64) - c = real.(b) + c = similar(b, Float64) + c .= real.(b) solve!(r, lu, c, transposeoptype; workspace) c .= imag.(b) solve!(i, lu, c, transposeoptype; workspace) diff --git a/test/solvers/issues.jl b/test/solvers/issues.jl index e48e603f..3674c2f9 100644 --- a/test/solvers/issues.jl +++ b/test/solvers/issues.jl @@ -229,4 +229,9 @@ end @test_throws ArgumentError sparsevec(CHOLMOD.Sparse(Symmetric(sparse([2.0 1; 1 2])))) end +@testset "ldiv! of a real UmfpackLU with a ComplexF32 vector" begin + A = sparse([2.0 1; 1 3]) + @test ldiv!(lu(A), ComplexF32[1, 2]) โ‰ˆ ComplexF32.(Matrix(A) \ [1, 2]) +end + end # module From 4448d27cbc7efb13ee55a5fdecad77eabf15e3e1 Mon Sep 17 00:00:00 2001 From: "Viral B. Shah" Date: Mon, 26 Oct 2026 13:58:55 +0000 Subject: [PATCH 09/30] Check that `circshift!` of a sparse matrix can hold the result before writing (#1115) MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit `circshift!(O, X, (r, c))` with a nonzero column shift resized `O`'s buffers and rewrote its pattern before converting `X`'s values to `eltype(O)` or checking that `indtype(O)` counts `nnz(X)` entries. A failed conversion left `O` changed, and a too-narrow index type left it corrupt. **Before** ```julia julia> O = sparse([1 0 0; 0 2 0; 0 0 3]); julia> circshift!(O, sparse([1.5 0 0; 0 0 0; 0 0 0]), (0, 1)) ERROR: InexactError: Int64(1.5) julia> Matrix(O) 3ร—3 Matrix{Int64}: 0 1 0 0 0 0 0 0 0 julia> O8 = SparseMatrixCSC{Float64,Int8}(sparse(1.0I, 12, 12)); julia> circshift!(O8, sparse(ones(12, 12)), (0, 1)) ERROR: InexactError: trunc(Int8, 145) julia> Matrix(O8) ERROR: AssertionError: _goodbuffers(S) ``` **After** ```julia julia> circshift!(O, sparse([1.5 0 0; 0 0 0; 0 0 0]), (0, 1)) ERROR: InexactError: Int64(1.5) julia> Matrix(O) 3ร—3 Matrix{Int64}: 1 0 0 0 2 0 0 0 3 julia> circshift!(O8, sparse(ones(12, 12)), (0, 1)) ERROR: ArgumentError: the index type Int8 cannot hold 144 stored entries; use a larger index type julia> O8 == sparse(1.0I, 12, 12) true ``` The column-shift branch now calls `_checkconvertible(eltype(O), view(getnzval(X), 1:nz))` and `_checknnzfits(indtype(O), nz)` before its first write, as `copyto!` does. The other paths were already safe: a zero column shift and an empty source go through `copy!`, a fixed-pattern destination through `_copyto_fixed!`, and the sparse-vector `circshift!` through `copy!`, which validates first. ๐Ÿค– Generated with [Claude Code](https://claude.com/claude-code) Co-authored-by: Claude Opus 5.5 (cherry picked from commit dc082aba6477430a30f67cb7cf87e7599c3e8455) --- src/sparsematrix.jl | 2 ++ test/linalg.jl | 7 +++++++ 2 files changed, 9 insertions(+) diff --git a/src/sparsematrix.jl b/src/sparsematrix.jl index 9fc496df..102db346 100644 --- a/src/sparsematrix.jl +++ b/src/sparsematrix.jl @@ -2292,6 +2292,8 @@ function circshift!(O::AbstractSparseMatrixCSC, X::AbstractSparseMatrixCSC, (r,c Base.mightalias(getnzval(O), getnzval(X))) && throw(ArgumentError(LazyString("cannot shift the columns of a ", nameof(typeof(X)), " into a destination that shares its buffers; use `circshift` or a separate destination"))) + _checkconvertible(eltype(O), view(getnzval(X), 1:nz)) + _checknnzfits(indtype(O), nz) ##### readjust output resize!(getcolptr(O), size(X, 2) + 1) resize!(getrowval(O), nz) diff --git a/test/linalg.jl b/test/linalg.jl index d39ca267..7c4cf213 100644 --- a/test/linalg.jl +++ b/test/linalg.jl @@ -59,6 +59,13 @@ include("testhelpers.jl") @test mismatch(O, [1 0 2; 0 3 0; 4 0 5.0]) === nothing end @test mismatch(circshift!(S, S, (1, 0)), [4 0 5.0; 1 0 2; 0 3 0]) === nothing + # a value or a count that the destination cannot hold is rejected before it is written to + O = sparse([1 0 0; 0 2 0; 0 0 3]) + @test_throws InexactError circshift!(O, sparse([1.5 0 0; 0 0 0; 0 0 0]), (0, 1)) + @test mismatch(O, [1 0 0; 0 2 0; 0 0 3]) === nothing + O8 = SparseMatrixCSC{Float64,Int8}(sparse(1.0I, 12, 12)) + @test_throws ArgumentError circshift!(O8, sparse(ones(12, 12)), (0, 1)) + @test mismatch(O8, Matrix(1.0I, 12, 12); Ti=Int8) === nothing end end From 87d4caea2a54353b86b8306665bb49be53f97396 Mon Sep 17 00:00:00 2001 From: "Viral B. Shah" Date: Mon, 26 Oct 2026 14:29:01 +0000 Subject: [PATCH 10/30] Throw an `ArgumentError` from `cholesky` and `ldlt` for an index type CHOLMOD cannot take (#1117) MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit `cholesky` and `ldlt` of a sparse matrix whose index type does not promote with `Int` to one CHOLMOD takes, such as `UInt64`, throw a bare `TypeError` from inside the `CHOLMOD.Sparse` conversion, where `lu` names the matrix type. `cholesky!` and `ldlt!` reject index types like `Int16` that `cholesky` and `ldlt` accept. **Before** ```julia julia> A = sparse([4.0 1 0; 1 4 1; 0 1 4]); julia> cholesky(SparseMatrixCSC{Float64,UInt64}(A)) ERROR: TypeError: in Sparse, in Ti, expected Ti<:Union{Int32, Int64}, got Type{UInt64} julia> cholesky(SparseMatrixCSC{Float64,Int16}(A)); # works julia> cholesky!(cholesky(A), SparseMatrixCSC{Float64,Int16}(A)) ERROR: TypeError: in Sparse, in Ti, expected Ti<:Union{Int32, Int64}, got Type{Int16} ``` **After** ```julia julia> cholesky(SparseMatrixCSC{Float64,UInt64}(A)) ERROR: ArgumentError: matrix type SparseMatrixCSC{Float64, UInt64} not supported: CHOLMOD needs an index type of Int32 or Int64. Try converting its index type to Int first. julia> cholesky!(cholesky(A), SparseMatrixCSC{Float64,Int16}(A)) SparseArrays.CHOLMOD.Factor{Float64, Int64} ... ``` The six `CHOLMOD.Sparse` conversions from a `SparseMatrixCSC` or a `Symmetric`/`Hermitian` view of one computed `Ti <: ITypes ? Ti : promote_type(Ti, Int)` and passed the result to `Sparse{Tv, Ti}` unchecked. They now call a helper, `cholmod_itype`, that does the same promotion and throws the `ArgumentError` when the result is not in `ITypes`; the message lists `ITypes`, which is `Int32` alone on a 32-bit build. `cholesky!(F, A)` and `ldlt!(F, A)` used to convert `A` to `Sparse{eltype(F), Ti}` with `A`'s own index type, which `factorize_p!` then converted again to the factor's. They now convert straight to the factor's index type, so they accept any integer index type, `UInt64` included, as long as the indices fit. This touches dispatch: those two methods now take the index type from `F` rather than from `A`. ๐Ÿค– Generated with [Claude Code](https://claude.com/claude-code) Co-authored-by: Claude Opus 5.5 (cherry picked from commit 31a10e816f1877270f5e0c4f2d88ba499adf4aca) --- src/solvers/cholmod.jl | 51 +++++++++++++++++++++++++----------------- test/solvers/issues.jl | 7 ++++++ 2 files changed, 38 insertions(+), 20 deletions(-) diff --git a/src/solvers/cholmod.jl b/src/solvers/cholmod.jl index 977da253..2e7d2be8 100644 --- a/src/solvers/cholmod.jl +++ b/src/solvers/cholmod.jl @@ -1166,8 +1166,19 @@ function Sparse{Tv, Ti}(A::SparseMatrixCSC{<:Any}, stype::Integer) where {Tv<:VT check_sparse(o) return o end + +# CHOLMOD takes Int32 or Int64 indices; other index types are promoted with Int +function cholmod_itype(A, ::Type{Ti}) where {Ti} + Ti <: ITypes && return Ti + T = promote_type(Ti, Int) + T <: ITypes || throw(ArgumentError(string("matrix type ", typeof(A), " not supported: ", + "CHOLMOD needs an index type of ", join(Base.uniontypes(ITypes), " or "), ". ", + "Try converting its index type to Int first."))) + return T +end + Sparse{Tv}(A::SparseMatrixCSC{<:Any, Ti}, stype::Integer) where {Tv<:VTypes, Ti} = - Sparse{Tv, Ti <: ITypes ? Ti : promote_type(Ti, Int)}(A, stype) + Sparse{Tv, cholmod_itype(A, Ti)}(A, stype) # handle promotion function Sparse(A::SparseMatrixCSC{Tv}, stype::Integer) where {Tv} T = promote_type(Tv, Float64) @@ -1191,7 +1202,7 @@ function Sparse{Tv, Ti}(A::SparseMatrixCSC) where {Tv<:VTypes, Ti<:ITypes} o end Sparse{Tv}(A::SparseMatrixCSC{<:Any, Ti}) where {Tv<:VTypes, Ti} = - Sparse{Tv, Ti <: ITypes ? Ti : promote_type(Ti, Int)}(A) + Sparse{Tv, cholmod_itype(A, Ti)}(A) function Sparse(A::SparseMatrixCSC{Tv}) where {Tv} T = promote_type(Tv, Float64) return Sparse{T}(A) @@ -1206,19 +1217,19 @@ Sparse{Tv, Ti}(A::Hermitian{<:Any, <:SparseMatrixCSC}) where {Tv<:VTypes, Ti<:IT Sparse{Tv, Ti}(A.data, A.uplo == 'L' ? -1 : 1) Sparse(A::Symmetric{Tv, SparseMatrixCSC{Tv,Ti}}) where {Tv<:Real, Ti} = - Sparse{promote_type(Tv, Float64), Ti <: ITypes ? Ti : promote_type(Ti, Int)}( + Sparse{promote_type(Tv, Float64), cholmod_itype(A, Ti)}( A.data, A.uplo == 'L' ? -1 : 1 ) Sparse(A::Symmetric{Tv, SparseMatrixCSC{Tv,Ti}}) where {Tv<:Union{Float16, Float32}, Ti} = - Sparse{Float32, Ti <: ITypes ? Ti : promote_type(Ti, Int)}(A.data, A.uplo == 'L' ? -1 : 1) + Sparse{Float32, cholmod_itype(A, Ti)}(A.data, A.uplo == 'L' ? -1 : 1) Sparse(A::Hermitian{Tv, SparseMatrixCSC{Tv,Ti}}) where {Tv, Ti} = - Sparse{promote_type(Tv, Float64), Ti <: ITypes ? Ti : promote_type(Ti, Int)}( + Sparse{promote_type(Tv, Float64), cholmod_itype(A, Ti)}( A.data, A.uplo == 'L' ? -1 : 1 ) Sparse(A::Hermitian{Tv, SparseMatrixCSC{Tv,Ti}}) where {Tv<:Union{Float16, Float32, ComplexF32, ComplexF16}, Ti} = - Sparse{promote_type(Float32, Tv), Ti <: ITypes ? Ti : promote_type(Ti, Int)}( + Sparse{promote_type(Float32, Tv), cholmod_itype(A, Ti)}( A.data, A.uplo == 'L' ? -1 : 1 ) @@ -1738,13 +1749,13 @@ true Input matrices not of those element types will be converted to these types as appropriate. """ -cholesky!(F::Factor, A::Union{SparseMatrixCSC{T, Ti}, - SparseMatrixCSC{Complex{T}, Ti}, - Symmetric{T,<:SparseMatrixCSC{T, Ti}}, - Hermitian{Complex{T},<:SparseMatrixCSC{Complex{T}, Ti}}, - Hermitian{T,<:SparseMatrixCSC{T, Ti}}}; - shift = 0.0, check::Bool = true) where {T<:Real, Ti} = - cholesky!(F, Sparse{eltype(F), Ti}(A); shift = shift, check = check) +cholesky!(F::Factor{Tv, Ti}, A::Union{SparseMatrixCSC{T}, + SparseMatrixCSC{Complex{T}}, + Symmetric{T,<:SparseMatrixCSC{T}}, + Hermitian{Complex{T},<:SparseMatrixCSC{Complex{T}}}, + Hermitian{T,<:SparseMatrixCSC{T}}}; + shift = 0.0, check::Bool = true) where {Tv, Ti, T<:Real} = + cholesky!(F, Sparse{Tv, Ti}(A); shift = shift, check = check) function cholesky(A::Sparse; shift::Real=0.0, check::Bool = true, perm::Union{Nothing,AbstractVector{<:Integer}}=nothing) @@ -1928,13 +1939,13 @@ true Input matrices not of those element types will be converted to these types as appropriate. """ -ldlt!(F::Factor, A::Union{SparseMatrixCSC{T, Ti}, - SparseMatrixCSC{Complex{T}, Ti}, - Symmetric{T, <:SparseMatrixCSC{T, Ti}}, - Hermitian{Complex{T}, <:SparseMatrixCSC{Complex{T}, Ti}}, - Hermitian{T, <:SparseMatrixCSC{T, Ti}}}; - shift = 0.0, check::Bool = true) where {T<:Real, Ti} = - ldlt!(F, Sparse{eltype(F), Ti}(A), shift = shift, check = check) +ldlt!(F::Factor{Tv, Ti}, A::Union{SparseMatrixCSC{T}, + SparseMatrixCSC{Complex{T}}, + Symmetric{T, <:SparseMatrixCSC{T}}, + Hermitian{Complex{T}, <:SparseMatrixCSC{Complex{T}}}, + Hermitian{T, <:SparseMatrixCSC{T}}}; + shift = 0.0, check::Bool = true) where {Tv, Ti, T<:Real} = + ldlt!(F, Sparse{Tv, Ti}(A), shift = shift, check = check) function ldlt(A::Sparse; shift::Real=0.0, check::Bool = true, perm::Union{Nothing,AbstractVector{<:Integer}}=nothing) diff --git a/test/solvers/issues.jl b/test/solvers/issues.jl index 3674c2f9..588c8acb 100644 --- a/test/solvers/issues.jl +++ b/test/solvers/issues.jl @@ -207,6 +207,13 @@ end @test F \ b โ‰ˆ Matrix(A) \ b end +@testset "cholesky and cholesky! with an index type other than CHOLMOD's" begin + A = sparse([4.0 1 0; 1 4 1; 0 1 4]) + @test_throws ArgumentError cholesky(SparseMatrixCSC{Float64,UInt64}(A)) + F = cholesky(A) + @test cholesky!(F, SparseMatrixCSC{Float64,Int16}(2A)) \ [1.0, 2.0, 3.0] โ‰ˆ (2A) \ [1.0, 2.0, 3.0] +end + @testset "Issues #26367 and #26368" begin A = sparse([0.0 1 0 0; 0 0 0 0]) F = qr(A) From 6cba848aaba4db2f3c0e8e122456f6dbe5c9b30c Mon Sep 17 00:00:00 2001 From: "Viral B. Shah" Date: Mon, 26 Oct 2026 14:34:57 +0000 Subject: [PATCH 11/30] Return a sparse matrix from `diff` of a sparse view, transpose or adjoint (#1118) MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit `diff` of a view, `transpose` or `adjoint` of a sparse matrix reaches the dense `AbstractArray` method and returns a `Matrix`, while `diff(A; dims)` itself is sparse. **Before** ```julia julia> A = sparse([1.0 0 2; 0 3 0; 4 0 5]); julia> typeof(diff(view(A, :, 1:2); dims=1)), typeof(diff(view(A, :, [1, 3]); dims=2)), typeof(diff(view(A, 1:2, :); dims=1)) (Matrix{Float64}, Matrix{Float64}, Matrix{Float64}) julia> typeof(diff(transpose(A); dims=1)), typeof(diff(A'; dims=2)) (Matrix{Float64}, Matrix{Float64}) ``` **After** ```julia julia> typeof(diff(view(A, :, 1:2); dims=1)), typeof(diff(view(A, :, [1, 3]); dims=2)), typeof(diff(view(A, 1:2, :); dims=1)) (SparseMatrixCSC{Float64, Int64}, SparseMatrixCSC{Float64, Int64}, SparseMatrixCSC{Float64, Int64}) julia> typeof(diff(transpose(A); dims=1)), typeof(diff(A'; dims=2)) (SparseMatrixCSC{Float64, Int64}, SparseMatrixCSC{Float64, Int64}) ``` The `diff` method was defined on `AbstractSparseMatrixCSC` only. It now takes `SparseMatrixCSCOrSubArray`: `sparse_diff1` and `sparse_diff2` read each column through `getnzrange` against the parent's `getrowval`/`getnzval` instead of `nzrange` and `getcolptr`, so they run directly on a view of all rows and any column subset (where `getcolptr` is undefined and `nzrange` is O(`col`)), and any other view goes through `_compressed`, its O(nnz) sparse copy. An `adjoint` or `transpose` of one of these takes `diff` of the parent along the other dimension and materializes the wrapper with `copy`, which is also O(nnz). For a `SparseMatrixCSC` the kernels do the same work as before. Vector views were already sparse. This touches dispatch. Potential breakage: `diff` of a sparse view, `transpose` or `adjoint` returned a dense `Matrix` and now returns a `SparseMatrixCSC`, so code that relied on the dense result type (writing into it with `setindex!` expecting dense storage, or dispatching on `Matrix`) sees a different type. ๐Ÿค– Generated with [Claude Code](https://claude.com/claude-code) Co-authored-by: Claude Opus 5.5 (cherry picked from commit e4a42b4ab8c831f39860b610b3a48c4d233e66dd) --- src/linalg.jl | 29 +++++++++++++---------------- test/linalg.jl | 2 ++ 2 files changed, 15 insertions(+), 16 deletions(-) diff --git a/src/linalg.jl b/src/linalg.jl index efd664e5..f4feb7c2 100644 --- a/src/linalg.jl +++ b/src/linalg.jl @@ -725,7 +725,7 @@ function _sparse_diff1_column!(rowval, nzval, numnz, rowval_a, nzval_a, rng, m) return numnz end -function sparse_diff1(S::AbstractSparseMatrixCSC{Tv,Ti}) where {Tv,Ti} +function sparse_diff1(S::SparseMatrixCSCOrColumnSubset{Tv,Ti}) where {Tv,Ti} m,n = size(S) m > 1 || return SparseMatrixCSC(0, n, fill(one(Ti),n+1), Ti[], Tv[]) colptr = Vector{Ti}(undef, n+1) @@ -735,7 +735,7 @@ function sparse_diff1(S::AbstractSparseMatrixCSC{Tv,Ti}) where {Tv,Ti} numnz = 0 @inbounds colptr[1] = 1 @inbounds for col = 1 : n - numnz = _sparse_diff1_column!(rowval, nzval, numnz, getrowval(S), getnzval(S), nzrange(S, col), m) + numnz = _sparse_diff1_column!(rowval, nzval, numnz, getrowval(S), getnzval(S), getnzrange(S, col), m) colptr[col+1] = numnz+1 end deleteat!(rowval, numnz+1:length(rowval)) @@ -743,7 +743,7 @@ function sparse_diff1(S::AbstractSparseMatrixCSC{Tv,Ti}) where {Tv,Ti} return SparseMatrixCSC(m-1, n, colptr, rowval, nzval) end -function sparse_diff2(a::AbstractSparseMatrixCSC{Tv,Ti}) where {Tv,Ti} +function sparse_diff2(a::SparseMatrixCSCOrColumnSubset{Tv,Ti}) where {Tv,Ti} m,n = size(a) colptr = Vector{Ti}(undef, max(n,1)) numnz = 2 * nnz(a) # upper bound; will shrink later @@ -752,7 +752,6 @@ function sparse_diff2(a::AbstractSparseMatrixCSC{Tv,Ti}) where {Tv,Ti} z = zero(Tv) - colptr_a = getcolptr(a) rowval_a = getrowval(a) nzval_a = getnzval(a) @@ -762,28 +761,21 @@ function sparse_diff2(a::AbstractSparseMatrixCSC{Tv,Ti}) where {Tv,Ti} n == 0 && return SparseMatrixCSC(m, n, colptr, rowval, nzval) - startA = colptr_a[1] - stopA = colptr_a[2] - - rA = startA : stopA - 1 + rA = getnzrange(a, 1) rowvalA = rowval_a[rA] nzvalA = nzval_a[rA] - lA = stopA - startA + lA = length(rA) end @inbounds for col = 1:n-1 - startB, stopB = startA, stopA - startA = colptr_a[col+1] - stopA = colptr_a[col+2] - rowvalB = rowvalA nzvalB = nzvalA lB = lA - rA = startA : stopA - 1 + rA = getnzrange(a, col+1) rowvalA = rowval_a[rA] nzvalA = nzval_a[rA] - lA = stopA - startA + lA = length(rA) ptrB = 1 ptrA = 1 @@ -834,10 +826,15 @@ function sparse_diff2(a::AbstractSparseMatrixCSC{Tv,Ti}) where {Tv,Ti} return SparseMatrixCSC(m, n-1, colptr, rowval, nzval) end -function diff(a::AbstractSparseMatrixCSC; dims::Integer) +function diff(a::SparseMatrixCSCOrSubArray; dims::Integer) 1 <= dims <= 2 || throw(ArgumentError("dimension $dims out of range (1:2)")) + a = _compressed(a) return dims == 1 ? sparse_diff1(a) : sparse_diff2(a) end +function diff(a::AdjOrTrans{<:Any,<:SparseMatrixCSCOrSubArray}; dims::Integer) + 1 <= dims <= 2 || throw(ArgumentError("dimension $dims out of range (1:2)")) + return copy(LinearAlgebra.wrapperop(a)(diff(parent(a); dims = 3 - dims))) +end function _sparse_diff(x, dims) dims == 1 || throw(ArgumentError("dimension $dims out of range (1:1)")) diff --git a/test/linalg.jl b/test/linalg.jl index 7c4cf213..d206f5eb 100644 --- a/test/linalg.jl +++ b/test/linalg.jl @@ -555,6 +555,8 @@ end @test (dims == 1 ? 3 : 4) โˆ‰ view(rowvals(D), nzrange(D, dims == 1 ? 2 : 1)) @static if COMPREHENSIVE @test diff(sparse(Int32.(1:7), Int32.(1:7), one(T)); dims) isa SparseMatrixCSC{T,Int32} + @test mismatch(diff(view(A, :, [1, 2, 4]); dims), diff(M[:, [1, 2, 4]]; dims)) === nothing + @test mismatch(diff(transpose(A); dims), diff(transpose(M); dims)) === nothing end end for dims in (0, 3, 7) From f74be420c695f0a9ad17cd7005d1006ddddf66fc Mon Sep 17 00:00:00 2001 From: "Viral B. Shah" Date: Mon, 26 Oct 2026 12:48:55 -0400 Subject: [PATCH 12/30] Honour `T` and the documented `rfn` forms in `sprand` with a value function (#1109) MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit `sprand` with a value function `rfn` mishandles `T`. Without an RNG, `sprand(m, n, p, rfn, T)` ignores `T`. Without `T`, it calls `rfn` extra times just to find the eltype. The vector `sprand([rng], n, p, rfn, T)` rejects the documented `rfn(k)` and `rfn(rng, k)`: it only accepts an undocumented `rfn(rng, T, k)`. **Before** ```julia julia> calls = Ref(0); rfn = k -> (calls[] += 1; rand(Int8, k)); julia> eltype(sprand(3, 3, 0.5, rfn, Float64)), calls[] (Int8, 2) julia> calls[] = 0; eltype(sprand(3, 3, 0.5, rfn)), calls[] (Int8, 3) julia> typeof(sprand(10, 0.5, k -> rand(Int8, k), Float64)) ERROR: MethodError: no method matching (::var"#5#6")(::Random.TaskLocalRNG, ::Type{Float64}, ::Int64) ``` **After** ```julia julia> eltype(sprand(3, 3, 0.5, rfn, Float64)), calls[] (Float64, 1) julia> calls[] = 0; eltype(sprand(3, 3, 0.5, rfn)), calls[] (Int8, 1) julia> typeof(sprand(10, 0.5, k -> rand(Int8, k), Float64)) SparseVector{Float64, Int64} ``` The matrix method without an RNG now passes `T` on. Without `T`, the matrix methods take the eltype from the values `rfn` returns, instead of a throwaway call `rfn(1)` that also drew from the RNG. The vector methods with `T` now call the documented `rfn(k)` or `rfn(rng, k)` and convert to `T`, as the matrix methods do. They keep accepting `rfn(rng, T, k)`, checked with `applicable`. The docstring now says what `T` after `rfn` does. Changed stream: seeded `sprand([rng], m, n, p, rfn)` calls with a custom `rfn` and no `T` draw less from the RNG, so their values change. Every other form passes `T` itself and is unchanged: `sprand([rng], [T], m, n, p)`, `sprandn` and the vector forms, `sprand(Xoshiro(42), 50, 50, 0.3)` among them, which is bitwise identical to `main`. ๐Ÿค– Generated with [Claude Code](https://claude.com/claude-code) --------- Co-authored-by: Claude Opus 5.5 (cherry picked from commit daa0ac09ad883a58caf633b57425d9a933b5365e) --- src/constructors.jl | 18 +++++++++++++++--- src/sparsevector.jl | 11 ++++++++--- test/constructors.jl | 1 + test/sparsevector.jl | 8 ++++++++ 4 files changed, 32 insertions(+), 6 deletions(-) diff --git a/src/constructors.jl b/src/constructors.jl index fbfc69e5..1625213a 100644 --- a/src/constructors.jl +++ b/src/constructors.jl @@ -697,6 +697,7 @@ You can sample nonzero values from a different distribution by passing a custom `rfn` function instead of `rand`. This should be a function `rfn(k)` that returns an array of `k` random numbers sampled from the desired distribution; alternatively, if `rng` is supplied, it should instead be a function `rfn(rng, k)`. +Passing `T` after `rfn` converts the values that `rfn` returns to `T`. # Examples ```jldoctest; setup = :(using Random; Random.seed!(1234)) @@ -711,16 +712,27 @@ julia> sprand(Float64, 3, 0.75) [2] = 0.49425 ``` """ -function sprand(r::AbstractRNG, m::Integer, n::Integer, density::AbstractFloat, rfn::Function, ::Type{T}=eltype(rfn(r, 1))) where T +function sprand(r::AbstractRNG, m::Integer, n::Integer, density::AbstractFloat, rfn::Function, ::Type{T}) where T m, n = Int(m), Int(n) + I, V = _sprandvalues(r, m, n, density, rfn) + return sparse_sortedlinearindices!(I, convert(Vector{T}, V), m, n) +end +function sprand(r::AbstractRNG, m::Integer, n::Integer, density::AbstractFloat, rfn::Function) + m, n = Int(m), Int(n) + I, V = _sprandvalues(r, m, n, density, rfn) + return sparse_sortedlinearindices!(I, convert(Vector{eltype(V)}, V), m, n) +end +function _sprandvalues(r::AbstractRNG, m::Int, n::Int, density::AbstractFloat, rfn::F) where {F} (m < 0 || n < 0) && throw(ArgumentError("invalid Array dimensions")) 0 <= density <= 1 || throw(ArgumentError("$density not in [0,1]")) I = randsubseq(r, 1:(m*n), density) - return sparse_sortedlinearindices!(I, convert(Vector{T}, rfn(r,length(I))), m, n) + return I, rfn(r, length(I)) end -sprand(m::Integer, n::Integer, density::AbstractFloat, rfn::Function, ::Type{T} = eltype(rfn(1))) where {T} = +sprand(m::Integer, n::Integer, density::AbstractFloat, rfn::Function) = sprand(default_rng(), m, n, density, (r, i) -> rfn(i)) +sprand(m::Integer, n::Integer, density::AbstractFloat, rfn::Function, ::Type{T}) where {T} = + sprand(default_rng(), m, n, density, (r, i) -> rfn(i), T) truebools(r::AbstractRNG, n::Integer) = fill(true, n) diff --git a/src/sparsevector.jl b/src/sparsevector.jl index 5d8f997a..893b7f0f 100644 --- a/src/sparsevector.jl +++ b/src/sparsevector.jl @@ -768,11 +768,16 @@ end ### Rand Construction -sprand(n::Integer, p::AbstractFloat, rfn::Function, ::Type{T}) where {T} = sprand(default_rng(), n, p, rfn, T) +# `rfn(rng, T, k)` predates the documented `rfn(k)` and `rfn(rng, k)` and is still accepted +function sprand(n::Integer, p::AbstractFloat, rfn::Function, ::Type{T}) where {T} + r = default_rng() + applicable(rfn, r, T, 0) && return sprand(r, n, p, rfn, T) + return sprand(r, n, p, (r, i) -> rfn(i), T) +end function sprand(r::AbstractRNG, n::Integer, p::AbstractFloat, rfn::Function, ::Type{T}) where T I = randsubseq(r, 1:convert(Int, n), p) - V = rfn(r, T, length(I)) - SparseVector(n, I, V) + V = applicable(rfn, r, T, length(I)) ? rfn(r, T, length(I)) : rfn(r, length(I)) + SparseVector(n, I, convert(Vector{T}, V)) end sprand(n::Integer, p::AbstractFloat, rfn::Function) = sprand(default_rng(), n, p, (r, i) -> rfn(i)) diff --git a/test/constructors.jl b/test/constructors.jl index 1cc66d3d..3293a597 100644 --- a/test/constructors.jl +++ b/test/constructors.jl @@ -290,6 +290,7 @@ end @static if COMPREHENSIVE m = sprand(Float32, 10, 10, 0.1) @test eltype(m) == Float32 + @test eltype(sprand(10, 10, 0.5, k -> rand(Int8, k), Float64)) == Float64 end m = sprand(Float64, 10, 10, 0.1) @test eltype(m) == Float64 diff --git a/test/sparsevector.jl b/test/sparsevector.jl index 303195a4..8628d1a9 100644 --- a/test/sparsevector.jl +++ b/test/sparsevector.jl @@ -242,6 +242,14 @@ end @test xv == xr @test sprand(20, 1, 0.5, rfn1) isa SparseMatrixCSC{Int8,Int} @test sprand(MersenneTwister(5), 20, 1, 0.5, rfn2) isa SparseMatrixCSC{Int8,Int} + # and with `T`, which the older `rfn(rng, T, k)` still accepts + @test sprand(20, 0.5, rfn1, Float64) isa SparseVector{Float64,Int} + @test sprand(MersenneTwister(5), 20, 0.5, rfn2, Float64) isa SparseVector{Float64,Int} + @test sprand(MersenneTwister(5), 20, 0.5, (r, T, k) -> rand(r, T, k), Float64) isa SparseVector{Float64,Int} + # without `T`, `rfn` is called once, for the values + calls = Ref(0) + sprand(MersenneTwister(5), 4, 4, 0.5, (r, k) -> (calls[] += 1; rand(r, Int8, k))) + @test calls[] == 1 end end end From a4039e08920d4c6485f13d633902cad642398c2a Mon Sep 17 00:00:00 2001 From: "Viral B. Shah" Date: Mon, 26 Oct 2026 12:49:50 -0400 Subject: [PATCH 13/30] Return an `n x 1` matrix from `SparseMatrixCSC` of a sparse vector view (#1113) MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit `SparseMatrixCSC(view(A, :, j))` returns a `SparseVector`, so a constructor of `SparseMatrixCSC` returns something that is not a `SparseMatrixCSC`. The same holds for `SparseMatrixCSC{Tv}` and for a view of a sparse vector, while `SparseMatrixCSC{Tv,Ti}` of the same view, and `SparseMatrixCSC` of a `SparseVector`, return an `n x 1` matrix. **Before** ```julia using SparseArrays A = sparse([1.0 0; 0 1]) for C in (SparseMatrixCSC, SparseMatrixCSC{Float64}, SparseMatrixCSC{Float64,Int}) r = C(view(A, :, 1)); println(C, " => ", typeof(r), " ", size(r)) end ``` ``` SparseMatrixCSC => SparseVector{Float64, Int64} (2,) SparseMatrixCSC{Float64} => SparseVector{Float64, Int64} (2,) SparseMatrixCSC{Float64, Int64} => SparseMatrixCSC{Float64, Int64} (2, 1) ``` **After** ``` SparseMatrixCSC => SparseMatrixCSC{Float64, Int64} (2, 1) SparseMatrixCSC{Float64} => SparseMatrixCSC{Float64, Int64} (2, 1) SparseMatrixCSC{Float64, Int64} => SparseMatrixCSC{Float64, Int64} (2, 1) ``` The wrapper constructors in `sparseconvert.jl` return `_sparsem(A)` directly when the eltype matches, and `_sparsem` of a 1-d view is the `SparseVector` that indexing the parent gives. Two methods on `SubArray{<:Any,1}` now pass that vector on to `SparseMatrixCSC`, as the typed `{Tv,Ti}` constructor already does. `sparse(view(A, :, j))` is unchanged and still returns a `SparseVector`. This touches dispatch. **Potential breakage:** this changes the result type of code that works today. A caller relying on `SparseMatrixCSC(view(A, :, j))` being a `SparseVector`, for example indexing it with one index as a vector or passing it where a vector is required, now gets an `n x 1` matrix. `sparse(view(A, :, j))` or `A[:, j]` gives the vector. Dense has no constructor from a vector: `Matrix(v)` and `Matrix(view(M, :, 1))` throw a `MethodError`. A sparse vector has long converted to an `n x 1` matrix, through `SparseMatrixCSC(x)` and `convert(SparseMatrixCSC, x)`, so a view of one now does the same. Matching dense exactly would break those. A view of a dense vector, which dense cannot convert either, threw a `TypeError` from an internal `invoke`. It now throws the `MethodError` that `Matrix` gives: **Before** ```julia julia> SparseMatrixCSC(view([1.0, 2.0], :)) ERROR: TypeError: in invoke: argument type error, expected Tuple{AbstractMatrix}, got a value of type Tuple{SubArray{Float64, 1, Vector{Float64}, Tuple{Base.Slice{Base.OneTo{Int64}}}, true}} ``` **After** ```julia julia> SparseMatrixCSC(view([1.0, 2.0], :)) ERROR: MethodError: no method matching SparseMatrixCSC(::SubArray{Float64, 1, Vector{Float64}, Tuple{Base.Slice{Base.OneTo{Int64}}}, true}) ``` ๐Ÿค– Generated with [Claude Code](https://claude.com/claude-code) --------- Co-authored-by: Claude Opus 5.5 (cherry picked from commit 1a566031e5b1149f723c098537b6bca5cc8a7a47) --- src/sparseconvert.jl | 5 +++++ test/constructors.jl | 3 +++ 2 files changed, 8 insertions(+) diff --git a/src/sparseconvert.jl b/src/sparseconvert.jl index 3f05adcd..2678d4bf 100644 --- a/src/sparseconvert.jl +++ b/src/sparseconvert.jl @@ -11,6 +11,9 @@ for wr in (Symmetric, Hermitian, Transpose, Adjoint, @eval SparseMatrixCSC{Tv}(A::$wr) where Tv = SparseMatrixCSC{Tv}(_sparsem(A)) @eval SparseMatrixCSC{Tv,Ti}(A::$wr) where {Tv,Ti} = SparseMatrixCSC{Tv,Ti}(_sparsem(A)) end +# a vector view converts to an `n x 1` matrix, as a sparse vector does +SparseMatrixCSC(A::SubArray{<:Any,1}) = SparseMatrixCSC(_sparsem(A)) +SparseMatrixCSC{Tv}(A::SubArray{Tv,1}) where Tv = SparseMatrixCSC(_sparsem(A)) """ iswrsparse(::S) @@ -66,6 +69,8 @@ function _sparsem(@nospecialize A::AbstractArray{Tv}) where Tv A end else + # a dense vector has no matrix constructor, as `Matrix(v)` has none + A isa AbstractMatrix || throw(MethodError(SparseMatrixCSC, (A,))) # explicitly call abstract matrix fallback using getindex(A,...) invoke(SparseMatrixCSC{Tv,Int}, Tuple{AbstractMatrix}, A) end diff --git a/test/constructors.jl b/test/constructors.jl index 3293a597..90e150bf 100644 --- a/test/constructors.jl +++ b/test/constructors.jl @@ -30,6 +30,9 @@ end @static if COMPREHENSIVE @test SparseMatrixCSC{ComplexF64}(Array(a)')::SparseMatrixCSC{ComplexF64,Int} == a' @test Array(SparseMatrixCSC{eltype(a), Int8}(a)) == Array(a) + @test mismatch(SparseMatrixCSC(view(a, :, 1)), Array(a)[:, 1:1]) === nothing + # a view of a dense vector has no matrix constructor, as for `Matrix` + @test_throws MethodError SparseMatrixCSC(view([1.0, 2.0], :)) end @test collect(a) == a @test SparseMatrixCSC(a)::typeof(a) == a From 4a80966b3ed8ebe4cf95dcecc9873c6eef8e15d7 Mon Sep 17 00:00:00 2001 From: "Viral B. Shah" Date: Mon, 26 Oct 2026 12:54:14 -0400 Subject: [PATCH 14/30] Check the matrix size in `ldlt!` and `cholesky!` of a `CHOLMOD.Factor` before refactorizing (#1122) MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit `ldlt!(F, A)` on an LLt `CHOLMOD.Factor` converts `F` to a simplicial LDLt before CHOLMOD checks the size of `A`, so a call with a matrix of the wrong size throws and still changes the kind of `F`. **Before** ```julia julia> A = sparse([4.0 1 0; 1 4 1; 0 1 4]); F = cholesky(A); julia> propertynames(F) (:L, :U, :PtL, :UP, :p, :ptr) julia> ldlt!(F, sparse(1.0I, 4, 4)) ERROR: SparseArrays.CHOLMOD.CHOLMODException("A and L dimensions do not match") julia> propertynames(F) (:L, :U, :PtL, :UP, :D, :LD, :DU, :PtLD, :DUP, :p, :ptr) julia> cholesky!(cholesky(A), sparse(1.0I, 4, 4)) ERROR: SparseArrays.CHOLMOD.CHOLMODException("A and L dimensions do not match") ``` **After** ```julia julia> ldlt!(F, sparse(1.0I, 4, 4)) ERROR: DimensionMismatch: matrix has size (4, 4), but the factorization has size (3, 3) julia> propertynames(F) (:L, :U, :PtL, :UP, :p, :ptr) julia> cholesky!(cholesky(A), sparse(1.0I, 4, 4)) ERROR: DimensionMismatch: matrix has size (4, 4), but the factorization has size (3, 3) ``` `ldlt!` calls `change_factor!` before `factorize_p!`, and only the latter reaches CHOLMOD's dimension check. A new `check_refactor_size` compares `size(A)` with `size(F)` under the factor's lock, before anything writes to `F`, in both `ldlt!` and `cholesky!`. `cholesky!` already left `F` unchanged, and gets the check so that the two throw the same error. The eltype is already fixed by dispatch and the stype by `check_hermitian_stype`, and a non-square matrix fails that check, so size is the only precondition missing. Potential breakage: a size mismatch now throws `DimensionMismatch` instead of `CHOLMODException`, from both `ldlt!` and `cholesky!`. ๐Ÿค– Generated with [Claude Code](https://claude.com/claude-code) Co-authored-by: Claude Opus 5.5 (cherry picked from commit 1a0d484c2f84eeded7f28f20f24e6d58c4f47797) --- src/solvers/cholmod.jl | 7 +++++++ test/solvers/issues.jl | 8 ++++++++ 2 files changed, 15 insertions(+) diff --git a/src/solvers/cholmod.jl b/src/solvers/cholmod.jl index 2e7d2be8..faae83ab 100644 --- a/src/solvers/cholmod.jl +++ b/src/solvers/cholmod.jl @@ -1671,6 +1671,11 @@ function check_hermitian_stype(A::Sparse) throw(ArgumentError("sparse matrix is not symmetric/Hermitian")) return nothing end +function check_refactor_size(F::Factor, A::Sparse) + size(A) == size(F) || throw(DimensionMismatch( + "matrix has size $(size(A)), but the factorization has size $(size(F))")) + return nothing +end ## Compute that symbolic factorization only function symbolic(A::Sparse{<:VTypes, Ti}; @@ -1710,6 +1715,7 @@ function cholesky!(F::Factor{Tv}, A::Sparse{Tv}; shift::Real=0.0, check::Bool = true) where Tv check_hermitian_stype(A) @lock F.lock begin + check_refactor_size(F, A) @cholmod_param final_ll = true begin factorize_p!(A, shift, F) end @@ -1901,6 +1907,7 @@ function ldlt!(F::Factor{Tv}, A::Sparse{Tv}; shift::Real=0.0, check::Bool = true) where Tv check_hermitian_stype(A) @lock F.lock begin + check_refactor_size(F, A) change_factor!(F, false, false, true, false) factorize_p!(A, shift, F) check && (issuccess(F) || throw(factorization_exception(F))) diff --git a/test/solvers/issues.jl b/test/solvers/issues.jl index 588c8acb..4d218bc3 100644 --- a/test/solvers/issues.jl +++ b/test/solvers/issues.jl @@ -214,6 +214,14 @@ end @test cholesky!(F, SparseMatrixCSC{Float64,Int16}(2A)) \ [1.0, 2.0, 3.0] โ‰ˆ (2A) \ [1.0, 2.0, 3.0] end +@testset "ldlt! with a matrix of another size leaves an LLt factor unchanged" begin + A = sparse([4.0 1 0; 1 4 1; 0 1 4]) + F = cholesky(A) + @test_throws DimensionMismatch ldlt!(F, sparse(1.0I, 4, 4)) + @test :D โˆ‰ propertynames(F) + @test F \ ones(3) โ‰ˆ A \ ones(3) +end + @testset "Issues #26367 and #26368" begin A = sparse([0.0 1 0 0; 0 0 0 0]) F = qr(A) From 89acc69e9181546954dcec33e433569f9c9cd663 Mon Sep 17 00:00:00 2001 From: "Viral B. Shah" Date: Mon, 26 Oct 2026 12:55:04 -0400 Subject: [PATCH 15/30] Check the buffers in the outer `SparseMatrixCSC` constructor before shortening them (#1123) MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit The outer `SparseMatrixCSC(m, n, colptr, rowval, nzval)` constructor shortens `rowval` and `nzval` that are longer than `m*n`, and it did that before validating them. A call that then threw had already truncated the caller's vectors. **Before** ```julia julia> rv = [1, 1]; nz = [1.0, 2.0]; julia> SparseMatrixCSC(1, 1, [1, 3], rv, nz) ERROR: ArgumentError: Invalid buffers for SparseMatrixCSC construction n=1, colptr=2-element Vector{Int64}, rowval=1-element Vector{Int64}, nzval=1-element Vector{Float64} julia> rv, nz ([1], [1.0]) julia> rv = [1, 1]; nz = [1.0, 2.0]; julia> SparseMatrixCSC(1, 1, Int32[1, 2], rv, nz) ERROR: MethodError: no method matching SparseMatrixCSC{Float64, Int64}(::Int64, ::Int64, ::Vector{Int32}, ::Vector{Int64}, ::Vector{Float64}) julia> rv, nz ([1], [1.0]) ``` **After** ```julia julia> SparseMatrixCSC(1, 1, [1, 3], rv, nz) ERROR: ArgumentError: Invalid buffers for SparseMatrixCSC construction n=1, colptr=2-element Vector{Int64}, rowval=2-element Vector{Int64}, nzval=2-element Vector{Float64} julia> SparseMatrixCSC(1, 1, Int32[1, 2], rv, nz) ERROR: ArgumentError: colptr and rowval must have the same element type, got Int32 and Int64 julia> rv, nz ([1, 1], [1.0, 2.0]) ``` The constructor now checks, before the first `resize!`, that the inner constructor will accept the shortened buffers: `colptr` has `n + 1` entries and `colptr[end] - 1` equals both shortened lengths, and `colptr` and `rowval` have the same element type. A valid call is unchanged: the buffers are still shortened in place and shared with the result. The one visible change is that mismatched index element types now throw an `ArgumentError` naming both types instead of a `MethodError` from the inner constructor; no such call ever succeeded. Left out: `FixedSparseCSC`'s outer constructors do not resize, so they need nothing. ๐Ÿค– Generated with [Claude Code](https://claude.com/claude-code) Co-authored-by: Claude Opus 5.5 (cherry picked from commit b86329c4770e8789e055384116b4315a269ea693) --- src/sparsematrix.jl | 10 ++++++++-- test/constructors.jl | 5 +++++ 2 files changed, 13 insertions(+), 2 deletions(-) diff --git a/src/sparsematrix.jl b/src/sparsematrix.jl index 102db346..d0de6acc 100644 --- a/src/sparsematrix.jl +++ b/src/sparsematrix.jl @@ -56,13 +56,19 @@ end function SparseMatrixCSC(m::Integer, n::Integer, colptr::Vector, rowval::Vector, nzval::Vector) Tv = eltype(nzval) Ti = promote_type(eltype(colptr), eltype(rowval)) + eltype(colptr) === eltype(rowval) || + throw(ArgumentError("colptr and rowval must have the same element type, got $(eltype(colptr)) and $(eltype(rowval))")) sparse_check_Ti(m, n, Ti) sparse_check(n, colptr, rowval, nzval) # silently shorten rowval and nzval to usable index positions. maxlen = abs(widemul(m, n)) isbitstype(Ti) && (maxlen = min(maxlen, typemax(Ti) - 1)) - length(rowval) > maxlen && resize!(rowval, maxlen) - length(nzval) > maxlen && resize!(nzval, maxlen) + if length(rowval) > maxlen || length(nzval) > maxlen + length(colptr) == n + 1 && colptr[end] - 1 == min(length(rowval), maxlen) == min(length(nzval), maxlen) || + throw(ArgumentError("Invalid buffers for SparseMatrixCSC construction n=$n, colptr=$(summary(colptr)), rowval=$(summary(rowval)), nzval=$(summary(nzval))")) + length(rowval) > maxlen && resize!(rowval, maxlen) + length(nzval) > maxlen && resize!(nzval, maxlen) + end SparseMatrixCSC{Tv,Ti}(m, n, colptr, rowval, nzval) end diff --git a/test/constructors.jl b/test/constructors.jl index 90e150bf..9f02a728 100644 --- a/test/constructors.jl +++ b/test/constructors.jl @@ -83,6 +83,11 @@ end @test isequal(sparse(big.([1,1,1,2,2,3,4,5]),big.([1,2,3,2,3,3,4,5]),big.([1,2,4,3,5,6,7,8]), 6, 6), SparseMatrixCSC(6, 6, big.([1,2,4,7,8,9,9]), big.([1,1,2,1,2,3,4,5]), big.([1,2,3,4,5,6,7,8]))) @test sparse(Any[1,2,3], Any[1,2,3], Any[1,1,1]) == sparse([1,2,3], [1,2,3], [1,1,1]) + # an invalid call does not shorten the caller's vectors + rv, nz = [1, 1], [1.0, 2.0] + @test_throws ArgumentError SparseMatrixCSC(1, 1, [1, 3], rv, nz) + @test_throws ArgumentError SparseMatrixCSC(1, 1, Int32[1, 2], rv, nz) + @test rv == [1, 1] && nz == [1.0, 2.0] end # with combine @test sparse([1, 1, 2, 2, 2], [1, 2, 1, 2, 2], 1.0, 2, 2, +) == sparse([1, 1, 2, 2], [1, 2, 1, 2], [1.0, 1.0, 1.0, 2.0], 2, 2) From 67707e0dbabce73b359df783fef0aad1cfb69cb1 Mon Sep 17 00:00:00 2001 From: "Viral B. Shah" Date: Mon, 26 Oct 2026 16:50:30 -0400 Subject: [PATCH 16/30] Leave a valid array when the predicate of `fkeep!`, `droptol!` or `dropzeros!` throws (#1120) MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit The vector `fkeep!`, and the matrix `droptol!` and `dropzeros!`, compact the storage in place with no guard. If the predicate throws partway through, the array is left invalid. Matrix `fkeep!` already handles this through `_fkeep_guarded!`. The predicate can throw on ordinary inputs: `abs(x) > tol` is non-boolean for a `missing` value, and `iszero` throws for an `Any` eltype holding a `String`. ```julia using SparseArrays; using SparseArrays: getcolptr, getrowval, getnzval, nonzeroinds x = sparsevec([1, 2, 3], [-1.0, 1.0, 2.0]) try SparseArrays.fkeep!((i, v) -> i == 3 ? error("boom") : v > 0, x) catch e; println(e) end @show nonzeroinds(x) nonzeros(x) y = sparsevec([1, 2, 3], Any[0, 1, "a"]) try dropzeros!(y) catch e; println(typeof(e)) end @show nonzeroinds(y) nonzeros(y) A = sparse([1, 2], [1, 2], Union{Missing,Float64}[0.1, missing], 2, 2) try droptol!(A, 0.5) catch e; println(typeof(e)) end @show getcolptr(A) getrowval(A) getnzval(A) B = sparse([1, 1, 2], [1, 2, 2], Any[0, 1, "a"], 2, 2) try dropzeros!(B) catch e; println(typeof(e)) end @show getcolptr(B) getrowval(B) getnzval(B) ``` **Before**: the vectors are left with a repeated index. The matrices are left with a stale column pointer, so in `A` the entry `0.1` moves from column 1 to column 2. ``` ErrorException("boom") nonzeroinds(x) = [2, 2, 3] nonzeros(x) = [1.0, 1.0, 2.0] MethodError nonzeroinds(y) = [2, 2, 3] nonzeros(y) = Any[1, 1, "a"] TypeError getcolptr(A) = [1, 1, 3] getrowval(A) = [1, 2] getnzval(A) = Union{Missing, Float64}[0.1, missing] MethodError getcolptr(B) = [1, 1, 4] getrowval(B) = [1, 1, 2] getnzval(B) = Any[1, 1, "a"] ``` **After**: the same errors are thrown, and every array is left valid. ``` ErrorException("boom") nonzeroinds(x) = [2, 3] nonzeros(x) = [1.0, 2.0] MethodError nonzeroinds(y) = [2, 3] nonzeros(y) = Any[1, "a"] TypeError getcolptr(A) = [1, 1, 2] getrowval(A) = [2] getnzval(A) = Union{Missing, Float64}[missing] MethodError getcolptr(B) = [1, 1, 3] getrowval(B) = [1, 2] getnzval(B) = Any[1, "a"] ``` The compaction writes each kept entry at a write position that trails the read position. When the predicate throws, the entries between those two positions are already copied forward but the storage has not been trimmed. For a matrix, the column pointers after the current column are stale too. The vector kernel now records `(read position, write position)` per entry, the same way the matrix kernel records its progress. A new vector method of `_fkeep_guarded!` then does what the matrix method does: it shifts the unvisited tail down, keeping the entry that threw, trims the storage, and rethrows. So the entries visited before the throw are filtered and the rest are kept. The matrix `droptol!` and `dropzeros!` now call `_fkeep_guarded!` as well. `tril!` and `triu!` stay on the unguarded kernel, because their index predicates cannot throw. Recording the progress costs one store per entry. Minimum time of `dropzeros!` over 100 runs, on a vector and on a 10^5ร—10^5 matrix, each with 10^6 stored entries: | case | `main` | branch | |---|---|---| | vector, half stored zeros | 0.510 ms | 0.512 ms | | vector, no stored zeros | 0.595 ms | 0.594 ms | | matrix, half stored zeros | 0.707 ms | 0.627 ms | | matrix, no stored zeros | 0.792 ms | 0.737 ms | The difference is within run-to-run noise: across three rounds, the matrix case on `main` ranged from 0.689 to 0.716 ms. ๐Ÿค– Generated with [Claude Code](https://claude.com/claude-code) Co-authored-by: Claude Opus 5.5 (cherry picked from commit bfa76ee9674734aaad8808433bab3b4db69b75c0) --- src/sparsematrix.jl | 8 ++++---- src/sparsevector.jl | 25 +++++++++++++++++++++++-- test/sparsematrix.jl | 3 +++ test/sparsevector.jl | 6 ++++++ 4 files changed, 36 insertions(+), 6 deletions(-) diff --git a/src/sparsematrix.jl b/src/sparsematrix.jl index d0de6acc..12787a2e 100644 --- a/src/sparsematrix.jl +++ b/src/sparsematrix.jl @@ -1544,8 +1544,8 @@ Base.sort(A::AbstractSparseMatrixCSC; kws...) = ## fkeep! and children tril!, triu!, droptol!, dropzeros[!] -# The predicates of `tril!`, `droptol!` and the other children do not throw, and call this -# with no `progress`. `fkeep!` passes one for a caller's `f`: see `_fkeep_guarded!`. +# The predicates of `tril!` and `triu!` do not throw, and call this with no `progress`. Those +# of `fkeep!`, `droptol!` and `dropzeros!` can, for some eltypes: see `_fkeep_guarded!`. function _fkeep!(f::F, A::AbstractSparseMatrixCSC, progress::P=nothing) where {F<:Function,P} An = size(A, 2) Acolptr = getcolptr(A) @@ -1689,7 +1689,7 @@ julia> droptol!(A, 1e-6) ``` """ droptol!(A::AbstractSparseMatrixCSC, tol) = - (_checkpatternwritable(A, :droptol!); _fkeep!((i, j, x) -> abs(x) > tol, A)) + (_checkpatternwritable(A, :droptol!); _fkeep_guarded!((i, j, x) -> abs(x) > tol, A)) """ dropzeros!(A::AbstractSparseMatrixCSC) @@ -1698,7 +1698,7 @@ Remove the stored zeros from `A`, in place, and return `A`. See also [`dropzeros`](@ref), which returns a copy, and [`droptol!`](@ref). """ -dropzeros!(A::AbstractSparseMatrixCSC) = (_checkpatternwritable(A, :dropzeros!); _fkeep!((i, j, x) -> _isnotzero(x), A)) +dropzeros!(A::AbstractSparseMatrixCSC) = (_checkpatternwritable(A, :dropzeros!); _fkeep_guarded!((i, j, x) -> _isnotzero(x), A)) """ dropzeros(A::AbstractSparseMatrixCSC;) diff --git a/src/sparsevector.jl b/src/sparsevector.jl index 893b7f0f..aae2ad4e 100644 --- a/src/sparsevector.jl +++ b/src/sparsevector.jl @@ -1894,8 +1894,7 @@ end Base.sort(x::AbstractCompressedVector; kws...) = sort!(_writablecopy(x); kws...) -function fkeep!(f, x::AbstractCompressedVector{Tv}) where Tv - _checkpatternwritable(x, :fkeep!) +function _fkeep!(f::F, x::AbstractCompressedVector, progress::Base.RefValue{Tuple{Int,Int}}) where F nzind = nonzeroinds(x) nzval = nonzeros(x) @@ -1903,6 +1902,7 @@ function fkeep!(f, x::AbstractCompressedVector{Tv}) where Tv @inbounds for xk in 1:nnz(x) xi = nzind[xk] xv = nzval[xk] + progress[] = (xk, x_writepos) # If this element should be kept, rewrite in new position if f(xi, xv) if x_writepos != xk @@ -1920,7 +1920,28 @@ function fkeep!(f, x::AbstractCompressedVector{Tv}) where Tv return x end +# As for a matrix: when `f` throws at entry `xk`, keep it and every later entry, and rethrow. +function _fkeep_guarded!(f::F, x::AbstractCompressedVector) where F + progress = Ref((1, 1)) + try + _fkeep!(f, x, progress) + catch + xk, x_writepos = progress[] + nzind, nzval = nonzeroinds(x), nonzeros(x) + shift = xk - x_writepos + if shift != 0 + x_nnz = length(nzind) + copyto!(nzind, x_writepos, nzind, xk, x_nnz - xk + 1) + copyto!(nzval, x_writepos, nzval, xk, x_nnz - xk + 1) + resize!(nzind, x_nnz - shift) + resize!(nzval, x_nnz - shift) + end + rethrow() + end + return x +end +fkeep!(f::F, x::AbstractCompressedVector) where {F} = (_checkpatternwritable(x, :fkeep!); _fkeep_guarded!(f, x)) """ droptol!(x::AbstractCompressedVector, tol) diff --git a/test/sparsematrix.jl b/test/sparsematrix.jl index 51bad629..0e470d02 100644 --- a/test/sparsematrix.jl +++ b/test/sparsematrix.jl @@ -1283,6 +1283,9 @@ end A = sparse(D) @test_throws ErrorException SparseArrays.fkeep!((i, j, x) -> error("stop"), A) @test mismatch(A, D) === nothing + A = sparse([1, 2], [1, 2], Union{Missing,Float64}[0.1, missing], 2, 2) + @test_throws TypeError droptol!(A, 0.5) + @test isequal(A, sparse([2], [2], Union{Missing,Float64}[missing], 2, 2)) && nnz(A) == 1 end @testset "tril! and triu! with a diagonal out of range" begin diff --git a/test/sparsevector.jl b/test/sparsevector.jl index 8628d1a9..b424b7c2 100644 --- a/test/sparsevector.jl +++ b/test/sparsevector.jl @@ -1438,6 +1438,12 @@ end SparseArrays.fkeep!(f_drop, xdrop) @test exact_equal(xdrop, SparseVector(7, [1, 3, 4, 7], [3., -1., 1., 3.])) @test_deprecated SparseArrays.fkeep!(xdrop, f_drop) + @static if COMPREHENSIVE + # the entries visited before the throw are filtered and the others are kept + x = sparsevec([1, 2, 3], [-1.0, 1.0, 2.0]) + @test_throws ErrorException SparseArrays.fkeep!((i, v) -> i == 3 ? error("stop") : v > 0, x) + @test mismatch(x, [0.0, 1.0, 2.0]) === nothing && nnz(x) == 2 + end end @static if COMPREHENSIVE From d9950c7eda716b3dfd3c8ce775c22b4e64e5e103 Mon Sep 17 00:00:00 2001 From: "Viral B. Shah" Date: Mon, 26 Oct 2026 16:59:09 -0400 Subject: [PATCH 17/30] Give sparse `\` and `qr(A) \ b` the eltype of dense `\` (#1124) MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit The eltype of sparse `A \ B` depended on the shape of `A`: a tall or wide `Float32` matrix with a `Float64` right-hand side gave `Float32`, where a square one and dense `\` give `Float64`, and an `Int` tall matrix with a `Float32` right-hand side gave `Float64` where dense gives `Float32`. So `\` was not inferable, and the least squares solve rounded `b` to `Float32`. The same rounding happened at the factorization level in `qr(A) \ b` and `qr(A)' \ b`, and `transpose(A) \ B` of a complex Hermitian `ComplexF32` matrix with a matrix `B` gave `ComplexF64`. ```julia using SparseArrays, LinearAlgebra for A in (Float32[1 0; 0 1; 1 1], Float32[2 0; 1 1], Float32[1 0 1; 0 1 1]) b = collect(1.0:size(A, 1)) println(size(A), ": sparse ", eltype(sparse(A) \ b), ", dense ", eltype(A \ b)) end println(Base.infer_return_type(\, Tuple{SparseMatrixCSC{Float32,Int}, Vector{Float64}})) A = sparse(Float32[1 0; 0 1; 1 1]); b = [1.0, 2.0, 3.0 + 1e-7] x = A \ b println(norm(Matrix{Float64}(A)' * (Matrix{Float64}(A) * x - b))) # normal-equations residual ``` **Before** ``` (3, 2): sparse Float32, dense Float64 (2, 2): sparse Float64, dense Float64 (2, 3): sparse Float32, dense Float64 Union{Vector{Float32}, Vector{Float64}} 3.6563210698207683e-7 ``` **After** ``` (3, 2): sparse Float64, dense Float64 (2, 2): sparse Float64, dense Float64 (2, 3): sparse Float64, dense Float64 Vector{Float64} 1.1102230246251565e-15 ``` ```julia using SparseArrays, LinearAlgebra A = sparse(ComplexF32[2 1im; -1im 2]); B = Float32[1 3; 2 4] println("transpose(A) \\ B: sparse ", eltype(transpose(A) \ B), ", dense ", eltype(transpose(Matrix(A)) \ B)) println(Base.infer_return_type(\, Tuple{typeof(transpose(A)), Matrix{Float32}})) T = sparse(Float32[1 0; 0 1; 1 1]); b = [1.0, 2.0, 3.0 + 1e-7] println("qr(A) \\ b: sparse ", eltype(qr(T) \ b), ", dense ", eltype(qr(Matrix(T)) \ b)) println("qr(A)' \\ b: sparse ", eltype(qr(T)' \ b[1:2]), ", dense ", eltype(qr(Matrix(T))' \ b[1:2])) ``` **Before** ``` transpose(A) \ B: sparse ComplexF64, dense ComplexF32 Union{Matrix{ComplexF32}, Matrix{ComplexF64}} qr(A) \ b: sparse Float32, dense Float64 qr(A)' \ b: sparse Float32, dense Float64 ``` **After** ``` transpose(A) \ B: sparse ComplexF32, dense ComplexF32 Matrix{ComplexF32} qr(A) \ b: sparse Float64, dense Float64 qr(A)' \ b: sparse Float64, dense Float64 ``` Only the LU branch of `_sparse_factorize` converted its solution to `typeof(one(eltype(A)) \ one(eltype(B)))`. The QR and LQ branches, including the `qr(A')` branch of the adjoint and transpose methods, returned the eltype of the factor, which SPQR gives the eltype of the matrix. They now convert the same way. Converting the result alone would keep the `Float32` solve, since `qr` of a `Float32` matrix factorizes in double precision and rounds the factor back, and its solve converts `B` to the factor's eltype. So before the QR or LQ factorization, a matrix of a narrower eltype than the solution is converted to the wider one. That costs nothing extra, because `qr` makes that copy anyway. A matrix and right-hand side of the same eltype factorize as before. The conversion is skipped for an abstract eltype, which has no `one`, so `tall \ Any[...]` still works. At the factorization level, `qr(A) \ b` and `qr(A)' \ b` (and `transpose(qr(A)) \ b`, which goes through the adjoint) converted `b` to the eltype of the factor. When dense `\` gives a wider solution, they now convert the factorization to that precision with the existing `QRSparse{T}(F)` constructor, as dense `\` does with `Factorization{T}(F)`. This copies the factor's arrays exactly and does not refactorize, so the accuracy is that of the `Float32` factorization, as for dense. A real factorization stays real for a complex right-hand side and keeps the reinterpreting solve. The factor is not converted when the target is not a BLAS float, such as `BigFloat`, or when `b` has an abstract eltype; those cases behave as before. `A \ b` does not take this path, because its factor already has the solution's precision. The transposed complex Hermitian `\` method was defined for a `Vector` only, so a matrix right-hand side went to LinearAlgebra's generic `\`. It now takes `VecOrMat`, and the matrix goes through the sparse Hermitian solve as the vector does. This touches dispatch. Potential breakage: code that works today gets a different result eltype. A tall or wide `Float32` matrix with a `Float64` right-hand side returned `Vector{Float32}` from `A \ b`, `qr(A) \ b` and `qr(A)' \ b`, and now returns `Vector{Float64}`, with the solve done in the promoted precision. `transpose(A) \ B` of a complex Hermitian `ComplexF32` matrix with a matrix `B` returned `Matrix{ComplexF64}` and now returns `Matrix{ComplexF32}`. Left out: at the factorization level, `lu` and `cholesky` still differ from dense. `lu` of a `Float32` or `ComplexF32` matrix is an `UmfpackLU{Float64}` or `UmfpackLU{ComplexF64}`, so `lu(A) \ b` returns double precision where dense keeps single. The factor does not record the original eltype, so this cannot be fixed by converting the solution. `cholesky(A) \ b` of a `Float32` matrix keeps a single-precision CHOLMOD factor and returns `Float32` for a `Float64` `b`, where dense returns `Float64`. Converting a CHOLMOD factor is not the same small change. `A \ b` is not affected by either. #999 changes the sparse right-hand side `\` beside this code. Its rectangular branch calls the dense right-hand side `\`, so it gets this fix, and the diffs touch different lines. ๐Ÿค– Generated with [Claude Code](https://claude.com/claude-code) --------- Co-authored-by: Claude Opus 5.5 (cherry picked from commit 5d9ed81b122dba4d2a48def8d4d42e70994a7c65) --- src/linalg.jl | 40 +++++++++++++++++++++++++--------------- src/solvers/spqr.jl | 14 ++++++++++++-- test/solvers/solvers.jl | 9 ++++++++- test/solvers/spqr.jl | 2 ++ 4 files changed, 47 insertions(+), 18 deletions(-) diff --git a/src/linalg.jl b/src/linalg.jl index f4feb7c2..03807a5b 100644 --- a/src/linalg.jl +++ b/src/linalg.jl @@ -554,7 +554,7 @@ function _dot(x::AbstractSparseVector, A::SparseMatrixCSCOrColumnSubset, y::Abst end ## end of symmetric/Hermitian -\(A::Transpose{<:Complex,<:Hermitian{<:Complex,<:AbstractSparseMatrixCSC}}, B::Vector) = copy(A) \ B +\(A::Transpose{<:Complex,<:Hermitian{<:Complex,<:AbstractSparseMatrixCSC}}, B::VecOrMat) = copy(A) \ B ## diagonal scaling @@ -1436,10 +1436,11 @@ end _hermitian_solve(A::AbstractSparseMatrixCSC) = !(eltype(A) <: Union{Integer, Complex{<:Integer}}) && ishermitian(A) -# `\` and `factorize` choose the same method from the structure of `A`; `f` (or `flu`, for -# the LU factorization) is applied in each branch so that the result stays inferrable. A -# Hermitian matrix is only wrapped, because `\` and `factorize` finish it differently. -function _sparse_factorize(f, A::AbstractSparseMatrixCSC, flu = f) +# `\` and `factorize` choose the same method from the structure of `A`; `f` (or `fsolver`, +# for the LU, QR and LQ factorizations) is applied in each branch so that the result stays +# inferrable, and the QR and LQ factorizations are of `fqr(A)`. A Hermitian matrix is only +# wrapped, because `\` and `factorize` finish it differently. +function _sparse_factorize(f, A::AbstractSparseMatrixCSC, fsolver = f, fqr = identity) m, n = size(A) if m == n if istril(A) @@ -1449,12 +1450,12 @@ function _sparse_factorize(f, A::AbstractSparseMatrixCSC, flu = f) elseif _hermitian_solve(A) return f(Hermitian(A)) end - return flu(lu(A)) + return fsolver(lu(A)) elseif m > n - return f(qr(A)) + return fsolver(qr(fqr(A))) else # A is wide, so the LQ factorization gives the minimum-norm solution - return f(lq(A)) + return fsolver(lq(fqr(A))) end end @@ -1462,15 +1463,24 @@ end _solve_factor(F) = F _solve_factor(D::Diagonal) = Diagonal(Vector(diag(D))) _sparse_solve(op, F, B) = \(op(_solve_factor(F)), B) -# the sparse LU works in double precision, so its solution is converted to the eltype that -# dense `\` gives -_sparse_lusolve(op, F, A, B) = - convert(AbstractArray{typeof(one(eltype(A)) \ one(eltype(B)))}, _sparse_solve(op, F, B)) +# The eltype that dense `\` gives. It is `nothing` for an abstract eltype, which has no `one` +# and whose solution is returned as the solver gives it. +_ldiv_eltype(A, B) = isconcretetype(eltype(A)) && isconcretetype(eltype(B)) ? + typeof(one(eltype(A)) \ one(eltype(B))) : nothing +# The sparse LU works in double precision and the QR and LQ in the precision of the matrix, +# so their solution is converted to the eltype that dense `\` gives, and for the QR and LQ a +# matrix of a narrower eltype than that is converted before it is factorized. +_sparse_solversolve(op, F, A, B) = _convert_solution(_ldiv_eltype(A, B), _sparse_solve(op, F, B)) +_convert_solution(T, X) = convert(AbstractArray{T}, X) +_convert_solution(::Nothing, X) = X +_qr_matrix(A, B) = _qr_matrix(A, _ldiv_eltype(A, B)) +_qr_matrix(A, T::Type) = convert(AbstractMatrix{promote_type(float(eltype(A)), real(T))}, A) +_qr_matrix(A, ::Nothing) = A function \(A::AbstractSparseMatrixCSC, B::AbstractVecOrMat) require_one_based_indexing(A, B) return _sparse_factorize(F -> _sparse_solve(identity, F, B), A, - F -> _sparse_lusolve(identity, F, A, B)) + F -> _sparse_solversolve(identity, F, A, B), A -> _qr_matrix(A, B)) end for (xformtype, xformop) in ((:Adjoint, :adjoint), (:Transpose, :transpose)) @eval function \(xformA::($xformtype){<:Any,<:AbstractSparseMatrixCSC}, B::AbstractVecOrMat) @@ -1479,12 +1489,12 @@ for (xformtype, xformop) in ((:Adjoint, :adjoint), (:Transpose, :transpose)) m, n = size(A) if m < n # A' is tall, so the least squares solve needs a factorization of A' itself - return _sparse_solve(identity, qr($xformop(A)), B) + return _sparse_solversolve(identity, qr($xformop(_qr_matrix(A, B))), A, B) end # Otherwise the transformed choice for A serves. For a wide A' that is the # factorization of A, which gives the minimum-norm solution. return _sparse_factorize(F -> _sparse_solve($xformop, F, B), A, - F -> _sparse_lusolve($xformop, F, A, B)) + F -> _sparse_solversolve($xformop, F, A, B), A -> _qr_matrix(A, B)) end end diff --git a/src/solvers/spqr.jl b/src/solvers/spqr.jl index cbfd69ef..6c565685 100644 --- a/src/solvers/spqr.jl +++ b/src/solvers/spqr.jl @@ -654,6 +654,15 @@ LinearAlgebra.rank(F::AdjointQRSparse) = rank(parent(F)) _rhs_eltype(::Type{T}, ::Type{<:Complex}) where {T<:LinearAlgebra.BlasReal} = Complex{T} _rhs_eltype(::Type{T}, ::Type) where {T} = T _rhs_array(F, B::AbstractVecOrMat) = convert(Array{_rhs_eltype(eltype(F), eltype(B))}, B) +# Rather than round the right-hand side to the precision of the factorization, a solve whose +# solution dense `\` gives a wider eltype converts the factorization to that precision, as +# dense `\` does. A real factorization stays real for a complex right-hand side. +_solve_eltype(F, B) = isconcretetype(eltype(B)) ? + promote_type(eltype(F), real(typeof(zero(eltype(B)) / oneunit(eltype(F))))) : eltype(F) +_promoted_factor(F::QRSparse, B) = _promoted_factor(F, _solve_eltype(F, B)) +_promoted_factor(F::QRSparse{T}, ::Type{T}) where {T} = F +_promoted_factor(F::QRSparse, ::Type{T}) where {T<:LinearAlgebra.BlasFloat} = QRSparse{T}(F) +_promoted_factor(F::QRSparse, ::Type) = F ## Two helper methods _ret_size(F::Union{QRSparse,AdjointQRSparse}, b::AbstractVector) = (size(F, 2),) @@ -733,7 +742,7 @@ julia> qr(A)\\fill(1.0, 4) 0.0 ``` """ -(\)(F::QRSparse, B::AbstractVecOrMat) = F\_rhs_array(F, B) +(\)(F::QRSparse, B::AbstractVecOrMat) = (G = _promoted_factor(F, B); G\_rhs_array(G, B)) function LinearAlgebra.ldiv!(X::StridedVecOrMat{T}, F::QRSparse{T}, B::StridedVecOrMat{T}; workspace::Union{Nothing,SpqrWS{T}} = nothing) where {T} @@ -842,7 +851,8 @@ julia> A'x 0.0 ``` """ -(\)(Fadj::AdjointQRSparse, B::AbstractVecOrMat) = Fadj\_rhs_array(Fadj, B) +(\)(Fadj::AdjointQRSparse, B::AbstractVecOrMat) = + (G = adjoint(_promoted_factor(parent(Fadj), B)); G\_rhs_array(G, B)) function LinearAlgebra.ldiv!(X::StridedVecOrMat{T}, Fadj::AdjointQRSparse{T}, B::StridedVecOrMat{T}; workspace::Union{Nothing,SpqrWS{T}} = nothing) where {T} diff --git a/test/solvers/solvers.jl b/test/solvers/solvers.jl index 7f763e38..21cbb8d1 100644 --- a/test/solvers/solvers.jl +++ b/test/solvers/solvers.jl @@ -129,7 +129,14 @@ end # the adjoint solve converts its result to the eltype the plain solve gives x = square' \ Float32.(b) @test x isa Vector{Float32} && x โ‰ˆ Matrix(square') \ b - # only the LU solve converts its result, so a least squares solve takes any right-hand side + # the least squares solve gives the eltype of dense `\`, and keeps the precision of `b` + x = @inferred Float32.(tall) \ b + @test x isa Vector{Float64} && x โ‰ˆ Matrix(tall) \ b + # so does the solve of a transposed complex Hermitian matrix with a matrix right-hand side + Ht = transpose(sparse(ComplexF32[4 im; -im 3])) + @test @inferred(Ht \ Float32[1 2; 3 4]) isa Matrix{ComplexF32} + # the solution for an abstract eltype is not converted, so a least squares solve takes + # any right-hand side @test tall \ Any[1.0, 2, 3] โ‰ˆ Matrix(tall) \ b end # when the Cholesky factorization of a Hermitian matrix fails, `factorize` returns diff --git a/test/solvers/spqr.jl b/test/solvers/spqr.jl index 2417d4a4..8adbd6ca 100644 --- a/test/solvers/spqr.jl +++ b/test/solvers/spqr.jl @@ -254,6 +254,8 @@ end b, B = collect(1.0:m), reshape(collect(1.0:3m), 3, m) eltyA <: Real && @test F.Q * b โ‰ˆ F.Q * eltyA.(b) eltyA <: Real && @test B * F.Q' โ‰ˆ eltyA.(B) * F.Q' + # a solve converts F to the precision of the solution, as dense `\` does + @test @inferred(F \ b) isa Vector{widen(eltyA)} end end From e9eb8a7931fec8d4db7f21aace35ef7280344506 Mon Sep 17 00:00:00 2001 From: "Viral B. Shah" Date: Mon, 26 Oct 2026 17:00:53 -0400 Subject: [PATCH 18/30] Throw before densifying `b` when a triangular `ldiv!` solution does not convert to its eltype (#1125) MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit `ldiv!` of a strided triangular matrix and a sparse-vector `b` stores the whole active range of `b` before it solves, so a solution that does not convert to `eltype(b)` throws with `b` already densified. **Before** ```julia julia> b = sparsevec([2], [1], 2); julia> ldiv!(UpperTriangular([2.0 1; 0 2]), b) ERROR: InexactError: Int64(0.5) julia> (SparseArrays.nonzeroinds(b), nonzeros(b)) ([1, 2], [0, 1]) ``` **After** ```julia julia> ldiv!(UpperTriangular([2.0 1; 0 2]), b) ERROR: InexactError: Int64(-0.25) julia> (SparseArrays.nonzeroinds(b), nonzeros(b)) ([2], [1]) julia> ldiv!(UpperTriangular([2.0 2; 0 2]), sparsevec([2], [2], 2)) 2-element SparseVector{Int64, Int64} with 2 stored entries: [1] = -1 [2] = 1 ``` The eltype computation of `\` moves into `_trisolve_eltype`. When `_storesafely(eltype(b), _trisolve_eltype(xA, b))` fails, as for an integer `b` with a float triangle or a real `b` with a complex one, `ldiv!` solves with `\`, runs `_checkconvertible` on the active range, and only then densifies `b` and copies the solution in. A solution whose values do convert is still stored, as with a dense `b`. The common case, where the solution's eltype is safe to store, keeps the in-place block solve. The `InexactError` now names the first value that fails, not the last one solved. Left out: the literal `0` that `_densify!` writes into new positions stays, because `_densezero` fails for an `Any` eltype, which the in-place solve supports. ๐Ÿค– Generated with [Claude Code](https://claude.com/claude-code) Co-authored-by: Claude Opus 5.5 (cherry picked from commit 8b27c32a4ffe17e76e40c4da5fa415721b05360a) --- src/sparsevector.jl | 18 +++++++++++++++--- test/triangular.jl | 4 ++++ 2 files changed, 19 insertions(+), 3 deletions(-) diff --git a/src/sparsevector.jl b/src/sparsevector.jl index aae2ad4e..a5d4fe42 100644 --- a/src/sparsevector.jl +++ b/src/sparsevector.jl @@ -1795,11 +1795,15 @@ function _checktrisolve(xA::UpperOrLowerTriangular, b::AbstractVector) return nothing end -function \(xA::StridedTriangularMaybeAdjOrTrans{TA}, b::AbstractCompressedVector{Tb}) where {TA<:Number,Tb<:Number} - _checktrisolve(xA, b) - TAb = xA isa UnitUpperOrUnitLowerTriangular ? +# the eltype of the solution of a triangular solve: a unit triangle does not divide +_trisolve_eltype(xA::UpperOrLowerTriangular{TA}, ::AbstractVector{Tb}) where {TA<:Number,Tb<:Number} = + xA isa UnitUpperOrUnitLowerTriangular ? typeof(zero(TA)*zero(Tb) + zero(TA)*zero(Tb)) : typeof((zero(TA)*zero(Tb) + zero(TA)*zero(Tb))/one(TA)) + +function \(xA::StridedTriangularMaybeAdjOrTrans{TA}, b::AbstractCompressedVector{Tb}) where {TA<:Number,Tb<:Number} + _checktrisolve(xA, b) + TAb = _trisolve_eltype(xA, b) r = convert(Array{TAb}, b) if nnz(b) != 0 active = _activerange(xA, b) @@ -1812,6 +1816,14 @@ function ldiv!(xA::StridedTriangularMaybeAdjOrTrans, b::AbstractCompressedVector _checktrisolve(xA, b) nnz(b) == 0 && return b active = _activerange(xA, b) + # a solution that may not convert to the eltype of `b` is solved out of place and checked, + # so that a failed conversion throws before `b` takes the dense pattern + if eltype(xA) <: Number && eltype(b) <: Number && !_storesafely(eltype(b), _trisolve_eltype(xA, b)) + x = view(xA \ b, active) + _checkconvertible(eltype(b), x) + copyto!(view(nonzeros(b), _densify!(b, active)), x) + return b + end # store the whole active section up front rather than entry by entry during the solve, # which also lets the solve run on a dense vector and reach tuned BLAS stored = _densify!(b, active) diff --git a/test/triangular.jl b/test/triangular.jl index 7a93b7c4..de29ab7c 100644 --- a/test/triangular.jl +++ b/test/triangular.jl @@ -730,6 +730,10 @@ end b = SparseVector(3, [1, 3], Any[1.0, 2.0]) x = ldiv!(L, copy(b)) @test nonzeroinds(x) == 1:3 && x โ‰ˆ L \ [1.0, 0.0, 2.0] + # a solution that does not convert to the eltype of `b` throws before `b` is densified + b = sparsevec([2], [1], 2) + @test_throws InexactError ldiv!(UpperTriangular([2.0 1; 0 2]), b) + @test nonzeroinds(b) == [2] && nonzeros(b) == [1] end @testset "scalar promotion" for eltypemat in eltypes From 88e4cf47771b9d21847530df0fadde3c2d8bf14f Mon Sep 17 00:00:00 2001 From: "Viral B. Shah" Date: Mon, 26 Oct 2026 17:01:57 -0400 Subject: [PATCH 19/30] Fall back to `lu` in `factorize` of a sparse Hermitian matrix whose LDLt fails (#1116) MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit `factorize` of a sparse indefinite Hermitian matrix threw `ZeroPivotException` when its LDLt needs pivoting. `\` on the same matrix works, and so does dense `factorize`. `cond(A, 1)` threw too, because it goes through `factorize`. **Before** ```julia julia> A = sparse([0.0 1; 1 0]); b = [1.0, 2.0]; julia> A \ b 2-element Vector{Float64}: 2.0 1.0 julia> factorize(A) \ b ERROR: ZeroPivotException: factorization encountered one or more zero pivots. Consider switching to a pivoted LU factorization. julia> cond(A, 1) ERROR: ZeroPivotException: factorization encountered one or more zero pivots. Consider switching to a pivoted LU factorization. ``` **After** ```julia julia> typeof(factorize(A)) SparseArrays.UMFPACK.UmfpackLU{Float64, Int64} julia> factorize(A) \ b 2-element Vector{Float64}: 2.0 1.0 julia> cond(A, 1) 1.0 ``` The Hermitian method tries `cholesky(A; check = false)` first. If that fails, it called `ldlt!(F, A)`, which throws when CHOLMOD's unpivoted LDLt meets a zero pivot. It now calls `ldlt!(F, A; check = false)`, and when that also fails it returns `lu(A)`, as the Hermitian `\` does (`cholmod.jl`). So `factorize` now succeeds wherever `\` does. A singular matrix still throws, now a `SingularException` from `lu`. The inferred return type of `factorize(::Hermitian{T,SparseMatrixCSC{T,Int}})` widens from `CHOLMOD.Factor{T,Int}` to `Union{CHOLMOD.Factor{T,Int}, UMFPACK.UmfpackLU{T,Int}}` (`UmfpackLU{Float64,Int}` for `Float32`). `factorize(::SparseMatrixCSC)` already returned a seven-member union that includes both. ๐Ÿค– Generated with [Claude Code](https://claude.com/claude-code) Co-authored-by: Claude Opus 5.5 (cherry picked from commit 6ab13eeb0bb965b23708da42aeed73848b2f8252) --- src/linalg.jl | 11 +++++------ test/solvers/issues.jl | 5 +++++ 2 files changed, 10 insertions(+), 6 deletions(-) diff --git a/src/linalg.jl b/src/linalg.jl index 03807a5b..b3987759 100644 --- a/src/linalg.jl +++ b/src/linalg.jl @@ -1504,12 +1504,11 @@ factorize(A::AbstractSparseMatrixCSC) = _sparse_factorize(_factorize_choice, A) function factorize(A::RealHermSymComplexHerm{<:Union{Float32,Float64},<:AbstractSparseMatrixCSC}) F = cholesky(A; check = false) - if LinearAlgebra.issuccess(F) - return F - else - ldlt!(F, A) - return F - end + LinearAlgebra.issuccess(F) && return F + ldlt!(F, A; check = false) + LinearAlgebra.issuccess(F) && return F + # an indefinite matrix may need pivoting, which only the LU does, as in `\` + return lu(A) end eigen(A::AbstractSparseMatrixCSC) = diff --git a/test/solvers/issues.jl b/test/solvers/issues.jl index 4d218bc3..51ca3b60 100644 --- a/test/solvers/issues.jl +++ b/test/solvers/issues.jl @@ -249,4 +249,9 @@ end @test ldiv!(lu(A), ComplexF32[1, 2]) โ‰ˆ ComplexF32.(Matrix(A) \ [1, 2]) end +@testset "factorize of an indefinite Hermitian matrix whose LDLt needs pivoting" begin + b = [1.0, 2.0] + @test factorize(sparse([0.0 1; 1 0])) \ b โ‰ˆ [0.0 1; 1 0] \ b +end + end # module From d4aca46f3a029325b24736579f3ee7c0dd94a1e8 Mon Sep 17 00:00:00 2001 From: "Viral B. Shah" Date: Mon, 26 Oct 2026 17:09:05 -0400 Subject: [PATCH 20/30] Convert scaled values before writing the destination of sparse `mul!`, `ldiv!` and `rdiv!` (#1121) MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Scaling a sparse matrix into a destination, by a number (`mul!(C, A, b)`, `mul!(C, b, A)`) or by a `Diagonal` (`ldiv!(C, D, A)`, and `D \ A` and `A / D` through it), gave `C` the pattern of `A` before converting the values to `eltype(C)`. A value that does not convert left `C` with the new pattern and uninitialized values. The in-place `Diagonal` forms (`rmul!`, `lmul!`, `rdiv!`, `ldiv!`, also on a sparse vector) stopped partway, with some entries scaled. ```julia using SparseArrays, LinearAlgebra A = sparse([2 0; 0 3]); D = Diagonal([2, 2]) C = spzeros(Int, 2, 2); try mul!(C, A, 0.5) catch e; println(e) end; @show nnz(C) C = spzeros(Int, 2, 2); try ldiv!(C, D, A) catch e; println(e) end; @show nnz(C) B = copy(A); try rdiv!(B, D) catch e; println(e) end; @show nonzeros(B) ``` **Before** ``` InexactError(:Int64, (Int64, 1.5)) nnz(C) = 2 InexactError(:Int64, (Int64, 1.5)) nnz(C) = 2 InexactError(:Int64, (Int64, 1.5)) nonzeros(B) = [1, 3] ``` **After** ``` InexactError(:Int64, (Int64, 1.5)) nnz(C) = 0 InexactError(:Int64, (Int64, 1.5)) nnz(C) = 0 InexactError(:Int64, (Int64, 1.5)) nonzeros(B) = [2, 3] ``` The kernels now check with `_storesafely`, as the `Diagonal` `mul!` already does, whether every value of the product's eltype converts to `eltype(C)`. When it does, nothing changes: the values are written in place, with no allocation. Otherwise they are written to a new vector of `eltype(C)`, and `C` takes that vector and the pattern of `A` only after every value has converted. Successful calls return the same results as before. Left out: in-place scaling by a number (`rmul!(A, b)`, `lmul!`, `rdiv!`, `ldiv!`) uses the dense kernels of Base and LinearAlgebra on the stored values, which stop partway just as they do for a dense array. An `op` that throws for a reason other than conversion still leaves `C` changed, as in the `Diagonal` `mul!`. ๐Ÿค– Generated with [Claude Code](https://claude.com/claude-code) Co-authored-by: Claude Opus 5.5 (cherry picked from commit 66a7b5cc929902031e4aade72b9c11ed6e1c0f7f) --- src/linalg.jl | 33 +++++++++++++++++++++++++-------- src/matmul.jl | 12 ++++++------ test/linalg.jl | 9 +++++++++ test/matmul.jl | 9 +++++++++ 4 files changed, 49 insertions(+), 14 deletions(-) diff --git a/src/linalg.jl b/src/linalg.jl index b3987759..4d59b8da 100644 --- a/src/linalg.jl +++ b/src/linalg.jl @@ -561,39 +561,43 @@ end # Scale the stored entries of `A` into `C`, which either is `A` or takes its pattern: # `C[i, j] = op(A[i, j], d[j])` columnwise, `C[i, j] = op(d[i], A[i, j])` rowwise. The # argument order of `op` is kept for non-commutative eltypes. The callers check that -# `length(d)` matches the scaled dimension and that `size(C) == size(A)`. +# `length(d)` matches the scaled dimension and that `size(C) == size(A)`. When a value might +# not convert to `eltype(C)`, the values go to a new vector that `C` takes once all are written. function _scalecols!(C::AbstractSparseMatrixCSC, A::AbstractSparseMatrixCSC, d::AbstractVector, op::F) where {F} - copyinds!(C, A) - Cnzval = getnzval(C) Anzval = getnzval(A) - resize!(Cnzval, length(Anzval)) + direct = _storesafely(eltype(C), promote_op(op, eltype(A), eltype(d))) + Cnzval = direct ? _takepattern!(C, A) : Vector{eltype(C)}(undef, nnz(A)) @inbounds for col in axes(A, 2) dcol = d[col] for p in nzrange(A, col) Cnzval[p] = op(Anzval[p], dcol) end end + direct || _storescaled!(C, A, Cnzval) return C end function _scalerows!(C::AbstractSparseMatrixCSC, A::AbstractSparseMatrixCSC, d::AbstractVector, op::F) where {F} - copyinds!(C, A) - Cnzval = getnzval(C) Anzval = getnzval(A) Arowval = getrowval(A) - resize!(Cnzval, length(Anzval)) + direct = _storesafely(eltype(C), promote_op(op, eltype(d), eltype(A))) + Cnzval = direct ? _takepattern!(C, A) : Vector{eltype(C)}(undef, nnz(A)) @inbounds for col in axes(A, 2), p in nzrange(A, col) Cnzval[p] = op(d[Arowval[p]], Anzval[p]) end + direct || _storescaled!(C, A, Cnzval) return C end function _scalerows!(x::AbstractSparseVector, d::AbstractVector, op::F) where {F} xnzval = nonzeros(x) xnzind = nonzeroinds(x) + direct = _storesafely(eltype(x), promote_op(op, eltype(d), eltype(x))) + vals = direct ? xnzval : similar(xnzval) @inbounds for p in eachindex(xnzind, xnzval) - xnzval[p] = op(d[xnzind[p]], xnzval[p]) + vals[p] = op(d[xnzind[p]], xnzval[p]) end + direct || copyto!(xnzval, vals) return x end @@ -1350,6 +1354,19 @@ function copyinds!(C::AbstractSparseMatrixCSC, A::AbstractSparseMatrixCSC; copy_ end end +# Give `C` the pattern of `A`, and return the values of `C`, of the length of those of `A`. +function _takepattern!(C::AbstractSparseMatrixCSC, A::AbstractSparseMatrixCSC) + copyinds!(C, A) + return resize!(getnzval(C), length(getnzval(A))) +end +# Store `nzval` in `C` with the pattern of `A`, converting it first so that a value that does +# not convert to `eltype(C)` throws with `C` unchanged. +function _storescaled!(C::AbstractSparseMatrixCSC, A::AbstractSparseMatrixCSC, nzval::AbstractVector) + vals = convert(Vector{eltype(C)}, nzval) + copyto!(_takepattern!(C, A), vals) + return C +end + """ rowcheck_index(A::AbstractSparseMatrixCSC, row::Integer, col::Integer) diff --git a/src/matmul.jl b/src/matmul.jl index 50de6da2..52d8a371 100644 --- a/src/matmul.jl +++ b/src/matmul.jl @@ -1022,18 +1022,18 @@ _lmulkeepszeros(b::Number, A) = (z = _implicit_zero(A); z isa UnknownZero || _is function mul!(C::AbstractSparseMatrixCSC, A::AbstractSparseMatrixCSC, b::Number) size(A)==size(C) || throw(DimensionMismatch("A has size $(size(A)) but C has size $(size(C))")) (_is_fixed(C) || !_rmulkeepszeros(A, b)) && return broadcast!(*, C, A, b) - copyinds!(C, A) - resize!(getnzval(C), length(getnzval(A))) - mul!(getnzval(C), getnzval(A), b) + _storesafely(eltype(C), promote_op(*, eltype(A), typeof(b))) || + return _storescaled!(C, A, view(getnzval(A), 1:nnz(A)) * b) + mul!(_takepattern!(C, A), getnzval(A), b) C end function mul!(C::AbstractSparseMatrixCSC, b::Number, A::AbstractSparseMatrixCSC) size(A)==size(C) || throw(DimensionMismatch("A has size $(size(A)) but C has size $(size(C))")) (_is_fixed(C) || !_lmulkeepszeros(b, A)) && return broadcast!(*, C, b, A) - copyinds!(C, A) - resize!(getnzval(C), length(getnzval(A))) - mul!(getnzval(C), b, getnzval(A)) + _storesafely(eltype(C), promote_op(*, typeof(b), eltype(A))) || + return _storescaled!(C, A, b * view(getnzval(A), 1:nnz(A))) + mul!(_takepattern!(C, A), b, getnzval(A)) C end diff --git a/test/linalg.jl b/test/linalg.jl index d206f5eb..3c3adc0e 100644 --- a/test/linalg.jl +++ b/test/linalg.jl @@ -351,6 +351,15 @@ end C = spzeros(ComplexF64, 3, 3) @test_throws DimensionMismatch ldiv!(C, Dl, A) @test nnz(C) == 0 + @static if COMPREHENSIVE + # a quotient that does not convert to the eltype of the destination throws with it unchanged + A0 = sparse([2 0; 0 3]); x0 = sparsevec([1, 2], [2, 3]); D2 = Diagonal([2, 2]) + for (f!, X0) in ((C -> ldiv!(C, D2, A0), sparse([1 0; 0 0])), (B -> rdiv!(B, D2), A0), (y -> ldiv!(D2, y), x0)) + X = copy(X0) + @test_throws InexactError f!(X) + @test mismatch(X, X0) === nothing + end + end end end diff --git a/test/matmul.jl b/test/matmul.jl index 4d458f9d..dee8035b 100644 --- a/test/matmul.jl +++ b/test/matmul.jl @@ -573,6 +573,15 @@ end @test mismatch(mul!(sC, sA, 0.5), 0.5 * dA) === nothing @test mismatch(lmul!(0.5, copy(sA)), 0.5 * dA) === nothing @test mul!(sC, 0.5, sA) == mul!(sC, sA, 0.5) + @static if COMPREHENSIVE + # a product that does not convert to the eltype of `C` throws before `C` takes the pattern of `A` + A0 = sparse([1 0; 0 3]); C0 = sparse([2 0; 0 0]) + for f! in (C -> mul!(C, A0, 0.5), C -> mul!(C, 0.5, A0)) + C = copy(C0) + @test_throws InexactError f!(C) + @test mismatch(C, C0) === nothing + end + end @testset "inverse scaling with mul!" begin bi = inv.(b) From 84c3606e5905cfbce7abbab03f92d0272fc51047 Mon Sep 17 00:00:00 2001 From: "Viral B. Shah" Date: Mon, 26 Oct 2026 18:08:09 -0400 Subject: [PATCH 21/30] Return `Inf` from sparse `cond` of a singular matrix instead of throwing (#1126) MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit `cond(A, 1)` and `cond(A, Inf)` of a singular sparse matrix throw `SingularException`. Dense `cond` returns `Inf` for these. ```julia using SparseArrays, LinearAlgebra for A in ([1.0 2; 3 6], [1.0 0; 2 0]), p in (1, Inf) println("cond(sparse($A), $p) = ", try cond(sparse(A), p) catch e; sprint(showerror, e) end) end ``` **Before** ``` cond(sparse([1.0 2.0; 3.0 6.0]), 1) = SingularException: matrix is singular; factorization failed. Zero pivot found at index 0 cond(sparse([1.0 2.0; 3.0 6.0]), Inf) = SingularException: matrix is singular; factorization failed. Zero pivot found at index 0 cond(sparse([1.0 0.0; 2.0 0.0]), 1) = SingularException: matrix is singular; factorization failed. Zero pivot found at index 2 cond(sparse([1.0 0.0; 2.0 0.0]), Inf) = SingularException: matrix is singular; factorization failed. Zero pivot found at index 2 ``` **After** ``` cond(sparse([1.0 2.0; 3.0 6.0]), 1) = Inf cond(sparse([1.0 2.0; 3.0 6.0]), Inf) = Inf cond(sparse([1.0 0.0; 2.0 0.0]), 1) = Inf cond(sparse([1.0 0.0; 2.0 0.0]), Inf) = Inf ``` `cond` estimates the norm of the inverse with `opnormestinv`, which solves with `factorize(A)`. A singular matrix makes `lu` throw. A triangular or diagonal matrix gets no factorization object, so its solve throws. `issuccess` on `lu(A; check = false)` would only cover the `lu` branch, so `cond` catches `SingularException` instead, as dense `cond` does. It then returns `convert(float(real(eltype(A))), Inf)`, the same value and type as dense. A `Float32` matrix gives `Inf32`. The inferred return types are unchanged. A singular Hermitian matrix still throws `ZeroPivotException` from the LDLt in `factorize`. That exception is not caught, because on `main` a nonsingular indefinite matrix such as `[0.0 1; 1 0]` throws it too. #1116 makes that path fall back to `lu`, which throws `SingularException`. With both PRs a singular Hermitian matrix also gives `Inf`; I checked this on #1116's branch with this change applied. This PR does not depend on #1116. ๐Ÿค– Generated with [Claude Code](https://claude.com/claude-code) Co-authored-by: Claude Opus 5.5 (cherry picked from commit 4c878618d95b6b6994e7f2529bcc1c80bb5a3db9) --- src/linalg.jl | 18 ++++++++++-------- test/solvers/issues.jl | 6 ++++++ 2 files changed, 16 insertions(+), 8 deletions(-) diff --git a/src/linalg.jl b/src/linalg.jl index 4d59b8da..446f071d 100644 --- a/src/linalg.jl +++ b/src/linalg.jl @@ -990,14 +990,16 @@ small enough to store dense. function cond(A::AbstractSparseMatrixCSC, p::Real=2) # as for a dense matrix; the estimator of the inverse's norm needs at least one column (p == 1 || p == Inf) && checksquare(A) == 0 && return zero(real(float(eltype(A)))) - if p == 1 - normAinv = opnormestinv(A) - normA = opnorm(A, 1) - return normA * normAinv - elseif p == Inf - normAinv = opnormestinv(copy(A')) - normA = opnorm(A, Inf) - return normA * normAinv + if p == 1 || p == Inf + # a singular matrix has an infinite condition number, as for a dense matrix; a + # triangular or diagonal `factorize(A)` reports it only in the solve + normAinv = try + p == 1 ? opnormestinv(A) : opnormestinv(copy(A')) + catch e + e isa SingularException || rethrow() + return convert(float(real(eltype(A))), Inf) + end + return opnorm(A, p) * normAinv elseif p == 2 throw(ArgumentError("only 1- and Inf-norm condition numbers are implemented for sparse matrices, for 2-norm try cond(Array(A), 2) instead")) else diff --git a/test/solvers/issues.jl b/test/solvers/issues.jl index 51ca3b60..60e094f9 100644 --- a/test/solvers/issues.jl +++ b/test/solvers/issues.jl @@ -254,4 +254,10 @@ end @test factorize(sparse([0.0 1; 1 0])) \ b โ‰ˆ [0.0 1; 1 0] \ b end +@testset "cond of a singular matrix is Inf" begin + @test cond(sparse([1.0 2; 3 6]), 1) === Inf + # the triangular solve, not the factorization, finds this one singular + @test cond(sparse([1.0 0; 2 0]), Inf) === Inf +end + end # module From e97706c5647abd71c7441340f0721083cb4fc242 Mon Sep 17 00:00:00 2001 From: "Viral B. Shah" Date: Mon, 26 Oct 2026 18:51:02 -0400 Subject: [PATCH 22/30] Handle sparse matrices of 2^63 or more entries in reductions, `sprand` and linear indexing (#1128) MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit `widelength`, the length that the sparse kernels compare with `nnz`, multiplied the dimensions as `Int64`. For a matrix of 2^63 entries or more it wrapped, so the reductions and `sprand` gave wrong results or threw, and linear indexing failed with a misleading `BoundsError`. **Before** ```julia julia> B = spzeros(2^62, 4); # 2^64 entries, none stored julia> any(iszero, B), all(isone, B), count(iszero, B; init = big(0)) (false, true, 0) julia> sum(x -> 1.0, B) ERROR: ArgumentError: reducing over an empty collection is not allowed; consider supplying `init` to the reducer julia> findmax(B) ERROR: ArgumentError: array slices must be non-empty julia> nnz(sprand(Xoshiro(1), 2^62, 4, 1e-17)) 0 julia> B[[1]] = [1.0] ERROR: BoundsError: attempt to access 4611686018427387904ร—4 SparseMatrixCSC{Float64, Int64} with 0 stored entries at index [1] ``` **After** ```julia julia> any(iszero, B), all(isone, B), count(iszero, B; init = big(0)) (true, false, 18446744073709551616) julia> sum(x -> 1.0, B) 1.8446744073709552e19 julia> findmax(B) (0.0, CartesianIndex(1, 1)) julia> nnz(sprand(Xoshiro(1), 2^62, 4, 1e-17)) 169 julia> B[[1]] = [1.0] ERROR: OverflowError: a 4611686018427387904ร—4 sparse matrix has more entries than an `Int` can index linearly; index it by row and column ``` `widelength` is now `prod(widen, size(x))`, an `Int128` on a 64-bit build and `Int64`, as before, on a 32-bit one. `count`, the two-array `mapreduce` and `_inorder` use it in place of `length` or `prod(size(A))`, and the counters of unstored runs in the ordered reductions are widened the same way. An `Int128` count would widen `z * n` for `Complex{Int}` and reach a different `^` method for an eltype that defines `^` for an `Int` power only, so `_addzeros` and `_powzeros` split it into `Int` factors and powers. For a matrix that fits an `Int` they make the same calls as before, and `prod` stays type-stable. When `m*n` overflows, `sprand` draws the rows column by column, which keeps the distribution; every size that fits draws the same entries as before. Linear indexing into such a matrix is not supported. `getindex` and `setindex!` with a linear index, range or array, including the fixed-pattern `setindex!`, throw an `OverflowError` before writing anything; indexing by row and column works. This touches dispatch: scalar `A[i] = v` gets a method on `AbstractSparseMatrixCSC`, as Base's fallback checks the index against the wrapping `length`. On a 10โดร—10โด matrix of density 10โปยณ, `main` against this branch: | | `main` | branch | |---|---|---| | `sum` / `prod` / `count` (Float64) | 10.1 / 13.7 / 11.6 ฮผs | 10.1 / 13.7 / 11.6 ฮผs | | ordered `mapreduce` (Float64 / Int) | 101.6 / 80.5 ฮผs | 101.0 / 77.1 ฮผs | | `mapreduce(*, -, A, A)` (Float64 / Int) | 420.8 / 135.2 ฮผs | 423.5 / 133.1 ฮผs | | `findmax` (Float64) | 221.4 ฮผs | 221.5 ฮผs | | `A[I] = x`, `A[I]` (10โด indices) | 448 / 152 ฮผs | 450 / 152 ฮผs | | `A[i] = v`, `A[i]` | 8.2 / 6.7 ns | 9.0 / 7.0 ns | | `sprand` | 580 ฮผs | 580 ฮผs | The result types and values of 18 reductions over 10 eltypes, and the entries `sprand` draws in each of its forms, match `main`. Left out: Base's own `length(A)` still wraps, and so do the checks that use it or `m*n`, such as `isempty`, `Base.checkbounds` and the dims reductions' `nnz(A) == m*n`. ๐Ÿค– Generated with [Claude Code](https://claude.com/claude-code) --------- Co-authored-by: Claude Opus 5.5 (cherry picked from commit 2bef1fd8b9e387500af770c0c1a31fec9ec4d56c) --- src/SparseArrays.jl | 6 +++++- src/constructors.jl | 25 ++++++++++++++++++++++- src/indexing.jl | 17 ++++++++++++++-- src/reductions.jl | 40 ++++++++++++++++++++++++++----------- src/sparsematrix.jl | 5 +++-- test/issues.jl | 14 +++++++++++++ test/trim/src/SparseTrim.jl | 7 +++++++ 7 files changed, 96 insertions(+), 18 deletions(-) diff --git a/src/SparseArrays.jl b/src/SparseArrays.jl index 04c8cf3b..e101d332 100644 --- a/src/SparseArrays.jl +++ b/src/SparseArrays.jl @@ -380,7 +380,11 @@ julia> collect(SparseArrays.iternz(x)) """ function iternz end -widelength(x::AbstractSparseArray) = prod(Int64.(size(x))) +widelength(x::AbstractSparseArray) = prod(widen, size(x)) +# Linear indexing compares an index with `length(A)`, which wraps for 2^63 entries or more +_checklinearindexing(A::AbstractSparseMatrixCSC) = widelength(A) <= typemax(Int) || + throw(OverflowError(string("a ", size(A, 1), "ร—", size(A, 2), " sparse matrix has more ", + "entries than an `Int` can index linearly; index it by row and column"))) const _restore_scalar_indexing = Expr[] diff --git a/src/constructors.jl b/src/constructors.jl index 1625213a..bf4452bc 100644 --- a/src/constructors.jl +++ b/src/constructors.jl @@ -714,20 +714,43 @@ julia> sprand(Float64, 3, 0.75) """ function sprand(r::AbstractRNG, m::Integer, n::Integer, density::AbstractFloat, rfn::Function, ::Type{T}) where T m, n = Int(m), Int(n) + if _sprandoverflows(m, n, density) + colptr, rowval, V = _sprand_bycolumn(r, m, n, density, rfn) + return SparseMatrixCSC(m, n, colptr, rowval, convert(Vector{T}, V)) + end I, V = _sprandvalues(r, m, n, density, rfn) return sparse_sortedlinearindices!(I, convert(Vector{T}, V), m, n) end function sprand(r::AbstractRNG, m::Integer, n::Integer, density::AbstractFloat, rfn::Function) m, n = Int(m), Int(n) + if _sprandoverflows(m, n, density) + colptr, rowval, V = _sprand_bycolumn(r, m, n, density, rfn) + return SparseMatrixCSC(m, n, colptr, rowval, convert(Vector{eltype(V)}, V)) + end I, V = _sprandvalues(r, m, n, density, rfn) return sparse_sortedlinearindices!(I, convert(Vector{eltype(V)}, V), m, n) end -function _sprandvalues(r::AbstractRNG, m::Int, n::Int, density::AbstractFloat, rfn::F) where {F} +# the arguments are checked first, since negative dimensions can overflow too +function _sprandoverflows(m::Int, n::Int, density::AbstractFloat) (m < 0 || n < 0) && throw(ArgumentError("invalid Array dimensions")) 0 <= density <= 1 || throw(ArgumentError("$density not in [0,1]")) + return Base.mul_with_overflow(m, n)[2] +end +function _sprandvalues(r::AbstractRNG, m::Int, n::Int, density::AbstractFloat, rfn::F) where {F} I = randsubseq(r, 1:(m*n), density) return I, rfn(r, length(I)) end +# an `Int` cannot index the `m*n` entries, so the stored rows are drawn column by column +function _sprand_bycolumn(r::AbstractRNG, m::Int, n::Int, density::AbstractFloat, rfn::F) where {F} + colptr = Vector{Int}(undef, n + 1) + colptr[1] = 1 + rowval, rows = Int[], Int[] + for col in 1:n + append!(rowval, randsubseq!(r, rows, 1:m, density)) + colptr[col + 1] = length(rowval) + 1 + end + return colptr, rowval, rfn(r, length(rowval)) +end sprand(m::Integer, n::Integer, density::AbstractFloat, rfn::Function) = sprand(default_rng(), m, n, density, (r, i) -> rfn(i)) diff --git a/src/indexing.jl b/src/indexing.jl index 0838852c..c0cda6d3 100644 --- a/src/indexing.jl +++ b/src/indexing.jl @@ -35,8 +35,10 @@ _integer_indices(I::AbstractArray{CartesianIndex{1}}) = map(first โˆ˜ Tuple, I) # Index types no method below matches (`CartesianIndex`, custom indices) are lowered and # dispatched again; indices that are already lowered fall back to Base. -@propagate_inbounds getindex(A::AbstractSparseMatrixCSC{Tv,Ti}, I...) where {Tv,Ti} = - _getindex_lowered(A, I, _lower_indices(A, I...)) +@propagate_inbounds function getindex(A::AbstractSparseMatrixCSC{Tv,Ti}, I...) where {Tv,Ti} + I isa Tuple{Integer} && _checklinearindexing(A) + return _getindex_lowered(A, I, _lower_indices(A, I...)) +end @propagate_inbounds getindex(x::AbstractSparseVector{Tv,Ti}, I...) where {Tv,Ti} = _getindex_lowered(x, I, _lower_indices(x, I...)) @propagate_inbounds _getindex_lowered(A, I, L) = A[L...] @@ -503,6 +505,7 @@ end function getindex(A::AbstractSparseMatrixCSC{Tv,Ti}, I::AbstractArray) where {Tv,Ti} require_one_based_indexing(A, I) + eltype(I) <: Integer && _checklinearindexing(A) ndims(I) == 2 || return _getindex_dense(A, I) szA = size(A) nA = szA[1]*szA[2] @@ -596,6 +599,13 @@ _checkscalarindex(::Integer) = nothing # dispatch helper for #29034 @RCI setindex!(A::AbstractSparseMatrixCSC, _v, _i::Integer, _j::Integer) = _setindex_scalar!(A, _v, _i, _j) +@RCI function setindex!(A::AbstractSparseMatrixCSC, _v, _i::Integer) + _checkscalarindex(_i) + _checklinearindexing(A) + checkbounds(A, _i) + col, row = divrem(Int(_i) - 1, size(A, 1)) + return _setindex_scalar!(A, _v, row + 1, col + 1) +end function _setindex_scalar!(A::AbstractSparseMatrixCSC{Tv,Ti}, _v, _i::Integer, _j::Integer) where {Tv,Ti<:Integer} _checkscalarindex(_i); _checkscalarindex(_j) @@ -1124,6 +1134,7 @@ end function setindex!(A::AbstractSparseMatrixCSC, x::AbstractArray, Ix::AbstractVector{<:Integer}) require_one_based_indexing(A, x, Ix) + _checklinearindexing(A) # to_indices does not compare a Bool mask's length with A's Ix isa AbstractVector{Bool} && checkbounds(A, Ix) (I,) = Base.ensure_indexable(to_indices(A, (Ix,))) @@ -1564,6 +1575,7 @@ getindex(A::AbstractSparseMatrixCSC, ::Colon) = A[1:end] function getindex(A::AbstractSparseMatrixCSC{Tv}, I::AbstractUnitRange) where Tv require_one_based_indexing(A, I) + _checklinearindexing(A) checkbounds(A, I) szA = size(A) nA = szA[1]*szA[2] @@ -1606,6 +1618,7 @@ end function getindex(A::AbstractSparseMatrixCSC{Tv,Ti}, I::AbstractVector) where {Tv,Ti} I isa Base.LogicalIndex && return _logical_index(A, I.mask) require_one_based_indexing(A, I) + eltype(I) <: Integer && _checklinearindexing(A) @boundscheck checkbounds(A, I) szA = size(A) nA = szA[1]*szA[2] diff --git a/src/reductions.jl b/src/reductions.jl index 1db9702c..8f0d86a6 100644 --- a/src/reductions.jl +++ b/src/reductions.jl @@ -95,13 +95,13 @@ function _mapreduce_ordered(f::F, op::G, A::SparseMatrixCSCOrColumnSubset{T}) wh v, dropout = _foldstored(f, op, z, nzvalview(A)) dropout && return v v = Base._InitialValue() - unstored = Int64(0) # spans columns, so it can exceed typemax(Int32) + unstored = widen(0) # spans columns, so it can exceed typemax(Int) @inbounds for col in axes(A, 2) prev = 0 for k in getnzrange(A, col) row = Int(rowval[k]) v = rf(v, _afterunstored(op, z, idem, unstored + (row - prev - 1), f(nzval[k]))) - unstored, prev = Int64(0), row + unstored, prev = widen(0), row end unstored += m - prev end @@ -187,14 +187,30 @@ end _mapreducezeros(f::F, op::Union{typeof(Base.add_sum),typeof(+)}, ::Type{T}, nzeros::Integer, v0) where {F,T} = nzeros == 0 ? Base.reduce_first(op, v0) : op(_addzeros(_unstored(op, f, T, v0), nzeros), v0) _mapreducezeros(f::F, op::Union{typeof(Base.mul_prod),typeof(*)},::Type{T}, nzeros::Integer, v0) where {F,T} = - nzeros == 0 ? Base.reduce_first(op, v0) : op(_unstored(op, f, T, v0)^nzeros, v0) + nzeros == 0 ? Base.reduce_first(op, v0) : op(_powzeros(_unstored(op, f, T, v0), nzeros), v0) +# `z` to the power of `nzeros`. An eltype may define `^` for an `Int` power only, so a count +# from `widelength` is split into `Int` powers, as `_addzeros` splits it into factors. +_powzeros(z, nzeros) = z^nzeros +function _powzeros(z, nzeros::Int128) + nzeros <= typemax(Int) && return z^Int(nzeros) + q, r = divrem(nzeros, typemax(Int)) + return (z^typemax(Int))^Int(q) * z^Int(r) +end _unstored(op, f, ::Type{T}, v0) where T = _widen(Base.reduce_first(op, f(zero(T))), v0) _widen(z, v0) = _widen(z, typeof(v0)) _widen(z, ::Type{S}) where S = convert(promote_type(typeof(z), S), z) # `nzeros` copies of `z` added up, in the type `+` keeps for them: a small integer wraps -# rather than widening to the count's `Int`, as when dense input is summed with `+` +# rather than widening to the count's `Int`, as when dense input is summed with `+`. A count +# from `widelength` is an `Int128`, which would widen `Complex{Int}` and the like, so it is +# split into `Int` factors. _addzeros(z, nzeros) = z * nzeros +function _addzeros(z, nzeros::Int128) + nzeros <= typemax(Int) && return z * Int(nzeros) + q, r = divrem(nzeros - 1, typemax(Int)) + return z * Int(q) * typemax(Int) + z * Int(r + 1) +end _addzeros(z::Base.BitInteger, nzeros) = z * (nzeros % typeof(z)) +_addzeros(z::Base.BitInteger, nzeros::Int128) = z * (nzeros % typeof(z)) _mapreducezeros(f::F, op::Union{typeof(min),typeof(max)}, ::Type{T}, nzeros::Integer, v0) where {F,T} = nzeros == 0 ? v0 : op(v0, f(zero(T))) _mapreducezeros(f::Base.ExtremaMap, op::typeof(Base._extrema_rf), ::Type{T}, nzeros::Integer, v0) where {T} = @@ -220,7 +236,7 @@ function Base._mapreduce(f::F, op::Union{typeof(Base.mul_prod),typeof(*)}, ::Bas elseif !_commutes(f, op, T) _mapreduce_ordered(f, op, A) else - v = Base.reduce_first(op, f(zero(T)))^nzeros + v = _powzeros(Base.reduce_first(op, f(zero(T))), nzeros) (nnzA == 0 || (_iszero(v) && Base.promote_op(f, T) <: Integer)) && return v r = v*Base._mapreduce(f, op, nzvalview(A)) # A zero ends the product unless a stored `Inf` or `NaN` makes it `NaN`. The product @@ -266,7 +282,7 @@ end # A reduction along a dimension, an empty array and a shape mismatch are left to Base's method. for T in (:SparseVectorOrView, :SparseMatrixCSCOrColumnSubset) @eval function Base.mapreduce(f::F, op::G, A::$T, B::$T; dims=:, init=Base._InitialValue(), kw...) where {F,G} - (dims === (:) && isempty(kw) && size(A) == size(B) && !isempty(A)) || + (dims === (:) && isempty(kw) && size(A) == size(B) && !iszero(widelength(A))) || return reduce(op, map(f, A, B); dims, init, kw...) return _commutes(f, op, eltype(A), eltype(B)) ? _mapreduce_unordered(f, op, init, A, B) : _mapreduce_ordered(f, op, init, A, B) @@ -280,7 +296,7 @@ function _mapreduce_ordered(f::F, op::G, init, A::AbstractArray{Ta}, B::Abstract Ai, Av, Bi, Bv = getrowval(A), getnzval(A), getrowval(B), getnzval(B) m = size(A, 1) v = init - unstored = Int64(0) # positions that neither array stores, since the last one folded + unstored = widen(0) # positions that neither array stores, since the last one folded @inbounds for col in axes(A, 2) ra, rb = getnzrange(A, col), getnzrange(B, col) ia, ib = first(ra), first(rb) @@ -290,18 +306,18 @@ function _mapreduce_ordered(f::F, op::G, init, A::AbstractArray{Ta}, B::Abstract row = Int(min(Ai[ia], Bi[ib])) x, ia, ib = _merge_step(f, Ai, Av, ia, Bi, Bv, ib) v = rf(_foldunstored(f, rf, v, Ta, Tb, unstored + (row - prev - 1)), x) - unstored, prev = Int64(0), row + unstored, prev = widen(0), row end while ia <= ea row = Int(Ai[ia]) v = rf(_foldunstored(f, rf, v, Ta, Tb, unstored + (row - prev - 1)), f(Av[ia], zero(Tb))) - unstored, prev = Int64(0), row + unstored, prev = widen(0), row ia += 1 end while ib <= eb row = Int(Bi[ib]) v = rf(_foldunstored(f, rf, v, Ta, Tb, unstored + (row - prev - 1)), f(zero(Ta), Bv[ib])) - unstored, prev = Int64(0), row + unstored, prev = widen(0), row ib += 1 end unstored += m - prev @@ -365,7 +381,7 @@ function _mapreduce_unordered(f::F, op::G, init, A::AbstractArray{Ta}, B::Abstra ib == eb && (v = rf(v, f(zero(Ta), Bv[ib]))) end w isa Base._InitialValue || (v = v isa Base._InitialValue ? w : op(v, w)) - nzeros = prod(Int64.(size(A))) - stored + nzeros = prod(widen, size(A)) - stored if v isa Base._InitialValue # nothing stored and no `init`: seed with one unstored position v = rf(v, f(zero(Ta), zero(Tb))) nzeros -= 1 @@ -627,7 +643,7 @@ end # The kernels above fold each slice's stored entries first and its unstored ones last. That # matches element order only when `op` commutes or every entry is stored; otherwise the # kernels below fold the unstored entries where they lie. -_inorder(f, op, A::AbstractArray{T}) where {T} = _commutes(f, op, T) || nnz(A) == length(A) +_inorder(f, op, A::AbstractArray{T}) where {T} = _commutes(f, op, T) || nnz(A) == widelength(A) function _mapreducerows_ordered!(f::F, op::G, R::AbstractArray, A::SparseMatrixCSCOrColumnSubset{T}) where {F,G,T} require_one_based_indexing(A, R) rowval, nzval = getrowval(A), getnzval(A) diff --git a/src/sparsematrix.jl b/src/sparsematrix.jl index 12787a2e..010e0dee 100644 --- a/src/sparsematrix.jl +++ b/src/sparsematrix.jl @@ -319,7 +319,8 @@ function _storedinds(S::SparseMatrixCSCColumnSubset) end return inds end -widelength(S::SparseMatrixCSCColumnSubset) = prod(Int64.(size(S))) +widelength(S::SparseMatrixCSCColumnSubset) = prod(widen, size(S)) +widelength(x::Union{SparseColumnView,SparseVectorView}) = length(x) # `A` itself where it exposes compressed storage through `nzrange`, `rowvals` and `nonzeros`, # otherwise its sparse copy _compressed(A) = A @@ -353,7 +354,7 @@ nnz(S::SparseMatrixCSCView) = length(_storedinds(S)) nnz1(S) = @inbounds sum(length.(getnzrange.(Ref(S), axes(S, 2)))) function Base._simple_count(pred, S::SparseMatrixCSCOrColumnSubset, init::T) where T - init + T(count(pred, nzvalview(S)) + pred(zero(eltype(S)))*(prod(size(S)) - nnz(S))) + init + T(count(pred, nzvalview(S)) + pred(zero(eltype(S)))*(widelength(S) - nnz(S))) end # the entries of a lazy adjoint or transpose are those of its parent, each under the wrapper Base._simple_count(pred, S::AdjOrTrans{<:Any,<:SparseMatrixCSCOrColumnSubset}, init) = diff --git a/test/issues.jl b/test/issues.jl index eb1972d7..647e16e6 100644 --- a/test/issues.jl +++ b/test/issues.jl @@ -968,6 +968,20 @@ end @inferred hcat(fixture(Float64, 2, 1), Matrix(1.0I, 2, 2)) end +@testset "a matrix of 2^64 entries does not overflow its length" begin + B = spzeros(2^62, 4) + @test any(iszero, B) && !all(isone, B) + @test sum(x -> 1.0, B) == 2.0^64 && prod(x -> 2.0, B) == Inf + @test count(iszero, B; init = big(0)) == big(2)^64 + @test findmax(B) == (0.0, CartesianIndex(1, 1)) + @test mapreduce(*, +, B, B) == 0.0 + @test_throws OverflowError B[[1]] + @test_throws OverflowError B[[1]] = [1.0] + @test_throws OverflowError B[1] = 1.0 + @test nnz(B) == 0 + @test nnz(sprand(Xoshiro(1), 2^62, 4, 1e-17)) > 0 +end + if isdefined(Docs, :undocumented_names) # new in Julia 1.11 @testset "docstrings (issue julia#52725)" begin @test isempty(Docs.undocumented_names(SparseArrays)) diff --git a/test/trim/src/SparseTrim.jl b/test/trim/src/SparseTrim.jl index 7d305668..18dd884f 100644 --- a/test/trim/src/SparseTrim.jl +++ b/test/trim/src/SparseTrim.jl @@ -91,6 +91,11 @@ function indexing() B[:, 2] .= 1.0 bad += check(B[:, 2] == [1.0, 1.0, 1.0]) bad += check(Array(view(A, :, 3)) == [0.0, 0.0, 5.0]) + bad += check(A[3] == 4.0 && A[[3, 5]] == sparsevec([1, 2], [4.0, 3.0], 2) && A[2:3] == sparsevec([2], [4.0], 2)) + C = copy(A) + C[2] = 1.0 + C[[4]] = [7.0] + bad += check(C[2, 1] == 1.0 && C[1, 2] == 7.0 && nnz(C) == 6) return bad end @@ -113,6 +118,8 @@ function algebra() bad += check(A - B == sparse([1, 2, 3, 3, 1], [1, 2, 1, 3, 3], [2.0, 3.0, 5.0, 5.0, -1.0], 3, 3)) bad += check(A * B == sparse([1, 3, 3], [3, 1, 3], [2.0, -5.0, 4.0], 3, 3)) bad += check(Matrix(A * D) == [4.0 0.0 0.0; 0.0 9.0 0.0; 28.0 0.0 25.0]) + bad += check(mapreduce(*, +, A, A) == 54.0 && prod(x -> x + 1.0, A) == 360.0) + bad += check(mapreduce(identity, (a, b) -> a + b, A) == 14.0) bad += check(sum(A) == 14.0 && sum(A; dims = 1) == [6.0 3.0 5.0]) bad += check(sum(A; dims = 2) == reshape([2.0, 3.0, 9.0], 3, 1)) bad += check(maximum(A) == 5.0 && minimum(B) == -1.0 && prod(A) == 0.0) From 61d114683e33aafc82f6ff9f8f4b5f4c13535c74 Mon Sep 17 00:00:00 2001 From: "Viral B. Shah" Date: Mon, 26 Oct 2026 20:10:56 -0400 Subject: [PATCH 23/30] Document that sparse `cond` of a singular matrix returns `Inf` (#1130) MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit #1126 made sparse `cond` return `Inf` for a singular matrix, but its docstring does not say so, nor that rounding can hide the singularity. Raised in the post-merge review of #1126. The material is the behaviour on `main`: ```julia julia> cond(sparse([1.0 2; 3 6]), 1) Inf julia> A = [1.0 2 3; 4 5 6; 7 8 9]; # singular julia> cond(sparse(A), 1), cond(A, 1) (1.0808639105689187e17, Inf) ``` The UMFPACK LU of the second matrix ends with a pivot that rounding leaves nonzero, while LAPACK's partial pivoting meets an exact zero, so the docstring says that a finite but very large result can still mean a singular matrix, and that sparse and dense can disagree on which case they report. Written by Claude Code. ๐Ÿค– Generated with [Claude Code](https://claude.com/claude-code) Co-authored-by: Claude Opus 5.5 (cherry picked from commit 5b98f3c435b465f281b69dcdf0f79c7e6bcd7663) --- src/linalg.jl | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/src/linalg.jl b/src/linalg.jl index 446f071d..ffb19d1b 100644 --- a/src/linalg.jl +++ b/src/linalg.jl @@ -986,6 +986,11 @@ estimator of Higham and Tisseur, from a factorization of `A` by `factorize` and forming the inverse, so the result can fall below the exact condition number. `p = 2`, the default, throws an `ArgumentError`, so `cond(A)` does too: use `cond(Array(A))` for a matrix small enough to store dense. + +A singular matrix gives `Inf` when its factorization meets an exactly zero pivot. Rounding +can hide that zero, and the result is then finite but very large. The sparse LU pivots in +another order than the dense one, so `cond(A, p)` can be finite where `cond(Array(A), p)` +is `Inf`. """ function cond(A::AbstractSparseMatrixCSC, p::Real=2) # as for a dense matrix; the estimator of the inverse's norm needs at least one column From 629baa6d6916a9aee748d6014b67aff44a66aa61 Mon Sep 17 00:00:00 2001 From: "Viral B. Shah" Date: Mon, 26 Oct 2026 20:22:59 -0400 Subject: [PATCH 24/30] Return the same float type from sparse `cond` for singular and nonsingular matrices (#1131) MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Sparse `cond` of a `Float32` or `ComplexF32` matrix returned `Float64` when the matrix went through the LU and was nonsingular, but `Float32` when it was singular (since #1126) or triangular, diagonal or positive definite. Dense `cond` returns `Float32` in every case. Raised in the post-merge review of #1126. **Before** ```julia julia> cond(sparse(Float32[1 2; 3 4]), 1) 20.99999999999999 julia> cond(sparse(Float32[1 2; 3 6]), 1) Inf32 julia> Base.return_types(cond, (SparseMatrixCSC{Float32,Int}, Int)) 1-element Vector{Any}: Union{Float32, Float64} ``` **After** ```julia julia> cond(sparse(Float32[1 2; 3 4]), 1) 21.0f0 julia> cond(sparse(Float32[1 2; 3 6]), 1) Inf32 julia> Base.return_types(cond, (SparseMatrixCSC{Float32,Int}, Int)) 1-element Vector{Any}: Float32 ``` UMFPACK works in double precision, so the estimate of the inverse's norm is a `Float64` for any eltype that takes the LU (a nonsymmetric matrix, or an indefinite Hermitian one where `factorize` falls back to `lu`), and `opnorm(A, p) * normAinv` promoted to it. `cond` now computes `T = float(real(eltype(A)))` once and returns `zero(T)`, `convert(T, Inf)` or `convert(T, opnorm(A, p) * normAinv)`. A sweep against dense `cond` over `Int`, `Float32`, `Float64`, `ComplexF32`, `ComplexF64`, `Rational{Int}`, `BigFloat` and `Complex{BigFloat}`, on general, singular, triangular, diagonal and Hermitian matrices, gives the dense result type in every case that runs; only the `Float32` and `ComplexF32` LU cases changed. Potential breakage: `cond` of a nonsingular `Float32` or `ComplexF32` sparse matrix that is not triangular, diagonal or positive definite now returns `Float32` instead of `Float64`, as dense `cond` does. Written by Claude Code. ๐Ÿค– Generated with [Claude Code](https://claude.com/claude-code) --------- Co-authored-by: Claude Opus 5.5 (cherry picked from commit eb474418bef53775befb309b11b6ece3e5f6d8b9) --- src/linalg.jl | 36 +++++++++++++++++++++++++----------- test/solvers/issues.jl | 7 ++++++- 2 files changed, 31 insertions(+), 12 deletions(-) diff --git a/src/linalg.jl b/src/linalg.jl index ffb19d1b..57a58bb2 100644 --- a/src/linalg.jl +++ b/src/linalg.jl @@ -993,18 +993,21 @@ another order than the dense one, so `cond(A, p)` can be finite where `cond(Arra is `Inf`. """ function cond(A::AbstractSparseMatrixCSC, p::Real=2) + # the type of dense `cond`, whichever precision the factorization of `A` works in + T = float(real(eltype(A))) # as for a dense matrix; the estimator of the inverse's norm needs at least one column - (p == 1 || p == Inf) && checksquare(A) == 0 && return zero(real(float(eltype(A)))) + (p == 1 || p == Inf) && checksquare(A) == 0 && return zero(T) if p == 1 || p == Inf - # a singular matrix has an infinite condition number, as for a dense matrix; a - # triangular or diagonal `factorize(A)` reports it only in the solve - normAinv = try - p == 1 ? opnormestinv(A) : opnormestinv(copy(A')) + B = p == 1 ? A : copy(A') + # a singular matrix has an infinite condition number, as for a dense matrix + F = try + factorize(B) catch e e isa SingularException || rethrow() - return convert(float(real(eltype(A))), Inf) + return convert(T, Inf) end - return opnorm(A, p) * normAinv + _issingularfactor(F) && return convert(T, Inf) + return convert(T, opnorm(A, p) * _opnormestinv(B, F, min(2, size(B, 1)))) elseif p == 2 throw(ArgumentError("only 1- and Inf-norm condition numbers are implemented for sparse matrices, for 2-norm try cond(Array(A), 2) instead")) else @@ -1012,11 +1015,22 @@ function cond(A::AbstractSparseMatrixCSC, p::Real=2) end end -function opnormestinv(A::AbstractSparseMatrixCSC{T}, t::Integer = min(2,maximum(size(A)))) where T +# `factorize` throws for a singular LU or Hermitian factorization, but leaves a triangular +# or diagonal one to the solve +_issingularfactor(F) = false +_issingularfactor(D::Diagonal) = any(iszero, D.diag) +_issingularfactor(L::LowerTriangular{<:Any,<:AbstractSparseMatrixCSC}) = !iszero(_trisingular(false, parent(L))) +_issingularfactor(U::UpperTriangular{<:Any,<:AbstractSparseMatrixCSC}) = !iszero(_trisingular(true, parent(U))) + +function opnormestinv(A::AbstractSparseMatrixCSC, t::Integer = min(2,maximum(size(A)))) + checksquare(A) + return _opnormestinv(A, factorize(A), t) +end + +# the estimate for the factorization `F` of `A` +function _opnormestinv(A::AbstractSparseMatrixCSC{T}, F, t::Integer) where T maxiter = 5 - # Check the input - n = checksquare(A) - F = factorize(A) + n = size(A, 1) if t <= 0 throw(ArgumentError("number of blocks must be a positive integer")) end diff --git a/test/solvers/issues.jl b/test/solvers/issues.jl index 60e094f9..2eff6dfa 100644 --- a/test/solvers/issues.jl +++ b/test/solvers/issues.jl @@ -256,8 +256,13 @@ end @testset "cond of a singular matrix is Inf" begin @test cond(sparse([1.0 2; 3 6]), 1) === Inf - # the triangular solve, not the factorization, finds this one singular + # `factorize` leaves the zero on the diagonal of a triangular or diagonal factor @test cond(sparse([1.0 0; 2 0]), Inf) === Inf + @test cond(sparse([1.0 0; 0 0]), 1) === Inf + # Hermitian, so `factorize` tries Cholesky and LDLt before the LU finds it singular + @test cond(sparse([1.0 1; 1 1]), 1) === Inf + # the LU works in double precision, but the result is of the matrix's precision + @test cond(sparse(Float32[1 2; 3 4]), 1) isa Float32 end end # module From 376ea564fb9627dcade26dc45085c9aac545f4f1 Mon Sep 17 00:00:00 2001 From: "Viral B. Shah" Date: Mon, 26 Oct 2026 23:20:37 -0400 Subject: [PATCH 25/30] Make whole views of a sparse vector or column behave like the parent for any index type (#1119) MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit A whole view of a sparse vector, `view(x, :)`, or of a sparse column, `view(A, :, j)`, should behave like the vector it views. Two things stopped it. `view(x, :)` of a `SparseVector{Tv,Int32}` is not a `SparseVectorView`, so it misses every O(nnz) method on that alias: `sum`, `maximum`, `norm` and `dot` visit all `n` elements, `rmul!` and `lmul!` take the dense fallback, `kron` is dense and `findnz` throws. And for any index type, `SparseVector(v)` of a whole view, which `sparse`, `convert`, `vcat`, `hcat` and `kron` all go through, reads it as a dense vector: O(n), stored zeros dropped, `Int` indices. ```julia using SparseArrays x = sparsevec(Int32[2, 4, 5], [1.0, 0.0, 3.0], 6); v = view(x, :) A = sparse(Int32[1, 3], Int32[1, 1], [1.0, 0.0], 4, 2); c = view(A, :, 1) ``` **Before** ``` SparseVector(v) = sparsevec([2, 5], [1.0, 3.0], 6) sparse(v) = sparsevec([2, 5], [1.0, 3.0], 6) vcat(v, v) = sparsevec([2, 5, 8, 11], [1.0, 3.0, 1.0, 3.0], 12) typeof(kron(v, x)) = Vector{Float64} typeof(SparseVector(c)) = SparseVector{Float64, Int64} findnz(v): MethodError: no method matching findnz(::SubArray{Float64, 1, SparseVector{Float64, Int32}, Tuple{Base.Slice{Base.OneTo{Int32}}}, false}) ``` **After** ``` SparseVector(v) = sparsevec(Int32[2, 4, 5], [1.0, 0.0, 3.0], 6) sparse(v) = sparsevec(Int32[2, 4, 5], [1.0, 0.0, 3.0], 6) vcat(v, v) = sparsevec(Int32[2, 4, 5, 8, 10, 11], [1.0, 0.0, 3.0, 1.0, 0.0, 3.0], 12) kron(v, x) = sparsevec(Int32[8, 10, 11, 20, 22, 23, 26, 28, 29], [1.0, 0.0, 3.0, 0.0, 0.0, 0.0, 3.0, 0.0, 9.0], 36) typeof(SparseVector(c)) = SparseVector{Float64, Int32} findnz(v) = (Int32[2, 4, 5], [1.0, 0.0, 3.0]) ``` For `v32 = view(sparsevec(Int32[1, n], [1.0, 2.0], n), :)` and `v64` likewise, with `n = 10^7`, the minimum of three runs and the `@allocated` of a warm call: | `f(v)` | `v32`, main | `v32`, branch | `v64`, main | `v64`, branch | |---------------------|-------------------|------------------|-------------------|------------------| | `sum` | 0.020 s, 0 B | 2.0e-8 s, 0 B | 3.0e-8 s, 0 B | 3.0e-8 s, 0 B | | `maximum` | 0.047 s, 0 B | 2.9e-8 s, 0 B | 3.0e-8 s, 0 B | 2.0e-8 s, 0 B | | `dot(v, v)` | 0.035 s, 0 B | 3.0e-8 s, 0 B | 2.1e-8 s, 0 B | 2.0e-8 s, 0 B | | `count(!iszero, v)` | 0.021 s, 0 B | 2.0e-8 s, 0 B | 0.021 s, 0 B | 2.0e-8 s, 0 B | | `findnz` | `MethodError` | 3.0e-8 s, 144 B | 3.0e-8 s, 160 B | 3.0e-8 s, 160 B | | `SparseVector` | 0.027 s, 19 MiB | 5.0e-8 s, 144 B | 0.023 s, 19 MiB | 4.0e-8 s, 160 B | | `vcat(v, v)` | 0.057 s, 38 MiB | 3.1e-7 s, 880 B | 0.052 s, 38 MiB | 2.6e-7 s, 960 B | A `SparseVector` stores its length as `Ti`, so `axes(x32)` is `(OneTo{Int32}(n),)` and `view(x32, :)` slices a `Slice{OneTo{Int32}}`, while `SparseVectorView` named `Slice{OneTo{Int}}` only. A matrix stores its size as `Int`, so `SparseColumnView` already matched for every index type. The fix widens `SparseVectorView` to `<:Tuple{Base.Slice{<:Base.OneTo}}`; no method on the alias assumes `Int` axes, since they read the parent through `nonzeros`, `nonzeroinds` and `length`. The `SparseVector` copy constructors now take `SparseVectorOrView` rather than `AbstractCompressedVector`, so a whole view or column view is copied from its stored entries with the parent's index type, and `sparsevec` of one copies it through them. `count` on a `SparseVector` already counted `f` over the stored values plus `f` of one zero times the number of unstored entries; that method (`_simple_count`) now also covers the whole and range views, and takes the zero from `_densezero`. This touches dispatch. Potential breakage, in each case towards what the parent already gives: - For a whole view or column view of any index type, `SparseVector(v)`, `sparse(v)`, `convert(SparseVector, v)` and sparse `vcat`, `hcat`, `hvcat` and `kron` that take one keep its stored zeros, so `nnz` of the result can be larger than before. - They also keep the parent's index type instead of `Int`. - For a whole view `v` of a `SparseVector` with a non-`Int` index type: `nonzeros(v)` is the parent's `Vector` rather than a `SubArray` of it; `nonzeroinds(v)` is the parent's own index vector rather than a fresh copy, so writing to it changes the parent; `findfirst`, `findlast`, `findnext` and `findprev` return `Ti` rather than `Int`; `diff(v)` keeps a stored zero; `kron(v, x)` and `kron(v, A)` are sparse rather than dense; and `A * v` with an `Int`-indexed `A` is a `SparseVector{Tv,Int}` rather than a `SparseVector{Tv,Ti}`. `length` keeps returning `Ti`, as #176 (JuliaLang/julia#45937) decided. Range views such as `view(x, 2:4)` still convert through the dense path, and `spdiagm(0 => view(x, :))` still stores every diagonal entry; both are left for separate PRs. ๐Ÿค– Generated with [Claude Code](https://claude.com/claude-code) --------- Co-authored-by: Claude Opus 5.5 (cherry picked from commit 61f75e00a0872438be41f43df04d4a8ea0795af9) --- src/SparseArrays.jl | 5 +++-- src/sparsevector.jl | 11 ++++++----- test/sparsevector.jl | 15 ++++++++++++--- 3 files changed, 21 insertions(+), 10 deletions(-) diff --git a/src/SparseArrays.jl b/src/SparseArrays.jl index e101d332..4a43ebd4 100644 --- a/src/SparseArrays.jl +++ b/src/SparseArrays.jl @@ -160,9 +160,10 @@ const SparseMatrixCSCSubArray{Tv,Ti} = SubArray{Tv,2,<:AbstractSparseMatrixCSC{T const SparseMatrixCSCOrSubArray{Tv,Ti} = Union{AbstractSparseMatrixCSC{Tv,Ti}, SparseMatrixCSCSubArray{Tv,Ti}} # Whole-column views of sparse matrices and whole views of sparse vectors share the -# sparse vector interface. +# sparse vector interface. The length of a `SparseVector` has its index type, so a whole +# view of one slices a `OneTo{Ti}`. const SparseColumnView{Tv,Ti} = SubArray{Tv,1,<:AbstractSparseMatrixCSC{Tv,Ti},Tuple{Base.Slice{Base.OneTo{Int}},Int},false} -const SparseVectorView{Tv,Ti} = SubArray{Tv,1,<:AbstractSparseVector{Tv,Ti},Tuple{Base.Slice{Base.OneTo{Int}}},false} +const SparseVectorView{Tv,Ti} = SubArray{Tv,1,<:AbstractSparseVector{Tv,Ti},<:Tuple{Base.Slice{<:Base.OneTo}},false} const SparseVectorOrView{Tv,Ti} = Union{AbstractCompressedVector{Tv,Ti}, SparseColumnView{Tv,Ti}, SparseVectorView{Tv,Ti}} const AdjOrTransSparseVectorOrView{Tv,Ti} = AdjOrTrans{Tv, <:SparseVectorOrView{Tv,Ti}} # view of a unit range of a sparse vector's indices diff --git a/src/sparsevector.jl b/src/sparsevector.jl index a5d4fe42..781eba6e 100644 --- a/src/sparsevector.jl +++ b/src/sparsevector.jl @@ -63,8 +63,8 @@ length(x::FixedSparseVector) = getfield(x, :n) size(x::SparseVector) = (getfield(x, :n),) size(x::FixedSparseVector) = (getfield(x, :n),) -function Base._simple_count(f, x::AbstractCompressedVector, init::T) where T - init + T(count(f, nonzeros(x)) + f(zero(eltype(x)))*(length(x) - nnz(x))) +function Base._simple_count(f, x::Union{SparseVectorOrView,SparseVectorPartialView}, init::T) where T + init + T(count(f, nonzeros(x)) + f(_densezero(x))*(length(x) - nnz(x))) end # implement the nnz - nzrange - nonzeros - rowvals interface for sparse vectors @@ -555,6 +555,7 @@ sparsevec(a::AbstractVector{T}) where {T} = SparseVector{T, Int}(a) sparsevec(a::AbstractArray) = sparsevec(vec(a)) sparsevec(a::AbstractSparseMatrixCSC) = a[:] sparsevec(a::AbstractSparseVector) = vec(a) +sparsevec(a::Union{SparseColumnView,SparseVectorView}) = SparseVector(a) sparse(a::AbstractVector) = sparsevec(a) function _dense2indval!(nzind::Vector{Ti}, nzval::Vector{Tv}, s::AbstractArray) where {Tv,Ti} @@ -595,9 +596,9 @@ SparseVector{Tv}(s::AbstractVector) where {Tv} = SparseVector{Tv,Int}(s) SparseVector(s::AbstractVector{Tv}) where {Tv} = SparseVector{Tv,Int}(s) # copy-constructors -SparseVector(s::AbstractCompressedVector{Tv,Ti}) where {Tv,Ti} = SparseVector{Tv,Ti}(s) -SparseVector{Tv}(s::AbstractCompressedVector{<:Any,Ti}) where {Tv,Ti} = SparseVector{Tv,Ti}(s) -function SparseVector{Tv,Ti}(s::AbstractCompressedVector) where {Tv,Ti} +SparseVector(s::SparseVectorOrView{Tv,Ti}) where {Tv,Ti} = SparseVector{Tv,Ti}(s) +SparseVector{Tv}(s::SparseVectorOrView{<:Any,Ti}) where {Tv,Ti} = SparseVector{Tv,Ti}(s) +function SparseVector{Tv,Ti}(s::SparseVectorOrView) where {Tv,Ti} copyind = Vector{Ti}(nonzeroinds(s)) copynz = Vector{Tv}(nonzeros(s)) SparseVector{Tv,Ti}(length(s), copyind, copynz) diff --git a/test/sparsevector.jl b/test/sparsevector.jl index b424b7c2..deaead50 100644 --- a/test/sparsevector.jl +++ b/test/sparsevector.jl @@ -1805,9 +1805,9 @@ end @test ((1 + im) * v)::T == (1 + im) * d @test copy(v)::T == d end - # the full-column and full-vector views scale in place through the parent's storage - # (`view(x, :)` is a `SparseVectorView` only when the parent's axes are `Int`-sized) - for v in (view(A, :, 1), (Ti === Int ? (view(x, :),) : ())...), a in (1.0 + im,) + # the full-column and full-vector views scale in place through the parent's storage, + # also when the parent's axes are not `Int` + for v in (view(A, :, 1), view(x, :), view(SparseVector{ComplexF64,Int128}(x), :)), a in (1.0 + im,) @test which(*, (typeof(v), typeof(a))).module === SparseArrays @test which(*, (typeof(a), typeof(v))).module === SparseArrays @test which(/, (typeof(v), typeof(a))).module === SparseArrays @@ -1831,6 +1831,15 @@ end @test @allocated(f(v)) < 10^4 end end + # `count` calls `f` on the stored entries and on one zero + calls = Ref(0) + @test count(z -> (calls[] += 1; !iszero(z)), view(SparseVector{ComplexF64,Int32}(y), :)) == 2 && calls[] == 3 + # a whole view converts with the parent's index type and stored zeros + z = SparseVector{ComplexF64,Int32}(x) + C = SparseMatrixCSC{ComplexF64,Int32}(A) + for (v, u) in ((view(z, :), z), (view(C, :, 1), C[:, 1])), f in (SparseVector, sparse, w -> vcat(w, w), w -> kron(w, z)) + @test exact_equal(f(v), f(u)) + end end end From 2fe7e0a766bb88a74bfba4d1184db8763aa506d4 Mon Sep 17 00:00:00 2001 From: "Viral B. Shah" Date: Mon, 26 Oct 2026 23:30:41 -0400 Subject: [PATCH 26/30] Build `spdiagm` from the stored entries of a sparse vector view (#1133) MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit `spdiagm` places only the stored entries of a sparse vector, but a view of one, `view(x, :)`, `view(x, 2:9)` or `view(A, :, j)`, is not an `AbstractSparseVector`, so it takes the dense branch and stores every element of the diagonal, zeros included. ```julia using SparseArrays x = sparsevec([2, 4, 5], [1.0, 0.0, 3.0], 6) spdiagm(0 => view(x, :)) ``` **Before** ``` 6ร—6 SparseMatrixCSC{Float64, Int64} with 6 stored entries: 0.0 โ‹… โ‹… โ‹… โ‹… โ‹… โ‹… 1.0 โ‹… โ‹… โ‹… โ‹… โ‹… โ‹… 0.0 โ‹… โ‹… โ‹… โ‹… โ‹… โ‹… 0.0 โ‹… โ‹… โ‹… โ‹… โ‹… โ‹… 3.0 โ‹… โ‹… โ‹… โ‹… โ‹… โ‹… 0.0 ``` **After** ``` 6ร—6 SparseMatrixCSC{Float64, Int64} with 3 stored entries: โ‹… โ‹… โ‹… โ‹… โ‹… โ‹… โ‹… 1.0 โ‹… โ‹… โ‹… โ‹… โ‹… โ‹… โ‹… โ‹… โ‹… โ‹… โ‹… โ‹… โ‹… 0.0 โ‹… โ‹… โ‹… โ‹… โ‹… โ‹… 3.0 โ‹… โ‹… โ‹… โ‹… โ‹… โ‹… โ‹… ``` For `x = sparsevec([1, n], [1.0, 2.0], n)` with `n = 10^6`, `spdiagm(0 => v)`, the minimum of five runs: | `v` | main | branch | |----------------|------------------------------|----------------------------| | `x` | 0.0068 s, 22.9 MiB, nnz 2 | 0.0067 s, 22.9 MiB, nnz 2 | | `view(x, :)` | 0.029 s, 76.3 MiB, nnz 10^6 | 0.0066 s, 22.9 MiB, nnz 2 | | `view(x, 1:n)` | 0.029 s, 76.3 MiB, nnz 10^6 | 0.0071 s, 22.9 MiB, nnz 2 | The helpers `_nnz`, `_indices` and `_nzvals` that `spdiagm_internal` reads each diagonal through had a sparse method for `AbstractSparseVector` only. They now take `SparseQVecOperand`, the existing alias for sparse vectors and their whole-column, whole-vector and range views, all of which provide `nnz`, `nonzeroinds` and `nonzeros`. The time that remains is the column pointers of the `n`ร—`n` result. Potential breakage: `spdiagm` of a sparse vector view no longer stores the zeros between the view's stored entries, so `nnz` of the result is that of the view. ๐Ÿค– Generated with [Claude Code](https://claude.com/claude-code) Co-authored-by: Claude Opus 5.5 (cherry picked from commit aee932add0f289a6ca32e702b6c963ad41596052) --- src/SparseArrays.jl | 3 ++- src/constructors.jl | 6 +++--- test/linalg.jl | 7 +++++++ 3 files changed, 12 insertions(+), 4 deletions(-) diff --git a/src/SparseArrays.jl b/src/SparseArrays.jl index 4a43ebd4..4f319880 100644 --- a/src/SparseArrays.jl +++ b/src/SparseArrays.jl @@ -209,7 +209,8 @@ const LinAlgLeftQs = Union{HessenbergQ,QRCompactWYQ,QRPackedQ} # transpose of the operand. Sparse matrices and their views, with or without an adjoint or # transpose, and a transposed sparse vector as a one-row matrix; and sparse vectors and their # views. An adjoint sparse vector is left to LinearAlgebra's generic (Q' * u')', which keeps -# the result an Adjoint as for a dense vector. +# the result an Adjoint as for a dense vector. `spdiagm` reads the vector operands by their +# stored entries. const SparseQMatOperand = Union{SparseMatrixCSCOrView, AdjOrTrans{<:Any,<:SparseMatrixCSCOrView}, Transpose{<:Any,<:SparseVectorOrView}} const SparseQVecOperand = Union{AbstractSparseVector, SparseVectorOrView, SparseVectorPartialView} diff --git a/src/constructors.jl b/src/constructors.jl index bf4452bc..bf3beb05 100644 --- a/src/constructors.jl +++ b/src/constructors.jl @@ -948,10 +948,10 @@ sparse(s::UniformScaling, m::Integer, n::Integer) = sparse(s, Dims((m, n))) sparse(::Type{Tv}, s::UniformScaling, m::Integer, n::Integer) where {Tv} = SparseMatrixCSC{Tv}(s, Dims((m, n))) sparse(::Type{Tv}, ::Type{Ti}, s::UniformScaling, m::Integer, n::Integer) where {Tv, Ti} = SparseMatrixCSC{Tv, Ti}(s, Dims((m, n))) -_nnz(v::AbstractSparseVector) = nnz(v) +_nnz(v::SparseQVecOperand) = nnz(v) _nnz(v::AbstractVector) = length(v) -function _indices(v::AbstractSparseVector, row, col) +function _indices(v::SparseQVecOperand, row, col) ix = nonzeroinds(v) return (row .+ ix, col .+ ix) end @@ -960,7 +960,7 @@ function _indices(v::AbstractVector, row, col) return (row+1:row+veclen, col+1:col+veclen) end -_nzvals(v::AbstractSparseVector) = nonzeros(v) +_nzvals(v::SparseQVecOperand) = nonzeros(v) _nzvals(v::AbstractVector) = v # Promoted element type of the diagonals, mirroring `Base.promote_eltypeof` diff --git a/test/linalg.jl b/test/linalg.jl index 3c3adc0e..24bb3f90 100644 --- a/test/linalg.jl +++ b/test/linalg.jl @@ -718,6 +718,13 @@ end x = fixturevec(Float64, 10); y = ones(Int, 9) @test spdiagm(0 => x, 1 => y)::SparseMatrixCSC{Float64,Int} == Bidiagonal(Vector(x), ones(9), :U) @test nnz(spdiagm(0 => x, 1 => y)) == length(y) + nnz(x) + @static if COMPREHENSIVE + # so do views of a sparse vector or column + A = fixture(Float64, 10, 3) + for v in (view(x, :), view(x, 2:9), view(A, :, 1)) + @test exact_equal(spdiagm(-1 => v), spdiagm(-1 => copy(v))) + end + end end @static if COMPREHENSIVE From d06853ba1cc14de2cca3bccab86fefdc893a28c9 Mon Sep 17 00:00:00 2001 From: "Viral B. Shah" Date: Mon, 26 Oct 2026 23:42:44 -0400 Subject: [PATCH 27/30] Take the index type of a mixed sparse and dense concatenation from its sparse arguments only (#1127) MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit The index type of a concatenation that mixes sparse and dense arguments depended on where the dense ones were. `hcat`, `vcat` and `cat` counted a dense first argument as `Int`, and `hvcat` counted every dense block as `Int`. So the same blocks in another order gave another index type, or threw. The all-sparse methods already promote over every argument. **Before** ```julia julia> n = 2; A32 = spzeros(Float64, Int32, 1, n); A64 = spzeros(Float64, Int64, 3_000_000_000, n); d = ones(1, n); julia> v32 = spzeros(Float64, Int32, 1); v64 = spzeros(3_000_000_000); julia> vcat(A32, A64, d) ERROR: ArgumentError: number of rows (m = 3000000002) does not fit in Ti = Int32 julia> vcat(v32, v64, ones(1)) ERROR: InexactError: trunc(Int32, 3000000002) julia> indtype(vcat(A32, d)), indtype(vcat(d, A32)) (Int32, Int64) julia> indtype([A32; d]), indtype([A32 A32; d d]) (Int32, Int64) julia> indtype([A32' I]), indtype([I A32']) (Int32, Int64) ``` **After** ```julia julia> vcat(A32, A64, d) 3000000002ร—2 SparseMatrixCSC{Float64, Int64} with 2 stored entries julia> vcat(v32, v64, ones(1)) 3000000002-element SparseVector{Float64, Int64} with 1 stored entry: [3000000002] = 1.0 julia> indtype(vcat(A32, d)), indtype(vcat(d, A32)) (Int32, Int32) julia> indtype([A32; d]), indtype([A32 A32; d d]) (Int32, Int32) julia> indtype([A32' I]), indtype([I A32']) (Int32, Int32) ``` The index type was read after `_makesparse` had turned every dense array into an `Int`-indexed sparse matrix. It is now read from the original arguments by one rule, `_catidxtype`, which `hcat`, `vcat`, `cat`, `hvcat` and the `sparse_*` functions share. The rule promotes over the sparse arguments, wherever they are, and gives `Int` when none is sparse. A dense array, a number, a structured matrix or `I` adds nothing. The size of the result is checked against the index type, and that check now covers vector results too, which used to throw an `InexactError`. `I` used to become `sparse(J, n, n)`, which has `Int` indices, and now becomes a `Diagonal`. The Concatenation section of the manual states the rule. Potential breakage, all of it a change of index type in code that works today: - A dense block or a number in first place no longer widens the result to `Int`. `vcat(d, A32)`, `[1.0 v32'; v32 A32]` and `[I A32]` now give `Int32`, as `vcat(A32, d)` already does. A result too large for that index type now throws, in either order. - `hvcat` with dense blocks no longer counts them as `Int`, so `[A32 A32; d d]` gives `Int32`. - Sparse arguments of different index types now give the promoted type, not that of the first: `vcat(A8, A32, d)` gives `Int32`, as `vcat(A8, A32)` already does. Written by Claude Code. ๐Ÿค– Generated with [Claude Code](https://claude.com/claude-code) --------- Co-authored-by: Claude Opus 5.5 (cherry picked from commit b575c39d0cfac1998fcba1fe0911a01898b968dc) --- docs/src/index.md | 5 ++ src/concatenation.jl | 116 +++++++++++++++++++++--------------------- test/concatenation.jl | 27 +++++++--- 3 files changed, 83 insertions(+), 65 deletions(-) diff --git a/docs/src/index.md b/docs/src/index.md index f76c3c7c..a2d99dbc 100644 --- a/docs/src/index.md +++ b/docs/src/index.md @@ -214,6 +214,11 @@ sparse. A dense argument contributes its nonzeros, and a `Diagonal` or `I` its d [`sparse_hcat`](@ref), [`sparse_vcat`](@ref) and [`sparse_hvcat`](@ref) return a sparse array even when no argument is sparse. +The index type of the result is promoted over the sparse arguments, wherever they are, +and is `Int` when none is sparse. A dense argument, a number or `I` adds none, so with a +dense `D`, `[A32; D]` and `[D; A32]` both have the `Int32` indices of `A32`. A result too +large for that index type throws an `ArgumentError`. + #### [Reductions](@id man-sparse-reductions) `sum(A)` and the other reductions over all entries visit the stored entries only. Along a diff --git a/src/concatenation.jl b/src/concatenation.jl index 8646554b..73d441a6 100644 --- a/src/concatenation.jl +++ b/src/concatenation.jl @@ -292,16 +292,15 @@ end # Concatenations involving un/annotated sparse/special matrices/vectors should yield sparse arrays -# the output array type is determined by the first element of the to be concatenated objects -# if this is a Number, the output would be dense by the fallback abstractarray.jl code (see cat_similar) -# so make sure that if that happens, the "array" is sparse (if more sparse arrays are involved, of course) +# `_sparse` makes the first argument of a concatenation that has to be sparse, with a +# number standing for a one-element vector _sparse(x::Number) = sparsevec([1], [x], 1) _sparse(A) = _makesparse(A) _makesparse(x::Number) = x _makesparse(x::AbstractVector) = convert(SparseVector, x)::SparseVector _makesparse(x::AbstractMatrix) = convert(SparseMatrixCSC, x)::SparseMatrixCSC # a `UniformScaling` has no size of its own: `LinearAlgebra._hcat`/`_vcat`/`_hvcat` size it -# from its neighbours and then call `promote_to_arrays_` below to make it sparse +# from its neighbours and then call `promote_to_arrays_` below to make it a `Diagonal` _makesparse(J::UniformScaling) = J anysparse() = false anysparse(X) = X isa AbstractArray && issparse(X) @@ -321,13 +320,13 @@ _densesparse(x) = anysparse(x) ? Array(x) : x const _SparseVecConcatGroup = Union{Vector, AbstractSparseVector} function hcat(X::_SparseVecConcatGroup...) if _concatsparse(X...) - return cat(map(sparse, X)...; dims=Val(2)) + return _sparse_cat_t(Val(2), promote_eltype(X...), X...) end return cat(map(_densesparse, X)...; dims=Val(2)) end function vcat(X::_SparseVecConcatGroup...) if _concatsparse(X...) - return cat(map(sparse, X)...; dims=Val(1)) + return _sparse_cat_t(Val(1), promote_eltype(X...), X...) end return cat(map(_densesparse, X)...; dims=Val(1)) end @@ -339,22 +338,30 @@ end # ambiguity is introduced for arrays that are not sparse (#431). const _SparseConcatGroup = Union{AbstractVecOrMat,Number} -# Base's `_cat_t` takes the output type from its first argument, so with a leading number -# it would build a dense array. Choose the destination from the first array instead, and -# keep the number as is so that it fills its block the way it does in dense -# concatenation (#383). A leading number still widens the index type to at least `Int`, -# as the one-element sparse vector standing in for it used to. `X` has already been -# through `_makesparse`. -_catleader(X1::AbstractArray, X...) = X1 -_catleader(X1::Number, X...) = _catleader(X...) -_catleader(X1::Number) = _sparse(X1) -_catdest(::Type{T}, shape, X1::AbstractArray, X...) where {T} = similar(X1, T, shape) -function _catdest(::Type{T}, shape, X1::Number, X...) where {T} - A = _catleader(X1, X...) - return similar(A, T, promote_type(Int, indtype(A)), shape) -end -_catidxtype(X1::AbstractArray, X...) = indtype(X1) -_catidxtype(X1::Number, X...) = promote_type(Int, indtype(_catleader(X1, X...))) +# The index type of a concatenation is promoted over its sparse arguments, wherever they +# are, and is `Int` when none is sparse. A dense array, a number, a structured matrix or a +# `UniformScaling` has no index type of its own, so it adds none, and the result does not +# depend on the order of the arguments. The size of the result is checked against it. +function _catidxtype(X::Vararg{Any,N}) where {N} + Ti = _sparseidxtype(X...) + return Ti === Union{} ? Int : Ti +end +_sparseidxtype() = Union{} +_sparseidxtype(x, X::Vararg{Any,N}) where {N} = + promote_type(anysparse(x) ? _wrappedidxtype(x) : Union{}, _sparseidxtype(X...)) +_wrappedidxtype(x::AbstractSparseArray) = indtype(x) +function _wrappedidxtype(x::AbstractArray) + p = parent(x) + return p === x ? indtype(_makesparse(x)) : _wrappedidxtype(p) +end +# The destination of a vector result is made here and not by Base's `_cat_t`, which takes +# its type from the first argument and so would build a dense array after a leading +# number. A number is kept as is, so that it fills its block as in dense concatenation (#383). +function _catvector(::Type{T}, n, Xin::Tuple) where {T} + Ti = _catidxtype(Xin...) + sparse_check_Ti(n, 1, Ti) + return spzeros(T, Ti, n) +end # The compiled method serves every value of `dims`, so `dims2cat(dims)` has an unknown # length, and `similar` of a sparse array cannot be resolved for a shape of unknown length, # which `juliac --trim` rejects. A sparse result has one or two dimensions, so each gets a @@ -371,10 +378,11 @@ end # Base's `__cat` writes each block into the destination with `setindex!`, which does not # store a zero, so a matrix result is assembled from the stored entries of the blocks. It # keeps the stored zeros of a sparse block, as `hcat` and `vcat` of sparse matrices do. -function _sparse_cat_t_shaped(catdims::Tuple{Vararg{Bool}}, ::Type{T}, X::Vararg{Any,N}) where {T,N} +function _sparse_cat_t_shaped(catdims::Tuple{Vararg{Bool}}, ::Type{T}, Xin::Vararg{Any,N}) where {T,N} + X = map(_makesparse, Xin) shape = Base.cat_size_shape(catdims, X...) - length(shape) == 1 && return Base.__cat(_catdest(T, shape, X...), shape, catdims, X...) - Ti = _catidxtype(X...) + length(shape) == 1 && return Base.__cat(_catvector(T, shape[1], Xin), shape, catdims, X...) + Ti = _catidxtype(Xin...) sparse_check_Ti(shape[1], shape[2], Ti) vertical, horizontal = catdims[1], length(catdims) == 2 && catdims[2] # a number fills its block, which spans the result in a dimension not concatenated @@ -397,9 +405,6 @@ end _catdense(x::AbstractArray) = Array(x) _catdense(x) = x _sparse_typed_hcat(::Type{T}, X...) where {T} = _sparse_cat_t(Val(2), T, X...) -# with only vectors the result is a vector, which Base's `typed_vcat` builds by assigning -# whole sparse vectors, which keeps their stored zeros -_sparse_typed_vcat(::Type{T}, X::AbstractVector...) where {T} = Base.typed_vcat(T, X...) _sparse_typed_vcat(::Type{T}, X...) where {T} = _sparse_cat_t(Val(1), T, X...) # `Vararg{_SparseConcatGroup,N}` makes Julia compile `cat_internal` for each argument @@ -409,51 +414,42 @@ _sparse_typed_vcat(::Type{T}, X...) where {T} = _sparse_cat_t(Val(1), T, X...) Base.@constprop :aggressive function cat_internal(dims, X1::_SparseConcatGroup, X::Vararg{_SparseConcatGroup,N}) where {N} T = promote_eltype(X1, X...) if _concatsparse(X1, X...) - return _sparse_cat_t(dims, T, _makesparse(X1), map(_makesparse, X)...) + return _sparse_cat_t(dims, T, X1, X...) end return Base._cat_t(dims, T, _densesparse(X1), map(_densesparse, X)...) end function hcat_internal(X1::_SparseConcatGroup, X::_SparseConcatGroup...) T = promote_eltype(X1, X...) if _concatsparse(X1, X...) - return _sparse_typed_hcat(T, _makesparse(X1), map(_makesparse, X)...) + return _sparse_typed_hcat(T, X1, X...) end return Base.typed_hcat(T, _densesparse(X1), map(_densesparse, X)...) end function vcat_internal(X1::_SparseConcatGroup, X::_SparseConcatGroup...) T = promote_eltype(X1, X...) if _concatsparse(X1, X...) - return _sparse_typed_vcat(T, _makesparse(X1), map(_makesparse, X)...) + return _sparse_typed_vcat(T, X1, X...) end return Base.typed_vcat(T, _densesparse(X1), map(_densesparse, X)...) end # `Vararg{_SparseConcatGroup,N}` for the same reason as in `cat_internal` function hvcat_internal(rows::Tuple{Vararg{Int}}, X1::_SparseConcatGroup, X::Vararg{_SparseConcatGroup,N}) where {N} if _concatsparse(X1, X...) - return _sparse_hvcat(rows, _makesparse(X1), map(_makesparse, X)...) + return _sparse_hvcat(rows, X1, X...) end return Base.typed_hvcat(Base.promote_eltypeof(X1, X...), rows, _densesparse(X1), map(_densesparse, X)...) end # `_hvcat_csc` reads the blocks by index. Splitting them into a tuple per block row would # give tuples whose length depends on the value of `rows`, so their types could not be # inferred and `juliac --trim` could not resolve the calls on them. -function _sparse_hvcat(rows::Tuple{Vararg{Int}}, X::Vararg{AbstractSparseMatrixCSC,N}) where {N} - return _hvcat_csc(promote_eltype(X...), promote_idxtype(X...), rows, X...) -end # A vector is an `nร—1` block and a number a `1ร—1` block, as in dense `hvcat`. A number is # stored unless scalar `setindex!` would leave it implicit, and a sparse vector keeps its -# stored zeros, as sparse matrix blocks do. A leading number widens the index type to at -# least `Int`, as it does for `hcat` and `vcat` (#383). -function _sparse_hvcat(rows::Tuple{Vararg{Int}}, X::Vararg{Any,N}) where {N} - Tv = promote_eltype(X...) - Ti = _hvcat_idxtype(X...) - return _hvcat_csc(Tv, Ti, rows, map(x -> _hvcat_block(Tv, Ti, x), X)...) -end -_hvcat_idxtype(X1::Number, X...) = promote_type(Int, _blocks_idxtype(X...)) -_hvcat_idxtype(X...) = _blocks_idxtype(X...) -_blocks_idxtype() = Union{} -_blocks_idxtype(x::Number, X...) = _blocks_idxtype(X...) -_blocks_idxtype(x::AbstractSparseArray, X...) = promote_type(indtype(x), _blocks_idxtype(X...)) +# stored zeros, as sparse matrix blocks do. +function _sparse_hvcat(rows::Tuple{Vararg{Int}}, Xin::Vararg{Any,N}) where {N} + Tv = promote_eltype(Xin...) + Ti = _catidxtype(Xin...) + return _hvcat_csc(Tv, Ti, rows, map(x -> _hvcat_block(Tv, Ti, _makesparse(x)), Xin)...) +end _hvcat_block(::Type, ::Type, B::AbstractSparseMatrixCSC) = B function _hvcat_block(::Type, ::Type, x::AbstractSparseVector{<:Any,Ti}) where {Ti} return SparseMatrixCSC(length(x), 1, Ti[1, nnz(x) + 1], nonzeroinds(x), nonzeros(x)) @@ -597,9 +593,17 @@ function Base.repeat(v::AbstractSparseVector, m, n) end -# make sure UniformScaling objects are converted to sparse matrices for concatenation +# A `UniformScaling` in a sparse concatenation becomes a `Diagonal`, which the concatenation +# makes sparse. As a structured block, not a sparse one, it does not widen the index type. promote_to_array_type(A::Tuple{Vararg{Union{_SparseConcatGroup,UniformScaling}}}) = _concatsparse(A...) ? SparseMatrixCSC : Matrix -promote_to_arrays_(n::Int, ::Type{SparseMatrixCSC}, J::UniformScaling) = sparse(J, n, n) +promote_to_arrays_(n::Int, ::Type{SparseMatrixCSC}, J::UniformScaling) = Diagonal(fill(J.ฮป, n)) +# `LinearAlgebra` sizes the `UniformScaling`s and concatenates the result, which is sparse +# only if an argument is. When none is, the first array or number is made sparse, with the +# `Int` indices that the index type of the result has then anyway. +_withsparsefirst(X::Vararg{Any,N}) where {N} = anysparse(X...) ? X : _sparsefirst(X...) +_sparsefirst() = () +_sparsefirst(J::UniformScaling, X::Vararg{Any,N}) where {N} = (J, _sparsefirst(X...)...) +_sparsefirst(x, X::Vararg{Any,N}) where {N} = (_sparse(x), X...) """ sparse_hcat(A...) @@ -611,9 +615,9 @@ Concatenate along dimension 2. Return a SparseMatrixCSC object. the concatenation with specialized "sparse" matrix types from LinearAlgebra.jl automatically yielded sparse output even in the absence of any SparseArray argument. """ -sparse_hcat(Xin::Union{AbstractVecOrMat,Number}...) = _sparse_cat_t(Val(2), promote_eltype(Xin...), map(_makesparse, Xin)...) +sparse_hcat(Xin::Union{AbstractVecOrMat,Number}...) = _sparse_cat_t(Val(2), promote_eltype(Xin...), Xin...) function sparse_hcat(X::Union{AbstractVecOrMat,UniformScaling,Number}...) - LinearAlgebra._hcat(_sparse(first(X)), map(_makesparse, Base.tail(X))...; array_type = SparseMatrixCSC) + LinearAlgebra._hcat(_withsparsefirst(X...)...; array_type = SparseMatrixCSC) end """ @@ -626,9 +630,9 @@ Concatenate along dimension 1. Return a SparseMatrixCSC object. the concatenation with specialized "sparse" matrix types from LinearAlgebra.jl automatically yielded sparse output even in the absence of any SparseArray argument. """ -sparse_vcat(Xin::Union{AbstractVecOrMat,Number}...) = _sparse_cat_t(Val(1), promote_eltype(Xin...), map(_makesparse, Xin)...) +sparse_vcat(Xin::Union{AbstractVecOrMat,Number}...) = _sparse_cat_t(Val(1), promote_eltype(Xin...), Xin...) function sparse_vcat(X::Union{AbstractVecOrMat,UniformScaling,Number}...) - LinearAlgebra._vcat(_sparse(first(X)), map(_makesparse, Base.tail(X))...; array_type = SparseMatrixCSC) + LinearAlgebra._vcat(_withsparsefirst(X...)...; array_type = SparseMatrixCSC) end """ @@ -643,9 +647,7 @@ arguments to concatenate in each block row. the concatenation with specialized "sparse" matrix types from LinearAlgebra.jl automatically yielded sparse output even in the absence of any SparseArray argument. """ -function sparse_hvcat(rows::Tuple{Vararg{Int}}, Xin::Union{AbstractVecOrMat,Number}...) - hvcat(rows, _sparse(first(Xin)), map(_makesparse, Base.tail(Xin))...) -end +sparse_hvcat(rows::Tuple{Vararg{Int}}, Xin::Union{AbstractVecOrMat,Number}...) = _sparse_hvcat(rows, Xin...) function sparse_hvcat(rows::Tuple{Vararg{Int}}, X::Union{AbstractVecOrMat,UniformScaling,Number}...) - LinearAlgebra._hvcat(rows, _sparse(first(X)), map(_makesparse, Base.tail(X))...; array_type = SparseMatrixCSC) + LinearAlgebra._hvcat(rows, _withsparsefirst(X...)...; array_type = SparseMatrixCSC) end diff --git a/test/concatenation.jl b/test/concatenation.jl index ef1cf303..ddcdeef8 100644 --- a/test/concatenation.jl +++ b/test/concatenation.jl @@ -91,7 +91,7 @@ include("testhelpers.jl") A32 = SparseMatrixCSC{Float64,Int32}(A) v32 = SparseVector{Float64,Int32}(vz) @test [A32 v32; v32' 1.0] isa SparseMatrixCSC{Float64,Int32} - @test [1.0 v32'; v32 A32] isa SparseMatrixCSC{Float64,Int} # a leading number widens + @test [1.0 v32'; v32 A32] isa SparseMatrixCSC{Float64,Int32} end @test_throws DimensionMismatch [A vz; 1 2] @static if COMPREHENSIVE @@ -142,6 +142,16 @@ include("testhelpers.jl") sz33 = spzeros(3, 3) de33 = Matrix(1.0I, 3, 3) @test [se33 de33; sz33 se33] == Array([se33 se33; sz33 se33 ]) + # the index type is promoted over the sparse blocks, not taken from the first + A8, A = spzeros(Float64, Int8, 1, 2), spzeros(200, 2) + @test mismatch(vcat(A8, A, ones(1, 2)), [zeros(201, 2); ones(1, 2)], Ti=Int) === nothing + v8, v = spzeros(Float64, Int8, 1), spzeros(200) + @test mismatch(vcat(v8, v, ones(1)), [zeros(201); 1], Ti=Int) === nothing + # a dense block or `I` adds no index type, wherever it is + @test mismatch(hcat(ones(1), v8), [1.0 0.0], Ti=Int8) === nothing + @test mismatch(hcat(I, A8'), [Matrix(1.0I, 2, 2) zeros(2)], Ti=Int8) === nothing + @test mismatch([A8' ones(2); A8' ones(2)], [zeros(2) ones(2); zeros(2) ones(2)], Ti=Int8) === nothing + @test mismatch(sparse_hcat(ones(2), A8', I), [ones(2) zeros(2) Matrix(1.0I, 2, 2)], Ti=Int8) === nothing end end @@ -393,9 +403,9 @@ end @static if COMPREHENSIVE B = sparse(I, 4, 4) @test mismatch(@inferred(hvcat(rows, S, B, B, S)), [dS I; I dS]; Ti=Int) === nothing - # a dense block converts to `Int` indices + # a dense block adds no index type S32 = SparseMatrixCSC{Float64,Int32}(S) - @test mismatch(@inferred(lit22(S32, A)), [dS A; A A]; Ti=Int) === nothing + @test mismatch(@inferred(lit22(S32, A)), [dS A; A A]; Ti=Int32) === nothing end end @@ -584,11 +594,12 @@ end @test mismatch(sparse_vcat(1, 2), [1, 2]) === nothing @static if COMPREHENSIVE @test sparse_hcat(1, 2)::SparseMatrixCSC == [1 2] - # a leading number widens a narrow index type as it did before - v8 = sparsevec(Int8[127], [2], 127); A8 = sparse(Int8[127], Int8[1], [2], 127, 2) - @test vcat(1, v8)::SparseVector{Int,Int} == vcat(1, Array(v8)) - @test sparse_vcat(1, v8)::SparseVector{Int,Int} == vcat(1, Array(v8)) - @test vcat(1, A8)::SparseMatrixCSC{Int,Int} == vcat(1, Array(A8)) + # a leading number adds no index type, and the size of the result is checked against it + v8 = sparsevec(Int8[126], [2], 126); A8 = sparse(Int8[126], Int8[1], [2], 126, 2) + @test vcat(1, v8)::SparseVector{Int,Int8} == vcat(1, Array(v8)) + @test sparse_vcat(1, v8)::SparseVector{Int,Int8} == vcat(1, Array(v8)) + @test vcat(1, A8)::SparseMatrixCSC{Int,Int8} == vcat(1, Array(A8)) + @test_throws ArgumentError vcat(1, 2, v8) # a vector whose length is not an `Int` stacked on a matrix w8 = sparsevec(Int8[1], [2.0], 2) @test mismatch(vcat(w8, ones(1, 1)), vcat(Array(w8), ones(1, 1)), Ti=Int8) === nothing From aba7095e21d3196b09459423b547ea4be4087c3d Mon Sep 17 00:00:00 2001 From: "Viral B. Shah" Date: Mon, 26 Oct 2026 23:56:39 -0400 Subject: [PATCH 28/30] Multiply sparse matrices by range views of a sparse vector with the sparse kernel (#1135) MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit A sparse matrix times a range view of a sparse vector, `A * view(x, 2:6)`, has no sparse method: whole views have one for each sparse left factor, but a range view falls through to LinearAlgebra's generic matrix-vector product, which writes into `similar(view)` element by element. With a triangular `A` that is quadratic. ```julia using SparseArrays, LinearAlgebra x = sparsevec([2, 4, 5], [1.0, 0.0, 3.0], 7); v = view(x, 2:6) U = UpperTriangular(sparse(1.0I, 5, 5)) U * v ``` **Before** ``` sparsevec([1, 4], [1.0, 3.0], 5) ``` **After**, as `U * view(x[2:6], :)` ``` sparsevec([1, 3, 4], [1.0, 0.0, 3.0], 5) ``` For `w = view(sparsevec([1, n], [1.0, 2.0], n), 1:n)` with `n = 10^5` and `A` the identity plus one entry, the minimum of three runs and the `@allocated` of a warm call: | `f(w)` | main | branch | |--------------------------|------------------|---------------------| | `A * w` | 0.0015 s, 272 B | 1.5e-5 s, 1.6 MiB | | `UpperTriangular(A) * w` | 20 s, 240 B | 0.0005 s, 2.4 MiB | | `Symmetric(A) * w` | 0.0018 s, 336 B | 0.0013 s, 3.9 MiB | The methods that send a whole view to `spmatmulv`, for a `SparseMatrixCSC`, a view of one, a triangular, symmetric or Hermitian wrapper and the eager adjoint of a triangular one, now also take `SparseVectorPartialView`, and `_spmatmulv_operand` gives `spmatmulv` the sparse vector a range view reads as, copied from its stored entries by `_sparsevector`. A whole view is passed on as before. The allocations are those of `spmatmulv`, as for a whole view. This touches dispatch. Potential breakage: the product is that of a whole view of `x[2:6]`, which stores an entry wherever products of stored entries contribute, so it can keep zeros the generic product dropped, and its index type is the promotion of both factors' rather than the vector's. Left out: a dense or `Diagonal` matrix times a sparse vector view, whole or range, returns a sparse vector from the generic product, where the same matrix times `x` returns a `Vector` or keeps the stored zeros; and `A' * w` takes the generic product for whole and range views alike. ๐Ÿค– Generated with [Claude Code](https://claude.com/claude-code) Co-authored-by: Claude Opus 5.5 (cherry picked from commit 832997f473e7b1cb106e9d068d841089da217b10) --- src/matmul.jl | 19 +++++++++++-------- test/matmul.jl | 9 +++++++++ 2 files changed, 20 insertions(+), 8 deletions(-) diff --git a/src/matmul.jl b/src/matmul.jl index 52d8a371..8c8ba86c 100644 --- a/src/matmul.jl +++ b/src/matmul.jl @@ -362,12 +362,15 @@ end # Sparse matrix multiplication as described in [Gustavson, 1978]: # http://dl.acm.org/citation.cfm?id=355796 -# spmatmul handles compressed vectors and whole-column/whole-vector views; other -# AbstractSparseVectors take the generic product. Plain CSC times a compressed vector is -# defined with the sparse vector products below. -*(A::SparseTriangular, B::SparseVectorOrView) = spmatmulv(A, B) -*(A::SparseMatrixCSCView, B::SparseVectorOrView) = spmatmulv(A, B) -*(A::AbstractSparseMatrixCSC, B::Union{SparseColumnView,SparseVectorView}) = spmatmulv(A, B) +# spmatmul handles compressed vectors and whole-column/whole-vector views, and a range view +# through the sparse vector it reads as; other AbstractSparseVectors take the generic +# product. Plain CSC times a compressed vector is defined with the sparse vector products +# below. +_spmatmulv_operand(x::Union{AbstractCompressedVector,SparseColumnView,SparseVectorView}) = x +_spmatmulv_operand(x::SparseVectorPartialView) = _sparsevector(x) +*(A::SparseTriangular, B::Union{SparseVectorOrView,SparseVectorPartialView}) = spmatmulv(A, _spmatmulv_operand(B)) +*(A::SparseMatrixCSCView, B::Union{SparseVectorOrView,SparseVectorPartialView}) = spmatmulv(A, _spmatmulv_operand(B)) +*(A::AbstractSparseMatrixCSC, B::Union{SparseColumnView,SparseVectorPartialView}) = spmatmulv(A, _spmatmulv_operand(B)) """ *(A::AbstractSparseMatrixCSC, B::AbstractSparseMatrixCSC) @@ -405,7 +408,7 @@ julia> A * B *(A::SparseMatrixCSCSymmHerm, B::Union{SparseOrTri,AdjOrTrans{<:Any,<:AbstractSparseMatrixCSC}}) = sparse(A) * B *(A::Union{SparseOrTri,AdjOrTrans{<:Any,<:AbstractSparseMatrixCSC}}, B::SparseMatrixCSCSymmHerm) = A * sparse(B) *(A::SparseMatrixCSCSymmHerm, B::SparseMatrixCSCSymmHerm) = sparse(A) * sparse(B) -*(A::SparseMatrixCSCSymmHerm, x::SparseVectorOrView) = sparse(A) * x +*(A::SparseMatrixCSCSymmHerm, x::Union{SparseVectorOrView,SparseVectorPartialView}) = sparse(A) * x # so is the eagerly made adjoint/transpose of a sparse triangular matrix for W in (:UpperTriangular, :LowerTriangular, :UnitUpperTriangular, :UnitLowerTriangular) @eval _sptriangular(A::$W) = $W(copy(parent(A))) @@ -415,7 +418,7 @@ end *(A::Union{SparseOrTri,AdjOrTrans{<:Any,<:AbstractSparseMatrixCSC},SparseMatrixCSCSymmHerm}, B::SparseAdjOrTransTriangular) = A * _sptriangular(B) *(A::SparseAdjOrTransTriangular, B::SparseAdjOrTransTriangular) = _sptriangular(A) * _sptriangular(B) -*(A::SparseAdjOrTransTriangular, x::SparseVectorOrView) = _sptriangular(A) * x +*(A::SparseAdjOrTransTriangular, x::Union{SparseVectorOrView,SparseVectorPartialView}) = _sptriangular(A) * x (*)(Da::Diagonal, A::Union{SparseMatrixCSCOrView, AdjOrTrans{<:Any,<:AbstractSparseMatrixCSC}}, Db::Diagonal) = Da * (A * Db) function (*)(Da::Diagonal, A::SparseMatrixCSC, Db::Diagonal) diff --git a/test/matmul.jl b/test/matmul.jl index dee8035b..5721a67a 100644 --- a/test/matmul.jl +++ b/test/matmul.jl @@ -972,6 +972,15 @@ end @test all(nonzeros(y) .!= 0.0) @test mismatch(y, Af * xf; Ti=Int, approx=true) === nothing @test mismatch((A * view(x, :))::SparseVector{Float64,Int}, Af * xf; approx=true) === nothing + @static if COMPREHENSIVE + # a range view multiplies as a whole view of the sparse vector it reads as + w = view(x, 3:11) + S = A[:, 3:11] + for M in (S, view(A, :, 3:11), UpperTriangular(S), Symmetric(S), UpperTriangular(S)') + @test which(*, (typeof(M), typeof(w))).module === SparseArrays + @test exact_equal(M * w, M * view(x[3:11], :)) + end + end y = *(transpose(A), x2) @test isa(y, SparseVector{Float64,Int}) From 18c8b0487b8a994ee4b2db355a5aae96396210a3 Mon Sep 17 00:00:00 2001 From: "Viral B. Shah" Date: Mon, 26 Oct 2026 23:57:57 -0400 Subject: [PATCH 29/30] Convert range views of a sparse vector from their stored entries (#1134) MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit A range view of a sparse vector, `view(x, 2:6)`, is converted to a `SparseVector` by reading it as a dense vector, so `SparseVector`, `sparse`, `sparsevec`, `convert` and sparse `vcat`, `hcat` and `hvcat` take O(n), drop stored zeros and give `Int` indices, and `kron` of a range view is dense. ```julia using SparseArrays x = sparsevec(Int32[2, 4, 5], [1.0, 0.0, 3.0], 7); v = view(x, 2:6) ``` **Before** ``` SparseVector(v) = sparsevec([1, 4], [1.0, 3.0], 5) vcat(v, v) = sparsevec([1, 4, 6, 9], [1.0, 3.0, 1.0, 3.0], 10) typeof(kron(v, x)) = Vector{Float64} ``` **After** ``` SparseVector(v) = sparsevec(Int32[1, 3, 4], [1.0, 0.0, 3.0], 5) vcat(v, v) = sparsevec(Int32[1, 3, 4, 6, 8, 9], [1.0, 0.0, 3.0, 1.0, 0.0, 3.0], 10) typeof(kron(v, x)) = SparseVector{Float64, Int32} ``` For `w = view(sparsevec([1, n], [1.0, 2.0], n), 1:n)` with `n = 10^7` and `y = sparsevec([1], [1.0], 2)`, the minimum of three runs and the `@allocated` of a warm call: | `f(w)` | main | branch | |-------------------|--------------------|--------------------| | `SparseVector(w)` | 0.030 s, 19.1 MiB | 2.3e-7 s, 240 B | | `vcat(w, w)` | 0.056 s, 38.2 MiB | 3.8e-7 s, 1.1 KiB | | `kron(w, y)` | 0.067 s, 152.6 MiB | 3.8e-7 s, 976 B | The `SparseVector` copy constructors and `sparsevec`, which #1119 widened to whole views, also take `SparseVectorPartialView`, whose `nonzeroinds` and `nonzeros` give its stored entries relative to the range. `convert` and the concatenations reach them through `_makesparse`. `kron` takes range views in `_SparseKronGroup` and in its vector methods, so they go through the same conversion. The `_sparsevector` helper that sparse `map` uses for views now calls `SparseVector` instead of repeating it. This touches dispatch. Potential breakage: for a range view, these conversions and concatenations keep its stored zeros and the parent's index type, and `kron` of one is sparse rather than dense, as for `x[2:6]`. Products with a range view are in #1135. ๐Ÿค– Generated with [Claude Code](https://claude.com/claude-code) Co-authored-by: Claude Opus 5.5 (cherry picked from commit 3d92185cd737d93a9cf0bea9c825cfb9f8e79ef5) --- src/linalg.jl | 12 +++++++----- src/sparsevector.jl | 14 ++++++-------- test/sparsevector.jl | 5 +++-- 3 files changed, 16 insertions(+), 15 deletions(-) diff --git a/src/linalg.jl b/src/linalg.jl index 57a58bb2..4417f927 100644 --- a/src/linalg.jl +++ b/src/linalg.jl @@ -1196,13 +1196,13 @@ end # Hermitian or triangular wrapper of a sparse vector or matrix or its adjoint or transpose. # Views stay out of the wrappers: each such wrapped view in these signatures makes loading # the package much slower, as Julia compares it with LinearAlgebra's triangular methods. -const _SparseKronGroup = Union{SparseVectorOrView, SparseMatrixCSCSubArray, - AdjOrTrans{<:Any,<:Union{SparseVectorOrView, SparseMatrixCSCSubArray}}, +const _SparseKronGroup = Union{SparseVectorOrView, SparseVectorPartialView, SparseMatrixCSCSubArray, + AdjOrTrans{<:Any,<:Union{SparseVectorOrView, SparseVectorPartialView, SparseMatrixCSCSubArray}}, SparseVecOrMatMaybeAdjOrTrans, HermOrSym{<:Any,<:SparseVecOrMatMaybeAdjOrTrans}, UpperOrLowerTriangular{<:Any,<:SparseVecOrMatMaybeAdjOrTrans}} _kronmatrix(A) = convert(SparseMatrixCSC, A) -_kronmatrix(x::Union{SparseColumnView,SparseVectorView}) = convert(SparseMatrixCSC, SparseVector(x)) +_kronmatrix(x::Union{SparseColumnView,SparseVectorView,SparseVectorPartialView}) = convert(SparseMatrixCSC, SparseVector(x)) # A number scales a view of a sparse matrix into a dense matrix unless the view selects # whole columns. Scaling the sparse matrix it converts to gives what the number gives with # a sparse matrix that is not a view. @@ -1340,8 +1340,10 @@ kron(A::_SparseKronGroup, B::_SparseKronGroup) = kron(_kronmatrix(A), _kronmatrix(B)) kron(A::_SparseKronGroup, B::_DenseKronGroup) = kron(A, sparse(B)) kron(A::_DenseKronGroup, B::_SparseKronGroup) = kron(sparse(A), B) -kron(x::SparseVectorOrView, y::SparseVectorOrView) = kron(SparseVector(x), SparseVector(y)) -kron(A::SparseVectorOrView, B::AdjOrTrans{<:Any,<:SparseVectorOrView}) = A .* B +kron(x::Union{SparseVectorOrView,SparseVectorPartialView}, y::Union{SparseVectorOrView,SparseVectorPartialView}) = + kron(SparseVector(x), SparseVector(y)) +kron(A::Union{SparseVectorOrView,SparseVectorPartialView}, + B::AdjOrTrans{<:Any,<:Union{SparseVectorOrView,SparseVectorPartialView}}) = A .* B # disambiguation kron(A::AbstractCompressedVector, B::AdjOrTrans{<:Any,<:AbstractCompressedVector}) = A .* B kron(a::Number, b::_SparseKronGroup) = a * _kronscaled(b) diff --git a/src/sparsevector.jl b/src/sparsevector.jl index 781eba6e..378c2d59 100644 --- a/src/sparsevector.jl +++ b/src/sparsevector.jl @@ -163,11 +163,9 @@ rowvals(x::SparseVectorOrView) = nonzeroinds(x) indtype(x::SparseColumnView) = indtype(parent(x)) indtype(x::Union{SparseVectorView, SparseVectorPartialView}) = indtype(parent(x)) -# the sparse vector a view reads as, built from its stored entries: `copy` of a column view -# widens the index type to `Int` and `SparseVector(view)` scans every element +# the sparse vector a view reads as _sparsevector(x::AbstractCompressedVector) = x -_sparsevector(x::Union{SparseColumnView,SparseVectorView,SparseVectorPartialView}) = - SparseVector(length(x), convert(Vector{indtype(x)}, nonzeroinds(x)), convert(Vector{eltype(x)}, nonzeros(x))) +_sparsevector(x::Union{SparseColumnView,SparseVectorView,SparseVectorPartialView}) = SparseVector(x) function Base.sizehint!(v::SparseVector, newlen::Integer) @@ -555,7 +553,7 @@ sparsevec(a::AbstractVector{T}) where {T} = SparseVector{T, Int}(a) sparsevec(a::AbstractArray) = sparsevec(vec(a)) sparsevec(a::AbstractSparseMatrixCSC) = a[:] sparsevec(a::AbstractSparseVector) = vec(a) -sparsevec(a::Union{SparseColumnView,SparseVectorView}) = SparseVector(a) +sparsevec(a::Union{SparseColumnView,SparseVectorView,SparseVectorPartialView}) = SparseVector(a) sparse(a::AbstractVector) = sparsevec(a) function _dense2indval!(nzind::Vector{Ti}, nzval::Vector{Tv}, s::AbstractArray) where {Tv,Ti} @@ -596,9 +594,9 @@ SparseVector{Tv}(s::AbstractVector) where {Tv} = SparseVector{Tv,Int}(s) SparseVector(s::AbstractVector{Tv}) where {Tv} = SparseVector{Tv,Int}(s) # copy-constructors -SparseVector(s::SparseVectorOrView{Tv,Ti}) where {Tv,Ti} = SparseVector{Tv,Ti}(s) -SparseVector{Tv}(s::SparseVectorOrView{<:Any,Ti}) where {Tv,Ti} = SparseVector{Tv,Ti}(s) -function SparseVector{Tv,Ti}(s::SparseVectorOrView) where {Tv,Ti} +SparseVector(s::Union{SparseVectorOrView{Tv,Ti},SparseVectorPartialView{Tv,Ti}}) where {Tv,Ti} = SparseVector{Tv,Ti}(s) +SparseVector{Tv}(s::Union{SparseVectorOrView{<:Any,Ti},SparseVectorPartialView{<:Any,Ti}}) where {Tv,Ti} = SparseVector{Tv,Ti}(s) +function SparseVector{Tv,Ti}(s::Union{SparseVectorOrView,SparseVectorPartialView}) where {Tv,Ti} copyind = Vector{Ti}(nonzeroinds(s)) copynz = Vector{Tv}(nonzeros(s)) SparseVector{Tv,Ti}(length(s), copyind, copynz) diff --git a/test/sparsevector.jl b/test/sparsevector.jl index deaead50..7c1be97e 100644 --- a/test/sparsevector.jl +++ b/test/sparsevector.jl @@ -1834,10 +1834,11 @@ end # `count` calls `f` on the stored entries and on one zero calls = Ref(0) @test count(z -> (calls[] += 1; !iszero(z)), view(SparseVector{ComplexF64,Int32}(y), :)) == 2 && calls[] == 3 - # a whole view converts with the parent's index type and stored zeros + # a view converts with the parent's index type and stored zeros z = SparseVector{ComplexF64,Int32}(x) C = SparseMatrixCSC{ComplexF64,Int32}(A) - for (v, u) in ((view(z, :), z), (view(C, :, 1), C[:, 1])), f in (SparseVector, sparse, w -> vcat(w, w), w -> kron(w, z)) + for (v, u) in ((view(z, :), z), (view(C, :, 1), C[:, 1]), (view(z, 2:6), z[2:6])), + f in (SparseVector, sparse, w -> vcat(w, w), w -> kron(w, z), w -> kron(w, C), w -> kron(w, w')) @test exact_equal(f(v), f(u)) end end From 42a8f69af2058c0275b8bf740b784fd176708ee7 Mon Sep 17 00:00:00 2001 From: "Viral B. Shah" Date: Tue, 27 Oct 2026 00:49:38 -0400 Subject: [PATCH 30/30] Recover the first-call compile time that #1128 added to `sprand` and `widelength` (#1137) MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit The TTFX check on #1136 flags the 1.14 backports: the cold run of `sprand(200, 200, 0.05); kron(S, S); S * rand(200)` went from 633 to 687 ms and the warm one from 274 to 289 ms. All of it comes from #1128. That PR raised the methods compiled by this workload from 465 to 514, and none of the other 28 backports changed it measurably. **Before** (`main`) ```julia julia> sprand(Xoshiro(1), 2^62, 4, 1e-17) 4611686018427387904ร—4 SparseMatrixCSC{Float64, Int64} with 169 stored entries ``` **After** ```julia julia> sprand(Xoshiro(1), 2^62, 4, 1e-17) ERROR: OverflowError: a 4611686018427387904ร—4 sparse matrix has more entries than `sprand` can draw from ``` First run of the TTFX workload in a fresh process, median of 16 runs interleaved, Julia 1.15.0-DEV.142 on macOS arm64. A cold run starts with an empty JIT object cache, as the check does: | | cold run | warm run | methods compiled | |---|---|---|---| | before #1128 (4c87861) | 282.8 ms | 116.3 ms | 465 | | `main` | 302.0 ms | 121.6 ms | 514 | | this PR | 282.9 ms | 116.2 ms | 450 | `widelength(x) = prod(widen, size(x))` went through `mapreduce`, `MappingRF` and `ComposedFunction`, about 35 methods on the first call of `kron`, `sprand` or a broadcast. It is now `*(map(widen, size(x))...)`: the same `Int128` product, from one call to `*`. `sizehint!` converts that bound to `Int`, which Base's `sizehint!` does anyway, so it is not compiled again for an `Int128` size. Inference and codegen of `sprand` covered the branch that draws column by column for a matrix of 2^63 or more entries, even though no realistic call takes it. That branch is removed: `sprand` throws an `OverflowError` for such a size. Before #1128 it returned a matrix with no stored entries. Should be backported to 1.14 so that #1136 picks it up. Written by Claude Code. ๐Ÿค– Generated with [Claude Code](https://claude.com/claude-code) Co-authored-by: Viral B. Shah Co-authored-by: Claude Opus 5.5 (cherry picked from commit 9bfad2792a57f2199bfc9024c6ef51d50259f454) --- src/SparseArrays.jl | 2 +- src/constructors.jl | 28 ++++------------------------ src/sparsematrix.jl | 4 ++-- test/issues.jl | 2 +- 4 files changed, 8 insertions(+), 28 deletions(-) diff --git a/src/SparseArrays.jl b/src/SparseArrays.jl index 4f319880..20110d14 100644 --- a/src/SparseArrays.jl +++ b/src/SparseArrays.jl @@ -382,7 +382,7 @@ julia> collect(SparseArrays.iternz(x)) """ function iternz end -widelength(x::AbstractSparseArray) = prod(widen, size(x)) +widelength(x::AbstractSparseArray) = *(map(widen, size(x))...) # Linear indexing compares an index with `length(A)`, which wraps for 2^63 entries or more _checklinearindexing(A::AbstractSparseMatrixCSC) = widelength(A) <= typemax(Int) || throw(OverflowError(string("a ", size(A, 1), "ร—", size(A, 2), " sparse matrix has more ", diff --git a/src/constructors.jl b/src/constructors.jl index bf3beb05..9eebafe4 100644 --- a/src/constructors.jl +++ b/src/constructors.jl @@ -714,43 +714,23 @@ julia> sprand(Float64, 3, 0.75) """ function sprand(r::AbstractRNG, m::Integer, n::Integer, density::AbstractFloat, rfn::Function, ::Type{T}) where T m, n = Int(m), Int(n) - if _sprandoverflows(m, n, density) - colptr, rowval, V = _sprand_bycolumn(r, m, n, density, rfn) - return SparseMatrixCSC(m, n, colptr, rowval, convert(Vector{T}, V)) - end I, V = _sprandvalues(r, m, n, density, rfn) return sparse_sortedlinearindices!(I, convert(Vector{T}, V), m, n) end function sprand(r::AbstractRNG, m::Integer, n::Integer, density::AbstractFloat, rfn::Function) m, n = Int(m), Int(n) - if _sprandoverflows(m, n, density) - colptr, rowval, V = _sprand_bycolumn(r, m, n, density, rfn) - return SparseMatrixCSC(m, n, colptr, rowval, convert(Vector{eltype(V)}, V)) - end I, V = _sprandvalues(r, m, n, density, rfn) return sparse_sortedlinearindices!(I, convert(Vector{eltype(V)}, V), m, n) end -# the arguments are checked first, since negative dimensions can overflow too -function _sprandoverflows(m::Int, n::Int, density::AbstractFloat) +function _sprandvalues(r::AbstractRNG, m::Int, n::Int, density::AbstractFloat, rfn::F) where {F} (m < 0 || n < 0) && throw(ArgumentError("invalid Array dimensions")) 0 <= density <= 1 || throw(ArgumentError("$density not in [0,1]")) - return Base.mul_with_overflow(m, n)[2] -end -function _sprandvalues(r::AbstractRNG, m::Int, n::Int, density::AbstractFloat, rfn::F) where {F} + # the entries are drawn as linear indices, which wrap for 2^63 entries or more + Base.mul_with_overflow(m, n)[2] && throw(OverflowError(string("a ", m, "ร—", n, + " sparse matrix has more entries than `sprand` can draw from"))) I = randsubseq(r, 1:(m*n), density) return I, rfn(r, length(I)) end -# an `Int` cannot index the `m*n` entries, so the stored rows are drawn column by column -function _sprand_bycolumn(r::AbstractRNG, m::Int, n::Int, density::AbstractFloat, rfn::F) where {F} - colptr = Vector{Int}(undef, n + 1) - colptr[1] = 1 - rowval, rows = Int[], Int[] - for col in 1:n - append!(rowval, randsubseq!(r, rows, 1:m, density)) - colptr[col + 1] = length(rowval) + 1 - end - return colptr, rowval, rfn(r, length(rowval)) -end sprand(m::Integer, n::Integer, density::AbstractFloat, rfn::Function) = sprand(default_rng(), m, n, density, (r, i) -> rfn(i)) diff --git a/src/sparsematrix.jl b/src/sparsematrix.jl index 010e0dee..02a782af 100644 --- a/src/sparsematrix.jl +++ b/src/sparsematrix.jl @@ -319,7 +319,7 @@ function _storedinds(S::SparseMatrixCSCColumnSubset) end return inds end -widelength(S::SparseMatrixCSCColumnSubset) = prod(widen, size(S)) +widelength(S::SparseMatrixCSCColumnSubset) = *(map(widen, size(S))...) widelength(x::Union{SparseColumnView,SparseVectorView}) = length(x) # `A` itself where it exposes compressed storage through `nzrange`, `rowvals` and `nonzeros`, # otherwise its sparse copy @@ -984,7 +984,7 @@ similar(S::AbstractSparseMatrixCSC, ::Type{TvNew}, ::Type{TiNew}, m::Integer, n: _sparsesimilar(S, TvNew, TiNew, (m, n)) function Base.sizehint!(S::SparseMatrixCSC, n::Integer) - nhint = min(n, widelength(S)) + nhint = Int(min(n, widelength(S))) sizehint!(getrowval(S), nhint) sizehint!(getnzval(S), nhint) return S diff --git a/test/issues.jl b/test/issues.jl index 647e16e6..b765ce0f 100644 --- a/test/issues.jl +++ b/test/issues.jl @@ -979,7 +979,7 @@ end @test_throws OverflowError B[[1]] = [1.0] @test_throws OverflowError B[1] = 1.0 @test nnz(B) == 0 - @test nnz(sprand(Xoshiro(1), 2^62, 4, 1e-17)) > 0 + @test_throws OverflowError sprand(Xoshiro(1), 2^62, 4, 1e-17) end if isdefined(Docs, :undocumented_names) # new in Julia 1.11