function  [nz,ctime,rnorm] = runmethods(A,b)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% function  [nz,ctime,rnorm] = runmethods(A,b)
% This function runs several methods to solve Ax=b.
%
% Outputs:  The vectors of outputs return performance.
% For method j: 
%
%  nz(j)      storage used by the method
%  ctime(j)   time taken by the method
%  rnorm(j)   norm of relative residual returned by the method
%                  ||b-Ax|| / ||b||
%
% Note that, in academic style, we throw away the computed solution x
% and just return performance data for the methods.
%
% The methods:
%  j=1  Cholesky on the original matrix
%  j=2  Cholesky with reverse Cuthill-McKee reordering (symrcm)
%  j=3  Cholesky with minimum degree reordering (symmmd)
%  j=4  Cholesky with approximate minimum degree reordering (symamd)
%  j=5  Cholesky with eigenpartioning reordering (specnd})
%
% Dianne O'Leary 06-2005
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

normb = norm(b);
nz    = zeros(1,5);
ctime = zeros(1,5);
rnorm = zeros(1,5);

n = length(b);
nz(1) = 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;
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;
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;
rnorm(3) = norm(b-A*x3)/normb;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Solving using the Cholesky factorization
%       with approximate minimum degree ordering
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

tic
p = symamd(A);
R = chol(A(p,p));
x4 = R \ (R' \ b(p));
x4(p) = x4;
nz(4) = 3*nnz(R) + 3*nnz(A) + 2*n;
ctime(4) = toc;
rnorm(4) = norm(b-A*x4)/normb;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Solving using the Cholesky factorization
%       with eigenvector partitioning
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if (n < 4000)  % having trouble with out-of-memory
   tic
   p = specnd(A);
   R = chol(A(p,p));
   x5 = R \ (R' \ b(p));
   x5(p) = x5;
   nz(5) = 3*nnz(R) + 3*nnz(A) + 2*n;
   ctime(5) = toc;
   rnorm(5) = norm(b-A*x5)/normb;
end


