#! /usr/bin/env python
#
def dqrlss ( a, lda, m, n, kr, b, jpvt, qraux ):
#*****************************************************************************80
#
## DQRLSS solves a linear system in a least squares sense.
#
# Discussion:
#
# DQRLSS must be preceded by a call to DQRANK.
#
# The system is to be solved is
# A * X = B
# where
# A is an M by N matrix with rank KR, as determined by DQRANK,
# B is a given M-vector,
# X is the N-vector to be computed.
#
# A solution X, with at most KR nonzero components, is found which
# minimizes the 2-norm of the residual (A*X-B).
#
# Licensing:
#
# This code is distributed under the GNU LGPL license.
#
# Modified:
#
# 31 August 2016
#
# Author:
#
# Original FORTRAN77 version by Dongarra, Moler, Bunch, Stewart.
# Python version by John Burkardt.
#
# Parameters:
#
# Input, real A(LDA,N), the QR factorization information
# from DQRANK. 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, integer KR, the rank of the matrix, as estimated
# by DQRANK.
#
# Input, real B(M), the right hand side of the linear system.
#
# Output, real X(N), a least squares solution to the
# linear system.
#
# Output, real R(M), the residual, B - A*X. RSD may
# overwite B.
#
# Input, integer JPVT(N), the pivot information from DQRANK.
# Columns JPVT(1), ..., JPVT(KR) of the original matrix are linearly
# independent to within the tolerance TOL and the remaining columns
# are linearly dependent.
#
# Input, real QRAUX(N), auxiliary information from DQRANK
# defining the QR factorization.
#
import numpy as np
from dqrsl import dqrsl
#
# Solve the reduced system of rank KR.
#
if ( 0 < kr ):
job = 110
qy, qty, x, r, ab, info = dqrsl ( a, lda, m, kr, qraux, b, job );
#
# Reverse the pivoting to recover the original variable ordering.
#
jpvt[0:n] = - jpvt[0:n]
for j in range ( kr, n ):
x = np.append ( x, 0.0 )
for j in range ( 0, n ):
if ( jpvt[j] <= 0 ):
k = - jpvt[j]
jpvt[j] = k
while ( k != j + 1 ):
t = x[j]
x[j] = x[k-1]
x[k-1] = t
jpvt[k-1] = - jpvt[k-1]
k = jpvt[k-1]
return x, r
def dqrlss_test ( ):
#*****************************************************************************80
#
## DQRLSS_TEST tests DQRLSS.
#
# Licensing:
#
# This code is distributed under the GNU LGPL license.
#
# Modified:
#
# 31 August 2016
#
# Author:
#
# John Burkardt
#
import numpy as np
import platform
from dqrank import dqrank
from r8vec_print import r8vec_print
print ( '' )
print ( 'DQRLSS_TEST' )
print ( ' Python version: %s' % ( platform.python_version ( ) ) )
print ( ' DQRLSS solves a rectangular linear system A*x = b' )
print ( ' with possibly non-full rank, in the least squares sense.' )
print ( ' Compare a tabulated solution X1 to the computed result X2.' )
#
# Set the matrix A.
#
m = 6;
n = 6;
a = np.zeros ( [ m, n ] )
b = np.zeros ( m )
x1 = np.zeros ( n )
for j in range ( 0, n ):
x1[j] = float ( j + 1 )
for i in range ( 0, m ):
a[i,j] = 1.0 / float ( i + j + 1 )
b[i] = b[i] + a[i,j] * x1[j]
print ( '' )
print ( ' The matrix A:' )
print ( '' )
for i in range ( 0, m ):
for j in range ( 0, n ):
print ( ' %14.6g' % ( a[i,j] ) ),
print ( '' )
r8vec_print ( n, x1, ' Solution x1:' )
r1 = b - np.dot ( a, x1 )
r8vec_print ( m, r1, ' Residual b-A*x1:' )
lda = m
tol = 0.000001
#
# Factor the matrix.
#
a_qr = a.copy ( )
kr, jpvt, qraux, a_qr = dqrank ( a_qr, lda, m, n, tol )
print ( '' )
print ( ' Matrix is of order m=%d by n=%d' % ( m, n ) )
print ( ' Condition number tolerance is %g' % ( tol ) )
print ( ' DQRANK estimates rank as %d' % ( kr ) )
#
# Solve the linear least squares problem.
#
x2, r = dqrlss ( a_qr, lda, m, n, kr, b, jpvt, qraux )
r8vec_print ( n, x2, ' Computed solution x2:' )
r2 = b - np.dot ( a, x2 )
r8vec_print ( m, r2, ' Residual b-A*x2:' )
#
# Terminate.
#
print ( '' )
print ( 'DQRLSS_TEST' )
print ( ' Normal end of execution.' )
return
if ( __name__ == '__main__' ):
from timestamp import timestamp
timestamp ( )
dqrlss_test ( )
timestamp ( )