function [Q,R,P, time ] =  qrnp(A,nb,trep,rcond1,lda, nx)
%
%   qrnp computes the QR decomposition of A without column interchanges
%   using LAPACK 3.0 routine DGEQRF.
%
%   R = qrnp(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'.  P = I since no pivoting is done.
% 
%   [Q, R, P] = qrnp(A) returns Householder transformations in Q,
%   a triangular matrix R and permutation matrix P=I such that A * P is
%   triangularized by the Householder transformations stored in Q.
%   Note that A * P is not equal to Q*R since the Householder 
%   transformations are not multiplied out so that A * P = Q*R
%   is true.
%
%   For the call qrnp(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 qrnp(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 LAPACK 3.0's DGEQRF.
%
%   [R ,time] = qrnp(A,nb); or [Q, R, P,time] = qrnp(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 qrnp (with nothing to the left of qrnp) will return the 
%   factorization time in the standard Matlab parameter ans. Therfore after
%   "qrnp(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, qrnp(A,nb); toc, time = ans,"
%   should report the elapsed time (from toc) close to the variable time. 
%
%   For the call
%       [R,time] = qrnp(A,nb,trep); OR qrnp(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 qrnp(A, nb, trep, rcond1) rcond1 is not used and is
%   included so that the arguments are similar to those of other calls. 
%
%   For the call qrnp(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.

%   Written by L. Foster. San Jose State University. Fortran code
%   from LAPACK 3.0.

%   Additional features and comments:
%
%   For the call qrnp(A,nb,trep,rcond1,lda, nx) the algorithm 
%   uses a block size of nb and a cross over point of nx.  When the
%   min(m,n) <= nx unblocked code is used.  Also unblocked code is used
%   in the factorization when either the number of rows or columns 
%   remaining to be factored becomes smaller than nx. The default
%   size is nx = 128 (chosen when there are less than 6 input parameters 
%   or if nx = 0) which corresponds to LAPACK 3.0's DGEQP3. 
%   The choice nx = nb also appears to be a reasonable value.  

%   Note that Q is not explictly formed in order not to slow down qrnp. 
%   One could write a routine (or use a LAPACK routine) that multiplies
%   by a Q written in factored form but I have not done so.
%   Even though Q is not constructed, after R = qrnp(A) 
%   R'*R should be numerically equal to A'*A.  Also the least squares
%   problem min || b - A x || can be solved by x = R \ (R' \ (A'*b ));
%   w = R \ (R' \ (A'*(b-A*x))); x = x + w.  Usually this solution is
%   accurate.  See Bjork, Numerical Solution for Least Squares Problems, 
%   SIAM, 1996, pp 250

%   Also even though Q is not constructed, after [Q, R, P] = qrnp(A)
%   P*R'*R*P' should be numerically equal to A'*A.  Also the least squares
%   problem min || b - A x || can be solved by : 
%   x = P * (R \ (R' \ (P' * (A'*b ))));
%   x = P * (R \ (R' \ (P' * (A'*(b-A*x) )))); x = x + w.
%   Usually this solution is accurate.

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

if ( nargin == 1 )
   nb = 0;         %reset below
   trep = 0.0;
   rcond1 = 0.0;
   lda = m;
   nx = 0;         %reset below
end
if ( nargin == 2 )
   trep = 0.0;
   rcond1 = 0.0;
   lda = m;
   nx = 0;         %reset below
end
if (nargin == 3 )
   rcond1 = 0.0;
   lda = m;
   nx = 0;         %reset below
end
if (nargin == 4 )
   lda = m;
   nx = 0;         %reset below
end
if (nargin == 5 )
   nx = 0;         %reset below
end

if ( lda < m )
    lda = m;
end
if ( nb == 0 )    % if nb == 0 choose the default nb
%  the following choice of nb matches the choice of LAPACK 3.0's ILAENV
%     (and the logic of DGEQP3)
   if ( minmn <= 128 ) 
      nb = 1;
   else
      nb = 32;
   end
end
if ( nx == 0 )    % if nx == 0 choose default nx
%  the following choice of nx matches the choice of LAPACK 3.0's ILAENV
   nx = 128;
%   nx = nb;
end   

ldatemp= zeros(lda,1);   % used to get the lda into the mex files
nbtemp = zeros(nb, 1);   % used to get the nb into the mex files
nxtemp = zeros(nx, 1);   % used to get the nx into the mex files
  
[B,jp,time2,jpt,work,work2] = ...
          qrnpf(A,nbtemp,nxtemp,trep,ldatemp);
if (nargout ~= 0 )          
   B = B(1:m,1:n); 
   jp = round(jp);
   minmn = min(m,n);
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 )
    Q = tril(B,-1);
    Q(1:minmn,1:minmn) = Q(1:minmn,1:minmn)+diag(work2(1:minmn));
    R = triu(B);
    P = eye(n,n);
    P = P(:,jp);
elseif ( nargout == 4 )
    Q = tril(B,-1);
    Q(1:minmn,1:minmn) = Q(1:minmn,1:minmn)+diag(work2(1:minmn));
    R = triu(B);
    P = eye(n,n);
    P = P(:,jp);
    time = time2;
else
    disp('must have 1, 2, 3, or 4 output arguments')
end  
