function [L,U,P,Q, time, comparisons, element_switches, ...
    flops_unbl, extra_flops, row_cols_searched, row_col_switches] = ...
    lucp(A,nb,trep,tol,lda)
%
%   lucp computes an LU decomposition of A with complete pivoting using 
%   a modification of LAPACK 3.0 routine DGETRF.
%
%   [L,U] = lucp(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] = lucp(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.
%
%   For the call lucp(A) or lucp(A,nb) the algorithm uses a block
%   size nb = 1. The value of the input parameter nb is ignored.
%   The parameter is in the  call to keep consistency with the
%   format for lurp, lupp, etc.
%
%   [L,U,time] = lucp(A,nb); or [L,U,P,Q,time] = lucp(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 lucp (with nothing to the left of lucp) will return the 
%   factorization time in the standard Matlab parameter ans. Therfore after
%   "lucp(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, lucp(A,nb); toc, time = ans,"
%   should report the elapsed time (from toc) close to the variable time. 
%
%   For the call
%       [L,U,time] = lucp(A,nb,trep); OR lucp(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 lucp(A,nb,trep,tol) tolerance pivoting will be used.
%   Only pivot when the magnitude of diagonal entry < (magnitude
%   of largest entry in lower right corner ) / tol.  The largest
%   magnitude entry of L and of D^(-1) U, where D is the diagonal of U
%   will be <= tol. Note if tol < 1 a value of tol = 1 is used.
%
%   For the call lucp(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 ] =  ...
%       lucp(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 is a modification of code in LAPACK 3.0.

[m,n]=size(A);

if ( nargin == 1 )
   nb = 1;
   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

nb =1;    % always have a block size 1 for lucp

ldatemp = zeros(lda,1);
  
[B,ip,jp,time2,rncomp,rnswch,rneflp,rnsear,rnswrc,ipt,jpt,work] = ...
          lucpf(A,nb,trep,tol,ldatemp);
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