% P1

% Examples...
spdIndices = [1 2 3 4 35 36 39 45 47 50 67 75 219 244 341 349 355 362 ...
              413 761 791  817 888 889 936 1184  1207 ...
              1211 1214 1283 1288 1289 1313 1314 1330 ...
              1347 1406 1409 1437 1439 1406 1606 1847 1909  1912 2213 2375];
nRep = 10;

disp('                               Solve/Factorization                           ')
disp(' Index    n   nnz(A)/n^2    chol   chol+amd  chol+symrcm    chol/amd  chol/symrcm')
disp('----------------------------------------------------------------------------------')
for i=spdIndices
    Problem = UFget(i);
    % Get the matrix, which will be in sparse format..
    A = Problem.A;
    [n,n]= size(A);
    rho = nnz(A)/n^2;
    b = randn(n,1);
    xe = A\b;
    
    % Cholesky on A..
    tic
    for k=1:nRep, G = chol(A,'lower'); end
    alfa0 = toc;
    tic
    for k=1:nRep,  x0 = G'\(G\b); end
    beta0 = toc;
    q0 = beta0/alfa0;
    
    % Cholesky with approximate minimum degree pre-ordering..
    tic
    for k=1:nRep, [G1,flag,p] = chol(A,'lower','vector'); end
    alfa1 = toc;
    tic
    for k=1:nRep, x1 = b(p); x1 = G1'\(G1\x1); x1(p) = x1; end
    beta1 = toc;
    q1 = beta1/alfa1;
    
    % Cholesky with Cuthill--McKee preordering...
    tic
    for k=1:nRep, q = symrcm(A); G2 = chol(A(q,q),'lower'); end
    alfa2 = toc;
    tic
    for k=1:nRep, x2 = G2'\(G2\b(q)); x2(q) = x2; end
    beta2 = toc;
    q2 = beta2/alfa2;
    r1 = (alfa0+beta0)/(alfa1+beta1);
    r2 = (alfa0+beta0)/(alfa2+beta2);
    fprintf('%5d  %5d   %8.6f  %8.3f %9.3f %9.3f   %12.3f  %7.3f\n',i,n,rho,q0,q1,q2,r1,r2)
end

%-1 if you did not use the permutations during the solve process