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