  function P2Grade()
% Checks out NearestKP
% 10 points
% Have to be smart in using sparse to set up atilde.
% Lost of zero rows. uOpt is tridiagonal
%3/10 - 7/10 for poor implementations that render bad tYOU/tCVL

clc
% Small example for correctness...
n1 = 20;
n2 = 50;
density = .1;
A = sparse([],[],[],n1*n2,n1*n2,0);
for i=1:n1
    idx = 1+(i-1)*n2:i*n2;
    A(idx,idx) = sprand(n2,n2,density)+speye(n2);
    if i>1
        A(idx-n2,idx) = sprand(n2,n2,density);
    end
    if i<n1
        A(idx+n2,idx) = sprand(n2,n2,density);
    end
end
[BYOU,CYOU] = NearestKP(A,n1,n2);
[BCVL,CCVL] = NearestKPCVL(A,n1,n2);
fprintf('n1 = %3d\nn2 = %3d\n',n1,n2)
fprintf('||BYOU x CYOU -  BCVL x CCVL || = %10.3e\n\n',norm(kron(BYOU,CYOU)-kron(BCVL,CCVL),1))

% Large example for efficiency...
n1 = 100;
n2 = 100;
density = .04;
A = sparse([],[],[],n1*n2,n1*n2,0);
for i=1:n1
    idx = 1+(i-1)*n2:i*n2;
    A(idx,idx) = sprand(n2,n2,density)+speye(n2);
    if i>1
        A(idx-n2,idx) = sprand(n2,n2,density);
    end
    if i<n1
        A(idx+n2,idx) = sprand(n2,n2,density);
    end
end

tic
[BYOU,CYOU] = NearestKP(A,n1,n2);
tYOU = toc;
tic
[BCVL,CCVL] = NearestKPCVL(A,n1,n2);
tCVL = toc;
fprintf('n1 = %3d\nn2 = %3d\n',n1,n2)
fprintf('tYOU/tCVL = %6.3f\n',tYOU/tCVL)

function [B,C] = NearestKPCVL(A,n1,n2)
% B (n1-by-n1) and C (n2-by-n2) minimize  and are both
% in full format.
%
%   AUTHOR: Andrew Loeb, 17/12/2015
%
%   INPUT:
%
%   A - n1*n2 x n1*n2
%           sparse n1-by-n1 block tridiagonal matrix with n2-by-n2 blocks
%
%   OUTPUT:
%   
%   B - n1 x n1
%           component of minimizing norm(A-kron(B,C),’fro’)
%   C - n2 x n2
%           component of minimizing norm(A-kron(B,C),’fro’)
%
%   NOTES:
%           B and C are in full format
%
%   SEE ALSO:
%   

% Initialize variables
I = [];
J = [];
V = [];

% We will construct Abar with a single call to sparse. We need the indices
% for the rows and columns of nonzero entries.
% Examine each tridiagonal block
for i = 1:n1
    for j = max(1,i-1):min(n1,i+1)
        % Pull off n2 x n2 block 
        Aij = A((i-1)*n2 + 1:i*n2, (j-1)*n2 + 1:j*n2);
        % Find the coordinates and nonzero elements
        [row, col, v] = find(Aij);
        % Each submatrix is mapped to a single row of Abar
        I = [I; (i + (j-1)*n1)*ones(length(v),1)];
        % The nonzero elements of Aij are mapped along the row
        J = [J; sub2ind([n2, n2], row, col)];
        % Store the nonzero elements
        V = [V; v];
    end
end

% Construct Abar
Abar = sparse(I, J, V, n1^2, n2^2);

% Find the largest singular value and its associated singular vectors
[U,S,V] = svds(Abar,1);

% Build B and C in full format
B = sqrt(S)*reshape(U,n1,n1);
C = sqrt(S)*reshape(V,n2,n2);

  function [B,C] = NearestKPCVL1(A,n1,n2)
% A (in sparse format) is an n1-by-n1 block tridiagonal matrix with n2-by-n2 blocks.
% B (n1-by-n1) and C (n2-by-n2) minimize norm(A-kron(B,C),'fro') and are both 
% in full format.

% Build Atilde
Atilde = sparse([],[],[],n1*n1,n2*n2,0);
k = 0;
for j=1:n1
    jdx = 1+(j-1)*n2:j*n2;
    for i=1:n1
       idx = 1+(i-1)*n2:i*n2;
       k = k+1;
       if abs(i-j)<=1
           Atilde(k,:) = reshape(A(idx,jdx),1,n2*n2);
       end
    end
end
[u1,sig1,v1] = svds(Atilde,1);
B = sqrt(sig1)*reshape(sparse(u1),n1,n1);
C = sqrt(sig1)*reshape(v1,n2,n2);


  function [B,C] = NearestKP2(A,n1,n2)
% A (in sparse format) is an n1-by-n1 block tridiagonal matrix with n2-by-n2 blocks.
% B (n1-by-n1) and C (n2-by-n2) minimize norm(A-kron(B,C),'fro') and are both 
% in full format.

% Build Atilde
Atilde = zeros(n1*n1,n2*n2);
k = 0;
for j=1:n1
    jdx = 1+(j-1)*n2:j*n2;
    for i=1:n1
       idx = 1+(i-1)*n2:i*n2;
       k = k+1;
       if abs(i-j)<=1
           Atilde(k,:) = reshape(A(idx,jdx),1,n2*n2);
       end
    end
end
[U,S,V] = svd(Atilde);
B = sqrt(S(1,1))*reshape(U(:,1),n1,n1);
C = sqrt(S(1,1))*reshape(V(:,1),n2,n2);