Hello everyone.
To avoid the XY problem I will first describe my actual use case. I am
writing a library for tensors to be used in continuum mechanics. These are
in dimensions 1,2,3 and typically of order 1,2,4.
Since tensor expressions are usually chained together like (A ⊗ B) * C ⋅ a
I want to make these stack allocated so I can get good performance without
making a bunch of mutating functions and having to preallocate everything
etc.
The way you stack allocate vector like objects in julia right now seems to
almost exclusively be with an immutable of a tuple. So my data structure
right now looks similar to:
immutable FourthOrderTensor{T} <: AbstractArray{T, 4}
data::NTuple{81, T}
end
Base.size(::FourthOrderTensor) = (3,3,3,3)
Base.linearindexing{T}(::Type{FourthOrderTensor{T}}) = Base.LinearFast()
Base.getindex(S::FourthOrderTensor, i::Int) = S.data[i]
If you'd like, you can think of the NTuple{81, T} as actually represent a
9x9 matrix.
The largest operations I will need to do is a double contraction between
two fourth order tensors which is equivalent to a 9x9 matrix multiplied
with another 9x9 matrix.
My actual question is, is there a good way to do this without generating a
huge amount of code?
I have looked at for example the FixedSizeArrays package but there they
seem to simply unroll every operation. For a 9x9 * 9x9 matrix
multiplication this unrolling creates a huge amount of code and doesn't
feel right.
I guess what I want is pretty
much https://github.com/JuliaLang/julia/issues/11902 but is there any
workaround as of now?
Thanks!
// Kristoffer
