kusano 7d535a
function [L,U,prow,pcol] = superlu(A,psparse)
kusano 7d535a
% SUPERLU : Supernodal LU factorization
kusano 7d535a
% 
kusano 7d535a
%  Executive summary:
kusano 7d535a
%
kusano 7d535a
%  [L,U,p] = superlu(A)          is like [L,U,P] = lu(A), but faster.
kusano 7d535a
%  [L,U,prow,pcol] = superlu(A)  preorders the columns of A by min degree,
kusano 7d535a
%                                    yielding A(prow,pcol) = L*U.
kusano 7d535a
%
kusano 7d535a
%  Details and options:
kusano 7d535a
%
kusano 7d535a
%  With one input and two or three outputs, SUPERLU has the same effect as LU,
kusano 7d535a
%  except that the pivoting permutation is returned as a vector, not a matrix:
kusano 7d535a
%
kusano 7d535a
%  [L,U,p] = superlu(A) returns unit lower triangular L, upper triangular U,
kusano 7d535a
%            and permutation vector p with A(p,:) = L*U.
kusano 7d535a
%  [L,U] = superlu(A) returns permuted triangular L and upper triangular U
kusano 7d535a
%            with A = L*U.
kusano 7d535a
%
kusano 7d535a
%  With a second input, the columns of A are permuted before factoring:
kusano 7d535a
%
kusano 7d535a
%  [L,U,prow] = superlu(A,psparse) returns triangular L and U and permutation 
kusano 7d535a
%            prow with A(prow,psparse) = L*U.
kusano 7d535a
%  [L,U] = superlu(A,psparse) returns permuted triangular L and triangular U 
kusano 7d535a
%            with A(:,psparse) = L*U.
kusano 7d535a
%  Here psparse will normally be a user-supplied permutation matrix or vector
kusano 7d535a
%  to be applied to the columns of A for sparsity.  COLAMD is one way to get
kusano 7d535a
%  such a permutation; see below to make SUPERLU compute it automatically.
kusano 7d535a
%  (If psparse is a permutation matrix, the matrix factored is A*psparse'.)
kusano 7d535a
%
kusano 7d535a
%  With a fourth output, a column permutation is computed and applied:
kusano 7d535a
%
kusano 7d535a
%  [L,U,prow,pcol] = superlu(A,psparse)  returns triangular L and U and
kusano 7d535a
%            permutations prow and pcol with A(prow,pcol) = L*U.
kusano 7d535a
%            Here psparse is a user-supplied column permutation for sparsity,
kusano 7d535a
%            and the matrix factored is A(:,psparse) (or A*psparse' if the
kusano 7d535a
%            input is a permutation matrix).  Output pcol is a permutation
kusano 7d535a
%            that first performs psparse, then postorders the etree of the 
kusano 7d535a
%            column intersection graph of A.  The postorder does not affect 
kusano 7d535a
%            sparsity, but makes supernodes in L consecutive.
kusano 7d535a
%  [L,U,prow,pcol] = superlu(A,0) is the same as ... = superlu(A,I); it does
kusano 7d535a
%            not permute for sparsity but it does postorder the etree.
kusano 7d535a
%  [L,U,prow,pcol] = superlu(A) is the same as ... = superlu(A,colamd(A));
kusano 7d535a
%            it uses column minimum degree to permute columns for sparsity,
kusano 7d535a
%            then postorders the etree and factors.
kusano 7d535a
%
kusano 7d535a
%  This m-file calls the mex-file MEXSUPERLU to do the work.
kusano 7d535a
%
kusano 7d535a
% See also COLAMD, LUSOLVE.
kusano 7d535a
%
kusano 7d535a
% John Gilbert, 6 April 1995.
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
% Note on permutation matrices and vectors:
kusano 7d535a
%
kusano 7d535a
% Names beginning with p are permutation vectors; 
kusano 7d535a
% names beginning with P are the corresponding permutation matrices.
kusano 7d535a
%
kusano 7d535a
% Thus  A(pfoo,pbar) is the same as Pfoo*A*Pbar'.
kusano 7d535a
%
kusano 7d535a
% We don't actually form any permutation matrices except Psparse,
kusano 7d535a
% but matrix notation is easier to untangle in the comments.
kusano 7d535a
kusano 7d535a
kusano 7d535a
[m,n] = size(A);
kusano 7d535a
if m ~= n 
kusano 7d535a
    error('matrix must be square.'); 
kusano 7d535a
end;
kusano 7d535a
if n == 0
kusano 7d535a
    L = []; U = []; prow = []; pcol = [];
kusano 7d535a
    return;
kusano 7d535a
end;
kusano 7d535a
kusano 7d535a
% If necessary, compute the column sparsity permutation.
kusano 7d535a
if nargin < 2 
kusano 7d535a
    if nargout >= 4
kusano 7d535a
        psparse = colamd(A);
kusano 7d535a
    else
kusano 7d535a
        psparse = 1:n;
kusano 7d535a
    end;
kusano 7d535a
end;
kusano 7d535a
if max(size(psparse)) <= 1
kusano 7d535a
    psparse = 1:n;
kusano 7d535a
end;
kusano 7d535a
kusano 7d535a
% Compute the permutation-matrix version of psparse,
kusano 7d535a
% which fits the internal data structures of mexsuperlu.
kusano 7d535a
if min(size(psparse)) == 1
kusano 7d535a
    Psparse = sparse(1:n,psparse,1,n,n);
kusano 7d535a
else
kusano 7d535a
    Psparse = psparse;
kusano 7d535a
    psparse = Psparse*[1:n]';
kusano 7d535a
end;
kusano 7d535a
kusano 7d535a
% Make sure the matrices are sparse.
kusano 7d535a
if ~issparse(A)
kusano 7d535a
    A = sparse(A);
kusano 7d535a
end;
kusano 7d535a
if ~issparse(Psparse)
kusano 7d535a
    Psparse = sparse(Psparse);
kusano 7d535a
end;
kusano 7d535a
kusano 7d535a
kusano 7d535a
% The output permutations from the mex-file are dense permutation vectors.
kusano 7d535a
[L,U,prowInv,pcolInv] = mexsuperlu(A,Psparse);
kusano 7d535a
prow = zeros(1,n);
kusano 7d535a
prow(prowInv) = 1:n;
kusano 7d535a
pcol = zeros(1,n);
kusano 7d535a
pcol(pcolInv) = 1:n;
kusano 7d535a
kusano 7d535a
% We now have
kusano 7d535a
%
kusano 7d535a
%    Prow*A*Psparse'*Post' = L*U   (1)
kusano 7d535a
%    Pcol' = Psparse'*Post'         
kusano 7d535a
%
kusano 7d535a
% (though we actually have the vectors, not the matrices, and
kusano 7d535a
% we haven't computed Post explicitly from Pcol and Psparse).
kusano 7d535a
% Now we figure out what the user really wanted returned,
kusano 7d535a
% and rewrite (1) accordingly.
kusano 7d535a
kusano 7d535a
if nargout == 4  
kusano 7d535a
    % Return both row and column permutations.  
kusano 7d535a
    % This is what we've already got.
kusano 7d535a
    % (1) becomes  Prow*A*Pcol' = L*U. 
kusano 7d535a
elseif nargout == 3
kusano 7d535a
    % Return row permutation only.  Fold the postorder perm 
kusano 7d535a
    % but not the sparsity perm into it, and into L and U.
kusano 7d535a
    % This preserves triangularity of L and U.
kusano 7d535a
    % (1) becomes (Post'*Prow) * A * Psparse' = (Post'*L*Post) * (Post'*U*Post).
kusano 7d535a
    postInv = pcolInv(psparse);
kusano 7d535a
    prow = prow(postInv);
kusano 7d535a
    L = L(postInv,postInv);
kusano 7d535a
    U = U(postInv,postInv);
kusano 7d535a
else % nargout <= 2
kusano 7d535a
    % Return no permutations.  Fold the postorder perm but
kusano 7d535a
    % not the sparsity perm into L and U.  Fold the pivoting
kusano 7d535a
    % perm into L, which destroys its triangularity.
kusano 7d535a
    % (1) becomes A * Psparse' = (Prow'*L*Post) * (Post'*U*Post).
kusano 7d535a
    postInv = pcolInv(psparse);
kusano 7d535a
    L = L(prowInv,postInv);
kusano 7d535a
    U = U(postInv,postInv);
kusano 7d535a
end;