  function P1Grade()
% Checks out NewDiag

% Generate a random example...
n = ceil(100*rand());
A = randn(n,n);
v = randperm(n);
k = floor(20*rand());
% Solve the problem by explicitly forming P and B...
P = eye(n,n);
P = P(v,:);
B = (P^k)'*A*P^k;
d0 = diag(B);
% Compare d0 with the vector returned by NewDiag..
d = NewDiag(A,v,k);
clc
fprintf('n = %2d, k = %2d, sum(abs(d-d0)) = %5.0f\n',n,k,sum(d-d0))

  function d = NewDiagCVL(A,v,k)
% A is nxn, v is a permutation of 1:n, and k is a nonnegative integer.
% d = diag(B) where B = {P^k]'*A*P^k and P = I(v,:), I = eye(n,n).
d = diag(A);
for i=1:k
    d(v) = d;
end

% Key fact: If B = P'*A*P, then diag(B) = P'*diag(A)

% -3 if you move O(n^2) data around
% -4 if you explicitly form P and have things like P*A or P^k*A
