function [itns,timest,nzerosR] = runmethods(A,b,tol,cuts,maxit,toprint)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% function [itns,timest,nzerosR] = runmethods(A,b,tol,cuts,maxit,toprint)
%
% This function runs a set of methods to solve the linear system
% Ax = b and then throws away the solution, returning only information
% on the cost of the various methods.
%
% The methods are:
%    1. Cholesky factorization
%    2. Gauss-Seidel
%    3. Conjugate gradients
%    4. Preconditioned Conjugate gradients, with preconditioner
%       symmetric Gauss-Seidel
%    5-(4+length(cuts)): Preconditioned Conjugate gradients, with 
%       preconditioner ICCG with various drop tolerances
% 
% Input:
%        A:    Array containing the matrix for the problem.
%        b:    Array containing the right-hand side for the problem.
%      tol:    The iteration stops with success when the residual
%              norm ||b - A x|| has been reduced by a factor tol.
%     cuts:    Parameters to use in the ICCG preconditioner for CG.
%    maxit:    The maximum number of iterations allowed for each method.
%  toprint:    = 1 if results should be displayed,
%              = 0 otherwise.
% Output:
%
%     itns:   Array containing the number of iterations for each method.
%   timest:   Array containing the time taken by each method.
%  nzerosR:   Array containing the number of nonzeros in any matrix
%             stored in addition to A.
%
% Dianne P. O'Leary 04/2006
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


n = length(b);
nnzA = nnz(A);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Try Cholesky factorization.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

tic
R = chol(A);
x = R\ (R' \ b);
timest(1) = toc;
nzerosR(1) = nnz(R);
itns(1) = 0;
[ir,jc,s] = find(R);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Try Gauss-Seidel, reducing iterations if problem is too big,
% since it will not converge anyway.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

tic
x = zeros(n,1);
if (n < 4000)
   [x,itns(2)] = gauss_seidel(A,b,x,tol,maxit);
else
   [x,itns(2)] = gauss_seidel(A,b,x,tol,maxit/10);
end
timest(2) = toc;
timest(2) = toc;
nzerosR(2) = 0;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Try CG with no preconditioner.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

tic
M = speye(n);
[x,flag,relres,itns(3),resvec] = pcg(A,b,tol,maxit,M,M');
timest(3) = toc;
nzerosR(3) = 0;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Try CG with Symmetric Gauss-Seidel preconditioner.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

tic
M = tril(A);
[x,flag,relres,itns(4),resvec] = pcg(A,b,tol,maxit,M,M');
timest(4) = toc;
nzerosR(4) = nnz(M);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Try CG with Incomplete Cholesky Preconditioners.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

j = 4;

for cut=cuts,
   j = j + 1;
   tic
   R = cholinc(A, cut );
   [x,flag,relres,itns(j),resvec] = pcg(A,b,tol,maxit,R',R);
   timest(j) = toc;
   nzerosR(j) = nnz(R);
end

storageest = nnzA + nzerosR;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Print a summary of results, if requested.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if (toprint=='yes')
  disp(sprintf('            n   cut   nnz(R)       #iter   time   storage_est'))
  disp(sprintf('CHOL %8d       %8d %8d %8d   %7d    %d', ...
       n,nzerosR(1),itns(1),fix(timest(1)*1000),storageest(1)))
  disp(sprintf('GS   %8d       %8d %8d %8d   %7d    %d', ...
       n,nzerosR(2),itns(2),fix(timest(2)*1000),storageest(2)))
  disp(sprintf('CG   %8d       %8d %8d %8d   %7d    %d', ...
       n,nzerosR(3),itns(3),fix(timest(3)*1000),storageest(3)))
  disp(sprintf('CGGS %8d       %8d %8d %8d   %7d    %d', ...
       n,nzerosR(4),itns(4),fix(timest(4)*1000),storageest(3)))
  for kk=1:length(cuts)
   disp(sprintf('ICCG %8d %5.2f %8d %8d  %7d      %d', ...
       n,cuts(kk),nzerosR(kk+4),itns(kk+4),fix(timest(kk+4)*1000), ...
       storageest(kk+4)))
  end
end