#! /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 ( )