Skip to content

Commit 4852432

Browse files
KristofferCKristofferC
committed
improve performance by explicitly instantiating functions for N = 1, 2 (#5)
1 parent 1c3301a commit 4852432

File tree

6 files changed

+237
-40
lines changed

6 files changed

+237
-40
lines changed

‎README.md

Lines changed: 3 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -2,12 +2,12 @@
22

33
[![Build Status](https://travis-ci.org/KristofferC/BlockArrays.jl.svg?branch=master)](https://travis-ci.org/KristofferC/BlockArrays.jl) [![codecov](https://codecov.io/gh/KristofferC/BlockArrays.jl/branch/master/graph/badge.svg)](https://codecov.io/gh/KristofferC/BlockArrays.jl)
44

5-
**Note:** Currently, a quite new build of julia master is needed to use this package.
6-
75
A block array is a partition of an array into blocks or subarrays, see [wikipedia](https://en.wikipedia.org/wiki/Block_matrix) for a good description. This package introduces the abstract type `AbstractBlockArray` for arrays that exhibit this block structure. Currently, two concrete types are implemented that have very similar API but differs in how the blocks are stored in memory. The type `BlockArray` stores each block contiguously while the type `PseudoBlockArray` stores the full matrix contiguously. Which one to use depends on your use case, `BlockArray` supports fast non copying extraction and insertion of blocks while `PseudoBlockArray` supports directly using the underlying full matrix in for example a linear solver. Both these types follow the `AbstractArray` interface and should work in arbitrary dimensions for arbitrary block types as long as the block type itself satisfies the `AbstractArray` interface.
86

97
This README will first provide an overview over the `BlockArray` type and then later discuss the few differences between `BlockArrays` and `PseudoBlockArrays`.
108

9+
**Note:** Currently, a quite new build of julia master is needed to use this package.
10+
1111
### Creating uninitialized `BlockArrays`.
1212

1313
A `BlockArray` can be created with the blocks left uninitialized using the `BlockArray(block_type, block_sizes...)` function.
@@ -150,7 +150,7 @@ The underlying array is accessed with `full` just like for `BlockArray`.
150150
## TODO
151151

152152
- Linear algebra stuff
153-
- Investigate performance
153+
- Investigate performance (for example that bounds check removal work properly)
154154

155155
## Author
156156

‎src/BlockArrays.jl

Lines changed: 12 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -2,19 +2,28 @@ __precompile__()
22

33
module BlockArrays
44

5-
export AbstractBlockArray, BlockArray, PseudoBlockArray, Block, getblock, getblock!, setblock!
5+
export Block, getblock, getblock!, setblock!, nblocks
6+
7+
export AbstractBlockArray, AbstractBlockMatrix, AbstractBlockVector, AbstractBlockVecOrMat
8+
export BlockArray, BlockMatrix, BlockVector, BlockVecOrMat
9+
export PseudoBlockArray, PseudoBlockMatrix, PseudoBlockVector, PseudoBlockVecOrMat
610

711
import Base: @propagate_inbounds
812
using Base.Cartesian
913

1014
abstract AbstractBlockArray{T, N, R} <: AbstractArray{T, N}
11-
1215
typealias AbstractBlockMatrix{T, R} AbstractBlockArray{T, 2, R}
1316
typealias AbstractBlockVector{T, R} AbstractBlockArray{T, 1, R}
1417
typealias AbstractBlockVecOrMat{T, R} Union{AbstractBlockMatrix{T, R}, AbstractBlockVector{T, R}}
1518

16-
nblocks(block_array::AbstractBlockArray) = nblocks(block_array.block_sizes)
19+
"""
20+
nblocks(block_array[, i])
21+
22+
The number of blocks in a block array, optionally in dimension i.
23+
"""
1724
nblocks(block_array::AbstractBlockArray, i::Int) = nblocks(block_array.block_sizes, i)
25+
nblocks(block_array::AbstractBlockArray) = nblocks(block_array.block_sizes)
26+
1827

1928
Base.similar{T,N}(block_array::AbstractBlockArray{T,N}) = similar(block_array, T)
2029
Base.size(arr::AbstractBlockArray) = map(sum, arr.block_sizes.sizes)

‎src/block.jl

Lines changed: 14 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -1,3 +1,6 @@
1+
# Note: Functions surrounded by a comment blocks are there because `Vararg` is sitll allocating.
2+
# When Vararg is fast enough, they can simply be removed
3+
14
"""
25
A `Block` is simply a wrapper around a set of indixes so that it can be used to dispatch on. By
36
indexing a `AbstractBlockArray` with a `Block` the the returned object value be that block.
@@ -8,5 +11,15 @@ end
811

912
Block{N}(n::Vararg{Int, N}) = Block{N}(n)
1013

11-
@propagate_inbounds Base.setindex!{T,N,R}(block_arr::AbstractBlockArray{T,N,R}, v, block::Block{N}) = setblock!(block_arr, v, block.n...)
14+
# vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv #
15+
@propagate_inbounds Base.getindex{T}(block_arr::AbstractBlockArray{T,1}, block::Block{1}) = getblock(block_arr, block.n[1])
16+
@propagate_inbounds Base.getindex{T}(block_arr::AbstractBlockArray{T,2}, block::Block{2}) = getblock(block_arr, block.n[1], block.n[2])
17+
# ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ #
1218
@propagate_inbounds Base.getindex{T,N}(block_arr::AbstractBlockArray{T,N}, block::Block{N}) = getblock(block_arr, block.n...)
19+
20+
# vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv #
21+
@propagate_inbounds Base.setindex!{T}(block_arr::AbstractBlockArray{T,1}, v, block::Block{1}) = setblock!(block_arr, v, block.n[1])
22+
@propagate_inbounds Base.setindex!{T}(block_arr::AbstractBlockArray{T,2}, v, block::Block{2}) = setblock!(block_arr, v, block.n[1], block.n[2])
23+
# ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ #
24+
@propagate_inbounds Base.setindex!{T,N}(block_arr::AbstractBlockArray{T,N}, v, block::Block{N}) = setblock!(block_arr, v, block.n...)
25+

‎src/blockarray.jl

Lines changed: 62 additions & 8 deletions
Original file line numberDiff line numberDiff line change
@@ -1,3 +1,6 @@
1+
# Note: Functions surrounded by a comment blocks are there because `Vararg` is sitll allocating.
2+
# When Vararg is fast enough, they can simply be removed
3+
14
##############
25
# BlockArray #
36
##############
@@ -33,37 +36,88 @@ end
3336
# Indexing #
3437
############
3538

39+
# vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv #
40+
function Base.getindex{T, N}(block_arr::BlockArray{T, N}, i::Int)
41+
@boundscheck checkbounds(block_arr, i)
42+
@inbounds v = block_arr[global2blockindex(block_arr.block_sizes, (i,))]
43+
return v
44+
end
45+
46+
function Base.getindex{T, N}(block_arr::BlockArray{T, N}, i::Int, j::Int)
47+
@boundscheck checkbounds(block_arr, i, j)
48+
@inbounds v = block_arr[global2blockindex(block_arr.block_sizes, (i,j))]
49+
return v
50+
end
51+
# ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ #
52+
3653
function Base.getindex{T, N}(block_arr::BlockArray{T, N}, i::Vararg{Int, N})
3754
@boundscheck checkbounds(block_arr, i...)
3855
@inbounds v = block_arr[global2blockindex(block_arr.block_sizes, i)]
3956
return v
4057
end
4158

59+
60+
# vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv #
61+
function Base.setindex!{T, N}(block_arr::BlockArray{T, N}, v, i::Int)
62+
@boundscheck checkbounds(block_arr, i)
63+
@inbounds block_arr[global2blockindex(block_arr.block_sizes, (i,))] = v
64+
return block_arr
65+
end
66+
67+
68+
function Base.setindex!{T, N}(block_arr::BlockArray{T, N}, v, i::Int, j::Int)
69+
@boundscheck checkbounds(block_arr, i, j)
70+
@inbounds block_arr[global2blockindex(block_arr.block_sizes, (i,j))] = v
71+
return block_arr
72+
end
73+
# ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ #
74+
75+
4276
function Base.setindex!{T, N}(block_arr::BlockArray{T, N}, v, i::Vararg{Int, N})
4377
@boundscheck checkbounds(block_arr, i...)
4478
@inbounds block_arr[global2blockindex(block_arr.block_sizes, i)] = v
4579
return block_arr
4680
end
4781

82+
# vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv #
83+
@propagate_inbounds getblock{T}(block_arr::BlockArray{T,1}, block_i::Int) = block_arr.blocks[block_i]
84+
@propagate_inbounds getblock{T}(block_arr::BlockArray{T,2}, block_i::Int, block_j::Int) = block_arr.blocks[block_i, block_j]
85+
# ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ #
4886
@propagate_inbounds getblock{T,N}(block_arr::BlockArray{T,N}, block::Vararg{Int, N}) = block_arr.blocks[block...]
4987

50-
function setblock!{T, N, R}(block_arr::BlockArray{T, N, R}, v, block::Vararg{Int, N})
51-
@boundscheck begin # TODO: Check if this eliminates the boundscheck with @inbounds
52-
for i in 1:N
53-
if size(v, i) != block_arr.block_sizes[i, block[i]]
54-
throw(ArgumentError(string("attempt to assign a $(size(v)) array to a ", ntuple(i -> block_arr.block_sizes[i][block[i]], Val{N}), " block")))
55-
end
88+
function _check_setblock!{T,N}(block_arr::BlockArray{T, N}, v, block::NTuple{N, Int})
89+
for i in 1:N
90+
if size(v, i) != block_arr.block_sizes[i, block[i]]
91+
throw(ArgumentError(string("attempt to assign a $(size(v)) array to a ", ntuple(i -> block_arr.block_sizes[i][block[i]], Val{N}), " block")))
5692
end
5793
end
94+
end
95+
96+
# vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv #
97+
function setblock!{T}(block_arr::BlockArray{T, 1}, v, block_i::Int)
98+
@boundscheck _check_setblock!(block_arr, v, (block_i,))
99+
@inbounds block_arr.blocks[block_i] = v
100+
return block_arr
101+
end
102+
103+
function setblock!{T}(block_arr::BlockArray{T, 2}, v, block_i::Int, block_j::Int)
104+
@boundscheck _check_setblock!(block_arr, v, (block_i, block_j))
105+
@inbounds block_arr.blocks[block_i, block_j] = v
106+
return block_arr
107+
end
108+
# ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ #
109+
110+
function setblock!{T, N}(block_arr::BlockArray{T, N}, v, block::Vararg{Int, N})
111+
@boundscheck _check_setblock!(block_arr, v, block)
58112
@inbounds block_arr.blocks[block...] = v
59113
return block_arr
60114
end
61115

62-
@propagate_inbounds function Base.getindex{T,N}(block_array::BlockArray{T, N}, block_index::BlockIndex)
116+
function Base.getindex{T,N}(block_array::BlockArray{T, N}, block_index::BlockIndex{N})
63117
return getblock(block_array, block_index.I...)[block_index.α...]
64118
end
65119

66-
@propagate_inbounds function Base.setindex!{T,N}(block_array::BlockArray{T, N}, v, block_index::BlockIndex)
120+
function Base.setindex!{T,N}(block_array::BlockArray{T, N}, v, block_index::BlockIndex{N})
67121
getblock(block_array, block_index.I...)[block_index.α...] = v
68122
end
69123

‎src/pseudo_blockarray.jl

Lines changed: 126 additions & 20 deletions
Original file line numberDiff line numberDiff line change
@@ -1,3 +1,6 @@
1+
# Note: Functions surrounded by a comment blocks are there because `Vararg` is sitll allocating.
2+
# When Vararg is fast enough, they can simply be removed
3+
14
####################
25
# PseudoBlockArray #
36
####################
@@ -39,36 +42,107 @@ end
3942
# Indexing #
4043
############
4144

45+
# vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv #
46+
function Base.getindex{T, N}(block_arr::PseudoBlockArray{T, N}, i::Int)
47+
@boundscheck checkbounds(block_arr, )
48+
@inbounds v = block_arr.blocks[i]
49+
return v
50+
end
51+
52+
function Base.getindex{T, N}(block_arr::PseudoBlockArray{T, N}, i::Int, j::Int)
53+
@boundscheck checkbounds(block_arr, i, j)
54+
@inbounds v = block_arr.blocks[i, j]
55+
return v
56+
end
57+
# ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ #
58+
4259
function Base.getindex{T, N}(block_arr::PseudoBlockArray{T, N}, i::Vararg{Int, N})
4360
@boundscheck checkbounds(block_arr, i...)
4461
@inbounds v = block_arr.blocks[i...]
4562
return v
4663
end
4764

65+
# vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv #
66+
function Base.setindex!{T, N}(block_arr::PseudoBlockArray{T, N}, v, i::Int)
67+
@boundscheck checkbounds(block_arr, i)
68+
@inbounds block_arr.blocks[i] = v
69+
return block_arr
70+
end
71+
72+
function Base.setindex!{T, N}(block_arr::PseudoBlockArray{T, N}, v, i::Int, j::Int)
73+
@boundscheck checkbounds(block_arr, i, j)
74+
@inbounds block_arr.blocks[i, j] = v
75+
return block_arr
76+
end
77+
# ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ #
78+
4879
function Base.setindex!{T, N}(block_arr::PseudoBlockArray{T, N}, v, i::Vararg{Int, N})
4980
@boundscheck checkbounds(block_arr, i...)
5081
@inbounds block_arr.blocks[i...] = v
5182
return block_arr
5283
end
5384

85+
# vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv #
86+
function getblock{T,N}(block_arr::PseudoBlockArray{T,N}, block_i::Int)
87+
range = globalrange(block_arr.block_sizes, (block_i,))
88+
return block_arr.blocks[range[1]]
89+
end
90+
91+
function getblock{T,N}(block_arr::PseudoBlockArray{T,N}, block_i::Int, block_j::Int)
92+
range = globalrange(block_arr.block_sizes, (block_i, block_j))
93+
return block_arr.blocks[range[1], range[2]]
94+
end
95+
# ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ #
96+
5497
function getblock{T,N}(block_arr::PseudoBlockArray{T,N}, block::Vararg{Int, N})
5598
range = globalrange(block_arr.block_sizes, block)
5699
return block_arr.blocks[range...]
57100
end
58101

59-
# TODO: Can this be written efficiently without a generated function?
60-
@generated function getblock!{T,N}(x, block_arr::PseudoBlockArray{T,N}, block::Vararg{Int, N})
61-
return quote
102+
function _check_getblock!{T, N}(blockrange, x, block_arr::PseudoBlockArray{T,N}, block::NTuple{N, Int})
103+
for i in 1:N
104+
if size(x, i) != length(blockrange[i])
105+
throw(ArgumentError(string("attempt to assign a ", ntuple(i -> block_arr.block_sizes[i, block[i]], Val{N}), " block to a $(size(x)) array")))
106+
end
107+
end
108+
end
62109

63-
blockrange = globalrange(block_arr.block_sizes, block)
110+
# vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv #
111+
@propagate_inbounds function getblock!{T}(x, block_arr::PseudoBlockArray{T,1}, block_i::Int)
112+
blockrange = globalrange(block_arr.block_sizes, (block_i,))
113+
@boundscheck _check_getblock!(blockrange, x, block_arr, (block_i,))
64114

65-
@boundscheck begin
66-
for i in 1:N
67-
if size(x, i) != length(blockrange[i])
68-
throw(ArgumentError(string("attempt to assign a ", ntuple(i -> block_arr.block_sizes[i, block[i]], Val{N}), " block to a $(size(x)) array")))
69-
end
70-
end
115+
arr = block_arr.blocks
116+
k_1 = 1
117+
@inbounds for i in blockrange[1]
118+
x[k_1] = arr[i]
119+
k_1 += 1
120+
end
121+
return x
122+
end
123+
124+
@propagate_inbounds function getblock!{T}(x, block_arr::PseudoBlockArray{T,2}, block_i::Int, block_j::Int)
125+
blockrange = globalrange(block_arr.block_sizes, (block_i,block_j))
126+
@boundscheck _check_getblock!(blockrange, x, block_arr, (block_i, block_j))
127+
128+
arr = block_arr.blocks
129+
k_2 = 1
130+
@inbounds for j in blockrange[2]
131+
k_1 = 1
132+
for i in blockrange[1]
133+
x[k_1, k_2] = arr[i, j]
134+
k_1 += 1
71135
end
136+
k_2 += 1
137+
end
138+
return x
139+
end
140+
# ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ #
141+
142+
@generated function getblock!{T,N}(x, block_arr::PseudoBlockArray{T,N}, block::Vararg{Int, N})
143+
return quote
144+
blockrange = globalrange(block_arr.block_sizes, block)
145+
@boundscheck _check_getblock!(blockrange, x, block_arr, block)
72146

73147
arr = block_arr.blocks
74148
@nexprs $N d -> k_d = 1
@@ -81,23 +155,55 @@ end
81155
end
82156
end
83157

84-
# TODO: Can this be written efficiently without a generated function?
85-
@generated function setblock!{T, N, R}(block_arr::PseudoBlockArray{T, N, R}, v, block::Vararg{Int, N})
86-
return quote
87-
@boundscheck begin # TODO: Check if this eliminates the boundscheck with @inbounds
88-
for i in 1:N
89-
if size(v, i) != block_arr.block_sizes[i, block[i]]
90-
throw(ArgumentError(string("attempt to assign a $(size(v)) array to a ", ntuple(i -> block_arr.block_sizes[i][block[i]], Val{N}), " block")))
91-
end
92-
end
158+
function _check_setblock!{T, N}(blockrange, x, block_arr::PseudoBlockArray{T,N}, block::NTuple{N, Int})
159+
for i in 1:N
160+
if size(x, i) != block_arr.block_sizes[i, block[i]]
161+
throw(ArgumentError(string("attempt to assign a $(size(x)) array to a ", ntuple(i -> block_arr.block_sizes[i][block[i]], Val{N}), " block")))
93162
end
163+
end
164+
end
165+
166+
# vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv #
167+
@propagate_inbounds function setblock!{T}(block_arr::PseudoBlockArray{T,1}, x, block_i::Int)
168+
blockrange = globalrange(block_arr.block_sizes, (block_i,))
169+
@boundscheck _check_setblock!(blockrange, x, block_arr, (block_i,))
170+
171+
arr = block_arr.blocks
172+
k_1 = 1
173+
@inbounds for i in blockrange[1]
174+
arr[i] = x[k_1]
175+
k_1 += 1
176+
end
177+
return x
178+
end
94179

180+
@propagate_inbounds function setblock!{T}(block_arr::PseudoBlockArray{T,2}, x, block_i::Int, block_j::Int)
181+
blockrange = globalrange(block_arr.block_sizes, (block_i,block_j))
182+
@boundscheck _check_setblock!(blockrange, x, block_arr, (block_i, block_j))
183+
184+
arr = block_arr.blocks
185+
k_2 = 1
186+
@inbounds for j in blockrange[2]
187+
k_1 = 1
188+
for i in blockrange[1]
189+
arr[i, j] = x[k_1, k_2]
190+
k_1 += 1
191+
end
192+
k_2 += 1
193+
end
194+
return x
195+
end
196+
# ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ #
197+
198+
@generated function setblock!{T, N}(block_arr::PseudoBlockArray{T, N}, x, block::Vararg{Int, N})
199+
return quote
95200
blockrange = globalrange(block_arr.block_sizes, block)
201+
@boundscheck _check_setblock!(blockrange, x, block_arr, block)
96202
arr = block_arr.blocks
97203
@nexprs $N d -> k_d = 1
98204
@inbounds begin
99205
@nloops $N i (d->(blockrange[d])) (d-> k_{d-1}=1) (d-> k_d+=1) begin
100-
(@nref $N arr i) = (@nref $N v k)
206+
(@nref $N arr i) = (@nref $N x k)
101207
end
102208
end
103209
end

0 commit comments

Comments
 (0)