Replacing normal arrays with static arrays

New to Julia
1

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!

2

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
3

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
4

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
5

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.

6

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

7

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

3 Likes