#! /usr/bin/env python3 # 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 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 MIT 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. # # Input: # # real A(LDA,N), the matrix whose decomposition is to be computed. # # integer LDA, the leading dimension of A, which must # be at least M. # # integer M, the number of rows of A. # # integer N, the number of columns of A. # # 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. # # 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. # # real QRAUX(N), will contain extra information defining # the QR factorization. # # real A(LDA,N), 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. # import numpy as np 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(). # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 31 August 2016 # # Author: # # John Burkardt # import numpy as np import platform 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 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 MIT 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 # # Input: # # real A(LDA,P), the N by P matrix whose decomposition is to # be computed. # # integer LDA, the leading dimension of the array A. LDA must # be at least N. # # integer N, the number of rows of the matrix A. # # integer P, the number of columns of the matrix A. # # 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. # # 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. # # real QRAUX(P), contains further information required # to recover the orthogonal part of the decomposition. # # 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(). # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 30 August 2016 # # Author: # # John Burkardt # import numpy as np import platform 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 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 MIT license. # # Modified: # # 31 August 2016 # # Author: # # Original FORTRAN77 version by Dongarra, Moler, Bunch, Stewart. # Python version by John Burkardt. # # 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. # # integer LDA, the leading dimension of A, which must # be at least M. # # integer M, the number of rows of A. # # integer N, the number of columns of A. # # integer KR, the rank of the matrix, as estimated # by dqrank. # # real B(M), the right hand side of the linear system. # # 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. # # real QRAUX(N), auxiliary information from dqrank # defining the QR factorization. # # Output: # # real X(N), a least squares solution to the # linear system. # # real R(M), the residual, B - A*X. RSD may # overwite B. # import numpy as np # # 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 MIT license. # # Modified: # # 31 August 2016 # # Author: # # John Burkardt # import numpy as np import platform 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 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 MIT 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 # # Input: # # real A(LDA,P), contains the output of dqrdc. # # integer LDA, the leading dimension of the array A. # # integer N, the number of rows of the matrix AK. It must # have the same value as N in dqrdc. # # 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. # # real QRAUX(P), the auxiliary output from dqrdc. # # real Y(N), a vector to be manipulated by dqrsl. # # 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. # # real QTY(N), contains Q' * Y, if requested. # # 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. # # 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. # # 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. # # 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(). # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 30 August 2016 # # Author: # # John Burkardt # import numpy as np import platform 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 def drotg ( a, b ): #*****************************************************************************80 # ## drotg() constructs a Givens plane rotation. # # Discussion: # # Given values A and B, this routine computes # # SIGMA = sign ( A ) if abs ( A ) > abs ( B ) # = sign ( B ) if abs ( A ) <= abs ( B ) # # R = SIGMA * ( A * A + B * B ) # # C = A / R if R is not 0 # = 1 if R is 0 # # S = B / R if R is not 0, # 0 if R is 0. # # The computed numbers then satisfy the equation # # ( C S ) ( A ) = ( R ) # ( -S C ) ( B ) = ( 0 ) # # The routine also computes # # Z = S if abs ( A ) > abs ( B ), # = 1 / C if abs ( A ) <= abs ( B ) and C is not 0, # = 1 if C is 0. # # The single value Z encodes C and S, and hence the rotation: # # If Z = 1, set C = 0 and S = 1 # If abs ( Z ) < 1, set C = sqrt ( 1 - Z * Z ) and S = Z # if abs ( Z ) > 1, set C = 1/ Z and S = sqrt ( 1 - C * C ) # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 06 September 2016 # # Author: # # Python version by John Burkardt. # # Reference: # # Jack Dongarra, Cleve Moler, Jim Bunch and Pete Stewart, # LINPACK User's Guide, # SIAM, 1979. # # Charles Lawson, Richard Hanson, David Kincaid, Fred Krogh, # Basic Linear Algebra Subprograms for Fortran Usage, # Algorithm 539, # ACM Transactions on Mathematical Software, # Volume 5, Number 3, September 1979, pages 308-323. # # Input: # # real A, B, the values A and B. # # Output: # # real C, S, the cosine and sine of the Givens rotation. # # real R, Z, the values R and Z. # import numpy as np if ( abs ( b ) < abs ( a ) ): roe = a else: roe = b scale = abs ( a ) + abs ( b ) if ( scale == 0.0 ): c = 1.0 s = 0.0 r = 0.0 else: r = scale * np.sqrt ( ( a / scale ) ** 2 + ( b / scale ) ** 2 ) if ( roe < 0.0 ): r = - r c = a / r s = b / r if ( 0.0 < abs ( c ) and abs ( c ) <= s ): z = 1.0 / c else: z = s return c, s, r, z def drotg_test ( ): #*****************************************************************************80 # ## drotg_test() tests drotg(). # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 06 September 2016 # # Author: # # John Burkardt # import numpy as np import platform test_num = 5 print ( '' ) print ( 'drotg_test' ) print ( ' Python version: %s' % ( platform.python_version ( ) ) ) print ( ' drotg generates a real Givens rotation' ) print ( ' ( C S ) * ( A ) = ( R )' ) print ( ' ( -S C ) ( B ) ( 0 )' ) print ( '' ) for test in range ( 0, test_num ): a = np.random.rand ( ) b = np.random.rand ( ) c, s, r, z = drotg ( a, b ) print ( '' ) print ( ' A = %g B = %g' % ( a, b ) ) print ( ' C = %g S = %g' % ( c, s ) ) print ( ' R = %g Z = %g' % ( r, z ) ) print ( ' C*A+S*B = %g' % ( c * a + s * b ) ) print ( ' -S*A+C*B = %g' % ( -s * a + c * b ) ) # # Terminate. # print ( '' ) print ( 'drotg_test' ) print ( ' Normal end of execution.' ) return def drot ( n, x, incx, y, incy, c, s ): #*****************************************************************************80 # ## drot() applies a plane rotation. # # Discussion: # # Note that the FORTRAN version of this function allowed users to pass in # X and Y data that was noncontiguous, (such as rows of a FORTRAN matrix). # The rotated data overwrote the input data, and so it might therefore # also be noncontiguous. # # This function does not assume that the output overwrites the input, # and treats the output vectors as new items of length exactly N. # # Note, moreover, that Python does NOT allow a matrix to be treated as a # vector in quite the simple way that FORTRAN does. # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 06 September 2016 # # Author: # # Python version by John Burkardt. # # Reference: # # Jack Dongarra, Cleve Moler, Jim Bunch and Pete Stewart, # LINPACK User's Guide, # SIAM, 1979. # # Charles Lawson, Richard Hanson, David Kincaid, Fred Krogh, # Basic Linear Algebra Subprograms for Fortran Usage, # Algorithm 539, # ACM Transactions on Mathematical Software, # Volume 5, Number 3, September 1979, pages 308-323. # # Input: # # integer N, the number of entries in the vectors. # # real X(INCX*N), one of the vectors to be rotated. # # integer INCX, the increment between successive entries of X. # # real Y(INCX*N), one of the vectors to be rotated. # # integer INCY, the increment between successive elements of Y. # # real C, S, parameters (presumably the cosine and # sine of some angle) that define a plane rotation. # # Output: # # real XR(N), the rotated vector. # # real YR(N), the rotated vector. # import numpy as np xr = np.zeros ( n ) yr = np.zeros ( n ) if ( n <= 0 ): pass elif ( incx == 1 and incy == 1 ): xr[0:n] = c * x[0:n] + s * y[0:n] yr[0:n] = c * y[0:n] - s * x[0:n] else: if ( 0 <= incx ): ix = 0 else: ix = ( - n + 1 ) * incx if ( 0 <= incy ): iy = 0 else: iy = ( - n + 1 ) * incy for i in range ( 0, n ): xr[ix] = c * x[ix] + s * y[iy] yr[iy] = c * y[iy] - s * x[ix] ix = ix + incx iy = iy + incy return xr, yr def drot_test ( ): #*****************************************************************************80 # ## drot_test() tests drot(). # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 06 September 2016 # # Author: # # John Burkardt # import numpy as np import platform n = 6 x = np.zeros ( n ) y = np.zeros ( n ) for i in range ( 0, n ): x[i] = float ( i + 1 ) y[i] = x[i] * x[i] - 12.0 print ( '' ) print ( 'drot_test' ) print ( ' Python version: %s' % ( platform.python_version ( ) ) ) print ( ' drot carries out a Givens rotation.' ) print ( '' ) print ( ' Vectors X and Y' ) print ( '' ) for i in range ( 0, n ): print ( ' %6d %14.6g %14.6g' % ( i, x[i], y[i] ) ) c = 0.5 s = np.sqrt ( 1.0 - c * c ) xr, yr = drot ( n, x, 1, y, 1, c, s ) print ( '' ) print ( ' xr, yr = drot ( n, x, 1, y, 1, %g, %g )' % ( c, s ) ) print ( '' ) print ( ' Rotated vectors XR and YR' ) print ( '' ) for i in range ( 0, n ): print ( ' %6d %14.6g %14.6g' % ( i, xr[i], yr[i] ) ) t = np.arctan2 ( y[0], x[0] ) c = np.cos ( t ) s = np.sin ( t ) xr, yr = drot ( n, x, 1, y, 1, c, s ) print ( '' ) print ( ' xr, yr = drot ( n, x, 1, y, 1, %g, %g )' % ( c, s ) ) print ( '' ) print ( ' Rotated vectors XR and YR' ) print ( '' ) for i in range ( 0, n ): print ( ' %6d %14.6g %14.6g' % ( i, xr[i], yr[i] ) ) # # Terminate. # print ( '' ) print ( 'drot_test' ) print ( ' Normal end of execution.' ) return def dsvdc ( a, lda, m, n, ldu, ldv, job ): #*****************************************************************************80 # ## dsvdc() computes the singular value decomposition of a real rectangular matrix. # # Discussion: # # This routine reduces an M by N matrix A to diagonal form by orthogonal # transformations U and V. The diagonal elements S(I) are the singular # values of A. The columns of U are the corresponding left singular # vectors, and the columns of V the right singular vectors. # # The form of the singular value decomposition is then # # A(MxN) = U(MxM) * S(MxN) * V(NxN)' # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 01 September 2016 # # Author: # # Original FORTRAN77 version by Dongarra, Moler, Bunch, Stewart. # Python version by John Burkardt. # # Reference: # # Dongarra, Moler, Bunch and 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 # # Input: # # real A(LDA,N), the M by N matrix whose singular # value decomposition is to be computed. # # integer LDA, the leading dimension of the array A. # LDA must be at least N. # # integer M, the number of rows of the matrix. # # integer N, the number of columns of the matrix A. # # integer LDU, the leading dimension of the array U. # LDU must be at least M. # # integer LDV, the leading dimension of the array V. # LDV must be at least N. # # integer JOB, controls the computation of the singular # vectors. It has the decimal expansion AB with the following meaning: # A = 0, for not compute the left singular vectors. # A = 1, return the M left singular vectors in U. # A >= 2, return the first min(M,N) singular vectors in U. # B = 0, for not compute the right singular vectors. # B = 1, return the right singular vectors in V. # # Output: # # real A(LDA,N), the matrix has been destroyed. Depending on the # user's requests, the matrix may contain other useful information. # # real S(MM), where MM = max(M+1,N). The first # min(M,N) entries of S contain the singular values of A arranged in # descending order of magnitude. # # real E(MM), where MM = max(M+1,N), ordinarily contains zeros. # However see the discussion of INFO for exceptions. # # real U(LDU,K). If JOBA = 1 then K = M # if 2 <= JOBA, then K = min(M,N). U contains the M by M matrix of # left singular vectors. U is not referenced if JOBA = 0. If M <= N # or if JOBA = 2, then U may be identified with A in the subroutine call. # # real V(LDV,N), the N by N matrix of right singular # vectors. V is not referenced if JOB is 0. If N <= M, then V may be # identified with A in the subroutine call. # # integer INFO, status indicator. # The singular values (and their corresponding singular vectors) # S(INFO+1), S(INFO+2),...,S(MN) are correct. Here MN = min ( M, N ). # Thus if INFO is 0, all the singular values and their vectors are # correct. In any event, the matrix B = U' * A * V is the bidiagonal # matrix with the elements of S on its diagonal and the elements of E on # its superdiagonal. Thus the singular values of A and B are the same. # import numpy as np maxit = 30 mm = max ( m + 1, n ) s = np.zeros ( mm ) e = np.zeros ( mm ) joba = ( job // 10 ) if ( joba == 1 ): u = np.zeros ( [ ldu, m ] ) elif ( 1 <= joba ): u = np.zeros ( [ ldu, min ( m, n ) ] ) v = np.zeros ( [ ldv, n ] ) work = np.zeros ( m ) info = 0 # # Determine what is to be computed. # wantu = 0 wantv = 0 jobu = int ( ( job % 100 ) / 10 ) if ( 1 < jobu ): ncu = min ( m, n ) else: ncu = m if ( jobu != 0 ): wantu = 1 if ( ( job % 10 ) != 0 ): wantv = 1 # # Reduce A to bidiagonal form, storing the diagonal elements # in S and the super-diagonal elements in E. # info = 0 nct = min ( m - 1, n ) nrt = max ( 0, min ( m, n - 2 ) ) lu = max ( nct, nrt ) for l in range ( 0, lu ): # # Compute the transformation for the L-th column and # place the L-th diagonal in S(L). # if ( l <= nct - 1 ): t = 0.0 for i in range ( l, m ): t = t + a[i,l] ** 2 t = np.sqrt ( t ) s[l] = t if ( s[l] != 0.0 ): if ( a[l,l] < 0.0 ): s[l] = - s[l] for i in range ( l, m ): a[i,l] = a[i,l] / s[l] a[l,l] = 1.0 + a[l,l] s[l] = - s[l] for j in range ( l + 1, n ): # # Apply the transformation. # if ( l <= nct - 1 and s[l] != 0.0 ): t = 0.0 for i in range ( l, m ): t = t - a[i,l] * a[i,j] t = t / a[l,l] for i in range ( l, m ): a[i,j] = a[i,j] + t * a[i,l] # # Place the L-th row of A into E for the # subsequent calculation of the row transformation. # e[j] = a[l,j] # # Place the transformation in U for subsequent back multiplication. # if ( wantu and l <= nct - 1 ): for i in range ( l, m ): u[i,l] = a[i,l] # # Compute the L-th row transformation and place the # L-th superdiagonal in E(L). # if ( l <= nrt - 1 ): t = 0.0 for i in range ( l + 1, n ): t = t + e[i] ** 2 e[l] = np.sqrt ( t ) if ( e[l] != 0.0 ): if ( e[l+1] < 0.0 ): e[l] = - e[l] e[l+1:n] = e[l+1:n] / e[l] e[l+1] = 1.0 + e[l+1] e[l] = - e[l] # # Apply the transformation. # if ( l + 1 <= m - 1 and e[l] != 0.0 ): work[l+1:m] = 0.0 for j in range ( l + 1, n ): for i in range ( l + 1, m ): work[i] = work[i] + e[j] * a[i,j] for j in range ( l + 1, n ): for i in range ( l + 1, m ): a[i,j] = a[i,j] - e[j] / e[l+1] * work[i] # # Place the transformation in V for subsequent back multiplication. # if ( wantv ): v[l+1:n,l] = e[l+1:n] # # Set up the final bidiagonal matrix of order MN. # mn = min ( m + 1, n ) nctp1 = nct + 1 nrtp1 = nrt + 1 if ( nct < n ): s[nctp1-1] = a[nctp1-1,nctp1-1] if ( m < mn ): s[mn-1] = 0.0 if ( nrtp1 < mn ): e[nrtp1-1] = a[nrtp1-1,mn-1] e[mn-1] = 0.0 # # If required, generate U. # if ( wantu ): for j in range ( nctp1 - 1, ncu ): for i in range ( 0, m ): u[i,j] = 0.0 for j in range ( nctp1 - 1, ncu ): u[j,j] = 1.0 for l in range ( nct - 1, -1, -1 ): if ( s[l] != 0.0 ): for j in range ( l + 1, ncu ): t = 0.0 for i in range ( l, m ): t = t - u[i,l] * u[i,j] t = t / u[l,l] for i in range ( l, m ): u[i,j] = u[i,j] + t * u[i,l] for i in range ( 0, l ): u[i,l] = 0.0 for i in range ( l, m ): u[i,l] = - u[i,l] u[l,l] = 1.0 + u[l,l] else: for i in range ( 0, m ): u[i,l] = 0.0 u[l,l] = 1.0 # # If it is required, generate V. # if ( wantv ): for l in range ( n - 1, -1, -1 ): if ( l <= nrt - 1 and e[l] != 0.0 ): for j in range ( l + 1, n ): t = 0.0 for i in range ( l + 1, n ): t = t - v[i,l] * v[i,j] t = t / v[l+1,l] for i in range ( l + 1, n ): v[i,j] = v[i,j] + t * v[i,l] for i in range ( 0, n ): v[i,l] = 0.0 v[l,l] = 1.0 # # Main iteration loop for the singular values. # mm = mn iter = 0 while ( 0 < mn ): # # If too many iterations have been performed, set flag and return. # if ( maxit <= iter ): print ( '' ) print ( 'dsvdc - Warning!' ) print ( ' MAXIT = %d <= ITER = %d' % ( maxit, iter ) ) info = mn return a, s, e, u, v, info # # This section of the program inspects for # negligible elements in the S and E arrays. # # On completion the variables KASE and L are set as follows: # # KASE = 1 if S(MN) and E(L-1) are negligible and L < MN # KASE = 2 if S(L) is negligible and L < MN # KASE = 3 if E(L-1) is negligible, L < MN, and # S(L), ..., S(MN) are not negligible (QR step). # KASE = 4 if E(MN-1) is negligible (convergence). # for l in range ( mn - 2, -2, -1 ): if ( l == -1 ): break test = abs ( s[l] ) + abs ( s[l+1] ) ztest = test + abs ( e[l] ) if ( ztest == test ): e[l] = 0.0 break if ( l == mn - 2 ): kase = 4 else: for ls in range ( mn - 1, l - 1, -1 ): if ( ls == l ): break test = 0.0 if ( ls != mn - 1 ): test = test + abs ( e[ls] ) if ( ls != l + 1 ): test = test + abs ( e[ls-1] ) ztest = test + abs ( s[ls] ) if ( ztest == test ): s[ls] = 0.0 break if ( ls == l ): kase = 3 elif ( ls == mn - 1 ): kase = 1 else: kase = 2 l = ls l = l + 1 # # Deflate negligible S(MN). # if ( kase == 1 ): mm1 = mn - 1 f = e[mn-2] e[mn-2] = 0.0 for k in range ( mm1 - 1, l - 1, -1 ): t1 = s[k] cs, sn, t1, f = drotg ( t1, f ) s[k] = t1 if ( k != l ): f = - sn * e[k-1] e[k-1] = cs * e[k-1] if ( wantv ): v1, v2 = drot ( n, v[0:n,k], 1, v[0:n,mn-1], 1, cs, sn ) for i in range ( 0, n ): v[i,k] = v1[i] v[i,mn-1] = v2[i] # # Split at negligible S(L). # elif ( kase == 2 ): f = e[l-1] e[l-1] = 0.0 for k in range ( l, mn ): t1 = s[k] cs, sn, t1, f = drotg ( t1, f ) s[k] = t1 f = - sn * e[k] e[k] = cs * e[k] if ( wantu ): u1, u2 = drot ( m, u[0:m,k], 1, u[0:m,l-1], 1, cs, sn ) for i in range ( 0, m ): u[i,k] = u1[i] u[i,l-1] = u2[i] # # Perform one QR step. # elif ( kase == 3 ): # # Calculate the shift. # scale = max ( abs ( s[mn-1] ), \ max ( abs ( s[mn-2] ), \ max ( abs ( e[mn-2] ), \ max ( abs ( s[l] ), \ abs ( e[l] ) ) ) ) ) sm = s[mn-1] / scale smm1 = s[mn-2] / scale emm1 = e[mn-2] / scale sl = s[l] / scale el = e[l] / scale b = ( ( smm1 + sm ) * ( smm1 - sm ) + emm1 * emm1 ) / 2.0 c = sm * sm * emm1 * emm1 shift = 0.0 if ( b != 0.0 or c != 0.0 ): shift = np.sqrt ( b * b + c ) if ( b < 0.0 ): shift = - shift shift = c / ( b + shift ) f = ( sl + sm ) * ( sl - sm ) - shift g = sl * el # # Chase zeros. # mm1 = mn - 1 for k in range ( l, mm1 ): cs, sn, f, g = drotg ( f, g ) if ( k != l ): e[k-1] = f f = cs * s[k] + sn * e[k] e[k] = cs * e[k] - sn * s[k] g = sn * s[k+1] s[k+1] = cs * s[k+1] if ( wantv ): v1, v2 = drot ( n, v[0:n,k], 1, v[0:n,k+1], 1, cs, sn ) for i in range ( 0, n ): v[i,k] = v1[i] v[i,k+1] = v2[i] cs, sn, f, g = drotg ( f, g ) s[k] = f f = cs * e[k] + sn * s[k+1] s[k+1] = - sn * e[k] + cs * s[k+1] g = sn * e[k+1] e[k+1] = cs * e[k+1] if ( wantu and k < m ): u1, u2 = drot ( m, u[0:m,k], 1, u[0:m,k+1], 1, cs, sn ) for i in range ( 0, m ): u[i,k] = u1[i] u[i,k+1] = u2[i] e[mn-2] = f iter = iter + 1 # # Convergence. # elif ( kase == 4 ): # # Make the singular value nonnegative. # if ( s[l] < 0.0 ): s[l] = - s[l] if ( wantv ): for i in range ( 0, n ): v[i,l] = - v[i,l] # # Order the singular values. # while ( True ): if ( l + 1 == mm ): break if ( s[l+1] <= s[l] ): break t = s[l] s[l] = s[l+1] s[l+1] = t if ( wantv and l < n - 1 ): for i in range ( 0, n ): t = v[i,l] v[i,l] = v[i,l+1] v[i,l+1] = t if ( wantu and l < m - 1 ): for i in range ( 0, m ): t = u[i,l] u[i,l] = u[i,l+1] u[i,l+1] = t l = l + 1 iter = 0 mn = mn - 1 return a, s, e, u, v, info def dsvdc_test ( ): #*****************************************************************************80 # ## dsvdc_test() tests dsvdc(). # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 07 September 2016 # # Author: # # John Burkardt # import numpy as np import platform m = 6 n = 4 print ( '' ) print ( 'dsvdc_test' ) print ( ' Python version: %s' % ( platform.python_version ( ) ) ) print ( ' dsvdc computes the singular value decomposition' ) print ( ' for an MxN matrix A in general storage.' ) print ( ' A = U * S * V\'' ) print ( '' ) print ( ' Matrix rows M = %d' % ( m ) ) print ( ' Matrix columns N = %d' % ( n ) ) # # Set A. # a = np.random.rand ( m, n ) r8mat_print ( m, n, a, ' The matrix A:' ) # # Decompose the matrix. # print ( '' ) print ( ' Decompose the matrix.' ) job = 11 lda = m ldu = m ldv = n a, s, e, u, v, info = dsvdc ( a, lda, m, n, ldu, ldv, job ) if ( info != 0 ): print ( '' ) print ( 'dsvdc_test - Warning:' ) print ( ' dsvdc returned nonzero INFO = %d' % ( info ) ) # return r8vec_print ( min ( m, n ), s, ' Singular values S:' ) r8mat_print ( m, m, u, ' Left Singular Vector Matrix U:' ) r8mat_print ( n, n, v, ' Right Singular Vector Matrix V:' ) sigma = np.zeros ( [ m, n ] ) for i in range ( 0, min ( m, n ) ): sigma[i,i] = s[i] usv = np.dot ( u, np.dot ( sigma, v.transpose ( ) ) ) r8mat_print ( m, n, usv, ' The product U * S * V\' (should equal A):' ) # # Terminate. # print ( '' ) print ( 'dsvdc_test' ) print ( ' Normal end of execution.' ) return def lstsq_solve_test ( ): #*****************************************************************************80 # ## lstsq_solve_test() tests the numpy.linalg.lstsq() operator. # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 29 August 2016 # # Author: # # John Burkardt # import numpy as np import platform print ( '' ) print ( 'lstsq_solve_test' ) print ( ' Python version: %s' % ( platform.python_version ( ) ) ) print ( ' The function x=np.linalg.lstsq(A,b)\n' ) print ( ' solves a linear system A*x = b in the least squares sense.' ) print ( ' Compare a tabulated solution X1 to the LSTSQ result X2.' ) prob_num = p00_prob_num ( ) print ( '' ) print ( ' Number of problems = %d' % ( prob_num ) ) print ( '' ) print ( ' Index M N ||B|| ||X1 - X2||' ), print ( ' ||X1|| ||X2|| ||R1|| ||R2||' ) print ( '' ) for prob in range ( 1, prob_num + 1 ): # # Get problem size. # m = p00_m ( prob ) n = p00_n ( prob ) # # Retrieve problem data. # a = p00_a ( prob, m, n ) b = p00_b ( prob, m ) x1 = p00_x ( prob, n ) b_norm = np.linalg.norm ( b ) x1_norm = np.linalg.norm ( x1 ) r1 = np.dot ( a, x1 ) - b r1_norm = np.linalg.norm ( r1 ) # # Use NUMPY's LINALG.LSTSQ on the problem. # x2, r, rank, s = np.linalg.lstsq ( a, b, rcond = None ) x2_norm = np.linalg.norm ( x2 ) r2 = np.dot ( a, x2 ) - b r2_norm = np.linalg.norm ( r2 ) # # Compare tabulated and computed solutions. # x_diff_norm = np.linalg.norm ( x1 - x2 ) # # Report results for this problem. # print ( ' %5d %4d %4d' % ( prob, m, n ), end = '' ) print ( ' %12.4g %12.4g' % ( b_norm, x_diff_norm ), end = '' ) print ( ' %12.4g %12.4g' % ( x1_norm, x2_norm ), end = '' ) print ( ' %12.4g %12.4g' % ( r1_norm, r2_norm ), end = '' ) print ( '' ) # # Terminate. # print ( '' ) print ( 'lstsq_solve_test' ) print ( ' Normal end of execution.' ) return def normal_solve ( m, n, a, b ): #*****************************************************************************80 # ## normal_solve() solves a linear system using the normal equations. # # Discussion: # # Given a presumably rectangular MxN system of equations A*x=b, this routine # sets up the NxN system A'*A*x=A'b. Assuming N <= M, and that A has full # column rank, the system will be solvable, and the vector x that is returned # will minimize the Euclidean norm of the residual. # # One drawback to this approach is that the condition number of the linear # system A'*A is effectively the square of the condition number of A, # meaning that there is a substantial loss of accuracy. # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 29 August 2016 # # Author: # # John Burkardt # # Reference: # # David Kahaner, Cleve Moler, Steven Nash, # Numerical Methods and Software, # Prentice Hall, 1989, # ISBN: 0-13-627258-4, # LC: TA345.K34. # # Input: # # integer M, the number of rows of A. # # integer N, the number of columns of A. # It must be the case that N <= M. # # real A(M,N), the matrix. # The matrix must have full column rank. # # real B(M), the right hand side. # # Output: # # real X(N), the least squares solution. # # integer FLAG, # 0, no error was detected. # 1, an error occurred. # import numpy as np flag = False x = np.zeros ( n ) if ( m < n ): flag = True else: ata = np.dot ( a.transpose ( ), a ) atb = np.dot ( a.transpose ( ), b ) ata_c, flag = r8mat_cholesky_factor ( n, ata ) if ( not flag ): x = r8mat_cholesky_solve ( n, ata_c, atb ) return x, flag def normal_solve_test ( ): #*****************************************************************************80 # ## normal_solve_test() tests normal_solve(). # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 29 August 2016 # # Author: # # John Burkardt # import numpy as np import platform print ( '' ) print ( 'normal_solve_test' ) print ( ' Python version: %s' % ( platform.python_version ( ) ) ) print ( ' normal_solve is a function with a simple interface which' ) print ( ' solves a linear system A*x = b in the least squares sense.' ) print ( ' Compare a tabulated solution X1 to the normal_solve result X2.' ) print ( '' ) print ( ' normal_solve cannot be applied when N < M,' ) print ( ' or if the matrix does not have full column rank.' ) prob_num = p00_prob_num ( ) print ( '' ) print ( ' Number of problems = %d' % ( prob_num ) ) print ( '' ) print ( ' Index M N ||B|| ||X1 - X2||' ), print ( ' ||X1|| ||X2|| ||R1|| ||R2||' ) print ( '' ) for prob in range ( 1, prob_num + 1 ): # # Get problem size. # m = p00_m ( prob ) n = p00_n ( prob ) # # Retrieve problem data. # a = p00_a ( prob, m, n ) b = p00_b ( prob, m ) x1 = p00_x ( prob, n ) b_norm = np.linalg.norm ( b ) x1_norm = np.linalg.norm ( x1 ) r1 = np.dot ( a, x1 ) - b r1_norm = np.linalg.norm ( r1 ) # # Use normal_solve on the problem. # x2, flag = normal_solve ( m, n, a, b ) if ( flag ): print ( ' %5d %4d %4d' % ( prob, m, n ) ), print ( ' %12.4g ------------' % ( b_norm ) ), print ( ' %12.4g ------------' % ( x1_norm ) ), print ( ' %12.4g ------------' % ( r1_norm ) ) else: x2_norm = np.linalg.norm ( x2 ) r2 = np.dot ( a, x2 ) - b r2_norm = np.linalg.norm ( r2 ) # # Compare tabulated and computed solutions. # x_diff_norm = np.linalg.norm ( x1 - x2 ) # # Report results for this problem. # print ( ' %5d %4d %4d' % ( prob, m, n ) ), print ( ' %12.4g %12.4g' % ( b_norm, x_diff_norm ) ), print ( ' %12.4g %12.4g' % ( x1_norm, x2_norm ) ), print ( ' %12.4g %12.4g' % ( r1_norm, r2_norm ) ) # # Terminate. # print ( '' ) print ( 'normal_solve_test' ) print ( ' Normal end of execution.' ) return def qr_solve ( m, n, a, b ): #*****************************************************************************80 # ## qr_solve() solves a linear system in the least squares sense. # # Discussion: # # If the matrix A has full column rank, then the solution X should be the # unique vector that minimizes the Euclidean norm of the residual. # # If the matrix A does not have full column rank, then the solution is # not unique; the vector X will minimize the residual norm, but so will # various other vectors. # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 27 August 2016 # # Author: # # John Burkardt # # Reference: # # David Kahaner, Cleve Moler, Steven Nash, # Numerical Methods and Software, # Prentice Hall, 1989, # ISBN: 0-13-627258-4, # LC: TA345.K34. # # Input: # # integer M, the number of rows of A. # # integer N, the number of columns of A. # # real A(M,N), the matrix. # # real B(M), the right hand side. # # Output: # # real X(N), the least squares solution. # lda = m eps = 2.220446049250313E-016 tol = 0.0 for i in range ( 0, m ): for j in range ( 0, n ): tol = max ( tol, abs ( a[i,j] ) ) tol = eps / tol # # Factor the matrix. # kr, jpvt, qraux, a = dqrank ( a, lda, m, n, tol ) # # Solve the least-squares problem. # x, r = dqrlss ( a, lda, m, n, kr, b, jpvt, qraux ) return x def qr_solve_test ( ): #*****************************************************************************80 # ## qr_solve_test() tests qr_solve(). # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 26 August 2016 # # Author: # # John Burkardt # import numpy as np import platform print ( '' ) print ( 'qr_solve_test' ) print ( ' Python version: %s' % ( platform.python_version ( ) ) ) print ( ' qr_solve is a function with a simple interface which' ) print ( ' solves a linear system A*x = b in the least squares sense.' ) print ( ' Compare a tabulated solution X1 to the qr_solve result X2.' ) prob_num = p00_prob_num ( ) print ( '' ) print ( ' Number of problems = %d' % ( prob_num ) ) print ( '' ) print ( ' Index M N ||B|| ||X1 - X2||' ), print ( ' ||X1|| ||X2|| ||R1|| ||R2||' ) print ( '' ) for prob in range ( 1, prob_num + 1 ): # # Get problem size. # m = p00_m ( prob ) n = p00_n ( prob ) # # Retrieve problem data. # a = p00_a ( prob, m, n ) b = p00_b ( prob, m ) x1 = p00_x ( prob, n ) b_norm = np.linalg.norm ( b ) x1_norm = np.linalg.norm ( x1 ) r1 = np.dot ( a, x1 ) - b r1_norm = np.linalg.norm ( r1 ) # # Use qr_solve on the problem. # Since we want to compute residuals ourselves, we need # to keep the originals of A and B, so we make copies # to send to qr_solve. # a_qr = a.copy ( ) b_qr = b.copy ( ) x2 = qr_solve ( m, n, a_qr, b_qr ) x2_norm = np.linalg.norm ( x2 ) r2 = np.dot ( a, x2 ) - b r2_norm = np.linalg.norm ( r2 ) # # Compare tabulated and computed solutions. # x_diff_norm = np.linalg.norm ( x1 - x2 ) # # Report results for this problem. # print ( ' %5d %4d %4d' % ( prob, m, n ) ), print ( ' %12.4g %12.4g' % ( b_norm, x_diff_norm ) ), print ( ' %12.4g %12.4g' % ( x1_norm, x2_norm ) ), print ( ' %12.4g %12.4g' % ( r1_norm, r2_norm ) ) # # Terminate. # print ( '' ) print ( 'qr_solve_test' ) print ( ' Normal end of execution.' ) return def qr_solve_tests ( ): #*****************************************************************************80 # ## qr_solve_tests() tests qr_solve(). # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 09 September 2016 # # Author: # # John Burkardt # import platform print ( '' ) print ( 'qr_solve_testS' ) print ( ' Python version: %s' % ( platform.python_version ( ) ) ) print ( ' Test the qr_solve library.' ) lstsq_solve_test ( ) normal_solve_test ( ) qr_solve_test ( ) svd_solve_test ( ) # # Terminate. # print ( '' ) print ( 'qr_solve_tests' ) print ( ' Normal end of execution.' ) return def r8mat_cholesky_factor ( n, a ): #*****************************************************************************80 # ## r8mat_cholesky_factor() computes the Cholesky factor of a symmetric matrix. # # Discussion: # # The matrix must be symmetric and positive semidefinite. # # For a positive semidefinite symmetric matrix A, the Cholesky factorization # is a lower triangular matrix L such that: # # A = L * L' # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 04 October 2012 # # Author: # # John Burkardt # # Input: # # integer N, the number of rows and columns of the matrix A. # # real A(N,N), the matrix. # # Output: # # real C(N,N), the N by N lower triangular Cholesky factor. # # bool FLAG: # False, no error occurred. # True, the matrix is not positive definite. # import numpy as np flag = False c = np.zeros ( [ n, n ] ) for j in range ( 0, n ): for i in range ( 0, n ): c[i,j] = a[i,j] for j in range ( 0, n ): c[0:j,j] = 0.0 for i in range ( j, n ): sum2 = c[j,i] for k in range ( 0, j ): sum2 = sum2 - c[j,k] * c[i,k] if ( i == j ): if ( sum2 <= 0.0 ): flag = True return c, flag else: c[i,j] = np.sqrt ( sum2 ) else: if ( c[j,j] != 0.0 ): c[i,j] = sum2 / c[j,j] else: c[i,j] = 0.0 return c, flag def r8mat_cholesky_factor_test ( ): #*****************************************************************************80 # ## r8mat_cholesky_factor_test() tests r8mat_cholesky_factor(). # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 28 August 2016 # # Author: # # John Burkardt # import numpy as np import platform n = 5 print ( '' ) print ( 'r8mat_cholesky_factor_test' ) print ( ' Python version: %s' % ( platform.python_version ( ) ) ) print ( ' r8mat_cholesky_factor determines the' ) print ( ' lower triangular Cholesky factorization' ) print ( ' of a positive definite symmetric matrix,' ) a = np.zeros ( [ n, n ] ) for i in range ( 0, n ): for j in range ( 0, n ): if ( i == j ): a[i,j] = 2.0 elif ( j == i - 1 or j == i + 1 ): a[i,j] = -1.0 r8mat_print ( n, n, a, ' Matrix to be factored:' ) # # Compute a Cholesky factor. # l, flag = r8mat_cholesky_factor ( n, a ) r8mat_print ( n, n, l, ' Cholesky factor L:' ) d = np.dot ( l, l.transpose ( ) ) r8mat_print ( n, n, d, ' Product L * L\':' ) # # Terminate. # print ( '' ) print ( 'r8mat_cholesky_factor_test:' ) print ( ' Normal end of execution.' ) return def r8mat_cholesky_solve ( n, l, b ): #*****************************************************************************80 # ## r8mat_cholesky_solve() solves a Cholesky factored linear system A * x = b. # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 28 August 2016 # # Author: # # John Burkardt # # Input: # # integer N, the number of rows and columns of the matrix A. # # real L(N,N), the N by N Cholesky factor of the # system matrix A. # # real B(N), the right hand side of the linear system. # # Output: # # real X(N), the solution of the linear system. # # # Solve L * y = b. # x = r8mat_l_solve ( n, l, b ) # # Solve L' * x = y. # x = r8mat_lt_solve ( n, l, x ) return x def r8mat_cholesky_solve_test ( ): #*****************************************************************************80 # ## r8mat_cholesky_solve_test() tests r8mat_cholesky_solve(). # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 28 August 2016 # # Author: # # John Burkardt # import numpy as np import platform n = 5 print ( '' ) print ( 'r8mat_cholesky_solve_test' ) print ( ' Python version: %s' % ( platform.python_version ( ) ) ) print ( ' r8mat_cholesky_solve solves a linear system' ) print ( ' using the lower triangular Cholesky factorization,' ) print ( ' for a positive definite symmetric matrix.' ) a = np.zeros ( [ n, n ] ) for i in range ( 0, n ): for j in range ( 0, n ): if ( i == j ): a[i,j] = 2.0 elif ( j == i - 1 or j == i + 1 ): a[i,j] = -1.0 r8mat_print ( n, n, a, ' Matrix to be factored:' ) # # Compute a Cholesky factor. # l, flag = r8mat_cholesky_factor ( n, a ) r8mat_print ( n, n, l, ' Cholesky factor L:' ) d = np.dot ( l, l.transpose ( ) ) r8mat_print ( n, n, d, ' Product L * L\':' ) # # Solve a linear system. # b = np.zeros ( n ) b[n-1] = float ( n + 1 ) r8vec_print ( n, b, ' Right hand side b:' ) x = r8mat_cholesky_solve ( n, l, b ) r8vec_print ( n, x, ' Computed solution x:' ) # # Terminate. # print ( '' ) print ( 'r8mat_cholesky_solve_test:' ) print ( ' Normal end of execution.' ) return def r8mat_l_solve ( n, a, b ): #*****************************************************************************80 # ## r8mat_l_solve() solves a lower triangular linear system. # # Discussion: # # An R8MAT is an MxN array of R8's, stored by (I,J) -> [I+J*M]. # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 28 August 2016 # # Author: # # John Burkardt # # Input: # # integer N, the number of rows and columns of # the matrix A. # # real A(N,N), the N by N lower triangular matrix. # # real B(N), the right hand side of the linear system. # # Output: # # real X(N), the solution of the linear system. # import numpy as np x = np.zeros ( n ) # # Solve L * x = b. # for i in range ( 0, n ): x[i] = b[i] for j in range ( 0, i ): x[i] = x[i] - a[i,j] * x[j] x[i] = x[i] / a[i,i] return x def r8mat_l_solve_test ( ): #*****************************************************************************80 # ## r8mat_l_solve_test() tests r8mat_l_solve(). # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 28 August 2015 # # Author: # # John Burkardt # import numpy as np import platform n = 4 a = np.array ( [ \ [ 1.0, 0.0, 0.0, 0.0 ], \ [ 2.0, 3.0, 0.0, 0.0 ], \ [ 4.0, 5.0, 6.0, 0.0 ], \ [ 7.0, 8.0, 9.0, 10.0 ] ] ) b = np.array ( [ 1.0, 8.0, 32.0, 90.0 ] ) print ( '' ) print ( 'r8mat_l_solve_test' ) print ( ' Python version: %s' % ( platform.python_version ( ) ) ) print ( ' r8mat_l_solve solves a lower triangular system.' ) r8mat_print ( n, n, a, ' Input matrix A:' ) r8vec_print ( n, b, ' Right hand side b:' ) x = r8mat_l_solve ( n, a, b ) r8vec_print ( n, x, ' Computed solution x:' ) r = np.dot ( a, x ) - b rnorm = np.linalg.norm ( r ) print ( '' ) print ( ' Norm of A*x-b = %g' % ( rnorm ) ) # # Terminate. # print ( '' ) print ( 'r8mat_l_solve_test' ) print ( ' Normal end of execution.' ) return def r8mat_lt_solve ( n, a, b ): #*****************************************************************************80 # ## r8mat_lt_solve() solves a transposed lower triangular linear system. # # Discussion: # # An R8MAT is an MxN array of R8's, stored by (I,J) -> [I+J*M]. # # Given the lower triangular matrix A, the linear system to be solved is: # # A' * x = b # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 28 August 2015 # # Author: # # John Burkardt # # Input: # # integer N, the number of rows and columns # of the matrix. # # real A(N,N), the N by N lower triangular matrix. # # real B(N,1), the right hand side of the linear system. # # Output: # # real X(N,1), the solution of the linear system. # import numpy as np # # Solve U' * x = b. # x = np.zeros ( n ) for i in range ( n - 1, -1, -1 ): x[i] = b[i] for j in range ( i + 1, n ): x[i] = x[i] - a[j,i] * x[j] x[i] = x[i] / a[i,i] return x def r8mat_lt_solve_test ( ): #*****************************************************************************80 # ## r8mat_lt_solve_test() tests r8mat_lt_solve(). # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 28 August 2016 # # Author: # # John Burkardt # import numpy as np import platform n = 4 a = np.array ( [ \ [ 1.0, 0.0, 0.0, 0.0 ], \ [ 2.0, 3.0, 0.0, 0.0 ], \ [ 4.0, 5.0, 6.0, 0.0 ], \ [ 7.0, 8.0, 9.0, 10.0 ] ] ) b = np.array ( [ 45.0, 53.0, 54.0, 40.0 ] ) print ( '' ) print ( 'r8mat_lt_solve_test' ) print ( ' Python version: %s' % ( platform.python_version ( ) ) ) print ( ' r8mat_lt_solve solves a transposed lower triangular system.' ) r8mat_print ( n, n, a, ' Input matrix A:' ) r8vec_print ( n, b, ' Right hand side b:' ) x = r8mat_lt_solve ( n, a, b ) r8vec_print ( n, x, ' Computed solution x:' ) r = np.dot ( np.transpose ( a ), x ) - b rnorm = np.linalg.norm ( r ) print ( '' ) print ( ' Norm of A\'*x-b = %g' % ( rnorm ) ) # # Terminate. # print ( '' ) print ( 'r8mat_lt_solve_test' ) print ( ' Normal end of execution.' ) return def r8mat_print ( m, n, a, title ): #*****************************************************************************80 # ## r8mat_print() prints an R8MAT. # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 31 August 2014 # # Author: # # John Burkardt # # Input: # # integer M, the number of rows in A. # # integer N, the number of columns in A. # # real A(M,N), the matrix. # # string TITLE, a title. # r8mat_print_some ( m, n, a, 0, 0, m - 1, n - 1, title ) return def r8mat_print_test ( ): #*****************************************************************************80 # ## r8mat_print_test() tests r8mat_print(). # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 10 February 2015 # # Author: # # John Burkardt # import numpy as np import platform print ( '' ) print ( 'r8mat_print_test' ) print ( ' Python version: %s' % ( platform.python_version ( ) ) ) print ( ' r8mat_print prints an R8MAT.' ) m = 4 n = 6 v = np.array ( [ \ [ 11.0, 12.0, 13.0, 14.0, 15.0, 16.0 ], [ 21.0, 22.0, 23.0, 24.0, 25.0, 26.0 ], [ 31.0, 32.0, 33.0, 34.0, 35.0, 36.0 ], [ 41.0, 42.0, 43.0, 44.0, 45.0, 46.0 ] ], dtype = np.float64 ) r8mat_print ( m, n, v, ' Here is an R8MAT:' ) # # Terminate. # print ( '' ) print ( 'r8mat_print_test:' ) print ( ' Normal end of execution.' ) return def r8mat_print_some ( m, n, a, ilo, jlo, ihi, jhi, title ): #*****************************************************************************80 # ## r8mat_print_some() prints out a portion of an R8MAT. # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 10 February 2015 # # Author: # # John Burkardt # # Input: # # integer M, N, the number of rows and columns of the matrix. # # real A(M,N), an M by N matrix to be printed. # # integer ILO, JLO, the first row and column to print. # # integer IHI, JHI, the last row and column to print. # # string TITLE, a title. # incx = 5 print ( '' ) print ( title ) if ( m <= 0 or n <= 0 ): print ( '' ) print ( ' (None)' ) return for j2lo in range ( max ( jlo, 0 ), min ( jhi + 1, n ), incx ): j2hi = j2lo + incx - 1 j2hi = min ( j2hi, n ) j2hi = min ( j2hi, jhi ) print ( '' ) print ( ' Col: ', end = '' ) for j in range ( j2lo, j2hi + 1 ): print ( '%7d ' % ( j ), end = '' ) print ( '' ) print ( ' Row' ) i2lo = max ( ilo, 0 ) i2hi = min ( ihi, m ) for i in range ( i2lo, i2hi + 1 ): print ( '%7d :' % ( i ), end = '' ) for j in range ( j2lo, j2hi + 1 ): print ( '%12g ' % ( a[i,j] ), end = '' ) print ( '' ) return def r8mat_print_some_test ( ): #*****************************************************************************80 # ## r8mat_print_some_test() tests r8mat_print_some(). # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 31 October 2014 # # Author: # # John Burkardt # import numpy as np import platform print ( '' ) print ( 'r8mat_print_some_test' ) print ( ' Python version: %s' % ( platform.python_version ( ) ) ) print ( ' r8mat_print_some prints some of an R8MAT.' ) m = 4 n = 6 v = np.array ( [ \ [ 11.0, 12.0, 13.0, 14.0, 15.0, 16.0 ], [ 21.0, 22.0, 23.0, 24.0, 25.0, 26.0 ], [ 31.0, 32.0, 33.0, 34.0, 35.0, 36.0 ], [ 41.0, 42.0, 43.0, 44.0, 45.0, 46.0 ] ], dtype = np.float64 ) r8mat_print_some ( m, n, v, 0, 3, 2, 5, ' Here is an R8MAT:' ) # # Terminate. # print ( '' ) print ( 'r8mat_print_some_test:' ) print ( ' Normal end of execution.' ) return def r8_mop ( i ): #*****************************************************************************80 # ## r8_mop() returns the I-th power of -1 as an R8 value. # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 01 June 2013 # # Author: # # John Burkardt # # Input: # # integer I, the power of -1. # # Output: # # real VALUE, the I-th power of -1. # if ( ( i % 2 ) == 0 ): value = + 1.0 else: value = - 1.0 return value def r8_mop_test ( ): #*****************************************************************************80 # ## r8_mop_test() tests r8_mop(). # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 06 December 2014 # # Author: # # John Burkardt # import numpy as np print ( '' ) print ( 'r8_mop_test' ) print ( ' r8_mop evaluates (-1.0)^I4 as an R8.' ) print ( '' ) print ( ' I4 r8_mop(I4)' ) print ( '' ) i4_min = -100; i4_max = +100; for test in range ( 0, 10 ): i4 = np.random.random_integers ( i4_min, i4_max ) r8 = r8_mop ( i4 ) print ( ' %4d %4.1f' % ( i4, r8 ) ) # # Terminate. # print ( '' ) print ( 'r8_mop_test' ) print ( ' Normal end of execution.' ) return def r8_sign ( x ): #*****************************************************************************80 # ## r8_sign() returns the sign of an R8. # # Discussion: # # The value is +1 if the number is positive or zero, and it is -1 otherwise. # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 03 June 2013 # # Author: # # John Burkardt # # Input: # # real X, the number whose sign is desired. # # Output: # # real VALUE, the sign of X. # if ( x < 0.0 ): value = -1.0 else: value = +1.0 return value def r8_sign_test ( ): #*****************************************************************************80 # ## r8_sign_test() tests r8_sign(). # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 28 September 2014 # # Author: # # John Burkardt # import numpy as np import platform test_num = 5 r8_test = np.array ( [ -1.25, -0.25, 0.0, +0.5, +9.0 ] ) print ( '' ) print ( 'r8_sign_test' ) print ( ' Python version: %s' % ( platform.python_version ( ) ) ) print ( ' r8_sign returns the sign of an R8.' ) print ( '' ) print ( ' R8 r8_sign(R8)' ) print ( '' ) for test in range ( 0, test_num ): r8 = r8_test[test] s = r8_sign ( r8 ) print ( ' %8.4f %8.0f' % ( r8, s ) ) # # Terminate. # print ( '' ) print ( 'r8_sign_test' ) print ( ' Normal end of execution.' ) return def r8vec_print ( n, a, title ): #*****************************************************************************80 # ## r8vec_print() prints an R8VEC. # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 31 August 2014 # # Author: # # John Burkardt # # Input: # # integer N, the dimension of the vector. # # real A(N), the vector to be printed. # # string TITLE, a title. # print ( '' ) print ( title ) print ( '' ) for i in range ( 0, n ): print ( '%6d: %12g' % ( i, a[i] ) ) def r8vec_print_test ( ): #*****************************************************************************80 # ## r8vec_print_test() tests r8vec_print(). # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 29 October 2014 # # Author: # # John Burkardt # import numpy as np import platform print ( '' ) print ( 'r8vec_print_test' ) print ( ' Python version: %s' % ( platform.python_version ( ) ) ) print ( ' r8vec_print prints an R8VEC.' ) n = 4 v = np.array ( [ 123.456, 0.000005, -1.0E+06, 3.14159265 ], dtype = np.float64 ) r8vec_print ( n, v, ' Here is an R8VEC:' ) # # Terminate. # print ( '' ) print ( 'r8vec_print_test:' ) print ( ' Normal end of execution.' ) return def svd_solve ( m, n, a, b ): #*****************************************************************************80 # ## svd_solve() solves a linear system in the least squares sense. # # Discussion: # # The vector X returned by this routine should always minimize the # Euclidean norm of the residual ||A*x-b||. # # If the matrix A does not have full column rank, then there are multiple # vectors that attain the minimum residual. In that case, the vector # X returned by this routine is the unique such minimizer that has the # the minimum possible Euclidean norm, that is, ||A*x-b|| and ||x|| # are both minimized. # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 31 August 2016 # # Author: # # John Burkardt # # Reference: # # David Kahaner, Cleve Moler, Steven Nash, # Numerical Methods and Software, # Prentice Hall, 1989, # ISBN: 0-13-627258-4, # LC: TA345.K34. # # Input: # # integer M, the number of rows of A. # # integer N, the number of columns of A. # # real A(M,N), the matrix. # # real B(M), the right hand side. # # Output: # # real X(N), the least squares solution. # import numpy as np # # Get the SVD. # a_copy = a lda = m ldu = m ldv = n job = 11 a_copy, sdiag, e, u, v, info = dsvdc ( a_copy, lda, m, n, ldu, ldv, job ) if ( info != 0 ): print ( '' ) print ( 'svd_solve - Fatal error!' ) print ( ' The SVD could not be calculated.' ) print ( ' LINPACK routine dsvdc returned a nonzero' ) print ( ' value of the error flag, INFO = %d' % ( info ) ) raise Exception ( 'svd_solve - Fatal error!' ) ub = np.dot ( u.transpose ( ), b ) sub = np.zeros ( n ) # # For singular problems, there may be tiny but nonzero singular values # that should be ignored. This is a reasonable attempt to avoid such # problems, although in general, the user might wish to control the tolerance. # smax = max ( sdiag ) eps = 2.220446049250313E-016 if ( smax <= eps ): smax = 1.0 stol = eps * smax for i in range ( 0, n ): if ( i <= m ): if ( stol <= sdiag[i] ): sub[i] = ub[i] / sdiag[i] x = np.dot ( v, sub ) return x def svd_solve_test ( ): #*****************************************************************************80 # ## svd_solve_test() tests svd_solve(). # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 31 August 2016 # # Author: # # John Burkardt # import numpy as np import platform print ( '' ) print ( 'svd_solve_test' ) print ( ' Python version: %s' % ( platform.python_version ( ) ) ) print ( ' svd_solve is a function with a simple interface which' ) print ( ' solves a linear system A*x = b in the least squares sense' ) print ( ' using the singular value decomposition (SVD).' ) print ( ' Compare a tabulated solution X1 to the qr_solve result X2.' ) prob_num = p00_prob_num ( ) print ( '' ) print ( ' Number of problems = %d' % ( prob_num ) ) print ( '' ) print ( ' Index M N ||B|| ||X1 - X2||' ), print ( ' ||X1|| ||X2|| ||R1|| ||R2||' ) print ( '' ) for prob in range ( 1, prob_num + 1 ): # # Get problem size. # m = p00_m ( prob ) n = p00_n ( prob ) # # Retrieve problem data. # a = p00_a ( prob, m, n ) b = p00_b ( prob, m ) x1 = p00_x ( prob, n ) b_norm = np.linalg.norm ( b ) x1_norm = np.linalg.norm ( x1 ) r1 = np.dot ( a, x1 ) - b r1_norm = np.linalg.norm ( r1 ) # # Use svd_solve on the problem. # Since we want to compute residuals ourselves, we need # to keep the originals of A and B, so we make copies # to send to svd_solve. # a_qr = a.copy ( ) b_qr = b.copy ( ) x2 = svd_solve ( m, n, a_qr, b_qr ) x2_norm = np.linalg.norm ( x2 ) r2 = np.dot ( a, x2 ) - b r2_norm = np.linalg.norm ( r2 ) # # Compare tabulated and computed solutions. # x_diff_norm = np.linalg.norm ( x1 - x2 ) # # Report results for this problem. # print ( ' %5d %4d %4d' % ( prob, m, n ) ), print ( ' %12.4g %12.4g' % ( b_norm, x_diff_norm ) ), print ( ' %12.4g %12.4g' % ( x1_norm, x2_norm ) ), print ( ' %12.4g %12.4g' % ( r1_norm, r2_norm ) ) # # Terminate. # print ( '' ) print ( 'svd_solve_test' ) print ( ' Normal end of execution.' ) return def p00_a ( prob, m, n ): #*****************************************************************************80 # ## p00_a() returns the matrix A for any least squares problem. # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 25 August 2016 # # Author: # # John Burkardt # # Input: # # integer PROB, the problem index. # # integer M, the number of equations. # # integer N, the number of variables. # # Output: # # real A(M,N), the matrix. # if ( prob == 1 ): a = p01_a ( m, n ) elif ( prob == 2 ): a = p02_a ( m, n ) elif ( prob == 3 ): a = p03_a ( m, n ) elif ( prob == 4 ): a = p04_a ( m, n ) elif ( prob == 5 ): a = p05_a ( m, n ) elif ( prob == 6 ): a = p06_a ( m, n ) else: print ( '' ) print ( 'p00_a - Fatal error!' ) print ( ' Illegal value of PROB = %d' % ( prob ) ) raise Exception ( 'p00_a - Fatal error!' ) return a def p00_b ( prob, m ): #*****************************************************************************80 # ## p00_b() returns the right hand side B for any least squares problem. # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 25 August 2016 # # Author: # # John Burkardt # # Input: # # integer PROB, the problem index. # # integer M, the number of equations. # # Output: # # real B(M), the right hand side. # if ( prob == 1 ): b = p01_b ( m ) elif ( prob == 2 ): b = p02_b ( m ) elif ( prob == 3 ): b = p03_b ( m ) elif ( prob == 4 ): b = p04_b ( m ) elif ( prob == 5 ): b = p05_b ( m ) elif ( prob == 6 ): b = p06_b ( m ) else: print ( '' ) print ( 'p00_b - Fatal error!' ) print ( ' Illegal value of PROB = %d' % ( prob ) ) raise Exception ( 'p00_b - Fatal error!' ) return b def p00_m ( prob ): #*****************************************************************************80 # ## p00_m() returns the number of equations M for any least squares problem. # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 25 August 2016 # # Author: # # John Burkardt # # Input: # # integer PROB, the problem index. # # Output: # # integer M, the number of equations. # if ( prob == 1 ): m = p01_m ( ) elif ( prob == 2 ): m = p02_m ( ) elif ( prob == 3 ): m = p03_m ( ) elif ( prob == 4 ): m = p04_m ( ) elif ( prob == 5 ): m = p05_m ( ) elif ( prob == 6 ): m = p06_m ( ) else: print ( '' ) print ( 'p00_m - Fatal error!' ) print ( ' Illegal value of PROB = %d' % ( prob ) ) raise Exception ( 'p00_m - Fatal error!' ) return m def p00_n ( prob ): #*****************************************************************************80 # ## p00_n() returns the number of variables N for any least squares problem. # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 25 August 2016 # # Author: # # John Burkardt # # Input: # # integer PROB, the problem index. # # Output: # # integer N, the number of variables. # if ( prob == 1 ): n = p01_n ( ) elif ( prob == 2 ): n = p02_n ( ) elif ( prob == 3 ): n = p03_n ( ) elif ( prob == 4 ): n = p04_n ( ) elif ( prob == 5 ): n = p05_n ( ) elif ( prob == 6 ): n = p06_n ( ) else: print ( '' ) print ( 'p00_n - Fatal error!' ) print ( ' Illegal value of PROB = %d' % ( prob ) ) raise Exception ( 'p00_n - Fatal error!' ) return n def p00_prob_num ( ): #*****************************************************************************80 # ## p00_prob_num() returns the number of least squares problems. # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 25 August 2016 # # Author: # # John Burkardt # # Output: # # real prob_num, the number of problems. # prob_num = 6 return prob_num def p00_x ( prob, n ): #*****************************************************************************80 # ## p00_x() returns the least squares solution X for any least squares problem. # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 25 August 2016 # # Author: # # John Burkardt # # Input: # # integer PROB, the problem index. # # integer N, the number of variables. # # Output: # # real X(N), the least squares solution. # if ( prob == 1 ): x = p01_x ( n ) elif ( prob == 2 ): x = p02_x ( n ) elif ( prob == 3 ): x = p03_x ( n ) elif ( prob == 4 ): x = p04_x ( n ) elif ( prob == 5 ): x = p05_x ( n ) elif ( prob == 6 ): x = p06_x ( n ) else: print ( '' ) print ( 'p00_x - Fatal error!' ) print ( ' Illegal value of PROB = %d' % ( prob ) ) raise Exception ( 'p00_x - Fatal error!' ) return x def p01_a ( m, n ): #*****************************************************************************80 # ## p01_a() returns the matrix A for problem 1. # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 25 August 2016 # # Author: # # John Burkardt # # Input: # # integer M, the number of equations. # # integer N, the number of variables. # # Output: # # real A(M,N), the matrix. # import numpy as np a = np.zeros ( [ m, n ] ) for i in range ( 0, m ): a[i,0] = 1.0 for j in range ( 1, n ): a[i,j] = a[i,j-1] * float ( i + 1 ) return a def p01_b ( m ): #*****************************************************************************80 # ## p01_b() returns the right hand side B for problem 1. # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 25 August 2016 # # Author: # # John Burkardt # # Input: # # integer M, the number of equations. # # Output: # # real B(M), the right hand side. # import numpy as np b = np.array ( [ 1.0, 2.3, 4.6, 3.1, 1.2 ] ) return b def p01_m ( ): #*****************************************************************************80 # ## p01_m() returns the number of equations M for problem 1. # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 25 August 2016 # # Author: # # John Burkardt # # Output: # # integer M, the number of equations. # m = 5 return m def p01_n ( ): #*****************************************************************************80 # ## p01_n() returns the number of variables N for problem 1. # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 25 August 2016 # # Author: # # John Burkardt # # Output: # # integer N, the number of variables. # n = 3 return n def p01_x ( n ): #*****************************************************************************80 # ## p01_x() returns the least squares solution X for problem 1. # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 25 August 2016 # # Author: # # John Burkardt # # Input: # # integer N, the number of variables. # # Output: # # real X(N), the least squares solution. # import numpy as np x = np.array ( [ -3.0200000, 4.4914286, -0.72857143 ] ) return x def p02_a ( m, n ): #*****************************************************************************80 # ## p02_a() returns the matrix A for problem 2. # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 25 August 2016 # # Author: # # John Burkardt # # Reference: # # Cleve Moler, # Numerical Computing with MATLAB, # SIAM, 2004, # ISBN13: 978-0-898716-60-3, # LC: QA297.M625, # ebook: https://www.mathworks.com/moler/chapters.html # # Input: # # integer M, the number of equations. # # Output: # # integer N, the number of variables. # # real A(M,N), the matrix. # import numpy as np a = np.zeros ( [ m, n ] ) for i in range ( 0, m ): a[i,n-1] = 1.0 for j in range ( n - 2, -1, -1 ): a[i,j] = a[i,j+1] * float ( i ) / 5.0 return a def p02_b ( m ): #*****************************************************************************80 # ## p02_b() returns the right hand side B for problem 2. # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 25 August 2016 # # Author: # # John Burkardt # # Input: # # integer M, the number of equations. # # Output: # # real B(M), the right hand side. # import numpy as np b = np.array ( [ 150.697, 179.323, 203.212, 226.505, 249.633, 281.422 ] ) return b def p02_m ( ): #*****************************************************************************80 # ## p02_m() returns the number of equations M for problem 2. # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 25 August 2016 # # Author: # # John Burkardt # # Output: # # integer M, the number of equations. # m = 6 return m def p02_n ( ): #*****************************************************************************80 # ## p02_n() returns the number of variables N for problem 2. # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 25 August 2016 # # Author: # # John Burkardt # # Output: # # integer N, the number of variables. # n = 3 return n def p02_x ( n ): #*****************************************************************************80 # ## p02_x() returns the least squares solution X for problem 2. # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 25 August 2016 # # Author: # # John Burkardt # # Input: # # integer N, the number of variables. # # Output: # # real X(N), the least squares solution. # import numpy as np x = np.array ( [ 5.7013, 121.1341, 152.4745 ] ) return x def p03_a ( m, n ): #*****************************************************************************80 # ## p03_a() returns the matrix A for problem 3. # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 25 August 2016 # # Author: # # John Burkardt # # Reference: # # Cleve Moler, # Numerical Computing with MATLAB, # SIAM, 2004, # ISBN13: 978-0-898716-60-3, # LC: QA297.M625, # ebook: https://www.mathworks.com/moler/chapters.html # # Input: # # integer M, the number of equations. # # integer N, the number of variables. # # Output: # # real A(M,N), the matrix. # import numpy as np a = np.zeros ( [ m, n ] ) a = np.array ( [ \ [ 1.0, 2.0, 3.0 ], \ [ 4.0, 5.0, 6.0 ], \ [ 7.0, 8.0, 9.0 ], \ [ 10.0, 11.0, 12.0 ], \ [ 13.0, 14.0, 15.0 ] ] ) return a def p03_b ( m ): #*****************************************************************************80 # ## p03_b() returns the right hand side B for problem 3. # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 25 August 2016 # # Author: # # John Burkardt # # Input: # # integer M, the number of equations. # # Output: # # real B(M), the right hand side. # import numpy as np b = np.array ( [ 16.0, 17.0, 18.0, 19.0, 20.0 ] ) return b def p03_m ( ): #*****************************************************************************80 # ## p03_m() returns the number of equations M for problem 3. # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 25 August 2016 # # Author: # # John Burkardt # # Output: # # integer M, the number of equations. # m = 5 return m def p03_n ( ): #*****************************************************************************80 # ## p03_n() returns the number of variables N for problem 3. # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 25 August 2016 # # Author: # # John Burkardt # # Output: # # integer N, the number of variables. # n = 3 return n def p03_x ( n ): #*****************************************************************************80 # ## p03_x() returns the least squares solution X for problem 3. # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 25 August 2016 # # Author: # # John Burkardt # # Input: # # integer N, the number of variables. # # Output: # # real X(N), the least squares solution. # import numpy as np x = np.array ( [ -7.5555556, 0.1111111, 7.7777778 ] ) return x def p04_a ( m, n ): #*****************************************************************************80 # ## p04_a() returns the matrix A for problem 4. # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 25 August 2016 # # Author: # # John Burkardt # # Input: # # integer M, the number of equations. # # integer N, the number of variables. # # Output: # # real A(M,N), the matrix. # import numpy as np a = np.zeros ( [ m, n ] ) for j in range ( 0, n ): for i in range ( 0, m ): a[i,j] = float ( j + 1 ) ** i return a def p04_b ( m ): #*****************************************************************************80 # ## p04_b() returns the right hand side B for problem 4. # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 25 August 2016 # # Author: # # John Burkardt # # Input: # # integer M, the number of equations. # # Output: # # real B(M), the right hand side. # import numpy as np b = np.array ( [ 15.0, 55.0, 225.0 ] ) return b def p04_m ( ): #*****************************************************************************80 # ## p04_m() returns the number of equations M for problem 4. # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 25 August 2016 # # Author: # # John Burkardt # # Output: # # integer M, the number of equations. # m = 3 return m def p04_n ( ): #*****************************************************************************80 # ## p04_n() returns the number of variables N for problem 4. # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 25 August 2016 # # Author: # # John Burkardt # # Output: # # integer N, the number of variables. # n = 5 return n def p04_x ( n ): #*****************************************************************************80 # ## p04_x() returns the least squares solution X for problem 4. # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 25 August 2016 # # Author: # # John Burkardt # # Input: # # integer N, the number of variables. # # Output: # # real X(N), the least squares solution. # import numpy as np x = np.array ( [ 1.0, 2.0, 3.0, 4.0, 5.0 ] ) return x def p05_a ( m, n ): #*****************************************************************************80 # ## p05_a() returns the matrix A for problem 5. # # Discussion: # # A is the Hilbert matrix. # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 25 August 2016 # # Author: # # John Burkardt # # Input: # # integer M, the number of equations. # # integer N, the number of variables. # # Output: # # real A(M,N), the matrix. # import numpy as np a = np.zeros ( [ m, n ] ) for j in range ( 0, n ): for i in range ( 0, m ): a[i,j] = 1.0 / float ( i + j + 1 ) return a def p05_b ( m ): #*****************************************************************************80 # ## p05_b() returns the right hand side B for problem 5. # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 25 August 2016 # # Author: # # John Burkardt # # Input: # # integer M, the number of equations. # # Output: # # real B(M), the right hand side. # import numpy as np b = np.zeros ( m ) b[0] = 1.0 return b def p05_m ( ): #*****************************************************************************80 # ## p05_m() returns the number of equations M for problem 5. # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 25 August 2016 # # Author: # # John Burkardt # # Output: # # integer M, the number of equations. # m = 10 return m def p05_n ( ): #*****************************************************************************80 # ## p05_n() returns the number of variables N for problem 5. # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 25 August 2016 # # Author: # # John Burkardt # # Output: # # integer N, the number of variables. # n = 10 return n def p05_x ( n ): #*****************************************************************************80 # ## p05_x() returns the least squares solution X for problem 5. # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 25 August 2016 # # Author: # # John Burkardt # # Input: # # integer N, the number of variables. # # Output: # # real X(N), the least squares solution. # from scipy.special import comb import numpy as np x = np.zeros ( n ) for i in range ( 0, n ): x[i] = r8_mop ( i + 2 ) * float ( i + 1 ) \ * comb ( n + i, n - 1 ) * comb ( n, n - i - 1 ) return x def p06_a ( m, n ): #*****************************************************************************80 # ## p06_a() returns the matrix A for problem 6. # # Discussion: # # A is a symmetric, orthogonal matrix. # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 25 August 2016 # # Author: # # John Burkardt # # Input: # # integer M, the number of equations. # # integer N, the number of variables. # # Output: # # real A(M,N), the matrix. # import numpy as np a = np.zeros ( [ m, n ] ) for i in range ( 0, m ): for j in range ( 0, n ): angle = float ( i + 1 ) * float ( j + 1 ) * np.pi / float ( n + 1 ) a[i,j] = np.sin ( angle ) a[0:m,0:n] = a[0:m,0:n] * np.sqrt ( 2.0 / ( n + 1 ) ) return a def p06_b ( m ): #*****************************************************************************80 # ## p06_b() returns the right hand side B for problem 6. # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 25 August 2016 # # Author: # # John Burkardt # # Input: # # integer M, the number of equations. # # Output: # # real B(M), the right hand side. # import numpy as np b = np.zeros ( m ) b[0] = 1.0 return b def p06_m ( ): #*****************************************************************************80 # ## p06_m() returns the number of equations M for problem 6. # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 25 August 2016 # # Author: # # John Burkardt # # Output: # # integer M, the number of equations. # m = 10 return m def p06_n ( ): #*****************************************************************************80 # ## p06_n() returns the number of variables N for problem 6. # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 25 August 2016 # # Author: # # John Burkardt # # Output: # # integer N, the number of variables. # n = 10 return n def p06_x ( n ): #*****************************************************************************80 # ## p06_x() returns the least squares solution X for problem 6. # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 25 August 2016 # # Author: # # John Burkardt # # Input: # # integer N, the number of variables. # # Output: # # real X(N), the least squares solution. # import numpy as np x = np.zeros ( n ) for i in range ( 0, n ): angle = float ( i + 1 ) * np.pi / float ( n + 1 ) x[i] = np.sin ( angle ) x[0:n] = x[0:n] * np.sqrt ( 2.0 / float ( n + 1) ) return x def r8_mop ( i ): #*****************************************************************************80 # ## r8_mop() returns the I-th power of -1 as an R8 value. # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 01 June 2013 # # Author: # # John Burkardt # # Input: # # integer I, the power of -1. # # Output: # # real VALUE, the I-th power of -1. # if ( ( i % 2 ) == 0 ): value = + 1.0 else: value = - 1.0 return value def r8_mop_test ( ): #*****************************************************************************80 # ## r8_mop_test() tests r8_mop(). # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 06 December 2014 # # Author: # # John Burkardt # import numpy as np print ( '' ) print ( 'r8_mop_test' ) print ( ' Python version: %s' % ( platform.python_version ( ) ) ) print ( ' r8_mop evaluates (-1.0)^I4 as an R8.' ) print ( '' ) print ( ' I4 r8_mop(I4)' ) print ( '' ) i4_min = -100 i4_max = +100 for test in range ( 0, 10 ): i4 = np.random.random_integers ( i4_min, i4_max ) r8 = r8_mop ( i4 ) print ( ' %4d %4.1f' % ( i4, r8 ) ) # # Terminate. # print ( '' ) print ( 'r8_mop_test' ) print ( ' Normal end of execution.' ) return def test_ls_test ( ): #*****************************************************************************80 # ## test_ls_test() tests test_ls(). # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 25 August 2016 # # Author: # # John Burkardt # import platform print ( '' ) print ( 'test_ls_test' ) print ( ' Python version: %s' % ( platform.python_version ( ) ) ) print ( ' Test test_ls.' ) ls_data_test ( ) lstsq_test ( ) # # Terminate. # print ( '' ) print ( 'test_ls_test' ) print ( ' Normal end of execution.' ) return def lstsq_test ( ): #*****************************************************************************80 # ## lstsq_test() tests the NUMPY LINALG least squares solver on the data. # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 25 August 2016 # # Author: # # John Burkardt # import numpy as np import platform print ( '' ) print ( 'lstsq_test' ) print ( ' Python version: %s' % ( platform.python_version ( ) ) ) print ( ' LSTSQ is the NUMPY LINALG least squares solver.' ) prob_num = p00_prob_num ( ) print ( '' ) print ( ' Number of problems = %d' % ( prob_num ) ) print ( '' ) print ( ' Index M N ||B|| ||X1-X2|| ||X1|| ' ), print ( ' ||X2|| ||R1|| ||R2||' ) print ( '' ) for prob in range ( 1, prob_num ): m = p00_m ( prob ) n = p00_n ( prob ) a = p00_a ( prob, m, n ) b = p00_b ( prob, m ) x1 = p00_x ( prob, n ) r1 = np.dot ( a, x1 ) - b b_norm = np.linalg.norm ( b ) x1_norm = np.linalg.norm ( x1 ) r1_norm = np.linalg.norm ( r1 ) [ x2, resids, rank, s ] = np.linalg.lstsq ( a, b ) r2 = np.dot ( a, x2 ) - b x2_norm = np.linalg.norm ( x2 ) r2_norm = np.linalg.norm ( r2 ) x_diff_norm = np.linalg.norm ( x1 - x2 ) print ( ' %5d %4d %4d %12.4g %12.4g %12.4g %12.4g %12.4g %12.4g' \ % ( prob, m, n, b_norm, x_diff_norm, x1_norm, x2_norm, r1_norm, r2_norm ) ) return def ls_data_test ( ): #*****************************************************************************80 # ## ls_data_test() summarizes the test data. # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 25 August 2016 # # Author: # # John Burkardt # import numpy as np import platform print ( '' ) print ( 'ls_data_test' ) print ( ' Python version: %s' % ( platform.python_version ( ) ) ) print ( ' Get each least squares test and compute the maximum residual.' ) print ( ' The L2 norm of the residual MUST be no greater than' ) print ( ' the L2 norm of the right hand side, else 0 is a better solution.' ) prob_num = p00_prob_num ( ) print ( '' ) print ( ' Number of problems = %d' % ( prob_num ) ) print ( '' ) print ( ' Index M N ||B|| ||X|| ||R||' ) print ( '' ) for prob in range ( 1, prob_num + 1 ): m = p00_m ( prob ) n = p00_n ( prob ) a = p00_a ( prob, m, n ) b = p00_b ( prob, m ) x = p00_x ( prob, n ) r = np.dot ( a, x ) - b b_norm = np.linalg.norm ( b ) x_norm = np.linalg.norm ( x ) r_norm = np.linalg.norm ( r ) print ( ' %5d %4d %4d %12.4g %12.4g %12.4g' \ % ( prob, m, n, b_norm, x_norm, r_norm ) ) return def timestamp ( ): #*****************************************************************************80 # ## timestamp() prints the date as a timestamp. # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 06 April 2013 # # Author: # # John Burkardt # import time t = time.time ( ) print ( time.ctime ( t ) ) return None if ( __name__ == '__main__' ): timestamp ( ) qr_solve_tests ( ) timestamp ( )