Botched the title a bit. Oh well.
On Friday, February 26, 2016 at 9:17:54 AM UTC+1, Kristoffer Carlsson wrote:
>
> 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
>
