Performance boosts by using matrix calculations

[phbasler]

I was thinking about how to reduce our runtime for the minimization of the models.
A lot of computation time goes into the calculation of the mass matrix at a given VEV. There we have expressions like

$$ \Lambda_{ij} = \frac{1}{2} L^{ijkl} \omega_k \omega_l $$

Here we sum over k and l for each value of (i,j).
I was therefore thinking of something like this

Before the first set of (i,j) we define

$$ \Omega_{kl} = \omega_k \omega_l $$

as a matrix. Then we could define the matrix

$$ (\tilde{L}{ij}){kl} = \Lambda_{ijkl} $$

and write the loop as

$$ L^{ijkl} \omega_k \omega_l = \mathrm{Tr}( \tilde{L}{ij} \Omega^T ) = \mathrm{Tr}( \tilde{L}{ij} \Omega) $$

This would remove the nested loop over k and l and replace it one matrix multiplication and trace calculation which are optimized by eigen.

What do you think of it?

[phbasler]

According to eigens lazy evaluation eigen only calculates the diagonal elements of the matrix product if we call (A*B).diagonal()