Add cconvert for pointers to BlockSkylineMatrix blocks

.github/workflows/ci.yml

@@ -29,6 +29,7 @@ jobs:
           - 'min'
           - 'lts'
           - '1'
+          - 'nightly'
         os:
           - ubuntu-latest
           - macOS-latest

src/BlockSkylineMatrix.jl

@@ -6,7 +6,8 @@ checkbandwidths(N, M, l::AbstractVector{Int}, u::AbstractVector{Int}) =
 function bb_blockstarts(ax, l::AbstractVector{Int}, u::AbstractVector{Int})
     N,M = blocksize.(ax,1)
     L,U = maximum(l), maximum(u)
-    b_start = BandedMatrix{Int}(undef, (N, M), (L, U))
+    # zero-filled so in-band slots for blocks outside a column's own band read as 0
+    b_start = BandedMatrix(Zeros{Int}(N, M), (L, U))
     -L > U && return b_start
 
     checkbandwidths(N, M, l, u)
@@ -455,6 +456,14 @@ _parent_blocks(V::BlockBandedBlock)::Tuple{Int,Int} =
 MemoryLayout(::Type{<:BlockBandedBlock}) = ColumnMajor()
 Base.elsize(::Type{<:BlockSkylineMatrix{T,R}}) where {T,R} = Base.elsize(R)
 
+function Base.cconvert(::Type{Ptr{T}}, V::BlockBandedBlock{T}) where T
+    A = parent(V)
+    K,J = _parent_blocks(V)
+    Block(K) ∈ blockcolsupport(A, Block(J)) || throw(ArgumentError("Pointer is only defined when inside blockcolsupport"))
+    # a range, not a single index, so empty blocks starting past the end of data are allowed
+    Base.cconvert(Ptr{T}, view(A.data, blockstart(A,K,J):lastindex(A.data)))
+end
+
 function Base.unsafe_convert(::Type{Ptr{T}}, V::BlockBandedBlock{T}) where T
     A = parent(V)
     K,J = _parent_blocks(V)

src/linalg.jl

@@ -106,6 +106,13 @@ strides(V::SubBlockSkylineMatrix{<:Any,LL,UU,<:Union{BlockRange1,Block1},Block1}
     (1,parent(V).block_sizes.block_strides[Int(parentindices(V)[2].block)])
 
 
+function Base.cconvert(::Type{Ptr{T}}, V::SubBlockSkylineMatrix{T,LL,UU,<:Union{BlockRange1,Block1},Block1}) where {T,LL,UU}
+    A = parent(V)
+    KR = parentindices(V)[1].block
+    J = parentindices(V)[2].block
+    Base.cconvert(Ptr{T}, view(A, first(KR), J))
+end
+
 function unsafe_convert(::Type{Ptr{T}}, V::SubBlockSkylineMatrix{T,LL,UU,<:Union{BlockRange1,Block1},Block1}) where {T,LL,UU}
     A = parent(V)
     JR = parentindices(V)[2]
@@ -117,6 +124,16 @@ end
 strides(V::SubBlockSkylineMatrix{<:Any,LL,UU,<:BlockRange1,<:BlockIndexRange1}) where {LL,UU} =
     (1,parent(V).block_sizes.block_strides[Int(Block(parentindices(V)[2]))])
 
+function Base.cconvert(::Type{Ptr{T}}, V::SubBlockSkylineMatrix{T,LL,UU,<:BlockRange1,<:BlockIndexRange1}) where {T,LL,UU}
+    A = parent(V)
+    JR = parentindices(V)[2]
+    K = first(parentindices(V)[1].block)
+    J = Block(JR)
+    K ∈ blockcolsupport(A, J) || throw(ArgumentError("Pointer is only defined when inside blockcolsupport"))
+    s = blockstart(A, Int(K), Int(J)) + (first(JR.block.indices[1]) - 1)*blockstride(A, Int(J))
+    Base.cconvert(Ptr{T}, view(A.data, s:lastindex(A.data)))
+end
+
 function unsafe_convert(::Type{Ptr{T}}, V::SubBlockSkylineMatrix{T,LL,UU,<:BlockRange1,<:BlockIndexRange1}) where {T,LL,UU}
     A = parent(V)
     JR = parentindices(V)[2]
@@ -127,6 +144,17 @@ function unsafe_convert(::Type{Ptr{T}}, V::SubBlockSkylineMatrix{T,LL,UU,<:Block
     p + sizeof(T)*(JR.block.indices[1][1]-1)*stride(V,2)
 end
 
+function Base.cconvert(::Type{Ptr{T}}, V::SubBlockSkylineMatrix{T,LL,UU,BlockIndexRange1,BlockIndexRange1}) where {T,LL,UU}
+    A = parent(V)
+    K = parentindices(V)[1].block.block
+    kr = parentindices(V)[1].block.indices[1]
+    J = parentindices(V)[2].block.block
+    jr = parentindices(V)[2].block.indices[1]
+    K ∈ blockcolsupport(A, J) || throw(ArgumentError("Pointer is only defined when inside blockcolsupport"))
+    s = blockstart(A, Int(K), Int(J)) + (first(kr) - 1) + (first(jr) - 1)*blockstride(A, Int(J))
+    Base.cconvert(Ptr{T}, view(A.data, s:lastindex(A.data)))
+end
+
 function unsafe_convert(::Type{Ptr{T}}, V::SubBlockSkylineMatrix{T,LL,UU,BlockIndexRange1,BlockIndexRange1}) where {T,LL,UU}
     A = parent(V)
     JR = parentindices(V)[2]

test/runtests.jl

@@ -2,8 +2,47 @@ using BlockBandedMatrices
 using ParallelTestRunner
 
 const init_code = quote
-   using Test
-   using BlockBandedMatrices
+    using Test
+    using BlockBandedMatrices
+
+    function strided_ptr(f, a::AbstractArray{T}) where {T}
+        a_cconv = Base.cconvert(Ptr{T}, a)
+        GC.@preserve a_cconv begin
+            f(Base.unsafe_convert(Ptr{T}, a_cconv))
+        end
+    end
+
+    """
+        check_strided_get(a::AbstractArray{T,N})
+
+    Test that array `a` implements the strided array interface for reading.
+    Checks stride consistency and that `unsafe_load` matches `getindex`.
+    """
+    function check_strided_get(a::AbstractArray{T,N})::Nothing where {T, N}
+        if !isbitstype(eltype(a))
+            error("a doesn't have isbits elements")
+        end
+        # Putting strided_ptr before the loop means that strided_ptr shouldn't error for empty arrays
+        strided_ptr(a) do a_ptr
+            for d in 1:N
+                if stride(a, d) != strides(a)[d]
+                    error("stride(a, d) doesn't equal strides(a)[d] for dimension $(d)")
+                end
+            end
+            for i in CartesianIndices(a)
+                el_ptr = a_ptr
+                for d in 1:N
+                    stride_in_bytes = stride(a, d) * Base.elsize(typeof(a))
+                    first_idx = first(axes(a, d))
+                    el_ptr += (i[d] - first_idx) * stride_in_bytes
+                end
+                if unsafe_load(el_ptr) !== a[i]
+                    error("getindex and unsafe_load mismatch at index $(i)")
+                end
+            end
+        end
+        nothing
+    end
 end
 
 # Start with autodiscovered tests

test/test_blockskyline.jl

@@ -165,6 +165,14 @@ Random.seed!(0)
         s = split(sprint(show, "text/plain", B), '\n')[2]
         @test s == " 1  │  0  0  │  ⋅  ⋅  ⋅  │  ⋅  ⋅  ⋅  ⋅"
     end
+
+    @testset "block_starts" begin
+        # block (3,2) is in the band of block_starts but not stored, so must be 0
+        bs = BlockBandedMatrices.BlockSkylineSizes([2,3,4], [1,2,3], [2,0,0], [0,1,1])
+        @test Matrix(bs.block_starts) == [1 10  0;
+                                          3 12 20;
+                                          6  0 23]
+    end
 end
 
 end # module

test/test_linalg.jl

@@ -10,6 +10,7 @@ using Test
 
 import BandedMatrices: BandError, bandeddata
 import BlockBandedMatrices: _BandedBlockBandedMatrix
+import ..check_strided_get
 
 @testset "lmul!/rmul!" begin
     C = BandedBlockBandedMatrix{Float64}(undef, 1:2,1:2, (1,1), (1,1))
@@ -45,6 +46,7 @@ end
     @test stride(V,2) == 7
     @test unsafe_load(pointer(V)) == 46
     @test unsafe_load(pointer(V) + stride(V,2)*sizeof(Float64)) == 53
+    check_strided_get(V)
 
     x = randn(size(A,2))
     @test A*x == (similar(x) .= MulAdd(A,x)) ≈ Matrix(A)*x

test/test_triblockbanded.jl

@@ -12,6 +12,7 @@ import BlockBandedMatrices: MemoryLayout, TriangularLayout,
                             blockrowstop, blockcolstop, ColumnMajor
 
 import BlockArrays: blockisequal
+import ..check_strided_get
 
 @testset "triangular" begin
     @testset "triangular BandedBlockBandedMatrix mul" begin
@@ -155,6 +156,7 @@ import BlockArrays: blockisequal
         V = view(A, Block.(2:3), Block(3))
         @test unsafe_load(pointer(V)) == A[2,4]
         @test unsafe_load(pointer(V)+sizeof(Float64)*stride(V,2)) == A[2,5]
+        check_strided_get(V)
         @test MemoryLayout(typeof(V)) == ColumnMajor()
 
         @test size(V) == (5,3)
@@ -165,6 +167,7 @@ import BlockArrays: blockisequal
 
         V = view(A, Block.(1:3), Block(3)[2:3])
         @test_throws ArgumentError pointer(V)
+        @test_throws ArgumentError Base.cconvert(Ptr{Float64}, V)
         @test MemoryLayout(typeof(V)) == ColumnMajor()
 
         b = randn(size(A,1))
@@ -173,6 +176,7 @@ import BlockArrays: blockisequal
         V = view(A, Block.(2:3), Block(3)[2:3])
         @test unsafe_load(pointer(V)) == A[2,5]
         @test unsafe_load(pointer(V)+sizeof(Float64)*stride(V,2)) == A[2,6]
+        check_strided_get(V)
         @test MemoryLayout(typeof(V)) == ColumnMajor()
 
         @test size(V) == (5,2)
@@ -182,6 +186,7 @@ import BlockArrays: blockisequal
 
         V = view(A, Block.(1:3), Block(3)[2:3])
         @test_throws ArgumentError pointer(V)
+        @test_throws ArgumentError Base.cconvert(Ptr{Float64}, V)
         @test MemoryLayout(typeof(V)) == ColumnMajor()
 
         b = randn(size(A,1))
@@ -193,6 +198,8 @@ import BlockArrays: blockisequal
         V = view(A, Block(N),  Block(N))
         V_22 = view(A, Block(N)[1:N],  Block(N)[1:N])
         @test unsafe_load(pointer(V_22)) == V_22[1,1] == V[1,1]
+        check_strided_get(V_22)
+        check_strided_get(V)
         @test strides(V_22) == strides(V) == (1,9)
         b = randn(N)
         @test copyto!(similar(b) , MulAdd(V,b)) == copyto!(similar(b) , MulAdd(V_22,b)) ==