function proj1demo ( )

%% PROJ1DEMO demonstrates how to manipulate Kronecker products.
%
%  Discussion:
%
%    It solves the linear system  K f = g
%    where K = kron(A,B)
%    using the singular value decompositions of A and B.
%
%    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.
%
%  Author:
%
%    Dianne P O'Leary 09/02
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'PROJ1DEMO\n' );
  fprintf ( 1, '  This function shows how to work with Kronecker products and SVD''s.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Suppose K = kronecker_product ( A, B )\n' );
  fprintf ( 1, '  Suppose we want to solve K*f=g\n' );
  fprintf ( 1, '  using the singular value decomposition of K.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  In this example, we choose random values for\n' );
  fprintf ( 1, '  f, A and B, then compute K and g.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Instead of solving K*f=g by taking the SVD of K\n' );
  fprintf ( 1, '  we compute the SVD''s of A and B, which implicitly\n' );
  fprintf ( 1, '  define the SVD of K.  Then we can determine the solution f\n' );
  fprintf ( 1, '  and we never actually formed the SVD of K, and we really\n' );
  fprintf ( 1, '  never had to form K either.\n' );

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

  n = 3;

  fprintf ( 1, '\n' );
  fprintf ( 1, '  The matrices A and B will be %dx%d\n', n, n );

  A = rand(n,n);
  B = rand(n,n) + eye(n);
%
%  "F" is a 3x3 matrix, 
%  "f" is the corresponding vector of length 9.
%
  F = rand(n,n);
  f = reshape(F,n*n,1);

  fprintf ( 1, '\n' );
  fprintf ( 1, '  The matrix K will have order %dx%d\n', n*n, n*n );
%
%  "KRON" is a MATLAB command that forms the Kronecker product of two matrices.
%  Here, we go ahead and actually form K, because N is small. 
%  But for our real problem, K will be too big.  
%  The point we are making here is that we can solve the 
%  problem K*f=g without actually forming K.
%
  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.%
%                                                      %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

  G = B * F * A';

%
% Notice that it gives the same data, so the
% following norm prints as a very small number.
%  (To compare G and g, we make a copy "Gg" of g that looks like a matrix,
%  and look at ||G - Gg||.
%

  Gg = reshape(g,n,n);

  fprintf ( 1, '\n' );
  fprintf ( 1, '  There are two ways to form the right hand side g.\n' );
  fprintf ( 1, '  To show they are the same, here is the norm of their difference:\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  norm(G - g) = %e\n', 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.        %
%                                                      %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Here is the error between the exact solution F\n' );
  fprintf ( 1, '  and the solution we computed:\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  norm(F - Fc) = %e\n', norm ( F - Fc, 'fro' ) );

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

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Now let''s show that we really can compute the SVD of K.\n' );
  fprintf ( 1, '  Does the product U * S * V'' = K?\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Norm ( K - U*S*V'' ) = %e\n', ...
    norm(K - kron(Ua,Ub)*diag(reshape(S,n*n,1))*kron(Va,Vb)', 'fro') );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Is our U matrix orthogonal?\n' );
  fprintf ( 1, '  norm ( U*U''-I ) = %e\n', norm(kron(Ua,Ub)*kron(Ua,Ub)' - eye(n*n),'fro'));
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Is our V matrix orthogonal?\n' );
  fprintf ( 1, '  norm ( U*V''-I ) = %e\n', norm(kron(Va,Vb)*kron(Va,Vb)' - eye(n*n),'fro'));

  fprintf ( 1, '\n' );
  fprintf ( 1, 'PROJ1DEMO:\n' );
  fprintf ( 1, '  Normal end of execution.\n' );

  return
end
