%
%  This script will demonstrate some of the sparse matrix 
%  capabilities of Matlab
%
%
%  D.C. Sorensen
%  Last modified 27 Aug 04
%
   NZhistory = [];
%
%  read a matrix from matrix market
%

   echo on;
   A = mmread('blckhole.mtx');

   echo off;

   [n,m] = size(A);
   disp('order of A ... n ='), disp(n)
   disp('  ')
   disp(' Strike a key ')
   disp('  ')
   pause
   disp(' Matlab Data Structure ')
   disp(' for Sparse Matrices')

   disp( ' i , j ,  Aij ')
   
   A(1:5,1:3)
   
   disp('  ')
   disp(' Strike a key ')
   disp('  ')
   pause



   norm(A-A',1)
   disp(' norm(A-A'',1) ')
   disp('  % Use one or inf norm for sparse A ')

   disp('  ')
   disp(' Strike a key ')
   disp('  ')
   pause

   A = A + 1000*speye(n);
   disp(' A = A + 1000*speye(n);')
   disp(' % to make A pos def' ) 
   disp(' % speye = sparse I ')
   disp('  ')
   disp(' Strike a key ')
   disp('  ')
   pause
%
%  -----------------------natural ordering-------------------------------------
%

   spy(A)
   title('Natural Ordering','Fontsize',15)
   disp(' spy(A) ')
   disp('  ')
   disp(' Strike a key ')
   disp('  ')
   pause

   L = chol(A);
   disp(' L = chol(A); ')

   hold
   spy(L,'r')
   title('Cholesky Factor Fill in Red ','Fontsize',15)
   disp(' spy(L) ')
   hold

   NonzerosL = nnz(L)
   NZhistory = [NZhistory , NonzerosL];
   disp('  ')
   disp(' Strike a key ')
   disp('  ')
   pause



%  more on
%  disp(' help sparfun ')
%
%  help sparfun
%
%  more off

%
%  -----------------------symrmc ordering-------------------------------------
%
   p = symrcm(A);
   disp(' p = symrcm(A); ')
   spy(A(p,p))
   title('RCM Ordering','Fontsize',15)
   disp(' spy(A(p,p)) ')
   disp('  ')
   disp(' Strike a key ')
   disp('  ')
   pause

   L = chol(A(p,p));
   disp(' L = chol(A(p,p)); ')
   hold
   title('Cholesky Factor Fill in Red ','Fontsize',15)
   spy(L,'r')
   disp(' spy(L,"r") ')
   hold

   NonzerosL = nnz(L)
   NZhistory = [NZhistory , NonzerosL];
   disp('  ')
   disp(' Strike a key ')
   disp('  ')
   pause
%
%  -----------------------symmmd ordering-------------------------------------
%
   
   
   p = symmmd(A);
   disp(' p = symmmd(A); ')

   spy(A(p,p))
   title('Symmetric Min Degree ','Fontsize',15)
   disp(' spy(A(p,p)) ')
   disp('  ')
   disp(' Strike a key ')
   disp('  ')
   pause

   L = chol(A(p,p));
   disp(' L = chol(A(p,p)); ')
   hold
   spy(L,'r')
   disp(' spy(L,"r") ')
   title('Cholesky Factor Fill in Red ','Fontsize',15)
   hold

   NonzerosL = nnz(L)
   NZhistory = [NZhistory , NonzerosL];
   disp('  ')
   disp(' Strike a key ')
   disp('  ')
   pause




%
%  -----------------------symamd ordering-------------------------------------
%
%  more on
%  help symamd
%  more off

   p = symamd(A);
   disp(' p = symamd(A); ')

   spy(A(p,p))
   title('Approximate Min Degree ','Fontsize',15)
   disp(' spy(A(p,p)) ')
   disp('  ')
   disp(' Strike a key ')
   disp('  ')
   pause

   L = chol(A(p,p));
   disp(' L = chol(A(p,p)); ')
   hold
   spy(L,'r')
   title('Cholesky Factor Fill in Red ','Fontsize',15)
   disp(' spy(L,"r") ')
   hold

   NonzerosL = nnz(L)
   NZhistory = [NZhistory , NonzerosL];
   disp('  ')
   disp(' Strike a key ')
   disp('  ')
   pause



%
%  Display a comparision of fill in for various methods
%
   bar(NZhistory,.3)
   xlabel('  natural       symrcm      symmmd      symamd   ','Fontsize',20)
   ylabel(' Nonzeros in   L','fontsize',20)