- This is the wrong order for locality
- 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.