d.elem() = s1.elem(0,0) * s2.elem(0,0);
d.elem() = l.elem(0,0) * r.elem(0,0);
traceSpinMultiply(const PSpinMatrix<T1,N>& l, const PSpinMatrix<T2,N>& r){  typename BinaryReturn<PSpinMatrix<T1,N>, PSpinMatrix<T2,N>, FnTraceSpinMultiply>::Type_t  d;  // The traceSpin is eaten here  d.elem() = s1.elem(0,0) * s2.elem(0,0);  for(int k=1; k < N; ++k)    d.elem() += s1.elem(0,k) * s2.elem(k,0);  for(int j=1; j < N; ++j)    for(int k=0; k < N; ++k)      d.elem() += s1.elem(j,k) * s2.elem(k,j);  return d;}