OffsetArrays vs Arrays speed

New to Julia
1

Hello Julia wizards!

I’ve been exploring this language for about a week, and so far I’m very very intrigued. I have a simulation code I’m writing that would greatly benefit from using OffsetArrays, but everything I try results in a massive slowdown.

I’ve boiled down the pattern to a simple example, and all of the documentation I can understand suggests that the Array version and the OffsetArray version should be the same speed. The benchmark result is:

Arrays: 730.958 μs (0 allocations: 0 bytes)
OffsetArrays: 129.026 ms (9998470 allocations: 167.82 MiB)

My guess is all of those unexpected allocations are the problem, but I have no idea where they are coming from. Help, please?

using BenchmarkTools
using OffsetArrays

mutable struct MyArray
    a::Array{Float64, 1}
    b::Array{Float64, 1}

    MyArray() = new()
    MyArray(a, b) = new(a,b)
end

mutable struct MyOffsetArray
    a::OffsetArray{Float64, 1}
    b::OffsetArray{Float64, 1}

    MyOffsetArray() = new()
    MyOffsetArray(a, b) = new(a,b)
end

function myoperate(arg::Union{MyArray, MyOffsetArray})
    for i = 1:length(arg.b)
        arg.b[i] = arg.b[i] + 0.5 * (arg.a[i] + arg.a[i+1])
    end
end

function timeit()
    nx = 1000000

    a = rand(nx+1)
    b = rand(nx)

    array = MyArray(a, b)
    offset = MyOffsetArray(a, b)

    @btime myoperate($array)
    @btime myoperate($offset)
end

timeit()
2

The problem is that OffsetArray{Float64, 1} is an abstract type (since the offset is also part of the type). You probably want to use

mutable struct MyParametricArray{T<:AbstractArray}
    a::T
    b::T
    MyParametricArray() = new()
    MyParametricArray(a, b) = new(a,b)
end

which will specialize correctly.

7 Likes
3

Wow, thanks! That led me the right way (and it would have been a long time before I figured that out…)

Arrays: 734.166 μs (0 allocations: 0 bytes)
OffsetArrays: 128.470 ms (9998470 allocations: 167.82 MiB)
OffsetArrays (parametric): 1.276 ms (0 allocations: 0 bytes)

Final code:

using BenchmarkTools
using OffsetArrays

mutable struct MyArray
    a::Array{Float64, 1}
    b::Array{Float64, 1}
end

mutable struct MyOffsetArray
    a::OffsetArray{Float64, 1}
    b::OffsetArray{Float64, 1}
end

mutable struct MyParametricArray{T<:AbstractArray}
    a::T
    b::T
end

function myoperate(arg::Union{MyArray, MyOffsetArray, MyParametricArray})
    for i = 1:length(arg.b)
        arg.b[i] = arg.b[i] + 0.5 * (arg.a[i] + arg.a[i+1])
    end
end

function timeit()
    nx = 1000000

    a = rand(nx+1)
    b = rand(nx)

    ao = OffsetArrays.Origin(1)(a)
    bo = OffsetArrays.Origin(1)(b)

    array = MyArray(a, b)
    offset = MyOffsetArray(ao, bo)
    parameter = MyParametricArray(ao, bo)

    @btime myoperate($array)
    @btime myoperate($offset)
    @btime myoperate($parameter)
end

timeit()
4

This seems highly risky. If a and b are OffsetArrays, those will likely not be valid indices (that’s the point of OffsetArrays, really). Use eachindex(a, b) instead.

Also, consider whether you need your type to be mutable. Immutable types are generally more efficient, and you can still mutate mutable fields of an immutable type, so your array code would still work.

6 Likes
5

Yes, I use eachindex in the real code. I was just trying to boil the problem down to the minimal case here.

Thanks for the tip about mutable.