% Partial Solution to Hmwk 3, Challenge 1:
% Exploring various algorithms to solve a linear
% system of equations A x = b.
% DPO  04/2010

load msc10848
S = Problem.A;
clear  Problem
b = ones(10848,1);
myplot = 0;
myrun=0;

tic, L = chol(S)'; t(1) = toc;                             
if (myplot)
   figure
   subplot(1,2,1)
   spy(S), title('A Sparse Symmetric Matrix')   
   subplot(1,2,2)
   spy(L), title('Cholesky decomposition of S')               
   print -deps msc_orig.eps 
end
nzz(1) = nnz(L);
tic, xo = S \ b; soltime(1) = toc;
err(1) = 0;

if (myrun)
%  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);                                             
if (myplot)
   figure
   subplot(1,2,1)
   spy(S(p,p)), title('S(p,p) after Cuthill-McKee ordering')  
end


%  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;                       
nzz(2) = nnz(L);
if (myplot)
   subplot(1,2,2)
   spy(L), title('chol(S(p,p)) after Cuthill-McKee ordering')
   print -deps msc_cm.eps 
end
tic, p = symrcm(S); xc(p,1) = S(p,p) \ b(p); soltime(2) = toc;
err(2) = norm(xc-xo);

%  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);                                         
if (myplot)
   figure
   subplot(1,2,1)
   spy(S(q,q)), title('S(q,q) after column count ordering')
end

%  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;                      
nzz(3)=nnz(L);
if (myplot)
   subplot(1,2,2)
   spy(L), title('chol(S(q,q)) after column count ordering')
   print -deps msc_cc.eps 
end
tic, p = colperm(S); xc(p,1) = S(p,p) \ b(p); soltime(3) = toc;
err(3) = norm(xc-xo);

%  The SYMAMD command uses the approximate minimimum degree                 
%  algorithm (a powerful graph-theoretic technique) to produce  
%  large blocks of zeros in the matrix.                         
              
r = symamd(S);                                            
if (myplot)
   figure
   subplot(1,2,1)
   spy(S(r,r)), title('S(r,r) after minimum degree ordering')
end

%  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;                        
nzz(4)=nnz(L);
if (myplot)
   subplot(1,2,2)
   spy(L), title('chol(S(r,r)) after minimum degree ordering')
   print -deps msc_md.eps 
end
tic, p = symamd(S); xc(p,1) = S(p,p) \ b(p); soltime(4) = toc;
err(4) = norm(xc-xo);


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) ))
                                                             %   
end % if myrun
                                                               
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%  Iterative methods
%
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% No precond.

b = ones(10848,1);
tic
 [x_no,flag_no,relres_no,iter_no,resvec_no] =pcg(S,b,1.e-5,1000);
soltime(5) = toc;
err(5) = norm(x_no-xo);
disp(sprintf('No Precond    flag=%d, rel.resid=%e,iter=%d',...
              flag_no, relres_no, iter_no))

% Incomplete Cholesky-infinity factorization

tic
R = cholinc(S,'inf');
 [x_ci,flag_ci,relres_ci,iter_ci,resvec_ci] =pcg(S,b,1.e-5,1000,R',R);
soltime(6) = toc;
err(6) = norm(x_ci-xo);
disp(sprintf('IncCholInf    flag=%d, rel.resid=%e,iter=%d',...
              flag_ci, relres_ci, iter_ci))

% Incomplete Cholesky-infinity factorization after amd

tic
r = symamd(S);
Sr = S(r,r);
Rr = cholinc(Sr,'inf');
 [x_ni,flag_ni,relres_ni,iter_ni,resvec_ni] =pcg(Sr,b(r),1.e-5,1000,Rr',Rr);
x_ni = x_ni(r);
soltime(7) = toc;
err(7) = norm(x_ni-xo);
disp(sprintf('IncCholInfAMD flag=%d, rel.resid=%e,iter=%d',...
              flag_ni, relres_ni, iter_ni))


disp(        'RESULTS SUMMARY    Soln.Time  Error')

disp(sprintf('original:          %7.2f     %e ',  soltime(1),err(1) ))
disp(sprintf('Cuthill-McKee:     %7.2f     %e ',  soltime(2),err(2) ))
disp(sprintf('Column count:      %7.2f     %e ',  soltime(3),err(3) ))
disp(sprintf('Minimum deg:       %7.2f     %e ',  soltime(4),err(4) ))
disp(sprintf('PCG, R=I           %7.2f     %e ',  soltime(5),err(5) ))
disp(sprintf('PCG, R=cholinc     %7.2f     %e ',  soltime(6),err(6) ))
disp(sprintf('PCG, R=cholinc,amd %7.2f     %e ',  soltime(7),err(7) ))

%----------------------------------------------------------

% Incomplete Cholesky factorization after amd

disp('Study pcg convergence as a function of droptol')
disp('   for incomplete Cholesky factorization after amd')

for kk=0:2:16,
   tic
   r = symamd(S);
   Sr = S(r,r);
   options.droptol = (10^kk)*eps;
   options.rdiag = 1;
   Rr = cholinc(Sr,options);
    [x_ni,flag_ni,relres_ni,iter_ni,resvec_ni] =pcg(Sr,b(r),1.e-5,1000,Rr',Rr);
   x_ni = x_ni(r);
   soltime(7) = toc;
   err(7) = norm(x_ni-xo);
   disp(sprintf('droptol = %e, nnz(Cholesky factor) = %d',...
        options.droptol,nnz(Rr) ))
   disp(sprintf('              flag=%d, rel.resid=%e,iter=%d',...
                 flag_ni, relres_ni, iter_ni))
   
   disp(sprintf('solution time =  %7.2f, error = %e ',  soltime(7),err(7) ))
end
