#! /usr/bin/env python
#
def dqrdc ( a, lda, n, p, jpvt, job ):
#*****************************************************************************80
#
## DQRDC computes the QR factorization of a real rectangular matrix.
#
# Discussion:
#
# DQRDC uses Householder transformations.
#
# Column pivoting based on the 2-norms of the reduced columns may be
# performed at the user's option.
#
# Licensing:
#
# This code is distributed under the GNU LGPL license.
#
# Modified:
#
# 30 August 2016
#
# Author:
#
# Original FORTRAN77 version by Dongarra, Moler, Bunch, Stewart.
# Python version by John Burkardt.
#
# Reference:
#
# Dongarra, Moler, Bunch, Stewart,
# LINPACK User's Guide,
# SIAM, (Society for Industrial and Applied Mathematics),
# 3600 University City Science Center,
# Philadelphia, PA, 19104-2688.
# ISBN 0-89871-172-X
#
# Parameters:
#
# Input, real A(LDA,P), the N by P matrix whose decomposition is to
# be computed.
#
# Input, integer LDA, the leading dimension of the array A. LDA must
# be at least N.
#
# Input, integer N, the number of rows of the matrix A.
#
# Input, integer P, the number of columns of the matrix A.
#
# Input, integer JPVT(P), integers that control the selection of the pivot
# columns. The K-th column A(*,K) of A is placed in one of three classes
# according to the value of JPVT(K).
# > 0, then A(K) is an initial column.
# = 0, then A(K) is a free column.
# < 0, then A(K) is a final column.
# Before the decomposition is computed, initial columns are moved to
# the beginning of the array A and final columns to the end. Both
# initial and final columns are frozen in place during the computation
# and only free columns are moved. At the K-th stage of the
# reduction, if A(*,K) is occupied by a free column it is interchanged
# with the free column of largest reduced norm. JPVT is not referenced
# if JOB == 0.
#
# Input, integer JOB, initiates column pivoting.
# 0, no pivoting is done.
# nonzero, pivoting is done.
#
# Output, real A(LDA,P), contains in its upper triangle the upper
# triangular matrix R of the QR factorization. Below its diagonal
# A contains information from which the orthogonal part of the
# decomposition can be recovered. Note that if pivoting has been
# requested, the decomposition is not that of the original matrix A
# but that of A with its columns permuted as described by JPVT.
#
# Output, real QRAUX(P), contains further information required
# to recover the orthogonal part of the decomposition.
#
# Output, integer JPVT(P), JPVT(K) contains the index of the column of the
# original matrix that has been interchanged into the K-th column, if
# pivoting was requested.
#
import numpy as np
pl = 1
pu = 0
qraux = np.zeros ( p )
work = np.zeros ( p )
#
# If pivoting is requested, rearrange the columns.
#
if ( job != 0 ):
for j in range ( 0, p ):
swapj = ( 0 < jpvt[j] )
if ( jpvt[j] < 0 ):
jpvt[j] = - ( j + 1 )
else:
jpvt[j] = ( j + 1 )
if ( swapj ):
if ( j + 1 != pl ):
for i in range ( 0, n ):
t = a[i,pl-1]
a[i,pl-1] = a[i,j]
a[i,j] = t
jpvt[j] = jpvt[pl-1]
jpvt[pl-1] = j + 1
pl = pl + 1
pu = p
for j in range ( p - 1, -1, -1 ):
if ( jpvt[j] < 0 ):
jpvt[j] = - jpvt[j]
if ( j + 1 != pu ):
for i in range ( 0, n ):
t = a[i,pu-1]
a[i,pu-1] = a[i,j]
a[i,j] = t
jp = jpvt[pu-1]
jpvt[pu-1] = jpvt[j]
jpvt[j] = jp
pu = pu - 1
#
# Compute the norms of the free columns.
#
for j in range ( pl - 1, pu ):
t = 0.0
for i in range ( 0, n ):
t = t + a[i,j] ** 2
qraux[j] = np.sqrt ( t )
work[j] = qraux[j]
#
# Perform the Householder reduction of A.
#
lup = min ( n, p )
for l in range ( 0, lup ):
#
# Bring the column of largest norm into the pivot position.
#
if ( pl <= l + 1 and l + 1 < pu ):
maxnrm = 0.0
maxj = l
for j in range ( l, pu ):
if ( maxnrm < qraux[j] ):
maxnrm = qraux[j]
maxj = j
if ( maxj != l ):
for i in range ( 0, n ):
t = a[i,l]
a[i,l] = a[i,maxj]
a[i,maxj] = t
qraux[maxj] = qraux[l]
work[maxj] = work[l]
jp = jpvt[maxj]
jpvt[maxj] = jpvt[l]
jpvt[l] = jp
#
# Compute the Householder transformation for column L.
#
qraux[l] = 0.0
if ( l + 1 != n ):
t = 0.0
for i in range ( l, n ):
t = t + a[i,l] ** 2
nrmxl = np.sqrt ( t )
if ( nrmxl != 0.0 ):
if ( a[l,l] < 0.0 ):
nrmxl = - abs ( nrmxl )
elif ( 0.0 < a[l,l] ):
nrmxl = abs ( nrmxl )
for i in range ( l, n ):
a[i,l] = a[i,l] / nrmxl
a[l,l] = 1.0 + a[l,l]
#
# Apply the transformation to the remaining columns, updating the norms.
#
for j in range ( l + 1, p ):
t = 0.0
for i in range ( l, n ):
t = t - a[i,l] * a[i,j]
t = t / a[l,l]
for i in range ( l, n ):
a[i,j] = a[i,j] + t * a[i,l]
if ( pl <= j and j <= pu ):
if ( qraux[j] != 0.0 ):
tt = 1.0 - ( abs ( a[l,j] ) / qraux[j] ) ** 2
tt = max ( tt, 0.0 )
t = tt
tt = 1.0 + 0.05 * tt * ( qraux[j] / work[j] ) ** 2
if ( tt != 1.0 ):
qraux[j] = qraux[j] * np.sqrt ( t )
else:
t = 0.0
for i in range ( l + 1, n ):
t = t + a[i,j] ** 2
qraux[j] = np.sqrt ( t )
work[j] = qraux[j]
#
# Save the transformation.
#
qraux[l] = a[l,l]
a[l,l] = - nrmxl
return a, qraux, jpvt
def dqrdc_test ( ):
#*****************************************************************************80
#
## DQRDC_TEST tests DQRDC.
#
# Discussion:
#
# DQRDC computes the QR factorization.
#
# Licensing:
#
# This code is distributed under the GNU LGPL license.
#
# Modified:
#
# 30 August 2016
#
# Author:
#
# John Burkardt
#
import numpy as np
import platform
from dqrsl import dqrsl
n = 3
p = 3
lda = n
print ( '' )
print ( 'DQRDC_TEST' )
print ( ' Python version: %s' % ( platform.python_version ( ) ) )
print ( ' DQRDC computes the QR decomposition of a rectangular' )
print ( ' matrix, but does not return Q and R explicitly.' )
print ( '' )
print ( ' Show how Q and R can be recovered using SQRSL.' )
#
# Set the matrix A.
#
a = np.array \
( [ \
[ 1.0, 1.0, 0.0 ], \
[ 1.0, 0.0, 1.0 ], \
[ 0.0, 1.0, 1.0 ] \
] )
print ( '' )
print ( ' The original matrix A:' )
print ( '' )
for i in range ( 0, n ):
for j in range ( 0, p ):
print ( ' %14.6g' % ( a[i,j] ) ),
print ( '' )
#
# Decompose the matrix.
#
print ( '' )
print ( ' Decompose the matrix.' )
job = 0
ipvt = np.zeros ( p, dtype = np.int32 )
a, qraux, ipvt = dqrdc ( a, lda, n, p, ipvt, job )
#
# Print out what DQRDC has stored in A...
#
print ( '' )
print ( ' The packed matrix A which describes Q and R:' )
print ( '' )
for i in range ( 0, n ):
for j in range ( 0, p ):
print ( ' %14.6g' % ( a[i,j] ) ),
print ( '' )
#
# ...and in QRAUX.
#
print ( '' )
print ( ' The QRAUX vector, containing some additional' )
print ( ' information defining Q:' )
print ( '' )
for i in range ( 0, n ):
print ( ' %14.6g' % ( qraux[i] ) )
print ( '' )
#
# Print out the resulting R factor.
#
r = np.zeros ( [ n, p ] )
print ( '' )
print ( ' The R factor:' )
print ( '' )
for i in range ( 0, n ):
for j in range ( 0, p ):
if ( i <= j ):
r[i,j] = a[i,j]
print ( ' %14.6g' % ( r[i,j] ) ),
print ( '' )
#
# Call DQRSL to extract the information about the Q matrix.
# We do this, essentially, by asking DQRSL to tell us the
# value of Q*Y, where Y is a column of the identity matrix.
#
q = np.zeros ( [ n, n ] )
job = 10000
for j in range ( 0, n ):
#
# Set the vector Y.
#
y = np.zeros ( n )
y[j] = 1.0
#
# Ask DQRSL to tell us what Q*Y is.
#
qy, qty, b, rsd, xb, info = dqrsl ( a, lda, n, p, qraux, y, job )
if ( info != 0 ):
print ( ' Error! DQRSL returns INFO = %d' % ( info ) )
return
#
# Copy QY into the appropriate column of Q.
#
for i in range ( 0, n ):
q[i,j] = qy[i]
#
# Now print out the Q matrix we have extracted.
#
print ( '' )
print ( ' The Q factor:' )
print ( '' )
for i in range ( 0, n ):
for j in range ( 0, n ):
print ( ' %14.6g' % ( q[i,j] ) ),
print ( '' )
#
# Compute Q*R to verify that it equals A.
#
b = np.dot ( q, r )
#
# Print the result.
#
print ( '' )
print ( ' The product Q * R:' )
print ( '' )
for i in range ( 0, n ):
for j in range ( 0, p ):
print ( ' %14.6g' % ( b[i,j] ) ),
print ( '' )
#
# Terminate.
#
print ( '' )
print ( 'DQRDC_TEST' )
print ( ' Normal end of execution.' )
return
if ( __name__ == '__main__' ):
from timestamp import timestamp
timestamp ( )
dqrdc_test ( )
timestamp ( )