Inner repeat without allocations

I want to achieve the effect of the following code:

X = randn(1000, 1000) 
Y = zeros(2000, 2000)
Y .= repeat(X, inner=(2,2))

The first two lines are just initialization. I would like the third line to run without allocations. What would be the best way to do this? Ideally, I’d be able to do this with broadcasting / dot notation / array functions from the standard library, rather than with for loops over 1000 x 1000 arrays. This is because my arrays can get quite large and performance is critical, and I’d like to take advantage of whatever multithreading / other optimizations those functions do rather than be responsible for writing my own efficient for loop.

I’m not sure if there’s an allocation-free version of repeat (although there should be), but this approach will always allocate the full result of calling repeat, then mutate Y to contain that result. So you really need something different than this “pure” version of repeat you’re using.

Hmm, I think it may be the case that an inner repeat is challenging to do without allocations, even with an output array? An intermediary output array may be necessary.

One can think of an inner repeat as an outer repeat followed by a permutedims, and I noticed that the permutedims! function doesn’t work in place – it requires a source and destination. This might be because we’d always like to access memory in column-major order?

So, to turn a 1000 x 1000 array X into a 2000 x 2000 array Y with inner repetition, we might need a second 2000 x 2000 array Y2. Then, we could do something like:

Y2 = reshape(Y2, 1000, 1000, 2, 2)
for i in 1:2
    for j in 1:2
        Y2[:, :, i, j] .= X
    end
end
Y = reshape(Y, 2, 1000, 2, 1000)
permutedims!(Y, Y2, (4, 1, 3, 2))
Y = reshape(Y, 2000, 2000)

What are people’s thoughts on my line of thinking and solution?

(Wait, this doesn’t seem right, Y2 could easily just be a view into X)

You don’t need permutedims, just reshape:

julia> X = rand(Int8, 3, 4);

julia> Y = zeros(6, 8);

julia> reshape(Y, 2, 3, 2, 4) .= reshape(X, 1, 3, 1, 4);

julia> Y == repeat(X, inner=(2,2))
true

However, if you have large arrays and care about performance, then finding a way not to repeat things might be even better.

1 Like