function  [nz,cflag,ctime,iter,rnorm] = runmethods(A,b,tol,droptol)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% function  [nz,cflag,ctime,iter,rnorm] = runmethods(A,b,tol,droptol)
% This function runs nine algorithms to solve Ax=b.
% For the iterative methods, the number tol determines
%   the convergence tolerance for the relative residual.
% The parameter droptol is the drop tolerance for the incomplete
%   Cholesky preconditioner.
%
% Outputs:  The vectors of outputs return performance.
% For each algorithm. 
%
%  nz      storage used by the method
%  cflag   error flag returned from the method
%  ctime   time taken by the method
%  iter    number of iterations taken by the method
%  rnorm   norm of relative residual returned by the method
%
% Note that, in academic style, we throw away the computed solution x
% and just return performance data for the algorithms.
%
% Dianne O'Leary 04-2005
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

cflag = zeros(9,1);
ctime = zeros(9,1);
normb = norm(b);

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

% Solving using the Cholesky factorization

tic
R = chol(A);
x1 = R \ (R' \ b);
nz(1) = 3*nnz(R) + 3*nnz(A) + 2*n;
ctime(1) = toc;
iter(1) = 1;
rnorm(1) = norm(b-A*x1)/normb;

% Solving using the Cholesky factorization
%       with Reverse-Cuthill-McKee ordering

tic
p = symrcm(A);
R = chol(A(p,p));
x2 = R \ (R' \ b(p));
x2(p) = x2;
nz(2) = 3*nnz(R) + 3*nnz(A) + 2*n;
ctime(2) = toc;
iter(2) = 1;
rnorm(2) = norm(b-A*x2)/normb;

% Solving using the Cholesky factorization
%       with minimum degree ordering
% (Note that Matlab has declared symmmd obsolete,
%  because it is too expensive, and is replacing
%  it by symamd, which produces an approximate
%  minimum degree ordering at less cost.)

tic
p = symmmd(A);
R = chol(A(p,p));
x3 = R \ (R' \ b(p));
x3(p) = x3;
nz(3) = 3*nnz(R) + 3*nnz(A) + 2*n;
ctime(3) = toc;
iter(3) = 1;
rnorm(3) = norm(b-A*x3)/normb;

% Solving using Gmres, restart =  20

restart = 20;
tic
[x4,cflag(4),relres,itr] = gmres(A,b,restart,tol);
nz(4) = 3*nnz(A) + 2*n + restart*n + restart^2;
ctime(4) = toc;
iter(4) = (itr(1)-1)*restart + itr(2);
rnorm(4) = norm(b-A*x4)/normb;

% Solving using Gmres, restart = 100

restart = 100;
tic
[x5,cflag(5),relres,itr] = gmres(A,b,restart,tol);
nz(5) = 3*nnz(A) + 2*n + restart*n + restart^2;
ctime(5) = toc;
iter(5) = (itr(1)-1)*restart + itr(2);
rnorm(5) = norm(b-A*x5)/normb;

%  CG, no preconditioning
tic
[x6,cflag(6),relres,iter(6)] = pcg(A,b,tol,500);
nz(6) = 3*nnz(A) + 5*n;
ctime(6) = toc;
rnorm(6) = norm(b-A*x6)/normb;

%  CG, diagonal preconditioning
tic
M = spdiags(diag(A),0,n,n);
[x7,cflag(7),relres,iter(7)] = pcg(A,b,tol,500,M);
nz(7) = 3*nnz(A) + 5*n + n;
ctime(7) = toc;
rnorm(7) = norm(b-A*x7)/normb;

%  CG, incomplete Cholesky preconditioning
tic
R = cholinc(A,droptol);
[x8,cflag(8),relres,iter(8)] = pcg(A,b,tol,500,R',R);
nz(8) = 3*nnz(A) + 5*n + 3*nnz(R);
ctime(8) = toc;
rnorm(8) = norm(b-A*x8)/normb;

%  CG, incomplete Cholesky preconditioning with minimum degree ordering
tic
R = cholinc(A(p,p),droptol);
[x9,cflag(9),relres,iter(9)] = pcg(A(p,p),b(p),tol,500,R',R);
x9(p) = x9;
nz(9) = 3*nnz(A) + 5*n + 3*nnz(R);
ctime(9) = toc;
rnorm(9) = norm(b-A*x9)/normb;