function fill_in_demo(A)

% Demonstrates fill-in when Cholesky and incomplete Cholesky
% decompositions are applied to a sparse matrix. 
%
% As currently set up, incomplete Cholesky decomposition is done
% using MATLAB's ichol command with options structure ic_opts and
% fields set as follows:  
%
%     ic_opts.type = 'ict'; ic_opts.michol = 'on'; ic_opts.droptol = t;
%
% for 1 = 1.e-1, 1.e-2, 1.e-3. 
%
% Input: 
%  A = sparse matrix. 
%

if ~issparse(A), 
  fprintf('\n fill_in_demo: Recasting the input matrix in sparse storage mode.\n'); 
  A = sparse(A); 
end

fig = gcf; clf; hold on; 
LA = tril(A);
spy(LA); title(['Lower Triangle of the Original Matrix']);
fprintf('\n Spy plot of the lower triangle of the original matrix.');

LAones = spones(LA);

init_storage =  length(find(tril(LAones)));
fprintf('\n Number of nonzeros in the original lower triangle = %d \n', init_storage);
fprintf('\n Press any key to continue.\n\n');
pause;

% Perform Cholesky decomposition. 

L = chol(A,'lower'); 
Lones = spones(L);
figure(fig);
spy(Lones.*LAones,'g');spy(Lones-LAones>0,'r'); 
title({'Complete Cholesky Decomposition';...
       ['Green: existing matrix entries. '...
        'Red: additional filled-in entries.']})
%pause;
new_cc_storage = length(find(Lones-LAones)); 
xlabel(['Number of filled-in entries = ',num2str(new_cc_storage)]);   
hold off; 

fprintf('\n Spy plot after complete Cholesky decomposition.');
fprintf('\n Number of new nonzeros from fill-in = %d ', new_cc_storage);
fprintf('\n   => %5.2f%% increase in storage \n', 100*new_cc_storage/init_storage);

% Perform incomplete Cholesky decomposition with increasingly small
% drop tolerances. 

ic_opts.type = 'ict'; ic_opts.michol = 'on'; 

fprintf('\n Perform incomplete Cholesky decomposition using MATLAB''s ichol with '); 
fprintf('\n options structure ic_opts and fields as follows:');
fprintf('\n         ic_opts.type = ''%s'', ic_opts.michol = ''%s'' \n', ic_opts.type, ic_opts.michol);

for i=1:5

  fprintf('\n Press any key to continue.\n\n');
  pause;

  t = 10^(-i);
  ic_opts.droptol = t;
  
  L = ichol(A,ic_opts); 
  Lones = spones(L);

  figure(fig); clf; spy(LAones,'g'); hold on; 
  title({['ICT Decomposition, t = ', num2str(t),'.'];...
         ['Green: existing matrix entries. '...
        'Red: additional filled-in entries.']})
  spy(Lones.*LAones,'g');spy(Lones-LAones>0,'r'); hold off; 
  new_ic_storage = length(find(Lones-LAones>0)); 
  xlabel(['Number of filled-in entries = ',num2str(new_ic_storage)]);       
  hold off; 
  fprintf('\n Spy plot after ICT decomposition with t = %4.2e.',t);
  fprintf('\n Number of new nonzeros from fill-in = %d ', new_ic_storage );
  fprintf('\n   => %5.2f%% increase in storage ', 100*new_ic_storage/init_storage);
  fprintf('\n   => %5.2f%% of increased storage required by complete Cholesky \n', 100*new_ic_storage/new_cc_storage);

end

