Convert (N)D Array to (N-2)D Array of SMatrix elements

General Usage
, ,
1

I want to do matrix inversion on a stack of matrices, based on this post. I’m having difficulty recasting the arrays though. In my case, the array will come as data or an input to a function, but here is a minimal example of what I’m trying to do ( fill something like the GMat array with the data from A)

using StaticArrays
A = rand(3,3,5,10)
GMat = Array{SMatrix{3,3,Float64,9}}(undef, 5,10)
reinterpret(Array{SMatrix{3,3,Float64}},A)

How would I do this conversion? Sorry for a very simple question.

2

Here’s how I’d write it:

using StaticArrays

let A = rand(3,3,5,10)
    A2 = reshape(A, 9, 5, 10)
    reinterpret(reshape, SMatrix{3,3, Float64, 9}, A2)
end

Which gives

5×10 reinterpret(reshape, SMatrix{3, 3, Float64, 9}, ::Array{Float64, 3}) with eltype SMatrix{3, 3, Float64, 9}:
 [0.643427 0.688092 0.518653; 0.1005 0.787188 0.146908; 0.859216 0.355823 0.12054]       …  [0.747454 0.0155871 0.349355; 0.589688 0.681063 0.839895; 0.302963 0.0427983 0.0466762]
 [0.842058 0.147054 0.143648; 0.666544 0.909276 0.890523; 0.956899 0.66385 0.901945]        [0.263943 0.274761 0.654176; 0.281182 0.504185 0.68238; 0.273182 0.763132 0.467453]
 [0.758898 0.763019 0.874907; 0.0393675 0.934689 0.670948; 0.060328 0.0557292 0.257968]     [0.205689 0.0310269 0.542356; 0.505019 0.304702 0.582388; 0.547555 0.144255 0.62276]
 [0.352791 0.846856 0.362235; 0.107646 0.312424 0.801409; 0.728911 0.870063 0.53615]        [0.658811 0.642384 0.467456; 0.633928 0.225012 0.698295; 0.155892 0.937625 0.818803]
 [0.49822 0.00705119 0.537992; 0.662411 0.338279 0.176099; 0.950147 0.660809 0.482299]      [0.0713601 0.571688 0.406371; 0.531959 0.618212 0.487209; 0.0732696 0.500866 0.481095]
4 Likes
3

For those of us who are not as smart as Mason :sweat:, the TensorCast package is perhaps more within reach:

using TensorCast
A = rand(3,3,5,10)
@cast B[k,l]{i,j} := A[i,j,k,l]

NB:
the curly brackets are used for the static array slices of desired output: 5x10 Matrix of SMatrix{3,3}

4 Likes
4 
julia> gsm=(SMatrix{3,3}(sm) for sm in eachslice(A,dims=(3,4)));

julia> collect(gsm)
5×10 Matrix{SMatrix{3, 3, Float64, 9}}:

or

SMatrix{size(A)[1:2]...}.(eachslice(A,dims=(3,4)))
1 Like