function x = lusolve(A,b,Pcol)
% LUSOLVE : Solve linear systems by supernodal LU factorization.
% 
%  x = lusolve(A, b) returns the solution to the linear system A*x = b,
%      using a supernodal LU factorization that is faster than Matlab's 
%      builtin LU.  This m-file just calls a mex routine to do the work.
%
%  By default, A is preordered by column minimum degree before factorization.
%  Optionally, the user can supply a desired column ordering:
%
%  x = lusolve(A, b, pcol) uses pcol as a column permutation.  
%      It still returns x = A\b, but it factors A(:,pcol) (if pcol is a 
%      permutation vector) or A*Pcol (if Pcol is a permutation matrix).
%       
%  x = lusolve(A, b, 0) suppresses the default minimum degree ordering;
%      that is, it forces the identity permutation on columns.
%
%  See also SUPERLU.
%
% John Gilbert, 6 April 1995
% Copyright (c) 1995 by Xerox Corporation.  All rights reserved.
% HELP COPYRIGHT for complete copyright and licensing notice.


[m,n] = size(A);
if m ~= n 
    error('matrix must be square'); 
end;
[mb,nb] = size(b);
if mb ~= n 
    error('right-hand side must have same row dimension as matrix');
end;
if n == 0
    x = [];
    return;
end;

% As necessary, compute the column permutation, and
% convert it from a permutation vector to a permutation matrix
% to fit the internal data structures of mexlusolve.
if nargin < 3 
    Pcol = colamd(A);
end;
if isempty(Pcol) | Pcol == 0
    Pcol = speye(n); 
end;
if min(size(Pcol)) == 1
    Pcol = sparse(1:n,Pcol,1,n,n);
end;

% Make sure the matrices are sparse and the vector is full.
if ~issparse(A)
    A = sparse(A);
end;
if issparse(b)
    b = full(b);
end;
if ~issparse(Pcol)
    Pcol = sparse(Pcol);
end;

x = mexlusolve(A,b,Pcol);