  function ShowSparseLU(A,C,thresh)
% A is a length-p cell array where A{i} is an r-by-r symmetric positive
%    definite matrix,i=1:p. These are the diagonal blocks.
% C is a length-(p-1) cell array where C{i} is and r-by-r matrix, i=1:p-1.
%    These are the subdiagonal blocks.
% thresh satisfies 0<=thresh<=1.
% Sets up a sparse representation of the symskew matrix T defined by
% A and C and computes [L,U,p] = lu(T,thresh). 
% Produces two figures:
% Figure1: Displays the sparsity of L using spy and reports the value
%    of thresh and the quotient nz(L)/nz(T) in the title.
% Figure2: Displays the sparsity of U using spy and reports the value
%    of thresh and the quotient nz(U)/nz(T) in the title. 

p = length(A);
[r,r] = size(A{1});
n = p*r;

% Set up T in sparse format...
T = sparse(n,n);
for i=1:p
    idx = 1+(i-1)*r:i*r;
    T(idx,idx) = A{i};
    if i<p
        T(idx+r,idx) = C{i};
        T(idx,idx+r) = -C{i}';
    end
end

% Sparse LU...
[L,U,p] = lu(T,thresh);

% Number of nonzeros...
[i,j,Tvals] = find(T); nonZerosInT = length(Tvals);
[i,j,Lvals] = find(L); nonZerosInL = length(Lvals);
[i,j,Uvals] = find(U); nonZerosInU = length(Uvals);
Ratio1 = nonZerosInL/nonZerosInT;
Ratio2 = nonZerosInU/nonZerosInT;

% Display....
figure
spy(L)
title(sprintf('thresh = %5.2f   nz(L)/nz(T) = %5.2f',thresh,Ratio1))
figure
spy(U)
title(sprintf('thresh = %5.2f   nz(U)/nz(T) = %5.2f',thresh,Ratio2))

% Number of nonzeros in nz(U)/nz(T) increases from .5 tp .6 as thresh
% increase from 0 to 1.

%-1 for minor set up T mistakes

        