Fast `diag(A' * B * A)`

[Mason]

@mikmoore

Shouldn’t this be += and initialized to D = zeros(V, size(A,2))? Or maybe use a temporary variable Di = zero(V) to do the accumulation before saving to D[i], if the compiler won’t hoist that on its own…

As written, the loop only does “work” on the final pair of j,k. The compiler may not be smart enough to have caught on, so it may not have a performance impact, but I’m suspicious that this function manages to be so much faster than even B*A despite that it should be doing “more” calculations. The bulk of that is likely the small dimension k=10 giving a meaningful opportunity to beat BLAS.

Crap, yeah you’re right that I should have done += on zeros. The good news is that it does not result in a big performance loss. Updating my post yet again.