# include # include # include # include int main ( ); void daxpy ( int n, double da, double dx[], int incx, double dy[], int incy ); double ddot ( int n, double dx[], int incx, double dy[], int incy ); double dnrm2 ( int n, double x[], int incx ); void drot ( int n, double x[], int incx, double y[], int incy, double c, double s ); void drotg ( double *sa, double *sb, double *c, double *s ); void dscal ( int n, double sa, double x[], int incx ); int dsvdc ( double a[], int lda, int m, int n, double s[], double e[], double u[], int ldu, double v[], int ldv, double work[], int job ); void dswap ( int n, double x[], int incx, double y[], int incy ); int i4_max ( int i1, int i2 ); int i4_min ( int i1, int i2 ); double r8_sign ( double x ); void r8mat_print ( int m, int n, double a[], char *title ); void r8mat_print_some ( int m, int n, double a[], int ilo, int jlo, int ihi, int jhi, char *title ); double *r8mat_transpose_new ( int m, int n, double a[] ); double *r8mat_uniform_01_new ( int m, int n, int *seed ); void svd_truncated_u ( int m, int n, double a[], double un[], double sn[], double v[] ); void svd_truncated_u_test ( int m, int n ); void svd_truncated_v ( int m, int n, double a[], double u[], double sm[], double vm[] ); void svd_truncated_v_test ( int m, int n ); void timestamp ( ); /******************************************************************************/ int main ( ) /******************************************************************************/ /* Purpose: svd_truncated() analyzes the truncated SVD algorithm. Licensing: This code is distributed under the MIT license. Modified: 19 March 2012 Author: John Burkardt */ { int m; int n; timestamp ( ); printf ( "\n" ); printf ( "SVD_TRUNCATED\n" ); printf ( " C version\n" ); printf ( " Demonstrate the use of the truncated or economy-size\n" ); printf ( " Singular Value Decomposition (SVD) for cases where\n" ); printf ( " the sizes of M and N are very different.\n" ); m = 4; n = 3; svd_truncated_u_test ( m, n ); m = 3; n = 4; svd_truncated_v_test ( m, n ); printf ( "\n" ); printf ( "SVD_TRUNCATED\n" ); printf ( " Normal end of execution.\n" ); printf ( "\n" ); timestamp ( ); return 0; } /******************************************************************************/ void daxpy ( int n, double da, double dx[], int incx, double dy[], int incy ) /******************************************************************************/ /* Purpose: DAXPY computes constant times a vector plus a vector. Discussion: This routine uses unrolled loops for increments equal to one. Licensing: This code is distributed under the MIT license. Modified: 30 March 2007 Author: Original FORTRAN77 version by Charles Lawson, Richard Hanson, David Kincaid, Fred Krogh. C version by John Burkardt. Reference: Jack Dongarra, Cleve Moler, Jim Bunch, 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. Parameters: Input, int N, the number of elements in DX and DY. Input, double DA, the multiplier of DX. Input, double DX[*], the first vector. Input, int INCX, the increment between successive entries of DX. Input/output, double DY[*], the second vector. On output, DY[*] has been replaced by DY[*] + DA * DX[*]. Input, int INCY, the increment between successive entries of DY. */ { int i; int ix; int iy; int m; if ( n <= 0 ) { return; } if ( da == 0.0 ) { return; } /* Code for unequal increments or equal increments not equal to 1. */ if ( incx != 1 || incy != 1 ) { if ( 0 <= incx ) { ix = 0; } else { ix = ( - n + 1 ) * incx; } if ( 0 <= incy ) { iy = 0; } else { iy = ( - n + 1 ) * incy; } for ( i = 0; i < n; i++ ) { dy[iy] = dy[iy] + da * dx[ix]; ix = ix + incx; iy = iy + incy; } } /* Code for both increments equal to 1. */ else { m = n % 4; for ( i = 0; i < m; i++ ) { dy[i] = dy[i] + da * dx[i]; } for ( i = m; i < n; i = i + 4 ) { dy[i ] = dy[i ] + da * dx[i ]; dy[i+1] = dy[i+1] + da * dx[i+1]; dy[i+2] = dy[i+2] + da * dx[i+2]; dy[i+3] = dy[i+3] + da * dx[i+3]; } } return; } /******************************************************************************/ double ddot ( int n, double dx[], int incx, double dy[], int incy ) /******************************************************************************/ /* Purpose: DDOT forms the dot product of two vectors. Discussion: This routine uses unrolled loops for increments equal to one. Licensing: This code is distributed under the MIT license. Modified: 30 March 2007 Author: Original FORTRAN77 version by Charles Lawson, Richard Hanson, David Kincaid, Fred Krogh. C version by John Burkardt. Reference: Jack Dongarra, Cleve Moler, Jim Bunch, 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. Parameters: Input, int N, the number of entries in the vectors. Input, double DX[*], the first vector. Input, int INCX, the increment between successive entries in DX. Input, double DY[*], the second vector. Input, int INCY, the increment between successive entries in DY. Output, double DDOT, the sum of the product of the corresponding entries of DX and DY. */ { double dtemp; int i; int ix; int iy; int m; dtemp = 0.0; if ( n <= 0 ) { return dtemp; } /* Code for unequal increments or equal increments not equal to 1. */ if ( incx != 1 || incy != 1 ) { if ( 0 <= incx ) { ix = 0; } else { ix = ( - n + 1 ) * incx; } if ( 0 <= incy ) { iy = 0; } else { iy = ( - n + 1 ) * incy; } for ( i = 0; i < n; i++ ) { dtemp = dtemp + dx[ix] * dy[iy]; ix = ix + incx; iy = iy + incy; } } /* Code for both increments equal to 1. */ else { m = n % 5; for ( i = 0; i < m; i++ ) { dtemp = dtemp + dx[i] * dy[i]; } for ( i = m; i < n; i = i + 5 ) { dtemp = dtemp + dx[i ] * dy[i ] + dx[i+1] * dy[i+1] + dx[i+2] * dy[i+2] + dx[i+3] * dy[i+3] + dx[i+4] * dy[i+4]; } } return dtemp; } /******************************************************************************/ double dnrm2 ( int n, double x[], int incx ) /******************************************************************************/ /* Purpose: DNRM2 returns the euclidean norm of a vector. Discussion: DNRM2 ( X ) = sqrt ( X' * X ) Licensing: This code is distributed under the MIT license. Modified: 30 March 2007 Author: Original FORTRAN77 version by Charles Lawson, Richard Hanson, David Kincaid, Fred Krogh. C version by John Burkardt. Reference: Jack Dongarra, Cleve Moler, Jim Bunch, 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. Parameters: Input, int N, the number of entries in the vector. Input, double X[*], the vector whose norm is to be computed. Input, int INCX, the increment between successive entries of X. Output, double DNRM2, the Euclidean norm of X. */ { double absxi; int i; int ix; double norm; double scale; double ssq; if ( n < 1 || incx < 1 ) { norm = 0.0; } else if ( n == 1 ) { norm = fabs ( x[0] ); } else { scale = 0.0; ssq = 1.0; ix = 0; for ( i = 0; i < n; i++ ) { if ( x[ix] != 0.0 ) { absxi = fabs ( x[ix] ); if ( scale < absxi ) { ssq = 1.0 + ssq * ( scale / absxi ) * ( scale / absxi ); scale = absxi; } else { ssq = ssq + ( absxi / scale ) * ( absxi / scale ); } } ix = ix + incx; } norm = scale * sqrt ( ssq ); } return norm; } /******************************************************************************/ void drot ( int n, double x[], int incx, double y[], int incy, double c, double s ) /******************************************************************************/ /* Purpose: DROT applies a plane rotation. Licensing: This code is distributed under the MIT license. Modified: 30 March 2007 Author: Original FORTRAN77 version by Charles Lawson, Richard Hanson, David Kincaid, Fred Krogh. C version by John Burkardt. Reference: Jack Dongarra, Cleve Moler, Jim Bunch, 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. Parameters: Input, int N, the number of entries in the vectors. Input/output, double X[*], one of the vectors to be rotated. Input, int INCX, the increment between successive entries of X. Input/output, double Y[*], one of the vectors to be rotated. Input, int INCY, the increment between successive elements of Y. Input, double C, S, parameters (presumably the cosine and sine of some angle) that define a plane rotation. */ { int i; int ix; int iy; double stemp; if ( n <= 0 ) { } else if ( incx == 1 && incy == 1 ) { for ( i = 0; i < n; i++ ) { stemp = c * x[i] + s * y[i]; y[i] = c * y[i] - s * x[i]; x[i] = stemp; } } else { if ( 0 <= incx ) { ix = 0; } else { ix = ( - n + 1 ) * incx; } if ( 0 <= incy ) { iy = 0; } else { iy = ( - n + 1 ) * incy; } for ( i = 0; i < n; i++ ) { stemp = c * x[ix] + s * y[iy]; y[iy] = c * y[iy] - s * x[ix]; x[ix] = stemp; ix = ix + incx; iy = iy + incy; } } return; } /******************************************************************************/ void drotg ( double *sa, double *sb, double *c, double *s ) /******************************************************************************/ /* Purpose: 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: 30 March 2007 Author: Original FORTRAN77 version by Charles Lawson, Richard Hanson, David Kincaid, Fred Krogh. C version by John Burkardt. Reference: Jack Dongarra, Cleve Moler, Jim Bunch, 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. Parameters: Input/output, double *SA, *SB, On input, SA and SB are the values A and B. On output, SA is overwritten with R, and SB is overwritten with Z. Output, double *C, *S, the cosine and sine of the Givens rotation. */ { double r; double roe; double scale; double z; if ( fabs ( *sb ) < fabs ( *sa ) ) { roe = *sa; } else { roe = *sb; } scale = fabs ( *sa ) + fabs ( *sb ); if ( scale == 0.0 ) { *c = 1.0; *s = 0.0; r = 0.0; } else { r = scale * sqrt ( ( *sa / scale ) * ( *sa / scale ) + ( *sb / scale ) * ( *sb / scale ) ); r = r8_sign ( roe ) * r; *c = *sa / r; *s = *sb / r; } if ( 0.0 < fabs ( *c ) && fabs ( *c ) <= *s ) { z = 1.0 / *c; } else { z = *s; } *sa = r; *sb = z; return; } /******************************************************************************/ void dscal ( int n, double sa, double x[], int incx ) /******************************************************************************/ /* Purpose: DSCAL scales a vector by a constant. Licensing: This code is distributed under the MIT license. Modified: 30 March 2007 Author: Original FORTRAN77 version by Charles Lawson, Richard Hanson, David Kincaid, Fred Krogh. C version by John Burkardt. Reference: Jack Dongarra, Cleve Moler, Jim Bunch, 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. Parameters: Input, int N, the number of entries in the vector. Input, double SA, the multiplier. Input/output, double X[*], the vector to be scaled. Input, int INCX, the increment between successive entries of X. */ { int i; int ix; int m; if ( n <= 0 ) { } else if ( incx == 1 ) { m = n % 5; for ( i = 0; i < m; i++ ) { x[i] = sa * x[i]; } for ( i = m; i < n; i = i + 5 ) { x[i] = sa * x[i]; x[i+1] = sa * x[i+1]; x[i+2] = sa * x[i+2]; x[i+3] = sa * x[i+3]; x[i+4] = sa * x[i+4]; } } else { if ( 0 <= incx ) { ix = 0; } else { ix = ( - n + 1 ) * incx; } for ( i = 0; i < n; i++ ) { x[ix] = sa * x[ix]; ix = ix + incx; } } return; } /******************************************************************************/ int dsvdc ( double a[], int lda, int m, int n, double s[], double e[], double u[], int ldu, double v[], int ldv, double work[], int job ) /******************************************************************************/ /* Purpose: 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: 03 May 2007 Author: Original FORTRAN77 version by Jack Dongarra, Cleve Moler, Jim Bunch, Pete Stewart, C version by John Burkardt. Reference: Jack Dongarra, Cleve Moler, Jim Bunch and Pete 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/output, double A[LDA*N]. On input, the M by N matrix whose singular value decomposition is to be computed. On output, the matrix has been destroyed. Depending on the user's requests, the matrix may contain other useful information. Input, int LDA, the leading dimension of the array A. LDA must be at least M. Input, int M, the number of rows of the matrix. Input, int N, the number of columns of the matrix A. Output, double S[MM], where MM = min(M+1,N). The first min(M,N) entries of S contain the singular values of A arranged in descending order of magnitude. Output, double E[MM], where MM = min(M+1,N), ordinarily contains zeros. However see the discussion of INFO for exceptions. Output, double 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. Input, int LDU, the leading dimension of the array U. LDU must be at least M. Output, double 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. Input, int LDV, the leading dimension of the array V. LDV must be at least N. Workspace, double WORK[M]. Input, int JOB, controls the computation of the singular vectors. It has the decimal expansion AB with the following meaning: A = 0, do 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, do not compute the right singular vectors. B = 1, return the right singular vectors in V. Output, int *DSVDC, status indicator INFO. 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. */ { double b; double c; double cs; double el; double emm1; double f; double g; int i; int info; int iter; int j; int jobu; int k; int kase; int kk; int l; int ll; int lls; int ls; int lu; int maxit = 30; int mm; int mm1; int mn; int nct; int nctp1; int ncu; int nrt; int nrtp1; double scale; double shift; double sl; double sm; double smm1; double sn; double t; double t1; double test; int wantu; int wantv; double ztest; /* Determine what is to be computed. */ info = 0; wantu = 0; wantv = 0; jobu = ( job % 100 ) / 10; if ( 1 < jobu ) { ncu = i4_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. */ nct = i4_min ( m-1, n ); nrt = i4_max ( 0, i4_min ( m, n-2 ) ); lu = i4_max ( nct, nrt ); for ( l = 1; l <= lu; l++ ) { /* Compute the transformation for the L-th column and place the L-th diagonal in S(L). */ if ( l <= nct ) { s[l-1] = dnrm2 ( m-l+1, a+l-1+(l-1)*lda, 1 ); if ( s[l-1] != 0.0 ) { if ( a[l-1+(l-1)*lda] != 0.0 ) { s[l-1] = r8_sign ( a[l-1+(l-1)*lda] ) * fabs ( s[l-1] ); } dscal ( m-l+1, 1.0 / s[l-1], a+l-1+(l-1)*lda, 1 ); a[l-1+(l-1)*lda] = 1.0 + a[l-1+(l-1)*lda]; } s[l-1] = -s[l-1]; } for ( j = l+1; j <= n; j++ ) { /* Apply the transformation. */ if ( l <= nct && s[l-1] != 0.0 ) { t = - ddot ( m-l+1, a+l-1+(l-1)*lda, 1, a+l-1+(j-1)*lda, 1 ) / a[l-1+(l-1)*lda]; daxpy ( m-l+1, t, a+l-1+(l-1)*lda, 1, a+l-1+(j-1)*lda, 1 ); } /* Place the L-th row of A into E for the subsequent calculation of the row transformation. */ e[j-1] = a[l-1+(j-1)*lda]; } /* Place the transformation in U for subsequent back multiplication. */ if ( wantu && l <= nct ) { for ( i = l; i <= m; i++ ) { u[i-1+(l-1)*ldu] = a[i-1+(l-1)*lda]; } } if ( l <= nrt ) { /* Compute the L-th row transformation and place the L-th superdiagonal in E(L). */ e[l-1] = dnrm2 ( n-l, e+l, 1 ); if ( e[l-1] != 0.0 ) { if ( e[l] != 0.0 ) { e[l-1] = r8_sign ( e[l] ) * fabs ( e[l-1] ); } dscal ( n-l, 1.0 / e[l-1], e+l, 1 ); e[l] = 1.0 + e[l]; } e[l-1] = -e[l-1]; /* Apply the transformation. */ if ( l+1 <= m && e[l-1] != 0.0 ) { for ( j = l+1; j <= m; j++ ) { work[j-1] = 0.0; } for ( j = l+1; j <= n; j++ ) { daxpy ( m-l, e[j-1], a+l+(j-1)*lda, 1, work+l, 1 ); } for ( j = l+1; j <= n; j++ ) { daxpy ( m-l, -e[j-1]/e[l], work+l, 1, a+l+(j-1)*lda, 1 ); } } /* Place the transformation in V for subsequent back multiplication. */ if ( wantv ) { for ( j = l+1; j <= n; j++ ) { v[j-1+(l-1)*ldv] = e[j-1]; } } } } /* Set up the final bidiagonal matrix of order MN. */ mn = i4_min ( m + 1, n ); nctp1 = nct + 1; nrtp1 = nrt + 1; if ( nct < n ) { s[nctp1-1] = a[nctp1-1+(nctp1-1)*lda]; } if ( m < mn ) { s[mn-1] = 0.0; } if ( nrtp1 < mn ) { e[nrtp1-1] = a[nrtp1-1+(mn-1)*lda]; } e[mn-1] = 0.0; /* If required, generate U. */ if ( wantu ) { for ( i = 1; i <= m; i++ ) { for ( j = nctp1; j <= ncu; j++ ) { u[(i-1)+(j-1)*ldu] = 0.0; } } for ( j = nctp1; j <= ncu; j++ ) { u[j-1+(j-1)*ldu] = 1.0; } for ( ll = 1; ll <= nct; ll++ ) { l = nct - ll + 1; if ( s[l-1] != 0.0 ) { for ( j = l+1; j <= ncu; j++ ) { t = - ddot ( m-l+1, u+(l-1)+(l-1)*ldu, 1, u+(l-1)+(j-1)*ldu, 1 ) / u[l-1+(l-1)*ldu]; daxpy ( m-l+1, t, u+(l-1)+(l-1)*ldu, 1, u+(l-1)+(j-1)*ldu, 1 ); } dscal ( m-l+1, -1.0, u+(l-1)+(l-1)*ldu, 1 ); u[l-1+(l-1)*ldu] = 1.0 + u[l-1+(l-1)*ldu]; for ( i = 1; i <= l-1; i++ ) { u[i-1+(l-1)*ldu] = 0.0; } } else { for ( i = 1; i <= m; i++ ) { u[i-1+(l-1)*ldu] = 0.0; } u[l-1+(l-1)*ldu] = 1.0; } } } /* If it is required, generate V. */ if ( wantv ) { for ( ll = 1; ll <= n; ll++ ) { l = n - ll + 1; if ( l <= nrt && e[l-1] != 0.0 ) { for ( j = l+1; j <= n; j++ ) { t = - ddot ( n-l, v+l+(l-1)*ldv, 1, v+l+(j-1)*ldv, 1 ) / v[l+(l-1)*ldv]; daxpy ( n-l, t, v+l+(l-1)*ldv, 1, v+l+(j-1)*ldv, 1 ); } } for ( i = 1; i <= n; i++ ) { v[i-1+(l-1)*ldv] = 0.0; } v[l-1+(l-1)*ldv] = 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 ) { info = mn; return 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 ( ll = 1; ll <= mn; ll++ ) { l = mn - ll; if ( l == 0 ) { break; } test = fabs ( s[l-1] ) + fabs ( s[l] ); ztest = test + fabs ( e[l-1] ); if ( ztest == test ) { e[l-1] = 0.0; break; } } if ( l == mn - 1 ) { kase = 4; } else { for ( lls = l + 1; lls <= mn + 1; lls++ ) { ls = mn - lls + l + 1; if ( ls == l ) { break; } test = 0.0; if ( ls != mn ) { test = test + fabs ( e[ls-1] ); } if ( ls != l + 1 ) { test = test + fabs ( e[ls-2] ); } ztest = test + fabs ( s[ls-1] ); if ( ztest == test ) { s[ls-1] = 0.0; break; } } if ( ls == l ) { kase = 3; } else if ( ls == mn ) { 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 ( kk = 1; kk <= mm1; kk++ ) { k = mm1 - kk + l; t1 = s[k-1]; drotg ( &t1, &f, &cs, &sn ); s[k-1] = t1; if ( k != l ) { f = -sn * e[k-2]; e[k-2] = cs * e[k-2]; } if ( wantv ) { drot ( n, v+0+(k-1)*ldv, 1, v+0+(mn-1)*ldv, 1, cs, sn ); } } } /* Split at negligible S(L). */ else if ( kase == 2 ) { f = e[l-2]; e[l-2] = 0.0; for ( k = l; k <= mn; k++ ) { t1 = s[k-1]; drotg ( &t1, &f, &cs, &sn ); s[k-1] = t1; f = - sn * e[k-1]; e[k-1] = cs * e[k-1]; if ( wantu ) { drot ( m, u+0+(k-1)*ldu, 1, u+0+(l-2)*ldu, 1, cs, sn ); } } } /* Perform one QR step. */ else if ( kase == 3 ) { /* Calculate the shift. */ scale = fmax ( fabs ( s[mn-1] ), fmax ( fabs ( s[mn-2] ), fmax ( fabs ( e[mn-2] ), fmax ( fabs ( s[l-1] ), fabs ( e[l-1] ) ) ) ) ); sm = s[mn-1] / scale; smm1 = s[mn-2] / scale; emm1 = e[mn-2] / scale; sl = s[l-1] / scale; el = e[l-1] / scale; b = ( ( smm1 + sm ) * ( smm1 - sm ) + emm1 * emm1 ) / 2.0; c = ( sm * emm1 ) * ( sm * emm1 ); shift = 0.0; if ( b != 0.0 || c != 0.0 ) { shift = 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 = l; k <= mm1; k++ ) { drotg ( &f, &g, &cs, &sn ); if ( k != l ) { e[k-2] = f; } f = cs * s[k-1] + sn * e[k-1]; e[k-1] = cs * e[k-1] - sn * s[k-1]; g = sn * s[k]; s[k] = cs * s[k]; if ( wantv ) { drot ( n, v+0+(k-1)*ldv, 1, v+0+k*ldv, 1, cs, sn ); } drotg ( &f, &g, &cs, &sn ); s[k-1] = f; f = cs * e[k-1] + sn * s[k]; s[k] = -sn * e[k-1] + cs * s[k]; g = sn * e[k]; e[k] = cs * e[k]; if ( wantu && k < m ) { drot ( m, u+0+(k-1)*ldu, 1, u+0+k*ldu, 1, cs, sn ); } } e[mn-2] = f; iter = iter + 1; } /* Convergence. */ else if ( kase == 4 ) { /* Make the singular value nonnegative. */ if ( s[l-1] < 0.0 ) { s[l-1] = -s[l-1]; if ( wantv ) { dscal ( n, -1.0, v+0+(l-1)*ldv, 1 ); } } /* Order the singular value. */ for ( ; ; ) { if ( l == mm ) { break; } if ( s[l] <= s[l-1] ) { break; } t = s[l-1]; s[l-1] = s[l]; s[l] = t; if ( wantv && l < n ) { dswap ( n, v+0+(l-1)*ldv, 1, v+0+l*ldv, 1 ); } if ( wantu && l < m ) { dswap ( m, u+0+(l-1)*ldu, 1, u+0+l*ldu, 1 ); } l = l + 1; } iter = 0; mn = mn - 1; } } return info; } /******************************************************************************/ void dswap ( int n, double x[], int incx, double y[], int incy ) /******************************************************************************/ /* Purpose: DSWAP interchanges two vectors. Licensing: This code is distributed under the MIT license. Modified: 30 March 2007 Author: Original FORTRAN77 version by Charles Lawson, Richard Hanson, David Kincaid, Fred Krogh. C version by John Burkardt. Reference: Jack Dongarra, Cleve Moler, Jim Bunch, 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. Parameters: Input, int N, the number of entries in the vectors. Input/output, double X[*], one of the vectors to swap. Input, int INCX, the increment between successive entries of X. Input/output, double Y[*], one of the vectors to swap. Input, int INCY, the increment between successive elements of Y. */ { int i; int ix; int iy; int m; double temp; if ( n <= 0 ) { } else if ( incx == 1 && incy == 1 ) { m = n % 3; for ( i = 0; i < m; i++ ) { temp = x[i]; x[i] = y[i]; y[i] = temp; } for ( i = m; i < n; i = i + 3 ) { temp = x[i]; x[i] = y[i]; y[i] = temp; temp = x[i+1]; x[i+1] = y[i+1]; y[i+1] = temp; temp = x[i+2]; x[i+2] = y[i+2]; y[i+2] = temp; } } else { if ( 0 <= incx ) { ix = 0; } else { ix = ( - n + 1 ) * incx; } if ( 0 <= incy ) { iy = 0; } else { iy = ( - n + 1 ) * incy; } for ( i = 0; i < n; i++ ) { temp = x[ix]; x[ix] = y[iy]; y[iy] = temp; ix = ix + incx; iy = iy + incy; } } return; } /******************************************************************************/ int i4_max ( int i1, int i2 ) /******************************************************************************/ /* Purpose: I4_MAX returns the maximum of two I4's. Licensing: This code is distributed under the MIT license. Modified: 29 August 2006 Author: John Burkardt Parameters: Input, int I1, I2, are two integers to be compared. Output, int I4_MAX, the larger of I1 and I2. */ { int value; if ( i2 < i1 ) { value = i1; } else { value = i2; } return value; } /******************************************************************************/ int i4_min ( int i1, int i2 ) /******************************************************************************/ /* Purpose: I4_MIN returns the smaller of two I4's. Licensing: This code is distributed under the MIT license. Modified: 29 August 2006 Author: John Burkardt Parameters: Input, int I1, I2, two integers to be compared. Output, int I4_MIN, the smaller of I1 and I2. */ { int value; if ( i1 < i2 ) { value = i1; } else { value = i2; } return value; } /******************************************************************************/ double r8_sign ( double x ) /******************************************************************************/ /* Purpose: R8_SIGN returns the sign of a R8. Licensing: This code is distributed under the MIT license. Modified: 18 October 2004 Author: John Burkardt Parameters: Input, double X, the number whose sign is desired. Output, double R8_SIGN, the sign of X. */ { double value; if ( x < 0.0 ) { value = -1.0; } else { value = 1.0; } return value; } /******************************************************************************/ void r8mat_print ( int m, int n, double a[], char *title ) /******************************************************************************/ /* Purpose: R8MAT_PRINT prints an R8MAT. Discussion: An R8MAT is a doubly dimensioned array of R8 values, stored as a vector in column-major order. Entry A(I,J) is stored as A[I+J*M] Licensing: This code is distributed under the MIT license. Modified: 28 May 2008 Author: John Burkardt Parameters: Input, int M, the number of rows in A. Input, int N, the number of columns in A. Input, double A[M*N], the M by N matrix. Input, char *TITLE, a title. */ { r8mat_print_some ( m, n, a, 1, 1, m, n, title ); return; } /******************************************************************************/ void r8mat_print_some ( int m, int n, double a[], int ilo, int jlo, int ihi, int jhi, char *title ) /******************************************************************************/ /* Purpose: R8MAT_PRINT_SOME prints some of an R8MAT. Discussion: An R8MAT is a doubly dimensioned array of R8 values, stored as a vector in column-major order. Licensing: This code is distributed under the MIT license. Modified: 26 June 2013 Author: John Burkardt Parameters: Input, int M, the number of rows of the matrix. M must be positive. Input, int N, the number of columns of the matrix. N must be positive. Input, double A[M*N], the matrix. Input, int ILO, JLO, IHI, JHI, designate the first row and column, and the last row and column to be printed. Input, char *TITLE, a title. */ { # define INCX 5 int i; int i2hi; int i2lo; int j; int j2hi; int j2lo; fprintf ( stdout, "\n" ); fprintf ( stdout, "%s\n", title ); if ( m <= 0 || n <= 0 ) { fprintf ( stdout, "\n" ); fprintf ( stdout, " (None)\n" ); return; } /* Print the columns of the matrix, in strips of 5. */ for ( j2lo = jlo; j2lo <= jhi; j2lo = j2lo + INCX ) { j2hi = j2lo + INCX - 1; if ( n < j2hi ) { j2hi = n; } if ( jhi < j2hi ) { j2hi = jhi; } fprintf ( stdout, "\n" ); /* For each column J in the current range... Write the header. */ fprintf ( stdout, " Col: "); for ( j = j2lo; j <= j2hi; j++ ) { fprintf ( stdout, " %7d ", j - 1 ); } fprintf ( stdout, "\n" ); fprintf ( stdout, " Row\n" ); fprintf ( stdout, "\n" ); /* Determine the range of the rows in this strip. */ if ( 1 < ilo ) { i2lo = ilo; } else { i2lo = 1; } if ( m < ihi ) { i2hi = m; } else { i2hi = ihi; } for ( i = i2lo; i <= i2hi; i++ ) { /* Print out (up to) 5 entries in row I, that lie in the current strip. */ fprintf ( stdout, "%5d:", i - 1 ); for ( j = j2lo; j <= j2hi; j++ ) { fprintf ( stdout, " %14f", a[i-1+(j-1)*m] ); } fprintf ( stdout, "\n" ); } } return; # undef INCX } /******************************************************************************/ double *r8mat_transpose_new ( int m, int n, double a[] ) /******************************************************************************/ /* Purpose: R8MAT_TRANSPOSE_NEW returns the transpose of a matrix. Discussion: An R8MAT is a doubly dimensioned array of R8 values, stored as a vector in column-major order. Licensing: This code is distributed under the MIT license. Modified: 27 June 2008 Author: John Burkardt Parameters: Input, int M, N, the number of rows and columns of the matrix A. Input, double A[M*N], the matrix whose transpose is desired. Output, double R8MAT_TRANSPOSE_NEW[N*M], the transposed matrix. */ { double *b; int i; int j; b = ( double * ) malloc ( n * m * sizeof ( double ) ); for ( j = 0; j < n; j++ ) { for ( i = 0; i < m; i++ ) { b[j+i*n] = a[i+j*m]; } } return b; } /******************************************************************************/ double *r8mat_uniform_01_new ( int m, int n, int *seed ) /******************************************************************************/ /* Purpose: R8MAT_UNIFORM_01_NEW fills an R8MAT with pseudorandom values scaled to [0,1]. Discussion: An R8MAT is a doubly dimensioned array of R8 values, stored as a vector in column-major order. This routine implements the recursion seed = 16807 * seed mod ( 2^31 - 1 ) unif = seed / ( 2^31 - 1 ) The integer arithmetic never requires more than 32 bits, including a sign bit. Licensing: This code is distributed under the MIT license. Modified: 30 June 2009 Author: John Burkardt Reference: Paul Bratley, Bennett Fox, Linus Schrage, A Guide to Simulation, Springer Verlag, pages 201-202, 1983. Bennett Fox, Algorithm 647: Implementation and Relative Efficiency of Quasirandom Sequence Generators, ACM Transactions on Mathematical Software, Volume 12, Number 4, pages 362-376, 1986. Philip Lewis, Allen Goodman, James Miller, A Pseudo-Random Number Generator for the System/360, IBM Systems Journal, Volume 8, pages 136-143, 1969. Parameters: Input, int M, N, the number of rows and columns. Input/output, int *SEED, the "seed" value. Normally, this value should not be 0, otherwise the output value of SEED will still be 0, and R8_UNIFORM will be 0. On output, SEED has been updated. Output, double R8MAT_UNIFORM_01_NEW[M*N], a matrix of pseudorandom values. */ { int i; int j; int k; double *r; r = ( double * ) malloc ( m * n * sizeof ( double ) ); for ( j = 0; j < n; j++ ) { for ( i = 0; i < m; i++ ) { k = *seed / 127773; *seed = 16807 * ( *seed - k * 127773 ) - k * 2836; if ( *seed < 0 ) { *seed = *seed + 2147483647; } r[i+j*m] = ( double ) ( *seed ) * 4.656612875E-10; } } return r; } /******************************************************************************/ void svd_truncated_u ( int m, int n, double a[], double un[], double sn[], double v[] ) /******************************************************************************/ /* Purpose: SVD_TRUNCATED_U gets the truncated SVD when N <= M Discussion: A(mxn) = U(mxm) * S(mxn) * V(nxn)' = Un(mxn) * Sn(nxn) * V(nxn)' Licensing: This code is distributed under the MIT license. Modified: 05 July 2013 Author: John Burkardt Parameters: Input, int M, N, the number of rows and columns in the matrix A. Input, double A[M*N], the matrix whose singular value decomposition we are investigating. Output, double UN[M*N], SN[N*N], V[N*N], the factors that form the singular value decomposition of A. */ { double *a_copy; double *e; int i; int info; int j; int lda; int ldu; int ldv; int job; double *sdiag; double *work; /* The correct size of E and SDIAG is min ( m+1, n). */ a_copy = ( double * ) malloc ( m * n * sizeof ( double ) ); e = ( double * ) malloc ( ( m + n ) * sizeof ( double ) ); sdiag = ( double * ) malloc ( ( m + n ) * sizeof ( double ) ); work = ( double * ) malloc ( m * sizeof ( double ) ); /* Compute the eigenvalues and eigenvectors. */ job = 21; lda = m; ldu = m; ldv = n; /* The input matrix is destroyed by the routine. Since we need to keep it around, we only pass a copy to the routine. */ for ( j = 0; j < n; j++ ) { for ( i = 0; i < m; i++ ) { a_copy[i+j*m] = a[i+j*m]; } } info = dsvdc ( a_copy, lda, m, n, sdiag, e, un, ldu, v, ldv, work, job ); if ( info != 0 ) { printf ( "\n" ); printf ( "SVD_TRUNCATED_U - Failure!\n" ); printf ( " The SVD could not be calculated.\n" ); printf ( " LINPACK routine DSVDC returned a nonzero\n" ); printf ( " value of the error flag, INFO = %d\n", info ); return; } /* Make the NxN matrix S from the diagonal values in SDIAG. */ for ( j = 0; j < n; j++ ) { for ( i = 0; i < n; i++ ) { if ( i == j ) { sn[i+j*n] = sdiag[i]; } else { sn[i+j*n] = 0.0; } } } free ( a_copy ); free ( e ); free ( sdiag ); free ( work ); return; } /******************************************************************************/ void svd_truncated_u_test ( int m, int n ) /******************************************************************************/ /* Purpose: SVD_TRUNCATED_U_TEST tests SVD_TRUNCATED_U. Licensing: This code is distributed under the MIT license. Modified: 05 July 2013 Author: John Burkardt */ { double *a; double *a_save; double err; int i; int j; int k; int seed; double *sn; double *un; double *v; printf ( "\n" ); printf ( "SVD_TRUNCATED_U_TEST\n" ); printf ( " M = %d\n", m ); printf ( " N = %d\n", n ); a = ( double * ) malloc ( m * n * sizeof ( double ) ); un = ( double * ) malloc ( m * n * sizeof ( double ) ); sn = ( double * ) malloc ( n * n * sizeof ( double ) ); v = ( double * ) malloc ( n * n * sizeof ( double ) ); seed = 123456789; a_save = r8mat_uniform_01_new ( m, n, &seed ); r8mat_print ( m, n, a_save, " A:" ); for ( j = 0; j < n; j++ ) { for ( i = 0; i < m; i++ ) { a[i+j*m] = a_save[i+j*m]; } } svd_truncated_u ( m, n, a, un, sn, v ); r8mat_print ( m, n, un, " UN:" ); r8mat_print ( n, n, sn, " SN:" ); r8mat_print ( n, n, v, " V:" ); /* Check the factorization by computing A = U * S * V' */ for ( j = 0; j < n; j++) { for ( i = 0; i < m; i++ ) { a[i+j*m] = 0.0; for ( k = 0; k < n; k++ ) { a[i+j*m] = a[i+j*m] + un[i+k*m] * sn[k+k*n] * v[j+k*n]; } } } err = 0.0; for ( j = 0; j < n; j++ ) { for ( i = 0; i < m; i++ ) { err = fmax ( err, fabs ( a[i+j*m] - a_save[i+j*m] ) ); } } printf ( "\n" ); printf ( " Maximum error |A - U*S*V'| = %g\n", err ); r8mat_print ( m, n, a, " Recomputed A = U * S * V':" ); free ( a ); free ( a_save ); free ( sn ); free ( un ); free ( v ); return; } /******************************************************************************/ void svd_truncated_v ( int m, int n, double a[], double u[], double sm[], double vm[] ) /******************************************************************************/ /* Purpose: SVD_TRUNCATED_V gets the truncated SVD when M <= N. Discussion: A(mxn) = U(mxm) * S(mxn) * V(nxn)' = U(mxm) * Sm(mxm) * Vm(mxn)' Licensing: This code is distributed under the MIT license. Modified: 05 July 2013 Author: John Burkardt Parameters: Input, int M, N, the number of rows and columns in the matrix A. Input, double A[M*N], the matrix whose singular value decomposition we are investigating. Output, double U[M*M], SM[M*M], VM[M*N], the factors that form the singular value decomposition of A. */ { double *a2; int m2; int n2; /* Transpose the matrix! */ a2 = r8mat_transpose_new ( m, n, a ); m2 = n; n2 = m; svd_truncated_u ( m2, n2, a2, vm, sm, u ); free ( a2 ); return; } /******************************************************************************/ void svd_truncated_v_test ( int m, int n ) /******************************************************************************/ /* Purpose: SVD_TRUNCATED_V_TEST tests SVD_TRUNCATED_V. Licensing: This code is distributed under the MIT license. Modified: 05 July 2013 Author: John Burkardt */ { double *a; double *a_save; double err; int i; int j; int k; int seed; double *sm; double *u; double *vm; printf ( "\n" ); printf ( "SVD_TRUNCATED_V_TEST\n" ); printf ( " M = %d\n", m ); printf ( " N = %d\n", n ); a = ( double * ) malloc ( m * n * sizeof ( double ) ); u = ( double * ) malloc ( m * m * sizeof ( double ) ); sm = ( double * ) malloc ( m * m * sizeof ( double ) ); vm = ( double * ) malloc ( n * m * sizeof ( double ) ); seed = 123456789; a_save = r8mat_uniform_01_new ( m, n, &seed ); r8mat_print ( m, n, a_save, " A:" ); for ( j = 0; j < n; j++ ) { for ( i = 0; i < m; i++ ) { a[i+j*m] = a_save[i+j*m]; } } svd_truncated_v ( m, n, a, u, sm, vm ); r8mat_print ( m, m, u, " U:" ); r8mat_print ( m, m, sm, " SM:" ); r8mat_print ( n, m, vm, " VM:" ); /* Check the factorization by computing A = U * S * V' */ for ( j = 0; j < n; j++) { for ( i = 0; i < m; i++ ) { a[i+j*m] = 0.0; for ( k = 0; k < m; k++ ) { a[i+j*m] = a[i+j*m] + u[i+k*m] * sm[k+k*m] * vm[j+k*n]; } } } err = 0.0; for ( j = 0; j < n; j++ ) { for ( i = 0; i < m; i++ ) { err = fmax ( err, fabs ( a[i+j*m] - a_save[i+j*m] ) ); } } printf ( "\n" ); printf ( " Maximum error |A - U*S*V'| = %g\n", err ); r8mat_print ( m, n, a, " Recomputed A = U * S * V':" ); free ( a ); free ( a_save ); free ( sm ); free ( u ); free ( vm ); return; } /******************************************************************************/ void timestamp ( void ) /******************************************************************************/ /* Purpose: TIMESTAMP prints the current YMDHMS date as a time stamp. Example: 31 May 2001 09:45:54 AM Licensing: This code is distributed under the MIT license. Modified: 24 September 2003 Author: John Burkardt Parameters: None */ { # define TIME_SIZE 40 static char time_buffer[TIME_SIZE]; const struct tm *tm; time_t now; now = time ( NULL ); tm = localtime ( &now ); strftime ( time_buffer, TIME_SIZE, "%d %B %Y %I:%M:%S %p", tm ); fprintf ( stdout, "%s\n", time_buffer ); return; # undef TIME_SIZE }