function [L,U,P,Q, time, comparisons, element_switches, ...
    flops_unbl, extra_flops, row_cols_searched, row_col_switches] = ...
    lurp(A,nb,trep,tol,lda, ratio, lwork)
%
%   lurp computes an LU decomposition of A using rook pivoting.
%
%   [L,U] = lurp(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] = lurp(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. If the parameter nb < 0 (see below) P and Q
%   are permutation vectors so A(P,Q) = L*U
%
%   For the call lurp(A,nb) the algorithm uses a block size of nb. 
%   If nb <= 0 or the call is lurp(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 default values may not be optimal for your computer.
%   They correspond to the default values in LAPACK 3.0's DGETRF.
%
%   [L,U,time] = lurp(A,nb); or [L,U,P,Q,time] = lurp(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 lurp (with nothing to the left of lurp) will return the 
%   factorization time in the standard Matlab parameter ans. Therefore after
%   "lurp(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, lurp(A,nb); toc, time = ans,"
%   should report the elapsed time (from toc) close to the variable time. 
%
%   For the call
%       [L,U,time] = lurp(A,nb,trep); OR lurp(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 lurp(A,nb,trep,tol); threshold pivoting will be used.
%   if tol >= 1 do threshold pivoting on rows and columns.
%      Pivot only if:
%         magnitude of pivot entry < (magnitude of the largest
%         entry in the corresponding row or column ) / tol
%      The maximum magnitude entry of L and of D^(-1) U will be <= tol,
%      where D is the diagonal of U
%   if tol < -1 do threshold pivoting on rows only. Use 
%      standard pivoting (i.e. tol = 1 ) on columns and pivot 
%      in a row only if:
%         magnitude of pivot entry < (magnitude of the largest
%         entry in the corresponding row ) / | tol |
%      The maximum magnitude entry of L will be <= 1 and
%      the maximum magnitude entry of D^(-1)U will be <= | tol |
%   if -1 < tol <= 1 use tol = 1
%   if tol = -1 do not search any rows (partial pivoting)
%   if there are 3 or fewer input parameters use the default tol = 1
%
%   For the call lurp(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 default value m is chosen if lda < m or lurp has 4 or fewer
%   input parameters.
%
%   For the call lurp(A,nb,trep, tol, lda, ratio) the blocked algorithm 
%   will switch to block size 1 if the "excess" calculations due to the 
%   blocking algorithm becomes more than ratio times the calculations
%   with nb = 1. The default value 0.25 is a reasonable choice. A switch
%   to the unblocked algorithm should almost never occur with this ratio.
%   The default value is chosen if ratio <= 0 or lurp has 5 or fewer
%   input parameters.
%
%   For the call lurp(A, nb, trep, tol, lda, ratio, lwork) a workspace 
%   of size lwork is used in the algorithm. 
%   If lwork >= max (5, m + n + nb*n) this workspace is used to reduce
%      cache misses and speed up the algorithm for larger matrices
%      on higher performance processors.
%   If 0 < lwork < max(5,m+n) lwork is reset to its minimum value
%      max(5,m+n) and the algorithm does not use extra workspace
%      to reduce cache misses.
%   The default value max(5, m + n + nb*n) is chosen if lwork <= 0 or if 
%   lurp has 6 or fewer input parameters.
%
%   The call 
%       [L,U,P,Q, time, comparisons, element_switches, flops_unbl,...
%       extra_flops, row_cols_searched, row_col_switches ] =  ...
%       lurp(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. copr.

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

if ( nargin == 1 )
%  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
   trep = 0.0;
   tol = 1.0;
   lda = m;
   ratio = .25;
   lwork = max(5, m + n + nb*n);
end
if ( nargin == 2 )
   trep = 0.0;
   tol = 1.0;
   lda = m;
   ratio = .25;
   lwork = max(5, m + n + nb*n);
end
if (nargin == 3 )
   tol = 1.0;
   lda = m;
   ratio = 0.25;
   lwork = max(5, m + n + nb*n);
end
if (nargin == 4 )
   lda = m;
   ratio = 0.25;
   lwork = max(5, m + n + nb*n);
end
if (nargin == 5)
   ratio = 0.25;
   lwork = max(5, m + n + nb*n);
end
if (nargin == 6)
   lwork = max(5, m + n + nb*n);
end 

% set default values for zero parameters
fvec = 0;      % flag for storing permutations as vectors (default no)
if ( nb <= 0 )
%  the following choice of nb matches the choice of LAPACK 3.0's ILAENV
%     ( and the logic of LAPACK 3.0's DGETRF)
   if ( nb < 0 )
      fvec = 1;
   end
   if ( minmn <= 64 )
      nb = 1;
   else
      nb = 64;
   end
end
% tol and trep already ok
if ( lda < m ) 
   lda = m;
end
if ( ratio == 0 )
   ratio = 0.25;
end
if ( nargin <= 6 )
    lwork = max(5, m + n + nb*n );
else
   if ( lwork <= 0 )
      lwork = max(5, m + n + nb*n);
   end
   if ( lwork <= max(5, m+n) )
      lwork = max(5, m+n);
   end
end

ldatemp = zeros(lda,1);   % used to get lda into mex file
                          % in a form suitable for dimension statement

if ( nb > n) nb = n; end      % nb <= min(m,n)
if ( nb > m ) nb = m; end     % otherwise get error in dgemv parameter 6

nbtemp=zeros(nb,1);    % used to get nb into the mex file
                       %   in a form suitable for dimension statements

work=zeros(lwork,1);   % used to create space for work space                     

[L,U,P,Q,ip,jp,time2,rncomp,rnswch,rneflp,rnrcsr,rnrcsw,ipt,jpt,iwork] = ...
        lurpf(A,nb,trep, tol, ldatemp, ratio, nbtemp, work, nargout,fvec);
% note lurpf, with different values of tol, is called for both 
%      rook pivoting (lurp.m ) and threshold rook pivoting (lutp.m). 
%      lurpf calls the fortran version of rook pivoting
if (nargout ~= 0 )
    if (lda > m)
      if ( m >= n )
         L = L(1:m,1:n);
      end
    end
    if ( lda > n )
        if (n >= m )
            U = U(1:m,1:n);
        end
    end
end   


if ( nargout == 0 )
   L = time2;
end
if ( nargout == 3)
   P = time2;
end
if (nargout == 5 )
   time = time2;
end
if ( nargout >= 4 & fvec == 1 )
   % makre row vectors, consistent  Matlab's lu
   P = P';
   Q = Q';
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 = rnrcsr;
    row_col_switches = rnrcsw; 
elseif (nargout >= 6 )
    disp('must have 1, 2, 3, 4 , 5 or 11 output arguments')
end

