2 /* -- translated by f2c (version 19990503).
3 You must link the resulting object file with the libraries:
4 -lf2c -lm (in that order)
9 /* Table of constant values */
11 static doublecomplex c_b1 = {0.,0.};
12 static doublecomplex c_b2 = {1.,0.};
13 static integer c__6 = 6;
14 static integer c__0 = 0;
15 static integer c__2 = 2;
16 static integer c__1 = 1;
17 static integer c_n1 = -1;
19 /* Subroutine */ int zgesvd_(char *jobu, char *jobvt, integer *m, integer *n,
20 doublecomplex *a, integer *lda, doublereal *s, doublecomplex *u,
21 integer *ldu, doublecomplex *vt, integer *ldvt, doublecomplex *work,
22 integer *lwork, doublereal *rwork, integer *info)
24 /* System generated locals */
26 integer a_dim1, a_offset, u_dim1, u_offset, vt_dim1, vt_offset, i__1[2],
31 Subroutine */ int s_cat(char *, char **, integer *, integer *, ftnlen);
32 double sqrt(doublereal);
35 static doublecomplex cdum[1];
37 static doublereal anrm;
38 static integer ierr, itau, ncvt, nrvt, i__;
39 extern logical lsame_(char *, char *);
40 static integer chunk, minmn;
41 extern /* Subroutine */ int zgemm_(char *, char *, integer *, integer *,
42 integer *, doublecomplex *, doublecomplex *, integer *,
43 doublecomplex *, integer *, doublecomplex *, doublecomplex *,
45 static integer wrkbl, itaup, itauq, mnthr, iwork;
46 static logical wntua, wntva, wntun, wntuo, wntvn, wntvo, wntus, wntvs;
48 extern doublereal dlamch_(char *);
49 static integer ir, iu;
50 extern /* Subroutine */ int dlascl_(char *, integer *, integer *,
51 doublereal *, doublereal *, integer *, integer *, doublereal *,
52 integer *, integer *), xerbla_(char *, integer *),
53 zgebrd_(integer *, integer *, doublecomplex *, integer *,
54 doublereal *, doublereal *, doublecomplex *, doublecomplex *,
55 doublecomplex *, integer *, integer *);
56 extern integer ilaenv_(integer *, char *, char *, integer *, integer *,
57 integer *, integer *, ftnlen, ftnlen);
58 extern doublereal zlange_(char *, integer *, integer *, doublecomplex *,
59 integer *, doublereal *);
60 static doublereal bignum;
61 extern /* Subroutine */ int zgelqf_(integer *, integer *, doublecomplex *,
62 integer *, doublecomplex *, doublecomplex *, integer *, integer *
63 ), zlascl_(char *, integer *, integer *, doublereal *, doublereal
64 *, integer *, integer *, doublecomplex *, integer *, integer *), zgeqrf_(integer *, integer *, doublecomplex *, integer *,
65 doublecomplex *, doublecomplex *, integer *, integer *), zlacpy_(
66 char *, integer *, integer *, doublecomplex *, integer *,
67 doublecomplex *, integer *), zlaset_(char *, integer *,
68 integer *, doublecomplex *, doublecomplex *, doublecomplex *,
70 static integer ldwrkr;
71 extern /* Subroutine */ int zbdsqr_(char *, integer *, integer *, integer
72 *, integer *, doublereal *, doublereal *, doublecomplex *,
73 integer *, doublecomplex *, integer *, doublecomplex *, integer *,
74 doublereal *, integer *);
75 static integer minwrk, ldwrku, maxwrk;
76 extern /* Subroutine */ int zungbr_(char *, integer *, integer *, integer
77 *, doublecomplex *, integer *, doublecomplex *, doublecomplex *,
78 integer *, integer *);
79 static doublereal smlnum;
80 static integer irwork;
81 extern /* Subroutine */ int zunmbr_(char *, char *, char *, integer *,
82 integer *, integer *, doublecomplex *, integer *, doublecomplex *,
83 doublecomplex *, integer *, doublecomplex *, integer *, integer *
84 ), zunglq_(integer *, integer *, integer *
85 , doublecomplex *, integer *, doublecomplex *, doublecomplex *,
86 integer *, integer *);
87 static logical lquery, wntuas, wntvas;
88 extern /* Subroutine */ int zungqr_(integer *, integer *, integer *,
89 doublecomplex *, integer *, doublecomplex *, doublecomplex *,
90 integer *, integer *);
91 static integer blk, ncu;
92 static doublereal dum[1], eps;
96 #define a_subscr(a_1,a_2) (a_2)*a_dim1 + a_1
97 #define a_ref(a_1,a_2) a[a_subscr(a_1,a_2)]
98 #define u_subscr(a_1,a_2) (a_2)*u_dim1 + a_1
99 #define u_ref(a_1,a_2) u[u_subscr(a_1,a_2)]
100 #define vt_subscr(a_1,a_2) (a_2)*vt_dim1 + a_1
101 #define vt_ref(a_1,a_2) vt[vt_subscr(a_1,a_2)]
104 /* -- LAPACK driver routine (version 3.0) --
105 Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
106 Courant Institute, Argonne National Lab, and Rice University
113 ZGESVD computes the singular value decomposition (SVD) of a complex
114 M-by-N matrix A, optionally computing the left and/or right singular
115 vectors. The SVD is written
117 A = U * SIGMA * conjugate-transpose(V)
119 where SIGMA is an M-by-N matrix which is zero except for its
120 min(m,n) diagonal elements, U is an M-by-M unitary matrix, and
121 V is an N-by-N unitary matrix. The diagonal elements of SIGMA
122 are the singular values of A; they are real and non-negative, and
123 are returned in descending order. The first min(m,n) columns of
124 U and V are the left and right singular vectors of A.
126 Note that the routine returns V**H, not V.
131 JOBU (input) CHARACTER*1
132 Specifies options for computing all or part of the matrix U:
133 = 'A': all M columns of U are returned in array U:
134 = 'S': the first min(m,n) columns of U (the left singular
135 vectors) are returned in the array U;
136 = 'O': the first min(m,n) columns of U (the left singular
137 vectors) are overwritten on the array A;
138 = 'N': no columns of U (no left singular vectors) are
141 JOBVT (input) CHARACTER*1
142 Specifies options for computing all or part of the matrix
144 = 'A': all N rows of V**H are returned in the array VT;
145 = 'S': the first min(m,n) rows of V**H (the right singular
146 vectors) are returned in the array VT;
147 = 'O': the first min(m,n) rows of V**H (the right singular
148 vectors) are overwritten on the array A;
149 = 'N': no rows of V**H (no right singular vectors) are
152 JOBVT and JOBU cannot both be 'O'.
155 The number of rows of the input matrix A. M >= 0.
158 The number of columns of the input matrix A. N >= 0.
160 A (input/output) COMPLEX*16 array, dimension (LDA,N)
161 On entry, the M-by-N matrix A.
163 if JOBU = 'O', A is overwritten with the first min(m,n)
164 columns of U (the left singular vectors,
166 if JOBVT = 'O', A is overwritten with the first min(m,n)
167 rows of V**H (the right singular vectors,
169 if JOBU .ne. 'O' and JOBVT .ne. 'O', the contents of A
173 The leading dimension of the array A. LDA >= max(1,M).
175 S (output) DOUBLE PRECISION array, dimension (min(M,N))
176 The singular values of A, sorted so that S(i) >= S(i+1).
178 U (output) COMPLEX*16 array, dimension (LDU,UCOL)
179 (LDU,M) if JOBU = 'A' or (LDU,min(M,N)) if JOBU = 'S'.
180 If JOBU = 'A', U contains the M-by-M unitary matrix U;
181 if JOBU = 'S', U contains the first min(m,n) columns of U
182 (the left singular vectors, stored columnwise);
183 if JOBU = 'N' or 'O', U is not referenced.
186 The leading dimension of the array U. LDU >= 1; if
187 JOBU = 'S' or 'A', LDU >= M.
189 VT (output) COMPLEX*16 array, dimension (LDVT,N)
190 If JOBVT = 'A', VT contains the N-by-N unitary matrix
192 if JOBVT = 'S', VT contains the first min(m,n) rows of
193 V**H (the right singular vectors, stored rowwise);
194 if JOBVT = 'N' or 'O', VT is not referenced.
197 The leading dimension of the array VT. LDVT >= 1; if
198 JOBVT = 'A', LDVT >= N; if JOBVT = 'S', LDVT >= min(M,N).
200 WORK (workspace/output) COMPLEX*16 array, dimension (LWORK)
201 On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
203 LWORK (input) INTEGER
204 The dimension of the array WORK. LWORK >= 1.
205 LWORK >= 2*MIN(M,N)+MAX(M,N).
206 For good performance, LWORK should generally be larger.
208 If LWORK = -1, then a workspace query is assumed; the routine
209 only calculates the optimal size of the WORK array, returns
210 this value as the first entry of the WORK array, and no error
211 message related to LWORK is issued by XERBLA.
213 RWORK (workspace) DOUBLE PRECISION array, dimension (5*min(M,N))
214 On exit, if INFO > 0, RWORK(1:MIN(M,N)-1) contains the
215 unconverged superdiagonal elements of an upper bidiagonal
216 matrix B whose diagonal is in S (not necessarily sorted).
217 B satisfies A = U * B * VT, so it has the same singular
218 values as A, and singular vectors related by U and VT.
220 INFO (output) INTEGER
221 = 0: successful exit.
222 < 0: if INFO = -i, the i-th argument had an illegal value.
223 > 0: if ZBDSQR did not converge, INFO specifies how many
224 superdiagonals of an intermediate bidiagonal form B
225 did not converge to zero. See the description of RWORK
228 =====================================================================
231 Test the input arguments
233 Parameter adjustments */
235 a_offset = 1 + a_dim1 * 1;
239 u_offset = 1 + u_dim1 * 1;
242 vt_offset = 1 + vt_dim1 * 1;
250 /* Writing concatenation */
251 i__1[0] = 1, a__1[0] = jobu;
252 i__1[1] = 1, a__1[1] = jobvt;
253 s_cat(ch__1, a__1, i__1, &c__2, (ftnlen)2);
254 mnthr = ilaenv_(&c__6, "ZGESVD", ch__1, m, n, &c__0, &c__0, (ftnlen)6, (
256 wntua = lsame_(jobu, "A");
257 wntus = lsame_(jobu, "S");
258 wntuas = wntua || wntus;
259 wntuo = lsame_(jobu, "O");
260 wntun = lsame_(jobu, "N");
261 wntva = lsame_(jobvt, "A");
262 wntvs = lsame_(jobvt, "S");
263 wntvas = wntva || wntvs;
264 wntvo = lsame_(jobvt, "O");
265 wntvn = lsame_(jobvt, "N");
267 lquery = *lwork == -1;
269 if (! (wntua || wntus || wntuo || wntun)) {
271 } else if (! (wntva || wntvs || wntvo || wntvn) || wntvo && wntuo) {
277 } else if (*lda < max(1,*m)) {
279 } else if (*ldu < 1 || wntuas && *ldu < *m) {
281 } else if (*ldvt < 1 || wntva && *ldvt < *n || wntvs && *ldvt < minmn) {
286 (Note: Comments in the code beginning "Workspace:" describe the
287 minimal amount of workspace needed at that point in the code,
288 as well as the preferred amount for good performance.
289 CWorkspace refers to complex workspace, and RWorkspace to
290 real workspace. NB refers to the optimal block size for the
291 immediately following subroutine, as returned by ILAENV.) */
293 if (*info == 0 && (*lwork >= 1 || lquery) && *m > 0 && *n > 0) {
296 /* Space needed for ZBDSQR is BDSPAC = 5*N */
301 /* Path 1 (M much larger than N, JOBU='N') */
303 maxwrk = *n + *n * ilaenv_(&c__1, "ZGEQRF", " ", m, n, &
304 c_n1, &c_n1, (ftnlen)6, (ftnlen)1);
306 i__2 = maxwrk, i__3 = (*n << 1) + (*n << 1) * ilaenv_(&
307 c__1, "ZGEBRD", " ", n, n, &c_n1, &c_n1, (ftnlen)
309 maxwrk = max(i__2,i__3);
310 if (wntvo || wntvas) {
312 i__2 = maxwrk, i__3 = (*n << 1) + (*n - 1) * ilaenv_(&
313 c__1, "ZUNGBR", "P", n, n, n, &c_n1, (ftnlen)
315 maxwrk = max(i__2,i__3);
318 maxwrk = max(minwrk,maxwrk);
319 } else if (wntuo && wntvn) {
321 /* Path 2 (M much larger than N, JOBU='O', JOBVT='N') */
323 wrkbl = *n + *n * ilaenv_(&c__1, "ZGEQRF", " ", m, n, &
324 c_n1, &c_n1, (ftnlen)6, (ftnlen)1);
326 i__2 = wrkbl, i__3 = *n + *n * ilaenv_(&c__1, "ZUNGQR",
327 " ", m, n, n, &c_n1, (ftnlen)6, (ftnlen)1);
328 wrkbl = max(i__2,i__3);
330 i__2 = wrkbl, i__3 = (*n << 1) + (*n << 1) * ilaenv_(&
331 c__1, "ZGEBRD", " ", n, n, &c_n1, &c_n1, (ftnlen)
333 wrkbl = max(i__2,i__3);
335 i__2 = wrkbl, i__3 = (*n << 1) + *n * ilaenv_(&c__1,
336 "ZUNGBR", "Q", n, n, n, &c_n1, (ftnlen)6, (ftnlen)
338 wrkbl = max(i__2,i__3);
340 i__2 = *n * *n + wrkbl, i__3 = *n * *n + *m * *n;
341 maxwrk = max(i__2,i__3);
342 minwrk = (*n << 1) + *m;
343 maxwrk = max(minwrk,maxwrk);
344 } else if (wntuo && wntvas) {
346 /* Path 3 (M much larger than N, JOBU='O', JOBVT='S' or
349 wrkbl = *n + *n * ilaenv_(&c__1, "ZGEQRF", " ", m, n, &
350 c_n1, &c_n1, (ftnlen)6, (ftnlen)1);
352 i__2 = wrkbl, i__3 = *n + *n * ilaenv_(&c__1, "ZUNGQR",
353 " ", m, n, n, &c_n1, (ftnlen)6, (ftnlen)1);
354 wrkbl = max(i__2,i__3);
356 i__2 = wrkbl, i__3 = (*n << 1) + (*n << 1) * ilaenv_(&
357 c__1, "ZGEBRD", " ", n, n, &c_n1, &c_n1, (ftnlen)
359 wrkbl = max(i__2,i__3);
361 i__2 = wrkbl, i__3 = (*n << 1) + *n * ilaenv_(&c__1,
362 "ZUNGBR", "Q", n, n, n, &c_n1, (ftnlen)6, (ftnlen)
364 wrkbl = max(i__2,i__3);
366 i__2 = wrkbl, i__3 = (*n << 1) + (*n - 1) * ilaenv_(&c__1,
367 "ZUNGBR", "P", n, n, n, &c_n1, (ftnlen)6, (
369 wrkbl = max(i__2,i__3);
371 i__2 = *n * *n + wrkbl, i__3 = *n * *n + *m * *n;
372 maxwrk = max(i__2,i__3);
373 minwrk = (*n << 1) + *m;
374 maxwrk = max(minwrk,maxwrk);
375 } else if (wntus && wntvn) {
377 /* Path 4 (M much larger than N, JOBU='S', JOBVT='N') */
379 wrkbl = *n + *n * ilaenv_(&c__1, "ZGEQRF", " ", m, n, &
380 c_n1, &c_n1, (ftnlen)6, (ftnlen)1);
382 i__2 = wrkbl, i__3 = *n + *n * ilaenv_(&c__1, "ZUNGQR",
383 " ", m, n, n, &c_n1, (ftnlen)6, (ftnlen)1);
384 wrkbl = max(i__2,i__3);
386 i__2 = wrkbl, i__3 = (*n << 1) + (*n << 1) * ilaenv_(&
387 c__1, "ZGEBRD", " ", n, n, &c_n1, &c_n1, (ftnlen)
389 wrkbl = max(i__2,i__3);
391 i__2 = wrkbl, i__3 = (*n << 1) + *n * ilaenv_(&c__1,
392 "ZUNGBR", "Q", n, n, n, &c_n1, (ftnlen)6, (ftnlen)
394 wrkbl = max(i__2,i__3);
395 maxwrk = *n * *n + wrkbl;
396 minwrk = (*n << 1) + *m;
397 maxwrk = max(minwrk,maxwrk);
398 } else if (wntus && wntvo) {
400 /* Path 5 (M much larger than N, JOBU='S', JOBVT='O') */
402 wrkbl = *n + *n * ilaenv_(&c__1, "ZGEQRF", " ", m, n, &
403 c_n1, &c_n1, (ftnlen)6, (ftnlen)1);
405 i__2 = wrkbl, i__3 = *n + *n * ilaenv_(&c__1, "ZUNGQR",
406 " ", m, n, n, &c_n1, (ftnlen)6, (ftnlen)1);
407 wrkbl = max(i__2,i__3);
409 i__2 = wrkbl, i__3 = (*n << 1) + (*n << 1) * ilaenv_(&
410 c__1, "ZGEBRD", " ", n, n, &c_n1, &c_n1, (ftnlen)
412 wrkbl = max(i__2,i__3);
414 i__2 = wrkbl, i__3 = (*n << 1) + *n * ilaenv_(&c__1,
415 "ZUNGBR", "Q", n, n, n, &c_n1, (ftnlen)6, (ftnlen)
417 wrkbl = max(i__2,i__3);
419 i__2 = wrkbl, i__3 = (*n << 1) + (*n - 1) * ilaenv_(&c__1,
420 "ZUNGBR", "P", n, n, n, &c_n1, (ftnlen)6, (
422 wrkbl = max(i__2,i__3);
423 maxwrk = (*n << 1) * *n + wrkbl;
424 minwrk = (*n << 1) + *m;
425 maxwrk = max(minwrk,maxwrk);
426 } else if (wntus && wntvas) {
428 /* Path 6 (M much larger than N, JOBU='S', JOBVT='S' or
431 wrkbl = *n + *n * ilaenv_(&c__1, "ZGEQRF", " ", m, n, &
432 c_n1, &c_n1, (ftnlen)6, (ftnlen)1);
434 i__2 = wrkbl, i__3 = *n + *n * ilaenv_(&c__1, "ZUNGQR",
435 " ", m, n, n, &c_n1, (ftnlen)6, (ftnlen)1);
436 wrkbl = max(i__2,i__3);
438 i__2 = wrkbl, i__3 = (*n << 1) + (*n << 1) * ilaenv_(&
439 c__1, "ZGEBRD", " ", n, n, &c_n1, &c_n1, (ftnlen)
441 wrkbl = max(i__2,i__3);
443 i__2 = wrkbl, i__3 = (*n << 1) + *n * ilaenv_(&c__1,
444 "ZUNGBR", "Q", n, n, n, &c_n1, (ftnlen)6, (ftnlen)
446 wrkbl = max(i__2,i__3);
448 i__2 = wrkbl, i__3 = (*n << 1) + (*n - 1) * ilaenv_(&c__1,
449 "ZUNGBR", "P", n, n, n, &c_n1, (ftnlen)6, (
451 wrkbl = max(i__2,i__3);
452 maxwrk = *n * *n + wrkbl;
453 minwrk = (*n << 1) + *m;
454 maxwrk = max(minwrk,maxwrk);
455 } else if (wntua && wntvn) {
457 /* Path 7 (M much larger than N, JOBU='A', JOBVT='N') */
459 wrkbl = *n + *n * ilaenv_(&c__1, "ZGEQRF", " ", m, n, &
460 c_n1, &c_n1, (ftnlen)6, (ftnlen)1);
462 i__2 = wrkbl, i__3 = *n + *m * ilaenv_(&c__1, "ZUNGQR",
463 " ", m, m, n, &c_n1, (ftnlen)6, (ftnlen)1);
464 wrkbl = max(i__2,i__3);
466 i__2 = wrkbl, i__3 = (*n << 1) + (*n << 1) * ilaenv_(&
467 c__1, "ZGEBRD", " ", n, n, &c_n1, &c_n1, (ftnlen)
469 wrkbl = max(i__2,i__3);
471 i__2 = wrkbl, i__3 = (*n << 1) + *n * ilaenv_(&c__1,
472 "ZUNGBR", "Q", n, n, n, &c_n1, (ftnlen)6, (ftnlen)
474 wrkbl = max(i__2,i__3);
475 maxwrk = *n * *n + wrkbl;
476 minwrk = (*n << 1) + *m;
477 maxwrk = max(minwrk,maxwrk);
478 } else if (wntua && wntvo) {
480 /* Path 8 (M much larger than N, JOBU='A', JOBVT='O') */
482 wrkbl = *n + *n * ilaenv_(&c__1, "ZGEQRF", " ", m, n, &
483 c_n1, &c_n1, (ftnlen)6, (ftnlen)1);
485 i__2 = wrkbl, i__3 = *n + *m * ilaenv_(&c__1, "ZUNGQR",
486 " ", m, m, n, &c_n1, (ftnlen)6, (ftnlen)1);
487 wrkbl = max(i__2,i__3);
489 i__2 = wrkbl, i__3 = (*n << 1) + (*n << 1) * ilaenv_(&
490 c__1, "ZGEBRD", " ", n, n, &c_n1, &c_n1, (ftnlen)
492 wrkbl = max(i__2,i__3);
494 i__2 = wrkbl, i__3 = (*n << 1) + *n * ilaenv_(&c__1,
495 "ZUNGBR", "Q", n, n, n, &c_n1, (ftnlen)6, (ftnlen)
497 wrkbl = max(i__2,i__3);
499 i__2 = wrkbl, i__3 = (*n << 1) + (*n - 1) * ilaenv_(&c__1,
500 "ZUNGBR", "P", n, n, n, &c_n1, (ftnlen)6, (
502 wrkbl = max(i__2,i__3);
503 maxwrk = (*n << 1) * *n + wrkbl;
504 minwrk = (*n << 1) + *m;
505 maxwrk = max(minwrk,maxwrk);
506 } else if (wntua && wntvas) {
508 /* Path 9 (M much larger than N, JOBU='A', JOBVT='S' or
511 wrkbl = *n + *n * ilaenv_(&c__1, "ZGEQRF", " ", m, n, &
512 c_n1, &c_n1, (ftnlen)6, (ftnlen)1);
514 i__2 = wrkbl, i__3 = *n + *m * ilaenv_(&c__1, "ZUNGQR",
515 " ", m, m, n, &c_n1, (ftnlen)6, (ftnlen)1);
516 wrkbl = max(i__2,i__3);
518 i__2 = wrkbl, i__3 = (*n << 1) + (*n << 1) * ilaenv_(&
519 c__1, "ZGEBRD", " ", n, n, &c_n1, &c_n1, (ftnlen)
521 wrkbl = max(i__2,i__3);
523 i__2 = wrkbl, i__3 = (*n << 1) + *n * ilaenv_(&c__1,
524 "ZUNGBR", "Q", n, n, n, &c_n1, (ftnlen)6, (ftnlen)
526 wrkbl = max(i__2,i__3);
528 i__2 = wrkbl, i__3 = (*n << 1) + (*n - 1) * ilaenv_(&c__1,
529 "ZUNGBR", "P", n, n, n, &c_n1, (ftnlen)6, (
531 wrkbl = max(i__2,i__3);
532 maxwrk = *n * *n + wrkbl;
533 minwrk = (*n << 1) + *m;
534 maxwrk = max(minwrk,maxwrk);
538 /* Path 10 (M at least N, but not much larger) */
540 maxwrk = (*n << 1) + (*m + *n) * ilaenv_(&c__1, "ZGEBRD",
541 " ", m, n, &c_n1, &c_n1, (ftnlen)6, (ftnlen)1);
542 if (wntus || wntuo) {
544 i__2 = maxwrk, i__3 = (*n << 1) + *n * ilaenv_(&c__1,
545 "ZUNGBR", "Q", m, n, n, &c_n1, (ftnlen)6, (ftnlen)
547 maxwrk = max(i__2,i__3);
551 i__2 = maxwrk, i__3 = (*n << 1) + *m * ilaenv_(&c__1,
552 "ZUNGBR", "Q", m, m, n, &c_n1, (ftnlen)6, (ftnlen)
554 maxwrk = max(i__2,i__3);
558 i__2 = maxwrk, i__3 = (*n << 1) + (*n - 1) * ilaenv_(&
559 c__1, "ZUNGBR", "P", n, n, n, &c_n1, (ftnlen)6, (
561 maxwrk = max(i__2,i__3);
563 minwrk = (*n << 1) + *m;
564 maxwrk = max(minwrk,maxwrk);
568 /* Space needed for ZBDSQR is BDSPAC = 5*M */
573 /* Path 1t(N much larger than M, JOBVT='N') */
575 maxwrk = *m + *m * ilaenv_(&c__1, "ZGELQF", " ", m, n, &
576 c_n1, &c_n1, (ftnlen)6, (ftnlen)1);
578 i__2 = maxwrk, i__3 = (*m << 1) + (*m << 1) * ilaenv_(&
579 c__1, "ZGEBRD", " ", m, m, &c_n1, &c_n1, (ftnlen)
581 maxwrk = max(i__2,i__3);
582 if (wntuo || wntuas) {
584 i__2 = maxwrk, i__3 = (*m << 1) + *m * ilaenv_(&c__1,
585 "ZUNGBR", "Q", m, m, m, &c_n1, (ftnlen)6, (
587 maxwrk = max(i__2,i__3);
590 maxwrk = max(minwrk,maxwrk);
591 } else if (wntvo && wntun) {
593 /* Path 2t(N much larger than M, JOBU='N', JOBVT='O') */
595 wrkbl = *m + *m * ilaenv_(&c__1, "ZGELQF", " ", m, n, &
596 c_n1, &c_n1, (ftnlen)6, (ftnlen)1);
598 i__2 = wrkbl, i__3 = *m + *m * ilaenv_(&c__1, "ZUNGLQ",
599 " ", m, n, m, &c_n1, (ftnlen)6, (ftnlen)1);
600 wrkbl = max(i__2,i__3);
602 i__2 = wrkbl, i__3 = (*m << 1) + (*m << 1) * ilaenv_(&
603 c__1, "ZGEBRD", " ", m, m, &c_n1, &c_n1, (ftnlen)
605 wrkbl = max(i__2,i__3);
607 i__2 = wrkbl, i__3 = (*m << 1) + (*m - 1) * ilaenv_(&c__1,
608 "ZUNGBR", "P", m, m, m, &c_n1, (ftnlen)6, (
610 wrkbl = max(i__2,i__3);
612 i__2 = *m * *m + wrkbl, i__3 = *m * *m + *m * *n;
613 maxwrk = max(i__2,i__3);
614 minwrk = (*m << 1) + *n;
615 maxwrk = max(minwrk,maxwrk);
616 } else if (wntvo && wntuas) {
618 /* Path 3t(N much larger than M, JOBU='S' or 'A',
621 wrkbl = *m + *m * ilaenv_(&c__1, "ZGELQF", " ", m, n, &
622 c_n1, &c_n1, (ftnlen)6, (ftnlen)1);
624 i__2 = wrkbl, i__3 = *m + *m * ilaenv_(&c__1, "ZUNGLQ",
625 " ", m, n, m, &c_n1, (ftnlen)6, (ftnlen)1);
626 wrkbl = max(i__2,i__3);
628 i__2 = wrkbl, i__3 = (*m << 1) + (*m << 1) * ilaenv_(&
629 c__1, "ZGEBRD", " ", m, m, &c_n1, &c_n1, (ftnlen)
631 wrkbl = max(i__2,i__3);
633 i__2 = wrkbl, i__3 = (*m << 1) + (*m - 1) * ilaenv_(&c__1,
634 "ZUNGBR", "P", m, m, m, &c_n1, (ftnlen)6, (
636 wrkbl = max(i__2,i__3);
638 i__2 = wrkbl, i__3 = (*m << 1) + *m * ilaenv_(&c__1,
639 "ZUNGBR", "Q", m, m, m, &c_n1, (ftnlen)6, (ftnlen)
641 wrkbl = max(i__2,i__3);
643 i__2 = *m * *m + wrkbl, i__3 = *m * *m + *m * *n;
644 maxwrk = max(i__2,i__3);
645 minwrk = (*m << 1) + *n;
646 maxwrk = max(minwrk,maxwrk);
647 } else if (wntvs && wntun) {
649 /* Path 4t(N much larger than M, JOBU='N', JOBVT='S') */
651 wrkbl = *m + *m * ilaenv_(&c__1, "ZGELQF", " ", m, n, &
652 c_n1, &c_n1, (ftnlen)6, (ftnlen)1);
654 i__2 = wrkbl, i__3 = *m + *m * ilaenv_(&c__1, "ZUNGLQ",
655 " ", m, n, m, &c_n1, (ftnlen)6, (ftnlen)1);
656 wrkbl = max(i__2,i__3);
658 i__2 = wrkbl, i__3 = (*m << 1) + (*m << 1) * ilaenv_(&
659 c__1, "ZGEBRD", " ", m, m, &c_n1, &c_n1, (ftnlen)
661 wrkbl = max(i__2,i__3);
663 i__2 = wrkbl, i__3 = (*m << 1) + (*m - 1) * ilaenv_(&c__1,
664 "ZUNGBR", "P", m, m, m, &c_n1, (ftnlen)6, (
666 wrkbl = max(i__2,i__3);
667 maxwrk = *m * *m + wrkbl;
668 minwrk = (*m << 1) + *n;
669 maxwrk = max(minwrk,maxwrk);
670 } else if (wntvs && wntuo) {
672 /* Path 5t(N much larger than M, JOBU='O', JOBVT='S') */
674 wrkbl = *m + *m * ilaenv_(&c__1, "ZGELQF", " ", m, n, &
675 c_n1, &c_n1, (ftnlen)6, (ftnlen)1);
677 i__2 = wrkbl, i__3 = *m + *m * ilaenv_(&c__1, "ZUNGLQ",
678 " ", m, n, m, &c_n1, (ftnlen)6, (ftnlen)1);
679 wrkbl = max(i__2,i__3);
681 i__2 = wrkbl, i__3 = (*m << 1) + (*m << 1) * ilaenv_(&
682 c__1, "ZGEBRD", " ", m, m, &c_n1, &c_n1, (ftnlen)
684 wrkbl = max(i__2,i__3);
686 i__2 = wrkbl, i__3 = (*m << 1) + (*m - 1) * ilaenv_(&c__1,
687 "ZUNGBR", "P", m, m, m, &c_n1, (ftnlen)6, (
689 wrkbl = max(i__2,i__3);
691 i__2 = wrkbl, i__3 = (*m << 1) + *m * ilaenv_(&c__1,
692 "ZUNGBR", "Q", m, m, m, &c_n1, (ftnlen)6, (ftnlen)
694 wrkbl = max(i__2,i__3);
695 maxwrk = (*m << 1) * *m + wrkbl;
696 minwrk = (*m << 1) + *n;
697 maxwrk = max(minwrk,maxwrk);
698 } else if (wntvs && wntuas) {
700 /* Path 6t(N much larger than M, JOBU='S' or 'A',
703 wrkbl = *m + *m * ilaenv_(&c__1, "ZGELQF", " ", m, n, &
704 c_n1, &c_n1, (ftnlen)6, (ftnlen)1);
706 i__2 = wrkbl, i__3 = *m + *m * ilaenv_(&c__1, "ZUNGLQ",
707 " ", m, n, m, &c_n1, (ftnlen)6, (ftnlen)1);
708 wrkbl = max(i__2,i__3);
710 i__2 = wrkbl, i__3 = (*m << 1) + (*m << 1) * ilaenv_(&
711 c__1, "ZGEBRD", " ", m, m, &c_n1, &c_n1, (ftnlen)
713 wrkbl = max(i__2,i__3);
715 i__2 = wrkbl, i__3 = (*m << 1) + (*m - 1) * ilaenv_(&c__1,
716 "ZUNGBR", "P", m, m, m, &c_n1, (ftnlen)6, (
718 wrkbl = max(i__2,i__3);
720 i__2 = wrkbl, i__3 = (*m << 1) + *m * ilaenv_(&c__1,
721 "ZUNGBR", "Q", m, m, m, &c_n1, (ftnlen)6, (ftnlen)
723 wrkbl = max(i__2,i__3);
724 maxwrk = *m * *m + wrkbl;
725 minwrk = (*m << 1) + *n;
726 maxwrk = max(minwrk,maxwrk);
727 } else if (wntva && wntun) {
729 /* Path 7t(N much larger than M, JOBU='N', JOBVT='A') */
731 wrkbl = *m + *m * ilaenv_(&c__1, "ZGELQF", " ", m, n, &
732 c_n1, &c_n1, (ftnlen)6, (ftnlen)1);
734 i__2 = wrkbl, i__3 = *m + *n * ilaenv_(&c__1, "ZUNGLQ",
735 " ", n, n, m, &c_n1, (ftnlen)6, (ftnlen)1);
736 wrkbl = max(i__2,i__3);
738 i__2 = wrkbl, i__3 = (*m << 1) + (*m << 1) * ilaenv_(&
739 c__1, "ZGEBRD", " ", m, m, &c_n1, &c_n1, (ftnlen)
741 wrkbl = max(i__2,i__3);
743 i__2 = wrkbl, i__3 = (*m << 1) + (*m - 1) * ilaenv_(&c__1,
744 "ZUNGBR", "P", m, m, m, &c_n1, (ftnlen)6, (
746 wrkbl = max(i__2,i__3);
747 maxwrk = *m * *m + wrkbl;
748 minwrk = (*m << 1) + *n;
749 maxwrk = max(minwrk,maxwrk);
750 } else if (wntva && wntuo) {
752 /* Path 8t(N much larger than M, JOBU='O', JOBVT='A') */
754 wrkbl = *m + *m * ilaenv_(&c__1, "ZGELQF", " ", m, n, &
755 c_n1, &c_n1, (ftnlen)6, (ftnlen)1);
757 i__2 = wrkbl, i__3 = *m + *n * ilaenv_(&c__1, "ZUNGLQ",
758 " ", n, n, m, &c_n1, (ftnlen)6, (ftnlen)1);
759 wrkbl = max(i__2,i__3);
761 i__2 = wrkbl, i__3 = (*m << 1) + (*m << 1) * ilaenv_(&
762 c__1, "ZGEBRD", " ", m, m, &c_n1, &c_n1, (ftnlen)
764 wrkbl = max(i__2,i__3);
766 i__2 = wrkbl, i__3 = (*m << 1) + (*m - 1) * ilaenv_(&c__1,
767 "ZUNGBR", "P", m, m, m, &c_n1, (ftnlen)6, (
769 wrkbl = max(i__2,i__3);
771 i__2 = wrkbl, i__3 = (*m << 1) + *m * ilaenv_(&c__1,
772 "ZUNGBR", "Q", m, m, m, &c_n1, (ftnlen)6, (ftnlen)
774 wrkbl = max(i__2,i__3);
775 maxwrk = (*m << 1) * *m + wrkbl;
776 minwrk = (*m << 1) + *n;
777 maxwrk = max(minwrk,maxwrk);
778 } else if (wntva && wntuas) {
780 /* Path 9t(N much larger than M, JOBU='S' or 'A',
783 wrkbl = *m + *m * ilaenv_(&c__1, "ZGELQF", " ", m, n, &
784 c_n1, &c_n1, (ftnlen)6, (ftnlen)1);
786 i__2 = wrkbl, i__3 = *m + *n * ilaenv_(&c__1, "ZUNGLQ",
787 " ", n, n, m, &c_n1, (ftnlen)6, (ftnlen)1);
788 wrkbl = max(i__2,i__3);
790 i__2 = wrkbl, i__3 = (*m << 1) + (*m << 1) * ilaenv_(&
791 c__1, "ZGEBRD", " ", m, m, &c_n1, &c_n1, (ftnlen)
793 wrkbl = max(i__2,i__3);
795 i__2 = wrkbl, i__3 = (*m << 1) + (*m - 1) * ilaenv_(&c__1,
796 "ZUNGBR", "P", m, m, m, &c_n1, (ftnlen)6, (
798 wrkbl = max(i__2,i__3);
800 i__2 = wrkbl, i__3 = (*m << 1) + *m * ilaenv_(&c__1,
801 "ZUNGBR", "Q", m, m, m, &c_n1, (ftnlen)6, (ftnlen)
803 wrkbl = max(i__2,i__3);
804 maxwrk = *m * *m + wrkbl;
805 minwrk = (*m << 1) + *n;
806 maxwrk = max(minwrk,maxwrk);
810 /* Path 10t(N greater than M, but not much larger) */
812 maxwrk = (*m << 1) + (*m + *n) * ilaenv_(&c__1, "ZGEBRD",
813 " ", m, n, &c_n1, &c_n1, (ftnlen)6, (ftnlen)1);
814 if (wntvs || wntvo) {
816 i__2 = maxwrk, i__3 = (*m << 1) + *m * ilaenv_(&c__1,
817 "ZUNGBR", "P", m, n, m, &c_n1, (ftnlen)6, (ftnlen)
819 maxwrk = max(i__2,i__3);
823 i__2 = maxwrk, i__3 = (*m << 1) + *n * ilaenv_(&c__1,
824 "ZUNGBR", "P", n, n, m, &c_n1, (ftnlen)6, (ftnlen)
826 maxwrk = max(i__2,i__3);
830 i__2 = maxwrk, i__3 = (*m << 1) + (*m - 1) * ilaenv_(&
831 c__1, "ZUNGBR", "Q", m, m, m, &c_n1, (ftnlen)6, (
833 maxwrk = max(i__2,i__3);
835 minwrk = (*m << 1) + *n;
836 maxwrk = max(minwrk,maxwrk);
839 work[1].r = (doublereal) maxwrk, work[1].i = 0.;
842 if (*lwork < minwrk && ! lquery) {
847 xerbla_("ZGESVD", &i__2);
853 /* Quick return if possible */
855 if (*m == 0 || *n == 0) {
857 work[1].r = 1., work[1].i = 0.;
862 /* Get machine constants */
865 smlnum = sqrt(dlamch_("S")) / eps;
866 bignum = 1. / smlnum;
868 /* Scale A if max element outside range [SMLNUM,BIGNUM] */
870 anrm = zlange_("M", m, n, &a[a_offset], lda, dum);
872 if (anrm > 0. && anrm < smlnum) {
874 zlascl_("G", &c__0, &c__0, &anrm, &smlnum, m, n, &a[a_offset], lda, &
876 } else if (anrm > bignum) {
878 zlascl_("G", &c__0, &c__0, &anrm, &bignum, m, n, &a[a_offset], lda, &
884 /* A has at least as many rows as columns. If A has sufficiently
885 more rows than columns, first reduce using the QR
886 decomposition (if sufficient workspace available) */
892 /* Path 1 (M much larger than N, JOBU='N')
893 No left singular vectors to be computed */
899 (CWorkspace: need 2*N, prefer N+N*NB)
900 (RWorkspace: need 0) */
902 i__2 = *lwork - iwork + 1;
903 zgeqrf_(m, n, &a[a_offset], lda, &work[itau], &work[iwork], &
906 /* Zero out below R */
910 zlaset_("L", &i__2, &i__3, &c_b1, &c_b1, &a_ref(2, 1), lda);
916 /* Bidiagonalize R in A
917 (CWorkspace: need 3*N, prefer 2*N+2*N*NB)
918 (RWorkspace: need N) */
920 i__2 = *lwork - iwork + 1;
921 zgebrd_(n, n, &a[a_offset], lda, &s[1], &rwork[ie], &work[
922 itauq], &work[itaup], &work[iwork], &i__2, &ierr);
924 if (wntvo || wntvas) {
926 /* If right singular vectors desired, generate P'.
927 (CWorkspace: need 3*N-1, prefer 2*N+(N-1)*NB)
930 i__2 = *lwork - iwork + 1;
931 zungbr_("P", n, n, n, &a[a_offset], lda, &work[itaup], &
932 work[iwork], &i__2, &ierr);
937 /* Perform bidiagonal QR iteration, computing right
938 singular vectors of A in A if desired
940 (RWorkspace: need BDSPAC) */
942 zbdsqr_("U", n, &ncvt, &c__0, &c__0, &s[1], &rwork[ie], &a[
943 a_offset], lda, cdum, &c__1, cdum, &c__1, &rwork[
946 /* If right singular vectors desired in VT, copy them there */
949 zlacpy_("F", n, n, &a[a_offset], lda, &vt[vt_offset],
953 } else if (wntuo && wntvn) {
955 /* Path 2 (M much larger than N, JOBU='O', JOBVT='N')
956 N left singular vectors to be overwritten on A and
957 no right singular vectors to be computed */
959 if (*lwork >= *n * *n + *n * 3) {
961 /* Sufficient workspace for a fast algorithm */
965 i__2 = wrkbl, i__3 = *lda * *n;
966 if (*lwork >= max(i__2,i__3) + *lda * *n) {
968 /* WORK(IU) is LDA by N, WORK(IR) is LDA by N */
972 } else /* if(complicated condition) */ {
974 i__2 = wrkbl, i__3 = *lda * *n;
975 if (*lwork >= max(i__2,i__3) + *n * *n) {
977 /* WORK(IU) is LDA by N, WORK(IR) is N by N */
983 /* WORK(IU) is LDWRKU by N, WORK(IR) is N by N */
985 ldwrku = (*lwork - *n * *n) / *n;
989 itau = ir + ldwrkr * *n;
993 (CWorkspace: need N*N+2*N, prefer N*N+N+N*NB)
996 i__2 = *lwork - iwork + 1;
997 zgeqrf_(m, n, &a[a_offset], lda, &work[itau], &work[iwork]
1000 /* Copy R to WORK(IR) and zero out below it */
1002 zlacpy_("U", n, n, &a[a_offset], lda, &work[ir], &ldwrkr);
1005 zlaset_("L", &i__2, &i__3, &c_b1, &c_b1, &work[ir + 1], &
1009 (CWorkspace: need N*N+2*N, prefer N*N+N+N*NB)
1012 i__2 = *lwork - iwork + 1;
1013 zungqr_(m, n, n, &a[a_offset], lda, &work[itau], &work[
1014 iwork], &i__2, &ierr);
1020 /* Bidiagonalize R in WORK(IR)
1021 (CWorkspace: need N*N+3*N, prefer N*N+2*N+2*N*NB)
1022 (RWorkspace: need N) */
1024 i__2 = *lwork - iwork + 1;
1025 zgebrd_(n, n, &work[ir], &ldwrkr, &s[1], &rwork[ie], &
1026 work[itauq], &work[itaup], &work[iwork], &i__2, &
1029 /* Generate left vectors bidiagonalizing R
1030 (CWorkspace: need N*N+3*N, prefer N*N+2*N+N*NB)
1031 (RWorkspace: need 0) */
1033 i__2 = *lwork - iwork + 1;
1034 zungbr_("Q", n, n, n, &work[ir], &ldwrkr, &work[itauq], &
1035 work[iwork], &i__2, &ierr);
1038 /* Perform bidiagonal QR iteration, computing left
1039 singular vectors of R in WORK(IR)
1040 (CWorkspace: need N*N)
1041 (RWorkspace: need BDSPAC) */
1043 zbdsqr_("U", n, &c__0, n, &c__0, &s[1], &rwork[ie], cdum,
1044 &c__1, &work[ir], &ldwrkr, cdum, &c__1, &rwork[
1048 /* Multiply Q in A by left singular vectors of R in
1049 WORK(IR), storing result in WORK(IU) and copying to A
1050 (CWorkspace: need N*N+N, prefer N*N+M*N)
1055 for (i__ = 1; i__3 < 0 ? i__ >= i__2 : i__ <= i__2; i__ +=
1058 i__4 = *m - i__ + 1;
1059 chunk = min(i__4,ldwrku);
1060 zgemm_("N", "N", &chunk, n, n, &c_b2, &a_ref(i__, 1),
1061 lda, &work[ir], &ldwrkr, &c_b1, &work[iu], &
1063 zlacpy_("F", &chunk, n, &work[iu], &ldwrku, &a_ref(
1070 /* Insufficient workspace for a fast algorithm */
1078 (CWorkspace: need 2*N+M, prefer 2*N+(M+N)*NB)
1081 i__3 = *lwork - iwork + 1;
1082 zgebrd_(m, n, &a[a_offset], lda, &s[1], &rwork[ie], &work[
1083 itauq], &work[itaup], &work[iwork], &i__3, &ierr);
1085 /* Generate left vectors bidiagonalizing A
1086 (CWorkspace: need 3*N, prefer 2*N+N*NB)
1089 i__3 = *lwork - iwork + 1;
1090 zungbr_("Q", m, n, n, &a[a_offset], lda, &work[itauq], &
1091 work[iwork], &i__3, &ierr);
1094 /* Perform bidiagonal QR iteration, computing left
1095 singular vectors of A in A
1096 (CWorkspace: need 0)
1097 (RWorkspace: need BDSPAC) */
1099 zbdsqr_("U", n, &c__0, m, &c__0, &s[1], &rwork[ie], cdum,
1100 &c__1, &a[a_offset], lda, cdum, &c__1, &rwork[
1105 } else if (wntuo && wntvas) {
1107 /* Path 3 (M much larger than N, JOBU='O', JOBVT='S' or 'A')
1108 N left singular vectors to be overwritten on A and
1109 N right singular vectors to be computed in VT */
1111 if (*lwork >= *n * *n + *n * 3) {
1113 /* Sufficient workspace for a fast algorithm */
1117 i__3 = wrkbl, i__2 = *lda * *n;
1118 if (*lwork >= max(i__3,i__2) + *lda * *n) {
1120 /* WORK(IU) is LDA by N and WORK(IR) is LDA by N */
1124 } else /* if(complicated condition) */ {
1126 i__3 = wrkbl, i__2 = *lda * *n;
1127 if (*lwork >= max(i__3,i__2) + *n * *n) {
1129 /* WORK(IU) is LDA by N and WORK(IR) is N by N */
1135 /* WORK(IU) is LDWRKU by N and WORK(IR) is N by N */
1137 ldwrku = (*lwork - *n * *n) / *n;
1141 itau = ir + ldwrkr * *n;
1145 (CWorkspace: need N*N+2*N, prefer N*N+N+N*NB)
1148 i__3 = *lwork - iwork + 1;
1149 zgeqrf_(m, n, &a[a_offset], lda, &work[itau], &work[iwork]
1152 /* Copy R to VT, zeroing out below it */
1154 zlacpy_("U", n, n, &a[a_offset], lda, &vt[vt_offset],
1158 zlaset_("L", &i__3, &i__2, &c_b1, &c_b1, &vt_ref(2, 1),
1162 (CWorkspace: need N*N+2*N, prefer N*N+N+N*NB)
1165 i__3 = *lwork - iwork + 1;
1166 zungqr_(m, n, n, &a[a_offset], lda, &work[itau], &work[
1167 iwork], &i__3, &ierr);
1173 /* Bidiagonalize R in VT, copying result to WORK(IR)
1174 (CWorkspace: need N*N+3*N, prefer N*N+2*N+2*N*NB)
1175 (RWorkspace: need N) */
1177 i__3 = *lwork - iwork + 1;
1178 zgebrd_(n, n, &vt[vt_offset], ldvt, &s[1], &rwork[ie], &
1179 work[itauq], &work[itaup], &work[iwork], &i__3, &
1181 zlacpy_("L", n, n, &vt[vt_offset], ldvt, &work[ir], &
1184 /* Generate left vectors bidiagonalizing R in WORK(IR)
1185 (CWorkspace: need N*N+3*N, prefer N*N+2*N+N*NB)
1188 i__3 = *lwork - iwork + 1;
1189 zungbr_("Q", n, n, n, &work[ir], &ldwrkr, &work[itauq], &
1190 work[iwork], &i__3, &ierr);
1192 /* Generate right vectors bidiagonalizing R in VT
1193 (CWorkspace: need N*N+3*N-1, prefer N*N+2*N+(N-1)*NB)
1196 i__3 = *lwork - iwork + 1;
1197 zungbr_("P", n, n, n, &vt[vt_offset], ldvt, &work[itaup],
1198 &work[iwork], &i__3, &ierr);
1201 /* Perform bidiagonal QR iteration, computing left
1202 singular vectors of R in WORK(IR) and computing right
1203 singular vectors of R in VT
1204 (CWorkspace: need N*N)
1205 (RWorkspace: need BDSPAC) */
1207 zbdsqr_("U", n, n, n, &c__0, &s[1], &rwork[ie], &vt[
1208 vt_offset], ldvt, &work[ir], &ldwrkr, cdum, &c__1,
1209 &rwork[irwork], info);
1212 /* Multiply Q in A by left singular vectors of R in
1213 WORK(IR), storing result in WORK(IU) and copying to A
1214 (CWorkspace: need N*N+N, prefer N*N+M*N)
1219 for (i__ = 1; i__2 < 0 ? i__ >= i__3 : i__ <= i__3; i__ +=
1222 i__4 = *m - i__ + 1;
1223 chunk = min(i__4,ldwrku);
1224 zgemm_("N", "N", &chunk, n, n, &c_b2, &a_ref(i__, 1),
1225 lda, &work[ir], &ldwrkr, &c_b1, &work[iu], &
1227 zlacpy_("F", &chunk, n, &work[iu], &ldwrku, &a_ref(
1234 /* Insufficient workspace for a fast algorithm */
1240 (CWorkspace: need 2*N, prefer N+N*NB)
1243 i__2 = *lwork - iwork + 1;
1244 zgeqrf_(m, n, &a[a_offset], lda, &work[itau], &work[iwork]
1247 /* Copy R to VT, zeroing out below it */
1249 zlacpy_("U", n, n, &a[a_offset], lda, &vt[vt_offset],
1253 zlaset_("L", &i__2, &i__3, &c_b1, &c_b1, &vt_ref(2, 1),
1257 (CWorkspace: need 2*N, prefer N+N*NB)
1260 i__2 = *lwork - iwork + 1;
1261 zungqr_(m, n, n, &a[a_offset], lda, &work[itau], &work[
1262 iwork], &i__2, &ierr);
1268 /* Bidiagonalize R in VT
1269 (CWorkspace: need 3*N, prefer 2*N+2*N*NB)
1272 i__2 = *lwork - iwork + 1;
1273 zgebrd_(n, n, &vt[vt_offset], ldvt, &s[1], &rwork[ie], &
1274 work[itauq], &work[itaup], &work[iwork], &i__2, &
1277 /* Multiply Q in A by left vectors bidiagonalizing R
1278 (CWorkspace: need 2*N+M, prefer 2*N+M*NB)
1281 i__2 = *lwork - iwork + 1;
1282 zunmbr_("Q", "R", "N", m, n, n, &vt[vt_offset], ldvt, &
1283 work[itauq], &a[a_offset], lda, &work[iwork], &
1286 /* Generate right vectors bidiagonalizing R in VT
1287 (CWorkspace: need 3*N-1, prefer 2*N+(N-1)*NB)
1290 i__2 = *lwork - iwork + 1;
1291 zungbr_("P", n, n, n, &vt[vt_offset], ldvt, &work[itaup],
1292 &work[iwork], &i__2, &ierr);
1295 /* Perform bidiagonal QR iteration, computing left
1296 singular vectors of A in A and computing right
1297 singular vectors of A in VT
1299 (RWorkspace: need BDSPAC) */
1301 zbdsqr_("U", n, n, m, &c__0, &s[1], &rwork[ie], &vt[
1302 vt_offset], ldvt, &a[a_offset], lda, cdum, &c__1,
1303 &rwork[irwork], info);
1311 /* Path 4 (M much larger than N, JOBU='S', JOBVT='N')
1312 N left singular vectors to be computed in U and
1313 no right singular vectors to be computed */
1315 if (*lwork >= *n * *n + *n * 3) {
1317 /* Sufficient workspace for a fast algorithm */
1320 if (*lwork >= wrkbl + *lda * *n) {
1322 /* WORK(IR) is LDA by N */
1327 /* WORK(IR) is N by N */
1331 itau = ir + ldwrkr * *n;
1335 (CWorkspace: need N*N+2*N, prefer N*N+N+N*NB)
1338 i__2 = *lwork - iwork + 1;
1339 zgeqrf_(m, n, &a[a_offset], lda, &work[itau], &work[
1340 iwork], &i__2, &ierr);
1342 /* Copy R to WORK(IR), zeroing out below it */
1344 zlacpy_("U", n, n, &a[a_offset], lda, &work[ir], &
1348 zlaset_("L", &i__2, &i__3, &c_b1, &c_b1, &work[ir + 1]
1352 (CWorkspace: need N*N+2*N, prefer N*N+N+N*NB)
1355 i__2 = *lwork - iwork + 1;
1356 zungqr_(m, n, n, &a[a_offset], lda, &work[itau], &
1357 work[iwork], &i__2, &ierr);
1363 /* Bidiagonalize R in WORK(IR)
1364 (CWorkspace: need N*N+3*N, prefer N*N+2*N+2*N*NB)
1365 (RWorkspace: need N) */
1367 i__2 = *lwork - iwork + 1;
1368 zgebrd_(n, n, &work[ir], &ldwrkr, &s[1], &rwork[ie], &
1369 work[itauq], &work[itaup], &work[iwork], &
1372 /* Generate left vectors bidiagonalizing R in WORK(IR)
1373 (CWorkspace: need N*N+3*N, prefer N*N+2*N+N*NB)
1376 i__2 = *lwork - iwork + 1;
1377 zungbr_("Q", n, n, n, &work[ir], &ldwrkr, &work[itauq]
1378 , &work[iwork], &i__2, &ierr);
1381 /* Perform bidiagonal QR iteration, computing left
1382 singular vectors of R in WORK(IR)
1383 (CWorkspace: need N*N)
1384 (RWorkspace: need BDSPAC) */
1386 zbdsqr_("U", n, &c__0, n, &c__0, &s[1], &rwork[ie],
1387 cdum, &c__1, &work[ir], &ldwrkr, cdum, &c__1,
1388 &rwork[irwork], info);
1390 /* Multiply Q in A by left singular vectors of R in
1391 WORK(IR), storing result in U
1392 (CWorkspace: need N*N)
1395 zgemm_("N", "N", m, n, n, &c_b2, &a[a_offset], lda, &
1396 work[ir], &ldwrkr, &c_b1, &u[u_offset], ldu);
1400 /* Insufficient workspace for a fast algorithm */
1405 /* Compute A=Q*R, copying result to U
1406 (CWorkspace: need 2*N, prefer N+N*NB)
1409 i__2 = *lwork - iwork + 1;
1410 zgeqrf_(m, n, &a[a_offset], lda, &work[itau], &work[
1411 iwork], &i__2, &ierr);
1412 zlacpy_("L", m, n, &a[a_offset], lda, &u[u_offset],
1416 (CWorkspace: need 2*N, prefer N+N*NB)
1419 i__2 = *lwork - iwork + 1;
1420 zungqr_(m, n, n, &u[u_offset], ldu, &work[itau], &
1421 work[iwork], &i__2, &ierr);
1427 /* Zero out below R in A */
1431 zlaset_("L", &i__2, &i__3, &c_b1, &c_b1, &a_ref(2, 1),
1434 /* Bidiagonalize R in A
1435 (CWorkspace: need 3*N, prefer 2*N+2*N*NB)
1436 (RWorkspace: need N) */
1438 i__2 = *lwork - iwork + 1;
1439 zgebrd_(n, n, &a[a_offset], lda, &s[1], &rwork[ie], &
1440 work[itauq], &work[itaup], &work[iwork], &
1443 /* Multiply Q in U by left vectors bidiagonalizing R
1444 (CWorkspace: need 2*N+M, prefer 2*N+M*NB)
1447 i__2 = *lwork - iwork + 1;
1448 zunmbr_("Q", "R", "N", m, n, n, &a[a_offset], lda, &
1449 work[itauq], &u[u_offset], ldu, &work[iwork],
1454 /* Perform bidiagonal QR iteration, computing left
1455 singular vectors of A in U
1457 (RWorkspace: need BDSPAC) */
1459 zbdsqr_("U", n, &c__0, m, &c__0, &s[1], &rwork[ie],
1460 cdum, &c__1, &u[u_offset], ldu, cdum, &c__1, &
1461 rwork[irwork], info);
1467 /* Path 5 (M much larger than N, JOBU='S', JOBVT='O')
1468 N left singular vectors to be computed in U and
1469 N right singular vectors to be overwritten on A */
1471 if (*lwork >= (*n << 1) * *n + *n * 3) {
1473 /* Sufficient workspace for a fast algorithm */
1476 if (*lwork >= wrkbl + (*lda << 1) * *n) {
1478 /* WORK(IU) is LDA by N and WORK(IR) is LDA by N */
1481 ir = iu + ldwrku * *n;
1483 } else if (*lwork >= wrkbl + (*lda + *n) * *n) {
1485 /* WORK(IU) is LDA by N and WORK(IR) is N by N */
1488 ir = iu + ldwrku * *n;
1492 /* WORK(IU) is N by N and WORK(IR) is N by N */
1495 ir = iu + ldwrku * *n;
1498 itau = ir + ldwrkr * *n;
1502 (CWorkspace: need 2*N*N+2*N, prefer 2*N*N+N+N*NB)
1505 i__2 = *lwork - iwork + 1;
1506 zgeqrf_(m, n, &a[a_offset], lda, &work[itau], &work[
1507 iwork], &i__2, &ierr);
1509 /* Copy R to WORK(IU), zeroing out below it */
1511 zlacpy_("U", n, n, &a[a_offset], lda, &work[iu], &
1515 zlaset_("L", &i__2, &i__3, &c_b1, &c_b1, &work[iu + 1]
1519 (CWorkspace: need 2*N*N+2*N, prefer 2*N*N+N+N*NB)
1522 i__2 = *lwork - iwork + 1;
1523 zungqr_(m, n, n, &a[a_offset], lda, &work[itau], &
1524 work[iwork], &i__2, &ierr);
1530 /* Bidiagonalize R in WORK(IU), copying result to
1532 (CWorkspace: need 2*N*N+3*N,
1533 prefer 2*N*N+2*N+2*N*NB)
1534 (RWorkspace: need N) */
1536 i__2 = *lwork - iwork + 1;
1537 zgebrd_(n, n, &work[iu], &ldwrku, &s[1], &rwork[ie], &
1538 work[itauq], &work[itaup], &work[iwork], &
1540 zlacpy_("U", n, n, &work[iu], &ldwrku, &work[ir], &
1543 /* Generate left bidiagonalizing vectors in WORK(IU)
1544 (CWorkspace: need 2*N*N+3*N, prefer 2*N*N+2*N+N*NB)
1547 i__2 = *lwork - iwork + 1;
1548 zungbr_("Q", n, n, n, &work[iu], &ldwrku, &work[itauq]
1549 , &work[iwork], &i__2, &ierr);
1551 /* Generate right bidiagonalizing vectors in WORK(IR)
1552 (CWorkspace: need 2*N*N+3*N-1,
1553 prefer 2*N*N+2*N+(N-1)*NB)
1556 i__2 = *lwork - iwork + 1;
1557 zungbr_("P", n, n, n, &work[ir], &ldwrkr, &work[itaup]
1558 , &work[iwork], &i__2, &ierr);
1561 /* Perform bidiagonal QR iteration, computing left
1562 singular vectors of R in WORK(IU) and computing
1563 right singular vectors of R in WORK(IR)
1564 (CWorkspace: need 2*N*N)
1565 (RWorkspace: need BDSPAC) */
1567 zbdsqr_("U", n, n, n, &c__0, &s[1], &rwork[ie], &work[
1568 ir], &ldwrkr, &work[iu], &ldwrku, cdum, &c__1,
1569 &rwork[irwork], info);
1571 /* Multiply Q in A by left singular vectors of R in
1572 WORK(IU), storing result in U
1573 (CWorkspace: need N*N)
1576 zgemm_("N", "N", m, n, n, &c_b2, &a[a_offset], lda, &
1577 work[iu], &ldwrku, &c_b1, &u[u_offset], ldu);
1579 /* Copy right singular vectors of R to A
1580 (CWorkspace: need N*N)
1583 zlacpy_("F", n, n, &work[ir], &ldwrkr, &a[a_offset],
1588 /* Insufficient workspace for a fast algorithm */
1593 /* Compute A=Q*R, copying result to U
1594 (CWorkspace: need 2*N, prefer N+N*NB)
1597 i__2 = *lwork - iwork + 1;
1598 zgeqrf_(m, n, &a[a_offset], lda, &work[itau], &work[
1599 iwork], &i__2, &ierr);
1600 zlacpy_("L", m, n, &a[a_offset], lda, &u[u_offset],
1604 (CWorkspace: need 2*N, prefer N+N*NB)
1607 i__2 = *lwork - iwork + 1;
1608 zungqr_(m, n, n, &u[u_offset], ldu, &work[itau], &
1609 work[iwork], &i__2, &ierr);
1615 /* Zero out below R in A */
1619 zlaset_("L", &i__2, &i__3, &c_b1, &c_b1, &a_ref(2, 1),
1622 /* Bidiagonalize R in A
1623 (CWorkspace: need 3*N, prefer 2*N+2*N*NB)
1624 (RWorkspace: need N) */
1626 i__2 = *lwork - iwork + 1;
1627 zgebrd_(n, n, &a[a_offset], lda, &s[1], &rwork[ie], &
1628 work[itauq], &work[itaup], &work[iwork], &
1631 /* Multiply Q in U by left vectors bidiagonalizing R
1632 (CWorkspace: need 2*N+M, prefer 2*N+M*NB)
1635 i__2 = *lwork - iwork + 1;
1636 zunmbr_("Q", "R", "N", m, n, n, &a[a_offset], lda, &
1637 work[itauq], &u[u_offset], ldu, &work[iwork],
1641 /* Generate right vectors bidiagonalizing R in A
1642 (CWorkspace: need 3*N-1, prefer 2*N+(N-1)*NB)
1645 i__2 = *lwork - iwork + 1;
1646 zungbr_("P", n, n, n, &a[a_offset], lda, &work[itaup],
1647 &work[iwork], &i__2, &ierr);
1650 /* Perform bidiagonal QR iteration, computing left
1651 singular vectors of A in U and computing right
1652 singular vectors of A in A
1654 (RWorkspace: need BDSPAC) */
1656 zbdsqr_("U", n, n, m, &c__0, &s[1], &rwork[ie], &a[
1657 a_offset], lda, &u[u_offset], ldu, cdum, &
1658 c__1, &rwork[irwork], info);
1662 } else if (wntvas) {
1664 /* Path 6 (M much larger than N, JOBU='S', JOBVT='S'
1666 N left singular vectors to be computed in U and
1667 N right singular vectors to be computed in VT */
1669 if (*lwork >= *n * *n + *n * 3) {
1671 /* Sufficient workspace for a fast algorithm */
1674 if (*lwork >= wrkbl + *lda * *n) {
1676 /* WORK(IU) is LDA by N */
1681 /* WORK(IU) is N by N */
1685 itau = iu + ldwrku * *n;
1689 (CWorkspace: need N*N+2*N, prefer N*N+N+N*NB)
1692 i__2 = *lwork - iwork + 1;
1693 zgeqrf_(m, n, &a[a_offset], lda, &work[itau], &work[
1694 iwork], &i__2, &ierr);
1696 /* Copy R to WORK(IU), zeroing out below it */
1698 zlacpy_("U", n, n, &a[a_offset], lda, &work[iu], &
1702 zlaset_("L", &i__2, &i__3, &c_b1, &c_b1, &work[iu + 1]
1706 (CWorkspace: need N*N+2*N, prefer N*N+N+N*NB)
1709 i__2 = *lwork - iwork + 1;
1710 zungqr_(m, n, n, &a[a_offset], lda, &work[itau], &
1711 work[iwork], &i__2, &ierr);
1717 /* Bidiagonalize R in WORK(IU), copying result to VT
1718 (CWorkspace: need N*N+3*N, prefer N*N+2*N+2*N*NB)
1719 (RWorkspace: need N) */
1721 i__2 = *lwork - iwork + 1;
1722 zgebrd_(n, n, &work[iu], &ldwrku, &s[1], &rwork[ie], &
1723 work[itauq], &work[itaup], &work[iwork], &
1725 zlacpy_("U", n, n, &work[iu], &ldwrku, &vt[vt_offset],
1728 /* Generate left bidiagonalizing vectors in WORK(IU)
1729 (CWorkspace: need N*N+3*N, prefer N*N+2*N+N*NB)
1732 i__2 = *lwork - iwork + 1;
1733 zungbr_("Q", n, n, n, &work[iu], &ldwrku, &work[itauq]
1734 , &work[iwork], &i__2, &ierr);
1736 /* Generate right bidiagonalizing vectors in VT
1737 (CWorkspace: need N*N+3*N-1,
1738 prefer N*N+2*N+(N-1)*NB)
1741 i__2 = *lwork - iwork + 1;
1742 zungbr_("P", n, n, n, &vt[vt_offset], ldvt, &work[
1743 itaup], &work[iwork], &i__2, &ierr)
1747 /* Perform bidiagonal QR iteration, computing left
1748 singular vectors of R in WORK(IU) and computing
1749 right singular vectors of R in VT
1750 (CWorkspace: need N*N)
1751 (RWorkspace: need BDSPAC) */
1753 zbdsqr_("U", n, n, n, &c__0, &s[1], &rwork[ie], &vt[
1754 vt_offset], ldvt, &work[iu], &ldwrku, cdum, &
1755 c__1, &rwork[irwork], info);
1757 /* Multiply Q in A by left singular vectors of R in
1758 WORK(IU), storing result in U
1759 (CWorkspace: need N*N)
1762 zgemm_("N", "N", m, n, n, &c_b2, &a[a_offset], lda, &
1763 work[iu], &ldwrku, &c_b1, &u[u_offset], ldu);
1767 /* Insufficient workspace for a fast algorithm */
1772 /* Compute A=Q*R, copying result to U
1773 (CWorkspace: need 2*N, prefer N+N*NB)
1776 i__2 = *lwork - iwork + 1;
1777 zgeqrf_(m, n, &a[a_offset], lda, &work[itau], &work[
1778 iwork], &i__2, &ierr);
1779 zlacpy_("L", m, n, &a[a_offset], lda, &u[u_offset],
1783 (CWorkspace: need 2*N, prefer N+N*NB)
1786 i__2 = *lwork - iwork + 1;
1787 zungqr_(m, n, n, &u[u_offset], ldu, &work[itau], &
1788 work[iwork], &i__2, &ierr);
1790 /* Copy R to VT, zeroing out below it */
1792 zlacpy_("U", n, n, &a[a_offset], lda, &vt[vt_offset],
1796 zlaset_("L", &i__2, &i__3, &c_b1, &c_b1, &vt_ref(2, 1)
1803 /* Bidiagonalize R in VT
1804 (CWorkspace: need 3*N, prefer 2*N+2*N*NB)
1805 (RWorkspace: need N) */
1807 i__2 = *lwork - iwork + 1;
1808 zgebrd_(n, n, &vt[vt_offset], ldvt, &s[1], &rwork[ie],
1809 &work[itauq], &work[itaup], &work[iwork], &
1812 /* Multiply Q in U by left bidiagonalizing vectors
1814 (CWorkspace: need 2*N+M, prefer 2*N+M*NB)
1817 i__2 = *lwork - iwork + 1;
1818 zunmbr_("Q", "R", "N", m, n, n, &vt[vt_offset], ldvt,
1819 &work[itauq], &u[u_offset], ldu, &work[iwork],
1822 /* Generate right bidiagonalizing vectors in VT
1823 (CWorkspace: need 3*N-1, prefer 2*N+(N-1)*NB)
1826 i__2 = *lwork - iwork + 1;
1827 zungbr_("P", n, n, n, &vt[vt_offset], ldvt, &work[
1828 itaup], &work[iwork], &i__2, &ierr)
1832 /* Perform bidiagonal QR iteration, computing left
1833 singular vectors of A in U and computing right
1834 singular vectors of A in VT
1836 (RWorkspace: need BDSPAC) */
1838 zbdsqr_("U", n, n, m, &c__0, &s[1], &rwork[ie], &vt[
1839 vt_offset], ldvt, &u[u_offset], ldu, cdum, &
1840 c__1, &rwork[irwork], info);
1850 /* Path 7 (M much larger than N, JOBU='A', JOBVT='N')
1851 M left singular vectors to be computed in U and
1852 no right singular vectors to be computed
1855 i__2 = *n + *m, i__3 = *n * 3;
1856 if (*lwork >= *n * *n + max(i__2,i__3)) {
1858 /* Sufficient workspace for a fast algorithm */
1861 if (*lwork >= wrkbl + *lda * *n) {
1863 /* WORK(IR) is LDA by N */
1868 /* WORK(IR) is N by N */
1872 itau = ir + ldwrkr * *n;
1875 /* Compute A=Q*R, copying result to U
1876 (CWorkspace: need N*N+2*N, prefer N*N+N+N*NB)
1879 i__2 = *lwork - iwork + 1;
1880 zgeqrf_(m, n, &a[a_offset], lda, &work[itau], &work[
1881 iwork], &i__2, &ierr);
1882 zlacpy_("L", m, n, &a[a_offset], lda, &u[u_offset],
1885 /* Copy R to WORK(IR), zeroing out below it */
1887 zlacpy_("U", n, n, &a[a_offset], lda, &work[ir], &
1891 zlaset_("L", &i__2, &i__3, &c_b1, &c_b1, &work[ir + 1]
1895 (CWorkspace: need N*N+N+M, prefer N*N+N+M*NB)
1898 i__2 = *lwork - iwork + 1;
1899 zungqr_(m, m, n, &u[u_offset], ldu, &work[itau], &
1900 work[iwork], &i__2, &ierr);
1906 /* Bidiagonalize R in WORK(IR)
1907 (CWorkspace: need N*N+3*N, prefer N*N+2*N+2*N*NB)
1908 (RWorkspace: need N) */
1910 i__2 = *lwork - iwork + 1;
1911 zgebrd_(n, n, &work[ir], &ldwrkr, &s[1], &rwork[ie], &
1912 work[itauq], &work[itaup], &work[iwork], &
1915 /* Generate left bidiagonalizing vectors in WORK(IR)
1916 (CWorkspace: need N*N+3*N, prefer N*N+2*N+N*NB)
1919 i__2 = *lwork - iwork + 1;
1920 zungbr_("Q", n, n, n, &work[ir], &ldwrkr, &work[itauq]
1921 , &work[iwork], &i__2, &ierr);
1924 /* Perform bidiagonal QR iteration, computing left
1925 singular vectors of R in WORK(IR)
1926 (CWorkspace: need N*N)
1927 (RWorkspace: need BDSPAC) */
1929 zbdsqr_("U", n, &c__0, n, &c__0, &s[1], &rwork[ie],
1930 cdum, &c__1, &work[ir], &ldwrkr, cdum, &c__1,
1931 &rwork[irwork], info);
1933 /* Multiply Q in U by left singular vectors of R in
1934 WORK(IR), storing result in A
1935 (CWorkspace: need N*N)
1938 zgemm_("N", "N", m, n, n, &c_b2, &u[u_offset], ldu, &
1939 work[ir], &ldwrkr, &c_b1, &a[a_offset], lda);
1941 /* Copy left singular vectors of A from A to U */
1943 zlacpy_("F", m, n, &a[a_offset], lda, &u[u_offset],
1948 /* Insufficient workspace for a fast algorithm */
1953 /* Compute A=Q*R, copying result to U
1954 (CWorkspace: need 2*N, prefer N+N*NB)
1957 i__2 = *lwork - iwork + 1;
1958 zgeqrf_(m, n, &a[a_offset], lda, &work[itau], &work[
1959 iwork], &i__2, &ierr);
1960 zlacpy_("L", m, n, &a[a_offset], lda, &u[u_offset],
1964 (CWorkspace: need N+M, prefer N+M*NB)
1967 i__2 = *lwork - iwork + 1;
1968 zungqr_(m, m, n, &u[u_offset], ldu, &work[itau], &
1969 work[iwork], &i__2, &ierr);
1975 /* Zero out below R in A */
1979 zlaset_("L", &i__2, &i__3, &c_b1, &c_b1, &a_ref(2, 1),
1982 /* Bidiagonalize R in A
1983 (CWorkspace: need 3*N, prefer 2*N+2*N*NB)
1984 (RWorkspace: need N) */
1986 i__2 = *lwork - iwork + 1;
1987 zgebrd_(n, n, &a[a_offset], lda, &s[1], &rwork[ie], &
1988 work[itauq], &work[itaup], &work[iwork], &
1991 /* Multiply Q in U by left bidiagonalizing vectors
1993 (CWorkspace: need 2*N+M, prefer 2*N+M*NB)
1996 i__2 = *lwork - iwork + 1;
1997 zunmbr_("Q", "R", "N", m, n, n, &a[a_offset], lda, &
1998 work[itauq], &u[u_offset], ldu, &work[iwork],
2003 /* Perform bidiagonal QR iteration, computing left
2004 singular vectors of A in U
2006 (RWorkspace: need BDSPAC) */
2008 zbdsqr_("U", n, &c__0, m, &c__0, &s[1], &rwork[ie],
2009 cdum, &c__1, &u[u_offset], ldu, cdum, &c__1, &
2010 rwork[irwork], info);
2016 /* Path 8 (M much larger than N, JOBU='A', JOBVT='O')
2017 M left singular vectors to be computed in U and
2018 N right singular vectors to be overwritten on A
2021 i__2 = *n + *m, i__3 = *n * 3;
2022 if (*lwork >= (*n << 1) * *n + max(i__2,i__3)) {
2024 /* Sufficient workspace for a fast algorithm */
2027 if (*lwork >= wrkbl + (*lda << 1) * *n) {
2029 /* WORK(IU) is LDA by N and WORK(IR) is LDA by N */
2032 ir = iu + ldwrku * *n;
2034 } else if (*lwork >= wrkbl + (*lda + *n) * *n) {
2036 /* WORK(IU) is LDA by N and WORK(IR) is N by N */
2039 ir = iu + ldwrku * *n;
2043 /* WORK(IU) is N by N and WORK(IR) is N by N */
2046 ir = iu + ldwrku * *n;
2049 itau = ir + ldwrkr * *n;
2052 /* Compute A=Q*R, copying result to U
2053 (CWorkspace: need 2*N*N+2*N, prefer 2*N*N+N+N*NB)
2056 i__2 = *lwork - iwork + 1;
2057 zgeqrf_(m, n, &a[a_offset], lda, &work[itau], &work[
2058 iwork], &i__2, &ierr);
2059 zlacpy_("L", m, n, &a[a_offset], lda, &u[u_offset],
2063 (CWorkspace: need 2*N*N+N+M, prefer 2*N*N+N+M*NB)
2066 i__2 = *lwork - iwork + 1;
2067 zungqr_(m, m, n, &u[u_offset], ldu, &work[itau], &
2068 work[iwork], &i__2, &ierr);
2070 /* Copy R to WORK(IU), zeroing out below it */
2072 zlacpy_("U", n, n, &a[a_offset], lda, &work[iu], &
2076 zlaset_("L", &i__2, &i__3, &c_b1, &c_b1, &work[iu + 1]
2083 /* Bidiagonalize R in WORK(IU), copying result to
2085 (CWorkspace: need 2*N*N+3*N,
2086 prefer 2*N*N+2*N+2*N*NB)
2087 (RWorkspace: need N) */
2089 i__2 = *lwork - iwork + 1;
2090 zgebrd_(n, n, &work[iu], &ldwrku, &s[1], &rwork[ie], &
2091 work[itauq], &work[itaup], &work[iwork], &
2093 zlacpy_("U", n, n, &work[iu], &ldwrku, &work[ir], &
2096 /* Generate left bidiagonalizing vectors in WORK(IU)
2097 (CWorkspace: need 2*N*N+3*N, prefer 2*N*N+2*N+N*NB)
2100 i__2 = *lwork - iwork + 1;
2101 zungbr_("Q", n, n, n, &work[iu], &ldwrku, &work[itauq]
2102 , &work[iwork], &i__2, &ierr);
2104 /* Generate right bidiagonalizing vectors in WORK(IR)
2105 (CWorkspace: need 2*N*N+3*N-1,
2106 prefer 2*N*N+2*N+(N-1)*NB)
2109 i__2 = *lwork - iwork + 1;
2110 zungbr_("P", n, n, n, &work[ir], &ldwrkr, &work[itaup]
2111 , &work[iwork], &i__2, &ierr);
2114 /* Perform bidiagonal QR iteration, computing left
2115 singular vectors of R in WORK(IU) and computing
2116 right singular vectors of R in WORK(IR)
2117 (CWorkspace: need 2*N*N)
2118 (RWorkspace: need BDSPAC) */
2120 zbdsqr_("U", n, n, n, &c__0, &s[1], &rwork[ie], &work[
2121 ir], &ldwrkr, &work[iu], &ldwrku, cdum, &c__1,
2122 &rwork[irwork], info);
2124 /* Multiply Q in U by left singular vectors of R in
2125 WORK(IU), storing result in A
2126 (CWorkspace: need N*N)
2129 zgemm_("N", "N", m, n, n, &c_b2, &u[u_offset], ldu, &
2130 work[iu], &ldwrku, &c_b1, &a[a_offset], lda);
2132 /* Copy left singular vectors of A from A to U */
2134 zlacpy_("F", m, n, &a[a_offset], lda, &u[u_offset],
2137 /* Copy right singular vectors of R from WORK(IR) to A */
2139 zlacpy_("F", n, n, &work[ir], &ldwrkr, &a[a_offset],
2144 /* Insufficient workspace for a fast algorithm */
2149 /* Compute A=Q*R, copying result to U
2150 (CWorkspace: need 2*N, prefer N+N*NB)
2153 i__2 = *lwork - iwork + 1;
2154 zgeqrf_(m, n, &a[a_offset], lda, &work[itau], &work[
2155 iwork], &i__2, &ierr);
2156 zlacpy_("L", m, n, &a[a_offset], lda, &u[u_offset],
2160 (CWorkspace: need N+M, prefer N+M*NB)
2163 i__2 = *lwork - iwork + 1;
2164 zungqr_(m, m, n, &u[u_offset], ldu, &work[itau], &
2165 work[iwork], &i__2, &ierr);
2171 /* Zero out below R in A */
2175 zlaset_("L", &i__2, &i__3, &c_b1, &c_b1, &a_ref(2, 1),
2178 /* Bidiagonalize R in A
2179 (CWorkspace: need 3*N, prefer 2*N+2*N*NB)
2180 (RWorkspace: need N) */
2182 i__2 = *lwork - iwork + 1;
2183 zgebrd_(n, n, &a[a_offset], lda, &s[1], &rwork[ie], &
2184 work[itauq], &work[itaup], &work[iwork], &
2187 /* Multiply Q in U by left bidiagonalizing vectors
2189 (CWorkspace: need 2*N+M, prefer 2*N+M*NB)
2192 i__2 = *lwork - iwork + 1;
2193 zunmbr_("Q", "R", "N", m, n, n, &a[a_offset], lda, &
2194 work[itauq], &u[u_offset], ldu, &work[iwork],
2198 /* Generate right bidiagonalizing vectors in A
2199 (CWorkspace: need 3*N-1, prefer 2*N+(N-1)*NB)
2202 i__2 = *lwork - iwork + 1;
2203 zungbr_("P", n, n, n, &a[a_offset], lda, &work[itaup],
2204 &work[iwork], &i__2, &ierr);
2207 /* Perform bidiagonal QR iteration, computing left
2208 singular vectors of A in U and computing right
2209 singular vectors of A in A
2211 (RWorkspace: need BDSPAC) */
2213 zbdsqr_("U", n, n, m, &c__0, &s[1], &rwork[ie], &a[
2214 a_offset], lda, &u[u_offset], ldu, cdum, &
2215 c__1, &rwork[irwork], info);
2219 } else if (wntvas) {
2221 /* Path 9 (M much larger than N, JOBU='A', JOBVT='S'
2223 M left singular vectors to be computed in U and
2224 N right singular vectors to be computed in VT
2227 i__2 = *n + *m, i__3 = *n * 3;
2228 if (*lwork >= *n * *n + max(i__2,i__3)) {
2230 /* Sufficient workspace for a fast algorithm */
2233 if (*lwork >= wrkbl + *lda * *n) {
2235 /* WORK(IU) is LDA by N */
2240 /* WORK(IU) is N by N */
2244 itau = iu + ldwrku * *n;
2247 /* Compute A=Q*R, copying result to U
2248 (CWorkspace: need N*N+2*N, prefer N*N+N+N*NB)
2251 i__2 = *lwork - iwork + 1;
2252 zgeqrf_(m, n, &a[a_offset], lda, &work[itau], &work[
2253 iwork], &i__2, &ierr);
2254 zlacpy_("L", m, n, &a[a_offset], lda, &u[u_offset],
2258 (CWorkspace: need N*N+N+M, prefer N*N+N+M*NB)
2261 i__2 = *lwork - iwork + 1;
2262 zungqr_(m, m, n, &u[u_offset], ldu, &work[itau], &
2263 work[iwork], &i__2, &ierr);
2265 /* Copy R to WORK(IU), zeroing out below it */
2267 zlacpy_("U", n, n, &a[a_offset], lda, &work[iu], &
2271 zlaset_("L", &i__2, &i__3, &c_b1, &c_b1, &work[iu + 1]
2278 /* Bidiagonalize R in WORK(IU), copying result to VT
2279 (CWorkspace: need N*N+3*N, prefer N*N+2*N+2*N*NB)
2280 (RWorkspace: need N) */
2282 i__2 = *lwork - iwork + 1;
2283 zgebrd_(n, n, &work[iu], &ldwrku, &s[1], &rwork[ie], &
2284 work[itauq], &work[itaup], &work[iwork], &
2286 zlacpy_("U", n, n, &work[iu], &ldwrku, &vt[vt_offset],
2289 /* Generate left bidiagonalizing vectors in WORK(IU)
2290 (CWorkspace: need N*N+3*N, prefer N*N+2*N+N*NB)
2293 i__2 = *lwork - iwork + 1;
2294 zungbr_("Q", n, n, n, &work[iu], &ldwrku, &work[itauq]
2295 , &work[iwork], &i__2, &ierr);
2297 /* Generate right bidiagonalizing vectors in VT
2298 (CWorkspace: need N*N+3*N-1,
2299 prefer N*N+2*N+(N-1)*NB)
2300 (RWorkspace: need 0) */
2302 i__2 = *lwork - iwork + 1;
2303 zungbr_("P", n, n, n, &vt[vt_offset], ldvt, &work[
2304 itaup], &work[iwork], &i__2, &ierr)
2308 /* Perform bidiagonal QR iteration, computing left
2309 singular vectors of R in WORK(IU) and computing
2310 right singular vectors of R in VT
2311 (CWorkspace: need N*N)
2312 (RWorkspace: need BDSPAC) */
2314 zbdsqr_("U", n, n, n, &c__0, &s[1], &rwork[ie], &vt[
2315 vt_offset], ldvt, &work[iu], &ldwrku, cdum, &
2316 c__1, &rwork[irwork], info);
2318 /* Multiply Q in U by left singular vectors of R in
2319 WORK(IU), storing result in A
2320 (CWorkspace: need N*N)
2323 zgemm_("N", "N", m, n, n, &c_b2, &u[u_offset], ldu, &
2324 work[iu], &ldwrku, &c_b1, &a[a_offset], lda);
2326 /* Copy left singular vectors of A from A to U */
2328 zlacpy_("F", m, n, &a[a_offset], lda, &u[u_offset],
2333 /* Insufficient workspace for a fast algorithm */
2338 /* Compute A=Q*R, copying result to U
2339 (CWorkspace: need 2*N, prefer N+N*NB)
2342 i__2 = *lwork - iwork + 1;
2343 zgeqrf_(m, n, &a[a_offset], lda, &work[itau], &work[
2344 iwork], &i__2, &ierr);
2345 zlacpy_("L", m, n, &a[a_offset], lda, &u[u_offset],
2349 (CWorkspace: need N+M, prefer N+M*NB)
2352 i__2 = *lwork - iwork + 1;
2353 zungqr_(m, m, n, &u[u_offset], ldu, &work[itau], &
2354 work[iwork], &i__2, &ierr);
2356 /* Copy R from A to VT, zeroing out below it */
2358 zlacpy_("U", n, n, &a[a_offset], lda, &vt[vt_offset],
2362 zlaset_("L", &i__2, &i__3, &c_b1, &c_b1, &vt_ref(2, 1)
2369 /* Bidiagonalize R in VT
2370 (CWorkspace: need 3*N, prefer 2*N+2*N*NB)
2371 (RWorkspace: need N) */
2373 i__2 = *lwork - iwork + 1;
2374 zgebrd_(n, n, &vt[vt_offset], ldvt, &s[1], &rwork[ie],
2375 &work[itauq], &work[itaup], &work[iwork], &
2378 /* Multiply Q in U by left bidiagonalizing vectors
2380 (CWorkspace: need 2*N+M, prefer 2*N+M*NB)
2383 i__2 = *lwork - iwork + 1;
2384 zunmbr_("Q", "R", "N", m, n, n, &vt[vt_offset], ldvt,
2385 &work[itauq], &u[u_offset], ldu, &work[iwork],
2388 /* Generate right bidiagonalizing vectors in VT
2389 (CWorkspace: need 3*N-1, prefer 2*N+(N-1)*NB)
2392 i__2 = *lwork - iwork + 1;
2393 zungbr_("P", n, n, n, &vt[vt_offset], ldvt, &work[
2394 itaup], &work[iwork], &i__2, &ierr)
2398 /* Perform bidiagonal QR iteration, computing left
2399 singular vectors of A in U and computing right
2400 singular vectors of A in VT
2402 (RWorkspace: need BDSPAC) */
2404 zbdsqr_("U", n, n, m, &c__0, &s[1], &rwork[ie], &vt[
2405 vt_offset], ldvt, &u[u_offset], ldu, cdum, &
2406 c__1, &rwork[irwork], info);
2418 Path 10 (M at least N, but not much larger)
2419 Reduce to bidiagonal form without QR decomposition */
2427 (CWorkspace: need 2*N+M, prefer 2*N+(M+N)*NB)
2428 (RWorkspace: need N) */
2430 i__2 = *lwork - iwork + 1;
2431 zgebrd_(m, n, &a[a_offset], lda, &s[1], &rwork[ie], &work[itauq],
2432 &work[itaup], &work[iwork], &i__2, &ierr);
2435 /* If left singular vectors desired in U, copy result to U
2436 and generate left bidiagonalizing vectors in U
2437 (CWorkspace: need 2*N+NCU, prefer 2*N+NCU*NB)
2440 zlacpy_("L", m, n, &a[a_offset], lda, &u[u_offset], ldu);
2447 i__2 = *lwork - iwork + 1;
2448 zungbr_("Q", m, &ncu, n, &u[u_offset], ldu, &work[itauq], &
2449 work[iwork], &i__2, &ierr);
2453 /* If right singular vectors desired in VT, copy result to
2454 VT and generate right bidiagonalizing vectors in VT
2455 (CWorkspace: need 3*N-1, prefer 2*N+(N-1)*NB)
2458 zlacpy_("U", n, n, &a[a_offset], lda, &vt[vt_offset], ldvt);
2459 i__2 = *lwork - iwork + 1;
2460 zungbr_("P", n, n, n, &vt[vt_offset], ldvt, &work[itaup], &
2461 work[iwork], &i__2, &ierr);
2465 /* If left singular vectors desired in A, generate left
2466 bidiagonalizing vectors in A
2467 (CWorkspace: need 3*N, prefer 2*N+N*NB)
2470 i__2 = *lwork - iwork + 1;
2471 zungbr_("Q", m, n, n, &a[a_offset], lda, &work[itauq], &work[
2472 iwork], &i__2, &ierr);
2476 /* If right singular vectors desired in A, generate right
2477 bidiagonalizing vectors in A
2478 (CWorkspace: need 3*N-1, prefer 2*N+(N-1)*NB)
2481 i__2 = *lwork - iwork + 1;
2482 zungbr_("P", n, n, n, &a[a_offset], lda, &work[itaup], &work[
2483 iwork], &i__2, &ierr);
2486 if (wntuas || wntuo) {
2492 if (wntvas || wntvo) {
2498 if (! wntuo && ! wntvo) {
2500 /* Perform bidiagonal QR iteration, if desired, computing
2501 left singular vectors in U and computing right singular
2504 (RWorkspace: need BDSPAC) */
2506 zbdsqr_("U", n, &ncvt, &nru, &c__0, &s[1], &rwork[ie], &vt[
2507 vt_offset], ldvt, &u[u_offset], ldu, cdum, &c__1, &
2508 rwork[irwork], info);
2509 } else if (! wntuo && wntvo) {
2511 /* Perform bidiagonal QR iteration, if desired, computing
2512 left singular vectors in U and computing right singular
2515 (RWorkspace: need BDSPAC) */
2517 zbdsqr_("U", n, &ncvt, &nru, &c__0, &s[1], &rwork[ie], &a[
2518 a_offset], lda, &u[u_offset], ldu, cdum, &c__1, &
2519 rwork[irwork], info);
2522 /* Perform bidiagonal QR iteration, if desired, computing
2523 left singular vectors in A and computing right singular
2526 (RWorkspace: need BDSPAC) */
2528 zbdsqr_("U", n, &ncvt, &nru, &c__0, &s[1], &rwork[ie], &vt[
2529 vt_offset], ldvt, &a[a_offset], lda, cdum, &c__1, &
2530 rwork[irwork], info);
2537 /* A has more columns than rows. If A has sufficiently more
2538 columns than rows, first reduce using the LQ decomposition (if
2539 sufficient workspace available) */
2545 /* Path 1t(N much larger than M, JOBVT='N')
2546 No right singular vectors to be computed */
2552 (CWorkspace: need 2*M, prefer M+M*NB)
2555 i__2 = *lwork - iwork + 1;
2556 zgelqf_(m, n, &a[a_offset], lda, &work[itau], &work[iwork], &
2559 /* Zero out above L */
2563 zlaset_("U", &i__2, &i__3, &c_b1, &c_b1, &a_ref(1, 2), lda);
2569 /* Bidiagonalize L in A
2570 (CWorkspace: need 3*M, prefer 2*M+2*M*NB)
2571 (RWorkspace: need M) */
2573 i__2 = *lwork - iwork + 1;
2574 zgebrd_(m, m, &a[a_offset], lda, &s[1], &rwork[ie], &work[
2575 itauq], &work[itaup], &work[iwork], &i__2, &ierr);
2576 if (wntuo || wntuas) {
2578 /* If left singular vectors desired, generate Q
2579 (CWorkspace: need 3*M, prefer 2*M+M*NB)
2582 i__2 = *lwork - iwork + 1;
2583 zungbr_("Q", m, m, m, &a[a_offset], lda, &work[itauq], &
2584 work[iwork], &i__2, &ierr);
2588 if (wntuo || wntuas) {
2592 /* Perform bidiagonal QR iteration, computing left singular
2593 vectors of A in A if desired
2595 (RWorkspace: need BDSPAC) */
2597 zbdsqr_("U", m, &c__0, &nru, &c__0, &s[1], &rwork[ie], cdum, &
2598 c__1, &a[a_offset], lda, cdum, &c__1, &rwork[irwork],
2601 /* If left singular vectors desired in U, copy them there */
2604 zlacpy_("F", m, m, &a[a_offset], lda, &u[u_offset], ldu);
2607 } else if (wntvo && wntun) {
2609 /* Path 2t(N much larger than M, JOBU='N', JOBVT='O')
2610 M right singular vectors to be overwritten on A and
2611 no left singular vectors to be computed */
2613 if (*lwork >= *m * *m + *m * 3) {
2615 /* Sufficient workspace for a fast algorithm */
2619 i__2 = wrkbl, i__3 = *lda * *n;
2620 if (*lwork >= max(i__2,i__3) + *lda * *m) {
2622 /* WORK(IU) is LDA by N and WORK(IR) is LDA by M */
2627 } else /* if(complicated condition) */ {
2629 i__2 = wrkbl, i__3 = *lda * *n;
2630 if (*lwork >= max(i__2,i__3) + *m * *m) {
2632 /* WORK(IU) is LDA by N and WORK(IR) is M by M */
2639 /* WORK(IU) is M by CHUNK and WORK(IR) is M by M */
2642 chunk = (*lwork - *m * *m) / *m;
2646 itau = ir + ldwrkr * *m;
2650 (CWorkspace: need M*M+2*M, prefer M*M+M+M*NB)
2653 i__2 = *lwork - iwork + 1;
2654 zgelqf_(m, n, &a[a_offset], lda, &work[itau], &work[iwork]
2657 /* Copy L to WORK(IR) and zero out above it */
2659 zlacpy_("L", m, m, &a[a_offset], lda, &work[ir], &ldwrkr);
2662 zlaset_("U", &i__2, &i__3, &c_b1, &c_b1, &work[ir +
2666 (CWorkspace: need M*M+2*M, prefer M*M+M+M*NB)
2669 i__2 = *lwork - iwork + 1;
2670 zunglq_(m, n, m, &a[a_offset], lda, &work[itau], &work[
2671 iwork], &i__2, &ierr);
2677 /* Bidiagonalize L in WORK(IR)
2678 (CWorkspace: need M*M+3*M, prefer M*M+2*M+2*M*NB)
2679 (RWorkspace: need M) */
2681 i__2 = *lwork - iwork + 1;
2682 zgebrd_(m, m, &work[ir], &ldwrkr, &s[1], &rwork[ie], &
2683 work[itauq], &work[itaup], &work[iwork], &i__2, &
2686 /* Generate right vectors bidiagonalizing L
2687 (CWorkspace: need M*M+3*M-1, prefer M*M+2*M+(M-1)*NB)
2690 i__2 = *lwork - iwork + 1;
2691 zungbr_("P", m, m, m, &work[ir], &ldwrkr, &work[itaup], &
2692 work[iwork], &i__2, &ierr);
2695 /* Perform bidiagonal QR iteration, computing right
2696 singular vectors of L in WORK(IR)
2697 (CWorkspace: need M*M)
2698 (RWorkspace: need BDSPAC) */
2700 zbdsqr_("U", m, m, &c__0, &c__0, &s[1], &rwork[ie], &work[
2701 ir], &ldwrkr, cdum, &c__1, cdum, &c__1, &rwork[
2705 /* Multiply right singular vectors of L in WORK(IR) by Q
2706 in A, storing result in WORK(IU) and copying to A
2707 (CWorkspace: need M*M+M, prefer M*M+M*N)
2712 for (i__ = 1; i__3 < 0 ? i__ >= i__2 : i__ <= i__2; i__ +=
2715 i__4 = *n - i__ + 1;
2716 blk = min(i__4,chunk);
2717 zgemm_("N", "N", m, &blk, m, &c_b2, &work[ir], &
2718 ldwrkr, &a_ref(1, i__), lda, &c_b1, &work[iu],
2720 zlacpy_("F", m, &blk, &work[iu], &ldwrku, &a_ref(1,
2727 /* Insufficient workspace for a fast algorithm */
2735 (CWorkspace: need 2*M+N, prefer 2*M+(M+N)*NB)
2736 (RWorkspace: need M) */
2738 i__3 = *lwork - iwork + 1;
2739 zgebrd_(m, n, &a[a_offset], lda, &s[1], &rwork[ie], &work[
2740 itauq], &work[itaup], &work[iwork], &i__3, &ierr);
2742 /* Generate right vectors bidiagonalizing A
2743 (CWorkspace: need 3*M, prefer 2*M+M*NB)
2746 i__3 = *lwork - iwork + 1;
2747 zungbr_("P", m, n, m, &a[a_offset], lda, &work[itaup], &
2748 work[iwork], &i__3, &ierr);
2751 /* Perform bidiagonal QR iteration, computing right
2752 singular vectors of A in A
2754 (RWorkspace: need BDSPAC) */
2756 zbdsqr_("L", m, n, &c__0, &c__0, &s[1], &rwork[ie], &a[
2757 a_offset], lda, cdum, &c__1, cdum, &c__1, &rwork[
2762 } else if (wntvo && wntuas) {
2764 /* Path 3t(N much larger than M, JOBU='S' or 'A', JOBVT='O')
2765 M right singular vectors to be overwritten on A and
2766 M left singular vectors to be computed in U */
2768 if (*lwork >= *m * *m + *m * 3) {
2770 /* Sufficient workspace for a fast algorithm */
2774 i__3 = wrkbl, i__2 = *lda * *n;
2775 if (*lwork >= max(i__3,i__2) + *lda * *m) {
2777 /* WORK(IU) is LDA by N and WORK(IR) is LDA by M */
2782 } else /* if(complicated condition) */ {
2784 i__3 = wrkbl, i__2 = *lda * *n;
2785 if (*lwork >= max(i__3,i__2) + *m * *m) {
2787 /* WORK(IU) is LDA by N and WORK(IR) is M by M */
2794 /* WORK(IU) is M by CHUNK and WORK(IR) is M by M */
2797 chunk = (*lwork - *m * *m) / *m;
2801 itau = ir + ldwrkr * *m;
2805 (CWorkspace: need M*M+2*M, prefer M*M+M+M*NB)
2808 i__3 = *lwork - iwork + 1;
2809 zgelqf_(m, n, &a[a_offset], lda, &work[itau], &work[iwork]
2812 /* Copy L to U, zeroing about above it */
2814 zlacpy_("L", m, m, &a[a_offset], lda, &u[u_offset], ldu);
2817 zlaset_("U", &i__3, &i__2, &c_b1, &c_b1, &u_ref(1, 2),
2821 (CWorkspace: need M*M+2*M, prefer M*M+M+M*NB)
2824 i__3 = *lwork - iwork + 1;
2825 zunglq_(m, n, m, &a[a_offset], lda, &work[itau], &work[
2826 iwork], &i__3, &ierr);
2832 /* Bidiagonalize L in U, copying result to WORK(IR)
2833 (CWorkspace: need M*M+3*M, prefer M*M+2*M+2*M*NB)
2834 (RWorkspace: need M) */
2836 i__3 = *lwork - iwork + 1;
2837 zgebrd_(m, m, &u[u_offset], ldu, &s[1], &rwork[ie], &work[
2838 itauq], &work[itaup], &work[iwork], &i__3, &ierr);
2839 zlacpy_("U", m, m, &u[u_offset], ldu, &work[ir], &ldwrkr);
2841 /* Generate right vectors bidiagonalizing L in WORK(IR)
2842 (CWorkspace: need M*M+3*M-1, prefer M*M+2*M+(M-1)*NB)
2845 i__3 = *lwork - iwork + 1;
2846 zungbr_("P", m, m, m, &work[ir], &ldwrkr, &work[itaup], &
2847 work[iwork], &i__3, &ierr);
2849 /* Generate left vectors bidiagonalizing L in U
2850 (CWorkspace: need M*M+3*M, prefer M*M+2*M+M*NB)
2853 i__3 = *lwork - iwork + 1;
2854 zungbr_("Q", m, m, m, &u[u_offset], ldu, &work[itauq], &
2855 work[iwork], &i__3, &ierr);
2858 /* Perform bidiagonal QR iteration, computing left
2859 singular vectors of L in U, and computing right
2860 singular vectors of L in WORK(IR)
2861 (CWorkspace: need M*M)
2862 (RWorkspace: need BDSPAC) */
2864 zbdsqr_("U", m, m, m, &c__0, &s[1], &rwork[ie], &work[ir],
2865 &ldwrkr, &u[u_offset], ldu, cdum, &c__1, &rwork[
2869 /* Multiply right singular vectors of L in WORK(IR) by Q
2870 in A, storing result in WORK(IU) and copying to A
2871 (CWorkspace: need M*M+M, prefer M*M+M*N))
2876 for (i__ = 1; i__2 < 0 ? i__ >= i__3 : i__ <= i__3; i__ +=
2879 i__4 = *n - i__ + 1;
2880 blk = min(i__4,chunk);
2881 zgemm_("N", "N", m, &blk, m, &c_b2, &work[ir], &
2882 ldwrkr, &a_ref(1, i__), lda, &c_b1, &work[iu],
2884 zlacpy_("F", m, &blk, &work[iu], &ldwrku, &a_ref(1,
2891 /* Insufficient workspace for a fast algorithm */
2897 (CWorkspace: need 2*M, prefer M+M*NB)
2900 i__2 = *lwork - iwork + 1;
2901 zgelqf_(m, n, &a[a_offset], lda, &work[itau], &work[iwork]
2904 /* Copy L to U, zeroing out above it */
2906 zlacpy_("L", m, m, &a[a_offset], lda, &u[u_offset], ldu);
2909 zlaset_("U", &i__2, &i__3, &c_b1, &c_b1, &u_ref(1, 2),
2913 (CWorkspace: need 2*M, prefer M+M*NB)
2916 i__2 = *lwork - iwork + 1;
2917 zunglq_(m, n, m, &a[a_offset], lda, &work[itau], &work[
2918 iwork], &i__2, &ierr);
2924 /* Bidiagonalize L in U
2925 (CWorkspace: need 3*M, prefer 2*M+2*M*NB)
2926 (RWorkspace: need M) */
2928 i__2 = *lwork - iwork + 1;
2929 zgebrd_(m, m, &u[u_offset], ldu, &s[1], &rwork[ie], &work[
2930 itauq], &work[itaup], &work[iwork], &i__2, &ierr);
2932 /* Multiply right vectors bidiagonalizing L by Q in A
2933 (CWorkspace: need 2*M+N, prefer 2*M+N*NB)
2936 i__2 = *lwork - iwork + 1;
2937 zunmbr_("P", "L", "C", m, n, m, &u[u_offset], ldu, &work[
2938 itaup], &a[a_offset], lda, &work[iwork], &i__2, &
2941 /* Generate left vectors bidiagonalizing L in U
2942 (CWorkspace: need 3*M, prefer 2*M+M*NB)
2945 i__2 = *lwork - iwork + 1;
2946 zungbr_("Q", m, m, m, &u[u_offset], ldu, &work[itauq], &
2947 work[iwork], &i__2, &ierr);
2950 /* Perform bidiagonal QR iteration, computing left
2951 singular vectors of A in U and computing right
2952 singular vectors of A in A
2954 (RWorkspace: need BDSPAC) */
2956 zbdsqr_("U", m, n, m, &c__0, &s[1], &rwork[ie], &a[
2957 a_offset], lda, &u[u_offset], ldu, cdum, &c__1, &
2958 rwork[irwork], info);
2966 /* Path 4t(N much larger than M, JOBU='N', JOBVT='S')
2967 M right singular vectors to be computed in VT and
2968 no left singular vectors to be computed */
2970 if (*lwork >= *m * *m + *m * 3) {
2972 /* Sufficient workspace for a fast algorithm */
2975 if (*lwork >= wrkbl + *lda * *m) {
2977 /* WORK(IR) is LDA by M */
2982 /* WORK(IR) is M by M */
2986 itau = ir + ldwrkr * *m;
2990 (CWorkspace: need M*M+2*M, prefer M*M+M+M*NB)
2993 i__2 = *lwork - iwork + 1;
2994 zgelqf_(m, n, &a[a_offset], lda, &work[itau], &work[
2995 iwork], &i__2, &ierr);
2997 /* Copy L to WORK(IR), zeroing out above it */
2999 zlacpy_("L", m, m, &a[a_offset], lda, &work[ir], &
3003 zlaset_("U", &i__2, &i__3, &c_b1, &c_b1, &work[ir +
3007 (CWorkspace: need M*M+2*M, prefer M*M+M+M*NB)
3010 i__2 = *lwork - iwork + 1;
3011 zunglq_(m, n, m, &a[a_offset], lda, &work[itau], &
3012 work[iwork], &i__2, &ierr);
3018 /* Bidiagonalize L in WORK(IR)
3019 (CWorkspace: need M*M+3*M, prefer M*M+2*M+2*M*NB)
3020 (RWorkspace: need M) */
3022 i__2 = *lwork - iwork + 1;
3023 zgebrd_(m, m, &work[ir], &ldwrkr, &s[1], &rwork[ie], &
3024 work[itauq], &work[itaup], &work[iwork], &
3027 /* Generate right vectors bidiagonalizing L in
3029 (CWorkspace: need M*M+3*M, prefer M*M+2*M+(M-1)*NB)
3032 i__2 = *lwork - iwork + 1;
3033 zungbr_("P", m, m, m, &work[ir], &ldwrkr, &work[itaup]
3034 , &work[iwork], &i__2, &ierr);
3037 /* Perform bidiagonal QR iteration, computing right
3038 singular vectors of L in WORK(IR)
3039 (CWorkspace: need M*M)
3040 (RWorkspace: need BDSPAC) */
3042 zbdsqr_("U", m, m, &c__0, &c__0, &s[1], &rwork[ie], &
3043 work[ir], &ldwrkr, cdum, &c__1, cdum, &c__1, &
3044 rwork[irwork], info);
3046 /* Multiply right singular vectors of L in WORK(IR) by
3047 Q in A, storing result in VT
3048 (CWorkspace: need M*M)
3051 zgemm_("N", "N", m, n, m, &c_b2, &work[ir], &ldwrkr, &
3052 a[a_offset], lda, &c_b1, &vt[vt_offset], ldvt);
3056 /* Insufficient workspace for a fast algorithm */
3062 (CWorkspace: need 2*M, prefer M+M*NB)
3065 i__2 = *lwork - iwork + 1;
3066 zgelqf_(m, n, &a[a_offset], lda, &work[itau], &work[
3067 iwork], &i__2, &ierr);
3069 /* Copy result to VT */
3071 zlacpy_("U", m, n, &a[a_offset], lda, &vt[vt_offset],
3075 (CWorkspace: need 2*M, prefer M+M*NB)
3078 i__2 = *lwork - iwork + 1;
3079 zunglq_(m, n, m, &vt[vt_offset], ldvt, &work[itau], &
3080 work[iwork], &i__2, &ierr);
3086 /* Zero out above L in A */
3090 zlaset_("U", &i__2, &i__3, &c_b1, &c_b1, &a_ref(1, 2),
3093 /* Bidiagonalize L in A
3094 (CWorkspace: need 3*M, prefer 2*M+2*M*NB)
3095 (RWorkspace: need M) */
3097 i__2 = *lwork - iwork + 1;
3098 zgebrd_(m, m, &a[a_offset], lda, &s[1], &rwork[ie], &
3099 work[itauq], &work[itaup], &work[iwork], &
3102 /* Multiply right vectors bidiagonalizing L by Q in VT
3103 (CWorkspace: need 2*M+N, prefer 2*M+N*NB)
3106 i__2 = *lwork - iwork + 1;
3107 zunmbr_("P", "L", "C", m, n, m, &a[a_offset], lda, &
3108 work[itaup], &vt[vt_offset], ldvt, &work[
3109 iwork], &i__2, &ierr);
3112 /* Perform bidiagonal QR iteration, computing right
3113 singular vectors of A in VT
3115 (RWorkspace: need BDSPAC) */
3117 zbdsqr_("U", m, n, &c__0, &c__0, &s[1], &rwork[ie], &
3118 vt[vt_offset], ldvt, cdum, &c__1, cdum, &c__1,
3119 &rwork[irwork], info);
3125 /* Path 5t(N much larger than M, JOBU='O', JOBVT='S')
3126 M right singular vectors to be computed in VT and
3127 M left singular vectors to be overwritten on A */
3129 if (*lwork >= (*m << 1) * *m + *m * 3) {
3131 /* Sufficient workspace for a fast algorithm */
3134 if (*lwork >= wrkbl + (*lda << 1) * *m) {
3136 /* WORK(IU) is LDA by M and WORK(IR) is LDA by M */
3139 ir = iu + ldwrku * *m;
3141 } else if (*lwork >= wrkbl + (*lda + *m) * *m) {
3143 /* WORK(IU) is LDA by M and WORK(IR) is M by M */
3146 ir = iu + ldwrku * *m;
3150 /* WORK(IU) is M by M and WORK(IR) is M by M */
3153 ir = iu + ldwrku * *m;
3156 itau = ir + ldwrkr * *m;
3160 (CWorkspace: need 2*M*M+2*M, prefer 2*M*M+M+M*NB)
3163 i__2 = *lwork - iwork + 1;
3164 zgelqf_(m, n, &a[a_offset], lda, &work[itau], &work[
3165 iwork], &i__2, &ierr);
3167 /* Copy L to WORK(IU), zeroing out below it */
3169 zlacpy_("L", m, m, &a[a_offset], lda, &work[iu], &
3173 zlaset_("U", &i__2, &i__3, &c_b1, &c_b1, &work[iu +
3177 (CWorkspace: need 2*M*M+2*M, prefer 2*M*M+M+M*NB)
3180 i__2 = *lwork - iwork + 1;
3181 zunglq_(m, n, m, &a[a_offset], lda, &work[itau], &
3182 work[iwork], &i__2, &ierr);
3188 /* Bidiagonalize L in WORK(IU), copying result to
3190 (CWorkspace: need 2*M*M+3*M,
3191 prefer 2*M*M+2*M+2*M*NB)
3192 (RWorkspace: need M) */
3194 i__2 = *lwork - iwork + 1;
3195 zgebrd_(m, m, &work[iu], &ldwrku, &s[1], &rwork[ie], &
3196 work[itauq], &work[itaup], &work[iwork], &
3198 zlacpy_("L", m, m, &work[iu], &ldwrku, &work[ir], &
3201 /* Generate right bidiagonalizing vectors in WORK(IU)
3202 (CWorkspace: need 2*M*M+3*M-1,
3203 prefer 2*M*M+2*M+(M-1)*NB)
3206 i__2 = *lwork - iwork + 1;
3207 zungbr_("P", m, m, m, &work[iu], &ldwrku, &work[itaup]
3208 , &work[iwork], &i__2, &ierr);
3210 /* Generate left bidiagonalizing vectors in WORK(IR)
3211 (CWorkspace: need 2*M*M+3*M, prefer 2*M*M+2*M+M*NB)
3214 i__2 = *lwork - iwork + 1;
3215 zungbr_("Q", m, m, m, &work[ir], &ldwrkr, &work[itauq]
3216 , &work[iwork], &i__2, &ierr);
3219 /* Perform bidiagonal QR iteration, computing left
3220 singular vectors of L in WORK(IR) and computing
3221 right singular vectors of L in WORK(IU)
3222 (CWorkspace: need 2*M*M)
3223 (RWorkspace: need BDSPAC) */
3225 zbdsqr_("U", m, m, m, &c__0, &s[1], &rwork[ie], &work[
3226 iu], &ldwrku, &work[ir], &ldwrkr, cdum, &c__1,
3227 &rwork[irwork], info);
3229 /* Multiply right singular vectors of L in WORK(IU) by
3230 Q in A, storing result in VT
3231 (CWorkspace: need M*M)
3234 zgemm_("N", "N", m, n, m, &c_b2, &work[iu], &ldwrku, &
3235 a[a_offset], lda, &c_b1, &vt[vt_offset], ldvt);
3237 /* Copy left singular vectors of L to A
3238 (CWorkspace: need M*M)
3241 zlacpy_("F", m, m, &work[ir], &ldwrkr, &a[a_offset],
3246 /* Insufficient workspace for a fast algorithm */
3251 /* Compute A=L*Q, copying result to VT
3252 (CWorkspace: need 2*M, prefer M+M*NB)
3255 i__2 = *lwork - iwork + 1;
3256 zgelqf_(m, n, &a[a_offset], lda, &work[itau], &work[
3257 iwork], &i__2, &ierr);
3258 zlacpy_("U", m, n, &a[a_offset], lda, &vt[vt_offset],
3262 (CWorkspace: need 2*M, prefer M+M*NB)
3265 i__2 = *lwork - iwork + 1;
3266 zunglq_(m, n, m, &vt[vt_offset], ldvt, &work[itau], &
3267 work[iwork], &i__2, &ierr);
3273 /* Zero out above L in A */
3277 zlaset_("U", &i__2, &i__3, &c_b1, &c_b1, &a_ref(1, 2),
3280 /* Bidiagonalize L in A
3281 (CWorkspace: need 3*M, prefer 2*M+2*M*NB)
3282 (RWorkspace: need M) */
3284 i__2 = *lwork - iwork + 1;
3285 zgebrd_(m, m, &a[a_offset], lda, &s[1], &rwork[ie], &
3286 work[itauq], &work[itaup], &work[iwork], &
3289 /* Multiply right vectors bidiagonalizing L by Q in VT
3290 (CWorkspace: need 2*M+N, prefer 2*M+N*NB)
3293 i__2 = *lwork - iwork + 1;
3294 zunmbr_("P", "L", "C", m, n, m, &a[a_offset], lda, &
3295 work[itaup], &vt[vt_offset], ldvt, &work[
3296 iwork], &i__2, &ierr);
3298 /* Generate left bidiagonalizing vectors of L in A
3299 (CWorkspace: need 3*M, prefer 2*M+M*NB)
3302 i__2 = *lwork - iwork + 1;
3303 zungbr_("Q", m, m, m, &a[a_offset], lda, &work[itauq],
3304 &work[iwork], &i__2, &ierr);
3307 /* Perform bidiagonal QR iteration, computing left
3308 singular vectors of A in A and computing right
3309 singular vectors of A in VT
3311 (RWorkspace: need BDSPAC) */
3313 zbdsqr_("U", m, n, m, &c__0, &s[1], &rwork[ie], &vt[
3314 vt_offset], ldvt, &a[a_offset], lda, cdum, &
3315 c__1, &rwork[irwork], info);
3319 } else if (wntuas) {
3321 /* Path 6t(N much larger than M, JOBU='S' or 'A',
3323 M right singular vectors to be computed in VT and
3324 M left singular vectors to be computed in U */
3326 if (*lwork >= *m * *m + *m * 3) {
3328 /* Sufficient workspace for a fast algorithm */
3331 if (*lwork >= wrkbl + *lda * *m) {
3333 /* WORK(IU) is LDA by N */
3338 /* WORK(IU) is LDA by M */
3342 itau = iu + ldwrku * *m;
3346 (CWorkspace: need M*M+2*M, prefer M*M+M+M*NB)
3349 i__2 = *lwork - iwork + 1;
3350 zgelqf_(m, n, &a[a_offset], lda, &work[itau], &work[
3351 iwork], &i__2, &ierr);
3353 /* Copy L to WORK(IU), zeroing out above it */
3355 zlacpy_("L", m, m, &a[a_offset], lda, &work[iu], &
3359 zlaset_("U", &i__2, &i__3, &c_b1, &c_b1, &work[iu +
3363 (CWorkspace: need M*M+2*M, prefer M*M+M+M*NB)
3366 i__2 = *lwork - iwork + 1;
3367 zunglq_(m, n, m, &a[a_offset], lda, &work[itau], &
3368 work[iwork], &i__2, &ierr);
3374 /* Bidiagonalize L in WORK(IU), copying result to U
3375 (CWorkspace: need M*M+3*M, prefer M*M+2*M+2*M*NB)
3376 (RWorkspace: need M) */
3378 i__2 = *lwork - iwork + 1;
3379 zgebrd_(m, m, &work[iu], &ldwrku, &s[1], &rwork[ie], &
3380 work[itauq], &work[itaup], &work[iwork], &
3382 zlacpy_("L", m, m, &work[iu], &ldwrku, &u[u_offset],
3385 /* Generate right bidiagonalizing vectors in WORK(IU)
3386 (CWorkspace: need M*M+3*M-1,
3387 prefer M*M+2*M+(M-1)*NB)
3390 i__2 = *lwork - iwork + 1;
3391 zungbr_("P", m, m, m, &work[iu], &ldwrku, &work[itaup]
3392 , &work[iwork], &i__2, &ierr);
3394 /* Generate left bidiagonalizing vectors in U
3395 (CWorkspace: need M*M+3*M, prefer M*M+2*M+M*NB)
3398 i__2 = *lwork - iwork + 1;
3399 zungbr_("Q", m, m, m, &u[u_offset], ldu, &work[itauq],
3400 &work[iwork], &i__2, &ierr);
3403 /* Perform bidiagonal QR iteration, computing left
3404 singular vectors of L in U and computing right
3405 singular vectors of L in WORK(IU)
3406 (CWorkspace: need M*M)
3407 (RWorkspace: need BDSPAC) */
3409 zbdsqr_("U", m, m, m, &c__0, &s[1], &rwork[ie], &work[
3410 iu], &ldwrku, &u[u_offset], ldu, cdum, &c__1,
3411 &rwork[irwork], info);
3413 /* Multiply right singular vectors of L in WORK(IU) by
3414 Q in A, storing result in VT
3415 (CWorkspace: need M*M)
3418 zgemm_("N", "N", m, n, m, &c_b2, &work[iu], &ldwrku, &
3419 a[a_offset], lda, &c_b1, &vt[vt_offset], ldvt);
3423 /* Insufficient workspace for a fast algorithm */
3428 /* Compute A=L*Q, copying result to VT
3429 (CWorkspace: need 2*M, prefer M+M*NB)
3432 i__2 = *lwork - iwork + 1;
3433 zgelqf_(m, n, &a[a_offset], lda, &work[itau], &work[
3434 iwork], &i__2, &ierr);
3435 zlacpy_("U", m, n, &a[a_offset], lda, &vt[vt_offset],
3439 (CWorkspace: need 2*M, prefer M+M*NB)
3442 i__2 = *lwork - iwork + 1;
3443 zunglq_(m, n, m, &vt[vt_offset], ldvt, &work[itau], &
3444 work[iwork], &i__2, &ierr);
3446 /* Copy L to U, zeroing out above it */
3448 zlacpy_("L", m, m, &a[a_offset], lda, &u[u_offset],
3452 zlaset_("U", &i__2, &i__3, &c_b1, &c_b1, &u_ref(1, 2),
3459 /* Bidiagonalize L in U
3460 (CWorkspace: need 3*M, prefer 2*M+2*M*NB)
3461 (RWorkspace: need M) */
3463 i__2 = *lwork - iwork + 1;
3464 zgebrd_(m, m, &u[u_offset], ldu, &s[1], &rwork[ie], &
3465 work[itauq], &work[itaup], &work[iwork], &
3468 /* Multiply right bidiagonalizing vectors in U by Q
3470 (CWorkspace: need 2*M+N, prefer 2*M+N*NB)
3473 i__2 = *lwork - iwork + 1;
3474 zunmbr_("P", "L", "C", m, n, m, &u[u_offset], ldu, &
3475 work[itaup], &vt[vt_offset], ldvt, &work[
3476 iwork], &i__2, &ierr);
3478 /* Generate left bidiagonalizing vectors in U
3479 (CWorkspace: need 3*M, prefer 2*M+M*NB)
3482 i__2 = *lwork - iwork + 1;
3483 zungbr_("Q", m, m, m, &u[u_offset], ldu, &work[itauq],
3484 &work[iwork], &i__2, &ierr);
3487 /* Perform bidiagonal QR iteration, computing left
3488 singular vectors of A in U and computing right
3489 singular vectors of A in VT
3491 (RWorkspace: need BDSPAC) */
3493 zbdsqr_("U", m, n, m, &c__0, &s[1], &rwork[ie], &vt[
3494 vt_offset], ldvt, &u[u_offset], ldu, cdum, &
3495 c__1, &rwork[irwork], info);
3505 /* Path 7t(N much larger than M, JOBU='N', JOBVT='A')
3506 N right singular vectors to be computed in VT and
3507 no left singular vectors to be computed
3510 i__2 = *n + *m, i__3 = *m * 3;
3511 if (*lwork >= *m * *m + max(i__2,i__3)) {
3513 /* Sufficient workspace for a fast algorithm */
3516 if (*lwork >= wrkbl + *lda * *m) {
3518 /* WORK(IR) is LDA by M */
3523 /* WORK(IR) is M by M */
3527 itau = ir + ldwrkr * *m;
3530 /* Compute A=L*Q, copying result to VT
3531 (CWorkspace: need M*M+2*M, prefer M*M+M+M*NB)
3534 i__2 = *lwork - iwork + 1;
3535 zgelqf_(m, n, &a[a_offset], lda, &work[itau], &work[
3536 iwork], &i__2, &ierr);
3537 zlacpy_("U", m, n, &a[a_offset], lda, &vt[vt_offset],
3540 /* Copy L to WORK(IR), zeroing out above it */
3542 zlacpy_("L", m, m, &a[a_offset], lda, &work[ir], &
3546 zlaset_("U", &i__2, &i__3, &c_b1, &c_b1, &work[ir +
3550 (CWorkspace: need M*M+M+N, prefer M*M+M+N*NB)
3553 i__2 = *lwork - iwork + 1;
3554 zunglq_(n, n, m, &vt[vt_offset], ldvt, &work[itau], &
3555 work[iwork], &i__2, &ierr);
3561 /* Bidiagonalize L in WORK(IR)
3562 (CWorkspace: need M*M+3*M, prefer M*M+2*M+2*M*NB)
3563 (RWorkspace: need M) */
3565 i__2 = *lwork - iwork + 1;
3566 zgebrd_(m, m, &work[ir], &ldwrkr, &s[1], &rwork[ie], &
3567 work[itauq], &work[itaup], &work[iwork], &
3570 /* Generate right bidiagonalizing vectors in WORK(IR)
3571 (CWorkspace: need M*M+3*M-1,
3572 prefer M*M+2*M+(M-1)*NB)
3575 i__2 = *lwork - iwork + 1;
3576 zungbr_("P", m, m, m, &work[ir], &ldwrkr, &work[itaup]
3577 , &work[iwork], &i__2, &ierr);
3580 /* Perform bidiagonal QR iteration, computing right
3581 singular vectors of L in WORK(IR)
3582 (CWorkspace: need M*M)
3583 (RWorkspace: need BDSPAC) */
3585 zbdsqr_("U", m, m, &c__0, &c__0, &s[1], &rwork[ie], &
3586 work[ir], &ldwrkr, cdum, &c__1, cdum, &c__1, &
3587 rwork[irwork], info);
3589 /* Multiply right singular vectors of L in WORK(IR) by
3590 Q in VT, storing result in A
3591 (CWorkspace: need M*M)
3594 zgemm_("N", "N", m, n, m, &c_b2, &work[ir], &ldwrkr, &
3595 vt[vt_offset], ldvt, &c_b1, &a[a_offset], lda);
3597 /* Copy right singular vectors of A from A to VT */
3599 zlacpy_("F", m, n, &a[a_offset], lda, &vt[vt_offset],
3604 /* Insufficient workspace for a fast algorithm */
3609 /* Compute A=L*Q, copying result to VT
3610 (CWorkspace: need 2*M, prefer M+M*NB)
3613 i__2 = *lwork - iwork + 1;
3614 zgelqf_(m, n, &a[a_offset], lda, &work[itau], &work[
3615 iwork], &i__2, &ierr);
3616 zlacpy_("U", m, n, &a[a_offset], lda, &vt[vt_offset],
3620 (CWorkspace: need M+N, prefer M+N*NB)
3623 i__2 = *lwork - iwork + 1;
3624 zunglq_(n, n, m, &vt[vt_offset], ldvt, &work[itau], &
3625 work[iwork], &i__2, &ierr);
3631 /* Zero out above L in A */
3635 zlaset_("U", &i__2, &i__3, &c_b1, &c_b1, &a_ref(1, 2),
3638 /* Bidiagonalize L in A
3639 (CWorkspace: need 3*M, prefer 2*M+2*M*NB)
3640 (RWorkspace: need M) */
3642 i__2 = *lwork - iwork + 1;
3643 zgebrd_(m, m, &a[a_offset], lda, &s[1], &rwork[ie], &
3644 work[itauq], &work[itaup], &work[iwork], &
3647 /* Multiply right bidiagonalizing vectors in A by Q
3649 (CWorkspace: need 2*M+N, prefer 2*M+N*NB)
3652 i__2 = *lwork - iwork + 1;
3653 zunmbr_("P", "L", "C", m, n, m, &a[a_offset], lda, &
3654 work[itaup], &vt[vt_offset], ldvt, &work[
3655 iwork], &i__2, &ierr);
3658 /* Perform bidiagonal QR iteration, computing right
3659 singular vectors of A in VT
3661 (RWorkspace: need BDSPAC) */
3663 zbdsqr_("U", m, n, &c__0, &c__0, &s[1], &rwork[ie], &
3664 vt[vt_offset], ldvt, cdum, &c__1, cdum, &c__1,
3665 &rwork[irwork], info);
3671 /* Path 8t(N much larger than M, JOBU='O', JOBVT='A')
3672 N right singular vectors to be computed in VT and
3673 M left singular vectors to be overwritten on A
3676 i__2 = *n + *m, i__3 = *m * 3;
3677 if (*lwork >= (*m << 1) * *m + max(i__2,i__3)) {
3679 /* Sufficient workspace for a fast algorithm */
3682 if (*lwork >= wrkbl + (*lda << 1) * *m) {
3684 /* WORK(IU) is LDA by M and WORK(IR) is LDA by M */
3687 ir = iu + ldwrku * *m;
3689 } else if (*lwork >= wrkbl + (*lda + *m) * *m) {
3691 /* WORK(IU) is LDA by M and WORK(IR) is M by M */
3694 ir = iu + ldwrku * *m;
3698 /* WORK(IU) is M by M and WORK(IR) is M by M */
3701 ir = iu + ldwrku * *m;
3704 itau = ir + ldwrkr * *m;
3707 /* Compute A=L*Q, copying result to VT
3708 (CWorkspace: need 2*M*M+2*M, prefer 2*M*M+M+M*NB)
3711 i__2 = *lwork - iwork + 1;
3712 zgelqf_(m, n, &a[a_offset], lda, &work[itau], &work[
3713 iwork], &i__2, &ierr);
3714 zlacpy_("U", m, n, &a[a_offset], lda, &vt[vt_offset],
3718 (CWorkspace: need 2*M*M+M+N, prefer 2*M*M+M+N*NB)
3721 i__2 = *lwork - iwork + 1;
3722 zunglq_(n, n, m, &vt[vt_offset], ldvt, &work[itau], &
3723 work[iwork], &i__2, &ierr);
3725 /* Copy L to WORK(IU), zeroing out above it */
3727 zlacpy_("L", m, m, &a[a_offset], lda, &work[iu], &
3731 zlaset_("U", &i__2, &i__3, &c_b1, &c_b1, &work[iu +
3738 /* Bidiagonalize L in WORK(IU), copying result to
3740 (CWorkspace: need 2*M*M+3*M,
3741 prefer 2*M*M+2*M+2*M*NB)
3742 (RWorkspace: need M) */
3744 i__2 = *lwork - iwork + 1;
3745 zgebrd_(m, m, &work[iu], &ldwrku, &s[1], &rwork[ie], &
3746 work[itauq], &work[itaup], &work[iwork], &
3748 zlacpy_("L", m, m, &work[iu], &ldwrku, &work[ir], &
3751 /* Generate right bidiagonalizing vectors in WORK(IU)
3752 (CWorkspace: need 2*M*M+3*M-1,
3753 prefer 2*M*M+2*M+(M-1)*NB)
3756 i__2 = *lwork - iwork + 1;
3757 zungbr_("P", m, m, m, &work[iu], &ldwrku, &work[itaup]
3758 , &work[iwork], &i__2, &ierr);
3760 /* Generate left bidiagonalizing vectors in WORK(IR)
3761 (CWorkspace: need 2*M*M+3*M, prefer 2*M*M+2*M+M*NB)
3764 i__2 = *lwork - iwork + 1;
3765 zungbr_("Q", m, m, m, &work[ir], &ldwrkr, &work[itauq]
3766 , &work[iwork], &i__2, &ierr);
3769 /* Perform bidiagonal QR iteration, computing left
3770 singular vectors of L in WORK(IR) and computing
3771 right singular vectors of L in WORK(IU)
3772 (CWorkspace: need 2*M*M)
3773 (RWorkspace: need BDSPAC) */
3775 zbdsqr_("U", m, m, m, &c__0, &s[1], &rwork[ie], &work[
3776 iu], &ldwrku, &work[ir], &ldwrkr, cdum, &c__1,
3777 &rwork[irwork], info);
3779 /* Multiply right singular vectors of L in WORK(IU) by
3780 Q in VT, storing result in A
3781 (CWorkspace: need M*M)
3784 zgemm_("N", "N", m, n, m, &c_b2, &work[iu], &ldwrku, &
3785 vt[vt_offset], ldvt, &c_b1, &a[a_offset], lda);
3787 /* Copy right singular vectors of A from A to VT */
3789 zlacpy_("F", m, n, &a[a_offset], lda, &vt[vt_offset],
3792 /* Copy left singular vectors of A from WORK(IR) to A */
3794 zlacpy_("F", m, m, &work[ir], &ldwrkr, &a[a_offset],
3799 /* Insufficient workspace for a fast algorithm */
3804 /* Compute A=L*Q, copying result to VT
3805 (CWorkspace: need 2*M, prefer M+M*NB)
3808 i__2 = *lwork - iwork + 1;
3809 zgelqf_(m, n, &a[a_offset], lda, &work[itau], &work[
3810 iwork], &i__2, &ierr);
3811 zlacpy_("U", m, n, &a[a_offset], lda, &vt[vt_offset],
3815 (CWorkspace: need M+N, prefer M+N*NB)
3818 i__2 = *lwork - iwork + 1;
3819 zunglq_(n, n, m, &vt[vt_offset], ldvt, &work[itau], &
3820 work[iwork], &i__2, &ierr);
3826 /* Zero out above L in A */
3830 zlaset_("U", &i__2, &i__3, &c_b1, &c_b1, &a_ref(1, 2),
3833 /* Bidiagonalize L in A
3834 (CWorkspace: need 3*M, prefer 2*M+2*M*NB)
3835 (RWorkspace: need M) */
3837 i__2 = *lwork - iwork + 1;
3838 zgebrd_(m, m, &a[a_offset], lda, &s[1], &rwork[ie], &
3839 work[itauq], &work[itaup], &work[iwork], &
3842 /* Multiply right bidiagonalizing vectors in A by Q
3844 (CWorkspace: need 2*M+N, prefer 2*M+N*NB)
3847 i__2 = *lwork - iwork + 1;
3848 zunmbr_("P", "L", "C", m, n, m, &a[a_offset], lda, &
3849 work[itaup], &vt[vt_offset], ldvt, &work[
3850 iwork], &i__2, &ierr);
3852 /* Generate left bidiagonalizing vectors in A
3853 (CWorkspace: need 3*M, prefer 2*M+M*NB)
3856 i__2 = *lwork - iwork + 1;
3857 zungbr_("Q", m, m, m, &a[a_offset], lda, &work[itauq],
3858 &work[iwork], &i__2, &ierr);
3861 /* Perform bidiagonal QR iteration, computing left
3862 singular vectors of A in A and computing right
3863 singular vectors of A in VT
3865 (RWorkspace: need BDSPAC) */
3867 zbdsqr_("U", m, n, m, &c__0, &s[1], &rwork[ie], &vt[
3868 vt_offset], ldvt, &a[a_offset], lda, cdum, &
3869 c__1, &rwork[irwork], info);
3873 } else if (wntuas) {
3875 /* Path 9t(N much larger than M, JOBU='S' or 'A',
3877 N right singular vectors to be computed in VT and
3878 M left singular vectors to be computed in U
3881 i__2 = *n + *m, i__3 = *m * 3;
3882 if (*lwork >= *m * *m + max(i__2,i__3)) {
3884 /* Sufficient workspace for a fast algorithm */
3887 if (*lwork >= wrkbl + *lda * *m) {
3889 /* WORK(IU) is LDA by M */
3894 /* WORK(IU) is M by M */
3898 itau = iu + ldwrku * *m;
3901 /* Compute A=L*Q, copying result to VT
3902 (CWorkspace: need M*M+2*M, prefer M*M+M+M*NB)
3905 i__2 = *lwork - iwork + 1;
3906 zgelqf_(m, n, &a[a_offset], lda, &work[itau], &work[
3907 iwork], &i__2, &ierr);
3908 zlacpy_("U", m, n, &a[a_offset], lda, &vt[vt_offset],
3912 (CWorkspace: need M*M+M+N, prefer M*M+M+N*NB)
3915 i__2 = *lwork - iwork + 1;
3916 zunglq_(n, n, m, &vt[vt_offset], ldvt, &work[itau], &
3917 work[iwork], &i__2, &ierr);
3919 /* Copy L to WORK(IU), zeroing out above it */
3921 zlacpy_("L", m, m, &a[a_offset], lda, &work[iu], &
3925 zlaset_("U", &i__2, &i__3, &c_b1, &c_b1, &work[iu +
3932 /* Bidiagonalize L in WORK(IU), copying result to U
3933 (CWorkspace: need M*M+3*M, prefer M*M+2*M+2*M*NB)
3934 (RWorkspace: need M) */
3936 i__2 = *lwork - iwork + 1;
3937 zgebrd_(m, m, &work[iu], &ldwrku, &s[1], &rwork[ie], &
3938 work[itauq], &work[itaup], &work[iwork], &
3940 zlacpy_("L", m, m, &work[iu], &ldwrku, &u[u_offset],
3943 /* Generate right bidiagonalizing vectors in WORK(IU)
3944 (CWorkspace: need M*M+3*M, prefer M*M+2*M+(M-1)*NB)
3947 i__2 = *lwork - iwork + 1;
3948 zungbr_("P", m, m, m, &work[iu], &ldwrku, &work[itaup]
3949 , &work[iwork], &i__2, &ierr);
3951 /* Generate left bidiagonalizing vectors in U
3952 (CWorkspace: need M*M+3*M, prefer M*M+2*M+M*NB)
3955 i__2 = *lwork - iwork + 1;
3956 zungbr_("Q", m, m, m, &u[u_offset], ldu, &work[itauq],
3957 &work[iwork], &i__2, &ierr);
3960 /* Perform bidiagonal QR iteration, computing left
3961 singular vectors of L in U and computing right
3962 singular vectors of L in WORK(IU)
3963 (CWorkspace: need M*M)
3964 (RWorkspace: need BDSPAC) */
3966 zbdsqr_("U", m, m, m, &c__0, &s[1], &rwork[ie], &work[
3967 iu], &ldwrku, &u[u_offset], ldu, cdum, &c__1,
3968 &rwork[irwork], info);
3970 /* Multiply right singular vectors of L in WORK(IU) by
3971 Q in VT, storing result in A
3972 (CWorkspace: need M*M)
3975 zgemm_("N", "N", m, n, m, &c_b2, &work[iu], &ldwrku, &
3976 vt[vt_offset], ldvt, &c_b1, &a[a_offset], lda);
3978 /* Copy right singular vectors of A from A to VT */
3980 zlacpy_("F", m, n, &a[a_offset], lda, &vt[vt_offset],
3985 /* Insufficient workspace for a fast algorithm */
3990 /* Compute A=L*Q, copying result to VT
3991 (CWorkspace: need 2*M, prefer M+M*NB)
3994 i__2 = *lwork - iwork + 1;
3995 zgelqf_(m, n, &a[a_offset], lda, &work[itau], &work[
3996 iwork], &i__2, &ierr);
3997 zlacpy_("U", m, n, &a[a_offset], lda, &vt[vt_offset],
4001 (CWorkspace: need M+N, prefer M+N*NB)
4004 i__2 = *lwork - iwork + 1;
4005 zunglq_(n, n, m, &vt[vt_offset], ldvt, &work[itau], &
4006 work[iwork], &i__2, &ierr);
4008 /* Copy L to U, zeroing out above it */
4010 zlacpy_("L", m, m, &a[a_offset], lda, &u[u_offset],
4014 zlaset_("U", &i__2, &i__3, &c_b1, &c_b1, &u_ref(1, 2),
4021 /* Bidiagonalize L in U
4022 (CWorkspace: need 3*M, prefer 2*M+2*M*NB)
4023 (RWorkspace: need M) */
4025 i__2 = *lwork - iwork + 1;
4026 zgebrd_(m, m, &u[u_offset], ldu, &s[1], &rwork[ie], &
4027 work[itauq], &work[itaup], &work[iwork], &
4030 /* Multiply right bidiagonalizing vectors in U by Q
4032 (CWorkspace: need 2*M+N, prefer 2*M+N*NB)
4035 i__2 = *lwork - iwork + 1;
4036 zunmbr_("P", "L", "C", m, n, m, &u[u_offset], ldu, &
4037 work[itaup], &vt[vt_offset], ldvt, &work[
4038 iwork], &i__2, &ierr);
4040 /* Generate left bidiagonalizing vectors in U
4041 (CWorkspace: need 3*M, prefer 2*M+M*NB)
4044 i__2 = *lwork - iwork + 1;
4045 zungbr_("Q", m, m, m, &u[u_offset], ldu, &work[itauq],
4046 &work[iwork], &i__2, &ierr);
4049 /* Perform bidiagonal QR iteration, computing left
4050 singular vectors of A in U and computing right
4051 singular vectors of A in VT
4053 (RWorkspace: need BDSPAC) */
4055 zbdsqr_("U", m, n, m, &c__0, &s[1], &rwork[ie], &vt[
4056 vt_offset], ldvt, &u[u_offset], ldu, cdum, &
4057 c__1, &rwork[irwork], info);
4069 Path 10t(N greater than M, but not much larger)
4070 Reduce to bidiagonal form without LQ decomposition */
4078 (CWorkspace: need 2*M+N, prefer 2*M+(M+N)*NB)
4081 i__2 = *lwork - iwork + 1;
4082 zgebrd_(m, n, &a[a_offset], lda, &s[1], &rwork[ie], &work[itauq],
4083 &work[itaup], &work[iwork], &i__2, &ierr);
4086 /* If left singular vectors desired in U, copy result to U
4087 and generate left bidiagonalizing vectors in U
4088 (CWorkspace: need 3*M-1, prefer 2*M+(M-1)*NB)
4091 zlacpy_("L", m, m, &a[a_offset], lda, &u[u_offset], ldu);
4092 i__2 = *lwork - iwork + 1;
4093 zungbr_("Q", m, m, n, &u[u_offset], ldu, &work[itauq], &work[
4094 iwork], &i__2, &ierr);
4098 /* If right singular vectors desired in VT, copy result to
4099 VT and generate right bidiagonalizing vectors in VT
4100 (CWorkspace: need 2*M+NRVT, prefer 2*M+NRVT*NB)
4103 zlacpy_("U", m, n, &a[a_offset], lda, &vt[vt_offset], ldvt);
4110 i__2 = *lwork - iwork + 1;
4111 zungbr_("P", &nrvt, n, m, &vt[vt_offset], ldvt, &work[itaup],
4112 &work[iwork], &i__2, &ierr);
4116 /* If left singular vectors desired in A, generate left
4117 bidiagonalizing vectors in A
4118 (CWorkspace: need 3*M-1, prefer 2*M+(M-1)*NB)
4121 i__2 = *lwork - iwork + 1;
4122 zungbr_("Q", m, m, n, &a[a_offset], lda, &work[itauq], &work[
4123 iwork], &i__2, &ierr);
4127 /* If right singular vectors desired in A, generate right
4128 bidiagonalizing vectors in A
4129 (CWorkspace: need 3*M, prefer 2*M+M*NB)
4132 i__2 = *lwork - iwork + 1;
4133 zungbr_("P", m, n, m, &a[a_offset], lda, &work[itaup], &work[
4134 iwork], &i__2, &ierr);
4137 if (wntuas || wntuo) {
4143 if (wntvas || wntvo) {
4149 if (! wntuo && ! wntvo) {
4151 /* Perform bidiagonal QR iteration, if desired, computing
4152 left singular vectors in U and computing right singular
4155 (RWorkspace: need BDSPAC) */
4157 zbdsqr_("L", m, &ncvt, &nru, &c__0, &s[1], &rwork[ie], &vt[
4158 vt_offset], ldvt, &u[u_offset], ldu, cdum, &c__1, &
4159 rwork[irwork], info);
4160 } else if (! wntuo && wntvo) {
4162 /* Perform bidiagonal QR iteration, if desired, computing
4163 left singular vectors in U and computing right singular
4166 (RWorkspace: need BDSPAC) */
4168 zbdsqr_("L", m, &ncvt, &nru, &c__0, &s[1], &rwork[ie], &a[
4169 a_offset], lda, &u[u_offset], ldu, cdum, &c__1, &
4170 rwork[irwork], info);
4173 /* Perform bidiagonal QR iteration, if desired, computing
4174 left singular vectors in A and computing right singular
4177 (RWorkspace: need BDSPAC) */
4179 zbdsqr_("L", m, &ncvt, &nru, &c__0, &s[1], &rwork[ie], &vt[
4180 vt_offset], ldvt, &a[a_offset], lda, cdum, &c__1, &
4181 rwork[irwork], info);
4188 /* Undo scaling if necessary */
4191 if (anrm > bignum) {
4192 dlascl_("G", &c__0, &c__0, &bignum, &anrm, &minmn, &c__1, &s[1], &
4195 if (*info != 0 && anrm > bignum) {
4197 dlascl_("G", &c__0, &c__0, &bignum, &anrm, &i__2, &c__1, &rwork[
4198 ie], &minmn, &ierr);
4200 if (anrm < smlnum) {
4201 dlascl_("G", &c__0, &c__0, &smlnum, &anrm, &minmn, &c__1, &s[1], &
4204 if (*info != 0 && anrm < smlnum) {
4206 dlascl_("G", &c__0, &c__0, &smlnum, &anrm, &i__2, &c__1, &rwork[
4207 ie], &minmn, &ierr);
4211 /* Return optimal workspace in WORK(1) */
4213 work[1].r = (doublereal) maxwrk, work[1].i = 0.;