function varargout = lus(A, d)
% LUS	Factoring matrix into elementary reversible matrices.
%	A SERM or a LUS can be implemented by integer mapping.
%
%	[P, L, U, S] = LUS(A, d) or [P, LU, S] = LUS(A, d) 
%	factorizes matrix A into a product of a permutation 
%	matrix P, a unit lower triangular matrix L and an upper
%	triangular matrix U with diagonal d (or L and U are 
%	stored compactly into one matrix LU) and a unit 
%	single-row matrix S, so that
%		P*A = L*U*S.
%
%	MaxErr = LUS(A) returns the absolutely maximum magnitude 
%	of upper bound of transform error in order to find an
%	optimal SERM factorization.
%
%	Transpose of vector formed permutation matrix:
%		for i=1:n, Pt(P(i))=i; end;
%
%	See also the biult-in function LU.
%
%	For any commercial applications of this program, a 
%	written permission MUST be obtained from the authors.
%	email:	phao@cis.pku.edu.cn
%		phao@dcs.qmul.ac.uk
%
%	If this program is used for research, please cite:
%
%	* Pengwei Hao, "Customizable Triangular Factorizations 
%	  of Matrices", Linear Algebra and Its Applications, 
%	  Vol.382, pp.135-154, May 2004.
%
%	* Pengwei Hao and Qingyun Shi, "Matrix Factorization 
%	  for Reversible Integer Mapping", IEEE Transactions 
%	  on Signal Processing, Vol.49, No.10, pp.2314-2324, 
%	  Oct. 2001.
%
%	By Pengwei Hao and Lei Yang, 2000-2005

error( nargchk( 1, 2, nargin ) );

[m,n] = size(A);
if (m~=n) error('Square matrix is acceptable!'); end;
if (abs(det(A))<eps) error('A nonsingular matrix is acceptable!'); end;

if (nargin==1)
    d = ones(1,m);
%     more examples:
%     d = [1 1 1 1 1 1 1 1];
%     d = [-1 -1 -1 1 1 -1 1 1];
end;

%disp( sprintf( 'The Determinant is %d', det(A) ) );
A0 = A;

% Initialization
P = 1:n;
S = [];
% Pr = [0, 5, 4, 1, 3, 6, 2, 7]+1;
Pr = 1:n;
disp(Pr);

A = A(:,Pr);
A1 = A;
for i = 1:(n-1),

% pivoting for largest denominator
    k = i;
    c = abs(A(i,n));
    for i1 = (i+1):n,
        c1 = abs(A(i1,n));
        if (c1>c) c = c1; k = i1; end;
    end;
% pivoting for least single-row elements
    k = i;
    c = realmax;
    for i1 = i:n,
        if (abs(A(i1,i)-d(i))<eps) c1 = 0; % the error when abs(A(i1,n))<eps????
        elseif (abs(A(i1,n))<eps) c1 = realmax;
        else c1 = abs( (A(i1,i)-d(i))/A(i1,n) );
        end;
%        if (abs(c1-0.5)<abs(c-0.5)) c = c1; k = i1; end;
        if (c1<c) c = c1; k = i1; end;
    end;

% permutaion
    tp = P(i);    P(i) = P(k);    P(k) = tp;
    for j = 1:n,
        ta = A(i,j);    A(i,j) = A(k,j);    A(k,j) = ta;
    end;

% Normalization
    if (abs(A(i,i)-d(i))<eps) s = 0;
    else s = (A(i,i)-d(i))/A(i,n);
    end;
    S(i) = s;
    for i1 = 1:n,
        A(i1,i) = A(i1,i) - s*A(i1,n);
    end;

% Elimination
    for i1 = (i+1):n,
        for j = (i+1):n,
            A(i1,j) = A(i1,j) - A(i1,i)*A(i,j)/A(i,i);%need to be divided by A(i,i)
        end;
    end;
    
    for i1 = i+1:n,% just elimiate for L ,need to be divided by A(i,i)
        A(i1,i) = A(i1,i) /A(i,i);
    end;

end;% end of the outmost for i

S = eye(n)+[zeros(n-1,1);1]*[S,0];
disp( 'RPLUS Factorization :' ); disp( A ); disp( S(n,:) );
L = tril(A) - diag(diag(A)) + eye(n);%retain the L when Factorization,now only need extraction
U = triu(A);
disp( L );
disp( U );
disp( P );
disp( 'L*U :' );disp( L*U );

% Determinant Check
A1 = A1(P,:)*inv(S);
disp( 'A1' );disp(A1);
for i = 1:n,
  D(i) = det(A1(1:i,1:i));
end;
disp( 'The Minors of the leading principal submatrices are :' ); disp( sprintf(' %d', D) );  disp(' ');


% Return out-arguments
if (nargout >= 4)
   varargout = {P, L, U, S,D,A};
elseif (nargout == 3)
   varargout = {P, A, S};
elseif (nargout == 2)
   varargout = {P, S};
end;

