traceSpinMultiply(const PMatrix<T1,N,C>& l, const PMatrix<T2,N,C>& r)
traceSpinMultiply(const PScalar<T1>& l, const PMatrix<T2,N,C>& r)
traceSpinMultiply(const PMatrix<T1,N,C>& l, const PMatrix<T2,N,C>& r){  typename BinaryReturn<PMatrix<T1,N,C>, PMatrix<T2,N,C>, FnTraceSpinMultiply>::Type_t  d;  for(int i=0; i < N; ++i)    for(int j=0; j < N; ++j)    {      d.elem(i,j) = traceSpinMultiply(l.elem(i,0), r.elem(0,j));      for(int k=1; k < N; ++k)	d.elem(i,j) += traceSpinMultiply(l.elem(i,k), r.elem(k,j));    }  return d;}