From a4d33f8065a06f741ec8fcf12c23edafea46e3a4 Mon Sep 17 00:00:00 2001 From: Sheehan Olver Date: Tue, 13 Oct 2026 16:18:42 +0100 Subject: [PATCH 1/4] support sparse for LayoutArray and friends --- Project.toml | 2 +- ext/ArrayLayoutsSparseArraysExt.jl | 52 ++++++++++++++++++++++++++++-- src/ArrayLayouts.jl | 10 ++++++ test/test_layoutarray.jl | 50 ++++++++++++++++++++++++++++ 4 files changed, 111 insertions(+), 3 deletions(-) diff --git a/Project.toml b/Project.toml index 8e609277..6f723e62 100644 --- a/Project.toml +++ b/Project.toml @@ -1,6 +1,6 @@ name = "ArrayLayouts" uuid = "4c555306-a7a7-4459-81d9-ec55ddd5c99a" -version = "1.15.1" +version = "1.16" authors = ["Sheehan Olver "] [deps] diff --git a/ext/ArrayLayoutsSparseArraysExt.jl b/ext/ArrayLayoutsSparseArraysExt.jl index 32b543e7..88212e27 100644 --- a/ext/ArrayLayoutsSparseArraysExt.jl +++ b/ext/ArrayLayoutsSparseArraysExt.jl @@ -1,15 +1,63 @@ module ArrayLayoutsSparseArraysExt using ArrayLayouts -using ArrayLayouts: copyto!_layout, equals_layout +using ArrayLayouts: copyto!_layout, equals_layout, sparse_layout, LayoutVecOrMat, UpperOrLowerTriangular using SparseArrays # Specifying the full namespace is necessary because of https://github.com/JuliaLang/julia/issues/48533 # See https://github.com/JuliaStats/LogExpFunctions.jl/pull/63 import ArrayLayouts.LinearAlgebra -using ArrayLayouts.LinearAlgebra: Adjoint, Transpose, Symmetric, Hermitian, Diagonal, UpperTriangular, UnitUpperTriangular, LowerTriangular, UnitLowerTriangular +using ArrayLayouts.LinearAlgebra: Adjoint, Transpose, Symmetric, Hermitian, Diagonal, UpperTriangular, UnitUpperTriangular, LowerTriangular, UnitLowerTriangular, + AdjOrTrans, HermOrSym import Base: copyto!, == +# use colsupport to deduce the sparsity pattern +ArrayLayouts.sparse_layout(_, A) = sparse_colsupport(A) + +function sparse_colsupport(A::AbstractMatrix{T}) where T + Base.require_one_based_indexing(A) + m, n = size(A) + colptr = Vector{Int}(undef, n+1) + rowval = Int[] + nzval = T[] + colptr[1] = 1 + for j = 1:n + for k in colsupport(A, j) + v = A[k, j] + if !iszero(v) + push!(rowval, k) + push!(nzval, v) + end + end + colptr[j+1] = length(rowval) + 1 + end + SparseMatrixCSC(m, n, colptr, rowval, nzval) +end + +function sparse_colsupport(a::AbstractVector{T}) where T + Base.require_one_based_indexing(a) + nzind = Int[] + nzval = T[] + for k in colsupport(a, 1) + v = a[k] + if !iszero(v) + push!(nzind, k) + push!(nzval, v) + end + end + SparseVector(length(a), nzind, nzval) +end + +const LayoutSparseVectors = Union{LayoutVector, SubArray{<:Any,1,<:LayoutArray}} +const LayoutSparseMatrixParents = Union{LayoutMatrix, SubArray{<:Any,2,<:LayoutArray}, AdjOrTrans{<:Any,<:LayoutVecOrMat}} +const LayoutSparseMatrices = Union{LayoutSparseMatrixParents, + HermOrSym{<:Any,<:LayoutMatrix}, + UpperOrLowerTriangular{<:Any,<:LayoutSparseMatrixParents}, + Diagonal{<:Any,<:LayoutVector}} + +SparseArrays.sparse(a::LayoutSparseVectors) = sparse_layout(a) +SparseArrays.sparse(A::LayoutSparseMatrices) = sparse_layout(A) + # ambiguity from sparsematrix.jl copyto!(dest::LayoutMatrix, src::SparseArrays.AbstractSparseMatrixCSC) = copyto!_layout(dest, src) diff --git a/src/ArrayLayouts.jl b/src/ArrayLayouts.jl index 353ab9d6..e6eddf0f 100644 --- a/src/ArrayLayouts.jl +++ b/src/ArrayLayouts.jl @@ -329,6 +329,16 @@ equals_layout(A, B) = equals_layout(MemoryLayout(A), MemoryLayout(B), A, B) ==(A::SubArray{<:Any,<:Any,<:LayoutArray}, B::AbstractArray) = equals_layout(A, B) ==(A::AbstractArray, B::SubArray{<:Any,<:Any,<:LayoutArray}) = equals_layout(A, B) +""" + sparse_layout(layout, A) + +Layout-aware hook for `SparseArrays.sparse(A)`, called with `layout = MemoryLayout(A)`. +Overload this to specialise conversion to a sparse array on memory layouts. The default +(defined in the SparseArrays extension) uses `colsupport` to deduce the sparsity pattern. +""" +function sparse_layout end +sparse_layout(A) = sparse_layout(MemoryLayout(A), A) + _fill_lmul!(β, A::AbstractArray) = iszero(β) ? zero!(A) : lmul!(β, A) diff --git a/test/test_layoutarray.jl b/test/test_layoutarray.jl index 3449cb1d..ffa2cf6d 100644 --- a/test/test_layoutarray.jl +++ b/test/test_layoutarray.jl @@ -834,6 +834,56 @@ ArrayLayouts.equals_layout(::MyEqualsLayout, _, A, B) = false end end +struct MySparseLayout <: MemoryLayout end +struct MyBidiagonalMatrix <: LayoutMatrix{Float64} end +Base.size(::MyBidiagonalMatrix) = (4,4) +Base.getindex(::MyBidiagonalMatrix, k::Int, j::Int) = k == j ? Float64(k) : k+1 == j ? 1.0 : error("outside colsupport") +MemoryLayout(::Type{MyBidiagonalMatrix}) = MySparseLayout() +ArrayLayouts.colsupport(::MySparseLayout, A, j) = max(1,minimum(j)-1):min(size(A,1),maximum(j)) +struct MyZeroVector <: LayoutVector{Float64} end +Base.size(::MyZeroVector) = (3,) +Base.getindex(::MyZeroVector, k::Int) = error("outside colsupport") +MemoryLayout(::Type{MyZeroVector}) = ZerosLayout() +struct MySparseOverloadLayout <: MemoryLayout end +struct MySparseOverloadVector <: LayoutVector{Float64} end +Base.size(::MySparseOverloadVector) = (3,) +Base.getindex(::MySparseOverloadVector, k::Int) = Float64(k) +MemoryLayout(::Type{MySparseOverloadVector}) = MySparseOverloadLayout() +ArrayLayouts.sparse_layout(::MySparseOverloadLayout, a) = sparsevec([2], [-1.0], 3) + +@testset "sparse" begin + A = [1.0 0 2; 0 3 4; 5 0 0] + M = MyMatrix(A) + for B in (M, view(M, 1:2, 2:3), M', transpose(M), Symmetric(M), Hermitian(M, :L), + UpperTriangular(M), UnitLowerTriangular(M), LowerTriangular(M'), UpperTriangular(view(M, :, :))) + S = sparse(B) + @test S isa SparseMatrixCSC{Float64,Int} + @test S == B + @test nnz(S) == count(!iszero, B) + end + v = MyVector([1.0, 0, 2]) + for b in (v, view(v, 2:3), view(M, :, 2)) + s = sparse(b) + @test s isa SparseVector{Float64,Int} + @test s == b + @test nnz(s) == count(!iszero, b) + end + D = sparse(Diagonal(v)) + @test D isa SparseMatrixCSC{Float64,Int} + @test D == Diagonal([1.0, 0, 2]) + @test nnz(D) == 2 + @test sparse(v') == v' + + B = MyBidiagonalMatrix() + S = sparse(B) + @test S == Bidiagonal([1.0:4;], ones(3), :U) + @test nnz(S) == 7 + s = sparse(MyZeroVector()) + @test s == zeros(3) + @test nnz(s) == 0 + @test sparse(MySparseOverloadVector()) == [0, -1.0, 0] +end + @testset "* for infinite layouts" begin tup = InfSymTridiagonal(), InfTridiagonal(), InfBidiagonal('U'), InfBidiagonal('L'), From d61e68f28b93fb66282f42e2ababda8a38e8c1ed Mon Sep 17 00:00:00 2001 From: Sheehan Olver Date: Thu, 15 Oct 2026 15:19:41 +0100 Subject: [PATCH 2/4] v1.17 --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index de648e12..a81e20e6 100644 --- a/Project.toml +++ b/Project.toml @@ -1,6 +1,6 @@ name = "ArrayLayouts" uuid = "4c555306-a7a7-4459-81d9-ec55ddd5c99a" -version = "1.16.1" +version = "1.17" authors = ["Sheehan Olver "] [deps] From 6867c5279c6be1409ea0d232022770d045def7f6 Mon Sep 17 00:00:00 2001 From: Sheehan Olver Date: Sat, 24 Oct 2026 18:45:20 -0400 Subject: [PATCH 3/4] Fix _zeroeltype for conrtinuum arrays --- src/mul.jl | 5 +++-- test/test_muladd.jl | 3 +++ 2 files changed, 6 insertions(+), 2 deletions(-) diff --git a/src/mul.jl b/src/mul.jl index c4cd7680..a76e10f2 100644 --- a/src/mul.jl +++ b/src/mul.jl @@ -33,8 +33,9 @@ axes(M::Mul) = _mul_axes(axes(M.A), axes(M.B)) zeroeltype(M) = _zeroeltype(eltype(M), M) # allow special casing where we know more about zero _zeroeltype(::Type{T}, M) where T = zero(T) -# zero(T) is not defined for arrays like Vector{Float64} so use the entries to determine the size -_zeroeltype(::Type{T}, M::Mul) where T<:AbstractArray = convert(T, zero(first(M.A) * first(M.B))) +# zero(T) is not defined for arrays like Vector{Float64} so use the entries to determine the size. +# Prefer zero(T) when it exists (e.g. SVector) since first may not be defined (e.g. quasi-arrays) +_zeroeltype(::Type{T}, M::Mul) where T<:AbstractArray = applicable(zero, T) ? zero(T) : convert(T, zero(first(M.A) * first(M.B))) zeroeltype(M::Mul{<:Any,<:Any,<:SubArray}) = zeroeltype(Mul(parent(M.A), M.B)) function _getindex(::Type{Tuple{AA}}, M::Mul, (k,)::Tuple{AA}) where AA diff --git a/test/test_muladd.jl b/test/test_muladd.jl index 84a22c68..d2896723 100644 --- a/test/test_muladd.jl +++ b/test/test_muladd.jl @@ -99,6 +99,9 @@ Random.seed!(0) @test ArrayLayouts.Mul(X, reshape(c,2,1))[1,1] == X[1] + 2X[2] v = reshape([[1.0, 2.0], [3.0, 4.0]], 1, 2) @test ArrayLayouts.Mul(v, c)[1] == [7.0, 10.0] + # eltypes with zero(T) should not need any entries + S = Matrix{SVector{2,Float64}}(undef, 1, 0) + @test ArrayLayouts.Mul(S, Float64[])[1] === zero(SVector{2,Float64}) end @testset "Diagonal Fill" begin From 0e68ab29aefe0cde07aaec7d38026137f35bf8e5 Mon Sep 17 00:00:00 2001 From: Sheehan Olver Date: Sun, 25 Oct 2026 09:33:02 -0400 Subject: [PATCH 4/4] Update Project.toml --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index 0b47929d..7b9bea05 100644 --- a/Project.toml +++ b/Project.toml @@ -1,6 +1,6 @@ name = "ArrayLayouts" uuid = "4c555306-a7a7-4459-81d9-ec55ddd5c99a" -version = "1.17" +version = "1.16" authors = ["Sheehan Olver "] [deps]