function [L,U,P,Q,rank1]=lupan3(A,nb,rcond1,mu)

%
%   lupan3 finds a rank-revealing LU factorization of A using
%   (a modification of) Pan's algorithm 3
% 
%   [L,U,P,Q] = lupan3(A); returns unit lower triangular matrix L, upper
%   triangular matrix U,  permutation matrix P and permutation matrix
%   Q so that P*A*Q = L*U.
%
%   [L,U,P,Q]=lupan3(A,nb); uses a block size nb in the algorithm. 
%   If there is only one parameter or if nb is zero use the default
%   size nb = 32.
%
%   For the call [L,U,P,Q,k]=lupan3(A,nb,rcond1) then
%   the  columns and rows of A are permutated so that B11, the k by k
%   upper left corner of P*A*Q, is well conditioned in the following
%   sense:  If we partition 
%
%        L = [ L11    0  ]    and  U = [  U11    U12 ]
%            [ L21   L22 ]             [   0     U22 ]
%
%   then
%          (1) cond(B11) is smaller, approximately, than 1 / rcond1
%          (2) if k is increased by one this property is violated
%          (3) the smallest singular value of B11 is  on order of
%              the kth singular value of A. 
%          (4) L22*U22 is on order of the (k+1)st singular value
%              of A
%   The default value (if lupan3 has 2 or fewer parameters or rcond1=0)
%   is rcond1 = eps (relative machine precision).

%  Note:  Optional parameter mu is not used in this implementation
%
%   Reference: "On the existance and computation of rank-revealing
%   LU factorizations", Linear Algebra and its Applications 316 (2000)
%   p. 199-222 by C.-T. Pan.


if ( nargin == 1)
   rcond1 = eps;
   nb = 0;
   mu = 1.1;
elseif (nargin == 2)
   rcond1 = eps;
   mu = 1.1;
elseif (nargin == 3)
   mu = 1.1;
end

% Note: if nb is input to zero qrxp will choose the default (ILAENV)
%       value of 32

[m,n]=size(A);
minmn = min(m,n);


%Pan's Algorithm 3:


%Step 1:  Move linearly independent columns of A to the left.
% We will replace Pan's Algorithm 2 with a rank revealing 
% QR factorization. In exact arithmetic this should produce
% the same permutation P as Pan's LU based algorithm 3. In
% inexact arithmetic the QR based algorithm should, in fact,
% be more accurate that Pan's LU algorithm since the LU
% algorithm forms the normal equations.  This will lose information
% about small - those less than sqrt(eps)*norm(A) -singular values.
% A natural choice of a QR algorithm would be to use qryp since
% this is based on Pan's QR algorithm.  However the qryp 
% implementation in ACM Algorithm 782 appears to have a
% bug, except that nb = 1 seems ok (but slower).  Therefore
% we have used qrxp from ACM Algorithm 782.

trep=0;    % in qrxp don't do repetitions (which can help timing of qrxp)
lda = 0;   % use default lda = declared number of rows in A ( = m)
[Qt,R,Q, time,rank1,crcond,svlues] = qrxp(A,nb,trep,rcond1,lda);
k = rank1;
B = permright(A,Q);   % forms A*Q


% Step 2: Move the linearly independent rows of the first k columns
% of A*Q to the first k rows.
% Rather than Algorithm 1 described in Pan's paper we will use
% qrxp.  This is faster than a Matlab implementation of Pan's
% Algorithm 1 since qrxp calls Fortran code and qrxp serves
% the same purpose.

[Qt,R,P]=qrxp(B(:,1:k)',nb,trep,0);  
P = P';
B = permleft(B,P);    % so B = P * B;


% Step 3:  Construct the rank revealing LU factorization of
% B = P*A*Q with the k by k upper left corner well conditioned, if possible

B11 = B(1:k,1:k);
B12 = B(1:k,k+1:n);
B21 = B(k+1:m,1:k);
B22 = B(k+1:m,k+1:n);
% factor B11 and construct portions of L and U
[L11,U11,P11]=lu(B11);
U12 = L11 \ ( permleft(B12,P11) );    %U12 = L11 \ (P11*B12);
L21 = B21 / U11;
U22 = B22 - L21 * U12;    % schur complement
% factor U22
[L22,U22,P22]=lu(U22);
% update L21 to reflect P22 permutation
L21=permleft(L21,P22);
% update P to reflect P11 and P12 permutations
P=permleft(P,[P11,zeros(k,m-k);zeros(m-k,k),P22]);
% construct L and U
L = [L11, zeros(k,minmn-k) ; L21  L22 ];
U = [U11, U12 ; zeros(minmn-k,k), U22 ];

%--------------------------------------------

function B = permright(A,P)
% Form, efficiently, A*P where P is a permutation matrix
% Usage B = permright(A,P);
% Input:  A  - an m by n matrix
%         P  - an n by n permutation matrix 
% Output: B = A * P
if ( length(P) > 0 )
  [pv,pvt]=find(P==1);
  B=A(:,pv);
else
  B = A ;
end

%--------------------------------------------

function B = permleft(A,P)
% Form, efficiently, P*A where P is a permutation matrix
% Usage B = permleft(A,P);
% Input:  A  - an m by n matrix
%         P  - an m by m permutation matrix 
% Output: B = P * A
if (length(P) > 0 )
  [pv,pvt]=find(P'==1);
  B=A(pv,:);
else
  B = A;
end