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)