Skip to content

Commit f34598c

Browse files
Ryan Benninkjishnub
Ryan Bennink
authored andcommitted
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 93e60d3 commit f34598c

File tree

1 file changed

+26
-0
lines changed

1 file changed

+26
-0
lines changed

src/OffsetArrays.jl

Lines changed: 26 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -1,6 +1,8 @@
11
module OffsetArrays
22

33
using Base: Indices, tail, @propagate_inbounds
4+
import Base: (*), convert, promote_rule
5+
46
@static if !isdefined(Base, :IdentityUnitRange)
57
const IdentityUnitRange = Base.Slice
68
else
@@ -45,10 +47,12 @@ OffsetArray(A::AbstractArray{T,N}, offsets::Vararg{Int,N}) where {T,N} =
4547
OffsetArray(A, offsets)
4648
OffsetArray(A::AbstractArray{T,0}) where {T} = OffsetArray(A, ())
4749

50+
# Create an uninitialized OffsetArray with given element type
4851
const ArrayInitializer = Union{UndefInitializer, Missing, Nothing}
4952
OffsetArray{T,N}(init::ArrayInitializer, inds::Indices{N}) where {T,N} =
5053
OffsetArray(Array{T,N}(init, map(indexlength, inds)), map(indexoffset, inds))
5154
OffsetArray{T}(init::ArrayInitializer, inds::Indices{N}) where {T,N} = OffsetArray{T,N}(init, inds)
55+
# Same thing, but taking multiple args for offsets/indices
5256
OffsetArray{T,N}(init::ArrayInitializer, inds::Vararg{AbstractUnitRange,N}) where {T,N} = OffsetArray{T,N}(init, inds)
5357
OffsetArray{T}(init::ArrayInitializer, inds::Vararg{AbstractUnitRange,N}) where {T,N} = OffsetArray{T,N}(init, inds)
5458

@@ -403,6 +407,28 @@ end
403407
no_offset_view(A::OffsetArray) = no_offset_view(parent(A))
404408

405409

410+
# Quick hack for matrix multiplication.
411+
# Ideally, one would instead improve LinearAlgebra's support of custom indexing.
412+
function (*)(A::OffsetMatrix, B::OffsetMatrix)
413+
matmult_check_axes(A, B)
414+
C = OffsetArray(parent(A) * parent(B), (axes(A,1), axes(B,2)))
415+
end
416+
417+
function (*)(A::OffsetMatrix, B::OffsetVector)
418+
matmult_check_axes(A, B)
419+
C = OffsetArray(parent(A) * parent(B), axes(A,1))
420+
end
421+
matmult_check_axes(A, B) = axes(A, 2) == axes(B, 1) || error("axes(A,2) must equal axes(B,1)")
422+
423+
(*)(A::OffsetMatrix, B::AbstractMatrix) = A * OffsetArray(B)
424+
(*)(A::OffsetMatrix, B::AbstractVector) = A * OffsetArray(B)
425+
(*)(A::AbstractMatrix, B::OffsetArray) = OffsetArray(A) * B
426+
(*)(A::AbstractVector, B::OffsetArray) = OffsetArray(A) * B
427+
428+
# An alternative to the above four methods would be to use promote_rule, but it doesn't get invoked
429+
# 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)}
430+
431+
406432
####
407433
# work around for segfault in searchsorted*
408434
# https://github.com/JuliaLang/julia/issues/33977

0 commit comments

Comments
 (0)