1 /* ./src_f77/dgesdd.f -- translated by f2c (version 20030320).
2 You must link the resulting object file with the libraries:
3 -lf2c -lm (in that order)
4 */
5
6 #include <punc/vf2c.h>
7
8 /* Table of constant values */
9
10 static integer c__1 = 1;
11 static integer c_n1 = -1;
12 static integer c__0 = 0;
13 static doublereal c_b227 = 0.;
14 static doublereal c_b248 = 1.;
15
dgesdd_(char * jobz,integer * m,integer * n,doublereal * a,integer * lda,doublereal * s,doublereal * u,integer * ldu,doublereal * vt,integer * ldvt,doublereal * work,integer * lwork,integer * iwork,integer * info,ftnlen jobz_len)16 /* Subroutine */ int dgesdd_(char *jobz, integer *m, integer *n, doublereal *
17 a, integer *lda, doublereal *s, doublereal *u, integer *ldu,
18 doublereal *vt, integer *ldvt, doublereal *work, integer *lwork,
19 integer *iwork, integer *info, ftnlen jobz_len)
20 {
21 /* System generated locals */
22 integer a_dim1, a_offset, u_dim1, u_offset, vt_dim1, vt_offset, i__1,
23 i__2, i__3;
24
25 /* Builtin functions */
26 double sqrt(doublereal);
27
28 /* Local variables */
29 static integer i__, ie, il, ir, iu, blk;
30 static doublereal dum[1], eps;
31 static integer ivt, iscl;
32 static doublereal anrm;
33 static integer idum[1], ierr, itau;
34 extern /* Subroutine */ int dgemm_(char *, char *, integer *, integer *,
35 integer *, doublereal *, doublereal *, integer *, doublereal *,
36 integer *, doublereal *, doublereal *, integer *, ftnlen, ftnlen);
37 extern logical lsame_(char *, char *, ftnlen, ftnlen);
38 static integer chunk, minmn, wrkbl, itaup, itauq, mnthr;
39 static logical wntqa;
40 static integer nwork;
41 static logical wntqn, wntqo, wntqs;
42 extern /* Subroutine */ int dbdsdc_(char *, char *, integer *, doublereal
43 *, doublereal *, doublereal *, integer *, doublereal *, integer *,
44 doublereal *, integer *, doublereal *, integer *, integer *,
45 ftnlen, ftnlen), dgebrd_(integer *, integer *, doublereal *,
46 integer *, doublereal *, doublereal *, doublereal *, doublereal *,
47 doublereal *, integer *, integer *);
48 extern doublereal dlamch_(char *, ftnlen), dlange_(char *, integer *,
49 integer *, doublereal *, integer *, doublereal *, ftnlen);
50 static integer bdspac;
51 extern /* Subroutine */ int dgelqf_(integer *, integer *, doublereal *,
52 integer *, doublereal *, doublereal *, integer *, integer *),
53 dlascl_(char *, integer *, integer *, doublereal *, doublereal *,
54 integer *, integer *, doublereal *, integer *, integer *, ftnlen),
55 dgeqrf_(integer *, integer *, doublereal *, integer *,
56 doublereal *, doublereal *, integer *, integer *), dlacpy_(char *,
57 integer *, integer *, doublereal *, integer *, doublereal *,
58 integer *, ftnlen), dlaset_(char *, integer *, integer *,
59 doublereal *, doublereal *, doublereal *, integer *, ftnlen),
60 xerbla_(char *, integer *, ftnlen), dorgbr_(char *, integer *,
61 integer *, integer *, doublereal *, integer *, doublereal *,
62 doublereal *, integer *, integer *, ftnlen);
63 extern integer ilaenv_(integer *, char *, char *, integer *, integer *,
64 integer *, integer *, ftnlen, ftnlen);
65 static doublereal bignum;
66 extern /* Subroutine */ int dormbr_(char *, char *, char *, integer *,
67 integer *, integer *, doublereal *, integer *, doublereal *,
68 doublereal *, integer *, doublereal *, integer *, integer *,
69 ftnlen, ftnlen, ftnlen), dorglq_(integer *, integer *, integer *,
70 doublereal *, integer *, doublereal *, doublereal *, integer *,
71 integer *), dorgqr_(integer *, integer *, integer *, doublereal *,
72 integer *, doublereal *, doublereal *, integer *, integer *);
73 static integer ldwrkl, ldwrkr, minwrk, ldwrku, maxwrk, ldwkvt;
74 static doublereal smlnum;
75 static logical wntqas, lquery;
76
77
78 /* -- LAPACK driver routine (version 3.0) -- */
79 /* Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., */
80 /* Courant Institute, Argonne National Lab, and Rice University */
81 /* October 31, 1999 */
82
83 /* .. Scalar Arguments .. */
84 /* .. */
85 /* .. Array Arguments .. */
86 /* .. */
87
88 /* Purpose */
89 /* ======= */
90
91 /* DGESDD computes the singular value decomposition (SVD) of a real */
92 /* M-by-N matrix A, optionally computing the left and right singular */
93 /* vectors. If singular vectors are desired, it uses a */
94 /* divide-and-conquer algorithm. */
95
96 /* The SVD is written */
97
98 /* A = U * SIGMA * transpose(V) */
99
100 /* where SIGMA is an M-by-N matrix which is zero except for its */
101 /* min(m,n) diagonal elements, U is an M-by-M orthogonal matrix, and */
102 /* V is an N-by-N orthogonal matrix. The diagonal elements of SIGMA */
103 /* are the singular values of A; they are real and non-negative, and */
104 /* are returned in descending order. The first min(m,n) columns of */
105 /* U and V are the left and right singular vectors of A. */
106
107 /* Note that the routine returns VT = V**T, not V. */
108
109 /* The divide and conquer algorithm makes very mild assumptions about */
110 /* floating point arithmetic. It will work on machines with a guard */
111 /* digit in add/subtract, or on those binary machines without guard */
112 /* digits which subtract like the Cray X-MP, Cray Y-MP, Cray C-90, or */
113 /* Cray-2. It could conceivably fail on hexadecimal or decimal machines */
114 /* without guard digits, but we know of none. */
115
116 /* Arguments */
117 /* ========= */
118
119 /* JOBZ (input) CHARACTER*1 */
120 /* Specifies options for computing all or part of the matrix U: */
121 /* = 'A': all M columns of U and all N rows of V**T are */
122 /* returned in the arrays U and VT; */
123 /* = 'S': the first min(M,N) columns of U and the first */
124 /* min(M,N) rows of V**T are returned in the arrays U */
125 /* and VT; */
126 /* = 'O': If M >= N, the first N columns of U are overwritten */
127 /* on the array A and all rows of V**T are returned in */
128 /* the array VT; */
129 /* otherwise, all columns of U are returned in the */
130 /* array U and the first M rows of V**T are overwritten */
131 /* in the array VT; */
132 /* = 'N': no columns of U or rows of V**T are computed. */
133
134 /* M (input) INTEGER */
135 /* The number of rows of the input matrix A. M >= 0. */
136
137 /* N (input) INTEGER */
138 /* The number of columns of the input matrix A. N >= 0. */
139
140 /* A (input/output) DOUBLE PRECISION array, dimension (LDA,N) */
141 /* On entry, the M-by-N matrix A. */
142 /* On exit, */
143 /* if JOBZ = 'O', A is overwritten with the first N columns */
144 /* of U (the left singular vectors, stored */
145 /* columnwise) if M >= N; */
146 /* A is overwritten with the first M rows */
147 /* of V**T (the right singular vectors, stored */
148 /* rowwise) otherwise. */
149 /* if JOBZ .ne. 'O', the contents of A are destroyed. */
150
151 /* LDA (input) INTEGER */
152 /* The leading dimension of the array A. LDA >= max(1,M). */
153
154 /* S (output) DOUBLE PRECISION array, dimension (min(M,N)) */
155 /* The singular values of A, sorted so that S(i) >= S(i+1). */
156
157 /* U (output) DOUBLE PRECISION array, dimension (LDU,UCOL) */
158 /* UCOL = M if JOBZ = 'A' or JOBZ = 'O' and M < N; */
159 /* UCOL = min(M,N) if JOBZ = 'S'. */
160 /* If JOBZ = 'A' or JOBZ = 'O' and M < N, U contains the M-by-M */
161 /* orthogonal matrix U; */
162 /* if JOBZ = 'S', U contains the first min(M,N) columns of U */
163 /* (the left singular vectors, stored columnwise); */
164 /* if JOBZ = 'O' and M >= N, or JOBZ = 'N', U is not referenced. */
165
166 /* LDU (input) INTEGER */
167 /* The leading dimension of the array U. LDU >= 1; if */
168 /* JOBZ = 'S' or 'A' or JOBZ = 'O' and M < N, LDU >= M. */
169
170 /* VT (output) DOUBLE PRECISION array, dimension (LDVT,N) */
171 /* If JOBZ = 'A' or JOBZ = 'O' and M >= N, VT contains the */
172 /* N-by-N orthogonal matrix V**T; */
173 /* if JOBZ = 'S', VT contains the first min(M,N) rows of */
174 /* V**T (the right singular vectors, stored rowwise); */
175 /* if JOBZ = 'O' and M < N, or JOBZ = 'N', VT is not referenced. */
176
177 /* LDVT (input) INTEGER */
178 /* The leading dimension of the array VT. LDVT >= 1; if */
179 /* JOBZ = 'A' or JOBZ = 'O' and M >= N, LDVT >= N; */
180 /* if JOBZ = 'S', LDVT >= min(M,N). */
181
182 /* WORK (workspace/output) DOUBLE PRECISION array, dimension (LWORK) */
183 /* On exit, if INFO = 0, WORK(1) returns the optimal LWORK; */
184
185 /* LWORK (input) INTEGER */
186 /* The dimension of the array WORK. LWORK >= 1. */
187 /* If JOBZ = 'N', */
188 /* LWORK >= 3*min(M,N) + max(max(M,N),6*min(M,N)). */
189 /* If JOBZ = 'O', */
190 /* LWORK >= 3*min(M,N)*min(M,N) + */
191 /* max(max(M,N),5*min(M,N)*min(M,N)+4*min(M,N)). */
192 /* If JOBZ = 'S' or 'A' */
193 /* LWORK >= 3*min(M,N)*min(M,N) + */
194 /* max(max(M,N),4*min(M,N)*min(M,N)+4*min(M,N)). */
195 /* For good performance, LWORK should generally be larger. */
196 /* If LWORK < 0 but other input arguments are legal, WORK(1) */
197 /* returns the optimal LWORK. */
198
199 /* IWORK (workspace) INTEGER array, dimension (8*min(M,N)) */
200
201 /* INFO (output) INTEGER */
202 /* = 0: successful exit. */
203 /* < 0: if INFO = -i, the i-th argument had an illegal value. */
204 /* > 0: DBDSDC did not converge, updating process failed. */
205
206 /* Further Details */
207 /* =============== */
208
209 /* Based on contributions by */
210 /* Ming Gu and Huan Ren, Computer Science Division, University of */
211 /* California at Berkeley, USA */
212
213 /* ===================================================================== */
214
215 /* .. Parameters .. */
216 /* .. */
217 /* .. Local Scalars .. */
218 /* .. */
219 /* .. Local Arrays .. */
220 /* .. */
221 /* .. External Subroutines .. */
222 /* .. */
223 /* .. External Functions .. */
224 /* .. */
225 /* .. Intrinsic Functions .. */
226 /* .. */
227 /* .. Executable Statements .. */
228
229 /* Test the input arguments */
230
231 /* Parameter adjustments */
232 a_dim1 = *lda;
233 a_offset = 1 + a_dim1;
234 a -= a_offset;
235 --s;
236 u_dim1 = *ldu;
237 u_offset = 1 + u_dim1;
238 u -= u_offset;
239 vt_dim1 = *ldvt;
240 vt_offset = 1 + vt_dim1;
241 vt -= vt_offset;
242 --work;
243 --iwork;
244
245 /* Function Body */
246 *info = 0;
247 minmn = min(*m,*n);
248 mnthr = (integer) (minmn * 11. / 6.);
249 wntqa = lsame_(jobz, "A", (ftnlen)1, (ftnlen)1);
250 wntqs = lsame_(jobz, "S", (ftnlen)1, (ftnlen)1);
251 wntqas = wntqa || wntqs;
252 wntqo = lsame_(jobz, "O", (ftnlen)1, (ftnlen)1);
253 wntqn = lsame_(jobz, "N", (ftnlen)1, (ftnlen)1);
254 minwrk = 1;
255 maxwrk = 1;
256 lquery = *lwork == -1;
257
258 if (! (wntqa || wntqs || wntqo || wntqn)) {
259 *info = -1;
260 } else if (*m < 0) {
261 *info = -2;
262 } else if (*n < 0) {
263 *info = -3;
264 } else if (*lda < max(1,*m)) {
265 *info = -5;
266 } else if (*ldu < 1 || wntqas && *ldu < *m || wntqo && *m < *n && *ldu < *
267 m) {
268 *info = -8;
269 } else if (*ldvt < 1 || wntqa && *ldvt < *n || wntqs && *ldvt < minmn ||
270 wntqo && *m >= *n && *ldvt < *n) {
271 *info = -10;
272 }
273
274 /* Compute workspace */
275 /* (Note: Comments in the code beginning "Workspace:" describe the */
276 /* minimal amount of workspace needed at that point in the code, */
277 /* as well as the preferred amount for good performance. */
278 /* NB refers to the optimal block size for the immediately */
279 /* following subroutine, as returned by ILAENV.) */
280
281 if (*info == 0 && *m > 0 && *n > 0) {
282 if (*m >= *n) {
283
284 /* Compute space needed for DBDSDC */
285
286 if (wntqn) {
287 bdspac = *n * 7;
288 } else {
289 bdspac = *n * 3 * *n + (*n << 2);
290 }
291 if (*m >= mnthr) {
292 if (wntqn) {
293
294 /* Path 1 (M much larger than N, JOBZ='N') */
295
296 wrkbl = *n + *n * ilaenv_(&c__1, "DGEQRF", " ", m, n, &
297 c_n1, &c_n1, (ftnlen)6, (ftnlen)1);
298 /* Computing MAX */
299 i__1 = wrkbl, i__2 = *n * 3 + (*n << 1) * ilaenv_(&c__1,
300 "DGEBRD", " ", n, n, &c_n1, &c_n1, (ftnlen)6, (
301 ftnlen)1);
302 wrkbl = max(i__1,i__2);
303 /* Computing MAX */
304 i__1 = wrkbl, i__2 = bdspac + *n;
305 maxwrk = max(i__1,i__2);
306 minwrk = bdspac + *n;
307 } else if (wntqo) {
308
309 /* Path 2 (M much larger than N, JOBZ='O') */
310
311 wrkbl = *n + *n * ilaenv_(&c__1, "DGEQRF", " ", m, n, &
312 c_n1, &c_n1, (ftnlen)6, (ftnlen)1);
313 /* Computing MAX */
314 i__1 = wrkbl, i__2 = *n + *n * ilaenv_(&c__1, "DORGQR",
315 " ", m, n, n, &c_n1, (ftnlen)6, (ftnlen)1);
316 wrkbl = max(i__1,i__2);
317 /* Computing MAX */
318 i__1 = wrkbl, i__2 = *n * 3 + (*n << 1) * ilaenv_(&c__1,
319 "DGEBRD", " ", n, n, &c_n1, &c_n1, (ftnlen)6, (
320 ftnlen)1);
321 wrkbl = max(i__1,i__2);
322 /* Computing MAX */
323 i__1 = wrkbl, i__2 = *n * 3 + *n * ilaenv_(&c__1, "DORMBR"
324 , "QLN", n, n, n, &c_n1, (ftnlen)6, (ftnlen)3);
325 wrkbl = max(i__1,i__2);
326 /* Computing MAX */
327 i__1 = wrkbl, i__2 = *n * 3 + *n * ilaenv_(&c__1, "DORMBR"
328 , "PRT", n, n, n, &c_n1, (ftnlen)6, (ftnlen)3);
329 wrkbl = max(i__1,i__2);
330 /* Computing MAX */
331 i__1 = wrkbl, i__2 = bdspac + *n * 3;
332 wrkbl = max(i__1,i__2);
333 maxwrk = wrkbl + (*n << 1) * *n;
334 minwrk = bdspac + (*n << 1) * *n + *n * 3;
335 } else if (wntqs) {
336
337 /* Path 3 (M much larger than N, JOBZ='S') */
338
339 wrkbl = *n + *n * ilaenv_(&c__1, "DGEQRF", " ", m, n, &
340 c_n1, &c_n1, (ftnlen)6, (ftnlen)1);
341 /* Computing MAX */
342 i__1 = wrkbl, i__2 = *n + *n * ilaenv_(&c__1, "DORGQR",
343 " ", m, n, n, &c_n1, (ftnlen)6, (ftnlen)1);
344 wrkbl = max(i__1,i__2);
345 /* Computing MAX */
346 i__1 = wrkbl, i__2 = *n * 3 + (*n << 1) * ilaenv_(&c__1,
347 "DGEBRD", " ", n, n, &c_n1, &c_n1, (ftnlen)6, (
348 ftnlen)1);
349 wrkbl = max(i__1,i__2);
350 /* Computing MAX */
351 i__1 = wrkbl, i__2 = *n * 3 + *n * ilaenv_(&c__1, "DORMBR"
352 , "QLN", n, n, n, &c_n1, (ftnlen)6, (ftnlen)3);
353 wrkbl = max(i__1,i__2);
354 /* Computing MAX */
355 i__1 = wrkbl, i__2 = *n * 3 + *n * ilaenv_(&c__1, "DORMBR"
356 , "PRT", n, n, n, &c_n1, (ftnlen)6, (ftnlen)3);
357 wrkbl = max(i__1,i__2);
358 /* Computing MAX */
359 i__1 = wrkbl, i__2 = bdspac + *n * 3;
360 wrkbl = max(i__1,i__2);
361 maxwrk = wrkbl + *n * *n;
362 minwrk = bdspac + *n * *n + *n * 3;
363 } else if (wntqa) {
364
365 /* Path 4 (M much larger than N, JOBZ='A') */
366
367 wrkbl = *n + *n * ilaenv_(&c__1, "DGEQRF", " ", m, n, &
368 c_n1, &c_n1, (ftnlen)6, (ftnlen)1);
369 /* Computing MAX */
370 i__1 = wrkbl, i__2 = *n + *m * ilaenv_(&c__1, "DORGQR",
371 " ", m, m, n, &c_n1, (ftnlen)6, (ftnlen)1);
372 wrkbl = max(i__1,i__2);
373 /* Computing MAX */
374 i__1 = wrkbl, i__2 = *n * 3 + (*n << 1) * ilaenv_(&c__1,
375 "DGEBRD", " ", n, n, &c_n1, &c_n1, (ftnlen)6, (
376 ftnlen)1);
377 wrkbl = max(i__1,i__2);
378 /* Computing MAX */
379 i__1 = wrkbl, i__2 = *n * 3 + *n * ilaenv_(&c__1, "DORMBR"
380 , "QLN", n, n, n, &c_n1, (ftnlen)6, (ftnlen)3);
381 wrkbl = max(i__1,i__2);
382 /* Computing MAX */
383 i__1 = wrkbl, i__2 = *n * 3 + *n * ilaenv_(&c__1, "DORMBR"
384 , "PRT", n, n, n, &c_n1, (ftnlen)6, (ftnlen)3);
385 wrkbl = max(i__1,i__2);
386 /* Computing MAX */
387 i__1 = wrkbl, i__2 = bdspac + *n * 3;
388 wrkbl = max(i__1,i__2);
389 maxwrk = wrkbl + *n * *n;
390 minwrk = bdspac + *n * *n + *n * 3;
391 }
392 } else {
393
394 /* Path 5 (M at least N, but not much larger) */
395
396 wrkbl = *n * 3 + (*m + *n) * ilaenv_(&c__1, "DGEBRD", " ", m,
397 n, &c_n1, &c_n1, (ftnlen)6, (ftnlen)1);
398 if (wntqn) {
399 /* Computing MAX */
400 i__1 = wrkbl, i__2 = bdspac + *n * 3;
401 maxwrk = max(i__1,i__2);
402 minwrk = *n * 3 + max(*m,bdspac);
403 } else if (wntqo) {
404 /* Computing MAX */
405 i__1 = wrkbl, i__2 = *n * 3 + *n * ilaenv_(&c__1, "DORMBR"
406 , "QLN", m, n, n, &c_n1, (ftnlen)6, (ftnlen)3);
407 wrkbl = max(i__1,i__2);
408 /* Computing MAX */
409 i__1 = wrkbl, i__2 = *n * 3 + *n * ilaenv_(&c__1, "DORMBR"
410 , "PRT", n, n, n, &c_n1, (ftnlen)6, (ftnlen)3);
411 wrkbl = max(i__1,i__2);
412 /* Computing MAX */
413 i__1 = wrkbl, i__2 = bdspac + *n * 3;
414 wrkbl = max(i__1,i__2);
415 maxwrk = wrkbl + *m * *n;
416 /* Computing MAX */
417 i__1 = *m, i__2 = *n * *n + bdspac;
418 minwrk = *n * 3 + max(i__1,i__2);
419 } else if (wntqs) {
420 /* Computing MAX */
421 i__1 = wrkbl, i__2 = *n * 3 + *n * ilaenv_(&c__1, "DORMBR"
422 , "QLN", m, n, n, &c_n1, (ftnlen)6, (ftnlen)3);
423 wrkbl = max(i__1,i__2);
424 /* Computing MAX */
425 i__1 = wrkbl, i__2 = *n * 3 + *n * ilaenv_(&c__1, "DORMBR"
426 , "PRT", n, n, n, &c_n1, (ftnlen)6, (ftnlen)3);
427 wrkbl = max(i__1,i__2);
428 /* Computing MAX */
429 i__1 = wrkbl, i__2 = bdspac + *n * 3;
430 maxwrk = max(i__1,i__2);
431 minwrk = *n * 3 + max(*m,bdspac);
432 } else if (wntqa) {
433 /* Computing MAX */
434 i__1 = wrkbl, i__2 = *n * 3 + *m * ilaenv_(&c__1, "DORMBR"
435 , "QLN", m, m, n, &c_n1, (ftnlen)6, (ftnlen)3);
436 wrkbl = max(i__1,i__2);
437 /* Computing MAX */
438 i__1 = wrkbl, i__2 = *n * 3 + *n * ilaenv_(&c__1, "DORMBR"
439 , "PRT", n, n, n, &c_n1, (ftnlen)6, (ftnlen)3);
440 wrkbl = max(i__1,i__2);
441 /* Computing MAX */
442 i__1 = maxwrk, i__2 = bdspac + *n * 3;
443 maxwrk = max(i__1,i__2);
444 minwrk = *n * 3 + max(*m,bdspac);
445 }
446 }
447 } else {
448
449 /* Compute space needed for DBDSDC */
450
451 if (wntqn) {
452 bdspac = *m * 7;
453 } else {
454 bdspac = *m * 3 * *m + (*m << 2);
455 }
456 if (*n >= mnthr) {
457 if (wntqn) {
458
459 /* Path 1t (N much larger than M, JOBZ='N') */
460
461 wrkbl = *m + *m * ilaenv_(&c__1, "DGELQF", " ", m, n, &
462 c_n1, &c_n1, (ftnlen)6, (ftnlen)1);
463 /* Computing MAX */
464 i__1 = wrkbl, i__2 = *m * 3 + (*m << 1) * ilaenv_(&c__1,
465 "DGEBRD", " ", m, m, &c_n1, &c_n1, (ftnlen)6, (
466 ftnlen)1);
467 wrkbl = max(i__1,i__2);
468 /* Computing MAX */
469 i__1 = wrkbl, i__2 = bdspac + *m;
470 maxwrk = max(i__1,i__2);
471 minwrk = bdspac + *m;
472 } else if (wntqo) {
473
474 /* Path 2t (N much larger than M, JOBZ='O') */
475
476 wrkbl = *m + *m * ilaenv_(&c__1, "DGELQF", " ", m, n, &
477 c_n1, &c_n1, (ftnlen)6, (ftnlen)1);
478 /* Computing MAX */
479 i__1 = wrkbl, i__2 = *m + *m * ilaenv_(&c__1, "DORGLQ",
480 " ", m, n, m, &c_n1, (ftnlen)6, (ftnlen)1);
481 wrkbl = max(i__1,i__2);
482 /* Computing MAX */
483 i__1 = wrkbl, i__2 = *m * 3 + (*m << 1) * ilaenv_(&c__1,
484 "DGEBRD", " ", m, m, &c_n1, &c_n1, (ftnlen)6, (
485 ftnlen)1);
486 wrkbl = max(i__1,i__2);
487 /* Computing MAX */
488 i__1 = wrkbl, i__2 = *m * 3 + *m * ilaenv_(&c__1, "DORMBR"
489 , "QLN", m, m, m, &c_n1, (ftnlen)6, (ftnlen)3);
490 wrkbl = max(i__1,i__2);
491 /* Computing MAX */
492 i__1 = wrkbl, i__2 = *m * 3 + *m * ilaenv_(&c__1, "DORMBR"
493 , "PRT", m, m, m, &c_n1, (ftnlen)6, (ftnlen)3);
494 wrkbl = max(i__1,i__2);
495 /* Computing MAX */
496 i__1 = wrkbl, i__2 = bdspac + *m * 3;
497 wrkbl = max(i__1,i__2);
498 maxwrk = wrkbl + (*m << 1) * *m;
499 minwrk = bdspac + (*m << 1) * *m + *m * 3;
500 } else if (wntqs) {
501
502 /* Path 3t (N much larger than M, JOBZ='S') */
503
504 wrkbl = *m + *m * ilaenv_(&c__1, "DGELQF", " ", m, n, &
505 c_n1, &c_n1, (ftnlen)6, (ftnlen)1);
506 /* Computing MAX */
507 i__1 = wrkbl, i__2 = *m + *m * ilaenv_(&c__1, "DORGLQ",
508 " ", m, n, m, &c_n1, (ftnlen)6, (ftnlen)1);
509 wrkbl = max(i__1,i__2);
510 /* Computing MAX */
511 i__1 = wrkbl, i__2 = *m * 3 + (*m << 1) * ilaenv_(&c__1,
512 "DGEBRD", " ", m, m, &c_n1, &c_n1, (ftnlen)6, (
513 ftnlen)1);
514 wrkbl = max(i__1,i__2);
515 /* Computing MAX */
516 i__1 = wrkbl, i__2 = *m * 3 + *m * ilaenv_(&c__1, "DORMBR"
517 , "QLN", m, m, m, &c_n1, (ftnlen)6, (ftnlen)3);
518 wrkbl = max(i__1,i__2);
519 /* Computing MAX */
520 i__1 = wrkbl, i__2 = *m * 3 + *m * ilaenv_(&c__1, "DORMBR"
521 , "PRT", m, m, m, &c_n1, (ftnlen)6, (ftnlen)3);
522 wrkbl = max(i__1,i__2);
523 /* Computing MAX */
524 i__1 = wrkbl, i__2 = bdspac + *m * 3;
525 wrkbl = max(i__1,i__2);
526 maxwrk = wrkbl + *m * *m;
527 minwrk = bdspac + *m * *m + *m * 3;
528 } else if (wntqa) {
529
530 /* Path 4t (N much larger than M, JOBZ='A') */
531
532 wrkbl = *m + *m * ilaenv_(&c__1, "DGELQF", " ", m, n, &
533 c_n1, &c_n1, (ftnlen)6, (ftnlen)1);
534 /* Computing MAX */
535 i__1 = wrkbl, i__2 = *m + *n * ilaenv_(&c__1, "DORGLQ",
536 " ", n, n, m, &c_n1, (ftnlen)6, (ftnlen)1);
537 wrkbl = max(i__1,i__2);
538 /* Computing MAX */
539 i__1 = wrkbl, i__2 = *m * 3 + (*m << 1) * ilaenv_(&c__1,
540 "DGEBRD", " ", m, m, &c_n1, &c_n1, (ftnlen)6, (
541 ftnlen)1);
542 wrkbl = max(i__1,i__2);
543 /* Computing MAX */
544 i__1 = wrkbl, i__2 = *m * 3 + *m * ilaenv_(&c__1, "DORMBR"
545 , "QLN", m, m, m, &c_n1, (ftnlen)6, (ftnlen)3);
546 wrkbl = max(i__1,i__2);
547 /* Computing MAX */
548 i__1 = wrkbl, i__2 = *m * 3 + *m * ilaenv_(&c__1, "DORMBR"
549 , "PRT", m, m, m, &c_n1, (ftnlen)6, (ftnlen)3);
550 wrkbl = max(i__1,i__2);
551 /* Computing MAX */
552 i__1 = wrkbl, i__2 = bdspac + *m * 3;
553 wrkbl = max(i__1,i__2);
554 maxwrk = wrkbl + *m * *m;
555 minwrk = bdspac + *m * *m + *m * 3;
556 }
557 } else {
558
559 /* Path 5t (N greater than M, but not much larger) */
560
561 wrkbl = *m * 3 + (*m + *n) * ilaenv_(&c__1, "DGEBRD", " ", m,
562 n, &c_n1, &c_n1, (ftnlen)6, (ftnlen)1);
563 if (wntqn) {
564 /* Computing MAX */
565 i__1 = wrkbl, i__2 = bdspac + *m * 3;
566 maxwrk = max(i__1,i__2);
567 minwrk = *m * 3 + max(*n,bdspac);
568 } else if (wntqo) {
569 /* Computing MAX */
570 i__1 = wrkbl, i__2 = *m * 3 + *m * ilaenv_(&c__1, "DORMBR"
571 , "QLN", m, m, n, &c_n1, (ftnlen)6, (ftnlen)3);
572 wrkbl = max(i__1,i__2);
573 /* Computing MAX */
574 i__1 = wrkbl, i__2 = *m * 3 + *m * ilaenv_(&c__1, "DORMBR"
575 , "PRT", m, n, m, &c_n1, (ftnlen)6, (ftnlen)3);
576 wrkbl = max(i__1,i__2);
577 /* Computing MAX */
578 i__1 = wrkbl, i__2 = bdspac + *m * 3;
579 wrkbl = max(i__1,i__2);
580 maxwrk = wrkbl + *m * *n;
581 /* Computing MAX */
582 i__1 = *n, i__2 = *m * *m + bdspac;
583 minwrk = *m * 3 + max(i__1,i__2);
584 } else if (wntqs) {
585 /* Computing MAX */
586 i__1 = wrkbl, i__2 = *m * 3 + *m * ilaenv_(&c__1, "DORMBR"
587 , "QLN", m, m, n, &c_n1, (ftnlen)6, (ftnlen)3);
588 wrkbl = max(i__1,i__2);
589 /* Computing MAX */
590 i__1 = wrkbl, i__2 = *m * 3 + *m * ilaenv_(&c__1, "DORMBR"
591 , "PRT", m, n, m, &c_n1, (ftnlen)6, (ftnlen)3);
592 wrkbl = max(i__1,i__2);
593 /* Computing MAX */
594 i__1 = wrkbl, i__2 = bdspac + *m * 3;
595 maxwrk = max(i__1,i__2);
596 minwrk = *m * 3 + max(*n,bdspac);
597 } else if (wntqa) {
598 /* Computing MAX */
599 i__1 = wrkbl, i__2 = *m * 3 + *m * ilaenv_(&c__1, "DORMBR"
600 , "QLN", m, m, n, &c_n1, (ftnlen)6, (ftnlen)3);
601 wrkbl = max(i__1,i__2);
602 /* Computing MAX */
603 i__1 = wrkbl, i__2 = *m * 3 + *m * ilaenv_(&c__1, "DORMBR"
604 , "PRT", n, n, m, &c_n1, (ftnlen)6, (ftnlen)3);
605 wrkbl = max(i__1,i__2);
606 /* Computing MAX */
607 i__1 = wrkbl, i__2 = bdspac + *m * 3;
608 maxwrk = max(i__1,i__2);
609 minwrk = *m * 3 + max(*n,bdspac);
610 }
611 }
612 }
613 work[1] = (doublereal) maxwrk;
614 }
615
616 if (*lwork < minwrk && ! lquery) {
617 *info = -12;
618 }
619 if (*info != 0) {
620 i__1 = -(*info);
621 xerbla_("DGESDD", &i__1, (ftnlen)6);
622 return 0;
623 } else if (lquery) {
624 return 0;
625 }
626
627 /* Quick return if possible */
628
629 if (*m == 0 || *n == 0) {
630 if (*lwork >= 1) {
631 work[1] = 1.;
632 }
633 return 0;
634 }
635
636 /* Get machine constants */
637
638 eps = dlamch_("P", (ftnlen)1);
639 smlnum = sqrt(dlamch_("S", (ftnlen)1)) / eps;
640 bignum = 1. / smlnum;
641
642 /* Scale A if max element outside range [SMLNUM,BIGNUM] */
643
644 anrm = dlange_("M", m, n, &a[a_offset], lda, dum, (ftnlen)1);
645 iscl = 0;
646 if (anrm > 0. && anrm < smlnum) {
647 iscl = 1;
648 dlascl_("G", &c__0, &c__0, &anrm, &smlnum, m, n, &a[a_offset], lda, &
649 ierr, (ftnlen)1);
650 } else if (anrm > bignum) {
651 iscl = 1;
652 dlascl_("G", &c__0, &c__0, &anrm, &bignum, m, n, &a[a_offset], lda, &
653 ierr, (ftnlen)1);
654 }
655
656 if (*m >= *n) {
657
658 /* A has at least as many rows as columns. If A has sufficiently */
659 /* more rows than columns, first reduce using the QR */
660 /* decomposition (if sufficient workspace available) */
661
662 if (*m >= mnthr) {
663
664 if (wntqn) {
665
666 /* Path 1 (M much larger than N, JOBZ='N') */
667 /* No singular vectors to be computed */
668
669 itau = 1;
670 nwork = itau + *n;
671
672 /* Compute A=Q*R */
673 /* (Workspace: need 2*N, prefer N+N*NB) */
674
675 i__1 = *lwork - nwork + 1;
676 dgeqrf_(m, n, &a[a_offset], lda, &work[itau], &work[nwork], &
677 i__1, &ierr);
678
679 /* Zero out below R */
680
681 i__1 = *n - 1;
682 i__2 = *n - 1;
683 dlaset_("L", &i__1, &i__2, &c_b227, &c_b227, &a[a_dim1 + 2],
684 lda, (ftnlen)1);
685 ie = 1;
686 itauq = ie + *n;
687 itaup = itauq + *n;
688 nwork = itaup + *n;
689
690 /* Bidiagonalize R in A */
691 /* (Workspace: need 4*N, prefer 3*N+2*N*NB) */
692
693 i__1 = *lwork - nwork + 1;
694 dgebrd_(n, n, &a[a_offset], lda, &s[1], &work[ie], &work[
695 itauq], &work[itaup], &work[nwork], &i__1, &ierr);
696 nwork = ie + *n;
697
698 /* Perform bidiagonal SVD, computing singular values only */
699 /* (Workspace: need N+BDSPAC) */
700
701 dbdsdc_("U", "N", n, &s[1], &work[ie], dum, &c__1, dum, &c__1,
702 dum, idum, &work[nwork], &iwork[1], info, (ftnlen)1,
703 (ftnlen)1);
704
705 } else if (wntqo) {
706
707 /* Path 2 (M much larger than N, JOBZ = 'O') */
708 /* N left singular vectors to be overwritten on A and */
709 /* N right singular vectors to be computed in VT */
710
711 ir = 1;
712
713 /* WORK(IR) is LDWRKR by N */
714
715 if (*lwork >= *lda * *n + *n * *n + *n * 3 + bdspac) {
716 ldwrkr = *lda;
717 } else {
718 ldwrkr = (*lwork - *n * *n - *n * 3 - bdspac) / *n;
719 }
720 itau = ir + ldwrkr * *n;
721 nwork = itau + *n;
722
723 /* Compute A=Q*R */
724 /* (Workspace: need N*N+2*N, prefer N*N+N+N*NB) */
725
726 i__1 = *lwork - nwork + 1;
727 dgeqrf_(m, n, &a[a_offset], lda, &work[itau], &work[nwork], &
728 i__1, &ierr);
729
730 /* Copy R to WORK(IR), zeroing out below it */
731
732 dlacpy_("U", n, n, &a[a_offset], lda, &work[ir], &ldwrkr, (
733 ftnlen)1);
734 i__1 = *n - 1;
735 i__2 = *n - 1;
736 dlaset_("L", &i__1, &i__2, &c_b227, &c_b227, &work[ir + 1], &
737 ldwrkr, (ftnlen)1);
738
739 /* Generate Q in A */
740 /* (Workspace: need N*N+2*N, prefer N*N+N+N*NB) */
741
742 i__1 = *lwork - nwork + 1;
743 dorgqr_(m, n, n, &a[a_offset], lda, &work[itau], &work[nwork],
744 &i__1, &ierr);
745 ie = itau;
746 itauq = ie + *n;
747 itaup = itauq + *n;
748 nwork = itaup + *n;
749
750 /* Bidiagonalize R in VT, copying result to WORK(IR) */
751 /* (Workspace: need N*N+4*N, prefer N*N+3*N+2*N*NB) */
752
753 i__1 = *lwork - nwork + 1;
754 dgebrd_(n, n, &work[ir], &ldwrkr, &s[1], &work[ie], &work[
755 itauq], &work[itaup], &work[nwork], &i__1, &ierr);
756
757 /* WORK(IU) is N by N */
758
759 iu = nwork;
760 nwork = iu + *n * *n;
761
762 /* Perform bidiagonal SVD, computing left singular vectors */
763 /* of bidiagonal matrix in WORK(IU) and computing right */
764 /* singular vectors of bidiagonal matrix in VT */
765 /* (Workspace: need N+N*N+BDSPAC) */
766
767 dbdsdc_("U", "I", n, &s[1], &work[ie], &work[iu], n, &vt[
768 vt_offset], ldvt, dum, idum, &work[nwork], &iwork[1],
769 info, (ftnlen)1, (ftnlen)1);
770
771 /* Overwrite WORK(IU) by left singular vectors of R */
772 /* and VT by right singular vectors of R */
773 /* (Workspace: need 2*N*N+3*N, prefer 2*N*N+2*N+N*NB) */
774
775 i__1 = *lwork - nwork + 1;
776 dormbr_("Q", "L", "N", n, n, n, &work[ir], &ldwrkr, &work[
777 itauq], &work[iu], n, &work[nwork], &i__1, &ierr, (
778 ftnlen)1, (ftnlen)1, (ftnlen)1);
779 i__1 = *lwork - nwork + 1;
780 dormbr_("P", "R", "T", n, n, n, &work[ir], &ldwrkr, &work[
781 itaup], &vt[vt_offset], ldvt, &work[nwork], &i__1, &
782 ierr, (ftnlen)1, (ftnlen)1, (ftnlen)1);
783
784 /* Multiply Q in A by left singular vectors of R in */
785 /* WORK(IU), storing result in WORK(IR) and copying to A */
786 /* (Workspace: need 2*N*N, prefer N*N+M*N) */
787
788 i__1 = *m;
789 i__2 = ldwrkr;
790 for (i__ = 1; i__2 < 0 ? i__ >= i__1 : i__ <= i__1; i__ +=
791 i__2) {
792 /* Computing MIN */
793 i__3 = *m - i__ + 1;
794 chunk = min(i__3,ldwrkr);
795 dgemm_("N", "N", &chunk, n, n, &c_b248, &a[i__ + a_dim1],
796 lda, &work[iu], n, &c_b227, &work[ir], &ldwrkr, (
797 ftnlen)1, (ftnlen)1);
798 dlacpy_("F", &chunk, n, &work[ir], &ldwrkr, &a[i__ +
799 a_dim1], lda, (ftnlen)1);
800 /* L10: */
801 }
802
803 } else if (wntqs) {
804
805 /* Path 3 (M much larger than N, JOBZ='S') */
806 /* N left singular vectors to be computed in U and */
807 /* N right singular vectors to be computed in VT */
808
809 ir = 1;
810
811 /* WORK(IR) is N by N */
812
813 ldwrkr = *n;
814 itau = ir + ldwrkr * *n;
815 nwork = itau + *n;
816
817 /* Compute A=Q*R */
818 /* (Workspace: need N*N+2*N, prefer N*N+N+N*NB) */
819
820 i__2 = *lwork - nwork + 1;
821 dgeqrf_(m, n, &a[a_offset], lda, &work[itau], &work[nwork], &
822 i__2, &ierr);
823
824 /* Copy R to WORK(IR), zeroing out below it */
825
826 dlacpy_("U", n, n, &a[a_offset], lda, &work[ir], &ldwrkr, (
827 ftnlen)1);
828 i__2 = *n - 1;
829 i__1 = *n - 1;
830 dlaset_("L", &i__2, &i__1, &c_b227, &c_b227, &work[ir + 1], &
831 ldwrkr, (ftnlen)1);
832
833 /* Generate Q in A */
834 /* (Workspace: need N*N+2*N, prefer N*N+N+N*NB) */
835
836 i__2 = *lwork - nwork + 1;
837 dorgqr_(m, n, n, &a[a_offset], lda, &work[itau], &work[nwork],
838 &i__2, &ierr);
839 ie = itau;
840 itauq = ie + *n;
841 itaup = itauq + *n;
842 nwork = itaup + *n;
843
844 /* Bidiagonalize R in WORK(IR) */
845 /* (Workspace: need N*N+4*N, prefer N*N+3*N+2*N*NB) */
846
847 i__2 = *lwork - nwork + 1;
848 dgebrd_(n, n, &work[ir], &ldwrkr, &s[1], &work[ie], &work[
849 itauq], &work[itaup], &work[nwork], &i__2, &ierr);
850
851 /* Perform bidiagonal SVD, computing left singular vectors */
852 /* of bidiagoal matrix in U and computing right singular */
853 /* vectors of bidiagonal matrix in VT */
854 /* (Workspace: need N+BDSPAC) */
855
856 dbdsdc_("U", "I", n, &s[1], &work[ie], &u[u_offset], ldu, &vt[
857 vt_offset], ldvt, dum, idum, &work[nwork], &iwork[1],
858 info, (ftnlen)1, (ftnlen)1);
859
860 /* Overwrite U by left singular vectors of R and VT */
861 /* by right singular vectors of R */
862 /* (Workspace: need N*N+3*N, prefer N*N+2*N+N*NB) */
863
864 i__2 = *lwork - nwork + 1;
865 dormbr_("Q", "L", "N", n, n, n, &work[ir], &ldwrkr, &work[
866 itauq], &u[u_offset], ldu, &work[nwork], &i__2, &ierr,
867 (ftnlen)1, (ftnlen)1, (ftnlen)1);
868
869 i__2 = *lwork - nwork + 1;
870 dormbr_("P", "R", "T", n, n, n, &work[ir], &ldwrkr, &work[
871 itaup], &vt[vt_offset], ldvt, &work[nwork], &i__2, &
872 ierr, (ftnlen)1, (ftnlen)1, (ftnlen)1);
873
874 /* Multiply Q in A by left singular vectors of R in */
875 /* WORK(IR), storing result in U */
876 /* (Workspace: need N*N) */
877
878 dlacpy_("F", n, n, &u[u_offset], ldu, &work[ir], &ldwrkr, (
879 ftnlen)1);
880 dgemm_("N", "N", m, n, n, &c_b248, &a[a_offset], lda, &work[
881 ir], &ldwrkr, &c_b227, &u[u_offset], ldu, (ftnlen)1, (
882 ftnlen)1);
883
884 } else if (wntqa) {
885
886 /* Path 4 (M much larger than N, JOBZ='A') */
887 /* M left singular vectors to be computed in U and */
888 /* N right singular vectors to be computed in VT */
889
890 iu = 1;
891
892 /* WORK(IU) is N by N */
893
894 ldwrku = *n;
895 itau = iu + ldwrku * *n;
896 nwork = itau + *n;
897
898 /* Compute A=Q*R, copying result to U */
899 /* (Workspace: need N*N+2*N, prefer N*N+N+N*NB) */
900
901 i__2 = *lwork - nwork + 1;
902 dgeqrf_(m, n, &a[a_offset], lda, &work[itau], &work[nwork], &
903 i__2, &ierr);
904 dlacpy_("L", m, n, &a[a_offset], lda, &u[u_offset], ldu, (
905 ftnlen)1);
906
907 /* Generate Q in U */
908 /* (Workspace: need N*N+2*N, prefer N*N+N+N*NB) */
909 i__2 = *lwork - nwork + 1;
910 dorgqr_(m, m, n, &u[u_offset], ldu, &work[itau], &work[nwork],
911 &i__2, &ierr);
912
913 /* Produce R in A, zeroing out other entries */
914
915 i__2 = *n - 1;
916 i__1 = *n - 1;
917 dlaset_("L", &i__2, &i__1, &c_b227, &c_b227, &a[a_dim1 + 2],
918 lda, (ftnlen)1);
919 ie = itau;
920 itauq = ie + *n;
921 itaup = itauq + *n;
922 nwork = itaup + *n;
923
924 /* Bidiagonalize R in A */
925 /* (Workspace: need N*N+4*N, prefer N*N+3*N+2*N*NB) */
926
927 i__2 = *lwork - nwork + 1;
928 dgebrd_(n, n, &a[a_offset], lda, &s[1], &work[ie], &work[
929 itauq], &work[itaup], &work[nwork], &i__2, &ierr);
930
931 /* Perform bidiagonal SVD, computing left singular vectors */
932 /* of bidiagonal matrix in WORK(IU) and computing right */
933 /* singular vectors of bidiagonal matrix in VT */
934 /* (Workspace: need N+N*N+BDSPAC) */
935
936 dbdsdc_("U", "I", n, &s[1], &work[ie], &work[iu], n, &vt[
937 vt_offset], ldvt, dum, idum, &work[nwork], &iwork[1],
938 info, (ftnlen)1, (ftnlen)1);
939
940 /* Overwrite WORK(IU) by left singular vectors of R and VT */
941 /* by right singular vectors of R */
942 /* (Workspace: need N*N+3*N, prefer N*N+2*N+N*NB) */
943
944 i__2 = *lwork - nwork + 1;
945 dormbr_("Q", "L", "N", n, n, n, &a[a_offset], lda, &work[
946 itauq], &work[iu], &ldwrku, &work[nwork], &i__2, &
947 ierr, (ftnlen)1, (ftnlen)1, (ftnlen)1);
948 i__2 = *lwork - nwork + 1;
949 dormbr_("P", "R", "T", n, n, n, &a[a_offset], lda, &work[
950 itaup], &vt[vt_offset], ldvt, &work[nwork], &i__2, &
951 ierr, (ftnlen)1, (ftnlen)1, (ftnlen)1);
952
953 /* Multiply Q in U by left singular vectors of R in */
954 /* WORK(IU), storing result in A */
955 /* (Workspace: need N*N) */
956
957 dgemm_("N", "N", m, n, n, &c_b248, &u[u_offset], ldu, &work[
958 iu], &ldwrku, &c_b227, &a[a_offset], lda, (ftnlen)1, (
959 ftnlen)1);
960
961 /* Copy left singular vectors of A from A to U */
962
963 dlacpy_("F", m, n, &a[a_offset], lda, &u[u_offset], ldu, (
964 ftnlen)1);
965
966 }
967
968 } else {
969
970 /* M .LT. MNTHR */
971
972 /* Path 5 (M at least N, but not much larger) */
973 /* Reduce to bidiagonal form without QR decomposition */
974
975 ie = 1;
976 itauq = ie + *n;
977 itaup = itauq + *n;
978 nwork = itaup + *n;
979
980 /* Bidiagonalize A */
981 /* (Workspace: need 3*N+M, prefer 3*N+(M+N)*NB) */
982
983 i__2 = *lwork - nwork + 1;
984 dgebrd_(m, n, &a[a_offset], lda, &s[1], &work[ie], &work[itauq], &
985 work[itaup], &work[nwork], &i__2, &ierr);
986 if (wntqn) {
987
988 /* Perform bidiagonal SVD, only computing singular values */
989 /* (Workspace: need N+BDSPAC) */
990
991 dbdsdc_("U", "N", n, &s[1], &work[ie], dum, &c__1, dum, &c__1,
992 dum, idum, &work[nwork], &iwork[1], info, (ftnlen)1,
993 (ftnlen)1);
994 } else if (wntqo) {
995 iu = nwork;
996 if (*lwork >= *m * *n + *n * 3 + bdspac) {
997
998 /* WORK( IU ) is M by N */
999
1000 ldwrku = *m;
1001 nwork = iu + ldwrku * *n;
1002 dlaset_("F", m, n, &c_b227, &c_b227, &work[iu], &ldwrku, (
1003 ftnlen)1);
1004 } else {
1005
1006 /* WORK( IU ) is N by N */
1007
1008 ldwrku = *n;
1009 nwork = iu + ldwrku * *n;
1010
1011 /* WORK(IR) is LDWRKR by N */
1012
1013 ir = nwork;
1014 ldwrkr = (*lwork - *n * *n - *n * 3) / *n;
1015 }
1016 nwork = iu + ldwrku * *n;
1017
1018 /* Perform bidiagonal SVD, computing left singular vectors */
1019 /* of bidiagonal matrix in WORK(IU) and computing right */
1020 /* singular vectors of bidiagonal matrix in VT */
1021 /* (Workspace: need N+N*N+BDSPAC) */
1022
1023 dbdsdc_("U", "I", n, &s[1], &work[ie], &work[iu], &ldwrku, &
1024 vt[vt_offset], ldvt, dum, idum, &work[nwork], &iwork[
1025 1], info, (ftnlen)1, (ftnlen)1);
1026
1027 /* Overwrite VT by right singular vectors of A */
1028 /* (Workspace: need N*N+2*N, prefer N*N+N+N*NB) */
1029
1030 i__2 = *lwork - nwork + 1;
1031 dormbr_("P", "R", "T", n, n, n, &a[a_offset], lda, &work[
1032 itaup], &vt[vt_offset], ldvt, &work[nwork], &i__2, &
1033 ierr, (ftnlen)1, (ftnlen)1, (ftnlen)1);
1034
1035 if (*lwork >= *m * *n + *n * 3 + bdspac) {
1036
1037 /* Overwrite WORK(IU) by left singular vectors of A */
1038 /* (Workspace: need N*N+2*N, prefer N*N+N+N*NB) */
1039
1040 i__2 = *lwork - nwork + 1;
1041 dormbr_("Q", "L", "N", m, n, n, &a[a_offset], lda, &work[
1042 itauq], &work[iu], &ldwrku, &work[nwork], &i__2, &
1043 ierr, (ftnlen)1, (ftnlen)1, (ftnlen)1);
1044
1045 /* Copy left singular vectors of A from WORK(IU) to A */
1046
1047 dlacpy_("F", m, n, &work[iu], &ldwrku, &a[a_offset], lda,
1048 (ftnlen)1);
1049 } else {
1050
1051 /* Generate Q in A */
1052 /* (Workspace: need N*N+2*N, prefer N*N+N+N*NB) */
1053
1054 i__2 = *lwork - nwork + 1;
1055 dorgbr_("Q", m, n, n, &a[a_offset], lda, &work[itauq], &
1056 work[nwork], &i__2, &ierr, (ftnlen)1);
1057
1058 /* Multiply Q in A by left singular vectors of */
1059 /* bidiagonal matrix in WORK(IU), storing result in */
1060 /* WORK(IR) and copying to A */
1061 /* (Workspace: need 2*N*N, prefer N*N+M*N) */
1062
1063 i__2 = *m;
1064 i__1 = ldwrkr;
1065 for (i__ = 1; i__1 < 0 ? i__ >= i__2 : i__ <= i__2; i__ +=
1066 i__1) {
1067 /* Computing MIN */
1068 i__3 = *m - i__ + 1;
1069 chunk = min(i__3,ldwrkr);
1070 dgemm_("N", "N", &chunk, n, n, &c_b248, &a[i__ +
1071 a_dim1], lda, &work[iu], &ldwrku, &c_b227, &
1072 work[ir], &ldwrkr, (ftnlen)1, (ftnlen)1);
1073 dlacpy_("F", &chunk, n, &work[ir], &ldwrkr, &a[i__ +
1074 a_dim1], lda, (ftnlen)1);
1075 /* L20: */
1076 }
1077 }
1078
1079 } else if (wntqs) {
1080
1081 /* Perform bidiagonal SVD, computing left singular vectors */
1082 /* of bidiagonal matrix in U and computing right singular */
1083 /* vectors of bidiagonal matrix in VT */
1084 /* (Workspace: need N+BDSPAC) */
1085
1086 dlaset_("F", m, n, &c_b227, &c_b227, &u[u_offset], ldu, (
1087 ftnlen)1);
1088 dbdsdc_("U", "I", n, &s[1], &work[ie], &u[u_offset], ldu, &vt[
1089 vt_offset], ldvt, dum, idum, &work[nwork], &iwork[1],
1090 info, (ftnlen)1, (ftnlen)1);
1091
1092 /* Overwrite U by left singular vectors of A and VT */
1093 /* by right singular vectors of A */
1094 /* (Workspace: need 3*N, prefer 2*N+N*NB) */
1095
1096 i__1 = *lwork - nwork + 1;
1097 dormbr_("Q", "L", "N", m, n, n, &a[a_offset], lda, &work[
1098 itauq], &u[u_offset], ldu, &work[nwork], &i__1, &ierr,
1099 (ftnlen)1, (ftnlen)1, (ftnlen)1);
1100 i__1 = *lwork - nwork + 1;
1101 dormbr_("P", "R", "T", n, n, n, &a[a_offset], lda, &work[
1102 itaup], &vt[vt_offset], ldvt, &work[nwork], &i__1, &
1103 ierr, (ftnlen)1, (ftnlen)1, (ftnlen)1);
1104 } else if (wntqa) {
1105
1106 /* Perform bidiagonal SVD, computing left singular vectors */
1107 /* of bidiagonal matrix in U and computing right singular */
1108 /* vectors of bidiagonal matrix in VT */
1109 /* (Workspace: need N+BDSPAC) */
1110
1111 dlaset_("F", m, m, &c_b227, &c_b227, &u[u_offset], ldu, (
1112 ftnlen)1);
1113 dbdsdc_("U", "I", n, &s[1], &work[ie], &u[u_offset], ldu, &vt[
1114 vt_offset], ldvt, dum, idum, &work[nwork], &iwork[1],
1115 info, (ftnlen)1, (ftnlen)1);
1116
1117 /* Set the right corner of U to identity matrix */
1118
1119 i__1 = *m - *n;
1120 i__2 = *m - *n;
1121 dlaset_("F", &i__1, &i__2, &c_b227, &c_b248, &u[*n + 1 + (*n
1122 + 1) * u_dim1], ldu, (ftnlen)1);
1123
1124 /* Overwrite U by left singular vectors of A and VT */
1125 /* by right singular vectors of A */
1126 /* (Workspace: need N*N+2*N+M, prefer N*N+2*N+M*NB) */
1127
1128 i__1 = *lwork - nwork + 1;
1129 dormbr_("Q", "L", "N", m, m, n, &a[a_offset], lda, &work[
1130 itauq], &u[u_offset], ldu, &work[nwork], &i__1, &ierr,
1131 (ftnlen)1, (ftnlen)1, (ftnlen)1);
1132 i__1 = *lwork - nwork + 1;
1133 dormbr_("P", "R", "T", n, n, m, &a[a_offset], lda, &work[
1134 itaup], &vt[vt_offset], ldvt, &work[nwork], &i__1, &
1135 ierr, (ftnlen)1, (ftnlen)1, (ftnlen)1);
1136 }
1137
1138 }
1139
1140 } else {
1141
1142 /* A has more columns than rows. If A has sufficiently more */
1143 /* columns than rows, first reduce using the LQ decomposition (if */
1144 /* sufficient workspace available) */
1145
1146 if (*n >= mnthr) {
1147
1148 if (wntqn) {
1149
1150 /* Path 1t (N much larger than M, JOBZ='N') */
1151 /* No singular vectors to be computed */
1152
1153 itau = 1;
1154 nwork = itau + *m;
1155
1156 /* Compute A=L*Q */
1157 /* (Workspace: need 2*M, prefer M+M*NB) */
1158
1159 i__1 = *lwork - nwork + 1;
1160 dgelqf_(m, n, &a[a_offset], lda, &work[itau], &work[nwork], &
1161 i__1, &ierr);
1162
1163 /* Zero out above L */
1164
1165 i__1 = *m - 1;
1166 i__2 = *m - 1;
1167 dlaset_("U", &i__1, &i__2, &c_b227, &c_b227, &a[(a_dim1 << 1)
1168 + 1], lda, (ftnlen)1);
1169 ie = 1;
1170 itauq = ie + *m;
1171 itaup = itauq + *m;
1172 nwork = itaup + *m;
1173
1174 /* Bidiagonalize L in A */
1175 /* (Workspace: need 4*M, prefer 3*M+2*M*NB) */
1176
1177 i__1 = *lwork - nwork + 1;
1178 dgebrd_(m, m, &a[a_offset], lda, &s[1], &work[ie], &work[
1179 itauq], &work[itaup], &work[nwork], &i__1, &ierr);
1180 nwork = ie + *m;
1181
1182 /* Perform bidiagonal SVD, computing singular values only */
1183 /* (Workspace: need M+BDSPAC) */
1184
1185 dbdsdc_("U", "N", m, &s[1], &work[ie], dum, &c__1, dum, &c__1,
1186 dum, idum, &work[nwork], &iwork[1], info, (ftnlen)1,
1187 (ftnlen)1);
1188
1189 } else if (wntqo) {
1190
1191 /* Path 2t (N much larger than M, JOBZ='O') */
1192 /* M right singular vectors to be overwritten on A and */
1193 /* M left singular vectors to be computed in U */
1194
1195 ivt = 1;
1196
1197 /* IVT is M by M */
1198
1199 il = ivt + *m * *m;
1200 if (*lwork >= *m * *n + *m * *m + *m * 3 + bdspac) {
1201
1202 /* WORK(IL) is M by N */
1203
1204 ldwrkl = *m;
1205 chunk = *n;
1206 } else {
1207 ldwrkl = *m;
1208 chunk = (*lwork - *m * *m) / *m;
1209 }
1210 itau = il + ldwrkl * *m;
1211 nwork = itau + *m;
1212
1213 /* Compute A=L*Q */
1214 /* (Workspace: need M*M+2*M, prefer M*M+M+M*NB) */
1215
1216 i__1 = *lwork - nwork + 1;
1217 dgelqf_(m, n, &a[a_offset], lda, &work[itau], &work[nwork], &
1218 i__1, &ierr);
1219
1220 /* Copy L to WORK(IL), zeroing about above it */
1221
1222 dlacpy_("L", m, m, &a[a_offset], lda, &work[il], &ldwrkl, (
1223 ftnlen)1);
1224 i__1 = *m - 1;
1225 i__2 = *m - 1;
1226 dlaset_("U", &i__1, &i__2, &c_b227, &c_b227, &work[il +
1227 ldwrkl], &ldwrkl, (ftnlen)1);
1228
1229 /* Generate Q in A */
1230 /* (Workspace: need M*M+2*M, prefer M*M+M+M*NB) */
1231
1232 i__1 = *lwork - nwork + 1;
1233 dorglq_(m, n, m, &a[a_offset], lda, &work[itau], &work[nwork],
1234 &i__1, &ierr);
1235 ie = itau;
1236 itauq = ie + *m;
1237 itaup = itauq + *m;
1238 nwork = itaup + *m;
1239
1240 /* Bidiagonalize L in WORK(IL) */
1241 /* (Workspace: need M*M+4*M, prefer M*M+3*M+2*M*NB) */
1242
1243 i__1 = *lwork - nwork + 1;
1244 dgebrd_(m, m, &work[il], &ldwrkl, &s[1], &work[ie], &work[
1245 itauq], &work[itaup], &work[nwork], &i__1, &ierr);
1246
1247 /* Perform bidiagonal SVD, computing left singular vectors */
1248 /* of bidiagonal matrix in U, and computing right singular */
1249 /* vectors of bidiagonal matrix in WORK(IVT) */
1250 /* (Workspace: need M+M*M+BDSPAC) */
1251
1252 dbdsdc_("U", "I", m, &s[1], &work[ie], &u[u_offset], ldu, &
1253 work[ivt], m, dum, idum, &work[nwork], &iwork[1],
1254 info, (ftnlen)1, (ftnlen)1);
1255
1256 /* Overwrite U by left singular vectors of L and WORK(IVT) */
1257 /* by right singular vectors of L */
1258 /* (Workspace: need 2*M*M+3*M, prefer 2*M*M+2*M+M*NB) */
1259
1260 i__1 = *lwork - nwork + 1;
1261 dormbr_("Q", "L", "N", m, m, m, &work[il], &ldwrkl, &work[
1262 itauq], &u[u_offset], ldu, &work[nwork], &i__1, &ierr,
1263 (ftnlen)1, (ftnlen)1, (ftnlen)1);
1264 i__1 = *lwork - nwork + 1;
1265 dormbr_("P", "R", "T", m, m, m, &work[il], &ldwrkl, &work[
1266 itaup], &work[ivt], m, &work[nwork], &i__1, &ierr, (
1267 ftnlen)1, (ftnlen)1, (ftnlen)1);
1268
1269 /* Multiply right singular vectors of L in WORK(IVT) by Q */
1270 /* in A, storing result in WORK(IL) and copying to A */
1271 /* (Workspace: need 2*M*M, prefer M*M+M*N) */
1272
1273 i__1 = *n;
1274 i__2 = chunk;
1275 for (i__ = 1; i__2 < 0 ? i__ >= i__1 : i__ <= i__1; i__ +=
1276 i__2) {
1277 /* Computing MIN */
1278 i__3 = *n - i__ + 1;
1279 blk = min(i__3,chunk);
1280 dgemm_("N", "N", m, &blk, m, &c_b248, &work[ivt], m, &a[
1281 i__ * a_dim1 + 1], lda, &c_b227, &work[il], &
1282 ldwrkl, (ftnlen)1, (ftnlen)1);
1283 dlacpy_("F", m, &blk, &work[il], &ldwrkl, &a[i__ * a_dim1
1284 + 1], lda, (ftnlen)1);
1285 /* L30: */
1286 }
1287
1288 } else if (wntqs) {
1289
1290 /* Path 3t (N much larger than M, JOBZ='S') */
1291 /* M right singular vectors to be computed in VT and */
1292 /* M left singular vectors to be computed in U */
1293
1294 il = 1;
1295
1296 /* WORK(IL) is M by M */
1297
1298 ldwrkl = *m;
1299 itau = il + ldwrkl * *m;
1300 nwork = itau + *m;
1301
1302 /* Compute A=L*Q */
1303 /* (Workspace: need M*M+2*M, prefer M*M+M+M*NB) */
1304
1305 i__2 = *lwork - nwork + 1;
1306 dgelqf_(m, n, &a[a_offset], lda, &work[itau], &work[nwork], &
1307 i__2, &ierr);
1308
1309 /* Copy L to WORK(IL), zeroing out above it */
1310
1311 dlacpy_("L", m, m, &a[a_offset], lda, &work[il], &ldwrkl, (
1312 ftnlen)1);
1313 i__2 = *m - 1;
1314 i__1 = *m - 1;
1315 dlaset_("U", &i__2, &i__1, &c_b227, &c_b227, &work[il +
1316 ldwrkl], &ldwrkl, (ftnlen)1);
1317
1318 /* Generate Q in A */
1319 /* (Workspace: need M*M+2*M, prefer M*M+M+M*NB) */
1320
1321 i__2 = *lwork - nwork + 1;
1322 dorglq_(m, n, m, &a[a_offset], lda, &work[itau], &work[nwork],
1323 &i__2, &ierr);
1324 ie = itau;
1325 itauq = ie + *m;
1326 itaup = itauq + *m;
1327 nwork = itaup + *m;
1328
1329 /* Bidiagonalize L in WORK(IU), copying result to U */
1330 /* (Workspace: need M*M+4*M, prefer M*M+3*M+2*M*NB) */
1331
1332 i__2 = *lwork - nwork + 1;
1333 dgebrd_(m, m, &work[il], &ldwrkl, &s[1], &work[ie], &work[
1334 itauq], &work[itaup], &work[nwork], &i__2, &ierr);
1335
1336 /* Perform bidiagonal SVD, computing left singular vectors */
1337 /* of bidiagonal matrix in U and computing right singular */
1338 /* vectors of bidiagonal matrix in VT */
1339 /* (Workspace: need M+BDSPAC) */
1340
1341 dbdsdc_("U", "I", m, &s[1], &work[ie], &u[u_offset], ldu, &vt[
1342 vt_offset], ldvt, dum, idum, &work[nwork], &iwork[1],
1343 info, (ftnlen)1, (ftnlen)1);
1344
1345 /* Overwrite U by left singular vectors of L and VT */
1346 /* by right singular vectors of L */
1347 /* (Workspace: need M*M+3*M, prefer M*M+2*M+M*NB) */
1348
1349 i__2 = *lwork - nwork + 1;
1350 dormbr_("Q", "L", "N", m, m, m, &work[il], &ldwrkl, &work[
1351 itauq], &u[u_offset], ldu, &work[nwork], &i__2, &ierr,
1352 (ftnlen)1, (ftnlen)1, (ftnlen)1);
1353 i__2 = *lwork - nwork + 1;
1354 dormbr_("P", "R", "T", m, m, m, &work[il], &ldwrkl, &work[
1355 itaup], &vt[vt_offset], ldvt, &work[nwork], &i__2, &
1356 ierr, (ftnlen)1, (ftnlen)1, (ftnlen)1);
1357
1358 /* Multiply right singular vectors of L in WORK(IL) by */
1359 /* Q in A, storing result in VT */
1360 /* (Workspace: need M*M) */
1361
1362 dlacpy_("F", m, m, &vt[vt_offset], ldvt, &work[il], &ldwrkl, (
1363 ftnlen)1);
1364 dgemm_("N", "N", m, n, m, &c_b248, &work[il], &ldwrkl, &a[
1365 a_offset], lda, &c_b227, &vt[vt_offset], ldvt, (
1366 ftnlen)1, (ftnlen)1);
1367
1368 } else if (wntqa) {
1369
1370 /* Path 4t (N much larger than M, JOBZ='A') */
1371 /* N right singular vectors to be computed in VT and */
1372 /* M left singular vectors to be computed in U */
1373
1374 ivt = 1;
1375
1376 /* WORK(IVT) is M by M */
1377
1378 ldwkvt = *m;
1379 itau = ivt + ldwkvt * *m;
1380 nwork = itau + *m;
1381
1382 /* Compute A=L*Q, copying result to VT */
1383 /* (Workspace: need M*M+2*M, prefer M*M+M+M*NB) */
1384
1385 i__2 = *lwork - nwork + 1;
1386 dgelqf_(m, n, &a[a_offset], lda, &work[itau], &work[nwork], &
1387 i__2, &ierr);
1388 dlacpy_("U", m, n, &a[a_offset], lda, &vt[vt_offset], ldvt, (
1389 ftnlen)1);
1390
1391 /* Generate Q in VT */
1392 /* (Workspace: need M*M+2*M, prefer M*M+M+M*NB) */
1393
1394 i__2 = *lwork - nwork + 1;
1395 dorglq_(n, n, m, &vt[vt_offset], ldvt, &work[itau], &work[
1396 nwork], &i__2, &ierr);
1397
1398 /* Produce L in A, zeroing out other entries */
1399
1400 i__2 = *m - 1;
1401 i__1 = *m - 1;
1402 dlaset_("U", &i__2, &i__1, &c_b227, &c_b227, &a[(a_dim1 << 1)
1403 + 1], lda, (ftnlen)1);
1404 ie = itau;
1405 itauq = ie + *m;
1406 itaup = itauq + *m;
1407 nwork = itaup + *m;
1408
1409 /* Bidiagonalize L in A */
1410 /* (Workspace: need M*M+4*M, prefer M*M+3*M+2*M*NB) */
1411
1412 i__2 = *lwork - nwork + 1;
1413 dgebrd_(m, m, &a[a_offset], lda, &s[1], &work[ie], &work[
1414 itauq], &work[itaup], &work[nwork], &i__2, &ierr);
1415
1416 /* Perform bidiagonal SVD, computing left singular vectors */
1417 /* of bidiagonal matrix in U and computing right singular */
1418 /* vectors of bidiagonal matrix in WORK(IVT) */
1419 /* (Workspace: need M+M*M+BDSPAC) */
1420
1421 dbdsdc_("U", "I", m, &s[1], &work[ie], &u[u_offset], ldu, &
1422 work[ivt], &ldwkvt, dum, idum, &work[nwork], &iwork[1]
1423 , info, (ftnlen)1, (ftnlen)1);
1424
1425 /* Overwrite U by left singular vectors of L and WORK(IVT) */
1426 /* by right singular vectors of L */
1427 /* (Workspace: need M*M+3*M, prefer M*M+2*M+M*NB) */
1428
1429 i__2 = *lwork - nwork + 1;
1430 dormbr_("Q", "L", "N", m, m, m, &a[a_offset], lda, &work[
1431 itauq], &u[u_offset], ldu, &work[nwork], &i__2, &ierr,
1432 (ftnlen)1, (ftnlen)1, (ftnlen)1);
1433 i__2 = *lwork - nwork + 1;
1434 dormbr_("P", "R", "T", m, m, m, &a[a_offset], lda, &work[
1435 itaup], &work[ivt], &ldwkvt, &work[nwork], &i__2, &
1436 ierr, (ftnlen)1, (ftnlen)1, (ftnlen)1);
1437
1438 /* Multiply right singular vectors of L in WORK(IVT) by */
1439 /* Q in VT, storing result in A */
1440 /* (Workspace: need M*M) */
1441
1442 dgemm_("N", "N", m, n, m, &c_b248, &work[ivt], &ldwkvt, &vt[
1443 vt_offset], ldvt, &c_b227, &a[a_offset], lda, (ftnlen)
1444 1, (ftnlen)1);
1445
1446 /* Copy right singular vectors of A from A to VT */
1447
1448 dlacpy_("F", m, n, &a[a_offset], lda, &vt[vt_offset], ldvt, (
1449 ftnlen)1);
1450
1451 }
1452
1453 } else {
1454
1455 /* N .LT. MNTHR */
1456
1457 /* Path 5t (N greater than M, but not much larger) */
1458 /* Reduce to bidiagonal form without LQ decomposition */
1459
1460 ie = 1;
1461 itauq = ie + *m;
1462 itaup = itauq + *m;
1463 nwork = itaup + *m;
1464
1465 /* Bidiagonalize A */
1466 /* (Workspace: need 3*M+N, prefer 3*M+(M+N)*NB) */
1467
1468 i__2 = *lwork - nwork + 1;
1469 dgebrd_(m, n, &a[a_offset], lda, &s[1], &work[ie], &work[itauq], &
1470 work[itaup], &work[nwork], &i__2, &ierr);
1471 if (wntqn) {
1472
1473 /* Perform bidiagonal SVD, only computing singular values */
1474 /* (Workspace: need M+BDSPAC) */
1475
1476 dbdsdc_("L", "N", m, &s[1], &work[ie], dum, &c__1, dum, &c__1,
1477 dum, idum, &work[nwork], &iwork[1], info, (ftnlen)1,
1478 (ftnlen)1);
1479 } else if (wntqo) {
1480 ldwkvt = *m;
1481 ivt = nwork;
1482 if (*lwork >= *m * *n + *m * 3 + bdspac) {
1483
1484 /* WORK( IVT ) is M by N */
1485
1486 dlaset_("F", m, n, &c_b227, &c_b227, &work[ivt], &ldwkvt,
1487 (ftnlen)1);
1488 nwork = ivt + ldwkvt * *n;
1489 } else {
1490
1491 /* WORK( IVT ) is M by M */
1492
1493 nwork = ivt + ldwkvt * *m;
1494 il = nwork;
1495
1496 /* WORK(IL) is M by CHUNK */
1497
1498 chunk = (*lwork - *m * *m - *m * 3) / *m;
1499 }
1500
1501 /* Perform bidiagonal SVD, computing left singular vectors */
1502 /* of bidiagonal matrix in U and computing right singular */
1503 /* vectors of bidiagonal matrix in WORK(IVT) */
1504 /* (Workspace: need M*M+BDSPAC) */
1505
1506 dbdsdc_("L", "I", m, &s[1], &work[ie], &u[u_offset], ldu, &
1507 work[ivt], &ldwkvt, dum, idum, &work[nwork], &iwork[1]
1508 , info, (ftnlen)1, (ftnlen)1);
1509
1510 /* Overwrite U by left singular vectors of A */
1511 /* (Workspace: need M*M+2*M, prefer M*M+M+M*NB) */
1512
1513 i__2 = *lwork - nwork + 1;
1514 dormbr_("Q", "L", "N", m, m, n, &a[a_offset], lda, &work[
1515 itauq], &u[u_offset], ldu, &work[nwork], &i__2, &ierr,
1516 (ftnlen)1, (ftnlen)1, (ftnlen)1);
1517
1518 if (*lwork >= *m * *n + *m * 3 + bdspac) {
1519
1520 /* Overwrite WORK(IVT) by left singular vectors of A */
1521 /* (Workspace: need M*M+2*M, prefer M*M+M+M*NB) */
1522
1523 i__2 = *lwork - nwork + 1;
1524 dormbr_("P", "R", "T", m, n, m, &a[a_offset], lda, &work[
1525 itaup], &work[ivt], &ldwkvt, &work[nwork], &i__2,
1526 &ierr, (ftnlen)1, (ftnlen)1, (ftnlen)1);
1527
1528 /* Copy right singular vectors of A from WORK(IVT) to A */
1529
1530 dlacpy_("F", m, n, &work[ivt], &ldwkvt, &a[a_offset], lda,
1531 (ftnlen)1);
1532 } else {
1533
1534 /* Generate P**T in A */
1535 /* (Workspace: need M*M+2*M, prefer M*M+M+M*NB) */
1536
1537 i__2 = *lwork - nwork + 1;
1538 dorgbr_("P", m, n, m, &a[a_offset], lda, &work[itaup], &
1539 work[nwork], &i__2, &ierr, (ftnlen)1);
1540
1541 /* Multiply Q in A by right singular vectors of */
1542 /* bidiagonal matrix in WORK(IVT), storing result in */
1543 /* WORK(IL) and copying to A */
1544 /* (Workspace: need 2*M*M, prefer M*M+M*N) */
1545
1546 i__2 = *n;
1547 i__1 = chunk;
1548 for (i__ = 1; i__1 < 0 ? i__ >= i__2 : i__ <= i__2; i__ +=
1549 i__1) {
1550 /* Computing MIN */
1551 i__3 = *n - i__ + 1;
1552 blk = min(i__3,chunk);
1553 dgemm_("N", "N", m, &blk, m, &c_b248, &work[ivt], &
1554 ldwkvt, &a[i__ * a_dim1 + 1], lda, &c_b227, &
1555 work[il], m, (ftnlen)1, (ftnlen)1);
1556 dlacpy_("F", m, &blk, &work[il], m, &a[i__ * a_dim1 +
1557 1], lda, (ftnlen)1);
1558 /* L40: */
1559 }
1560 }
1561 } else if (wntqs) {
1562
1563 /* Perform bidiagonal SVD, computing left singular vectors */
1564 /* of bidiagonal matrix in U and computing right singular */
1565 /* vectors of bidiagonal matrix in VT */
1566 /* (Workspace: need M+BDSPAC) */
1567
1568 dlaset_("F", m, n, &c_b227, &c_b227, &vt[vt_offset], ldvt, (
1569 ftnlen)1);
1570 dbdsdc_("L", "I", m, &s[1], &work[ie], &u[u_offset], ldu, &vt[
1571 vt_offset], ldvt, dum, idum, &work[nwork], &iwork[1],
1572 info, (ftnlen)1, (ftnlen)1);
1573
1574 /* Overwrite U by left singular vectors of A and VT */
1575 /* by right singular vectors of A */
1576 /* (Workspace: need 3*M, prefer 2*M+M*NB) */
1577
1578 i__1 = *lwork - nwork + 1;
1579 dormbr_("Q", "L", "N", m, m, n, &a[a_offset], lda, &work[
1580 itauq], &u[u_offset], ldu, &work[nwork], &i__1, &ierr,
1581 (ftnlen)1, (ftnlen)1, (ftnlen)1);
1582 i__1 = *lwork - nwork + 1;
1583 dormbr_("P", "R", "T", m, n, m, &a[a_offset], lda, &work[
1584 itaup], &vt[vt_offset], ldvt, &work[nwork], &i__1, &
1585 ierr, (ftnlen)1, (ftnlen)1, (ftnlen)1);
1586 } else if (wntqa) {
1587
1588 /* Perform bidiagonal SVD, computing left singular vectors */
1589 /* of bidiagonal matrix in U and computing right singular */
1590 /* vectors of bidiagonal matrix in VT */
1591 /* (Workspace: need M+BDSPAC) */
1592
1593 dlaset_("F", n, n, &c_b227, &c_b227, &vt[vt_offset], ldvt, (
1594 ftnlen)1);
1595 dbdsdc_("L", "I", m, &s[1], &work[ie], &u[u_offset], ldu, &vt[
1596 vt_offset], ldvt, dum, idum, &work[nwork], &iwork[1],
1597 info, (ftnlen)1, (ftnlen)1);
1598
1599 /* Set the right corner of VT to identity matrix */
1600
1601 i__1 = *n - *m;
1602 i__2 = *n - *m;
1603 dlaset_("F", &i__1, &i__2, &c_b227, &c_b248, &vt[*m + 1 + (*m
1604 + 1) * vt_dim1], ldvt, (ftnlen)1);
1605
1606 /* Overwrite U by left singular vectors of A and VT */
1607 /* by right singular vectors of A */
1608 /* (Workspace: need 2*M+N, prefer 2*M+N*NB) */
1609
1610 i__1 = *lwork - nwork + 1;
1611 dormbr_("Q", "L", "N", m, m, n, &a[a_offset], lda, &work[
1612 itauq], &u[u_offset], ldu, &work[nwork], &i__1, &ierr,
1613 (ftnlen)1, (ftnlen)1, (ftnlen)1);
1614 i__1 = *lwork - nwork + 1;
1615 dormbr_("P", "R", "T", n, n, m, &a[a_offset], lda, &work[
1616 itaup], &vt[vt_offset], ldvt, &work[nwork], &i__1, &
1617 ierr, (ftnlen)1, (ftnlen)1, (ftnlen)1);
1618 }
1619
1620 }
1621
1622 }
1623
1624 /* Undo scaling if necessary */
1625
1626 if (iscl == 1) {
1627 if (anrm > bignum) {
1628 dlascl_("G", &c__0, &c__0, &bignum, &anrm, &minmn, &c__1, &s[1], &
1629 minmn, &ierr, (ftnlen)1);
1630 }
1631 if (anrm < smlnum) {
1632 dlascl_("G", &c__0, &c__0, &smlnum, &anrm, &minmn, &c__1, &s[1], &
1633 minmn, &ierr, (ftnlen)1);
1634 }
1635 }
1636
1637 /* Return optimal workspace in WORK(1) */
1638
1639 work[1] = (doublereal) maxwrk;
1640
1641 return 0;
1642
1643 /* End of DGESDD */
1644
1645 } /* dgesdd_ */
1646
1647