%  This Matlab file generates a sparse matrix (the 5-point
%  operator on an m x m grid) and reorders it using three
%  algorithms.

%  A better example is given by the Matlab sparse demo, from
%  which much of this code is taken.

%                                Dianne O'Leary December 1996

m=50;
N = m * m;
b = -ones(N,5);
d = [0 1, -1, m, -m];
b(:,1)=4*ones(N,1);
for i=1:m-1,
  b(i*m+1,2)=0;
  b(i*m,3)=0;
end
a = spdiags(b,d,N,N);
S = sparse(a);
pct=100;

%  Press the "Start" button to see a demonstration which   
%  shows that reordering the rows and columns of a         
%  sparse matrix S can affect the time and storage required
%  for a matrix operation such as factoring S into its     
%  Cholesky decomposition, S=L*L'.          

spy(S), title('A Sparse Symmetric Matrix')           
n = nnz(S);                                          

%  A SPY plot shows the nonzero elements in a matrix.             
%  This spy plot shows a SPARSE symmetric positive definite       
%  matrix derived from a portion of the Harwell-Boeing test matrix
%  "west0479", a matrix describing connections in a model         
%  of a diffraction column in a chemical plant.                   
                                                                
title('spy(S)')             
%  print spar50.ps 

%  Now we compute the Cholesky factor L, where S=L*L'.           
%  Notice that L contains MANY more nonzero elements than        
%  the unfactored S, because the computation of the Cholesky     
%  factorization creates "fill-in" nonzeros.  This slows down the
%  algorithm and increases storage cost.                         
                                                               
tic, L = chol(S)'; t(1) = toc;                             
spy(L), title('Cholesky decomposition of S')               
nzz(1) = nnz(L);
%  print spar50.ps -append

%  By reordering the rows and columns of a matrix, it may be
%  possible to reduce the amount of fill-in created by      
%  factorization, thereby reducing time and storage cost.   
%   
%  We will now try three different orderings supported by   
%  MATLAB.                                                  
%  The SYMRCM command uses the reverse Cuthill-McKee             
%  reordering algorithm to move all nonzero elements closer to   
%  the diagonal, reducing the "bandwidth" of the original matrix.
                                                               
p = symrcm(S);                                             
spy(S(p,p)), title('S(p,p) after Cuthill-McKee ordering')  
%  print spar50.ps -append

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

tic, L = chol(S(p,p))'; t(2) = toc;                       
spy(L), title('chol(S(p,p)) after Cuthill-McKee ordering')
nzz(2) = nnz(L);
%  print spar50.ps -append

%  The COLPERM command uses the column count reordering       
%  algorithm to move rows and columns with higher nonzero     
%  count towards the end of the matrix.                       
                                                            
q = colperm(S);                                         
spy(S(q,q)), title('S(q,q) after column count ordering')

%  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, L = chol(S(q,q))'; t(3) = toc;                      
spy(L), title('chol(S(q,q)) after column count ordering')
nzz(3)=nnz(L);
%  print spar50.ps -append

%  The SYMMMD command uses the minimimum degree                 
%  algorithm (a powerful graph-theoretic technique) to produce  
%  large blocks of zeros in the matrix.                         
                                                              
r = symmmd(S);                                            
spy(S(r,r)), title('S(r,r) after minimum degree ordering')
%  print spar50.ps -append

%  The blocks of zeros produced by the minimum degree            
%  algorithm are preserved during the Cholesky factorization.    
%  This can significantly reduce time and storage costs.         
                                                               
tic, L = chol(S(r,r))'; t(4) = toc;                        
spy(L), title('chol(S(r,r)) after minimum degree ordering')
nzz(4)=nnz(L);
%  print spar50.ps -append

disp(        'RESULTS SUMMARY    NonZeroes   Fact.Time')
                                                               
disp(sprintf('original:          %7.0f     %f ',  nzz(1),t(1) ))
disp(sprintf('Cuthill-McKee:     %7.0f     %f ',  nzz(2),t(2) ))
disp(sprintf('Column count:      %7.0f     %f ',  nzz(3),t(3) ))
disp(sprintf('Minimum deg:       %7.0f     %f ',  nzz(4),t(4) ))
                                                             %   
                                                               
%           5 x 5 grid:  N=25 

%     RESULTS SUMMARY    NonZeroes   Fact.Time
%     original:              129     0.001894 
%     Cuthill-McKee:         115     0.002322 
%     Column count:          146     0.002438 
%     Minimum deg:           116     0.002337 



%           50 x 50 grid:  N=2500

%     RESULTS SUMMARY    NonZeroes   Fact.Time
%     original:           125049     1.724076 
%     Cuthill-McKee:       87025     1.038575 
%     Column count:       339131     9.634029 
%     Minimum deg:         37297     0.424364