3 /* Table of constant values */
5 static integer c__1 = 1;
7 /* Subroutine */ int slasdq_(char *uplo, integer *sqre, integer *n, integer *
8 ncvt, integer *nru, integer *ncc, real *d__, real *e, real *vt,
9 integer *ldvt, real *u, integer *ldu, real *c__, integer *ldc, real *
12 /* System generated locals */
13 integer c_dim1, c_offset, u_dim1, u_offset, vt_dim1, vt_offset, i__1,
22 extern logical lsame_(char *, char *);
23 extern /* Subroutine */ int slasr_(char *, char *, char *, integer *,
24 integer *, real *, real *, real *, integer *);
26 extern /* Subroutine */ int sswap_(integer *, real *, integer *, real *,
27 integer *), xerbla_(char *, integer *), slartg_(real *,
28 real *, real *, real *, real *);
30 extern /* Subroutine */ int sbdsqr_(char *, integer *, integer *, integer
31 *, integer *, real *, real *, real *, integer *, real *, integer *
32 , real *, integer *, real *, integer *);
35 /* -- LAPACK auxiliary routine (version 3.1) -- */
36 /* Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */
39 /* .. Scalar Arguments .. */
41 /* .. Array Arguments .. */
47 /* SLASDQ computes the singular value decomposition (SVD) of a real */
48 /* (upper or lower) bidiagonal matrix with diagonal D and offdiagonal */
49 /* E, accumulating the transformations if desired. Letting B denote */
50 /* the input bidiagonal matrix, the algorithm computes orthogonal */
51 /* matrices Q and P such that B = Q * S * P' (P' denotes the transpose */
52 /* of P). The singular values S are overwritten on D. */
54 /* The input matrix U is changed to U * Q if desired. */
55 /* The input matrix VT is changed to P' * VT if desired. */
56 /* The input matrix C is changed to Q' * C if desired. */
58 /* See "Computing Small Singular Values of Bidiagonal Matrices With */
59 /* Guaranteed High Relative Accuracy," by J. Demmel and W. Kahan, */
60 /* LAPACK Working Note #3, for a detailed description of the algorithm. */
65 /* UPLO (input) CHARACTER*1 */
66 /* On entry, UPLO specifies whether the input bidiagonal matrix */
67 /* is upper or lower bidiagonal, and wether it is square are */
69 /* UPLO = 'U' or 'u' B is upper bidiagonal. */
70 /* UPLO = 'L' or 'l' B is lower bidiagonal. */
72 /* SQRE (input) INTEGER */
73 /* = 0: then the input matrix is N-by-N. */
74 /* = 1: then the input matrix is N-by-(N+1) if UPLU = 'U' and */
75 /* (N+1)-by-N if UPLU = 'L'. */
77 /* The bidiagonal matrix has */
78 /* N = NL + NR + 1 rows and */
79 /* M = N + SQRE >= N columns. */
81 /* N (input) INTEGER */
82 /* On entry, N specifies the number of rows and columns */
83 /* in the matrix. N must be at least 0. */
85 /* NCVT (input) INTEGER */
86 /* On entry, NCVT specifies the number of columns of */
87 /* the matrix VT. NCVT must be at least 0. */
89 /* NRU (input) INTEGER */
90 /* On entry, NRU specifies the number of rows of */
91 /* the matrix U. NRU must be at least 0. */
93 /* NCC (input) INTEGER */
94 /* On entry, NCC specifies the number of columns of */
95 /* the matrix C. NCC must be at least 0. */
97 /* D (input/output) REAL array, dimension (N) */
98 /* On entry, D contains the diagonal entries of the */
99 /* bidiagonal matrix whose SVD is desired. On normal exit, */
100 /* D contains the singular values in ascending order. */
102 /* E (input/output) REAL array. */
103 /* dimension is (N-1) if SQRE = 0 and N if SQRE = 1. */
104 /* On entry, the entries of E contain the offdiagonal entries */
105 /* of the bidiagonal matrix whose SVD is desired. On normal */
106 /* exit, E will contain 0. If the algorithm does not converge, */
107 /* D and E will contain the diagonal and superdiagonal entries */
108 /* of a bidiagonal matrix orthogonally equivalent to the one */
109 /* given as input. */
111 /* VT (input/output) REAL array, dimension (LDVT, NCVT) */
112 /* On entry, contains a matrix which on exit has been */
113 /* premultiplied by P', dimension N-by-NCVT if SQRE = 0 */
114 /* and (N+1)-by-NCVT if SQRE = 1 (not referenced if NCVT=0). */
116 /* LDVT (input) INTEGER */
117 /* On entry, LDVT specifies the leading dimension of VT as */
118 /* declared in the calling (sub) program. LDVT must be at */
119 /* least 1. If NCVT is nonzero LDVT must also be at least N. */
121 /* U (input/output) REAL array, dimension (LDU, N) */
122 /* On entry, contains a matrix which on exit has been */
123 /* postmultiplied by Q, dimension NRU-by-N if SQRE = 0 */
124 /* and NRU-by-(N+1) if SQRE = 1 (not referenced if NRU=0). */
126 /* LDU (input) INTEGER */
127 /* On entry, LDU specifies the leading dimension of U as */
128 /* declared in the calling (sub) program. LDU must be at */
129 /* least max( 1, NRU ) . */
131 /* C (input/output) REAL array, dimension (LDC, NCC) */
132 /* On entry, contains an N-by-NCC matrix which on exit */
133 /* has been premultiplied by Q' dimension N-by-NCC if SQRE = 0 */
134 /* and (N+1)-by-NCC if SQRE = 1 (not referenced if NCC=0). */
136 /* LDC (input) INTEGER */
137 /* On entry, LDC specifies the leading dimension of C as */
138 /* declared in the calling (sub) program. LDC must be at */
139 /* least 1. If NCC is nonzero, LDC must also be at least N. */
141 /* WORK (workspace) REAL array, dimension (4*N) */
142 /* Workspace. Only referenced if one of NCVT, NRU, or NCC is */
143 /* nonzero, and if N is at least 2. */
145 /* INFO (output) INTEGER */
146 /* On exit, a value of 0 indicates a successful exit. */
147 /* If INFO < 0, argument number -INFO is illegal. */
148 /* If INFO > 0, the algorithm did not converge, and INFO */
149 /* specifies how many superdiagonals did not converge. */
151 /* Further Details */
152 /* =============== */
154 /* Based on contributions by */
155 /* Ming Gu and Huan Ren, Computer Science Division, University of */
156 /* California at Berkeley, USA */
158 /* ===================================================================== */
160 /* .. Parameters .. */
162 /* .. Local Scalars .. */
164 /* .. External Subroutines .. */
166 /* .. External Functions .. */
168 /* .. Intrinsic Functions .. */
170 /* .. Executable Statements .. */
172 /* Test the input parameters. */
174 /* Parameter adjustments */
178 vt_offset = 1 + vt_dim1;
181 u_offset = 1 + u_dim1;
184 c_offset = 1 + c_dim1;
191 if (lsame_(uplo, "U")) {
194 if (lsame_(uplo, "L")) {
199 } else if (*sqre < 0 || *sqre > 1) {
203 } else if (*ncvt < 0) {
205 } else if (*nru < 0) {
207 } else if (*ncc < 0) {
209 } else if (*ncvt == 0 && *ldvt < 1 || *ncvt > 0 && *ldvt < max(1,*n)) {
211 } else if (*ldu < max(1,*nru)) {
213 } else if (*ncc == 0 && *ldc < 1 || *ncc > 0 && *ldc < max(1,*n)) {
218 xerbla_("SLASDQ", &i__1);
225 /* ROTATE is true if any singular vectors desired, false otherwise */
227 rotate = *ncvt > 0 || *nru > 0 || *ncc > 0;
231 /* If matrix non-square upper bidiagonal, rotate to be lower */
232 /* bidiagonal. The rotations are on the right. */
234 if (iuplo == 1 && sqre1 == 1) {
236 for (i__ = 1; i__ <= i__1; ++i__) {
237 slartg_(&d__[i__], &e[i__], &cs, &sn, &r__);
239 e[i__] = sn * d__[i__ + 1];
240 d__[i__ + 1] = cs * d__[i__ + 1];
247 slartg_(&d__[*n], &e[*n], &cs, &sn, &r__);
257 /* Update singular vectors if desired. */
260 slasr_("L", "V", "F", &np1, ncvt, &work[1], &work[np1], &vt[
265 /* If matrix lower bidiagonal, rotate to be upper bidiagonal */
266 /* by applying Givens rotations on the left. */
270 for (i__ = 1; i__ <= i__1; ++i__) {
271 slartg_(&d__[i__], &e[i__], &cs, &sn, &r__);
273 e[i__] = sn * d__[i__ + 1];
274 d__[i__ + 1] = cs * d__[i__ + 1];
282 /* If matrix (N+1)-by-N lower bidiagonal, one additional */
283 /* rotation is needed. */
286 slartg_(&d__[*n], &e[*n], &cs, &sn, &r__);
294 /* Update singular vectors if desired. */
298 slasr_("R", "V", "F", nru, n, &work[1], &work[np1], &u[
301 slasr_("R", "V", "F", nru, &np1, &work[1], &work[np1], &u[
307 slasr_("L", "V", "F", n, ncc, &work[1], &work[np1], &c__[
310 slasr_("L", "V", "F", &np1, ncc, &work[1], &work[np1], &c__[
316 /* Call SBDSQR to compute the SVD of the reduced real */
317 /* N-by-N upper bidiagonal matrix. */
319 sbdsqr_("U", n, ncvt, nru, ncc, &d__[1], &e[1], &vt[vt_offset], ldvt, &u[
320 u_offset], ldu, &c__[c_offset], ldc, &work[1], info);
322 /* Sort the singular values into ascending order (insertion sort on */
323 /* singular values, but only one transposition per singular vector) */
326 for (i__ = 1; i__ <= i__1; ++i__) {
328 /* Scan for smallest D(I). */
333 for (j = i__ + 1; j <= i__2; ++j) {
342 /* Swap singular values and vectors. */
344 d__[isub] = d__[i__];
347 sswap_(ncvt, &vt[isub + vt_dim1], ldvt, &vt[i__ + vt_dim1],
351 sswap_(nru, &u[isub * u_dim1 + 1], &c__1, &u[i__ * u_dim1 + 1]
355 sswap_(ncc, &c__[isub + c_dim1], ldc, &c__[i__ + c_dim1], ldc)