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