# include # include # include # include # include # include "svd_snowfall.h" /******************************************************************************/ char ch_cap ( char ch ) /******************************************************************************/ /* Purpose: CH_CAP capitalizes a single character. Discussion: This routine should be equivalent to the library "toupper" function. Licensing: This code is distributed under the MIT license. Modified: 19 July 1998 Author: John Burkardt Parameters: Input, char CH, the character to capitalize. Output, char CH_CAP, the capitalized character. */ { if ( 97 <= ch && ch <= 122 ) { ch = ch - 32; } return ch; } /******************************************************************************/ int ch_eqi ( char ch1, char ch2 ) /******************************************************************************/ /* Purpose: CH_EQI is TRUE (1) if two characters are equal, disregarding case. Licensing: This code is distributed under the MIT license. Modified: 13 June 2003 Author: John Burkardt Parameters: Input, char CH1, CH2, the characters to compare. Output, int CH_EQI, is TRUE (1) if the two characters are equal, disregarding case and FALSE (0) otherwise. */ { int value; if ( 97 <= ch1 && ch1 <= 122 ) { ch1 = ch1 - 32; } if ( 97 <= ch2 && ch2 <= 122 ) { ch2 = ch2 - 32; } if ( ch1 == ch2 ) { value = 1; } else { value = 0; } return value; } /******************************************************************************/ int ch_to_digit ( char ch ) /******************************************************************************/ /* Purpose: CH_TO_DIGIT returns the integer value of a base 10 digit. Example: CH DIGIT --- ----- '0' 0 '1' 1 ... ... '9' 9 ' ' 0 'X' -1 Licensing: This code is distributed under the MIT license. Modified: 13 June 2003 Author: John Burkardt Parameters: Input, char CH, the decimal digit, '0' through '9' or blank are legal. Output, int CH_TO_DIGIT, the corresponding integer value. If the character was 'illegal', then DIGIT is -1. */ { int digit; if ( '0' <= ch && ch <= '9' ) { digit = ch - '0'; } else if ( ch == ' ' ) { digit = 0; } else { digit = -1; } return digit; } /******************************************************************************/ 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 file_column_count ( char *input_filename ) /******************************************************************************/ /* Purpose: FILE_COLUMN_COUNT counts the number of columns in the first line of a file. Discussion: The file is assumed to be a simple text file. Most lines of the file is presumed to consist of COLUMN_NUM words, separated by spaces. There may also be some blank lines, and some comment lines, which have a "#" in column 1. The routine tries to find the first non-comment non-blank line and counts the number of words in that line. If all lines are blanks or comments, it goes back and tries to analyze a comment line. Licensing: This code is distributed under the MIT license. Modified: 13 June 2003 Author: John Burkardt Parameters: Input, char *INPUT_FILENAME, the name of the file. Output, int FILE_COLUMN_COUNT, the number of columns assumed to be in the file. */ { # define LINE_MAX 255 int column_num; char *error; FILE *input; int got_one; char line[LINE_MAX]; /* Open the file. */ input = fopen ( input_filename, "r" ); if ( !input ) { fprintf ( stderr, "\n" ); fprintf ( stderr, "FILE_COLUMN_COUNT - Fatal error!\n" ); fprintf ( stderr, " Could not open the input file: \"%s\"\n", input_filename ); exit ( 1 ); } /* Read one line, but skip blank lines and comment lines. */ got_one = 0; for ( ; ; ) { error = fgets ( line, LINE_MAX, input ); if ( !error ) { break; } if ( s_len_trim ( line ) == 0 ) { continue; } if ( line[0] == '#' ) { continue; } got_one = 1; break; } if ( got_one == 0 ) { fclose ( input ); input = fopen ( input_filename, "r" ); for ( ; ; ) { error = fgets ( line, LINE_MAX, input ); if ( !error ) { break; } if ( s_len_trim ( line ) == 0 ) { continue; } got_one = 1; break; } } fclose ( input ); if ( got_one == 0 ) { fprintf ( stderr, "\n" ); fprintf ( stderr, "FILE_COLUMN_COUNT - Warning!\n" ); fprintf ( stderr, " The file does not seem to contain any data.\n" ); exit ( 1 ); } column_num = s_word_count ( line ); return column_num; # undef LINE_MAX } /******************************************************************************/ int file_row_count ( char *input_filename ) /******************************************************************************/ /* Purpose: FILE_ROW_COUNT counts the number of row records in a file. Discussion: It does not count lines that are blank, or that begin with a comment symbol '#'. Licensing: This code is distributed under the MIT license. Modified: 13 June 2003 Author: John Burkardt Parameters: Input, char *INPUT_FILENAME, the name of the input file. Output, int FILE_ROW_COUNT, the number of rows found. */ { # define LINE_MAX 255 int comment_num; char *error; FILE *input; char line[LINE_MAX]; int record_num; int row_num; row_num = 0; comment_num = 0; record_num = 0; input = fopen ( input_filename, "r" ); if ( !input ) { fprintf ( stderr, "\n" ); fprintf ( stderr, "FILE_ROW_COUNT - Fatal error!\n" ); fprintf ( stderr, " Could not open the input file: \"%s\"\n", input_filename ); exit ( 1 ); } for ( ; ; ) { error = fgets ( line, LINE_MAX, input ); if ( !error ) { break; } record_num = record_num + 1; if ( line[0] == '#' ) { comment_num = comment_num + 1; continue; } if ( s_len_trim ( line ) == 0 ) { comment_num = comment_num + 1; continue; } row_num = row_num + 1; } fclose ( input ); return row_num; # undef LINE_MAX } /******************************************************************************/ 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 r8col_normalize_li ( int m, int n, double a[] ) /******************************************************************************/ /* Purpose: R8COL_NORMALIZE_LI normalizes an R8COL with the column infinity norm. Discussion: Each column is scaled so that the entry of maximum norm has the value 1. Licensing: This code is distributed under the MIT license. Modified: 08 February 2012 Author: John Burkardt Parameters: Input, int M, N, the number of rows and columns. Input/output, double A[M*N], the array to be normalized. */ { double c; int i; int j; for ( j = 0; j < n; j++ ) { c = a[0+j*m]; for ( i = 1; i < m; i++ ) { if ( fabs ( c ) < fabs ( a[i+j*m] ) ) { c = a[i+j*m]; } } if ( c != 0.0 ) { for ( i = 0; i < m; i++ ) { a[i+m*j] = a[i+m*j] / c; } } } return; } /******************************************************************************/ void r8col_reverse ( int m, int n, double a[] ) /******************************************************************************/ /* Purpose: R8COL_REVERSE reverses the order of columns in an R8COL. Discussion: To reverse the columns is to start with something like 11 12 13 14 15 21 22 23 24 25 31 32 33 34 35 41 42 43 44 45 51 52 53 54 55 and return 15 14 13 12 11 25 24 23 22 21 35 34 33 32 31 45 44 43 42 41 55 54 53 52 51 Licensing: This code is distributed under the MIT license. Modified: 06 May 2013 Author: John Burkardt Parameters: Input, int M, N, the number of rows and columns. Input/output, double A[M*N], the matrix whose columns are to be flipped. */ { int i; int j; double t; for ( i = 0; i < m; i++ ) { for ( j = 0; j < ( n / 2 ); j++ ) { t = a[i+ j *m]; a[i+ j *m] = a[i+(n-1-j)*m]; a[i+(n-1-j)*m] = t; } } return; } /******************************************************************************/ double *r8mat_data_read ( char *input_filename, int m, int n ) /******************************************************************************/ /* Purpose: R8MAT_DATA_READ reads the data from an R8MAT file. Discussion: An R8MAT is an array of R8's. The file is assumed to contain one record per line. Records beginning with the '#' character are comments, and are ignored. Blank lines are also ignored. Each line that is not ignored is assumed to contain exactly (or at least) M real numbers, representing the coordinates of a point. There are assumed to be exactly (or at least) N such records. Licensing: This code is distributed under the MIT license. Modified: 27 January 2005 Author: John Burkardt Parameters: Input, char *INPUT_FILENAME, the name of the input file. Input, int M, the number of spatial dimensions. Input, int N, the number of points. The program will stop reading data once N values have been read. Output, double R8MAT_DATA_READ[M*N], the data. */ { # define LINE_MAX 255 int error; char *got_string; FILE *input; int i; int j; char line[255]; double *table; double *x; input = fopen ( input_filename, "r" ); if ( !input ) { fprintf ( stderr, "\n" ); fprintf ( stderr, "R8MAT_DATA_READ - Fatal error!\n" ); fprintf ( stderr, " Could not open the input file: \"%s\"\n", input_filename ); exit ( 1 ); } table = ( double * ) malloc ( m * n * sizeof ( double ) ); x = ( double * ) malloc ( m * sizeof ( double ) ); j = 0; while ( j < n ) { got_string = fgets ( line, LINE_MAX, input ); if ( !got_string ) { break; } if ( line[0] == '#' || s_len_trim ( line ) == 0 ) { continue; } error = s_to_r8vec ( line, m, x ); if ( error == 1 ) { continue; } for ( i = 0; i < m; i++ ) { table[i+j*m] = x[i]; } j = j + 1; } fclose ( input ); free ( x ); return table; # undef LINE_MAX } /******************************************************************************/ void r8mat_header_read ( char *input_filename, int *m, int *n ) /******************************************************************************/ /* Purpose: R8MAT_HEADER_READ reads the header from an R8MAT file. Discussion: An R8MAT is an array of R8's. Licensing: This code is distributed under the MIT license. Modified: 04 June 2004 Author: John Burkardt Parameters: Input, char *INPUT_FILENAME, the name of the input file. Output, int *M, the number of spatial dimensions. Output, int *N, the number of points. */ { *m = file_column_count ( input_filename ); if ( *m <= 0 ) { fprintf ( stderr, "\n" ); fprintf ( stderr, "R8MAT_HEADER_READ - Fatal error!\n" ); fprintf ( stderr, " FILE_COLUMN_COUNT failed.\n" ); exit ( 1 ); } *n = file_row_count ( input_filename ); if ( *n <= 0 ) { fprintf ( stderr, "\n" ); fprintf ( stderr, "R8MAT_HEADER_READ - Fatal error!\n" ); fprintf ( stderr, " FILE_ROW_COUNT failed.\n" ); exit ( 1 ); } return; } /******************************************************************************/ double *r8mat_mm_new ( int n1, int n2, int n3, double a[], double b[] ) /******************************************************************************/ /* Purpose: R8MAT_MM_NEW multiplies two matrices. Discussion: An R8MAT is a doubly dimensioned array of R8 values, stored as a vector in column-major order. For this routine, the result is returned as the function value. Licensing: This code is distributed under the MIT license. Modified: 08 April 2009 Author: John Burkardt Parameters: Input, int N1, N2, N3, the order of the matrices. Input, double A[N1*N2], double B[N2*N3], the matrices to multiply. Output, double R8MAT_MM[N1*N3], the product matrix C = A * B. */ { double *c; int i; int j; int k; c = ( double * ) malloc ( n1 * n3 * sizeof ( double ) ); for ( i = 0; i < n1; i++ ) { for ( j = 0; j < n3; j++ ) { c[i+j*n1] = 0.0; for ( k = 0; k < n2; k++ ) { c[i+j*n1] = c[i+j*n1] + a[i+k*n1] * b[k+j*n2]; } } } return c; } /******************************************************************************/ double *r8mat_mmt_new ( int n1, int n2, int n3, double a[], double b[] ) /******************************************************************************/ /* Purpose: R8MAT_MMT_NEW computes C = A * B'. Discussion: An R8MAT is a doubly dimensioned array of R8 values, stored as a vector in column-major order. For this routine, the result is returned as the function value. Licensing: This code is distributed under the MIT license. Modified: 13 November 2012 Author: John Burkardt Parameters: Input, int N1, N2, N3, the order of the matrices. Input, double A[N1*N2], double B[N3*N2], the matrices to multiply. Output, double R8MAT_MMT[N1*N3], the product matrix C = A * B'. */ { double *c; int i; int j; int k; c = ( double * ) malloc ( n1 * n3 * sizeof ( double ) ); for ( i = 0; i < n1; i++ ) { for ( j = 0; j < n3; j++ ) { c[i+j*n1] = 0.0; for ( k = 0; k < n2; k++ ) { c[i+j*n1] = c[i+j*n1] + a[i+k*n1] * b[j+k*n3]; } } } return c; } /******************************************************************************/ 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: 20 August 2010 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; j2hi = i4_min ( j2hi, n ); j2hi = i4_min ( 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. */ i2lo = i4_max ( ilo, 1 ); i2hi = i4_min ( ihi, m ); 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 } /******************************************************************************/ void r8mat_svd_linpack ( int m, int n, double a[], double u[], double s[], double v[] ) /******************************************************************************/ /* Purpose: R8MAT_SVD_LINPACK gets the SVD of a matrix using a call to LINPACK. Discussion: The singular value decomposition of a real MxN matrix A has the form: A = U * S * V' where U is MxM orthogonal, S is MxN, and entirely zero except for the diagonal; V is NxN orthogonal. Moreover, the nonzero entries of S are positive, and appear in order, from largest magnitude to smallest. This routine calls the LINPACK routine DSVDC to compute the factorization. Licensing: This code is distributed under the MIT license. Modified: 19 June 2012 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], S[M*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 = 11; 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, u, ldu, v, ldv, work, job ); if ( info != 0 ) { printf ( "\n" ); printf ( "R8MAT_SVD_LINPACK - 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 MxN matrix S from the diagonal values in SDIAG. */ for ( j = 0; j < n; j++ ) { for ( i = 0; i < m; i++ ) { if ( i == j ) { s[i+j*m] = sdiag[i]; } else { s[i+j*m] = 0.0; } } } /* Note that we do NOT need to transpose the V that comes out of LINPACK! */ free ( a_copy ); free ( e ); free ( sdiag ); free ( work ); return; } /******************************************************************************/ double *r8mat_svd_low_rank ( int m, int n, int r, double u[], double s[], double v[] ) /******************************************************************************/ /* Purpose: R8MAT_SVD_LOW_RANK forms a rank R approximation using the SVD. Discussion: The singular value decomposition of a real MxN matrix A has the form: A = U * S * V' where U is MxM orthogonal, S is MxN, and entirely zero except for the diagonal; V is NxN orthogonal. For R <= min ( M, N ), we can form a rank R approximation to A: Ar = Ur * Sr * Vr' where Ur is MxR orthogonal, Sr is RxR, and entirely zero except for the diagonal; Vr is NxR orthogonal. Strictly speaking, we use RM = min ( R, M ), RN = min ( R, N ): Ur is M x RM, Sr is RM x RN, Vr is N x RN Licensing: This code is distributed under the MIT license. Modified: 09 May 2013 Author: John Burkardt Parameters: Input, int M, N, the number of rows and columns in the matrix A. Input, int R, the desired rank of the approximation. Input, double U[M*M], S[M*N], V[N*N], the factors that form the singular value decomposition of A. Output, double R8MAT_SVD_LOW_RANK[M*N], the rank R SVD approximation to A. */ { int i; int j; int rm; int rn; double *s1; double *sv1; double *u1; double *usv1; double *v1; rm = i4_min ( m, r ); rn = i4_min ( n, r ); u1 = ( double * ) malloc ( m * rm * sizeof ( double ) ); s1 = ( double * ) malloc ( rm * rn * sizeof ( double ) ); v1 = ( double * ) malloc ( n * rn * sizeof ( double ) ); for ( j = 0; j < rn; j++ ) { for ( i = 0; i < rm; i++ ) { s1[i+j*rm] = s[i+j*m]; } } for ( j = 0; j < rn; j++ ) { for ( i = 0; i < n; i++ ) { v1[i+j*n] = v[i+j*n]; } } for ( j = 0; j < rm; j++ ) { for ( i = 0; i < m; i++ ) { u1[i+j*m] = u[i+j*m]; } } sv1 = r8mat_mmt_new ( rm, rn, n, s1, v1 ); usv1 = r8mat_mm_new ( m, rm, n, u1, sv1 ); free ( s1 ); free ( sv1 ); free ( u1 ); free ( v1 ); return usv1; } /******************************************************************************/ void r8row_reverse ( int m, int n, double a[] ) /******************************************************************************/ /* Purpose: R8ROW_REVERSE reverses the order of the rows of an R8ROW. Discussion: To reverse the rows is to start with something like 11 12 13 14 15 21 22 23 24 25 31 32 33 34 35 41 42 43 44 45 51 52 53 54 55 and return 51 52 53 54 55 41 42 43 44 45 31 32 33 34 35 21 22 23 24 25 11 12 13 14 15 Licensing: This code is distributed under the MIT license. Modified: 06 May 2013 Author: John Burkardt Parameters: Input, int M, N, the number of rows and columns. Input/output, double A[M*N], the matrix whose rows are to be flipped. */ { int i; int j; double t; for ( j = 0; j < n; j++ ) { for ( i = 0; i < ( m / 2 ); i++ ) { t = a[ i+j*m]; a[ i+j*m] = a[m-1-i+j*m]; a[m-1-i+j*m] = t; } } return; } /******************************************************************************/ double *r8vec_cum0_new ( int n, double a[] ) /******************************************************************************/ /* Purpose: R8VEC_CUM0_NEW computes the cumulutive sums of an R8VEC. Discussion: An R8VEC is a vector of R8's. Input: A = { 1.0, 2.0, 3.0, 4.0 } Output: A_CUM = { 0.0, 1.0, 3.0, 6.0, 10.0 } Licensing: This code is distributed under the MIT license. Modified: 07 May 2012 Author: John Burkardt Parameters: Input, int N, the number of entries in the vector. Input, double A[N], the vector to be summed. Output, double R8VEC_CUM0_NEW[N+1], the cumulative sums. */ { double *a_cum; int i; a_cum = ( double * ) malloc ( ( n + 1 ) * sizeof ( double ) ); a_cum[0] = 0.0; for ( i = 1; i <= n; i++ ) { a_cum[i] = a_cum[i-1] + a[i-1]; } return a_cum; } /******************************************************************************/ void r8vec_print ( int n, double a[], char *title ) /******************************************************************************/ /* Purpose: R8VEC_PRINT prints an R8VEC. Discussion: An R8VEC is a vector of R8's. Licensing: This code is distributed under the MIT license. Modified: 08 April 2009 Author: John Burkardt Parameters: Input, int N, the number of components of the vector. Input, double A[N], the vector to be printed. Input, char *TITLE, a title. */ { int i; fprintf ( stdout, "\n" ); fprintf ( stdout, "%s\n", title ); fprintf ( stdout, "\n" ); for ( i = 0; i < n; i++ ) { fprintf ( stdout, " %8d: %14f\n", i, a[i] ); } return; } /******************************************************************************/ double r8vec_sum ( int n, double a[] ) /******************************************************************************/ /* Purpose: R8VEC_SUM returns the sum of an R8VEC. Discussion: An R8VEC is a vector of R8's. Licensing: This code is distributed under the MIT license. Modified: 26 August 2008 Author: John Burkardt Parameters: Input, int N, the number of entries in the vector. Input, double A[N], the vector. Output, double R8VEC_SUM, the sum of the vector. */ { int i; double value; value = 0.0; for ( i = 0; i < n; i++ ) { value = value + a[i]; } return value; } /******************************************************************************/ int s_len_trim ( char *s ) /******************************************************************************/ /* Purpose: S_LEN_TRIM returns the length of a string to the last nonblank. Licensing: This code is distributed under the MIT license. Modified: 26 April 2003 Author: John Burkardt Parameters: Input, char *S, a pointer to a string. Output, int S_LEN_TRIM, the length of the string to the last nonblank. If S_LEN_TRIM is 0, then the string is entirely blank. */ { int n; char *t; n = strlen ( s ); t = s + strlen ( s ) - 1; while ( 0 < n ) { if ( *t != ' ' ) { return n; } t--; n--; } return n; } /******************************************************************************/ double s_to_r8 ( char *s, int *lchar, int *error ) /******************************************************************************/ /* Purpose: S_TO_R8 reads an R8 value from a string. Discussion: We have had some trouble with input of the form 1.0E-312. For now, let's assume anything less than 1.0E-20 is zero. This routine will read as many characters as possible until it reaches the end of the string, or encounters a character which cannot be part of the real number. Legal input is: 1 blanks, 2 '+' or '-' sign, 2.5 spaces 3 integer part, 4 decimal point, 5 fraction part, 6 'E' or 'e' or 'D' or 'd', exponent marker, 7 exponent sign, 8 exponent integer part, 9 exponent decimal point, 10 exponent fraction part, 11 blanks, 12 final comma or semicolon. with most quantities optional. Example: S R '1' 1.0 ' 1 ' 1.0 '1A' 1.0 '12,34,56' 12.0 ' 34 7' 34.0 '-1E2ABCD' -100.0 '-1X2ABCD' -1.0 ' 2E-1' 0.2 '23.45' 23.45 '-4.2E+2' -420.0 '17d2' 1700.0 '-14e-2' -0.14 'e2' 100.0 '-12.73e-9.23' -12.73 * 10.0^(-9.23) Licensing: This code is distributed under the MIT license. Modified: 24 June 2005 Author: John Burkardt Parameters: Input, char *S, the string containing the data to be read. Reading will begin at position 1 and terminate at the end of the string, or when no more characters can be read to form a legal real. Blanks, commas, or other nonnumeric data will, in particular, cause the conversion to halt. Output, int *LCHAR, the number of characters read from the string to form the number, including any terminating characters such as a trailing comma or blanks. Output, int *ERROR, is TRUE (1) if an error occurred and FALSE (0) otherwise. Output, double S_TO_R8, the value that was read from the string. */ { char c; int ihave; int isgn; int iterm; int jbot; int jsgn; int jtop; int nchar; int ndig; double r; double rbot; double rexp; double rtop; char TAB = 9; nchar = s_len_trim ( s ); *error = 0; r = 0.0; *lchar = -1; isgn = 1; rtop = 0.0; rbot = 1.0; jsgn = 1; jtop = 0; jbot = 1; ihave = 1; iterm = 0; for ( ; ; ) { c = s[*lchar+1]; *lchar = *lchar + 1; /* Blank or TAB character. */ if ( c == ' ' || c == TAB ) { if ( ihave == 2 ) { } else if ( ihave == 6 || ihave == 7 ) { iterm = 1; } else if ( 1 < ihave ) { ihave = 11; } } /* Comma. */ else if ( c == ',' || c == ';' ) { if ( ihave != 1 ) { iterm = 1; ihave = 12; *lchar = *lchar + 1; } } /* Minus sign. */ else if ( c == '-' ) { if ( ihave == 1 ) { ihave = 2; isgn = -1; } else if ( ihave == 6 ) { ihave = 7; jsgn = -1; } else { iterm = 1; } } /* Plus sign. */ else if ( c == '+' ) { if ( ihave == 1 ) { ihave = 2; } else if ( ihave == 6 ) { ihave = 7; } else { iterm = 1; } } /* Decimal point. */ else if ( c == '.' ) { if ( ihave < 4 ) { ihave = 4; } else if ( 6 <= ihave && ihave <= 8 ) { ihave = 9; } else { iterm = 1; } } /* Exponent marker. */ else if ( ch_eqi ( c, 'E' ) || ch_eqi ( c, 'D' ) ) { if ( ihave < 6 ) { ihave = 6; } else { iterm = 1; } } /* Digit. */ else if ( ihave < 11 && '0' <= c && c <= '9' ) { if ( ihave <= 2 ) { ihave = 3; } else if ( ihave == 4 ) { ihave = 5; } else if ( ihave == 6 || ihave == 7 ) { ihave = 8; } else if ( ihave == 9 ) { ihave = 10; } ndig = ch_to_digit ( c ); if ( ihave == 3 ) { rtop = 10.0 * rtop + ( double ) ndig; } else if ( ihave == 5 ) { rtop = 10.0 * rtop + ( double ) ndig; rbot = 10.0 * rbot; } else if ( ihave == 8 ) { jtop = 10 * jtop + ndig; } else if ( ihave == 10 ) { jtop = 10 * jtop + ndig; jbot = 10 * jbot; } } /* Anything else is regarded as a terminator. */ else { iterm = 1; } /* If we haven't seen a terminator, and we haven't examined the entire string, go get the next character. */ if ( iterm == 1 || nchar <= *lchar + 1 ) { break; } } /* If we haven't seen a terminator, and we have examined the entire string, then we're done, and LCHAR is equal to NCHAR. */ if ( iterm != 1 && (*lchar) + 1 == nchar ) { *lchar = nchar; } /* Number seems to have terminated. Have we got a legal number? Not if we terminated in states 1, 2, 6 or 7! */ if ( ihave == 1 || ihave == 2 || ihave == 6 || ihave == 7 ) { *error = 1; return r; } /* Number seems OK. Form it. We have had some trouble with input of the form 1.0E-312. For now, let's assume anything less than 1.0E-20 is zero. */ if ( jtop == 0 ) { rexp = 1.0; } else { if ( jbot == 1 ) { if ( jsgn * jtop < -20 ) { rexp = 0.0; } else { rexp = pow ( ( double ) 10.0, ( double ) ( jsgn * jtop ) ); } } else { if ( jsgn * jtop < -20 * jbot ) { rexp = 0.0; } else { rexp = jsgn * jtop; rexp = rexp / jbot; rexp = pow ( ( double ) 10.0, ( double ) rexp ); } } } r = isgn * rexp * rtop / rbot; return r; } /******************************************************************************/ int s_to_r8vec ( char *s, int n, double rvec[] ) /******************************************************************************/ /* Purpose: S_TO_R8VEC reads an R8VEC from a string. Licensing: This code is distributed under the MIT license. Modified: 19 February 2001 Author: John Burkardt Parameters: Input, char *S, the string to be read. Input, int N, the number of values expected. Output, double RVEC[N], the values read from the string. Output, int S_TO_R8VEC, is TRUE (1) if an error occurred and FALSE (0) otherwise. */ { int error; int i; int lchar; error = 0; for ( i = 0; i < n; i++ ) { rvec[i] = s_to_r8 ( s, &lchar, &error ); if ( error ) { return error; } s = s + lchar; } return error; } /******************************************************************************/ int s_word_count ( char *s ) /******************************************************************************/ /* Purpose: S_WORD_COUNT counts the number of "words" in a string. Licensing: This code is distributed under the MIT license. Modified: 16 September 2015 Author: John Burkardt Parameters: Input, char *S, the string to be examined. Output, int S_WORD_COUNT, the number of "words" in the string. Words are presumed to be separated by one or more blanks. */ { int blank; int word_num; char *t; word_num = 0; blank = 1; t = s; while ( *t ) { if ( *t == ' ' || *t == '\n' ) { blank = 1; } else if ( blank ) { word_num = word_num + 1; blank = 0; } t++; } return word_num; } /******************************************************************************/ void timestamp ( ) /******************************************************************************/ /* 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 }