%
% Script to test iterative solvers on systems in MATRIX-MARKET format
% Requires the routine mmread 
%
% Default parameters are:
%       - Maximum number of matvecs: 2000 
%       - Tolerance: 1e-8
%
% Programmed by Martin van Gijzen
%
%   The software is distributed without any warranty.
%
%   Copyright (c) December 2008
%   Version August 5, 2010
%

clear all; close all;
[filename, pathname] = uigetfile('*.mtx','Choose a matrix');
[dum, name] = fileparts(filename);
matrix = [pathname, filename];
rhs  = [pathname, name, '_b.mtx'];
[A, n, m, entries, rep, field, symm] = mmread(matrix);
[b, n1, m1, entries, rep, field, symm] = mmread(rhs);
if ( n1 ~= n )
   error('Size of rhs-vector incompatible with matrix  dimension'); 
end
if ( m1 ~= 1 )
   b = b(:,1);
   warning('rhs contains more vectors, taking first'); 
end
disp(['The size of the matrix is ',num2str(n)]);
spy(A);

prec = (input('ILU preconditioning (y/n)? ','s') == 'y');
if ( prec )
   droptol = input('Drop tolerance for ILU precondioner droptol = ');
   [L U] =luinc(A,droptol);
end

tol = 1e-8;
max_mv = 2000;

choice = 1;
fig = figure;
hold on;
xlabel('Number of MATVECS')
ylabel('|r|/|b|')
title(['Convergence for ',filename] );
grid on;

c = ['b','g','r','c','m','y','k'];
k = 0;
methods = [];

while ( choice < 7 )

   choice = menu('Choose a solver: ', 'GMRES(m)','Bi-CGSTAB','CGS','Bi-CG','QMR','IDRS(s)','STOP' );
   if ( choice == 7 )
      close all;
      return
   end

   k = k + 1;
   if k > 7
      k = 1;
   end

   t = cputime;
   if ( choice == 1 )
      m = input('m = ');
      max_it = ceil(max_mv/m);
      if ( prec ) 
         [x,flag,relres,iter,resvec] = gmres(A,b,m,tol,max_it, L,U );
      else
         [x,flag,relres,iter,resvec] = gmres(A,b,m,tol,max_it);
      end
      iter = m*(iter(1)-1)+iter(2);
      solver = ['GMRES(',num2str(m),')'];
   elseif ( choice == 2 )
      max_it = ceil(max_mv/2);
      if ( prec ) 
         [x,flag,relres,iter,resvec] = bicgstab(A,b,tol,max_it, L,U);
      else
         [x,flag,relres,iter,resvec] = bicgstab(A,b,tol,max_it);
      end
      solver = 'Bi-CGSTAB';
   elseif ( choice == 3 )
      max_it = ceil(max_mv/2);
      if ( prec ) 
         [x,flag,relres,iter,resvec] = cgs(A,b,tol,max_it, L,U);
      else
         [x,flag,relres,iter,resvec] = cgs(A,b,tol,max_it);
      end
      solver = 'CGS';
   elseif ( choice == 4 )
      max_it = ceil(max_mv/2);
      if ( prec ) 
         [x,flag,relres,iter,resvec] = bicg(A,b,tol,max_it, L,U);
      else
         [x,flag,relres,iter,resvec] = bicg(A,b,tol,max_it);
      end
      solver = 'Bi-CG';
   elseif ( choice == 5 )
      max_it = ceil(max_mv/2);
      if ( prec ) 
         [x,flag,relres,iter,resvec] = qmr(A,b,tol,max_it, L,U);
      else
         [x,flag,relres,iter,resvec] = qmr(A,b,tol,max_it);
      end
      solver = 'QMR';
   elseif ( choice == 6 )
      s = input('s = ');
      max_it = max_mv;
      if ( prec )
         [x,flag,relres,iter,resvec] = idrs(A,b,s,tol,max_it,L,U);
      else
         [x,flag,relres,iter,resvec] = idrs(A,b,s,tol,max_it);
      end
      solver = ['IDR(',num2str(s),')'];
   end 

   time = cputime - t;

   true_err = norm(b-A*x)/norm(b);
   disp(['Results for ', solver]);
   disp(['|b-Ax|/|b| = ',num2str(true_err)]);
   disp(['Iterations: ',num2str(iter)]);
   disp(['CPU time: ',num2str(time),'s.']);
   disp(' ');

   figure(fig);
   if ( choice == 3 || choice == 4 || choice == 5 )
      xas = 2*[0:1:length(resvec)-1];
   else
      xas = [0:1:length(resvec)-1];
   end 
   resvec = log10(resvec/norm(resvec(1)));
   plot(xas,resvec,c(k));
   if isempty(methods)
      methods = solver;
   else
      methods = char(methods,solver);
   end
   legend(methods);

end