diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index 7bd483e..12efa78 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -29,6 +29,7 @@ jobs: - 'min' - 'lts' - '1' + - 'nightly' os: - ubuntu-latest - macOS-latest diff --git a/src/BlockSkylineMatrix.jl b/src/BlockSkylineMatrix.jl index ccbc7f2..6c2f42f 100644 --- a/src/BlockSkylineMatrix.jl +++ b/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) diff --git a/src/linalg.jl b/src/linalg.jl index b018ae2..9d7247a 100644 --- a/src/linalg.jl +++ b/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] diff --git a/test/runtests.jl b/test/runtests.jl index a017f93..b30af9a 100644 --- a/test/runtests.jl +++ b/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 diff --git a/test/test_blockskyline.jl b/test/test_blockskyline.jl index 5385648..79a3e1c 100644 --- a/test/test_blockskyline.jl +++ b/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 diff --git a/test/test_linalg.jl b/test/test_linalg.jl index 92f2380..09c6bef 100644 --- a/test/test_linalg.jl +++ b/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 diff --git a/test/test_triblockbanded.jl b/test/test_triblockbanded.jl index efda3a7..d40c43e 100644 --- a/test/test_triblockbanded.jl +++ b/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)) ==