Repository navigation
Backports for release 1.14 - #1136
Merged
Merged
Conversation
…roadcast such as `A .+ 1` (#1105) 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 <noreply@anthropic.com> (cherry picked from commit ed2edcc)
…ractSparseMatrixCSC` (#1011) 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 <ViralBShah@users.noreply.github.com> Co-authored-by: Claude Opus 5.5 <noreply@anthropic.com> (cherry picked from commit 48204b2)
…to a sparse vector (#1107) `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 <noreply@anthropic.com> (cherry picked from commit b6282ef)
`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 <noreply@anthropic.com>
(cherry picked from commit 71e5eed)
…index type (#1110) `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 <noreply@anthropic.com> (cherry picked from commit 8f7bd49)
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 <noreply@anthropic.com>
(cherry picked from commit 353601e)
… `F` (#1114) 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 <noreply@anthropic.com> (cherry picked from commit a5211ce)
) `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 <noreply@anthropic.com> (cherry picked from commit 8036bf3)
… writing (#1115) `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 <noreply@anthropic.com> (cherry picked from commit dc082ab)
… CHOLMOD cannot take (#1117) `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 <noreply@anthropic.com> (cherry picked from commit 31a10e8)
…oint (#1118) `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 <noreply@anthropic.com> (cherry picked from commit e4a42b4)
…nction (#1109) `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 <noreply@anthropic.com> (cherry picked from commit daa0ac0)
…ew (#1113) `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 <noreply@anthropic.com> (cherry picked from commit 1a56603)
…` before refactorizing (#1122) `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 <noreply@anthropic.com> (cherry picked from commit 1a0d484)
…hortening them (#1123) 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 <noreply@anthropic.com> (cherry picked from commit b86329c)
…opzeros!` throws (#1120) 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 <noreply@anthropic.com> (cherry picked from commit bfa76ee)
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 <noreply@anthropic.com>
(cherry picked from commit 5d9ed81)
…ot convert to its eltype (#1125) `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 <noreply@anthropic.com> (cherry picked from commit 8b27c32)
…DLt fails (#1116) `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 <noreply@anthropic.com> (cherry picked from commit 6ab13ee)
…, `ldiv!` and `rdiv!` (#1121) 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 <noreply@anthropic.com> (cherry picked from commit 66a7b5c)
…ing (#1126) `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 <noreply@anthropic.com> (cherry picked from commit 4c87861)
…` and linear indexing (#1128) `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 <noreply@anthropic.com> (cherry picked from commit 2bef1fd)
#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 <noreply@anthropic.com> (cherry picked from commit 5b98f3c)
…gular matrices (#1131) 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 <noreply@anthropic.com> (cherry picked from commit eb47441)
…for any index type (#1119) 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 <noreply@anthropic.com> (cherry picked from commit 61f75e0)
`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 <noreply@anthropic.com>
(cherry picked from commit aee932a)
…s sparse arguments only (#1127) 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 <noreply@anthropic.com> (cherry picked from commit b575c39)
…parse kernel (#1135) 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 <noreply@anthropic.com> (cherry picked from commit 832997f)
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 <noreply@anthropic.com>
(cherry picked from commit 3d92185)
TTFX: No robust change.
DetailsJulia 1.15.0-DEV.167 · Linux X64 · 4 blocks (ABBA)
Dependents (median change over blocks)
Dependents per block
A change is robust when every block clears the threshold (precompile 8%, load (cold) 10%, run (cold) 15%, load+run (warm) 15%), the samples of the two commits do not overlap, and the medians differ by at least the floor (precompile 0.25 s, load (cold) 0.1 s, run (cold) 0.05 s, load+run (warm) 0.05 s). Load and run are also judged with the GC off, and either counts. Over all packages whose base value is at least the floor, each metric's geometric mean must clear a threshold in every block (precompile 2%, load (cold) 2%, run (cold) 3%, load+run (warm) 3%). |
Codecov Report❌ Patch coverage is
Additional details and impacted files@@ Coverage Diff @@
## release-1.14 #1136 +/- ##
================================================
- Coverage 94.82% 94.67% -0.15%
================================================
Files 16 16
Lines 10253 10422 +169
================================================
+ Hits 9722 9867 +145
- Misses 531 555 +24 ☔ View full report in Codecov by Harness. 🚀 New features to boost your workflow:
|
ViralBShah
added a commit
that referenced
this pull request
Oct 27, 2026
…`widelength` (#1137) 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 <ViralBShah@users.noreply.github.com> Co-authored-by: Claude Opus 5.5 <noreply@anthropic.com>
…`widelength` (#1137) 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 <ViralBShah@users.noreply.github.com> Co-authored-by: Claude Opus 5.5 <noreply@anthropic.com> (cherry picked from commit 9bfad27)
ViralBShah
force-pushed
the
backports-release-1.14
branch
from
October 27, 2026 04:51
00c2e37 to
42a8f69
Compare
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
Backports of the merged PRs labelled
backport 1.14, cherry-picked with-xin the order they were merged onmain.release-1.14starts at b0aa455 (#1104), the commit JuliaLang/julia'srelease-1.14pins, and every PR merged onmainsince then is labelled, so this branch's tree is identical tomain's at 9bfad27.A .+ 1#1105rowvalsandnonzeroson a subtype ofAbstractSparseMatrixCSC#1011copyto!writes them into a sparse vector #1107CHOLMOD.Sparsehas one column insparsevec#1108vcat,hcatandrepeatfits the index type #1110sparse!before writing #1111lu!(F::UmfpackLU, S)writes toF#1114UmfpackLUagainst aComplexF32vector inldiv!#1112circshift!of a sparse matrix can hold the result before writing #1115ArgumentErrorfromcholeskyandldltfor an index type CHOLMOD cannot take #1117diffof a sparse view, transpose or adjoint #1118Tand the documentedrfnforms insprandwith a value function #1109n x 1matrix fromSparseMatrixCSCof a sparse vector view #1113ldlt!andcholesky!of aCHOLMOD.Factorbefore refactorizing #1122SparseMatrixCSCconstructor before shortening them #1123fkeep!,droptol!ordropzeros!throws #1120\andqr(A) \ bthe eltype of dense\#1124bwhen a triangularldiv!solution does not convert to its eltype #1125luinfactorizeof a sparse Hermitian matrix whose LDLt fails #1116mul!,ldiv!andrdiv!#1121Inffrom sparsecondof a singular matrix instead of throwing #1126sprandand linear indexing #1128condof a singular matrix returnsInf#1130condfor singular and nonsingular matrices #1131spdiagmfrom the stored entries of a sparse vector view #1133sprandandwidelength#1137#1137 removes the first-call compile time that #1128 added, which the TTFX check flagged on this PR.
#1132 merged into the branch of #1131 and is part of its commit.
Not bug fixes: #1011 adds a deprecation, #1109 changes the
sprandvalue stream for a value function, #1113 changes the result shape ofSparseMatrixCSCof a vector view, and #1119, #1133, #1134 and #1135 add methods for views. They are included because they are labelled and 1.14 is not released yet.🤖 Generated with Claude Code