function [Q,R,P, time, rank1, crcond, svlues] =  ...
                   qrxp2(A,nb,trep, rcond1,lda, job, nQ)
%   qrxp computes the QR decomposition of A with column interchanges
%   using DGEQPX in algorithm 782 from the collected algorithms of the ACM.
%
%   R = qrxp(A) stores an "psychologically upper triangular matrix"
%   (i.e. a product of upper triangular and permutation matrices) in R.
%   If A * P = Q * R1 where P is a permutation matrix and Q is 
%   orthogonal then R is R1 * P'.
% 
%   [Q, R, P] = qrxp(A) returns a triangular matrix R such A * P is
%   triangularized by the an orthogonal transformation.  Usually
%   (unless job = 1 or 2, see below) Q = zero and is 
%   used as a place keeper.
%
%   For the call qrxp(A,nb) the algorithm uses a block size of nb. 
%   A block size of nb > 1 will often be faster than nb = 1.
%   If nb = 0 or the call is qrxp(A) the default block size
%   is selected which depends on minmn = min(m,n):
%      for minmn <= 128, nb = 1
%      for 128 < minmn, nb = 32.
%   These values may not be optimal for your computer.
%   They correspond to the default values in DGEQPX.
%
%   [R ,time] = qrxp(A,nb); or [Q, R, P,time] = qrxp(A,nb); will
%   return the factorization time as well as R or Q, R and P. The time
%   returned is for the Fortran portion of the code.  A call
%   to qrxp (with nothing to the left of qrp) will return the factor-
%   ization time in the standard Matlab parameter ans.  Therefore after
%   "qrxp(A,nb); time = ans;" time will contain the factorization time.
%   For this last call the time reported should be close to the time
%   reported by tic and toc in Matlab so "tic, qrxp(A,nb); toc, time = ans,"
%   should report the elapsed time (from toc) close to the variable time. 
%
%   For the call
%       [R,time] = qrxp(A,nb,trep); OR  qrxp(A,nb,trep); time=ans;
%   the factorization will be done multiple times to improve the 
%   accuracy of the time determination according to:
%        0 < trep    do repeatedly until the elapsed time exceeds trep
%       -1 <= trep <= 0   do only once
%        trep < -1  do |trep| repetitions
%   For smaller matrices (less than about 500) one should have trep = 1
%   or greater to time accurately. The time returned is the time for one 
%   factorization (=total time of factorizations / number of repetitions).
%
%   For the call qrxp(A, nb, trep, rcond1) the  columns of A are 
%   permuted so that R11, the leading rank1 columns of R,
%   in the QR factorization of A*P, have 
%       estimated condition number of R11 <= 1 / rcond1 
%   and this condition is violated if any additional columns of A are
%   factored, increasing the size of R11.    The default value (if qrxp
%   has < 4 parameters or rcond1=0) is rcond1 = eps (relative machine
%   precision).
%
%   [Q, R, P,time, rank1] = qrxp(A,nb,trep,rcond1); will 
%   also return the calculated rank, rank1, corresponding to the
%   stopping criteria rcond1.
%
%   For the call qrxp(A, nb, trep, rcond1, lda) the matrix A is 
%   stored, in the Fortran code, as a partial filled array with lda rows.
%   If lda < m then lda is reset to m. The choice of lda can affect the 
%   speed of the algorithm. Some languages such as Fortran typically
%   use partially filled arrays.
%
%   The call qrxp(A,nb,trep, rcond1,lda, job) the output matrix 
%   Q contains an orthogonal matrix used in A*P=Q*R (for job =3)
%   or (for job = 2) the transpose of this matrix. If job = 1 (the default) 
%   Q is not formed.
%
%   The call qrxp(A,nb,trep,rcond1,lda, job, nQ) the output matrix 
%   Q is the first nQ rows of Q (for job =2 or job = 3) apparently
%   calculated by a different manner. If job = 1 (the default) Q is 
%   not formed.
%
%   [Q,R,P, time,rank1,crcond,svlues] = qrxp(A,nb,trep,rcond1,lda,job); 
%   also calculates and returns
%       crcond - 1 / crcond is estimated condition number of R11,
%       svlues - estimated singular values where     
%          SVLUES(1): largest singular value of R(1:rank1,1:rand1)
%          SVLUES(2): smallest singular value of R(1:rank1,1:rank1)
%          SVLUES(3): smallest singular value of R(1:rank1+1,1:rank1+1)
%          SVLUES(4): smallest singular value of R

%   Written by L. Foster, San Jose State University.  The Fortan code is
%   from "Algorithm 782: Codes for rank-revealing qr factorizations 
%   of dense matrices." ACM Transactions on Mathematical Software 24, 
%   p. 254-257, 1998, by Bischof, C. H. and Quintana-Orti, G.

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

if ( nargin == 1 )
   nb = 0;         
   trep = 0.0;
   rcond1 = 0.0;
   lda = m;
   job = 1;
   nQ = 1;
end
if ( nargin == 2 )
   trep = 0.0;
   rcond1 = 0.0;
   lda = m;
   job = 1;
   nQ = 1;
end
if (nargin == 3 )
   rcond1 = 0.0;
   lda = m;
   job = 1;
   nQ = 1;
end
if (nargin == 4 )
    lda = m;
    job = 1;
    nQ = 1;  
end
if (nargin == 5 )
   job = 1;
   nQ = 1;
end
if (nargin == 6)
   nQ = m;     
end  
if ( nQ > m )
   nQ = m;
end

if ( lda < m )
    lda = m;
end 

% Note if nb == 0 the value of nb is reset in qrxpb.f to be
%   the default value of nb selected by DGEQPX (nb=32 for minmn >= 128,
%   nb = 1 for minmn < 128).  The interface qrxpg5.f sets aside
%   enough memory

if (job == 1 )
    rowsQ = lda;
    colsQ = 1;
elseif (job == 2 )
    rowsQ = lda;
    colsQ = nQ;
elseif (job == 3 )
    rowsQ = nQ;
    colsQ = lda;
else
    disp('job must be 1, 2 or 3')
end    

%Q = zeros(rowsQ,colsQ);
Q = sparse(rowsQ,colsQ);
if (rowsQ >= colsQ )
   Q(1:colsQ,1:colsQ) = speye(colsQ,colsQ);%eye(colsQ,colsQ);
else
   Q(1:rowsQ,1:rowsQ) = speye(rowsQ,rowsQ);%eye(rowsQ,rowsQ);
   % Q, as input, could be changed to other matrices
end
svlues = zeros(4,1);

ldatemp= zeros(lda,1);   % used to get the lda into the mex files
if ( nb > 0 )
    nbtemp = zeros(nb, 1);   % used to get the nb into the mex files
else
    nbtemp = zeros(abs(nb), 2);   % used to get the nb into the mex files
end
  
[B, jp,time2,crcond,rank1,svlues,jpt, work] = ...
          qrxpf(A,nbtemp,rcond1,trep,ldatemp,job,Q);

if (nargout ~= 0 )          
   B = B(1:m,1:n); 
   jp = round(jp);
   minmn = min(m,n);
   if (job == 1 )
       Q=0;
   end
   if ( job == 2 )
       Q = Q(1:m,1:nQ);
       Q=Q';
   end
   if ( job == 3 )
       Q= Q(1:nQ,1:m);
   end   
end   

if (nargout == 0 )
    Q = time2;
elseif (nargout == 1 )
    Q = triu(B);
    rjp(jp(1:n)) = 1:n ;
    Q = Q(:,rjp);
elseif ( nargout == 2 )
    Q = triu(B);
    rjp(jp(1:n)) = 1:n ;
    Q = Q(:,rjp);
    R = time2;
elseif (nargout == 3 )
    R = triu(B);
    P = speye(n,n);
%   P = eye(n,n);
    P = P(:,jp);
elseif ( nargout == 4 )
    R = triu(B);
    P = speye(n,n);    
%    P = eye(n,n);
    P = P(:,jp);
    time = time2;
elseif (nargout <=  7)
    R = triu(B);
    P = speye(n,n);
%   P = eye(n,n);
    P = P(:,jp);
    time = time2;
else
    disp('must have 1, 2, 3, 4, 5, 6 or 7 output arguments')
end  
