Broadcast matrix valued function, output into existing array

  1. This is the wrong order for locality
  2. Working with lots of 2x2 arrays is an ideal case for StaticArrays — instead of a 3d array, use a 1d array of 2x2 SMatrix’s. This will also fix (1).

Using an array of SMatrix also makes it easier to use broadcasting. For example:

using StaticArrays

function matop(a, b)
    return a * @SMatrix[b a; a b];
end

N = 100;
aar = range(0,0.1,N);
bar = range(0,0.05,N);
car = matop.(aar, bar)

and then you can access elements as e.g. car[i][1,2]. This is much more efficient than your previous implementation (especially if you put it into a function to avoid globals), because it no longer heap-allocates an Array on every call to matop.

4 Likes