Similar MArray with different size

General Usage
, ,
1

I want to create an array similar to some other array, which could be an Array, a MArray, or some other kind of mutable array-like container. But I would also like to set its dimensions to m×n. I’ve tried similar and it doesn’t seem to work. Consider the following:

using StaticArrays
m = n = 2
println(typeof(similar([1, 2, 3])))
println(typeof(similar([1, 2, 3], m, n)))
println(typeof(similar(@MVector([1, 2, 3]))))
println(typeof(similar(@MVector([1, 2, 3]), m, n)))

which outputs

Vector{Int64}
Matrix{Int64}
MVector{3, Int64}
Matrix{Int64}

I would have hoped

typeof(similar(@MVector([1, 2, 3]), m, n))

to output

MMatrix{2, 2, Int64, 4}

Is this expectation unreasonable? To get the behavior I want, the only way I have found for now is to test explicitly:

x = @MVector([1, 2, 3])
typeof(isa(x, MArray) ? MMatrix{m, n, eltype(x)}(undef) : similar(x, m, n))

which outputs what I expect

MMatrix{2, 2, Int64, 4}

Of course, that just works for MArrays. Is there a better way?

2

I think reshape is what you want, see example here.

3

This is slightly unreasonable because in your similar variables n and m have values that are not statically known, so similar that returns MMatrix would be type unstable. StaticArrays.jl generally prefers to fall back to normal Arrays when the alternative is type instability, and IMO that choice makes sense. Type-unstable StaticArrays.jl code is generally slower than just using Array.

It depends on the context. Is this just for top level data allocation before actual performance-sensitive part is run? If so, then it can’t get much nicer than your explicit test. If you need to call similar in actually performance-critical part, then you likely need some additional adjustments to your code to ensure type stability.

4

You could write a similar method for MArray that takes Val arguments (for sizes that are known statically) pretty easily.

5

Yes, and there is already similar that takes Size which does pretty much the same thing.

2 Likes
6

Thanks everyone for the very good suggestions and warnings. I think I came up with a decent solution. As @stevengj pointed out, it is pretty easy to create a new similar function. As @mateuszbaran pointed out, we must be mindful of type stability. To test the speed gain, I just execute mul! 1000 times.

using StaticArrays, BenchmarkTools, LinearAlgebra

similar2(x, dims) = similar(x, dims)
similar2(x :: MArray, dims) = MArray{Tuple{dims...}, eltype(x), length(dims), prod(dims)}(undef)

function mul_many_times!(out,A,x, n)
    for _ in 1:n
        mul!(out,A,x)
    end
end

function test_mul(x)
    d = length(x)
    A = similar2(x,(d,d))
    out = similar2(x,d)
    out .= 1
    return mul_many_times!(out,A,x,1000)
end

x = [1,2,3,4]
xm = @MVector [1,2,3,4]

Then

@btime test_mul($x)
@btime test_mul($xm)

outputs

  12.200 μs (2 allocations: 288 bytes)
  3.375 μs (6 allocations: 384 bytes)
7

hmm, is this just used to generate a Matrix type from a Vector type? ArrayInterface docs mention a similar problem: Julia's Extended Linear Algebra Interface · ArrayInterface.jl

zeromat(u) = u .* u' .* false

function mul_many_times!(out,A,x, n)
    for _ in 1:n
        mul!(out,A,x)
    end
end

function test_mul(x)
    A = zeromat(x)
    out = similar(x)
    out .= 1
    return mul_many_times!(out,A,x,1000)
end

that gives:

julia> @btime test_mul($x)
  25.600 μs (2 allocations: 288 bytes)

julia> @btime Clapeyron.test_mul($xm)
  8.900 μs (2 allocations: 192 bytes)

so the additional allocations in the static case are gone.

1 Like
8

You’re welcome :slightly_smiling_face: .

Additional allocations are also gone if the MArray similar2 is marked as @inline.

2 Likes