Skip to content

Commit 945ef97

Browse files
author
Ryan Bennink
committed
Preliminary support for matrix multiplication.
- New methods for * for combinations of OffsetArray and AbstractArray Unfortunately this introduces method ambiguities with Base. - As an alternative, tried defining a new promote_rule, but somehow it doesn't get invoked. New convenience constructors - Wrap any (non-offset) array as an OffsetArray with standard indices - Remove layer of indirection when constructing with explicit offsets Cosmetic edits to constructor section for improved readability
1 parent e22f56e commit 945ef97

File tree

1 file changed

+36
-2
lines changed

1 file changed

+36
-2
lines changed

src/OffsetArrays.jl

Lines changed: 36 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -3,32 +3,43 @@ VERSION < v"0.7.0-beta2.199" && __precompile__()
33
module OffsetArrays
44

55
using Base: Indices, tail, @propagate_inbounds
6+
import Base: (*), convert, promote_rule
7+
68
@static if !isdefined(Base, :IdentityUnitRange)
79
const IdentityUnitRange = Base.Slice
810
else
911
using Base: IdentityUnitRange
1012
end
1113

12-
export OffsetArray, OffsetVector
14+
export OffsetArray, OffsetVector, OffsetMatrix
1315

1416
struct OffsetArray{T,N,AA<:AbstractArray} <: AbstractArray{T,N}
1517
parent::AA
1618
offsets::NTuple{N,Int}
1719
end
1820
OffsetVector{T,AA<:AbstractArray} = OffsetArray{T,1,AA}
21+
OffsetMatrix{T,AA<:AbstractArray} = OffsetArray{T,2,AA}
1922

23+
# Construct from existing array
2024
OffsetArray(A::AbstractArray{T,N}, offsets::NTuple{N,Int}) where {T,N} =
2125
OffsetArray{T,N,typeof(A)}(A, offsets)
2226
OffsetArray(A::AbstractArray{T,N}, offsets::Vararg{Int,N}) where {T,N} =
2327
OffsetArray(A, offsets)
28+
OffsetArray(A::AbstractArray{T,0}) where {T} = OffsetArray{T,0,typeof(A)}(A, ())
29+
OffsetArray(A::AbstractArray) = convert(OffsetArray, A)
30+
OffsetArray(A::OffsetArray) = A
31+
32+
# Convert from existing array
33+
convert(::Type{OffsetArray}, A::AbstractArray{T,N}) where {T,N} = OffsetArray(A, ntuple(x->0, N))
2434

35+
# Create an uninitialized OffsetArray with given element type
2536
const ArrayInitializer = Union{UndefInitializer, Missing, Nothing}
2637
OffsetArray{T,N}(init::ArrayInitializer, inds::Indices{N}) where {T,N} =
2738
OffsetArray{T,N,Array{T,N}}(Array{T,N}(init, map(indexlength, inds)), map(indexoffset, inds))
2839
OffsetArray{T}(init::ArrayInitializer, inds::Indices{N}) where {T,N} = OffsetArray{T,N}(init, inds)
40+
# Same thing, but taking multiple args for offsets/indices
2941
OffsetArray{T,N}(init::ArrayInitializer, inds::Vararg{AbstractUnitRange,N}) where {T,N} = OffsetArray{T,N}(init, inds)
3042
OffsetArray{T}(init::ArrayInitializer, inds::Vararg{AbstractUnitRange,N}) where {T,N} = OffsetArray{T,N}(init, inds)
31-
OffsetArray(A::AbstractArray{T,0}) where {T} = OffsetArray{T,0,typeof(A)}(A, ())
3243

3344
# OffsetVector constructors
3445
OffsetVector(A::AbstractVector, offset) = OffsetArray(A, offset)
@@ -59,6 +70,7 @@ OffsetArray(A::AbstractArray{T,N}, inds::Vararg{AbstractUnitRange,N}) where {T,N
5970
OffsetArray(A, inds)
6071

6172
# avoid a level of indirection when nesting OffsetArrays
73+
OffsetArray(A::OffsetArray, offsets::NTuple{N,Int}) where {N} = OffsetArray(parent(A), offsets)
6274
function OffsetArray(A::OffsetArray, inds::NTuple{N,AbstractUnitRange}) where {N}
6375
OffsetArray(parent(A), inds)
6476
end
@@ -283,6 +295,28 @@ end
283295
no_offset_view(A::OffsetArray) = no_offset_view(parent(A))
284296

285297

298+
# Quick hack for matrix multiplication.
299+
# Ideally, one would instead improve LinearAlgebra's support of custom indexing.
300+
function (*)(A::OffsetMatrix, B::OffsetMatrix)
301+
matmult_check_axes(A, B)
302+
C = OffsetArray(parent(A) * parent(B), (axes(A,1), axes(B,2)))
303+
end
304+
305+
function (*)(A::OffsetMatrix, B::OffsetVector)
306+
matmult_check_axes(A, B)
307+
C = OffsetArray(parent(A) * parent(B), axes(A,1))
308+
end
309+
matmult_check_axes(A, B) = axes(A, 2) == axes(B, 1) || error("axes(A,2) must equal axes(B,1)")
310+
311+
(*)(A::OffsetMatrix, B::AbstractMatrix) = A * OffsetArray(B)
312+
(*)(A::OffsetMatrix, B::AbstractVector) = A * OffsetArray(B)
313+
(*)(A::AbstractMatrix, B::OffsetArray) = OffsetArray(A) * B
314+
(*)(A::AbstractVector, B::OffsetArray) = OffsetArray(A) * B
315+
316+
# An alternative to the above four methods would be to use promote_rule, but it doesn't get invoked
317+
# promote_rule(::Type{A1}, ::Type{A2}) where A1<:AbstractArray{<:Any,N} where A2<:OffsetArray{<:Any,N,A3} where {N,A3} = OffsetArray{eltype(promote_type(A1, A3)), N, promote_type(A1, A3)}
318+
319+
286320
####
287321
# work around for segfault in searchsorted*
288322
# https://github.com/JuliaLang/julia/issues/33977

0 commit comments

Comments
 (0)