00001 /* slasdq.f -- translated by f2c (version 20061008). 00002 You must link the resulting object file with libf2c: 00003 on Microsoft Windows system, link with libf2c.lib; 00004 on Linux or Unix systems, link with .../path/to/libf2c.a -lm 00005 or, if you install libf2c.a in a standard place, with -lf2c -lm 00006 -- in that order, at the end of the command line, as in 00007 cc *.o -lf2c -lm 00008 Source for libf2c is in /netlib/f2c/libf2c.zip, e.g., 00009 00010 http://www.netlib.org/f2c/libf2c.zip 00011 */ 00012 00013 #include "f2c.h" 00014 #include "blaswrap.h" 00015 00016 /* Table of constant values */ 00017 00018 static integer c__1 = 1; 00019 00020 /* Subroutine */ int slasdq_(char *uplo, integer *sqre, integer *n, integer * 00021 ncvt, integer *nru, integer *ncc, real *d__, real *e, real *vt, 00022 integer *ldvt, real *u, integer *ldu, real *c__, integer *ldc, real * 00023 work, integer *info) 00024 { 00025 /* System generated locals */ 00026 integer c_dim1, c_offset, u_dim1, u_offset, vt_dim1, vt_offset, i__1, 00027 i__2; 00028 00029 /* Local variables */ 00030 integer i__, j; 00031 real r__, cs, sn; 00032 integer np1, isub; 00033 real smin; 00034 integer sqre1; 00035 extern logical lsame_(char *, char *); 00036 extern /* Subroutine */ int slasr_(char *, char *, char *, integer *, 00037 integer *, real *, real *, real *, integer *); 00038 integer iuplo; 00039 extern /* Subroutine */ int sswap_(integer *, real *, integer *, real *, 00040 integer *), xerbla_(char *, integer *), slartg_(real *, 00041 real *, real *, real *, real *); 00042 logical rotate; 00043 extern /* Subroutine */ int sbdsqr_(char *, integer *, integer *, integer 00044 *, integer *, real *, real *, real *, integer *, real *, integer * 00045 , real *, integer *, real *, integer *); 00046 00047 00048 /* -- LAPACK auxiliary routine (version 3.2) -- */ 00049 /* Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */ 00050 /* November 2006 */ 00051 00052 /* .. Scalar Arguments .. */ 00053 /* .. */ 00054 /* .. Array Arguments .. */ 00055 /* .. */ 00056 00057 /* Purpose */ 00058 /* ======= */ 00059 00060 /* SLASDQ computes the singular value decomposition (SVD) of a real */ 00061 /* (upper or lower) bidiagonal matrix with diagonal D and offdiagonal */ 00062 /* E, accumulating the transformations if desired. Letting B denote */ 00063 /* the input bidiagonal matrix, the algorithm computes orthogonal */ 00064 /* matrices Q and P such that B = Q * S * P' (P' denotes the transpose */ 00065 /* of P). The singular values S are overwritten on D. */ 00066 00067 /* The input matrix U is changed to U * Q if desired. */ 00068 /* The input matrix VT is changed to P' * VT if desired. */ 00069 /* The input matrix C is changed to Q' * C if desired. */ 00070 00071 /* See "Computing Small Singular Values of Bidiagonal Matrices With */ 00072 /* Guaranteed High Relative Accuracy," by J. Demmel and W. Kahan, */ 00073 /* LAPACK Working Note #3, for a detailed description of the algorithm. */ 00074 00075 /* Arguments */ 00076 /* ========= */ 00077 00078 /* UPLO (input) CHARACTER*1 */ 00079 /* On entry, UPLO specifies whether the input bidiagonal matrix */ 00080 /* is upper or lower bidiagonal, and wether it is square are */ 00081 /* not. */ 00082 /* UPLO = 'U' or 'u' B is upper bidiagonal. */ 00083 /* UPLO = 'L' or 'l' B is lower bidiagonal. */ 00084 00085 /* SQRE (input) INTEGER */ 00086 /* = 0: then the input matrix is N-by-N. */ 00087 /* = 1: then the input matrix is N-by-(N+1) if UPLU = 'U' and */ 00088 /* (N+1)-by-N if UPLU = 'L'. */ 00089 00090 /* The bidiagonal matrix has */ 00091 /* N = NL + NR + 1 rows and */ 00092 /* M = N + SQRE >= N columns. */ 00093 00094 /* N (input) INTEGER */ 00095 /* On entry, N specifies the number of rows and columns */ 00096 /* in the matrix. N must be at least 0. */ 00097 00098 /* NCVT (input) INTEGER */ 00099 /* On entry, NCVT specifies the number of columns of */ 00100 /* the matrix VT. NCVT must be at least 0. */ 00101 00102 /* NRU (input) INTEGER */ 00103 /* On entry, NRU specifies the number of rows of */ 00104 /* the matrix U. NRU must be at least 0. */ 00105 00106 /* NCC (input) INTEGER */ 00107 /* On entry, NCC specifies the number of columns of */ 00108 /* the matrix C. NCC must be at least 0. */ 00109 00110 /* D (input/output) REAL array, dimension (N) */ 00111 /* On entry, D contains the diagonal entries of the */ 00112 /* bidiagonal matrix whose SVD is desired. On normal exit, */ 00113 /* D contains the singular values in ascending order. */ 00114 00115 /* E (input/output) REAL array. */ 00116 /* dimension is (N-1) if SQRE = 0 and N if SQRE = 1. */ 00117 /* On entry, the entries of E contain the offdiagonal entries */ 00118 /* of the bidiagonal matrix whose SVD is desired. On normal */ 00119 /* exit, E will contain 0. If the algorithm does not converge, */ 00120 /* D and E will contain the diagonal and superdiagonal entries */ 00121 /* of a bidiagonal matrix orthogonally equivalent to the one */ 00122 /* given as input. */ 00123 00124 /* VT (input/output) REAL array, dimension (LDVT, NCVT) */ 00125 /* On entry, contains a matrix which on exit has been */ 00126 /* premultiplied by P', dimension N-by-NCVT if SQRE = 0 */ 00127 /* and (N+1)-by-NCVT if SQRE = 1 (not referenced if NCVT=0). */ 00128 00129 /* LDVT (input) INTEGER */ 00130 /* On entry, LDVT specifies the leading dimension of VT as */ 00131 /* declared in the calling (sub) program. LDVT must be at */ 00132 /* least 1. If NCVT is nonzero LDVT must also be at least N. */ 00133 00134 /* U (input/output) REAL array, dimension (LDU, N) */ 00135 /* On entry, contains a matrix which on exit has been */ 00136 /* postmultiplied by Q, dimension NRU-by-N if SQRE = 0 */ 00137 /* and NRU-by-(N+1) if SQRE = 1 (not referenced if NRU=0). */ 00138 00139 /* LDU (input) INTEGER */ 00140 /* On entry, LDU specifies the leading dimension of U as */ 00141 /* declared in the calling (sub) program. LDU must be at */ 00142 /* least max( 1, NRU ) . */ 00143 00144 /* C (input/output) REAL array, dimension (LDC, NCC) */ 00145 /* On entry, contains an N-by-NCC matrix which on exit */ 00146 /* has been premultiplied by Q' dimension N-by-NCC if SQRE = 0 */ 00147 /* and (N+1)-by-NCC if SQRE = 1 (not referenced if NCC=0). */ 00148 00149 /* LDC (input) INTEGER */ 00150 /* On entry, LDC specifies the leading dimension of C as */ 00151 /* declared in the calling (sub) program. LDC must be at */ 00152 /* least 1. If NCC is nonzero, LDC must also be at least N. */ 00153 00154 /* WORK (workspace) REAL array, dimension (4*N) */ 00155 /* Workspace. Only referenced if one of NCVT, NRU, or NCC is */ 00156 /* nonzero, and if N is at least 2. */ 00157 00158 /* INFO (output) INTEGER */ 00159 /* On exit, a value of 0 indicates a successful exit. */ 00160 /* If INFO < 0, argument number -INFO is illegal. */ 00161 /* If INFO > 0, the algorithm did not converge, and INFO */ 00162 /* specifies how many superdiagonals did not converge. */ 00163 00164 /* Further Details */ 00165 /* =============== */ 00166 00167 /* Based on contributions by */ 00168 /* Ming Gu and Huan Ren, Computer Science Division, University of */ 00169 /* California at Berkeley, USA */ 00170 00171 /* ===================================================================== */ 00172 00173 /* .. Parameters .. */ 00174 /* .. */ 00175 /* .. Local Scalars .. */ 00176 /* .. */ 00177 /* .. External Subroutines .. */ 00178 /* .. */ 00179 /* .. External Functions .. */ 00180 /* .. */ 00181 /* .. Intrinsic Functions .. */ 00182 /* .. */ 00183 /* .. Executable Statements .. */ 00184 00185 /* Test the input parameters. */ 00186 00187 /* Parameter adjustments */ 00188 --d__; 00189 --e; 00190 vt_dim1 = *ldvt; 00191 vt_offset = 1 + vt_dim1; 00192 vt -= vt_offset; 00193 u_dim1 = *ldu; 00194 u_offset = 1 + u_dim1; 00195 u -= u_offset; 00196 c_dim1 = *ldc; 00197 c_offset = 1 + c_dim1; 00198 c__ -= c_offset; 00199 --work; 00200 00201 /* Function Body */ 00202 *info = 0; 00203 iuplo = 0; 00204 if (lsame_(uplo, "U")) { 00205 iuplo = 1; 00206 } 00207 if (lsame_(uplo, "L")) { 00208 iuplo = 2; 00209 } 00210 if (iuplo == 0) { 00211 *info = -1; 00212 } else if (*sqre < 0 || *sqre > 1) { 00213 *info = -2; 00214 } else if (*n < 0) { 00215 *info = -3; 00216 } else if (*ncvt < 0) { 00217 *info = -4; 00218 } else if (*nru < 0) { 00219 *info = -5; 00220 } else if (*ncc < 0) { 00221 *info = -6; 00222 } else if (*ncvt == 0 && *ldvt < 1 || *ncvt > 0 && *ldvt < max(1,*n)) { 00223 *info = -10; 00224 } else if (*ldu < max(1,*nru)) { 00225 *info = -12; 00226 } else if (*ncc == 0 && *ldc < 1 || *ncc > 0 && *ldc < max(1,*n)) { 00227 *info = -14; 00228 } 00229 if (*info != 0) { 00230 i__1 = -(*info); 00231 xerbla_("SLASDQ", &i__1); 00232 return 0; 00233 } 00234 if (*n == 0) { 00235 return 0; 00236 } 00237 00238 /* ROTATE is true if any singular vectors desired, false otherwise */ 00239 00240 rotate = *ncvt > 0 || *nru > 0 || *ncc > 0; 00241 np1 = *n + 1; 00242 sqre1 = *sqre; 00243 00244 /* If matrix non-square upper bidiagonal, rotate to be lower */ 00245 /* bidiagonal. The rotations are on the right. */ 00246 00247 if (iuplo == 1 && sqre1 == 1) { 00248 i__1 = *n - 1; 00249 for (i__ = 1; i__ <= i__1; ++i__) { 00250 slartg_(&d__[i__], &e[i__], &cs, &sn, &r__); 00251 d__[i__] = r__; 00252 e[i__] = sn * d__[i__ + 1]; 00253 d__[i__ + 1] = cs * d__[i__ + 1]; 00254 if (rotate) { 00255 work[i__] = cs; 00256 work[*n + i__] = sn; 00257 } 00258 /* L10: */ 00259 } 00260 slartg_(&d__[*n], &e[*n], &cs, &sn, &r__); 00261 d__[*n] = r__; 00262 e[*n] = 0.f; 00263 if (rotate) { 00264 work[*n] = cs; 00265 work[*n + *n] = sn; 00266 } 00267 iuplo = 2; 00268 sqre1 = 0; 00269 00270 /* Update singular vectors if desired. */ 00271 00272 if (*ncvt > 0) { 00273 slasr_("L", "V", "F", &np1, ncvt, &work[1], &work[np1], &vt[ 00274 vt_offset], ldvt); 00275 } 00276 } 00277 00278 /* If matrix lower bidiagonal, rotate to be upper bidiagonal */ 00279 /* by applying Givens rotations on the left. */ 00280 00281 if (iuplo == 2) { 00282 i__1 = *n - 1; 00283 for (i__ = 1; i__ <= i__1; ++i__) { 00284 slartg_(&d__[i__], &e[i__], &cs, &sn, &r__); 00285 d__[i__] = r__; 00286 e[i__] = sn * d__[i__ + 1]; 00287 d__[i__ + 1] = cs * d__[i__ + 1]; 00288 if (rotate) { 00289 work[i__] = cs; 00290 work[*n + i__] = sn; 00291 } 00292 /* L20: */ 00293 } 00294 00295 /* If matrix (N+1)-by-N lower bidiagonal, one additional */ 00296 /* rotation is needed. */ 00297 00298 if (sqre1 == 1) { 00299 slartg_(&d__[*n], &e[*n], &cs, &sn, &r__); 00300 d__[*n] = r__; 00301 if (rotate) { 00302 work[*n] = cs; 00303 work[*n + *n] = sn; 00304 } 00305 } 00306 00307 /* Update singular vectors if desired. */ 00308 00309 if (*nru > 0) { 00310 if (sqre1 == 0) { 00311 slasr_("R", "V", "F", nru, n, &work[1], &work[np1], &u[ 00312 u_offset], ldu); 00313 } else { 00314 slasr_("R", "V", "F", nru, &np1, &work[1], &work[np1], &u[ 00315 u_offset], ldu); 00316 } 00317 } 00318 if (*ncc > 0) { 00319 if (sqre1 == 0) { 00320 slasr_("L", "V", "F", n, ncc, &work[1], &work[np1], &c__[ 00321 c_offset], ldc); 00322 } else { 00323 slasr_("L", "V", "F", &np1, ncc, &work[1], &work[np1], &c__[ 00324 c_offset], ldc); 00325 } 00326 } 00327 } 00328 00329 /* Call SBDSQR to compute the SVD of the reduced real */ 00330 /* N-by-N upper bidiagonal matrix. */ 00331 00332 sbdsqr_("U", n, ncvt, nru, ncc, &d__[1], &e[1], &vt[vt_offset], ldvt, &u[ 00333 u_offset], ldu, &c__[c_offset], ldc, &work[1], info); 00334 00335 /* Sort the singular values into ascending order (insertion sort on */ 00336 /* singular values, but only one transposition per singular vector) */ 00337 00338 i__1 = *n; 00339 for (i__ = 1; i__ <= i__1; ++i__) { 00340 00341 /* Scan for smallest D(I). */ 00342 00343 isub = i__; 00344 smin = d__[i__]; 00345 i__2 = *n; 00346 for (j = i__ + 1; j <= i__2; ++j) { 00347 if (d__[j] < smin) { 00348 isub = j; 00349 smin = d__[j]; 00350 } 00351 /* L30: */ 00352 } 00353 if (isub != i__) { 00354 00355 /* Swap singular values and vectors. */ 00356 00357 d__[isub] = d__[i__]; 00358 d__[i__] = smin; 00359 if (*ncvt > 0) { 00360 sswap_(ncvt, &vt[isub + vt_dim1], ldvt, &vt[i__ + vt_dim1], 00361 ldvt); 00362 } 00363 if (*nru > 0) { 00364 sswap_(nru, &u[isub * u_dim1 + 1], &c__1, &u[i__ * u_dim1 + 1] 00365 , &c__1); 00366 } 00367 if (*ncc > 0) { 00368 sswap_(ncc, &c__[isub + c_dim1], ldc, &c__[i__ + c_dim1], ldc) 00369 ; 00370 } 00371 } 00372 /* L40: */ 00373 } 00374 00375 return 0; 00376 00377 /* End of SLASDQ */ 00378 00379 } /* slasdq_ */