function failed = trylusolve(A,b);
% TRYLUSOLVE : Test the Matlab interface to lusolve.
%
% failed = trylusolve;  
%          This runs several tests on lusolve, using a matrix with the
%          structure of "smallmesh".  It returns a list of failed tests.
%
% failed = trylusolve(A); 
% failed = trylusolve(A,b); 
%          This times the solution of a system with the given matrix,
%          both with lusolve and with Matlab's builtin "\" operator.
%          Usually, "\" will be faster for triangular or symmetric postive
%          definite matrices, and lusolve will be faster otherwise.
%
% Copyright (c) 1995 by Xerox Corporation.  All rights reserved.
% HELP COPYRIGHT for complete copyright and licensing notice.

disp(' ');
disp(['Testing LUSOLVE on ' time]);
disp(' ');

format compact
resetrandoms;
if nargin == 0
    A = hbo('smallmesh');
end;
[n,m] = size(A);
if n~=m, error('matrix must be square'), end;

failed = [];
ntest = 0;
tol = 100 * eps * n^2;

if nargin >= 1 % Just try a solve with the given input matrix.
    ntest=ntest+1;
    v = spparms;
    spparms('spumoni',1);
    fprintf('Matrix size: %d by %d with %d nonzeros.\n',n,m,nnz(A));
    disp(' ');
    if nargin == 1
        b = rand(n,1);
    end;
    tic; x=A\b; t1=toc;
    fprintf('A\\b time = %d seconds.\n',t1);
    disp(' ');
    tic; x=lusolve(A,b); t2=toc;
    fprintf('LUSOLVE time = %d seconds.\n',t2);
    disp(' ');
    ratio = t1/t2;
    fprintf('Ratio = %d\n',ratio);
    residual_norm = norm(A*x-b,inf)/norm(A,inf)
    disp(' ');
    if residual_norm > tol
         disp('*** FAILED ***'); failed = [failed ntest];
    end;
    spparms(v);
    return;
end;

% With no inputs, try several tests on the small mesh.


A = sprandn(A);
p = randperm(n);
I = speye(n);
P = I(p,:);
b = rand(n,2);

ntest=ntest+1;
fprintf('Test %d: Input perm 0.\n',ntest);
x = lusolve(A,b,0);
residual_norm = norm(A*x-b,inf)
disp(' ');
if residual_norm > tol 
    disp('*** FAILED ***'), disp(' '), failed = [failed ntest]; 
end;

ntest=ntest+1;
fprintf('Test %d: Input perm given as vector.\n',ntest);
x = lusolve(A,b,p);
residual_norm = norm(A*x-b,inf)
disp(' ');
if residual_norm > tol
    disp('*** FAILED ***^G'), disp(' '), failed = [failed ntest];
end;

ntest=ntest+1;
fprintf('Test %d: Input perm given as matrix.\n',ntest);
x = lusolve(A,b,P);
residual_norm = norm(A*x-b,inf)
disp(' ');
if residual_norm > tol 
    disp('*** FAILED ***'), disp(' '), failed = [failed ntest]; 
end;

ntest=ntest+1;
fprintf('Test %d: No input perm (colamd done internally).\n',ntest);
x = lusolve(A,b);
residual_norm = norm(A*x-b,inf)
disp(' ');
if residual_norm > tol 
    disp('*** FAILED ***'), disp(' '), failed = [failed ntest]; 
end;

ntest=ntest+1;
fprintf('Test %d: Empty matrix.\n',ntest);
x = lusolve([],[]);
disp(' ');
% if max(size(x))
if length(x)
    x
    disp('*** FAILED ***^G'), disp(' ');
    failed = [failed ntest];
end;

ntest=ntest+1;
fprintf('Test %d: Timing versus \\ on the 3-element airfoil matrix.\n',ntest);
A = hbo('airfoil2');
n = max(size(A));
A = sprandn(A);
b = rand(n,1);
tic; x=A\b; t1=toc;
fprintf('A\\b time = %d seconds.\n',t1);
tic; x=lusolve(A,b); t2=toc;
fprintf('LUSOLVE time = %d seconds.\n',t2);
ratio = t1/t2;
fprintf('Ratio = %d\n',ratio);
disp(' ');
if ratio < 1, disp('*** FAILED ***'), disp(' '), failed = [failed ntest]; end;


nfailed = length(failed);
if nfailed
    fprintf('\n%d tests failed.\n\n',nfailed);
else
    fprintf('\nAll tests passed.\n\n',nfailed);
end;