function TestPermKron()
clc
rand('seed',0)
% Generate a small example that can be checked explicitly and displayed...
m = 3; n = 4; 
[x,u] = sort(rand(1,m)); Im = eye(m,m); P = Im(:,u)
[y,v] = sort(rand(1,n)); In = eye(n,n); Q = In(:,v)
R0 = kron(P,Q)
w = PermKron(u,v)
Imn = eye(m*n,m*n);
R = Imn(:,w);
err1 = norm(R - R0)

% A small example that can be checked explicitly and displayed...
m = 300; n = 400; 
[x,u] = sort(rand(1,m)); 
[y,v] = sort(rand(1,n)); 
wCVL = PermKronCVL(u,v);
wYou = PermKron(u,v);
err2 = norm(wCVL-wYou)



  function w = PermKronCVL(u,v)
% u is a permutation of 1:m that encodes an m-by-m permutation P
% v is a permutation of 1:n that encodes an n-by-n permutation Q
% w is a permutation of 1:m*n that encodes R = kron(P,Q)
m = length(u);
n = length(v);
w = [];
for k=1:m
    % In block row k, put a copy of Q in block column u(k).
    % This block starts in column s+1 where
    s = (u(k)-1)*n;
    % Update w so that it includes the indices associated with this block
    % row...
    w = [w s+v];
end


