%SPARSITY Demonstrate effect of sparsity orderings.

%	Copyright (c) 1984-93 by The MathWorks, Inc.

clf
clc
echo on

% This demonstration shows the effect of different orderings
% in a sparse matrix on the time and storage requirements for
% sparse matrix operations.
% We begin with a Harwell-Boeing sparse test matrix, west0479.
% This matrix describes the connections in a model of a
% diffraction column in a chemical plant.

load('west0479.mat')
A = west0479;

% A "spy" plot shows the locations of the nonzero elements.

spy(A)
title('spy(A)')

% Press any key to continue after pauses.
pause

clc

% Here are spy plots of A, A', A*A' and A'*A.

subplot(2,2,1), spy(A),    title('A'),     drawnow
subplot(2,2,2), spy(A'),   title('A'''),   drawnow
subplot(2,2,3), spy(A*A'), title('A*A'''), drawnow
subplot(2,2,4), spy(A'*A); title('A''*A'), drawnow

% The spy plots of A*A' and A'*A show "fill-in", the creation of
% nonzero elements in positions that are zero in the original matrix.
% Storage requirements for sparse matrices are proportional to
% the number of nonzero elements.

pause

clf
clc

% We continue with A*A', which is symmetric and positive definite.

S = A*A';

subplot(2,2,1), spy(S), title('S')

pause

clc

% The computation of the Cholesky factorization creates fill-in.
% Let's time the factorization and count the number of nonzeros
% in the triangular factor.

tic
R = chol(S);
toc
nz = nnz(R)

subplot(2,2,2), spy(R), title('chol(S)')

pause

clc

% The reverse Cuthill-McKee algorithm renumbers the rows and columns
% in an attempt to reduce the bandwidth or profile of the matrix.

p = symrcm(S);
subplot(2,2,3), spy(S(p,p)), title('S(p,p)')

pause

clc

% The fill-in is confined to the band, so the Cholesky factorization
% of the reordered matrix takes less time and less storage.

tic
R = chol(S(p,p));
toc
nz = nnz(R)

subplot(2,2,4), spy(R), title('chol(S(p,p))')

pause

clf
clc

% The column count ordering moves rows and columns with higher
% nonzero count towards the end of the matrix.

q = colperm(S);
subplot(2,2,1), spy(S(q,q)), title('S(q,q)')

pause

clc

% For this example, the column count ordering happens to reduce
% the time and storage for Cholesky factorization, but this
% behavior cannot be expected in general.

tic
R = chol(S(q,q));
toc
nz = nnz(R)

subplot(2,2,2), spy(R), title('chol(S(q,q))')

pause

clc

% Minimum degree reordering is a powerful, graph theoretic algorithm
% that produces large blocks of zeros in the matrix.

r = symmmd(S);
subplot(2,2,3), spy(S(r,r)), title('S(r,r)')

pause

clc

% These blocks are preserved during the Cholesky factorization.
% As a result, minimum degree requires the least time and storage.

tic
R = chol(S(r,r));
toc
nz = nnz(R)

subplot(2,2,4), spy(R), title('chol(S(r,r))')

pause

% End
echo off
disp('>>')
