function P2()
% Checks out ArrowLS
clc
% Example 1.
m = 1000; n = 100; p = 50; rho = .1;
A = {sprandn(m,n,rho), sprandn(m,n,rho), sprandn(m,n,rho), sprandn(m,n,rho)};
C = {sprandn(m,p,rho), sprandn(m,p,rho), sprandn(m,p,rho), sprandn(m,p,rho)};
d = {randn(m,1),randn(m,1),randn(m,1),randn(m,1)};

tic
xLSCVL = ArrowLSCVL(A,C,d);
tCVL = toc;
tic
xLSYOU = ArrowLS(A,C,d)
tYOU = toc;
tRatio = tYOU/tCVL;
err = norm(xLSYOU - xLSCVL);
fprintf('Example 1. norm(xLSYOU-xLSCVL) = %10.3e    tYou/tCVL = %6.3f \n',err,tRatio)

% Example 2.
% Same structure as example 1, only now A{i},:1) is dense. This will create
% a massive amount of fill-in in R
m = 1000; n = 100; p = 50; rho = .1;
A = {sprandn(m,n,rho), sprandn(m,n,rho), sprandn(m,n,rho), sprandn(m,n,rho)};
A{1}(:,1) = ones(m,1); A{2}(:,1) = ones(m,1); A{3}(:,1) = ones(m,1); A{4}(:,1) = ones(m,1);
C = {sprandn(m,p,rho), sprandn(m,p,rho), sprandn(m,p,rho), sprandn(m,p,rho)};
d = {randn(m,1),randn(m,1),randn(m,1),randn(m,1)};
tic
xLSCVL = ArrowLSCVL(A,C,d);
tCVL = toc;
tic
xLSYOU = ArrowLS(A,C,d);
tYOU = toc;
tRatio = tYOU/tCVL;
err = norm(xLSYOU - xLSCVL);
fprintf('Example 2. norm(xLSYOU-xLSCVL) = %10.3e    tYou/tCVL = %6.3f \n',err,tRatio)

  function xLS = ArrowLSCVL(A,C,d)
% Solves the 4-by-5 block arrow LS problem min{Atilde*x = b}
% where 
%   Atilde = A{1}   0     0     0    C{1}   b =   d{1}
%             0    A{2}   0     0    C{2}         d{2}
%             0     0    A{3}   0    C{3}         d{3}
%             0     0     0    A{4}  C{4}         d{4} 

% Set up the block matrix in sparse format...
[m,n] = size(A{1});
[m,p] = size(C{1});
% For each A-block and C-block we mus figure out where their entries are
% located in the big matrix..
i = []; j = []; s = [];
for k=1:4
   % record ijs information for block A{k}..
   [ik,jk,sk] = find(A{k});
   i = [i; ik+(k-1)*m];
   j = [j; jk+(k-1)*n];
   s = [s; sk];
   % record ijs information for block C{k}...
   [ik,jk,sk] = find(C{k});
   i = [i;ik+(k-1)*m];
   j = [j;jk+4*n];
   s = [s;sk];
end
% Now we can set up the matrix...
Atilde = sparse(i,j,s,4*m,4*n+p);
% Use colamd because QR of Atilde(:,piv) may have a much sparser R than
% QR of A..
piv = colamd(Atilde);
Atilde = Atilde(:,piv);
% Solve min norm(Atilde*x - b} using the method of seminormal equations.
% Afterwards, permute solution to get solution to the original problem...
b = [d{1};d{2};d{3};d{4}];
R = qr(Atilde,0);
xSN = R\(R'\(Atilde'*b));
r = b - Atilde*xSN;
e = R\(R'\(Atilde'*r));
xLS = xSN + e;
% Undo the column permutation...
xLS(piv) = xLS;





               