#! /usr/bin/env python
#
def dqrsl ( a, lda, n, k, qraux, y, job ):
#*****************************************************************************80
#
## DQRSL computes transformations, projections, and least squares solutions.
#
# Discussion:
#
# DQRSL requires the output of DQRDC.
#
# For K <= min(N,P), let AK be the matrix
#
# AK = ( A(JPVT(1)), A(JPVT(2)), ..., A(JPVT(K)) )
#
# formed from columns JPVT(1), ..., JPVT(K) of the original
# N by P matrix A that was input to DQRDC. If no pivoting was
# done, AK consists of the first K columns of A in their
# original order. DQRDC produces a factored orthogonal matrix Q
# and an upper triangular matrix R such that
#
# AK = Q * (R)
# (0)
#
# This information is contained in coded form in the arrays
# A and QRAUX.
#
# The parameters QY, QTY, B, RSD, and AB are not referenced
# if their computation is not requested and in this case
# can be replaced by dummy variables in the calling program.
# To save storage, the user may in some cases use the same
# array for different parameters in the calling sequence. A
# frequently occuring example is when one wishes to compute
# any of B, RSD, or AB and does not need Y or QTY. In this
# case one may identify Y, QTY, and one of B, RSD, or AB, while
# providing separate arrays for anything else that is to be
# computed.
#
# Thus the calling sequence
#
# call dqrsl ( a, lda, n, k, qraux, y, dum, y, b, y, dum, 110, info )
#
# will result in the computation of B and RSD, with RSD
# overwriting Y. More generally, each item in the following
# list contains groups of permissible identifications for
# a single calling sequence.
#
# 1. (Y,QTY,B) (RSD) (AB) (QY)
#
# 2. (Y,QTY,RSD) (B) (AB) (QY)
#
# 3. (Y,QTY,AB) (B) (RSD) (QY)
#
# 4. (Y,QY) (QTY,B) (RSD) (AB)
#
# 5. (Y,QY) (QTY,RSD) (B) (AB)
#
# 6. (Y,QY) (QTY,AB) (B) (RSD)
#
# In any group the value returned in the array allocated to
# the group corresponds to the last member of the group.
#
# 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), contains the output of DQRDC.
#
# Input, integer LDA, the leading dimension of the array A.
#
# Input, integer N, the number of rows of the matrix AK. It must
# have the same value as N in DQRDC.
#
# Input, integer K, the number of columns of the matrix AK. K
# must not be greater than min(N,P), where P is the same as in the
# calling sequence to DQRDC.
#
# Input, real QRAUX(P), the auxiliary output from DQRDC.
#
# Input, real Y(N), a vector to be manipulated by DQRSL.
#
# Input, integer JOB, specifies what is to be computed. JOB has
# the decimal expansion ABCDE, with the following meaning:
# if A /= 0, compute QY.
# if B /= 0, compute QTY.
# if C /= 0, compute QTY and B.
# if D /= 0, compute QTY and RSD.
# if E /= 0, compute QTY and AB.
# Note that a request to compute B, RSD, or AB automatically triggers
# the computation of QTY, for which an array must be provided in the
# calling sequence.
#
# Output, real QY(N), contains Q * Y, if requested.
#
# Output, real QTY(N), contains Q' * Y, if requested.
#
# Output, real B(K), the solution of the least squares problem
# minimize norm2 ( Y - AK * B),
# if its computation has been requested. Note that if pivoting was
# requested in DQRDC, the J-th component of B will be associated with
# column JPVT(J) of the original matrix A that was input into DQRDC.
#
# Output, real RSD(N), the least squares residual Y - AK * B,
# if its computation has been requested. RSD is also the orthogonal
# projection of Y onto the orthogonal complement of the column space
# of AK.
#
# Output, real AB(N), the least squares approximation Ak * B,
# if its computation has been requested. AB is also the orthogonal
# projection of Y onto the column space of A.
#
# Output, integer INFO, is zero unless the computation of B has
# been requested and R is exactly singular. In this case, INFO is the
# index of the first zero diagonal element of R, and B is left unaltered.
#
import numpy as np
qy = np.zeros ( n )
qty = np.zeros ( n )
b = np.zeros ( k )
rsd = np.zeros ( n )
ab = np.zeros ( n )
#
# Set info flag.
#
info = 0
#
# Determine what is to be computed.
#
cqy = int ( job / 10000 ) != 0
cqty = ( job % 10000 ) != 0
cb = int ( ( job % 1000 ) / 100 ) != 0
cr = int ( ( job % 100 ) / 10 ) != 0
cab = ( job % 10 ) != 0
ju = min ( k, n - 1 )
#
# Special action when N = 1.
#
if ( ju == 0 ):
qy[0] = y[0]
qty[0] = y[0]
ab[0] = y[0]
if ( a[0,0] == 0.0 ):
info = 1
else:
b[0] = y[0] / a[0]
rsd[0] = 0.0
return qy, qty, b, rsd, ab, info
#
# Set up to compute QY or QTY.
#
for i in range ( 0, n ):
qy[i] = y[i]
qty[i] = y[i]
#
# Compute QY.
#
if ( cqy ):
for j in range ( ju - 1, -1, -1 ):
if ( qraux[j] != 0.0 ):
ajj = a[j,j]
a[j,j] = qraux[j]
t = 0.0
for i in range ( j, n ):
t = t - a[i,j] * qy[i]
t = t / a[j,j]
for i in range ( j, n ):
qy[i] = qy[i] + t * a[i,j]
a[j,j] = ajj
#
# Compute Q'*Y.
#
if ( cqty ):
for j in range ( 0, ju ):
if ( qraux[j] != 0.0 ):
ajj = a[j,j]
a[j,j] = qraux[j]
t = 0.0
for i in range ( j, n ):
t = t - a[i,j] * qty[i]
t = t / a[j,j]
for i in range ( j, n ):
qty[i] = qty[i] + t * a[i,j]
a[j,j] = ajj
#
# Set up to compute B, RSD, or AB.
#
for i in range ( 0, k ):
b[i] = qty[i]
ab[i] = qty[i]
for i in range ( k, n ):
rsd[i] = qty[i]
#
# Compute B.
#
if ( cb ):
for j in range ( k - 1, -1, -1 ):
if ( a[j,j] == 0.0 ):
info = j + 1
break
b[j] = b[j] / a[j,j]
if ( 0 < j ):
t = - b[j]
for i in range ( 0, j ):
b[i] = b[i] + t * a[i,j]
#
# Compute RSD or AB as required.
#
if ( cr or cab ):
for j in range ( ju - 1, -1, -1 ):
if ( qraux[j] != 0.0 ):
ajj = a[j,j]
a[j,j] = qraux[j]
if ( cr ):
t = 0.0
for i in range ( j, n ):
t = t - a[i,j] * rsd[i]
t = t / a[j,j]
for i in range ( j, n ):
rsd[i] = rsd[i] + t * a[i,j]
if ( cab ):
t = 0.0
for i in range ( j, n ):
t = t - a[i,j] * ab[i]
t = t / a[j,j]
for i in range ( j, n ):
ab[i] = ab[i] + t * a[i,j]
a[j,j] = ajj
return qy, qty, b, rsd, ab, info
def dqrsl_test ( ):
#*****************************************************************************80
#
## DQRSL_TEST tests DQRSL.
#
# Discussion:
#
# DQRSL can solve a linear system that was factored by DQRDC.
#
# Licensing:
#
# This code is distributed under the GNU LGPL license.
#
# Modified:
#
# 30 August 2016
#
# Author:
#
# John Burkardt
#
import numpy as np
import platform
from dqrdc import dqrdc
n = 5
p = 3
lda = n
print ( '' )
print ( 'DQRSL_TEST' )
print ( ' Python version: %s' % ( platform.python_version ( ) ) )
print ( ' DQRSL solves a rectangular linear system A*x=b in the' )
print ( ' least squares sense after A has been factored by DQRDC.' )
#
# Set the matrix A.
#
a = np.array \
( [ \
[ 1.0, 1.0, 1.0 ], \
[ 1.0, 2.0, 4.0 ], \
[ 1.0, 3.0, 9.0 ], \
[ 1.0, 4.0, 16.0 ], \
[ 1.0, 5.0, 25.0 ] \
] )
print ( '' )
print ( ' The 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 ( n, dtype = np.int32 )
a, qraux, ipvt = dqrdc ( a, lda, n, p, ipvt, job )
#
# Call DQRSL to compute the least squares solution A*x=b.
#
job = 110
y = np.array \
( [ \
1.0, \
2.3, \
4.6, \
3.1, \
1.2 \
] )
b2 = np.array \
( [ \
-3.02, \
4.4914286, \
-0.72857143
] )
qy, qty, b, r, xb, info = dqrsl ( a, lda, n, p, qraux, y, job )
if ( info != 0 ):
print ( '' )
print ( 'DQRSL_TEST - Warning!' )
print ( ' DQRSL returns INFO = %d' % ( info ) )
return
print ( '' )
print ( ' X X(expected):' )
print ( '' )
for i in range ( 0, p ):
print ( ' %14.6g %14.6g' % ( b[i], b2[i] ) )
#
# Terminate.
#
print ( '' )
print ( 'DQRSL_TEST' )
print ( ' Normal end of execution.' )
return
if ( __name__ == '__main__' ):
from timestamp import timestamp
timestamp ( )
dqrsl_test ( )
timestamp ( )