#! /usr/bin/env python # 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 GNU LGPL 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 # # Parameters: # # Input, real A(LDA,N), the M by N matrix whose singular # value decomposition is to be computed. # # Input, integer LDA, the leading dimension of the array A. # LDA must be at least N. # # Input, integer M, the number of rows of the matrix. # # Input, integer N, the number of columns of the matrix A. # # Input, integer LDU, the leading dimension of the array U. # LDU must be at least M. # # Input, integer LDV, the leading dimension of the array V. # LDV must be at least N. # # Input, 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. # # Output, 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. # # Output, real E(MM), where MM = max(M+1,N), ordinarily contains zeros. # However see the discussion of INFO for exceptions. # # Output, 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. # # Output, 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. # # Output, 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 from drot import drot from drotg import drotg 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. # # Discussion: # # DSVDC computes the singular value decomposition: # # A = U * S * V' # # Licensing: # # This code is distributed under the GNU LGPL license. # # Modified: # # 07 September 2016 # # Author: # # John Burkardt # import numpy as np import platform from r8mat_print import r8mat_print from r8mat_uniform_01 import r8mat_uniform_01 from r8vec_print import r8vec_print 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. # seed = 123456789 a, seed = r8mat_uniform_01 ( m, n, seed ) 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 if ( __name__ == '__main__' ): from timestamp import timestamp timestamp ( ) dsvdc_test ( ) timestamp ( )