3 /* Table of constant values */
5 static integer c__1 = 1;
7 /* Subroutine */ int dlasdq_(char *uplo, integer *sqre, integer *n, integer *
8 ncvt, integer *nru, integer *ncc, doublereal *d__, doublereal *e,
9 doublereal *vt, integer *ldvt, doublereal *u, integer *ldu,
10 doublereal *c__, integer *ldc, doublereal *work, integer *info)
12 /* System generated locals */
13 integer c_dim1, c_offset, u_dim1, u_offset, vt_dim1, vt_offset, i__1,
18 doublereal r__, cs, sn;
22 extern logical lsame_(char *, char *);
23 extern /* Subroutine */ int dlasr_(char *, char *, char *, integer *,
24 integer *, doublereal *, doublereal *, doublereal *, integer *), dswap_(integer *, doublereal *, integer *
25 , doublereal *, integer *);
27 extern /* Subroutine */ int dlartg_(doublereal *, doublereal *,
28 doublereal *, doublereal *, doublereal *), xerbla_(char *,
29 integer *), dbdsqr_(char *, integer *, integer *, integer
30 *, integer *, doublereal *, doublereal *, doublereal *, integer *,
31 doublereal *, integer *, doublereal *, integer *, doublereal *,
36 /* -- LAPACK auxiliary routine (version 3.1) -- */
37 /* Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */
40 /* .. Scalar Arguments .. */
42 /* .. Array Arguments .. */
48 /* DLASDQ computes the singular value decomposition (SVD) of a real */
49 /* (upper or lower) bidiagonal matrix with diagonal D and offdiagonal */
50 /* E, accumulating the transformations if desired. Letting B denote */
51 /* the input bidiagonal matrix, the algorithm computes orthogonal */
52 /* matrices Q and P such that B = Q * S * P' (P' denotes the transpose */
53 /* of P). The singular values S are overwritten on D. */
55 /* The input matrix U is changed to U * Q if desired. */
56 /* The input matrix VT is changed to P' * VT if desired. */
57 /* The input matrix C is changed to Q' * C if desired. */
59 /* See "Computing Small Singular Values of Bidiagonal Matrices With */
60 /* Guaranteed High Relative Accuracy," by J. Demmel and W. Kahan, */
61 /* LAPACK Working Note #3, for a detailed description of the algorithm. */
66 /* UPLO (input) CHARACTER*1 */
67 /* On entry, UPLO specifies whether the input bidiagonal matrix */
68 /* is upper or lower bidiagonal, and wether it is square are */
70 /* UPLO = 'U' or 'u' B is upper bidiagonal. */
71 /* UPLO = 'L' or 'l' B is lower bidiagonal. */
73 /* SQRE (input) INTEGER */
74 /* = 0: then the input matrix is N-by-N. */
75 /* = 1: then the input matrix is N-by-(N+1) if UPLU = 'U' and */
76 /* (N+1)-by-N if UPLU = 'L'. */
78 /* The bidiagonal matrix has */
79 /* N = NL + NR + 1 rows and */
80 /* M = N + SQRE >= N columns. */
82 /* N (input) INTEGER */
83 /* On entry, N specifies the number of rows and columns */
84 /* in the matrix. N must be at least 0. */
86 /* NCVT (input) INTEGER */
87 /* On entry, NCVT specifies the number of columns of */
88 /* the matrix VT. NCVT must be at least 0. */
90 /* NRU (input) INTEGER */
91 /* On entry, NRU specifies the number of rows of */
92 /* the matrix U. NRU must be at least 0. */
94 /* NCC (input) INTEGER */
95 /* On entry, NCC specifies the number of columns of */
96 /* the matrix C. NCC must be at least 0. */
98 /* D (input/output) DOUBLE PRECISION array, dimension (N) */
99 /* On entry, D contains the diagonal entries of the */
100 /* bidiagonal matrix whose SVD is desired. On normal exit, */
101 /* D contains the singular values in ascending order. */
103 /* E (input/output) DOUBLE PRECISION array. */
104 /* dimension is (N-1) if SQRE = 0 and N if SQRE = 1. */
105 /* On entry, the entries of E contain the offdiagonal entries */
106 /* of the bidiagonal matrix whose SVD is desired. On normal */
107 /* exit, E will contain 0. If the algorithm does not converge, */
108 /* D and E will contain the diagonal and superdiagonal entries */
109 /* of a bidiagonal matrix orthogonally equivalent to the one */
110 /* given as input. */
112 /* VT (input/output) DOUBLE PRECISION array, dimension (LDVT, NCVT) */
113 /* On entry, contains a matrix which on exit has been */
114 /* premultiplied by P', dimension N-by-NCVT if SQRE = 0 */
115 /* and (N+1)-by-NCVT if SQRE = 1 (not referenced if NCVT=0). */
117 /* LDVT (input) INTEGER */
118 /* On entry, LDVT specifies the leading dimension of VT as */
119 /* declared in the calling (sub) program. LDVT must be at */
120 /* least 1. If NCVT is nonzero LDVT must also be at least N. */
122 /* U (input/output) DOUBLE PRECISION array, dimension (LDU, N) */
123 /* On entry, contains a matrix which on exit has been */
124 /* postmultiplied by Q, dimension NRU-by-N if SQRE = 0 */
125 /* and NRU-by-(N+1) if SQRE = 1 (not referenced if NRU=0). */
127 /* LDU (input) INTEGER */
128 /* On entry, LDU specifies the leading dimension of U as */
129 /* declared in the calling (sub) program. LDU must be at */
130 /* least max( 1, NRU ) . */
132 /* C (input/output) DOUBLE PRECISION array, dimension (LDC, NCC) */
133 /* On entry, contains an N-by-NCC matrix which on exit */
134 /* has been premultiplied by Q' dimension N-by-NCC if SQRE = 0 */
135 /* and (N+1)-by-NCC if SQRE = 1 (not referenced if NCC=0). */
137 /* LDC (input) INTEGER */
138 /* On entry, LDC specifies the leading dimension of C as */
139 /* declared in the calling (sub) program. LDC must be at */
140 /* least 1. If NCC is nonzero, LDC must also be at least N. */
142 /* WORK (workspace) DOUBLE PRECISION array, dimension (4*N) */
143 /* Workspace. Only referenced if one of NCVT, NRU, or NCC is */
144 /* nonzero, and if N is at least 2. */
146 /* INFO (output) INTEGER */
147 /* On exit, a value of 0 indicates a successful exit. */
148 /* If INFO < 0, argument number -INFO is illegal. */
149 /* If INFO > 0, the algorithm did not converge, and INFO */
150 /* specifies how many superdiagonals did not converge. */
152 /* Further Details */
153 /* =============== */
155 /* Based on contributions by */
156 /* Ming Gu and Huan Ren, Computer Science Division, University of */
157 /* California at Berkeley, USA */
159 /* ===================================================================== */
161 /* .. Parameters .. */
163 /* .. Local Scalars .. */
165 /* .. External Subroutines .. */
167 /* .. External Functions .. */
169 /* .. Intrinsic Functions .. */
171 /* .. Executable Statements .. */
173 /* Test the input parameters. */
175 /* Parameter adjustments */
179 vt_offset = 1 + vt_dim1;
182 u_offset = 1 + u_dim1;
185 c_offset = 1 + c_dim1;
192 if (lsame_(uplo, "U")) {
195 if (lsame_(uplo, "L")) {
200 } else if (*sqre < 0 || *sqre > 1) {
204 } else if (*ncvt < 0) {
206 } else if (*nru < 0) {
208 } else if (*ncc < 0) {
210 } else if (*ncvt == 0 && *ldvt < 1 || *ncvt > 0 && *ldvt < max(1,*n)) {
212 } else if (*ldu < max(1,*nru)) {
214 } else if (*ncc == 0 && *ldc < 1 || *ncc > 0 && *ldc < max(1,*n)) {
219 xerbla_("DLASDQ", &i__1);
226 /* ROTATE is true if any singular vectors desired, false otherwise */
228 rotate = *ncvt > 0 || *nru > 0 || *ncc > 0;
232 /* If matrix non-square upper bidiagonal, rotate to be lower */
233 /* bidiagonal. The rotations are on the right. */
235 if (iuplo == 1 && sqre1 == 1) {
237 for (i__ = 1; i__ <= i__1; ++i__) {
238 dlartg_(&d__[i__], &e[i__], &cs, &sn, &r__);
240 e[i__] = sn * d__[i__ + 1];
241 d__[i__ + 1] = cs * d__[i__ + 1];
248 dlartg_(&d__[*n], &e[*n], &cs, &sn, &r__);
258 /* Update singular vectors if desired. */
261 dlasr_("L", "V", "F", &np1, ncvt, &work[1], &work[np1], &vt[
266 /* If matrix lower bidiagonal, rotate to be upper bidiagonal */
267 /* by applying Givens rotations on the left. */
271 for (i__ = 1; i__ <= i__1; ++i__) {
272 dlartg_(&d__[i__], &e[i__], &cs, &sn, &r__);
274 e[i__] = sn * d__[i__ + 1];
275 d__[i__ + 1] = cs * d__[i__ + 1];
283 /* If matrix (N+1)-by-N lower bidiagonal, one additional */
284 /* rotation is needed. */
287 dlartg_(&d__[*n], &e[*n], &cs, &sn, &r__);
295 /* Update singular vectors if desired. */
299 dlasr_("R", "V", "F", nru, n, &work[1], &work[np1], &u[
302 dlasr_("R", "V", "F", nru, &np1, &work[1], &work[np1], &u[
308 dlasr_("L", "V", "F", n, ncc, &work[1], &work[np1], &c__[
311 dlasr_("L", "V", "F", &np1, ncc, &work[1], &work[np1], &c__[
317 /* Call DBDSQR to compute the SVD of the reduced real */
318 /* N-by-N upper bidiagonal matrix. */
320 dbdsqr_("U", n, ncvt, nru, ncc, &d__[1], &e[1], &vt[vt_offset], ldvt, &u[
321 u_offset], ldu, &c__[c_offset], ldc, &work[1], info);
323 /* Sort the singular values into ascending order (insertion sort on */
324 /* singular values, but only one transposition per singular vector) */
327 for (i__ = 1; i__ <= i__1; ++i__) {
329 /* Scan for smallest D(I). */
334 for (j = i__ + 1; j <= i__2; ++j) {
343 /* Swap singular values and vectors. */
345 d__[isub] = d__[i__];
348 dswap_(ncvt, &vt[isub + vt_dim1], ldvt, &vt[i__ + vt_dim1],
352 dswap_(nru, &u[isub * u_dim1 + 1], &c__1, &u[i__ * u_dim1 + 1]
356 dswap_(ncc, &c__[isub + c_dim1], ldc, &c__[i__ + c_dim1], ldc)