kusano 7d535a
function info = trytime(A);
kusano 7d535a
% TRYTIME : time SUPERLU with 3 outputs against Matlab's LU
kusano 7d535a
%
kusano 7d535a
% info = trytime(A);
kusano 7d535a
% normally info is ||P*A*x-L*U*x|| / ||A|| (for a random x);
kusano 7d535a
% but info is at least 10^6 if the factors are not triangular,
kusano 7d535a
% or the permutation isn't a permutation.
kusano 7d535a
% Copyright (c) 1995 by Xerox Corporation.  All rights reserved.
kusano 7d535a
% HELP COPYRIGHT for complete copyright and licensing notice.
kusano 7d535a
kusano 7d535a
kusano 7d535a
n = max(size(A));
kusano 7d535a
info = 0;
kusano 7d535a
kusano 7d535a
format compact
kusano 7d535a
disp('SUPERLU with 3 outputs:');
kusano 7d535a
kusano 7d535a
tic;[L,U,P]=lu(A);t1=toc;
kusano 7d535a
fprintf('LU time = %d seconds.\n',t1);
kusano 7d535a
tic;[l,u,pr]=superlu(A);t2=toc;
kusano 7d535a
fprintf('SUPERLU time = %d seconds.\n',t2);
kusano 7d535a
ratio = t1/t2;
kusano 7d535a
fprintf('Ratio = %d\n',ratio);
kusano 7d535a
kusano 7d535a
if any(any(triu(l,1)))
kusano 7d535a
    disp('L is *NOT* lower triangular.');
kusano 7d535a
    info = info + 10^6;
kusano 7d535a
else
kusano 7d535a
    disp('L is lower triangular.');
kusano 7d535a
end;
kusano 7d535a
if nnz(l) == nnz(l+l)
kusano 7d535a
    disp('L has no explicit zeros.');
kusano 7d535a
else
kusano 7d535a
    disp('L contains explicit zeros.');
kusano 7d535a
    info = info+10^6;
kusano 7d535a
end;
kusano 7d535a
if any(any(tril(u,-1)))
kusano 7d535a
    disp('U is *NOT* upper triangular.');
kusano 7d535a
    info = info + 10^6;
kusano 7d535a
else
kusano 7d535a
    disp('U is upper triangular.');
kusano 7d535a
end;
kusano 7d535a
if nnz(u) == nnz(u+u)
kusano 7d535a
    disp('U has no explicit zeros.');
kusano 7d535a
else
kusano 7d535a
    disp('U contains explicit zeros.');
kusano 7d535a
    info = info+10^6;
kusano 7d535a
end;
kusano 7d535a
if pr == [1:n]
kusano 7d535a
    disp('PROW is the identity permutation.');
kusano 7d535a
elseif isperm(pr)
kusano 7d535a
    disp('PROW is a non-identity permutation.');
kusano 7d535a
else 
kusano 7d535a
    disp('PROW is *NOT* a permutation.');
kusano 7d535a
    info = info + 10^6;
kusano 7d535a
end;
kusano 7d535a
kusano 7d535a
x = rand(n,1);
kusano 7d535a
rnorm = norm(A(pr,:)*x - l*(u*x),inf)/norm(A,inf);
kusano 7d535a
fprintf(1,'||A(PROW,:)*x -  L*U*x||/||A|| = %d\n', rnorm);
kusano 7d535a
kusano 7d535a
info = info + rnorm;
kusano 7d535a
disp(' ');