tion)</span>
  X = reshape(1:D1*D1,[D1 D1]); XT = X';              <span class="comment">% vectorized indices</span>
  I=0*X; I(i,i)=1; ii=X(I==1)'; I=0*X; I(k,k)=1; kk=X(I==1)';
  I=0*X; I(i,k)=1; ik=X(I==1)'; ki=XT(I==1)';

  Mdm(k,:)  = mdm*Mdm(i,:) + mds*Sdm(ii,:);                    <span class="comment">% chainrule</span>
  Mds(k,:)  = mdm*Mds(i,:) + mds*Sds(ii,:);
  Sdm(kk,:) = sdm*Mdm(i,:) + sds*Sdm(ii,:);
  Sds(kk,:) = sdm*Mds(i,:) + sds*Sds(ii,:);
  dCdm      = Cdm*Mdm(i,:) + Cds*Sdm(ii,:);
  dCds