Resize! and multidimensional arrays

General Usage
1

Background

We want to make use of adaptive mesh refinement (AMR) for our hyperbolic PDE solver package Trixi.jl. For us, the most natural data layout is given by multidimensional arrays. To use AMR, we need to resize! the last dimension of these arrays. Currently, multidimensional Arrays can not be resize!ed and such a feature will not be implemented in the foreseeable future because of performance considerations.

Alternatives

One option might be to use ElasticArrays.jl. However, needing to use ElasticArray doesn’t feel completely good to me, since we also want to experiment with other array types such as PaddedMatrices.jl, see https://github.com/trixi-framework/Trixi.jl/issues/166. Then, we would need to wrap an Array inside something from PaddedMatrices.jl inside something from ElasticArrays.jl - or the other way round? Is that possible at all? If we go further, I fear this looks like it could easily explode.

We could try to use resize!able Vectors which we somehow wrap inside Trixi.jl as our multidimensional Arrays. The most straightforward approach I can think of fails:

  julia> v = ones(6)
  6-element Array{Float64,1}:
  1.0
  1.0
  1.0
  1.0
  1.0
  1.0

  julia> function foo(v)
            vm = reshape(v, (2, length(v)Ă·2))
            vm[1, 2]
        end
  foo (generic function with 1 method)

  julia> foo(v)
  1.0

  julia> v[3] = 5
  5

  julia> foo(v)
  5.0

  julia> resize!(v, 8)
  ERROR: cannot resize array with shared data
  Stacktrace:
  [1] _growend! at ./array.jl:892 [inlined]
  [2] resize!(::Array{Float64,1}, ::Int64) at ./array.jl:1085
  [3] top-level scope at REPL[6]:1

It’s okay that resize! doesn’t work when the data is shared.

Question 1: Is there a way to be able to use resize! on a Vector when the shared data went out of scope?

Instead, it looks like we need to use something like

julia> function foo(v)
          vm = unsafe_wrap(Array{eltype(v), 2}, pointer(v), (2, length(v)Ă·2))
          vm[1, 2]
      end

That’s okay and works for us in the current approach. It will also allow us to wrap the Vector with more type information on the sizes, e.g. using HybridArrays.jl. However, we need to take care, e.g. if we need to copy a given Vector locally, wrap it as a multidimensional Array - we must use GC.@preserve.

Question 2: Is there a possibility to use something like that without having to GC.@preserve copies of Vectors which are only used indirectly via multidimensional wrapped versions afterwards?

Question 3: Is there a better alternative?

1 Like
2

Dumb question: why does your algorithm have to be in-place?

1 Like
3

That’s not a dump question at all. In the past, we’ve just created new Arrays when doing AMR and that was totally okay for us. However, we want to make use of OrdinaryDiffEq.jl etc. in the future. There, it’s only possible to resize! intermediate caches used for the time integration but not to set them to something completely new.

Of course, one could also argue that this is a restriction of OrdinaryDiffEq.jl and should be fixed there. On the other hand, we can reduce allocations and save runtime if we are able to re-use already allocated memory. For example, standard Vectors won’t allocate on every call to resize!, which makes this approach very efficient for us.

4

I still do not understand why. But in any case,

  1. If you need to resize a lot, just allocate a “large enough” array and keep track of the dimensions you use at the momene. You can wrap this in a struct and implement the AbstractArray interface very simply.

  2. If you need to resize very infrequently, just copying over the data to a freshly allocated container should be reasonably cheap.

1 Like
5

Thanks for investing your time to answer my questions. We will continue to discuss options internally and use some kind of wrappers of Vectors (as also done by ElasticArrays.jl).

1 Like