%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                                      %
% This sample code demonstrates how to manipulate      %
% Kronecker products of matrices without using a lot   %
% of storage.                                          %
%                                                      %
% It solves the linear system  K f = g                 %
% where K = kron(A,B)                                  %
% using the singular value decompositions of A and B.  %
%                                                      %
% Notation:  F is a square array.  f contains these    %
%                same entries, stacked by column.      %
%            G is a square array.  g contains these    %
%                same entries, stacked by column.      %
%                                                      %
% projdemo.m     Dianne P O'Leary 09/02                %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                                      %
% Initialize the data.  The problem is small enough    %
% so that you can print everything out and look at it. %
%                                                      %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

n = 3;
F = rand(n,n);
A = rand(n,n);
B = rand(n,n) + eye(n);


f = reshape(F,n*n,1);
K = kron(A,B);
g = K * f;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                                      %
% We could solve the problem by saying fcomp = K \ g,  %
% but instead we will take advantage of the Kronecker  %
% structure to avoid ever needing K or its factors.    %
%                                                      %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                                      %
% Here is another way to compute the right-hand side.  %
%                                                      %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

G = B * F * A';

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                                      %
% Notice that it gives the same data, so the           %
% following norm prints as a very small number.        %
%                                                      %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Gg = reshape(g,n,n);
disp('norm(G - Gg)')
disp(norm(G - Gg, 'fro'))

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                                      %
% Let's solve the problem using SVD's of A and B,      %
% as a model for understanding the data manipulation   %
% in your homework.                                    %
%                                                      %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


[Ua, Sa, Va] = svd(A);
[Ub, Sb, Vb] = svd(B);

Ghat = Ub'*G*Ua;

S = diag(Sb)*(diag(Sa))';

Fhat = Ghat ./ S;
Fc = Vb*Fhat*Va';

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                                      %
% Notice that it gives almost the exact answer, so the %
% following norm prints as a very small number.        %
%                                                      %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

disp('norm(F - Fc)')
disp(norm(F - Fc, 'fro'))

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                                      %
% For completeness, we verify that                     %
% kron(Ua,Ub) S kron(Va,Vb) is the SVD of K.           %
%                                                      %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


disp('norm(K-svd(K))')
disp(norm(K - kron(Ua,Ub)*diag(reshape(S,n*n,1))*kron(Va,Vb)', 'fro'))
disp('norm(UU^{T}-I)')
disp(norm(kron(Ua,Ub)*kron(Ua,Ub)' - eye(n*n),'fro'))
disp('norm(VV^{T}-I)')
disp(norm(kron(Va,Vb)*kron(Va,Vb)' - eye(n*n),'fro'))