Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

Preliminary support for matrix multiplication. #93

Open
wants to merge 1 commit into
base: master
Choose a base branch
from
Open

Preliminary support for matrix multiplication. #93

Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
38 changes: 36 additions & 2 deletions src/OffsetArrays.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -3,32 +3,43 @@ VERSION < v"0.7.0-beta2.199" && __precompile__()
module OffsetArrays

using Base: Indices, tail, @propagate_inbounds
import Base: (*), convert, promote_rule

@static if !isdefined(Base, :IdentityUnitRange)
const IdentityUnitRange = Base.Slice
else
using Base: IdentityUnitRange
end

export OffsetArray, OffsetVector
export OffsetArray, OffsetVector, OffsetMatrix

struct OffsetArray{T,N,AA<:AbstractArray} <: AbstractArray{T,N}
parent::AA
offsets::NTuple{N,Int}
end
OffsetVector{T,AA<:AbstractArray} = OffsetArray{T,1,AA}
OffsetMatrix{T,AA<:AbstractArray} = OffsetArray{T,2,AA}

# Construct from existing array
OffsetArray(A::AbstractArray{T,N}, offsets::NTuple{N,Int}) where {T,N} =
OffsetArray{T,N,typeof(A)}(A, offsets)
OffsetArray(A::AbstractArray{T,N}, offsets::Vararg{Int,N}) where {T,N} =
OffsetArray(A, offsets)
OffsetArray(A::AbstractArray{T,0}) where {T} = OffsetArray{T,0,typeof(A)}(A, ())
OffsetArray(A::AbstractArray) = convert(OffsetArray, A)
OffsetArray(A::OffsetArray) = A

# Convert from existing array
convert(::Type{OffsetArray}, A::AbstractArray{T,N}) where {T,N} = OffsetArray(A, ntuple(x->0, N))

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

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

# avoid a level of indirection when nesting OffsetArrays
OffsetArray(A::OffsetArray, offsets::NTuple{N,Int}) where {N} = OffsetArray(parent(A), offsets)
function OffsetArray(A::OffsetArray, inds::NTuple{N,AbstractUnitRange}) where {N}
OffsetArray(parent(A), inds)
end
Expand Down Expand Up @@ -283,6 +295,28 @@ end
no_offset_view(A::OffsetArray) = no_offset_view(parent(A))


# Quick hack for matrix multiplication.
# Ideally, one would instead improve LinearAlgebra's support of custom indexing.
function (*)(A::OffsetMatrix, B::OffsetMatrix)
matmult_check_axes(A, B)
C = OffsetArray(parent(A) * parent(B), (axes(A,1), axes(B,2)))
end

function (*)(A::OffsetMatrix, B::OffsetVector)
matmult_check_axes(A, B)
C = OffsetArray(parent(A) * parent(B), axes(A,1))
end
matmult_check_axes(A, B) = axes(A, 2) == axes(B, 1) || error("axes(A,2) must equal axes(B,1)")

(*)(A::OffsetMatrix, B::AbstractMatrix) = A * OffsetArray(B)
(*)(A::OffsetMatrix, B::AbstractVector) = A * OffsetArray(B)
(*)(A::AbstractMatrix, B::OffsetArray) = OffsetArray(A) * B
(*)(A::AbstractVector, B::OffsetArray) = OffsetArray(A) * B

# An alternative to the above four methods would be to use promote_rule, but it doesn't get invoked
# 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)}


####
# work around for segfault in searchsorted*
# https://github.com/JuliaLang/julia/issues/33977
Expand Down