function [ kr, jpvt, qraux, a ] = dqrank ( a, lda, m, n, tol )

%*****************************************************************************80
%
%% DQRANK QR-factors a rectangular matrix and estimates its rank.
%
%  Discussion:
%
%    This routine is used in conjunction with dqrlss to solve
%    overdetermined, underdetermined and singular linear systems
%    in a least squares sense.
%
%    DQRANK uses the LINPACK subroutine DQRDC to compute the QR
%    factorization, with column pivoting, of an M by N matrix A.
%    The numerical rank is determined using the tolerance TOL.
%
%    Note that on output, ABS ( A(1,1) ) / ABS ( A(KR,KR) ) is an estimate
%    of the condition number of the matrix of independent columns,
%    and of R.  This estimate will be <= 1/TOL.
%
%  Modified:
%
%    29 August 2016
%
%  Author:
%
%    Original FORTRAN77 version by Jack Dongarra, Cleve Moler, Jim Bunch, 
%    Pete Stewart.
%    MATLAB version by John Burkardt.
%
%  Reference:
%
%    Jack Dongarra, Cleve Moler, Jim Bunch, Pete Stewart,
%    LINPACK User's Guide,
%    SIAM, 1979,
%    ISBN13: 978-0-898711-72-1,
%    LC: QA214.L56.
%
%  Parameters:
%
%    Input/output, real A(LDA,N).  On input, the matrix whose
%    decomposition is to be computed.  On output, the information from DQRDC.
%    The triangular matrix R of the QR factorization is contained in the
%    upper triangle and information needed to recover the orthogonal
%    matrix Q is stored below the diagonal in A and in the vector QRAUX.
%
%    Input, integer LDA, the leading dimension of A, which must
%    be at least M.
%
%    Input, integer M, the number of rows of A.
%
%    Input, integer N, the number of columns of A.
%
%    Input, real TOL, a relative tolerance used to determine the
%    numerical rank.  The problem should be scaled so that all the elements
%    of A have roughly the same absolute accuracy, EPS.  Then a reasonable
%    value for TOL is roughly EPS divided by the magnitude of the largest
%    element.
%
%    Output, integer KR, the numerical rank.
%
%    Output, integer JPVT(N), the pivot information from DQRDC.
%    Columns JPVT(1), ..., JPVT(KR) of the original matrix are linearly
%    independent to within the tolerance TOL and the remaining columns
%    are linearly dependent.
%
%    Output, real QRAUX(N), will contain extra information defining
%    the QR factorization.
%
%    Output, real A(LDA,N), factorization information.
%
  jpvt = zeros ( n, 1 );
  job = 1;

  [ a, qraux, jpvt ] = dqrdc ( a, lda, m, n, jpvt, job );

  kr = 0;
  k = min ( m, n );

  for j = 1 : k
    if ( abs ( a(j,j) ) <= tol * abs ( a(1,1) ) )
      return
    end
    kr = j;
  end

  return
end

