function [B1,C1,B2,C2] = NearestKP(A,p,q,opts)
% A is a p-by-p block matrix with q-by-q blocks.
% Computes B1 and B2 (both p-by-p) and C1 and C2 (both q-by-q) so that
% norm(A - kron(B1,C1) - kron(B2,C2),'fro') is minimized.

% Set up the tilde matrix in sparse format...
Atilde = sparse(0,0);
mu = 0;
for j=1:p
cols = 1+(j-1)*q:j*q;
for i=1:p
rows = 1+(i-1)*q:i*q;
Aij = A(rows,cols);
v = Aij(:)';
Atilde = [Atilde;v];
end
end
% Compute its two biggest singular triples..,
[U,S,V] = svds(Atilde,2,'L',opts);
% Reshape to get the required matrices...
B1 = sqrt(S(1,1))*reshape(full(U(:,1)),p,p);
C1 = sqrt(S(1,1))*reshape(full(V(:,1)),q,q);
B2 = sqrt(S(2,2))*reshape(full(U(:,2)),p,p);
C2 = sqrt(S(2,2))*reshape(full(V(:,2)),q,q);