#! /usr/bin/env python
#
def 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.
#
#  Licensing:
#
#    This code is distributed under the GNU LGPL license.
#
#  Modified:
#
#    31 August 2016
#
#  Author:
#
#    Original FORTRAN77 version by Jack Dongarra, Cleve Moler, Jim Bunch, 
#    Pete Stewart.
#    Python 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.
#
  import numpy as np
  from dqrdc import dqrdc

  jpvt = np.zeros ( n, dtype = np.int32 )
  job = 1

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

  kr = 0
  k = min ( m, n )

  for j in range ( 0, k ):
    if ( abs ( a[j,j] ) <= tol * abs ( a[0,0] ) ):
      break
    kr = j + 1

  return kr, jpvt, qraux, a

def dqrank_test ( ):

#*****************************************************************************80
#
## DQRANK_TEST tests DQRANK.
#
#  Discussion:
#
#    DQRANK factors a rectangular matrix and estimates its rank.
#
#  Licensing:
#
#    This code is distributed under the GNU LGPL license.
#
#  Modified:
#
#    31 August 2016
#
#  Author:
#
#    John Burkardt
#
  import numpy as np
  import platform
  from r8mat_print import r8mat_print
  from r8vec_print import r8vec_print

  m = 4
  n = 4
  lda = n
  tol = 0.000001

  print ( '' )
  print ( 'DQRANK_TEST' )
  print ( '  Python version: %s' % ( platform.python_version ( ) ) )
  print ( '  DQRANK factors a rectangular matrix and estimates its rank.' )

  a = np.zeros ( [ m, n ] )

  for i in range ( 0, m ):
    for j in range ( 0, n ):
      if ( i == j ):
        if ( i == 0 or i == m - 1 ):
          a[i,j] = 1.0 + 0.000000001
        else:
          a[i,j] = 2.0
      elif ( i == j + 1 or i == j - 1 ):
        a[i,j] = -1.0

  r8mat_print ( m, n, a, '  Matrix A:' )
        
  kr, jpvt, qraux, a = dqrank ( a, lda, m, n, tol )

  print ( '' )
  print ( ' Matrix rank estimated to be %d' % ( kr ) )

  r8mat_print ( m, n, a, '  QR factorization information in A:' )

  r8vec_print ( m, qraux, '  QR factorization informaiton in QRAUX:' )
#
#  Terminate.
#
  print ( '' )
  print ( 'DQRANK_TEST' )
  print ( '  Normal end of execution.' )
  return

if ( __name__ == '__main__' ):
  from timestamp import timestamp
  timestamp ( )
  dqrank_test ( )
  timestamp ( )