  function P1Grade()
% Checks out Qnn
% 10 points
% You have to vectorize the computation of the subscript and value vectors
% that are passed to sparse. Double loop implementations are very slow
% and your score will be around 5/10 if tYOU/tCVL is >=10

% Small example for correctness...
n = 10;
Q = Qnn(n);
err1 = norm(Q'*Q - speye(n^2,n^2),1);
err2 = norm(Q*Q' - speye(n^2,n^2),1);
clc
fprintf('n = %1d\n',n)
fprintf('|| Q''Q - I|| = %10.3e\n',err1)
fprintf('|| QQ'' - I|| = %10.3e\n',err2)
fprintf('|| Q - QCVL || = %10.3e\n',norm(Q-QnnCVL(n),1))

% Large example for efficiency
n = 500;
tic
Q = Qnn(n);
tYOU = toc;
tic
Q = QnnCVL(n);
tCVL = toc;
fprintf('\nn = %1d\n',n)
fprintf('tYOU/tCVL = %6.3f\n',tYOU/tCVL)


  function Q = QnnCVL(n)
% n is a positive integer
% Q is the n^2-by-n^2 orthogonal matrix Q_{n,n} in sparse format.

N1 = n*(n+1)/2;  % The number of symmetric basis vectors
N2 = n*(n-1)/2;  % The number of skw-symmetric basis vectors
N = n*n;

% The symmetric rank-1 columns...
% Q(i1(k),j1(k)) = 1, k=1:n
i1 = zeros(1,n); i1(1)=1;
j1 = zeros(1,n);j1(1)=1;
for k=1:n-1
    i1(k+1)=i1(k)+n+1;
    j1(k+1)=j1(k)+n+1-k;   
end
s1=ones(1,n);

% The rank-2 columns.
% The symmetric ones have nonzeros at 
%         Q(i2(k),j2(k)) = s2(k), k=1:(2*N2)
% The skew-symmetric ones have nonzeros at 
%         Q(i2(k),j3(k)) = s3(k), k=1:(2*N2)
s2= ones(1,2*N2)/sqrt(2);
v = ones(1,n)/sqrt(2);
L = N1;
i2=zeros(1,2*N2);j2=i2;j3=i2;s3=i2;
m=0;
for k=1:n-1
    ran=m+(1:2*(n-k));
    colRange = j1(k)+1:j1(k+1)-1;
    i2(ran) = [i1(k)+1:i1(k)+n-k   i1(k)+(n:n:(n-k)*n) ];
    j2(ran) = [colRange            colRange] ; 
    F = L+1;
    L = F+(n-k-1);
    j3(ran) = [F:L F:L];
    s3(ran) = [ones(1,n-k) -ones(1,n-k)];
    m = m+2*(n-k);
end
Q = sparse([i1 i2 i2]',[j1 j2 j3]',[s1 s2 s3/sqrt(2) ]',N,N);



  function Q = QnnCVL1(n)
% n is a positive integer
% Q is the n^2-by-n^2 orthogonal matrix Q_{n,n} in sparse format.
% This is a super inefficient implementation. 

alpha = 1/sqrt(2);
I = speye(n,n);
% Allocate space. n columns have one nonzero, all the rest have two
% nonzeros...
Q = sparse([],[],[],n^2,n^2,2*n^2-n);
% Indices that point to the next columns to fill in Q...
p = 0;
r = n*(n+1)/2;
for i=1:n
    p = p+1;
    % Store vec(S_{ii}) in column p...
    ei = I(:,i);
    Q(:,p) = reshape(ei*ei',n^2,1);
    for j=i+1:n
        % Store vec(alpha*S_{ij}) in column p...
        p = p+1;
        ej = I(:,j);
        vij=ei*ej';
        % Q(:,p) = alpha*reshape(ei*ej'+ej*ei',n^2,1);
        Q(:,p) = alpha*reshape(vij+vij',n^2,1);
        % Store vec(alpha*T_{ij}) in column r...
        r = r+1;
        % Q(:,r) = alpha*reshape(ei*ej'-ej*ei',n^2,1);
        Q(:,r) = alpha*reshape(vij'-vij,n^2,1);
    end
end
