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