function pi = permute2(mA, mC, mD)
% pi = permute2(mA, mC, mD)
%
% Define a rank mA * (mC + mD) permutation matrix Pi so that
%   Pi ( A \kron [C;D] ) = [ A \kron C ; A \kron D ]
% Store Pi as a permutation vector pi, such that
%   Y = Pi X   <==> Y = X(pi, :)
%   Y = Pi' X  <==> Y(pi,:) = X

piv1 = zeros(mA * mC, 1);
piv2 = zeros(mA * mD, 1);
i1 = 0;
i2 = 0;
for k = 0 : (mA-1)
   piv1(i1 + (1:mC)) = k * (mC+mD) + (1:mC);
   i1 = i1 + mC;
   piv2(i2 + (1 : mD)) = k * (mC+mD) + mC + (1:mD);
   i2 = i2 + mD;
end

pi = [piv1; piv2];
