Hi all: looking at some piece of a code I am developing, I wanted to make sure that there is no type stability issues. So, in the process of reviewing of my code I found a type stability issue when using StaticArrays. Below, I provide three working codes that exemplify the issue.

The first one makes @inferred to throw an error about the type stability issue:

ERROR: return type MArray{Tuple{2,2},Float64,2,4} does not match inferred return type Any

Stacktrace:

[1] top-level scope at none:0

(I had to copy and paste this code in the REPL as running it as a file didnâ€™t write the error to the REPL.)

```
using StaticArrays
using Test
function foo(nelementnodes::Int, le::Float64, area::Float64, Ey::Float64)
xi = gausspoints1d()
w = gaussweights1d()
Ke = MMatrix{nelementnodes,nelementnodes}(zeros(nelementnodes,nelementnodes))
@inbounds for i = 1:10
dNdx = shapefunder1d(nelementnodes, xi[i])
Ke += dNdx * dNdx' * w[i]
end
Ke *= Ey * area * 2.0 / le
return Ke
end
function shapefunder1d(nelementnodes::Int, xi::Float64)
dNdxi = MVector{nelementnodes}(zeros(nelementnodes))
if nelementnodes == 2
dNdxi[1] = -0.5
dNdxi[2] = 0.5
elseif nelementnodes == 3
dNdxi[1] = -0.5 + xi
dNdxi[2] = -2.0 * xi
dNdxi[3] = 0.5 + xi
else
throw("shapefunder1d(): illegal value of nelementnodes")
end
return dNdxi
end
function gausspoints1d()
x = MVector{10}(zeros(10))
x[1] = - 0.973906528517171720077964012084
x[2] = - 0.865063366688984510732096688423
x[3] = - 0.679409568299024406234327365115
x[4] = - 0.433395394129247290799265943166
x[5] = - 0.148874338981631210884826001130
x[6] = 0.148874338981631210884826001130
x[7] = 0.433395394129247290799265943166
x[8] = 0.679409568299024406234327365115
x[9] = 0.865063366688984510732096688423
x[10] = 0.973906528517171720077964012084
end
function gaussweights1d()
w = MVector{10}(zeros(10))
w[1] = 0.0666713443086881375935688098933
w[2] = 0.149451349150580593145776339658
w[3] = 0.219086362515982043995534934228
w[4] = 0.269266719309996355091226921569
w[5] = 0.295524224714752870173892994651
w[6] = 0.295524224714752870173892994651
w[7] = 0.269266719309996355091226921569
w[8] = 0.219086362515982043995534934228
w[9] = 0.149451349150580593145776339658
w[10] = 0.0666713443086881375935688098933
end
@inferred foo(2, 400.0, 0.8, 1000000.0)
```

I believe the problem is in the branching inside the function shapefunder1d(). I suspect this because when I remove the â€śifâ€ť statement, the @inferred correctly infers the type:

2Ă—2 MArray{Tuple{2,2},Float64,2,4}:

666.713 -666.713

-666.713 666.713

This is the code with the â€śifâ€ť statement removed from the shapefunder1d() function:

```
using StaticArrays
using Test
function foo(nelementnodes::Int, le::Float64, area::Float64, Ey::Float64)
xi = gausspoints1d()
w = gaussweights1d()
Ke = MMatrix{2,2}(zeros(2,2))
@inbounds for i = 1:10
dNdx = shapefunder1d(nelementnodes, xi[i])
Ke += dNdx * dNdx' * w[i]
end
Ke *= Ey * area * 2.0 / le
return Ke
end
function shapefunder1d(nelementnodes::Int, xi::Float64)
dNdxi = MVector{2}(zeros(2))
dNdxi[1] = -0.5
dNdxi[2] = 0.5
return dNdxi
end
function gausspoints1d()
x = MVector{10}(zeros(10))
x[1] = - 0.973906528517171720077964012084
x[2] = - 0.865063366688984510732096688423
x[3] = - 0.679409568299024406234327365115
x[4] = - 0.433395394129247290799265943166
x[5] = - 0.148874338981631210884826001130
x[6] = 0.148874338981631210884826001130
x[7] = 0.433395394129247290799265943166
x[8] = 0.679409568299024406234327365115
x[9] = 0.865063366688984510732096688423
x[10] = 0.973906528517171720077964012084
end
function gaussweights1d()
w = MVector{10}(zeros(10))
w[1] = 0.0666713443086881375935688098933
w[2] = 0.149451349150580593145776339658
w[3] = 0.219086362515982043995534934228
w[4] = 0.269266719309996355091226921569
w[5] = 0.295524224714752870173892994651
w[6] = 0.295524224714752870173892994651
w[7] = 0.269266719309996355091226921569
w[8] = 0.219086362515982043995534934228
w[9] = 0.149451349150580593145776339658
w[10] = 0.0666713443086881375935688098933
end
@inferred foo(2, 400.0, 0.8, 1000000.0)
```

Finally, I tested the first code (the one with the â€śifâ€ť statement) but using standard Julia arrays. This has no type stability problems. The @inferred correctly infers the type:

2Ă—2 Array{Float64,2}:

666.713 -666.713

-666.713 666.713

The code using standard Julia arrays is:

```
using Test
function foo(nelementnodes::Int, le::Float64, area::Float64, Ey::Float64)
xi = gausspoints1d()
w = gaussweights1d()
Ke = zeros(nelementnodes,nelementnodes)
@inbounds for i = 1:10
dNdx = shapefunder1d(nelementnodes, xi[i])
Ke += dNdx * dNdx' * w[i]
end
Ke *= Ey * area * 2.0 / le
return Ke
end
function shapefunder1d(nelementnodes::Int, xi::Float64)
dNdxi = zeros(nelementnodes)
if nelementnodes == 2
dNdxi[1] = -0.5
dNdxi[2] = 0.5
elseif nelementnodes == 3
dNdxi[1] = -0.5 + xi
dNdxi[2] = -2.0 * xi
dNdxi[3] = 0.5 + xi
else
throw("shapefunder1d(): illegal value of nelementnodes")
end
return dNdxi
end
function gausspoints1d()
x = zeros(10)
x[1] = - 0.973906528517171720077964012084
x[2] = - 0.865063366688984510732096688423
x[3] = - 0.679409568299024406234327365115
x[4] = - 0.433395394129247290799265943166
x[5] = - 0.148874338981631210884826001130
x[6] = 0.148874338981631210884826001130
x[7] = 0.433395394129247290799265943166
x[8] = 0.679409568299024406234327365115
x[9] = 0.865063366688984510732096688423
x[10] = 0.973906528517171720077964012084
end
function gaussweights1d()
w = zeros(10)
w[1] = 0.0666713443086881375935688098933
w[2] = 0.149451349150580593145776339658
w[3] = 0.219086362515982043995534934228
w[4] = 0.269266719309996355091226921569
w[5] = 0.295524224714752870173892994651
w[6] = 0.295524224714752870173892994651
w[7] = 0.269266719309996355091226921569
w[8] = 0.219086362515982043995534934228
w[9] = 0.149451349150580593145776339658
w[10] = 0.0666713443086881375935688098933
end
@inferred foo(2, 400.0, 0.8, 1000000.0)
```

Any idea how to fix the type stability issue in the first code without removing the â€śifâ€ť statement as shapefunder1d() will have to compute arrays of different sizes depending on â€śnelementnodesâ€ť? Also, I would like to insist on using StaticArrays as compared to standard Julia arrays they perform better for small sized arrays.

Thanks.

EDIT: SOLUTION

Taking the advices of rdeits and Kristoffer.carlsson I modified the code that uses StaticArrays. No more type stability issues.

Solution 1:

```
using StaticArrays
using Test
function foo(::Val{nelementnodes}, le::Float64, area::Float64, Ey::Float64) where nelementnodes
xi = gausspoints1d()
w = gaussweights1d()
Ke = zeros(SMatrix{nelementnodes,nelementnodes})
@inbounds for i = 1:10
dNdx = shapefunder1d(Val(nelementnodes), xi[i])
Ke += dNdx * dNdx' * w[i]
end
Ke *= Ey * area * 2.0 / le
return Ke
end
function shapefunder1d(::Val{2}, xi::Float64)
return @SVector [-0.5, 0.5]
end
function shapefunder1d(::Val{3}, xi::Float64)
return @SVector [-0.5 + xi, -2.0 * xi, 0.5 + xi]
end
function gausspoints1d()
return @SVector [- 0.973906528517171720077964012084,
- 0.865063366688984510732096688423,
- 0.679409568299024406234327365115,
- 0.433395394129247290799265943166,
- 0.148874338981631210884826001130,
0.148874338981631210884826001130,
0.433395394129247290799265943166,
0.679409568299024406234327365115,
0.865063366688984510732096688423,
0.973906528517171720077964012084]
end
function gaussweights1d()
return @SVector [0.0666713443086881375935688098933,
0.149451349150580593145776339658,
0.219086362515982043995534934228,
0.269266719309996355091226921569,
0.295524224714752870173892994651,
0.295524224714752870173892994651,
0.269266719309996355091226921569,
0.219086362515982043995534934228,
0.149451349150580593145776339658,
0.0666713443086881375935688098933]
end
nelementnodes = 2
@inferred foo(Val(nelementnodes), 400.0, 0.8, 1000000.0)
nelementnodes = 3
@inferred foo(Val(nelementnodes), 400.0, 0.8, 1000000.0)
```

Solution 2:

```
using StaticArrays
using Test
function foo(Ke::MMatrix{2,2,Float64}, le::Float64, area::Float64, Ey::Float64)
xi = gausspoints1d()
w = gaussweights1d()
@inbounds for i = 1:10
dNdx = shapefunder1d(2, xi[i])
Ke .+= dNdx * dNdx' * w[i]
end
Ke .*= Ey * area * 2.0 / le
end
function foo(Ke::MMatrix{3,3,Float64}, le::Float64, area::Float64, Ey::Float64)
xi = gausspoints1d()
w = gaussweights1d()
@inbounds for i = 1:10
dNdx = shapefunder1d(3, xi[i])
Ke .+= dNdx * dNdx' * w[i]
end
Ke .*= Ey * area * 2.0 / le
end
function shapefunder1d(nelementnodes::Int, xi::Float64)
nelementnodes == 2 && return @SVector [-0.5, 0.5]
nelementnodes == 3 && return @SVector [-0.5 + xi, -2.0 * xi, 0.5 + xi]
end
function gausspoints1d()
return @SVector [- 0.973906528517171720077964012084,
- 0.865063366688984510732096688423,
- 0.679409568299024406234327365115,
- 0.433395394129247290799265943166,
- 0.148874338981631210884826001130,
0.148874338981631210884826001130,
0.433395394129247290799265943166,
0.679409568299024406234327365115,
0.865063366688984510732096688423,
0.973906528517171720077964012084]
end
function gaussweights1d()
return @SVector [0.0666713443086881375935688098933,
0.149451349150580593145776339658,
0.219086362515982043995534934228,
0.269266719309996355091226921569,
0.295524224714752870173892994651,
0.295524224714752870173892994651,
0.269266719309996355091226921569,
0.219086362515982043995534934228,
0.149451349150580593145776339658,
0.0666713443086881375935688098933]
end
nelementnodes = 2
Ke = zeros(MMatrix{nelementnodes,nelementnodes})
@inferred foo(Ke, 400.0, 0.8, 1000000.0)
nelementnodes = 3
Ke = zeros(MMatrix{nelementnodes,nelementnodes})
@inferred foo(Ke, 400.0, 0.8, 1000000.0)
```