function [L,U,P,Q, time, comparisons, element_switches, ...
    flops_unbl, extra_flops, row_cols_searched, row_col_switches] = ...
    lupp(A,nb,trep,tol,lda)    
%
%   lupp computes an LU decomposition of A with partial pivoting using
%   LAPACK 3.0 routine DGETRF.
%
%   [L,U] = lupp(A) stores an "psychologically upper triangular matrix"
%   (i.e. a product of lower triangular and permutation matrices) in U 
%   and a "psychologically lower triangular matrix" (i.e. a product
%   of lower triangular and permutation matrices) in L, so
%   that A = L*U. 
% 
%   [L,U,P,Q] = lupp(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.  Q = I for partial pivoting.
%
%   For the call lupp(A,nb) the algorithm uses a block size of nb. 
%   If nb = 0 or the call is lupp(A) the default block size
%   is selected which depends on minmn = min(m,n):
%      for minmn <= 64, nb = 1
%      for 64 < minmn, nb = 64.
%   These values may not be optimal for your computer.
%   They correspond to the default values in LAPACK's DGETRF.
%
%   [L,U,time] = lupp(A,nb); or [L,U,P,Q,time] = lupp(A,nb); will 
%   return the factorization time as well as L, U, P and Q. The time
%   returned is for the Fortran portion of the code. A call
%   to lupp (with nothing to the left of lupp) will return the 
%   factorization time in the standard Matlab parameter ans. Therfore after
%   "lupp(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, lupp(A,nb); toc, time = ans,"
%   should report the elapsed time (from toc) close to the variable time. 
%
%   For the call
%       [L,U,time] = lupp(A,nb,trep); OR lupp(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 lupp(A,nb,trep,tol) tolerance pivoting will be used.
%   Only pivot when the magnitude of diagonal entry < (magnitude
%   of largest entry below the diagonal ) / tol.  The largest
%   magnitude entry of L will be <= tol. Note if tol < 1  a value
%   of tol = 1 is used.
%
%   For the call lupp(A, nb, trep, tol, lda) the matrix A is stored,
%   in the Fortran code, as a partial filled array with lda rows. 
%   The choice of lda can affect the speed of the algorithm. Some 
%   languages such as Fortran typically use partially filled arrays.
%
%   The call 
%       [L,U,P,Q, time, comparisons, element_switches, flops_unbl,...
%       extra_flops, row_cols_searched, row_col_switches ] =  ...
%       lupp(A,nb);
%    also returns the number of comparison in the factorization,
%    the number of element switches, the flops in the unblocked
%    algorithm, the extra flops due to the blocked algorithm, 
%    the number of rows and columns searched, and the number of 
%    row and column switches.  These counts are approximate and are
%    better for larger m and n.

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

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

if ( nargin == 1 )
   nb = 0;    % reset below
   trep = 0.0;
   tol = 1;
   lda = m;
end
if ( nargin == 2 )
    trep = 0.0;
    tol = 1;
   lda = m;
end
if (nargin == 3 )
   tol = 1;
   lda = m;
end  
if ( nargin == 4 )
   lda = m;
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 LAPACK 3.0's DGETRF)
   if ( minmn <= 64 ) 
      nb = 1;
   else
      nb = 64;
   end
end

ldatemp = zeros(lda,1);
nbtemp = zeros(1,1);  % nbtemp not used in current version of lupp 
  
[B,ip,jp,time2,rncomp,rnswch,rneflp,rnsear,rnswrc,ipt,jpt,work] = ...
          luppf(A,nb,trep,tol,ldatemp,nbtemp);

if (nargout ~= 0 )          
   B = B(1:m,1:n); 
   ip = ip(1:m);
   ip = round(ip);
   jp = round(jp);
   minmn = min(m,n);
end   

if (nargout == 0 )
   L = time2;
elseif (nargout == 1 )
    L = B;
elseif ( nargout == 2 )
    L = tril(B,-1);
    U = triu(B);
    clear B ;            % to reduce memory demands
    L = L(:,1:minmn);
    U = U(1:minmn,:);
    L = L + eye(m,minmn);
    rip(ip(1:m)) = 1:m;
    L = L(rip,:);
    rjp(jp(1:n)) = 1:n ;
    U = U(:,rjp);
elseif (nargout == 3 )
    L = tril(B,-1);
    U = triu(B);
    clear B ;            % to reduce memory demands    
    L = L(:,1:minmn);
    U = U(1:minmn,:);
    L = L + eye(m,minmn);
    rip(ip(1:m)) = 1:m;
    L = L(rip,:);
    rjp(jp(1:n)) = 1:n ;
    U = U(:,rjp);
    P= time2;
elseif ( nargout == 4 | nargout == 5 | nargout == 11)
    L = tril(B,-1);
    U = triu(B);
    clear B ;            % to reduce memory demands
    L = L(:,1:minmn);
    U = U(1:minmn,:);
    L = L + eye(m,minmn);
    P = eye(m,m);
    P = P(ip,:);
    Q= eye(n,n);
    Q = Q(:,jp);
else
    disp('must have 1, 2, 3, 4 , 5 or 11 output arguments')
end  

if (nargout == 5 )
    time = time2;
end    

if (nargout == 11 )
    time = time2;
    comparisons = rncomp;
    element_switches = rnswch;
    flops_unbl=n*(2*m-n-1)/2+n*(n-1)*(3*m-n-1)/3;  % for m >= n
    if ( n > m ) flops_unbl = m*(m-1)/2+m*(m-1)*(3*n-m-1)/3; end;    
    extra_flops = rneflp;
    row_cols_searched = rnsear;
    row_col_switches = rnswrc; 
end