 # Replacing normal arrays with static arrays

I’m trying to perform a 3D rotation on a diagonal matrix, which seems to work with the code below.

``````using LinearAlgebra
using StaticArrays

function euler(φ, θ, ψ)

R = Array{Float64}(undef,(3,3))
cosφ = cos(φ); cosψ = cos(ψ); cosθ = cos(θ)
sinφ = sin(φ); sinψ = sin(ψ); sinθ = sin(θ)

R[1,1] = cosφ*cosθ*cosψ - sinφ*sinψ
R[1,2] = sinφ*cosθ*cosψ + cosφ*sinψ
R[1,3] = -sinθ*cosψ

R[2,1] = -cosφ*cosθ*sinψ - sinφ*cosψ
R[2,2] = -sinφ*cosθ*sinψ + cosφ*cosψ
R[2,3] = sinθ*sinψ

R[3,1] = cosφ*sinθ
R[3,2] = sinφ*sinθ
R[3,3] = cosθ

return R

end

R = euler(pi,0.0,0.0)
v = [1.0; 2.0 ; 3.0]
m = transpose(R) * Diagonal(v) * R
``````

It would seem a good case for static arrays, but I cannot seem to get the right syntax. Even basic things like

``````
function euler(φ, θ, ψ)

R = SMatrix{3,3,Float64,9}
cosφ = cos(φ); cosψ = cos(ψ); cosθ = cos(θ)
sinφ = sin(φ); sinψ = sin(ψ); sinθ = sin(θ)

R[1,1] = cosφ*cosθ*cosψ - sinφ*sinψ
R[1,2] = sinφ*cosθ*cosψ + cosφ*sinψ
R[1,3] = -sinθ*cosψ

R[2,1] = -cosφ*cosθ*sinψ - sinφ*cosψ
R[2,2] = -sinφ*cosθ*sinψ + cosφ*cosψ
R[2,3] = sinθ*sinψ

R[3,1] = cosφ*sinθ
R[3,2] = sinφ*sinθ
R[3,3] = cosθ

return R

end

``````

returns an error as `R[1,1]` seemingly cannot be indexed. I’m also unable to create a diagonal matrix, or an identity matrix (`eye(SMatrix{3,3})` suggested in the docs fails for me, with “eye undefined”).

Basically, I have no idea how to work with StaticArrays – I assumed it was going to be the exact same code and syntax as regular arrays, only they’d have a fixed immutable size, but clearly I got this wrong.

Any clarifications or suggestions welcome!

`SMatrix` is immutable, but you can use `MMatrix` like this:

``````julia> function euler(φ, θ, ψ)
R = MMatrix{3,3,Float64}(undef)
cosφ = cos(φ); cosψ = cos(ψ); cosθ = cos(θ)
sinφ = sin(φ); sinψ = sin(ψ); sinθ = sin(θ)

R[1,1] = cosφ*cosθ*cosψ - sinφ*sinψ
R[1,2] = sinφ*cosθ*cosψ + cosφ*sinψ
R[1,3] = -sinθ*cosψ

R[2,1] = -cosφ*cosθ*sinψ - sinφ*cosψ
R[2,2] = -sinφ*cosθ*sinψ + cosφ*cosψ
R[2,3] = sinθ*sinψ

R[3,1] = cosφ*sinθ
R[3,2] = sinφ*sinθ
R[3,3] = cosθ

return SMatrix(R)
end
``````
1 Like

A couple of things. Firstly, as @foobar_lv2 said, you cannot mutate an `SMatrix` after it has been constructed; but `R` above isn’t an `SMatrix`, it’s a `DataType`, it’s like writing `x = Int` instead of `x = 3`.

You can use an `MMatrix` that you mutate, but that doesn’t seem necessary here, and might be slower. You can use an `SMatrix`, but you must create it ‘in one go’:

``````R = SMatrix{3,3,Float64}(cosφ*cosθ*cosψ - sinφ*sinψ,
sinφ*cosθ*cosψ + cosφ*sinψ, -sinθ*cosψ,
-cosφ*cosθ*sinψ - sinφ*cosψ, -sinφ*cosθ*sinψ + cosφ*cosψ,
sinθ*sinψ, cosφ*sinθ, sinφ*sinθ, cosθ)
``````

If you want to preserve readability, you can try this:

``````S = SMatrix{3,3,Float64}(
#= [1,1] =#  cosφ*cosθ*cosψ - sinφ*sinψ,
#= [1,2] =#  sinφ*cosθ*cosψ + cosφ*sinψ,
#= [1,3] =#  -sinθ*cosψ,

#= [2,1] =#  -cosφ*cosθ*sinψ - sinφ*cosψ,
#= [2,2] =#  -sinφ*cosθ*sinψ + cosφ*cosψ,
#= [2,3] =#  sinθ*sinψ,

#= [3,1] =#  cosφ*sinθ,
#= [3,2] =#  sinφ*sinθ,
#= [3,3] =#  cosθ)
``````

(This works in the REPL but not in Juno, for some reason.)

2 Likes

The `eye` function has been removed from Julia, that’s not just a StaticArrays thing. It’s sort of a wasteful construct, to allocate a full matrix when you only have a constant along the diagonal. Instead you are encouraged to use the `UniformScaling` construct, `I`, from the LinearAlgebra stdlib.

For small StaticArrays, though, it might very occasionally be useful to actually instantiate one. Then you do:

``````using LinearAlgebra
R = MMatrix{3,3,Float64}(I)
``````

You will run into the same issue if you want to create an identity matrix for normal arrays. In that case use

``````R = Matrix{Float64}(I, 3, 3)
``````
1 Like

The “must create in one go” used to annoy me to no end, until I finally saw the light: The compiler is really good at eliding heap allocs. Nowadays I just make an `ref=Base.RefValue{T}()` and then `unsafe_store!` into its `pointer_from_objref`. Even Base is using that pattern in `ReinterpretArray`, and `MArray` is implemented the same way. That’s just the current interface for `alloca` / mutation of immutables / local vars that have their address taken.

But yes, even if both compile to the same, I agree that this specific example looks clearer when written in one go.

Ah, I didn’t see that you turned your `MMatrix` into an `SMatrix` at the end there.

Note that https://github.com/FugroRoames/Rotations.jl already has this functionality and is based on StaticArrays.

3 Likes