%  Code: compLinSlvrs
%
%  This code compares the nonsymmetric iterative solvers from 
%  Matlab:  gmres, gmres5, bicgstab, bicg, qmr, cgs
%
%  It is set up to solve the Sherman 2 problem from Harwell Bowing
%  collection (available from Matrix Market).  The banded 
%  matrix consiting of bands [-10:10] is used as a pre-conditioner.  
%
%  Assumes files sherman2.mtx  sherman2_rhs1.mtx 
%
%          from Matrix Market exist
%
%-------------------
%
%  D.C. Sorensen
%  30 Mar 01
%
%  Revised 27 Oct 04  DCS
%   Replaced flopcounts with timings - flops no longer available
%   Relabled and re-ordered plots
%

 A = mmread('sherman2.mtx'); 
 b = mmread('sherman2_rhs1.mtx');

%
%   Construct a preconditioner from the diagonals [-10:10] 
%
 d = [-10:10];
 B = spdiags(A,d); 
 [n,n] = size(A);
 M = spdiags(B,d,n,n);  
 spy(A)    
 title('Sparsity Pattern of A','fontsize',12)
 pause
 hold
 spy(M,'r')    
 title('Block Jacobi PreConditioner','fontsize',12)
 hold

 pause

 [L,U] = lu(M); 


 spy(L,'r')
 hold
 pause

 spy(U)
 title('LU Decomposition of PreConditioner','fontsize',12)
 hold

 pause
 flpcnt = [];
%
% Uncomment for unpreconditioned methods
%
   L = speye(n);
   U = speye(n);
%
 tic
 [x,flag,relres,iter,resvec] = gmres(A,b,5,1e-8,100,L,U); 
 flops = toc;
 flpcnt = [flpcnt, flops];
 resvec = resvec/norm(b);
 semilogy(resvec,'m','linewidth',1.5) 
 title('Gmres_5','fontsize',12)
 hold
 pause

 
%flops(0);
 tic
 [x,flag,relres,iter,resvec] = bicgstab(A,b,1e-8,100,L,U);
 flops = toc;
 flpcnt = [flpcnt, flops];
 resvec = resvec/norm(b);
 semilogy(resvec,'linewidth',1.5) 
 title('BiCGstab','fontsize',12)
 pause

 [x,flag,relres,iter,resvec] = bicg(A,b,1e-8,100,L,U);
 flops = toc;
 flpcnt = [flpcnt, flops];
 resvec = resvec/norm(b);
 semilogy(resvec,'g','linewidth',1.5) 
 title('BiCG','fontsize',12)
 pause

%flops(0);
 tic
 [x,flag,relres,iter,resvec] = cgs(A,b,1e-8,100,L,U);
 flops = toc;
 flpcnt = [flpcnt, flops];
 resvec = resvec/norm(b);
 semilogy(resvec,'r','linewidth',1.5) 
 title('CGS','fontsize',12)
 pause

%flops(0);
 tic
 [x,flag,relres,iter,resvec] = gmres(A,b,n,1e-8,100,L,U); 
 flops = toc;
 flpcnt = [flpcnt, flops];
 resvec = resvec/norm(b);
 semilogy(resvec,'k','linewidth',1.5) 
 title('Gmres','fontsize',12)
 pause

%flops(0);
 tic
 [x,flag,relres,iter,resvec] = qmr(A,b,1e-8,100,L,U);
 flops = toc;
 flpcnt = [flpcnt, flops];
 resvec = resvec/norm(b);
 semilogy(resvec, 'k:','linewidth',1.5) 
 title('QMR','fontsize',12)

 pause
 title('Linear System Residual for Various methods','fontsize',12)
 xlabel('Iteration Number k','fontsize',12)
 ylabel('Scaled Residual  || b - Ax_k ||/||b||  at step k','fontsize',12)

%flops(0);
 tic
 hold
legend('gmres5','bicgstab','bicg','cgs','gmres','qmr')
  
  method = ['    gmres5    ',   'bicgstab  ','  bicg   ','   cgs    ','  gmres  ','   qmr    '];
  disp(' ')
  disp('Timings (in secs) for the various methods')
  disp(' ')
  disp(method)
  disp(flpcnt)
  disp(' ')
  disp('Warning: Timings can be unreliable in Matlab')
%
% sort(flpcnt')