1 /* ../netlib/zgesvd.f -- translated by f2c (version 20100827). You must link the resulting object file with libf2c: on Microsoft Windows system, link with libf2c.lib;
2 on Linux or Unix systems, link with .../path/to/libf2c.a -lm or, if you install libf2c.a in a standard place, with -lf2c -lm -- in that order, at the end of the command line, as in cc *.o -lf2c -lm Source for libf2c is in /netlib/f2c/libf2c.zip, e.g., http://www.netlib.org/f2c/libf2c.zip */
3 #include "FLA_f2c.h" /* Table of constant values */
4 static doublecomplex c_b1 =
5 {
6 0.,0.
7 }
8 ;
9 static doublecomplex c_b2 =
10 {
11 1.,0.
12 }
13 ;
14 static integer c__6 = 6;
15 static integer c__0 = 0;
16 static integer c_n1 = -1;
17 static integer c__1 = 1;
18 /* > \brief <b> ZGESVD computes the singular value decomposition (SVD) for GE matrices</b> */
19 /* =========== DOCUMENTATION =========== */
20 /* Online html documentation available at */
21 /* http://www.netlib.org/lapack/explore-html/ */
22 /* > \htmlonly */
23 /* > Download ZGESVD + dependencies */
24 /* > <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zgesvd. f"> */
25 /* > [TGZ]</a> */
26 /* > <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zgesvd. f"> */
27 /* > [ZIP]</a> */
28 /* > <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zgesvd. f"> */
29 /* > [TXT]</a> */
30 /* > \endhtmlonly */
31 /* Definition: */
32 /* =========== */
33 /* SUBROUTINE ZGESVD( JOBU, JOBVT, M, N, A, LDA, S, U, LDU, VT, LDVT, */
34 /* WORK, LWORK, RWORK, INFO ) */
35 /* .. Scalar Arguments .. */
36 /* CHARACTER JOBU, JOBVT */
37 /* INTEGER INFO, LDA, LDU, LDVT, LWORK, M, N */
38 /* .. */
39 /* .. Array Arguments .. */
40 /* DOUBLE PRECISION RWORK( * ), S( * ) */
41 /* COMPLEX*16 A( LDA, * ), U( LDU, * ), VT( LDVT, * ), */
42 /* $ WORK( * ) */
43 /* .. */
44 /* > \par Purpose: */
45 /* ============= */
46 /* > */
47 /* > \verbatim */
48 /* > */
49 /* > ZGESVD computes the singular value decomposition (SVD) of a complex */
50 /* > M-by-N matrix A, optionally computing the left and/or right singular */
51 /* > vectors. The SVD is written */
52 /* > */
53 /* > A = U * SIGMA * conjugate-transpose(V) */
54 /* > */
55 /* > where SIGMA is an M-by-N matrix which is zero except for its */
56 /* > min(m,n) diagonal elements, U is an M-by-M unitary matrix, and */
57 /* > V is an N-by-N unitary matrix. The diagonal elements of SIGMA */
58 /* > are the singular values of A;
59 they are real and non-negative, and */
60 /* > are returned in descending order. The first min(m,n) columns of */
61 /* > U and V are the left and right singular vectors of A. */
62 /* > */
63 /* > Note that the routine returns V**H, not V. */
64 /* > \endverbatim */
65 /* Arguments: */
66 /* ========== */
67 /* > \param[in] JOBU */
68 /* > \verbatim */
69 /* > JOBU is CHARACTER*1 */
70 /* > Specifies options for computing all or part of the matrix U: */
71 /* > = 'A': all M columns of U are returned in array U: */
72 /* > = 'S': the first min(m,n) columns of U (the left singular */
73 /* > vectors) are returned in the array U;
74 */
75 /* > = 'O': the first min(m,n) columns of U (the left singular */
76 /* > vectors) are overwritten on the array A;
77 */
78 /* > = 'N': no columns of U (no left singular vectors) are */
79 /* > computed. */
80 /* > \endverbatim */
81 /* > */
82 /* > \param[in] JOBVT */
83 /* > \verbatim */
84 /* > JOBVT is CHARACTER*1 */
85 /* > Specifies options for computing all or part of the matrix */
86 /* > V**H: */
87 /* > = 'A': all N rows of V**H are returned in the array VT;
88 */
89 /* > = 'S': the first min(m,n) rows of V**H (the right singular */
90 /* > vectors) are returned in the array VT;
91 */
92 /* > = 'O': the first min(m,n) rows of V**H (the right singular */
93 /* > vectors) are overwritten on the array A;
94 */
95 /* > = 'N': no rows of V**H (no right singular vectors) are */
96 /* > computed. */
97 /* > */
98 /* > JOBVT and JOBU cannot both be 'O'. */
99 /* > \endverbatim */
100 /* > */
101 /* > \param[in] M */
102 /* > \verbatim */
103 /* > M is INTEGER */
104 /* > The number of rows of the input matrix A. M >= 0. */
105 /* > \endverbatim */
106 /* > */
107 /* > \param[in] N */
108 /* > \verbatim */
109 /* > N is INTEGER */
110 /* > The number of columns of the input matrix A. N >= 0. */
111 /* > \endverbatim */
112 /* > */
113 /* > \param[in,out] A */
114 /* > \verbatim */
115 /* > A is COMPLEX*16 array, dimension (LDA,N) */
116 /* > On entry, the M-by-N matrix A. */
117 /* > On exit, */
118 /* > if JOBU = 'O', A is overwritten with the first min(m,n) */
119 /* > columns of U (the left singular vectors, */
120 /* > stored columnwise);
121 */
122 /* > if JOBVT = 'O', A is overwritten with the first min(m,n) */
123 /* > rows of V**H (the right singular vectors, */
124 /* > stored rowwise);
125 */
126 /* > if JOBU .ne. 'O' and JOBVT .ne. 'O', the contents of A */
127 /* > are destroyed. */
128 /* > \endverbatim */
129 /* > */
130 /* > \param[in] LDA */
131 /* > \verbatim */
132 /* > LDA is INTEGER */
133 /* > The leading dimension of the array A. LDA >= max(1,M). */
134 /* > \endverbatim */
135 /* > */
136 /* > \param[out] S */
137 /* > \verbatim */
138 /* > S is DOUBLE PRECISION array, dimension (min(M,N)) */
139 /* > The singular values of A, sorted so that S(i) >= S(i+1). */
140 /* > \endverbatim */
141 /* > */
142 /* > \param[out] U */
143 /* > \verbatim */
144 /* > U is COMPLEX*16 array, dimension (LDU,UCOL) */
145 /* > (LDU,M) if JOBU = 'A' or (LDU,min(M,N)) if JOBU = 'S'. */
146 /* > If JOBU = 'A', U contains the M-by-M unitary matrix U;
147 */
148 /* > if JOBU = 'S', U contains the first min(m,n) columns of U */
149 /* > (the left singular vectors, stored columnwise);
150 */
151 /* > if JOBU = 'N' or 'O', U is not referenced. */
152 /* > \endverbatim */
153 /* > */
154 /* > \param[in] LDU */
155 /* > \verbatim */
156 /* > LDU is INTEGER */
157 /* > The leading dimension of the array U. LDU >= 1;
158 if */
159 /* > JOBU = 'S' or 'A', LDU >= M. */
160 /* > \endverbatim */
161 /* > */
162 /* > \param[out] VT */
163 /* > \verbatim */
164 /* > VT is COMPLEX*16 array, dimension (LDVT,N) */
165 /* > If JOBVT = 'A', VT contains the N-by-N unitary matrix */
166 /* > V**H;
167 */
168 /* > if JOBVT = 'S', VT contains the first min(m,n) rows of */
169 /* > V**H (the right singular vectors, stored rowwise);
170 */
171 /* > if JOBVT = 'N' or 'O', VT is not referenced. */
172 /* > \endverbatim */
173 /* > */
174 /* > \param[in] LDVT */
175 /* > \verbatim */
176 /* > LDVT is INTEGER */
177 /* > The leading dimension of the array VT. LDVT >= 1;
178 if */
179 /* > JOBVT = 'A', LDVT >= N;
180 if JOBVT = 'S', LDVT >= min(M,N). */
181 /* > \endverbatim */
182 /* > */
183 /* > \param[out] WORK */
184 /* > \verbatim */
185 /* > WORK is COMPLEX*16 array, dimension (MAX(1,LWORK)) */
186 /* > On exit, if INFO = 0, WORK(1) returns the optimal LWORK. */
187 /* > \endverbatim */
188 /* > */
189 /* > \param[in] LWORK */
190 /* > \verbatim */
191 /* > LWORK is INTEGER */
192 /* > The dimension of the array WORK. */
193 /* > LWORK >= MAX(1,2*MIN(M,N)+MAX(M,N)). */
194 /* > For good performance, LWORK should generally be larger. */
195 /* > */
196 /* > If LWORK = -1, then a workspace query is assumed;
197 the routine */
198 /* > only calculates the optimal size of the WORK array, returns */
199 /* > this value as the first entry of the WORK array, and no error */
200 /* > message related to LWORK is issued by XERBLA. */
201 /* > \endverbatim */
202 /* > */
203 /* > \param[out] RWORK */
204 /* > \verbatim */
205 /* > RWORK is DOUBLE PRECISION array, dimension (5*min(M,N)) */
206 /* > On exit, if INFO > 0, RWORK(1:MIN(M,N)-1) contains the */
207 /* > unconverged superdiagonal elements of an upper bidiagonal */
208 /* > matrix B whose diagonal is in S (not necessarily sorted). */
209 /* > B satisfies A = U * B * VT, so it has the same singular */
210 /* > values as A, and singular vectors related by U and VT. */
211 /* > \endverbatim */
212 /* > */
213 /* > \param[out] INFO */
214 /* > \verbatim */
215 /* > INFO is INTEGER */
216 /* > = 0: successful exit. */
217 /* > < 0: if INFO = -i, the i-th argument had an illegal value. */
218 /* > > 0: if ZBDSQR did not converge, INFO specifies how many */
219 /* > superdiagonals of an intermediate bidiagonal form B */
220 /* > did not converge to zero. See the description of RWORK */
221 /* > above for details. */
222 /* > \endverbatim */
223 /* Authors: */
224 /* ======== */
225 /* > \author Univ. of Tennessee */
226 /* > \author Univ. of California Berkeley */
227 /* > \author Univ. of Colorado Denver */
228 /* > \author NAG Ltd. */
229 /* > \date April 2012 */
230 /* > \ingroup complex16GEsing */
231 /* ===================================================================== */
232 /* Subroutine */
zgesvd_(char * jobu,char * jobvt,integer * m,integer * n,doublecomplex * a,integer * lda,doublereal * s,doublecomplex * u,integer * ldu,doublecomplex * vt,integer * ldvt,doublecomplex * work,integer * lwork,doublereal * rwork,integer * info)233 int zgesvd_(char *jobu, char *jobvt, integer *m, integer *n, doublecomplex *a, integer *lda, doublereal *s, doublecomplex *u, integer *ldu, doublecomplex *vt, integer *ldvt, doublecomplex *work, integer *lwork, doublereal *rwork, integer *info)
234 {
235 /* System generated locals */
236 integer a_dim1, a_offset, u_dim1, u_offset, vt_dim1, vt_offset, i__2, i__3, i__4;
237 char ch__1[2];
238 /* Builtin functions */
239 /* Subroutine */
240
241 double sqrt(doublereal);
242 /* Local variables */
243 integer i__, ie, ir, iu, blk, ncu;
244 doublereal dum[1], eps;
245 integer nru;
246 doublecomplex cdum[1];
247 integer iscl;
248 doublereal anrm;
249 integer ierr, itau, ncvt, nrvt, lwork_zgebrd__, lwork_zgelqf__, lwork_zgeqrf__;
250 extern logical lsame_(char *, char *);
251 integer chunk, minmn;
252 extern /* Subroutine */
253 int zgemm_(char *, char *, integer *, integer *, integer *, doublecomplex *, doublecomplex *, integer *, doublecomplex *, integer *, doublecomplex *, doublecomplex *, integer *);
254 integer wrkbl, itaup, itauq, mnthr, iwork;
255 logical wntua, wntva, wntun, wntuo, wntvn, wntvo, wntus, wntvs;
256 extern doublereal dlamch_(char *);
257 extern /* Subroutine */
258 int dlascl_(char *, integer *, integer *, doublereal *, doublereal *, integer *, integer *, doublereal *, integer *, integer *), xerbla_(char *, integer *), zgebrd_();
259 extern integer ilaenv_(integer *, char *, char *, integer *, integer *, integer *, integer *);
260 extern doublereal zlange_(char *, integer *, integer *, doublecomplex *, integer *, doublereal *);
261 doublereal bignum;
262 extern /* Subroutine */
263 int zgelqf_(), zlascl_(char *, integer *, integer *, doublereal *, doublereal *, integer *, integer *, doublecomplex *, integer *, integer *), zgeqrf_(), zlacpy_(char *, integer *, integer *, doublecomplex *, integer *, doublecomplex *, integer *), zlaset_(char *, integer *, integer *, doublecomplex *, doublecomplex *, doublecomplex *, integer *);
264 integer ldwrkr;
265 extern /* Subroutine */
266 int zbdsqr_(char *, integer *, integer *, integer *, integer *, doublereal *, doublereal *, doublecomplex *, integer *, doublecomplex *, integer *, doublecomplex *, integer *, doublereal *, integer *);
267 integer minwrk, ldwrku, maxwrk;
268 extern /* Subroutine */
269 int zungbr_();
270 doublereal smlnum;
271 integer irwork;
272 extern /* Subroutine */
273 int zunmbr_(char *, char *, char *, integer *, integer *, integer *, doublecomplex *, integer *, doublecomplex *, doublecomplex *, integer *, doublecomplex *, integer *, integer * ), zunglq_();
274 logical lquery, wntuas, wntvas;
275 extern /* Subroutine */
276 int zungqr_();
277 integer lwork_zungbr_p__, lwork_zungbr_q__, lwork_zunglq_m__, lwork_zunglq_n__, lwork_zungqr_m__, lwork_zungqr_n__;
278 /* -- LAPACK driver routine (version 3.4.1) -- */
279 /* -- LAPACK is a software package provided by Univ. of Tennessee, -- */
280 /* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- */
281 /* April 2012 */
282 /* .. Scalar Arguments .. */
283 /* .. */
284 /* .. Array Arguments .. */
285 /* .. */
286 /* ===================================================================== */
287 /* .. Parameters .. */
288 /* .. */
289 /* .. Local Scalars .. */
290 /* .. */
291 /* .. Local Arrays .. */
292 /* .. */
293 /* .. External Subroutines .. */
294 /* .. */
295 /* .. External Functions .. */
296 /* .. */
297 /* .. Intrinsic Functions .. */
298 /* .. */
299 /* .. Executable Statements .. */
300 /* Test the input arguments */
301 /* Parameter adjustments */
302 a_dim1 = *lda;
303 a_offset = 1 + a_dim1;
304 a -= a_offset;
305 --s;
306 u_dim1 = *ldu;
307 u_offset = 1 + u_dim1;
308 u -= u_offset;
309 vt_dim1 = *ldvt;
310 vt_offset = 1 + vt_dim1;
311 vt -= vt_offset;
312 --work;
313 --rwork;
314 /* Function Body */
315 *info = 0;
316 minmn = min(*m,*n);
317 wntua = lsame_(jobu, "A");
318 wntus = lsame_(jobu, "S");
319 wntuas = wntua || wntus;
320 wntuo = lsame_(jobu, "O");
321 wntun = lsame_(jobu, "N");
322 wntva = lsame_(jobvt, "A");
323 wntvs = lsame_(jobvt, "S");
324 wntvas = wntva || wntvs;
325 wntvo = lsame_(jobvt, "O");
326 wntvn = lsame_(jobvt, "N");
327 lquery = *lwork == -1;
328 if (! (wntua || wntus || wntuo || wntun))
329 {
330 *info = -1;
331 }
332 else if (! (wntva || wntvs || wntvo || wntvn) || wntvo && wntuo)
333 {
334 *info = -2;
335 }
336 else if (*m < 0)
337 {
338 *info = -3;
339 }
340 else if (*n < 0)
341 {
342 *info = -4;
343 }
344 else if (*lda < max(1,*m))
345 {
346 *info = -6;
347 }
348 else if (*ldu < 1 || wntuas && *ldu < *m)
349 {
350 *info = -9;
351 }
352 else if (*ldvt < 1 || wntva && *ldvt < *n || wntvs && *ldvt < minmn)
353 {
354 *info = -11;
355 }
356 /* Compute workspace */
357 /* (Note: Comments in the code beginning "Workspace:" describe the */
358 /* minimal amount of workspace needed at that point in the code, */
359 /* as well as the preferred amount for good performance. */
360 /* CWorkspace refers to complex workspace, and RWorkspace to */
361 /* real workspace. NB refers to the optimal block size for the */
362 /* immediately following subroutine, as returned by ILAENV.) */
363 if (*info == 0)
364 {
365 minwrk = 1;
366 maxwrk = 1;
367 if (*m >= *n && minmn > 0)
368 {
369 /* Space needed for ZBDSQR is BDSPAC = 5*N */
370 mnthr = ilaenv_(&c__6, "ZGESVD", ch__1, m, n, &c__0, &c__0);
371 /* Compute space needed for ZGEQRF */
372 zgeqrf_(m, n, &a[a_offset], lda, dum, dum, &c_n1, &ierr);
373 lwork_zgeqrf__ = (integer) dum[0];
374 /* Compute space needed for ZUNGQR */
375 zungqr_(m, n, n, &a[a_offset], lda, dum, dum, &c_n1, &ierr);
376 lwork_zungqr_n__ = (integer) dum[0];
377 zungqr_(m, m, n, &a[a_offset], lda, dum, dum, &c_n1, &ierr);
378 lwork_zungqr_m__ = (integer) dum[0];
379 /* Compute space needed for ZGEBRD */
380 zgebrd_(n, n, &a[a_offset], lda, &s[1], dum, dum, dum, dum, &c_n1, &ierr);
381 lwork_zgebrd__ = (integer) dum[0];
382 /* Compute space needed for ZUNGBR */
383 zungbr_("P", n, n, n, &a[a_offset], lda, dum, dum, &c_n1, &ierr);
384 lwork_zungbr_p__ = (integer) dum[0];
385 zungbr_("Q", n, n, n, &a[a_offset], lda, dum, dum, &c_n1, &ierr);
386 lwork_zungbr_q__ = (integer) dum[0];
387 if (*m >= mnthr)
388 {
389 if (wntun)
390 {
391 /* Path 1 (M much larger than N, JOBU='N') */
392 maxwrk = *n + lwork_zgeqrf__;
393 /* Computing MAX */
394 i__2 = maxwrk;
395 i__3 = (*n << 1) + lwork_zgebrd__; // , expr subst
396 maxwrk = max(i__2,i__3);
397 if (wntvo || wntvas)
398 {
399 /* Computing MAX */
400 i__2 = maxwrk;
401 i__3 = (*n << 1) + lwork_zungbr_p__; // , expr subst
402 maxwrk = max(i__2,i__3);
403 }
404 minwrk = *n * 3;
405 }
406 else if (wntuo && wntvn)
407 {
408 /* Path 2 (M much larger than N, JOBU='O', JOBVT='N') */
409 wrkbl = *n + lwork_zgeqrf__;
410 /* Computing MAX */
411 i__2 = wrkbl;
412 i__3 = *n + lwork_zungqr_n__; // , expr subst
413 wrkbl = max(i__2,i__3);
414 /* Computing MAX */
415 i__2 = wrkbl;
416 i__3 = (*n << 1) + lwork_zgebrd__; // , expr subst
417 wrkbl = max(i__2,i__3);
418 /* Computing MAX */
419 i__2 = wrkbl;
420 i__3 = (*n << 1) + lwork_zungbr_q__; // , expr subst
421 wrkbl = max(i__2,i__3);
422 /* Computing MAX */
423 i__2 = *n * *n + wrkbl;
424 i__3 = *n * *n + *m * *n; // , expr subst
425 maxwrk = max(i__2,i__3);
426 minwrk = (*n << 1) + *m;
427 }
428 else if (wntuo && wntvas)
429 {
430 /* Path 3 (M much larger than N, JOBU='O', JOBVT='S' or */
431 /* 'A') */
432 wrkbl = *n + lwork_zgeqrf__;
433 /* Computing MAX */
434 i__2 = wrkbl;
435 i__3 = *n + lwork_zungqr_n__; // , expr subst
436 wrkbl = max(i__2,i__3);
437 /* Computing MAX */
438 i__2 = wrkbl;
439 i__3 = (*n << 1) + lwork_zgebrd__; // , expr subst
440 wrkbl = max(i__2,i__3);
441 /* Computing MAX */
442 i__2 = wrkbl;
443 i__3 = (*n << 1) + lwork_zungbr_q__; // , expr subst
444 wrkbl = max(i__2,i__3);
445 /* Computing MAX */
446 i__2 = wrkbl;
447 i__3 = (*n << 1) + lwork_zungbr_p__; // , expr subst
448 wrkbl = max(i__2,i__3);
449 /* Computing MAX */
450 i__2 = *n * *n + wrkbl;
451 i__3 = *n * *n + *m * *n; // , expr subst
452 maxwrk = max(i__2,i__3);
453 minwrk = (*n << 1) + *m;
454 }
455 else if (wntus && wntvn)
456 {
457 /* Path 4 (M much larger than N, JOBU='S', JOBVT='N') */
458 wrkbl = *n + lwork_zgeqrf__;
459 /* Computing MAX */
460 i__2 = wrkbl;
461 i__3 = *n + lwork_zungqr_n__; // , expr subst
462 wrkbl = max(i__2,i__3);
463 /* Computing MAX */
464 i__2 = wrkbl;
465 i__3 = (*n << 1) + lwork_zgebrd__; // , expr subst
466 wrkbl = max(i__2,i__3);
467 /* Computing MAX */
468 i__2 = wrkbl;
469 i__3 = (*n << 1) + lwork_zungbr_q__; // , expr subst
470 wrkbl = max(i__2,i__3);
471 maxwrk = *n * *n + wrkbl;
472 minwrk = (*n << 1) + *m;
473 }
474 else if (wntus && wntvo)
475 {
476 /* Path 5 (M much larger than N, JOBU='S', JOBVT='O') */
477 wrkbl = *n + lwork_zgeqrf__;
478 /* Computing MAX */
479 i__2 = wrkbl;
480 i__3 = *n + lwork_zungqr_n__; // , expr subst
481 wrkbl = max(i__2,i__3);
482 /* Computing MAX */
483 i__2 = wrkbl;
484 i__3 = (*n << 1) + lwork_zgebrd__; // , expr subst
485 wrkbl = max(i__2,i__3);
486 /* Computing MAX */
487 i__2 = wrkbl;
488 i__3 = (*n << 1) + lwork_zungbr_q__; // , expr subst
489 wrkbl = max(i__2,i__3);
490 /* Computing MAX */
491 i__2 = wrkbl;
492 i__3 = (*n << 1) + lwork_zungbr_p__; // , expr subst
493 wrkbl = max(i__2,i__3);
494 maxwrk = (*n << 1) * *n + wrkbl;
495 minwrk = (*n << 1) + *m;
496 }
497 else if (wntus && wntvas)
498 {
499 /* Path 6 (M much larger than N, JOBU='S', JOBVT='S' or */
500 /* 'A') */
501 wrkbl = *n + lwork_zgeqrf__;
502 /* Computing MAX */
503 i__2 = wrkbl;
504 i__3 = *n + lwork_zungqr_n__; // , expr subst
505 wrkbl = max(i__2,i__3);
506 /* Computing MAX */
507 i__2 = wrkbl;
508 i__3 = (*n << 1) + lwork_zgebrd__; // , expr subst
509 wrkbl = max(i__2,i__3);
510 /* Computing MAX */
511 i__2 = wrkbl;
512 i__3 = (*n << 1) + lwork_zungbr_q__; // , expr subst
513 wrkbl = max(i__2,i__3);
514 /* Computing MAX */
515 i__2 = wrkbl;
516 i__3 = (*n << 1) + lwork_zungbr_p__; // , expr subst
517 wrkbl = max(i__2,i__3);
518 maxwrk = *n * *n + wrkbl;
519 minwrk = (*n << 1) + *m;
520 }
521 else if (wntua && wntvn)
522 {
523 /* Path 7 (M much larger than N, JOBU='A', JOBVT='N') */
524 wrkbl = *n + lwork_zgeqrf__;
525 /* Computing MAX */
526 i__2 = wrkbl;
527 i__3 = *n + lwork_zungqr_m__; // , expr subst
528 wrkbl = max(i__2,i__3);
529 /* Computing MAX */
530 i__2 = wrkbl;
531 i__3 = (*n << 1) + lwork_zgebrd__; // , expr subst
532 wrkbl = max(i__2,i__3);
533 /* Computing MAX */
534 i__2 = wrkbl;
535 i__3 = (*n << 1) + lwork_zungbr_q__; // , expr subst
536 wrkbl = max(i__2,i__3);
537 maxwrk = *n * *n + wrkbl;
538 minwrk = (*n << 1) + *m;
539 }
540 else if (wntua && wntvo)
541 {
542 /* Path 8 (M much larger than N, JOBU='A', JOBVT='O') */
543 wrkbl = *n + lwork_zgeqrf__;
544 /* Computing MAX */
545 i__2 = wrkbl;
546 i__3 = *n + lwork_zungqr_m__; // , expr subst
547 wrkbl = max(i__2,i__3);
548 /* Computing MAX */
549 i__2 = wrkbl;
550 i__3 = (*n << 1) + lwork_zgebrd__; // , expr subst
551 wrkbl = max(i__2,i__3);
552 /* Computing MAX */
553 i__2 = wrkbl;
554 i__3 = (*n << 1) + lwork_zungbr_q__; // , expr subst
555 wrkbl = max(i__2,i__3);
556 /* Computing MAX */
557 i__2 = wrkbl;
558 i__3 = (*n << 1) + lwork_zungbr_p__; // , expr subst
559 wrkbl = max(i__2,i__3);
560 maxwrk = (*n << 1) * *n + wrkbl;
561 minwrk = (*n << 1) + *m;
562 }
563 else if (wntua && wntvas)
564 {
565 /* Path 9 (M much larger than N, JOBU='A', JOBVT='S' or */
566 /* 'A') */
567 wrkbl = *n + lwork_zgeqrf__;
568 /* Computing MAX */
569 i__2 = wrkbl;
570 i__3 = *n + lwork_zungqr_m__; // , expr subst
571 wrkbl = max(i__2,i__3);
572 /* Computing MAX */
573 i__2 = wrkbl;
574 i__3 = (*n << 1) + lwork_zgebrd__; // , expr subst
575 wrkbl = max(i__2,i__3);
576 /* Computing MAX */
577 i__2 = wrkbl;
578 i__3 = (*n << 1) + lwork_zungbr_q__; // , expr subst
579 wrkbl = max(i__2,i__3);
580 /* Computing MAX */
581 i__2 = wrkbl;
582 i__3 = (*n << 1) + lwork_zungbr_p__; // , expr subst
583 wrkbl = max(i__2,i__3);
584 maxwrk = *n * *n + wrkbl;
585 minwrk = (*n << 1) + *m;
586 }
587 }
588 else
589 {
590 /* Path 10 (M at least N, but not much larger) */
591 zgebrd_(m, n, &a[a_offset], lda, &s[1], dum, dum, dum, dum, & c_n1, &ierr);
592 lwork_zgebrd__ = (integer) dum[0];
593 maxwrk = (*n << 1) + lwork_zgebrd__;
594 if (wntus || wntuo)
595 {
596 zungbr_("Q", m, n, n, &a[a_offset], lda, dum, dum, &c_n1, &ierr);
597 lwork_zungbr_q__ = (integer) dum[0];
598 /* Computing MAX */
599 i__2 = maxwrk;
600 i__3 = (*n << 1) + lwork_zungbr_q__; // , expr subst
601 maxwrk = max(i__2,i__3);
602 }
603 if (wntua)
604 {
605 zungbr_("Q", m, m, n, &a[a_offset], lda, dum, dum, &c_n1, &ierr);
606 lwork_zungbr_q__ = (integer) dum[0];
607 /* Computing MAX */
608 i__2 = maxwrk;
609 i__3 = (*n << 1) + lwork_zungbr_q__; // , expr subst
610 maxwrk = max(i__2,i__3);
611 }
612 if (! wntvn)
613 {
614 /* Computing MAX */
615 i__2 = maxwrk;
616 i__3 = (*n << 1) + lwork_zungbr_p__; // , expr subst
617 maxwrk = max(i__2,i__3);
618 minwrk = (*n << 1) + *m;
619 }
620 }
621 }
622 else if (minmn > 0)
623 {
624 /* Space needed for ZBDSQR is BDSPAC = 5*M */
625 mnthr = ilaenv_(&c__6, "ZGESVD", ch__1, m, n, &c__0, &c__0);
626 /* Compute space needed for ZGELQF */
627 zgelqf_(m, n, &a[a_offset], lda, dum, dum, &c_n1, &ierr);
628 lwork_zgelqf__ = (integer) dum[0];
629 /* Compute space needed for ZUNGLQ */
630 zunglq_(n, n, m, dum, n, dum, dum, &c_n1, &ierr);
631 lwork_zunglq_n__ = (integer) dum[0];
632 zunglq_(m, n, m, &a[a_offset], lda, dum, dum, &c_n1, &ierr);
633 lwork_zunglq_m__ = (integer) dum[0];
634 /* Compute space needed for ZGEBRD */
635 zgebrd_(m, m, &a[a_offset], lda, &s[1], dum, dum, dum, dum, &c_n1, &ierr);
636 lwork_zgebrd__ = (integer) dum[0];
637 /* Compute space needed for ZUNGBR P */
638 zungbr_("P", m, m, m, &a[a_offset], n, dum, dum, &c_n1, &ierr);
639 lwork_zungbr_p__ = (integer) dum[0];
640 /* Compute space needed for ZUNGBR Q */
641 zungbr_("Q", m, m, m, &a[a_offset], n, dum, dum, &c_n1, &ierr);
642 lwork_zungbr_q__ = (integer) dum[0];
643 if (*n >= mnthr)
644 {
645 if (wntvn)
646 {
647 /* Path 1t(N much larger than M, JOBVT='N') */
648 maxwrk = *m + lwork_zgelqf__;
649 /* Computing MAX */
650 i__2 = maxwrk;
651 i__3 = (*m << 1) + lwork_zgebrd__; // , expr subst
652 maxwrk = max(i__2,i__3);
653 if (wntuo || wntuas)
654 {
655 /* Computing MAX */
656 i__2 = maxwrk;
657 i__3 = (*m << 1) + lwork_zungbr_q__; // , expr subst
658 maxwrk = max(i__2,i__3);
659 }
660 minwrk = *m * 3;
661 }
662 else if (wntvo && wntun)
663 {
664 /* Path 2t(N much larger than M, JOBU='N', JOBVT='O') */
665 wrkbl = *m + lwork_zgelqf__;
666 /* Computing MAX */
667 i__2 = wrkbl;
668 i__3 = *m + lwork_zunglq_m__; // , expr subst
669 wrkbl = max(i__2,i__3);
670 /* Computing MAX */
671 i__2 = wrkbl;
672 i__3 = (*m << 1) + lwork_zgebrd__; // , expr subst
673 wrkbl = max(i__2,i__3);
674 /* Computing MAX */
675 i__2 = wrkbl;
676 i__3 = (*m << 1) + lwork_zungbr_p__; // , expr subst
677 wrkbl = max(i__2,i__3);
678 /* Computing MAX */
679 i__2 = *m * *m + wrkbl;
680 i__3 = *m * *m + *m * *n; // , expr subst
681 maxwrk = max(i__2,i__3);
682 minwrk = (*m << 1) + *n;
683 }
684 else if (wntvo && wntuas)
685 {
686 /* Path 3t(N much larger than M, JOBU='S' or 'A', */
687 /* JOBVT='O') */
688 wrkbl = *m + lwork_zgelqf__;
689 /* Computing MAX */
690 i__2 = wrkbl;
691 i__3 = *m + lwork_zunglq_m__; // , expr subst
692 wrkbl = max(i__2,i__3);
693 /* Computing MAX */
694 i__2 = wrkbl;
695 i__3 = (*m << 1) + lwork_zgebrd__; // , expr subst
696 wrkbl = max(i__2,i__3);
697 /* Computing MAX */
698 i__2 = wrkbl;
699 i__3 = (*m << 1) + lwork_zungbr_p__; // , expr subst
700 wrkbl = max(i__2,i__3);
701 /* Computing MAX */
702 i__2 = wrkbl;
703 i__3 = (*m << 1) + lwork_zungbr_q__; // , expr subst
704 wrkbl = max(i__2,i__3);
705 /* Computing MAX */
706 i__2 = *m * *m + wrkbl;
707 i__3 = *m * *m + *m * *n; // , expr subst
708 maxwrk = max(i__2,i__3);
709 minwrk = (*m << 1) + *n;
710 }
711 else if (wntvs && wntun)
712 {
713 /* Path 4t(N much larger than M, JOBU='N', JOBVT='S') */
714 wrkbl = *m + lwork_zgelqf__;
715 /* Computing MAX */
716 i__2 = wrkbl;
717 i__3 = *m + lwork_zunglq_m__; // , expr subst
718 wrkbl = max(i__2,i__3);
719 /* Computing MAX */
720 i__2 = wrkbl;
721 i__3 = (*m << 1) + lwork_zgebrd__; // , expr subst
722 wrkbl = max(i__2,i__3);
723 /* Computing MAX */
724 i__2 = wrkbl;
725 i__3 = (*m << 1) + lwork_zungbr_p__; // , expr subst
726 wrkbl = max(i__2,i__3);
727 maxwrk = *m * *m + wrkbl;
728 minwrk = (*m << 1) + *n;
729 }
730 else if (wntvs && wntuo)
731 {
732 /* Path 5t(N much larger than M, JOBU='O', JOBVT='S') */
733 wrkbl = *m + lwork_zgelqf__;
734 /* Computing MAX */
735 i__2 = wrkbl;
736 i__3 = *m + lwork_zunglq_m__; // , expr subst
737 wrkbl = max(i__2,i__3);
738 /* Computing MAX */
739 i__2 = wrkbl;
740 i__3 = (*m << 1) + lwork_zgebrd__; // , expr subst
741 wrkbl = max(i__2,i__3);
742 /* Computing MAX */
743 i__2 = wrkbl;
744 i__3 = (*m << 1) + lwork_zungbr_p__; // , expr subst
745 wrkbl = max(i__2,i__3);
746 /* Computing MAX */
747 i__2 = wrkbl;
748 i__3 = (*m << 1) + lwork_zungbr_q__; // , expr subst
749 wrkbl = max(i__2,i__3);
750 maxwrk = (*m << 1) * *m + wrkbl;
751 minwrk = (*m << 1) + *n;
752 }
753 else if (wntvs && wntuas)
754 {
755 /* Path 6t(N much larger than M, JOBU='S' or 'A', */
756 /* JOBVT='S') */
757 wrkbl = *m + lwork_zgelqf__;
758 /* Computing MAX */
759 i__2 = wrkbl;
760 i__3 = *m + lwork_zunglq_m__; // , expr subst
761 wrkbl = max(i__2,i__3);
762 /* Computing MAX */
763 i__2 = wrkbl;
764 i__3 = (*m << 1) + lwork_zgebrd__; // , expr subst
765 wrkbl = max(i__2,i__3);
766 /* Computing MAX */
767 i__2 = wrkbl;
768 i__3 = (*m << 1) + lwork_zungbr_p__; // , expr subst
769 wrkbl = max(i__2,i__3);
770 /* Computing MAX */
771 i__2 = wrkbl;
772 i__3 = (*m << 1) + lwork_zungbr_q__; // , expr subst
773 wrkbl = max(i__2,i__3);
774 maxwrk = *m * *m + wrkbl;
775 minwrk = (*m << 1) + *n;
776 }
777 else if (wntva && wntun)
778 {
779 /* Path 7t(N much larger than M, JOBU='N', JOBVT='A') */
780 wrkbl = *m + lwork_zgelqf__;
781 /* Computing MAX */
782 i__2 = wrkbl;
783 i__3 = *m + lwork_zunglq_n__; // , expr subst
784 wrkbl = max(i__2,i__3);
785 /* Computing MAX */
786 i__2 = wrkbl;
787 i__3 = (*m << 1) + lwork_zgebrd__; // , expr subst
788 wrkbl = max(i__2,i__3);
789 /* Computing MAX */
790 i__2 = wrkbl;
791 i__3 = (*m << 1) + lwork_zungbr_p__; // , expr subst
792 wrkbl = max(i__2,i__3);
793 maxwrk = *m * *m + wrkbl;
794 minwrk = (*m << 1) + *n;
795 }
796 else if (wntva && wntuo)
797 {
798 /* Path 8t(N much larger than M, JOBU='O', JOBVT='A') */
799 wrkbl = *m + lwork_zgelqf__;
800 /* Computing MAX */
801 i__2 = wrkbl;
802 i__3 = *m + lwork_zunglq_n__; // , expr subst
803 wrkbl = max(i__2,i__3);
804 /* Computing MAX */
805 i__2 = wrkbl;
806 i__3 = (*m << 1) + lwork_zgebrd__; // , expr subst
807 wrkbl = max(i__2,i__3);
808 /* Computing MAX */
809 i__2 = wrkbl;
810 i__3 = (*m << 1) + lwork_zungbr_p__; // , expr subst
811 wrkbl = max(i__2,i__3);
812 /* Computing MAX */
813 i__2 = wrkbl;
814 i__3 = (*m << 1) + lwork_zungbr_q__; // , expr subst
815 wrkbl = max(i__2,i__3);
816 maxwrk = (*m << 1) * *m + wrkbl;
817 minwrk = (*m << 1) + *n;
818 }
819 else if (wntva && wntuas)
820 {
821 /* Path 9t(N much larger than M, JOBU='S' or 'A', */
822 /* JOBVT='A') */
823 wrkbl = *m + lwork_zgelqf__;
824 /* Computing MAX */
825 i__2 = wrkbl;
826 i__3 = *m + lwork_zunglq_n__; // , expr subst
827 wrkbl = max(i__2,i__3);
828 /* Computing MAX */
829 i__2 = wrkbl;
830 i__3 = (*m << 1) + lwork_zgebrd__; // , expr subst
831 wrkbl = max(i__2,i__3);
832 /* Computing MAX */
833 i__2 = wrkbl;
834 i__3 = (*m << 1) + lwork_zungbr_p__; // , expr subst
835 wrkbl = max(i__2,i__3);
836 /* Computing MAX */
837 i__2 = wrkbl;
838 i__3 = (*m << 1) + lwork_zungbr_q__; // , expr subst
839 wrkbl = max(i__2,i__3);
840 maxwrk = *m * *m + wrkbl;
841 minwrk = (*m << 1) + *n;
842 }
843 }
844 else
845 {
846 /* Path 10t(N greater than M, but not much larger) */
847 zgebrd_(m, n, &a[a_offset], lda, &s[1], dum, dum, dum, dum, & c_n1, &ierr);
848 lwork_zgebrd__ = (integer) dum[0];
849 maxwrk = (*m << 1) + lwork_zgebrd__;
850 if (wntvs || wntvo)
851 {
852 /* Compute space needed for ZUNGBR P */
853 zungbr_("P", m, n, m, &a[a_offset], n, dum, dum, &c_n1, & ierr);
854 lwork_zungbr_p__ = (integer) dum[0];
855 /* Computing MAX */
856 i__2 = maxwrk;
857 i__3 = (*m << 1) + lwork_zungbr_p__; // , expr subst
858 maxwrk = max(i__2,i__3);
859 }
860 if (wntva)
861 {
862 zungbr_("P", n, n, m, &a[a_offset], n, dum, dum, &c_n1, & ierr);
863 lwork_zungbr_p__ = (integer) dum[0];
864 /* Computing MAX */
865 i__2 = maxwrk;
866 i__3 = (*m << 1) + lwork_zungbr_p__; // , expr subst
867 maxwrk = max(i__2,i__3);
868 }
869 if (! wntun)
870 {
871 /* Computing MAX */
872 i__2 = maxwrk;
873 i__3 = (*m << 1) + lwork_zungbr_q__; // , expr subst
874 maxwrk = max(i__2,i__3);
875 minwrk = (*m << 1) + *n;
876 }
877 }
878 }
879 maxwrk = max(maxwrk,minwrk);
880 work[1].r = (doublereal) maxwrk;
881 work[1].i = 0.; // , expr subst
882 if (*lwork < minwrk && ! lquery)
883 {
884 *info = -13;
885 }
886 }
887 if (*info != 0)
888 {
889 i__2 = -(*info);
890 xerbla_("ZGESVD", &i__2);
891 return 0;
892 }
893 else if (lquery)
894 {
895 return 0;
896 }
897 /* Quick return if possible */
898 if (*m == 0 || *n == 0)
899 {
900 return 0;
901 }
902 /* Get machine constants */
903 eps = dlamch_("P");
904 smlnum = sqrt(dlamch_("S")) / eps;
905 bignum = 1. / smlnum;
906 /* Scale A if max element outside range [SMLNUM,BIGNUM] */
907 anrm = zlange_("M", m, n, &a[a_offset], lda, dum);
908 iscl = 0;
909 if (anrm > 0. && anrm < smlnum)
910 {
911 iscl = 1;
912 zlascl_("G", &c__0, &c__0, &anrm, &smlnum, m, n, &a[a_offset], lda, & ierr);
913 }
914 else if (anrm > bignum)
915 {
916 iscl = 1;
917 zlascl_("G", &c__0, &c__0, &anrm, &bignum, m, n, &a[a_offset], lda, & ierr);
918 }
919 if (*m >= *n)
920 {
921 /* A has at least as many rows as columns. If A has sufficiently */
922 /* more rows than columns, first reduce using the QR */
923 /* decomposition (if sufficient workspace available) */
924 if (*m >= mnthr)
925 {
926 if (wntun)
927 {
928 /* Path 1 (M much larger than N, JOBU='N') */
929 /* No left singular vectors to be computed */
930 itau = 1;
931 iwork = itau + *n;
932 /* Compute A=Q*R */
933 /* (CWorkspace: need 2*N, prefer N+N*NB) */
934 /* (RWorkspace: need 0) */
935 i__2 = *lwork - iwork + 1;
936 zgeqrf_(m, n, &a[a_offset], lda, &work[itau], &work[iwork], & i__2, &ierr);
937 /* Zero out below R */
938 i__2 = *n - 1;
939 i__3 = *n - 1;
940 zlaset_("L", &i__2, &i__3, &c_b1, &c_b1, &a[a_dim1 + 2], lda);
941 ie = 1;
942 itauq = 1;
943 itaup = itauq + *n;
944 iwork = itaup + *n;
945 /* Bidiagonalize R in A */
946 /* (CWorkspace: need 3*N, prefer 2*N+2*N*NB) */
947 /* (RWorkspace: need N) */
948 i__2 = *lwork - iwork + 1;
949 zgebrd_(n, n, &a[a_offset], lda, &s[1], &rwork[ie], &work[ itauq], &work[itaup], &work[iwork], &i__2, &ierr);
950 ncvt = 0;
951 if (wntvo || wntvas)
952 {
953 /* If right singular vectors desired, generate P'. */
954 /* (CWorkspace: need 3*N-1, prefer 2*N+(N-1)*NB) */
955 /* (RWorkspace: 0) */
956 i__2 = *lwork - iwork + 1;
957 zungbr_("P", n, n, n, &a[a_offset], lda, &work[itaup], & work[iwork], &i__2, &ierr);
958 ncvt = *n;
959 }
960 irwork = ie + *n;
961 /* Perform bidiagonal QR iteration, computing right */
962 /* singular vectors of A in A if desired */
963 /* (CWorkspace: 0) */
964 /* (RWorkspace: need BDSPAC) */
965 zbdsqr_("U", n, &ncvt, &c__0, &c__0, &s[1], &rwork[ie], &a[ a_offset], lda, cdum, &c__1, cdum, &c__1, &rwork[ irwork], info);
966 /* If right singular vectors desired in VT, copy them there */
967 if (wntvas)
968 {
969 zlacpy_("F", n, n, &a[a_offset], lda, &vt[vt_offset], ldvt);
970 }
971 }
972 else if (wntuo && wntvn)
973 {
974 /* Path 2 (M much larger than N, JOBU='O', JOBVT='N') */
975 /* N left singular vectors to be overwritten on A and */
976 /* no right singular vectors to be computed */
977 if (*lwork >= *n * *n + *n * 3)
978 {
979 /* Sufficient workspace for a fast algorithm */
980 ir = 1;
981 /* Computing MAX */
982 i__2 = wrkbl;
983 i__3 = *lda * *n; // , expr subst
984 if (*lwork >= max(i__2,i__3) + *lda * *n)
985 {
986 /* WORK(IU) is LDA by N, WORK(IR) is LDA by N */
987 ldwrku = *lda;
988 ldwrkr = *lda;
989 }
990 else /* if(complicated condition) */
991 {
992 /* Computing MAX */
993 i__2 = wrkbl;
994 i__3 = *lda * *n; // , expr subst
995 if (*lwork >= max(i__2,i__3) + *n * *n)
996 {
997 /* WORK(IU) is LDA by N, WORK(IR) is N by N */
998 ldwrku = *lda;
999 ldwrkr = *n;
1000 }
1001 else
1002 {
1003 /* WORK(IU) is LDWRKU by N, WORK(IR) is N by N */
1004 ldwrku = (*lwork - *n * *n) / *n;
1005 ldwrkr = *n;
1006 }
1007 }
1008 itau = ir + ldwrkr * *n;
1009 iwork = itau + *n;
1010 /* Compute A=Q*R */
1011 /* (CWorkspace: need N*N+2*N, prefer N*N+N+N*NB) */
1012 /* (RWorkspace: 0) */
1013 i__2 = *lwork - iwork + 1;
1014 zgeqrf_(m, n, &a[a_offset], lda, &work[itau], &work[iwork] , &i__2, &ierr);
1015 /* Copy R to WORK(IR) and zero out below it */
1016 zlacpy_("U", n, n, &a[a_offset], lda, &work[ir], &ldwrkr);
1017 i__2 = *n - 1;
1018 i__3 = *n - 1;
1019 zlaset_("L", &i__2, &i__3, &c_b1, &c_b1, &work[ir + 1], & ldwrkr);
1020 /* Generate Q in A */
1021 /* (CWorkspace: need N*N+2*N, prefer N*N+N+N*NB) */
1022 /* (RWorkspace: 0) */
1023 i__2 = *lwork - iwork + 1;
1024 zungqr_(m, n, n, &a[a_offset], lda, &work[itau], &work[ iwork], &i__2, &ierr);
1025 ie = 1;
1026 itauq = itau;
1027 itaup = itauq + *n;
1028 iwork = itaup + *n;
1029 /* Bidiagonalize R in WORK(IR) */
1030 /* (CWorkspace: need N*N+3*N, prefer N*N+2*N+2*N*NB) */
1031 /* (RWorkspace: need N) */
1032 i__2 = *lwork - iwork + 1;
1033 zgebrd_(n, n, &work[ir], &ldwrkr, &s[1], &rwork[ie], & work[itauq], &work[itaup], &work[iwork], &i__2, & ierr);
1034 /* Generate left vectors bidiagonalizing R */
1035 /* (CWorkspace: need N*N+3*N, prefer N*N+2*N+N*NB) */
1036 /* (RWorkspace: need 0) */
1037 i__2 = *lwork - iwork + 1;
1038 zungbr_("Q", n, n, n, &work[ir], &ldwrkr, &work[itauq], & work[iwork], &i__2, &ierr);
1039 irwork = ie + *n;
1040 /* Perform bidiagonal QR iteration, computing left */
1041 /* singular vectors of R in WORK(IR) */
1042 /* (CWorkspace: need N*N) */
1043 /* (RWorkspace: need BDSPAC) */
1044 zbdsqr_("U", n, &c__0, n, &c__0, &s[1], &rwork[ie], cdum, &c__1, &work[ir], &ldwrkr, cdum, &c__1, &rwork[ irwork], info);
1045 iu = itauq;
1046 /* Multiply Q in A by left singular vectors of R in */
1047 /* WORK(IR), storing result in WORK(IU) and copying to A */
1048 /* (CWorkspace: need N*N+N, prefer N*N+M*N) */
1049 /* (RWorkspace: 0) */
1050 i__2 = *m;
1051 i__3 = ldwrku;
1052 for (i__ = 1;
1053 i__3 < 0 ? i__ >= i__2 : i__ <= i__2;
1054 i__ += i__3)
1055 {
1056 /* Computing MIN */
1057 i__4 = *m - i__ + 1;
1058 chunk = min(i__4,ldwrku);
1059 zgemm_("N", "N", &chunk, n, n, &c_b2, &a[i__ + a_dim1] , lda, &work[ir], &ldwrkr, &c_b1, &work[iu], & ldwrku);
1060 zlacpy_("F", &chunk, n, &work[iu], &ldwrku, &a[i__ + a_dim1], lda);
1061 /* L10: */
1062 }
1063 }
1064 else
1065 {
1066 /* Insufficient workspace for a fast algorithm */
1067 ie = 1;
1068 itauq = 1;
1069 itaup = itauq + *n;
1070 iwork = itaup + *n;
1071 /* Bidiagonalize A */
1072 /* (CWorkspace: need 2*N+M, prefer 2*N+(M+N)*NB) */
1073 /* (RWorkspace: N) */
1074 i__3 = *lwork - iwork + 1;
1075 zgebrd_(m, n, &a[a_offset], lda, &s[1], &rwork[ie], &work[ itauq], &work[itaup], &work[iwork], &i__3, &ierr);
1076 /* Generate left vectors bidiagonalizing A */
1077 /* (CWorkspace: need 3*N, prefer 2*N+N*NB) */
1078 /* (RWorkspace: 0) */
1079 i__3 = *lwork - iwork + 1;
1080 zungbr_("Q", m, n, n, &a[a_offset], lda, &work[itauq], & work[iwork], &i__3, &ierr);
1081 irwork = ie + *n;
1082 /* Perform bidiagonal QR iteration, computing left */
1083 /* singular vectors of A in A */
1084 /* (CWorkspace: need 0) */
1085 /* (RWorkspace: need BDSPAC) */
1086 zbdsqr_("U", n, &c__0, m, &c__0, &s[1], &rwork[ie], cdum, &c__1, &a[a_offset], lda, cdum, &c__1, &rwork[ irwork], info);
1087 }
1088 }
1089 else if (wntuo && wntvas)
1090 {
1091 /* Path 3 (M much larger than N, JOBU='O', JOBVT='S' or 'A') */
1092 /* N left singular vectors to be overwritten on A and */
1093 /* N right singular vectors to be computed in VT */
1094 if (*lwork >= *n * *n + *n * 3)
1095 {
1096 /* Sufficient workspace for a fast algorithm */
1097 ir = 1;
1098 /* Computing MAX */
1099 i__3 = wrkbl;
1100 i__2 = *lda * *n; // , expr subst
1101 if (*lwork >= max(i__3,i__2) + *lda * *n)
1102 {
1103 /* WORK(IU) is LDA by N and WORK(IR) is LDA by N */
1104 ldwrku = *lda;
1105 ldwrkr = *lda;
1106 }
1107 else /* if(complicated condition) */
1108 {
1109 /* Computing MAX */
1110 i__3 = wrkbl;
1111 i__2 = *lda * *n; // , expr subst
1112 if (*lwork >= max(i__3,i__2) + *n * *n)
1113 {
1114 /* WORK(IU) is LDA by N and WORK(IR) is N by N */
1115 ldwrku = *lda;
1116 ldwrkr = *n;
1117 }
1118 else
1119 {
1120 /* WORK(IU) is LDWRKU by N and WORK(IR) is N by N */
1121 ldwrku = (*lwork - *n * *n) / *n;
1122 ldwrkr = *n;
1123 }
1124 }
1125 itau = ir + ldwrkr * *n;
1126 iwork = itau + *n;
1127 /* Compute A=Q*R */
1128 /* (CWorkspace: need N*N+2*N, prefer N*N+N+N*NB) */
1129 /* (RWorkspace: 0) */
1130 i__3 = *lwork - iwork + 1;
1131 zgeqrf_(m, n, &a[a_offset], lda, &work[itau], &work[iwork] , &i__3, &ierr);
1132 /* Copy R to VT, zeroing out below it */
1133 zlacpy_("U", n, n, &a[a_offset], lda, &vt[vt_offset], ldvt);
1134 if (*n > 1)
1135 {
1136 i__3 = *n - 1;
1137 i__2 = *n - 1;
1138 zlaset_("L", &i__3, &i__2, &c_b1, &c_b1, &vt[vt_dim1 + 2], ldvt);
1139 }
1140 /* Generate Q in A */
1141 /* (CWorkspace: need N*N+2*N, prefer N*N+N+N*NB) */
1142 /* (RWorkspace: 0) */
1143 i__3 = *lwork - iwork + 1;
1144 zungqr_(m, n, n, &a[a_offset], lda, &work[itau], &work[ iwork], &i__3, &ierr);
1145 ie = 1;
1146 itauq = itau;
1147 itaup = itauq + *n;
1148 iwork = itaup + *n;
1149 /* Bidiagonalize R in VT, copying result to WORK(IR) */
1150 /* (CWorkspace: need N*N+3*N, prefer N*N+2*N+2*N*NB) */
1151 /* (RWorkspace: need N) */
1152 i__3 = *lwork - iwork + 1;
1153 zgebrd_(n, n, &vt[vt_offset], ldvt, &s[1], &rwork[ie], & work[itauq], &work[itaup], &work[iwork], &i__3, & ierr);
1154 zlacpy_("L", n, n, &vt[vt_offset], ldvt, &work[ir], & ldwrkr);
1155 /* Generate left vectors bidiagonalizing R in WORK(IR) */
1156 /* (CWorkspace: need N*N+3*N, prefer N*N+2*N+N*NB) */
1157 /* (RWorkspace: 0) */
1158 i__3 = *lwork - iwork + 1;
1159 zungbr_("Q", n, n, n, &work[ir], &ldwrkr, &work[itauq], & work[iwork], &i__3, &ierr);
1160 /* Generate right vectors bidiagonalizing R in VT */
1161 /* (CWorkspace: need N*N+3*N-1, prefer N*N+2*N+(N-1)*NB) */
1162 /* (RWorkspace: 0) */
1163 i__3 = *lwork - iwork + 1;
1164 zungbr_("P", n, n, n, &vt[vt_offset], ldvt, &work[itaup], &work[iwork], &i__3, &ierr);
1165 irwork = ie + *n;
1166 /* Perform bidiagonal QR iteration, computing left */
1167 /* singular vectors of R in WORK(IR) and computing right */
1168 /* singular vectors of R in VT */
1169 /* (CWorkspace: need N*N) */
1170 /* (RWorkspace: need BDSPAC) */
1171 zbdsqr_("U", n, n, n, &c__0, &s[1], &rwork[ie], &vt[ vt_offset], ldvt, &work[ir], &ldwrkr, cdum, &c__1, &rwork[irwork], info);
1172 iu = itauq;
1173 /* Multiply Q in A by left singular vectors of R in */
1174 /* WORK(IR), storing result in WORK(IU) and copying to A */
1175 /* (CWorkspace: need N*N+N, prefer N*N+M*N) */
1176 /* (RWorkspace: 0) */
1177 i__3 = *m;
1178 i__2 = ldwrku;
1179 for (i__ = 1;
1180 i__2 < 0 ? i__ >= i__3 : i__ <= i__3;
1181 i__ += i__2)
1182 {
1183 /* Computing MIN */
1184 i__4 = *m - i__ + 1;
1185 chunk = min(i__4,ldwrku);
1186 zgemm_("N", "N", &chunk, n, n, &c_b2, &a[i__ + a_dim1] , lda, &work[ir], &ldwrkr, &c_b1, &work[iu], & ldwrku);
1187 zlacpy_("F", &chunk, n, &work[iu], &ldwrku, &a[i__ + a_dim1], lda);
1188 /* L20: */
1189 }
1190 }
1191 else
1192 {
1193 /* Insufficient workspace for a fast algorithm */
1194 itau = 1;
1195 iwork = itau + *n;
1196 /* Compute A=Q*R */
1197 /* (CWorkspace: need 2*N, prefer N+N*NB) */
1198 /* (RWorkspace: 0) */
1199 i__2 = *lwork - iwork + 1;
1200 zgeqrf_(m, n, &a[a_offset], lda, &work[itau], &work[iwork] , &i__2, &ierr);
1201 /* Copy R to VT, zeroing out below it */
1202 zlacpy_("U", n, n, &a[a_offset], lda, &vt[vt_offset], ldvt);
1203 if (*n > 1)
1204 {
1205 i__2 = *n - 1;
1206 i__3 = *n - 1;
1207 zlaset_("L", &i__2, &i__3, &c_b1, &c_b1, &vt[vt_dim1 + 2], ldvt);
1208 }
1209 /* Generate Q in A */
1210 /* (CWorkspace: need 2*N, prefer N+N*NB) */
1211 /* (RWorkspace: 0) */
1212 i__2 = *lwork - iwork + 1;
1213 zungqr_(m, n, n, &a[a_offset], lda, &work[itau], &work[ iwork], &i__2, &ierr);
1214 ie = 1;
1215 itauq = itau;
1216 itaup = itauq + *n;
1217 iwork = itaup + *n;
1218 /* Bidiagonalize R in VT */
1219 /* (CWorkspace: need 3*N, prefer 2*N+2*N*NB) */
1220 /* (RWorkspace: N) */
1221 i__2 = *lwork - iwork + 1;
1222 zgebrd_(n, n, &vt[vt_offset], ldvt, &s[1], &rwork[ie], & work[itauq], &work[itaup], &work[iwork], &i__2, & ierr);
1223 /* Multiply Q in A by left vectors bidiagonalizing R */
1224 /* (CWorkspace: need 2*N+M, prefer 2*N+M*NB) */
1225 /* (RWorkspace: 0) */
1226 i__2 = *lwork - iwork + 1;
1227 zunmbr_("Q", "R", "N", m, n, n, &vt[vt_offset], ldvt, & work[itauq], &a[a_offset], lda, &work[iwork], & i__2, &ierr);
1228 /* Generate right vectors bidiagonalizing R in VT */
1229 /* (CWorkspace: need 3*N-1, prefer 2*N+(N-1)*NB) */
1230 /* (RWorkspace: 0) */
1231 i__2 = *lwork - iwork + 1;
1232 zungbr_("P", n, n, n, &vt[vt_offset], ldvt, &work[itaup], &work[iwork], &i__2, &ierr);
1233 irwork = ie + *n;
1234 /* Perform bidiagonal QR iteration, computing left */
1235 /* singular vectors of A in A and computing right */
1236 /* singular vectors of A in VT */
1237 /* (CWorkspace: 0) */
1238 /* (RWorkspace: need BDSPAC) */
1239 zbdsqr_("U", n, n, m, &c__0, &s[1], &rwork[ie], &vt[ vt_offset], ldvt, &a[a_offset], lda, cdum, &c__1, &rwork[irwork], info);
1240 }
1241 }
1242 else if (wntus)
1243 {
1244 if (wntvn)
1245 {
1246 /* Path 4 (M much larger than N, JOBU='S', JOBVT='N') */
1247 /* N left singular vectors to be computed in U and */
1248 /* no right singular vectors to be computed */
1249 if (*lwork >= *n * *n + *n * 3)
1250 {
1251 /* Sufficient workspace for a fast algorithm */
1252 ir = 1;
1253 if (*lwork >= wrkbl + *lda * *n)
1254 {
1255 /* WORK(IR) is LDA by N */
1256 ldwrkr = *lda;
1257 }
1258 else
1259 {
1260 /* WORK(IR) is N by N */
1261 ldwrkr = *n;
1262 }
1263 itau = ir + ldwrkr * *n;
1264 iwork = itau + *n;
1265 /* Compute A=Q*R */
1266 /* (CWorkspace: need N*N+2*N, prefer N*N+N+N*NB) */
1267 /* (RWorkspace: 0) */
1268 i__2 = *lwork - iwork + 1;
1269 zgeqrf_(m, n, &a[a_offset], lda, &work[itau], &work[ iwork], &i__2, &ierr);
1270 /* Copy R to WORK(IR), zeroing out below it */
1271 zlacpy_("U", n, n, &a[a_offset], lda, &work[ir], & ldwrkr);
1272 i__2 = *n - 1;
1273 i__3 = *n - 1;
1274 zlaset_("L", &i__2, &i__3, &c_b1, &c_b1, &work[ir + 1] , &ldwrkr);
1275 /* Generate Q in A */
1276 /* (CWorkspace: need N*N+2*N, prefer N*N+N+N*NB) */
1277 /* (RWorkspace: 0) */
1278 i__2 = *lwork - iwork + 1;
1279 zungqr_(m, n, n, &a[a_offset], lda, &work[itau], & work[iwork], &i__2, &ierr);
1280 ie = 1;
1281 itauq = itau;
1282 itaup = itauq + *n;
1283 iwork = itaup + *n;
1284 /* Bidiagonalize R in WORK(IR) */
1285 /* (CWorkspace: need N*N+3*N, prefer N*N+2*N+2*N*NB) */
1286 /* (RWorkspace: need N) */
1287 i__2 = *lwork - iwork + 1;
1288 zgebrd_(n, n, &work[ir], &ldwrkr, &s[1], &rwork[ie], & work[itauq], &work[itaup], &work[iwork], & i__2, &ierr);
1289 /* Generate left vectors bidiagonalizing R in WORK(IR) */
1290 /* (CWorkspace: need N*N+3*N, prefer N*N+2*N+N*NB) */
1291 /* (RWorkspace: 0) */
1292 i__2 = *lwork - iwork + 1;
1293 zungbr_("Q", n, n, n, &work[ir], &ldwrkr, &work[itauq] , &work[iwork], &i__2, &ierr);
1294 irwork = ie + *n;
1295 /* Perform bidiagonal QR iteration, computing left */
1296 /* singular vectors of R in WORK(IR) */
1297 /* (CWorkspace: need N*N) */
1298 /* (RWorkspace: need BDSPAC) */
1299 zbdsqr_("U", n, &c__0, n, &c__0, &s[1], &rwork[ie], cdum, &c__1, &work[ir], &ldwrkr, cdum, &c__1, &rwork[irwork], info);
1300 /* Multiply Q in A by left singular vectors of R in */
1301 /* WORK(IR), storing result in U */
1302 /* (CWorkspace: need N*N) */
1303 /* (RWorkspace: 0) */
1304 zgemm_("N", "N", m, n, n, &c_b2, &a[a_offset], lda, & work[ir], &ldwrkr, &c_b1, &u[u_offset], ldu);
1305 }
1306 else
1307 {
1308 /* Insufficient workspace for a fast algorithm */
1309 itau = 1;
1310 iwork = itau + *n;
1311 /* Compute A=Q*R, copying result to U */
1312 /* (CWorkspace: need 2*N, prefer N+N*NB) */
1313 /* (RWorkspace: 0) */
1314 i__2 = *lwork - iwork + 1;
1315 zgeqrf_(m, n, &a[a_offset], lda, &work[itau], &work[ iwork], &i__2, &ierr);
1316 zlacpy_("L", m, n, &a[a_offset], lda, &u[u_offset], ldu);
1317 /* Generate Q in U */
1318 /* (CWorkspace: need 2*N, prefer N+N*NB) */
1319 /* (RWorkspace: 0) */
1320 i__2 = *lwork - iwork + 1;
1321 zungqr_(m, n, n, &u[u_offset], ldu, &work[itau], & work[iwork], &i__2, &ierr);
1322 ie = 1;
1323 itauq = itau;
1324 itaup = itauq + *n;
1325 iwork = itaup + *n;
1326 /* Zero out below R in A */
1327 i__2 = *n - 1;
1328 i__3 = *n - 1;
1329 zlaset_("L", &i__2, &i__3, &c_b1, &c_b1, &a[a_dim1 + 2], lda);
1330 /* Bidiagonalize R in A */
1331 /* (CWorkspace: need 3*N, prefer 2*N+2*N*NB) */
1332 /* (RWorkspace: need N) */
1333 i__2 = *lwork - iwork + 1;
1334 zgebrd_(n, n, &a[a_offset], lda, &s[1], &rwork[ie], & work[itauq], &work[itaup], &work[iwork], & i__2, &ierr);
1335 /* Multiply Q in U by left vectors bidiagonalizing R */
1336 /* (CWorkspace: need 2*N+M, prefer 2*N+M*NB) */
1337 /* (RWorkspace: 0) */
1338 i__2 = *lwork - iwork + 1;
1339 zunmbr_("Q", "R", "N", m, n, n, &a[a_offset], lda, & work[itauq], &u[u_offset], ldu, &work[iwork], &i__2, &ierr) ;
1340 irwork = ie + *n;
1341 /* Perform bidiagonal QR iteration, computing left */
1342 /* singular vectors of A in U */
1343 /* (CWorkspace: 0) */
1344 /* (RWorkspace: need BDSPAC) */
1345 zbdsqr_("U", n, &c__0, m, &c__0, &s[1], &rwork[ie], cdum, &c__1, &u[u_offset], ldu, cdum, &c__1, & rwork[irwork], info);
1346 }
1347 }
1348 else if (wntvo)
1349 {
1350 /* Path 5 (M much larger than N, JOBU='S', JOBVT='O') */
1351 /* N left singular vectors to be computed in U and */
1352 /* N right singular vectors to be overwritten on A */
1353 if (*lwork >= (*n << 1) * *n + *n * 3)
1354 {
1355 /* Sufficient workspace for a fast algorithm */
1356 iu = 1;
1357 if (*lwork >= wrkbl + (*lda << 1) * *n)
1358 {
1359 /* WORK(IU) is LDA by N and WORK(IR) is LDA by N */
1360 ldwrku = *lda;
1361 ir = iu + ldwrku * *n;
1362 ldwrkr = *lda;
1363 }
1364 else if (*lwork >= wrkbl + (*lda + *n) * *n)
1365 {
1366 /* WORK(IU) is LDA by N and WORK(IR) is N by N */
1367 ldwrku = *lda;
1368 ir = iu + ldwrku * *n;
1369 ldwrkr = *n;
1370 }
1371 else
1372 {
1373 /* WORK(IU) is N by N and WORK(IR) is N by N */
1374 ldwrku = *n;
1375 ir = iu + ldwrku * *n;
1376 ldwrkr = *n;
1377 }
1378 itau = ir + ldwrkr * *n;
1379 iwork = itau + *n;
1380 /* Compute A=Q*R */
1381 /* (CWorkspace: need 2*N*N+2*N, prefer 2*N*N+N+N*NB) */
1382 /* (RWorkspace: 0) */
1383 i__2 = *lwork - iwork + 1;
1384 zgeqrf_(m, n, &a[a_offset], lda, &work[itau], &work[ iwork], &i__2, &ierr);
1385 /* Copy R to WORK(IU), zeroing out below it */
1386 zlacpy_("U", n, n, &a[a_offset], lda, &work[iu], & ldwrku);
1387 i__2 = *n - 1;
1388 i__3 = *n - 1;
1389 zlaset_("L", &i__2, &i__3, &c_b1, &c_b1, &work[iu + 1] , &ldwrku);
1390 /* Generate Q in A */
1391 /* (CWorkspace: need 2*N*N+2*N, prefer 2*N*N+N+N*NB) */
1392 /* (RWorkspace: 0) */
1393 i__2 = *lwork - iwork + 1;
1394 zungqr_(m, n, n, &a[a_offset], lda, &work[itau], & work[iwork], &i__2, &ierr);
1395 ie = 1;
1396 itauq = itau;
1397 itaup = itauq + *n;
1398 iwork = itaup + *n;
1399 /* Bidiagonalize R in WORK(IU), copying result to */
1400 /* WORK(IR) */
1401 /* (CWorkspace: need 2*N*N+3*N, */
1402 /* prefer 2*N*N+2*N+2*N*NB) */
1403 /* (RWorkspace: need N) */
1404 i__2 = *lwork - iwork + 1;
1405 zgebrd_(n, n, &work[iu], &ldwrku, &s[1], &rwork[ie], & work[itauq], &work[itaup], &work[iwork], & i__2, &ierr);
1406 zlacpy_("U", n, n, &work[iu], &ldwrku, &work[ir], & ldwrkr);
1407 /* Generate left bidiagonalizing vectors in WORK(IU) */
1408 /* (CWorkspace: need 2*N*N+3*N, prefer 2*N*N+2*N+N*NB) */
1409 /* (RWorkspace: 0) */
1410 i__2 = *lwork - iwork + 1;
1411 zungbr_("Q", n, n, n, &work[iu], &ldwrku, &work[itauq] , &work[iwork], &i__2, &ierr);
1412 /* Generate right bidiagonalizing vectors in WORK(IR) */
1413 /* (CWorkspace: need 2*N*N+3*N-1, */
1414 /* prefer 2*N*N+2*N+(N-1)*NB) */
1415 /* (RWorkspace: 0) */
1416 i__2 = *lwork - iwork + 1;
1417 zungbr_("P", n, n, n, &work[ir], &ldwrkr, &work[itaup] , &work[iwork], &i__2, &ierr);
1418 irwork = ie + *n;
1419 /* Perform bidiagonal QR iteration, computing left */
1420 /* singular vectors of R in WORK(IU) and computing */
1421 /* right singular vectors of R in WORK(IR) */
1422 /* (CWorkspace: need 2*N*N) */
1423 /* (RWorkspace: need BDSPAC) */
1424 zbdsqr_("U", n, n, n, &c__0, &s[1], &rwork[ie], &work[ ir], &ldwrkr, &work[iu], &ldwrku, cdum, &c__1, &rwork[irwork], info);
1425 /* Multiply Q in A by left singular vectors of R in */
1426 /* WORK(IU), storing result in U */
1427 /* (CWorkspace: need N*N) */
1428 /* (RWorkspace: 0) */
1429 zgemm_("N", "N", m, n, n, &c_b2, &a[a_offset], lda, & work[iu], &ldwrku, &c_b1, &u[u_offset], ldu);
1430 /* Copy right singular vectors of R to A */
1431 /* (CWorkspace: need N*N) */
1432 /* (RWorkspace: 0) */
1433 zlacpy_("F", n, n, &work[ir], &ldwrkr, &a[a_offset], lda);
1434 }
1435 else
1436 {
1437 /* Insufficient workspace for a fast algorithm */
1438 itau = 1;
1439 iwork = itau + *n;
1440 /* Compute A=Q*R, copying result to U */
1441 /* (CWorkspace: need 2*N, prefer N+N*NB) */
1442 /* (RWorkspace: 0) */
1443 i__2 = *lwork - iwork + 1;
1444 zgeqrf_(m, n, &a[a_offset], lda, &work[itau], &work[ iwork], &i__2, &ierr);
1445 zlacpy_("L", m, n, &a[a_offset], lda, &u[u_offset], ldu);
1446 /* Generate Q in U */
1447 /* (CWorkspace: need 2*N, prefer N+N*NB) */
1448 /* (RWorkspace: 0) */
1449 i__2 = *lwork - iwork + 1;
1450 zungqr_(m, n, n, &u[u_offset], ldu, &work[itau], & work[iwork], &i__2, &ierr);
1451 ie = 1;
1452 itauq = itau;
1453 itaup = itauq + *n;
1454 iwork = itaup + *n;
1455 /* Zero out below R in A */
1456 i__2 = *n - 1;
1457 i__3 = *n - 1;
1458 zlaset_("L", &i__2, &i__3, &c_b1, &c_b1, &a[a_dim1 + 2], lda);
1459 /* Bidiagonalize R in A */
1460 /* (CWorkspace: need 3*N, prefer 2*N+2*N*NB) */
1461 /* (RWorkspace: need N) */
1462 i__2 = *lwork - iwork + 1;
1463 zgebrd_(n, n, &a[a_offset], lda, &s[1], &rwork[ie], & work[itauq], &work[itaup], &work[iwork], & i__2, &ierr);
1464 /* Multiply Q in U by left vectors bidiagonalizing R */
1465 /* (CWorkspace: need 2*N+M, prefer 2*N+M*NB) */
1466 /* (RWorkspace: 0) */
1467 i__2 = *lwork - iwork + 1;
1468 zunmbr_("Q", "R", "N", m, n, n, &a[a_offset], lda, & work[itauq], &u[u_offset], ldu, &work[iwork], &i__2, &ierr) ;
1469 /* Generate right vectors bidiagonalizing R in A */
1470 /* (CWorkspace: need 3*N-1, prefer 2*N+(N-1)*NB) */
1471 /* (RWorkspace: 0) */
1472 i__2 = *lwork - iwork + 1;
1473 zungbr_("P", n, n, n, &a[a_offset], lda, &work[itaup], &work[iwork], &i__2, &ierr);
1474 irwork = ie + *n;
1475 /* Perform bidiagonal QR iteration, computing left */
1476 /* singular vectors of A in U and computing right */
1477 /* singular vectors of A in A */
1478 /* (CWorkspace: 0) */
1479 /* (RWorkspace: need BDSPAC) */
1480 zbdsqr_("U", n, n, m, &c__0, &s[1], &rwork[ie], &a[ a_offset], lda, &u[u_offset], ldu, cdum, & c__1, &rwork[irwork], info);
1481 }
1482 }
1483 else if (wntvas)
1484 {
1485 /* Path 6 (M much larger than N, JOBU='S', JOBVT='S' */
1486 /* or 'A') */
1487 /* N left singular vectors to be computed in U and */
1488 /* N right singular vectors to be computed in VT */
1489 if (*lwork >= *n * *n + *n * 3)
1490 {
1491 /* Sufficient workspace for a fast algorithm */
1492 iu = 1;
1493 if (*lwork >= wrkbl + *lda * *n)
1494 {
1495 /* WORK(IU) is LDA by N */
1496 ldwrku = *lda;
1497 }
1498 else
1499 {
1500 /* WORK(IU) is N by N */
1501 ldwrku = *n;
1502 }
1503 itau = iu + ldwrku * *n;
1504 iwork = itau + *n;
1505 /* Compute A=Q*R */
1506 /* (CWorkspace: need N*N+2*N, prefer N*N+N+N*NB) */
1507 /* (RWorkspace: 0) */
1508 i__2 = *lwork - iwork + 1;
1509 zgeqrf_(m, n, &a[a_offset], lda, &work[itau], &work[ iwork], &i__2, &ierr);
1510 /* Copy R to WORK(IU), zeroing out below it */
1511 zlacpy_("U", n, n, &a[a_offset], lda, &work[iu], & ldwrku);
1512 i__2 = *n - 1;
1513 i__3 = *n - 1;
1514 zlaset_("L", &i__2, &i__3, &c_b1, &c_b1, &work[iu + 1] , &ldwrku);
1515 /* Generate Q in A */
1516 /* (CWorkspace: need N*N+2*N, prefer N*N+N+N*NB) */
1517 /* (RWorkspace: 0) */
1518 i__2 = *lwork - iwork + 1;
1519 zungqr_(m, n, n, &a[a_offset], lda, &work[itau], & work[iwork], &i__2, &ierr);
1520 ie = 1;
1521 itauq = itau;
1522 itaup = itauq + *n;
1523 iwork = itaup + *n;
1524 /* Bidiagonalize R in WORK(IU), copying result to VT */
1525 /* (CWorkspace: need N*N+3*N, prefer N*N+2*N+2*N*NB) */
1526 /* (RWorkspace: need N) */
1527 i__2 = *lwork - iwork + 1;
1528 zgebrd_(n, n, &work[iu], &ldwrku, &s[1], &rwork[ie], & work[itauq], &work[itaup], &work[iwork], & i__2, &ierr);
1529 zlacpy_("U", n, n, &work[iu], &ldwrku, &vt[vt_offset], ldvt);
1530 /* Generate left bidiagonalizing vectors in WORK(IU) */
1531 /* (CWorkspace: need N*N+3*N, prefer N*N+2*N+N*NB) */
1532 /* (RWorkspace: 0) */
1533 i__2 = *lwork - iwork + 1;
1534 zungbr_("Q", n, n, n, &work[iu], &ldwrku, &work[itauq] , &work[iwork], &i__2, &ierr);
1535 /* Generate right bidiagonalizing vectors in VT */
1536 /* (CWorkspace: need N*N+3*N-1, */
1537 /* prefer N*N+2*N+(N-1)*NB) */
1538 /* (RWorkspace: 0) */
1539 i__2 = *lwork - iwork + 1;
1540 zungbr_("P", n, n, n, &vt[vt_offset], ldvt, &work[ itaup], &work[iwork], &i__2, &ierr) ;
1541 irwork = ie + *n;
1542 /* Perform bidiagonal QR iteration, computing left */
1543 /* singular vectors of R in WORK(IU) and computing */
1544 /* right singular vectors of R in VT */
1545 /* (CWorkspace: need N*N) */
1546 /* (RWorkspace: need BDSPAC) */
1547 zbdsqr_("U", n, n, n, &c__0, &s[1], &rwork[ie], &vt[ vt_offset], ldvt, &work[iu], &ldwrku, cdum, & c__1, &rwork[irwork], info);
1548 /* Multiply Q in A by left singular vectors of R in */
1549 /* WORK(IU), storing result in U */
1550 /* (CWorkspace: need N*N) */
1551 /* (RWorkspace: 0) */
1552 zgemm_("N", "N", m, n, n, &c_b2, &a[a_offset], lda, & work[iu], &ldwrku, &c_b1, &u[u_offset], ldu);
1553 }
1554 else
1555 {
1556 /* Insufficient workspace for a fast algorithm */
1557 itau = 1;
1558 iwork = itau + *n;
1559 /* Compute A=Q*R, copying result to U */
1560 /* (CWorkspace: need 2*N, prefer N+N*NB) */
1561 /* (RWorkspace: 0) */
1562 i__2 = *lwork - iwork + 1;
1563 zgeqrf_(m, n, &a[a_offset], lda, &work[itau], &work[ iwork], &i__2, &ierr);
1564 zlacpy_("L", m, n, &a[a_offset], lda, &u[u_offset], ldu);
1565 /* Generate Q in U */
1566 /* (CWorkspace: need 2*N, prefer N+N*NB) */
1567 /* (RWorkspace: 0) */
1568 i__2 = *lwork - iwork + 1;
1569 zungqr_(m, n, n, &u[u_offset], ldu, &work[itau], & work[iwork], &i__2, &ierr);
1570 /* Copy R to VT, zeroing out below it */
1571 zlacpy_("U", n, n, &a[a_offset], lda, &vt[vt_offset], ldvt);
1572 if (*n > 1)
1573 {
1574 i__2 = *n - 1;
1575 i__3 = *n - 1;
1576 zlaset_("L", &i__2, &i__3, &c_b1, &c_b1, &vt[ vt_dim1 + 2], ldvt);
1577 }
1578 ie = 1;
1579 itauq = itau;
1580 itaup = itauq + *n;
1581 iwork = itaup + *n;
1582 /* Bidiagonalize R in VT */
1583 /* (CWorkspace: need 3*N, prefer 2*N+2*N*NB) */
1584 /* (RWorkspace: need N) */
1585 i__2 = *lwork - iwork + 1;
1586 zgebrd_(n, n, &vt[vt_offset], ldvt, &s[1], &rwork[ie], &work[itauq], &work[itaup], &work[iwork], & i__2, &ierr);
1587 /* Multiply Q in U by left bidiagonalizing vectors */
1588 /* in VT */
1589 /* (CWorkspace: need 2*N+M, prefer 2*N+M*NB) */
1590 /* (RWorkspace: 0) */
1591 i__2 = *lwork - iwork + 1;
1592 zunmbr_("Q", "R", "N", m, n, n, &vt[vt_offset], ldvt, &work[itauq], &u[u_offset], ldu, &work[iwork], &i__2, &ierr);
1593 /* Generate right bidiagonalizing vectors in VT */
1594 /* (CWorkspace: need 3*N-1, prefer 2*N+(N-1)*NB) */
1595 /* (RWorkspace: 0) */
1596 i__2 = *lwork - iwork + 1;
1597 zungbr_("P", n, n, n, &vt[vt_offset], ldvt, &work[ itaup], &work[iwork], &i__2, &ierr) ;
1598 irwork = ie + *n;
1599 /* Perform bidiagonal QR iteration, computing left */
1600 /* singular vectors of A in U and computing right */
1601 /* singular vectors of A in VT */
1602 /* (CWorkspace: 0) */
1603 /* (RWorkspace: need BDSPAC) */
1604 zbdsqr_("U", n, n, m, &c__0, &s[1], &rwork[ie], &vt[ vt_offset], ldvt, &u[u_offset], ldu, cdum, & c__1, &rwork[irwork], info);
1605 }
1606 }
1607 }
1608 else if (wntua)
1609 {
1610 if (wntvn)
1611 {
1612 /* Path 7 (M much larger than N, JOBU='A', JOBVT='N') */
1613 /* M left singular vectors to be computed in U and */
1614 /* no right singular vectors to be computed */
1615 /* Computing MAX */
1616 i__2 = *n + *m;
1617 i__3 = *n * 3; // , expr subst
1618 if (*lwork >= *n * *n + max(i__2,i__3))
1619 {
1620 /* Sufficient workspace for a fast algorithm */
1621 ir = 1;
1622 if (*lwork >= wrkbl + *lda * *n)
1623 {
1624 /* WORK(IR) is LDA by N */
1625 ldwrkr = *lda;
1626 }
1627 else
1628 {
1629 /* WORK(IR) is N by N */
1630 ldwrkr = *n;
1631 }
1632 itau = ir + ldwrkr * *n;
1633 iwork = itau + *n;
1634 /* Compute A=Q*R, copying result to U */
1635 /* (CWorkspace: need N*N+2*N, prefer N*N+N+N*NB) */
1636 /* (RWorkspace: 0) */
1637 i__2 = *lwork - iwork + 1;
1638 zgeqrf_(m, n, &a[a_offset], lda, &work[itau], &work[ iwork], &i__2, &ierr);
1639 zlacpy_("L", m, n, &a[a_offset], lda, &u[u_offset], ldu);
1640 /* Copy R to WORK(IR), zeroing out below it */
1641 zlacpy_("U", n, n, &a[a_offset], lda, &work[ir], & ldwrkr);
1642 i__2 = *n - 1;
1643 i__3 = *n - 1;
1644 zlaset_("L", &i__2, &i__3, &c_b1, &c_b1, &work[ir + 1] , &ldwrkr);
1645 /* Generate Q in U */
1646 /* (CWorkspace: need N*N+N+M, prefer N*N+N+M*NB) */
1647 /* (RWorkspace: 0) */
1648 i__2 = *lwork - iwork + 1;
1649 zungqr_(m, m, n, &u[u_offset], ldu, &work[itau], & work[iwork], &i__2, &ierr);
1650 ie = 1;
1651 itauq = itau;
1652 itaup = itauq + *n;
1653 iwork = itaup + *n;
1654 /* Bidiagonalize R in WORK(IR) */
1655 /* (CWorkspace: need N*N+3*N, prefer N*N+2*N+2*N*NB) */
1656 /* (RWorkspace: need N) */
1657 i__2 = *lwork - iwork + 1;
1658 zgebrd_(n, n, &work[ir], &ldwrkr, &s[1], &rwork[ie], & work[itauq], &work[itaup], &work[iwork], & i__2, &ierr);
1659 /* Generate left bidiagonalizing vectors in WORK(IR) */
1660 /* (CWorkspace: need N*N+3*N, prefer N*N+2*N+N*NB) */
1661 /* (RWorkspace: 0) */
1662 i__2 = *lwork - iwork + 1;
1663 zungbr_("Q", n, n, n, &work[ir], &ldwrkr, &work[itauq] , &work[iwork], &i__2, &ierr);
1664 irwork = ie + *n;
1665 /* Perform bidiagonal QR iteration, computing left */
1666 /* singular vectors of R in WORK(IR) */
1667 /* (CWorkspace: need N*N) */
1668 /* (RWorkspace: need BDSPAC) */
1669 zbdsqr_("U", n, &c__0, n, &c__0, &s[1], &rwork[ie], cdum, &c__1, &work[ir], &ldwrkr, cdum, &c__1, &rwork[irwork], info);
1670 /* Multiply Q in U by left singular vectors of R in */
1671 /* WORK(IR), storing result in A */
1672 /* (CWorkspace: need N*N) */
1673 /* (RWorkspace: 0) */
1674 zgemm_("N", "N", m, n, n, &c_b2, &u[u_offset], ldu, & work[ir], &ldwrkr, &c_b1, &a[a_offset], lda);
1675 /* Copy left singular vectors of A from A to U */
1676 zlacpy_("F", m, n, &a[a_offset], lda, &u[u_offset], ldu);
1677 }
1678 else
1679 {
1680 /* Insufficient workspace for a fast algorithm */
1681 itau = 1;
1682 iwork = itau + *n;
1683 /* Compute A=Q*R, copying result to U */
1684 /* (CWorkspace: need 2*N, prefer N+N*NB) */
1685 /* (RWorkspace: 0) */
1686 i__2 = *lwork - iwork + 1;
1687 zgeqrf_(m, n, &a[a_offset], lda, &work[itau], &work[ iwork], &i__2, &ierr);
1688 zlacpy_("L", m, n, &a[a_offset], lda, &u[u_offset], ldu);
1689 /* Generate Q in U */
1690 /* (CWorkspace: need N+M, prefer N+M*NB) */
1691 /* (RWorkspace: 0) */
1692 i__2 = *lwork - iwork + 1;
1693 zungqr_(m, m, n, &u[u_offset], ldu, &work[itau], & work[iwork], &i__2, &ierr);
1694 ie = 1;
1695 itauq = itau;
1696 itaup = itauq + *n;
1697 iwork = itaup + *n;
1698 /* Zero out below R in A */
1699 i__2 = *n - 1;
1700 i__3 = *n - 1;
1701 zlaset_("L", &i__2, &i__3, &c_b1, &c_b1, &a[a_dim1 + 2], lda);
1702 /* Bidiagonalize R in A */
1703 /* (CWorkspace: need 3*N, prefer 2*N+2*N*NB) */
1704 /* (RWorkspace: need N) */
1705 i__2 = *lwork - iwork + 1;
1706 zgebrd_(n, n, &a[a_offset], lda, &s[1], &rwork[ie], & work[itauq], &work[itaup], &work[iwork], & i__2, &ierr);
1707 /* Multiply Q in U by left bidiagonalizing vectors */
1708 /* in A */
1709 /* (CWorkspace: need 2*N+M, prefer 2*N+M*NB) */
1710 /* (RWorkspace: 0) */
1711 i__2 = *lwork - iwork + 1;
1712 zunmbr_("Q", "R", "N", m, n, n, &a[a_offset], lda, & work[itauq], &u[u_offset], ldu, &work[iwork], &i__2, &ierr) ;
1713 irwork = ie + *n;
1714 /* Perform bidiagonal QR iteration, computing left */
1715 /* singular vectors of A in U */
1716 /* (CWorkspace: 0) */
1717 /* (RWorkspace: need BDSPAC) */
1718 zbdsqr_("U", n, &c__0, m, &c__0, &s[1], &rwork[ie], cdum, &c__1, &u[u_offset], ldu, cdum, &c__1, & rwork[irwork], info);
1719 }
1720 }
1721 else if (wntvo)
1722 {
1723 /* Path 8 (M much larger than N, JOBU='A', JOBVT='O') */
1724 /* M left singular vectors to be computed in U and */
1725 /* N right singular vectors to be overwritten on A */
1726 /* Computing MAX */
1727 i__2 = *n + *m;
1728 i__3 = *n * 3; // , expr subst
1729 if (*lwork >= (*n << 1) * *n + max(i__2,i__3))
1730 {
1731 /* Sufficient workspace for a fast algorithm */
1732 iu = 1;
1733 if (*lwork >= wrkbl + (*lda << 1) * *n)
1734 {
1735 /* WORK(IU) is LDA by N and WORK(IR) is LDA by N */
1736 ldwrku = *lda;
1737 ir = iu + ldwrku * *n;
1738 ldwrkr = *lda;
1739 }
1740 else if (*lwork >= wrkbl + (*lda + *n) * *n)
1741 {
1742 /* WORK(IU) is LDA by N and WORK(IR) is N by N */
1743 ldwrku = *lda;
1744 ir = iu + ldwrku * *n;
1745 ldwrkr = *n;
1746 }
1747 else
1748 {
1749 /* WORK(IU) is N by N and WORK(IR) is N by N */
1750 ldwrku = *n;
1751 ir = iu + ldwrku * *n;
1752 ldwrkr = *n;
1753 }
1754 itau = ir + ldwrkr * *n;
1755 iwork = itau + *n;
1756 /* Compute A=Q*R, copying result to U */
1757 /* (CWorkspace: need 2*N*N+2*N, prefer 2*N*N+N+N*NB) */
1758 /* (RWorkspace: 0) */
1759 i__2 = *lwork - iwork + 1;
1760 zgeqrf_(m, n, &a[a_offset], lda, &work[itau], &work[ iwork], &i__2, &ierr);
1761 zlacpy_("L", m, n, &a[a_offset], lda, &u[u_offset], ldu);
1762 /* Generate Q in U */
1763 /* (CWorkspace: need 2*N*N+N+M, prefer 2*N*N+N+M*NB) */
1764 /* (RWorkspace: 0) */
1765 i__2 = *lwork - iwork + 1;
1766 zungqr_(m, m, n, &u[u_offset], ldu, &work[itau], & work[iwork], &i__2, &ierr);
1767 /* Copy R to WORK(IU), zeroing out below it */
1768 zlacpy_("U", n, n, &a[a_offset], lda, &work[iu], & ldwrku);
1769 i__2 = *n - 1;
1770 i__3 = *n - 1;
1771 zlaset_("L", &i__2, &i__3, &c_b1, &c_b1, &work[iu + 1] , &ldwrku);
1772 ie = 1;
1773 itauq = itau;
1774 itaup = itauq + *n;
1775 iwork = itaup + *n;
1776 /* Bidiagonalize R in WORK(IU), copying result to */
1777 /* WORK(IR) */
1778 /* (CWorkspace: need 2*N*N+3*N, */
1779 /* prefer 2*N*N+2*N+2*N*NB) */
1780 /* (RWorkspace: need N) */
1781 i__2 = *lwork - iwork + 1;
1782 zgebrd_(n, n, &work[iu], &ldwrku, &s[1], &rwork[ie], & work[itauq], &work[itaup], &work[iwork], & i__2, &ierr);
1783 zlacpy_("U", n, n, &work[iu], &ldwrku, &work[ir], & ldwrkr);
1784 /* Generate left bidiagonalizing vectors in WORK(IU) */
1785 /* (CWorkspace: need 2*N*N+3*N, prefer 2*N*N+2*N+N*NB) */
1786 /* (RWorkspace: 0) */
1787 i__2 = *lwork - iwork + 1;
1788 zungbr_("Q", n, n, n, &work[iu], &ldwrku, &work[itauq] , &work[iwork], &i__2, &ierr);
1789 /* Generate right bidiagonalizing vectors in WORK(IR) */
1790 /* (CWorkspace: need 2*N*N+3*N-1, */
1791 /* prefer 2*N*N+2*N+(N-1)*NB) */
1792 /* (RWorkspace: 0) */
1793 i__2 = *lwork - iwork + 1;
1794 zungbr_("P", n, n, n, &work[ir], &ldwrkr, &work[itaup] , &work[iwork], &i__2, &ierr);
1795 irwork = ie + *n;
1796 /* Perform bidiagonal QR iteration, computing left */
1797 /* singular vectors of R in WORK(IU) and computing */
1798 /* right singular vectors of R in WORK(IR) */
1799 /* (CWorkspace: need 2*N*N) */
1800 /* (RWorkspace: need BDSPAC) */
1801 zbdsqr_("U", n, n, n, &c__0, &s[1], &rwork[ie], &work[ ir], &ldwrkr, &work[iu], &ldwrku, cdum, &c__1, &rwork[irwork], info);
1802 /* Multiply Q in U by left singular vectors of R in */
1803 /* WORK(IU), storing result in A */
1804 /* (CWorkspace: need N*N) */
1805 /* (RWorkspace: 0) */
1806 zgemm_("N", "N", m, n, n, &c_b2, &u[u_offset], ldu, & work[iu], &ldwrku, &c_b1, &a[a_offset], lda);
1807 /* Copy left singular vectors of A from A to U */
1808 zlacpy_("F", m, n, &a[a_offset], lda, &u[u_offset], ldu);
1809 /* Copy right singular vectors of R from WORK(IR) to A */
1810 zlacpy_("F", n, n, &work[ir], &ldwrkr, &a[a_offset], lda);
1811 }
1812 else
1813 {
1814 /* Insufficient workspace for a fast algorithm */
1815 itau = 1;
1816 iwork = itau + *n;
1817 /* Compute A=Q*R, copying result to U */
1818 /* (CWorkspace: need 2*N, prefer N+N*NB) */
1819 /* (RWorkspace: 0) */
1820 i__2 = *lwork - iwork + 1;
1821 zgeqrf_(m, n, &a[a_offset], lda, &work[itau], &work[ iwork], &i__2, &ierr);
1822 zlacpy_("L", m, n, &a[a_offset], lda, &u[u_offset], ldu);
1823 /* Generate Q in U */
1824 /* (CWorkspace: need N+M, prefer N+M*NB) */
1825 /* (RWorkspace: 0) */
1826 i__2 = *lwork - iwork + 1;
1827 zungqr_(m, m, n, &u[u_offset], ldu, &work[itau], & work[iwork], &i__2, &ierr);
1828 ie = 1;
1829 itauq = itau;
1830 itaup = itauq + *n;
1831 iwork = itaup + *n;
1832 /* Zero out below R in A */
1833 i__2 = *n - 1;
1834 i__3 = *n - 1;
1835 zlaset_("L", &i__2, &i__3, &c_b1, &c_b1, &a[a_dim1 + 2], lda);
1836 /* Bidiagonalize R in A */
1837 /* (CWorkspace: need 3*N, prefer 2*N+2*N*NB) */
1838 /* (RWorkspace: need N) */
1839 i__2 = *lwork - iwork + 1;
1840 zgebrd_(n, n, &a[a_offset], lda, &s[1], &rwork[ie], & work[itauq], &work[itaup], &work[iwork], & i__2, &ierr);
1841 /* Multiply Q in U by left bidiagonalizing vectors */
1842 /* in A */
1843 /* (CWorkspace: need 2*N+M, prefer 2*N+M*NB) */
1844 /* (RWorkspace: 0) */
1845 i__2 = *lwork - iwork + 1;
1846 zunmbr_("Q", "R", "N", m, n, n, &a[a_offset], lda, & work[itauq], &u[u_offset], ldu, &work[iwork], &i__2, &ierr) ;
1847 /* Generate right bidiagonalizing vectors in A */
1848 /* (CWorkspace: need 3*N-1, prefer 2*N+(N-1)*NB) */
1849 /* (RWorkspace: 0) */
1850 i__2 = *lwork - iwork + 1;
1851 zungbr_("P", n, n, n, &a[a_offset], lda, &work[itaup], &work[iwork], &i__2, &ierr);
1852 irwork = ie + *n;
1853 /* Perform bidiagonal QR iteration, computing left */
1854 /* singular vectors of A in U and computing right */
1855 /* singular vectors of A in A */
1856 /* (CWorkspace: 0) */
1857 /* (RWorkspace: need BDSPAC) */
1858 zbdsqr_("U", n, n, m, &c__0, &s[1], &rwork[ie], &a[ a_offset], lda, &u[u_offset], ldu, cdum, & c__1, &rwork[irwork], info);
1859 }
1860 }
1861 else if (wntvas)
1862 {
1863 /* Path 9 (M much larger than N, JOBU='A', JOBVT='S' */
1864 /* or 'A') */
1865 /* M left singular vectors to be computed in U and */
1866 /* N right singular vectors to be computed in VT */
1867 /* Computing MAX */
1868 i__2 = *n + *m;
1869 i__3 = *n * 3; // , expr subst
1870 if (*lwork >= *n * *n + max(i__2,i__3))
1871 {
1872 /* Sufficient workspace for a fast algorithm */
1873 iu = 1;
1874 if (*lwork >= wrkbl + *lda * *n)
1875 {
1876 /* WORK(IU) is LDA by N */
1877 ldwrku = *lda;
1878 }
1879 else
1880 {
1881 /* WORK(IU) is N by N */
1882 ldwrku = *n;
1883 }
1884 itau = iu + ldwrku * *n;
1885 iwork = itau + *n;
1886 /* Compute A=Q*R, copying result to U */
1887 /* (CWorkspace: need N*N+2*N, prefer N*N+N+N*NB) */
1888 /* (RWorkspace: 0) */
1889 i__2 = *lwork - iwork + 1;
1890 zgeqrf_(m, n, &a[a_offset], lda, &work[itau], &work[ iwork], &i__2, &ierr);
1891 zlacpy_("L", m, n, &a[a_offset], lda, &u[u_offset], ldu);
1892 /* Generate Q in U */
1893 /* (CWorkspace: need N*N+N+M, prefer N*N+N+M*NB) */
1894 /* (RWorkspace: 0) */
1895 i__2 = *lwork - iwork + 1;
1896 zungqr_(m, m, n, &u[u_offset], ldu, &work[itau], & work[iwork], &i__2, &ierr);
1897 /* Copy R to WORK(IU), zeroing out below it */
1898 zlacpy_("U", n, n, &a[a_offset], lda, &work[iu], & ldwrku);
1899 i__2 = *n - 1;
1900 i__3 = *n - 1;
1901 zlaset_("L", &i__2, &i__3, &c_b1, &c_b1, &work[iu + 1] , &ldwrku);
1902 ie = 1;
1903 itauq = itau;
1904 itaup = itauq + *n;
1905 iwork = itaup + *n;
1906 /* Bidiagonalize R in WORK(IU), copying result to VT */
1907 /* (CWorkspace: need N*N+3*N, prefer N*N+2*N+2*N*NB) */
1908 /* (RWorkspace: need N) */
1909 i__2 = *lwork - iwork + 1;
1910 zgebrd_(n, n, &work[iu], &ldwrku, &s[1], &rwork[ie], & work[itauq], &work[itaup], &work[iwork], & i__2, &ierr);
1911 zlacpy_("U", n, n, &work[iu], &ldwrku, &vt[vt_offset], ldvt);
1912 /* Generate left bidiagonalizing vectors in WORK(IU) */
1913 /* (CWorkspace: need N*N+3*N, prefer N*N+2*N+N*NB) */
1914 /* (RWorkspace: 0) */
1915 i__2 = *lwork - iwork + 1;
1916 zungbr_("Q", n, n, n, &work[iu], &ldwrku, &work[itauq] , &work[iwork], &i__2, &ierr);
1917 /* Generate right bidiagonalizing vectors in VT */
1918 /* (CWorkspace: need N*N+3*N-1, */
1919 /* prefer N*N+2*N+(N-1)*NB) */
1920 /* (RWorkspace: need 0) */
1921 i__2 = *lwork - iwork + 1;
1922 zungbr_("P", n, n, n, &vt[vt_offset], ldvt, &work[ itaup], &work[iwork], &i__2, &ierr) ;
1923 irwork = ie + *n;
1924 /* Perform bidiagonal QR iteration, computing left */
1925 /* singular vectors of R in WORK(IU) and computing */
1926 /* right singular vectors of R in VT */
1927 /* (CWorkspace: need N*N) */
1928 /* (RWorkspace: need BDSPAC) */
1929 zbdsqr_("U", n, n, n, &c__0, &s[1], &rwork[ie], &vt[ vt_offset], ldvt, &work[iu], &ldwrku, cdum, & c__1, &rwork[irwork], info);
1930 /* Multiply Q in U by left singular vectors of R in */
1931 /* WORK(IU), storing result in A */
1932 /* (CWorkspace: need N*N) */
1933 /* (RWorkspace: 0) */
1934 zgemm_("N", "N", m, n, n, &c_b2, &u[u_offset], ldu, & work[iu], &ldwrku, &c_b1, &a[a_offset], lda);
1935 /* Copy left singular vectors of A from A to U */
1936 zlacpy_("F", m, n, &a[a_offset], lda, &u[u_offset], ldu);
1937 }
1938 else
1939 {
1940 /* Insufficient workspace for a fast algorithm */
1941 itau = 1;
1942 iwork = itau + *n;
1943 /* Compute A=Q*R, copying result to U */
1944 /* (CWorkspace: need 2*N, prefer N+N*NB) */
1945 /* (RWorkspace: 0) */
1946 i__2 = *lwork - iwork + 1;
1947 zgeqrf_(m, n, &a[a_offset], lda, &work[itau], &work[ iwork], &i__2, &ierr);
1948 zlacpy_("L", m, n, &a[a_offset], lda, &u[u_offset], ldu);
1949 /* Generate Q in U */
1950 /* (CWorkspace: need N+M, prefer N+M*NB) */
1951 /* (RWorkspace: 0) */
1952 i__2 = *lwork - iwork + 1;
1953 zungqr_(m, m, n, &u[u_offset], ldu, &work[itau], & work[iwork], &i__2, &ierr);
1954 /* Copy R from A to VT, zeroing out below it */
1955 zlacpy_("U", n, n, &a[a_offset], lda, &vt[vt_offset], ldvt);
1956 if (*n > 1)
1957 {
1958 i__2 = *n - 1;
1959 i__3 = *n - 1;
1960 zlaset_("L", &i__2, &i__3, &c_b1, &c_b1, &vt[ vt_dim1 + 2], ldvt);
1961 }
1962 ie = 1;
1963 itauq = itau;
1964 itaup = itauq + *n;
1965 iwork = itaup + *n;
1966 /* Bidiagonalize R in VT */
1967 /* (CWorkspace: need 3*N, prefer 2*N+2*N*NB) */
1968 /* (RWorkspace: need N) */
1969 i__2 = *lwork - iwork + 1;
1970 zgebrd_(n, n, &vt[vt_offset], ldvt, &s[1], &rwork[ie], &work[itauq], &work[itaup], &work[iwork], & i__2, &ierr);
1971 /* Multiply Q in U by left bidiagonalizing vectors */
1972 /* in VT */
1973 /* (CWorkspace: need 2*N+M, prefer 2*N+M*NB) */
1974 /* (RWorkspace: 0) */
1975 i__2 = *lwork - iwork + 1;
1976 zunmbr_("Q", "R", "N", m, n, n, &vt[vt_offset], ldvt, &work[itauq], &u[u_offset], ldu, &work[iwork], &i__2, &ierr);
1977 /* Generate right bidiagonalizing vectors in VT */
1978 /* (CWorkspace: need 3*N-1, prefer 2*N+(N-1)*NB) */
1979 /* (RWorkspace: 0) */
1980 i__2 = *lwork - iwork + 1;
1981 zungbr_("P", n, n, n, &vt[vt_offset], ldvt, &work[ itaup], &work[iwork], &i__2, &ierr) ;
1982 irwork = ie + *n;
1983 /* Perform bidiagonal QR iteration, computing left */
1984 /* singular vectors of A in U and computing right */
1985 /* singular vectors of A in VT */
1986 /* (CWorkspace: 0) */
1987 /* (RWorkspace: need BDSPAC) */
1988 zbdsqr_("U", n, n, m, &c__0, &s[1], &rwork[ie], &vt[ vt_offset], ldvt, &u[u_offset], ldu, cdum, & c__1, &rwork[irwork], info);
1989 }
1990 }
1991 }
1992 }
1993 else
1994 {
1995 /* M .LT. MNTHR */
1996 /* Path 10 (M at least N, but not much larger) */
1997 /* Reduce to bidiagonal form without QR decomposition */
1998 ie = 1;
1999 itauq = 1;
2000 itaup = itauq + *n;
2001 iwork = itaup + *n;
2002 /* Bidiagonalize A */
2003 /* (CWorkspace: need 2*N+M, prefer 2*N+(M+N)*NB) */
2004 /* (RWorkspace: need N) */
2005 i__2 = *lwork - iwork + 1;
2006 zgebrd_(m, n, &a[a_offset], lda, &s[1], &rwork[ie], &work[itauq], &work[itaup], &work[iwork], &i__2, &ierr);
2007 if (wntuas)
2008 {
2009 /* If left singular vectors desired in U, copy result to U */
2010 /* and generate left bidiagonalizing vectors in U */
2011 /* (CWorkspace: need 2*N+NCU, prefer 2*N+NCU*NB) */
2012 /* (RWorkspace: 0) */
2013 zlacpy_("L", m, n, &a[a_offset], lda, &u[u_offset], ldu);
2014 if (wntus)
2015 {
2016 ncu = *n;
2017 }
2018 if (wntua)
2019 {
2020 ncu = *m;
2021 }
2022 i__2 = *lwork - iwork + 1;
2023 zungbr_("Q", m, &ncu, n, &u[u_offset], ldu, &work[itauq], & work[iwork], &i__2, &ierr);
2024 }
2025 if (wntvas)
2026 {
2027 /* If right singular vectors desired in VT, copy result to */
2028 /* VT and generate right bidiagonalizing vectors in VT */
2029 /* (CWorkspace: need 3*N-1, prefer 2*N+(N-1)*NB) */
2030 /* (RWorkspace: 0) */
2031 zlacpy_("U", n, n, &a[a_offset], lda, &vt[vt_offset], ldvt);
2032 i__2 = *lwork - iwork + 1;
2033 zungbr_("P", n, n, n, &vt[vt_offset], ldvt, &work[itaup], & work[iwork], &i__2, &ierr);
2034 }
2035 if (wntuo)
2036 {
2037 /* If left singular vectors desired in A, generate left */
2038 /* bidiagonalizing vectors in A */
2039 /* (CWorkspace: need 3*N, prefer 2*N+N*NB) */
2040 /* (RWorkspace: 0) */
2041 i__2 = *lwork - iwork + 1;
2042 zungbr_("Q", m, n, n, &a[a_offset], lda, &work[itauq], &work[ iwork], &i__2, &ierr);
2043 }
2044 if (wntvo)
2045 {
2046 /* If right singular vectors desired in A, generate right */
2047 /* bidiagonalizing vectors in A */
2048 /* (CWorkspace: need 3*N-1, prefer 2*N+(N-1)*NB) */
2049 /* (RWorkspace: 0) */
2050 i__2 = *lwork - iwork + 1;
2051 zungbr_("P", n, n, n, &a[a_offset], lda, &work[itaup], &work[ iwork], &i__2, &ierr);
2052 }
2053 irwork = ie + *n;
2054 if (wntuas || wntuo)
2055 {
2056 nru = *m;
2057 }
2058 if (wntun)
2059 {
2060 nru = 0;
2061 }
2062 if (wntvas || wntvo)
2063 {
2064 ncvt = *n;
2065 }
2066 if (wntvn)
2067 {
2068 ncvt = 0;
2069 }
2070 if (! wntuo && ! wntvo)
2071 {
2072 /* Perform bidiagonal QR iteration, if desired, computing */
2073 /* left singular vectors in U and computing right singular */
2074 /* vectors in VT */
2075 /* (CWorkspace: 0) */
2076 /* (RWorkspace: need BDSPAC) */
2077 zbdsqr_("U", n, &ncvt, &nru, &c__0, &s[1], &rwork[ie], &vt[ vt_offset], ldvt, &u[u_offset], ldu, cdum, &c__1, & rwork[irwork], info);
2078 }
2079 else if (! wntuo && wntvo)
2080 {
2081 /* Perform bidiagonal QR iteration, if desired, computing */
2082 /* left singular vectors in U and computing right singular */
2083 /* vectors in A */
2084 /* (CWorkspace: 0) */
2085 /* (RWorkspace: need BDSPAC) */
2086 zbdsqr_("U", n, &ncvt, &nru, &c__0, &s[1], &rwork[ie], &a[ a_offset], lda, &u[u_offset], ldu, cdum, &c__1, & rwork[irwork], info);
2087 }
2088 else
2089 {
2090 /* Perform bidiagonal QR iteration, if desired, computing */
2091 /* left singular vectors in A and computing right singular */
2092 /* vectors in VT */
2093 /* (CWorkspace: 0) */
2094 /* (RWorkspace: need BDSPAC) */
2095 zbdsqr_("U", n, &ncvt, &nru, &c__0, &s[1], &rwork[ie], &vt[ vt_offset], ldvt, &a[a_offset], lda, cdum, &c__1, & rwork[irwork], info);
2096 }
2097 }
2098 }
2099 else
2100 {
2101 /* A has more columns than rows. If A has sufficiently more */
2102 /* columns than rows, first reduce using the LQ decomposition (if */
2103 /* sufficient workspace available) */
2104 if (*n >= mnthr)
2105 {
2106 if (wntvn)
2107 {
2108 /* Path 1t(N much larger than M, JOBVT='N') */
2109 /* No right singular vectors to be computed */
2110 itau = 1;
2111 iwork = itau + *m;
2112 /* Compute A=L*Q */
2113 /* (CWorkspace: need 2*M, prefer M+M*NB) */
2114 /* (RWorkspace: 0) */
2115 i__2 = *lwork - iwork + 1;
2116 zgelqf_(m, n, &a[a_offset], lda, &work[itau], &work[iwork], & i__2, &ierr);
2117 /* Zero out above L */
2118 i__2 = *m - 1;
2119 i__3 = *m - 1;
2120 zlaset_("U", &i__2, &i__3, &c_b1, &c_b1, &a[(a_dim1 << 1) + 1] , lda);
2121 ie = 1;
2122 itauq = 1;
2123 itaup = itauq + *m;
2124 iwork = itaup + *m;
2125 /* Bidiagonalize L in A */
2126 /* (CWorkspace: need 3*M, prefer 2*M+2*M*NB) */
2127 /* (RWorkspace: need M) */
2128 i__2 = *lwork - iwork + 1;
2129 zgebrd_(m, m, &a[a_offset], lda, &s[1], &rwork[ie], &work[ itauq], &work[itaup], &work[iwork], &i__2, &ierr);
2130 if (wntuo || wntuas)
2131 {
2132 /* If left singular vectors desired, generate Q */
2133 /* (CWorkspace: need 3*M, prefer 2*M+M*NB) */
2134 /* (RWorkspace: 0) */
2135 i__2 = *lwork - iwork + 1;
2136 zungbr_("Q", m, m, m, &a[a_offset], lda, &work[itauq], & work[iwork], &i__2, &ierr);
2137 }
2138 irwork = ie + *m;
2139 nru = 0;
2140 if (wntuo || wntuas)
2141 {
2142 nru = *m;
2143 }
2144 /* Perform bidiagonal QR iteration, computing left singular */
2145 /* vectors of A in A if desired */
2146 /* (CWorkspace: 0) */
2147 /* (RWorkspace: need BDSPAC) */
2148 zbdsqr_("U", m, &c__0, &nru, &c__0, &s[1], &rwork[ie], cdum, & c__1, &a[a_offset], lda, cdum, &c__1, &rwork[irwork], info);
2149 /* If left singular vectors desired in U, copy them there */
2150 if (wntuas)
2151 {
2152 zlacpy_("F", m, m, &a[a_offset], lda, &u[u_offset], ldu);
2153 }
2154 }
2155 else if (wntvo && wntun)
2156 {
2157 /* Path 2t(N much larger than M, JOBU='N', JOBVT='O') */
2158 /* M right singular vectors to be overwritten on A and */
2159 /* no left singular vectors to be computed */
2160 if (*lwork >= *m * *m + *m * 3)
2161 {
2162 /* Sufficient workspace for a fast algorithm */
2163 ir = 1;
2164 /* Computing MAX */
2165 i__2 = wrkbl;
2166 i__3 = *lda * *n; // , expr subst
2167 if (*lwork >= max(i__2,i__3) + *lda * *m)
2168 {
2169 /* WORK(IU) is LDA by N and WORK(IR) is LDA by M */
2170 ldwrku = *lda;
2171 chunk = *n;
2172 ldwrkr = *lda;
2173 }
2174 else /* if(complicated condition) */
2175 {
2176 /* Computing MAX */
2177 i__2 = wrkbl;
2178 i__3 = *lda * *n; // , expr subst
2179 if (*lwork >= max(i__2,i__3) + *m * *m)
2180 {
2181 /* WORK(IU) is LDA by N and WORK(IR) is M by M */
2182 ldwrku = *lda;
2183 chunk = *n;
2184 ldwrkr = *m;
2185 }
2186 else
2187 {
2188 /* WORK(IU) is M by CHUNK and WORK(IR) is M by M */
2189 ldwrku = *m;
2190 chunk = (*lwork - *m * *m) / *m;
2191 ldwrkr = *m;
2192 }
2193 }
2194 itau = ir + ldwrkr * *m;
2195 iwork = itau + *m;
2196 /* Compute A=L*Q */
2197 /* (CWorkspace: need M*M+2*M, prefer M*M+M+M*NB) */
2198 /* (RWorkspace: 0) */
2199 i__2 = *lwork - iwork + 1;
2200 zgelqf_(m, n, &a[a_offset], lda, &work[itau], &work[iwork] , &i__2, &ierr);
2201 /* Copy L to WORK(IR) and zero out above it */
2202 zlacpy_("L", m, m, &a[a_offset], lda, &work[ir], &ldwrkr);
2203 i__2 = *m - 1;
2204 i__3 = *m - 1;
2205 zlaset_("U", &i__2, &i__3, &c_b1, &c_b1, &work[ir + ldwrkr], &ldwrkr);
2206 /* Generate Q in A */
2207 /* (CWorkspace: need M*M+2*M, prefer M*M+M+M*NB) */
2208 /* (RWorkspace: 0) */
2209 i__2 = *lwork - iwork + 1;
2210 zunglq_(m, n, m, &a[a_offset], lda, &work[itau], &work[ iwork], &i__2, &ierr);
2211 ie = 1;
2212 itauq = itau;
2213 itaup = itauq + *m;
2214 iwork = itaup + *m;
2215 /* Bidiagonalize L in WORK(IR) */
2216 /* (CWorkspace: need M*M+3*M, prefer M*M+2*M+2*M*NB) */
2217 /* (RWorkspace: need M) */
2218 i__2 = *lwork - iwork + 1;
2219 zgebrd_(m, m, &work[ir], &ldwrkr, &s[1], &rwork[ie], & work[itauq], &work[itaup], &work[iwork], &i__2, & ierr);
2220 /* Generate right vectors bidiagonalizing L */
2221 /* (CWorkspace: need M*M+3*M-1, prefer M*M+2*M+(M-1)*NB) */
2222 /* (RWorkspace: 0) */
2223 i__2 = *lwork - iwork + 1;
2224 zungbr_("P", m, m, m, &work[ir], &ldwrkr, &work[itaup], & work[iwork], &i__2, &ierr);
2225 irwork = ie + *m;
2226 /* Perform bidiagonal QR iteration, computing right */
2227 /* singular vectors of L in WORK(IR) */
2228 /* (CWorkspace: need M*M) */
2229 /* (RWorkspace: need BDSPAC) */
2230 zbdsqr_("U", m, m, &c__0, &c__0, &s[1], &rwork[ie], &work[ ir], &ldwrkr, cdum, &c__1, cdum, &c__1, &rwork[ irwork], info);
2231 iu = itauq;
2232 /* Multiply right singular vectors of L in WORK(IR) by Q */
2233 /* in A, storing result in WORK(IU) and copying to A */
2234 /* (CWorkspace: need M*M+M, prefer M*M+M*N) */
2235 /* (RWorkspace: 0) */
2236 i__2 = *n;
2237 i__3 = chunk;
2238 for (i__ = 1;
2239 i__3 < 0 ? i__ >= i__2 : i__ <= i__2;
2240 i__ += i__3)
2241 {
2242 /* Computing MIN */
2243 i__4 = *n - i__ + 1;
2244 blk = min(i__4,chunk);
2245 zgemm_("N", "N", m, &blk, m, &c_b2, &work[ir], & ldwrkr, &a[i__ * a_dim1 + 1], lda, &c_b1, & work[iu], &ldwrku);
2246 zlacpy_("F", m, &blk, &work[iu], &ldwrku, &a[i__ * a_dim1 + 1], lda);
2247 /* L30: */
2248 }
2249 }
2250 else
2251 {
2252 /* Insufficient workspace for a fast algorithm */
2253 ie = 1;
2254 itauq = 1;
2255 itaup = itauq + *m;
2256 iwork = itaup + *m;
2257 /* Bidiagonalize A */
2258 /* (CWorkspace: need 2*M+N, prefer 2*M+(M+N)*NB) */
2259 /* (RWorkspace: need M) */
2260 i__3 = *lwork - iwork + 1;
2261 zgebrd_(m, n, &a[a_offset], lda, &s[1], &rwork[ie], &work[ itauq], &work[itaup], &work[iwork], &i__3, &ierr);
2262 /* Generate right vectors bidiagonalizing A */
2263 /* (CWorkspace: need 3*M, prefer 2*M+M*NB) */
2264 /* (RWorkspace: 0) */
2265 i__3 = *lwork - iwork + 1;
2266 zungbr_("P", m, n, m, &a[a_offset], lda, &work[itaup], & work[iwork], &i__3, &ierr);
2267 irwork = ie + *m;
2268 /* Perform bidiagonal QR iteration, computing right */
2269 /* singular vectors of A in A */
2270 /* (CWorkspace: 0) */
2271 /* (RWorkspace: need BDSPAC) */
2272 zbdsqr_("L", m, n, &c__0, &c__0, &s[1], &rwork[ie], &a[ a_offset], lda, cdum, &c__1, cdum, &c__1, &rwork[ irwork], info);
2273 }
2274 }
2275 else if (wntvo && wntuas)
2276 {
2277 /* Path 3t(N much larger than M, JOBU='S' or 'A', JOBVT='O') */
2278 /* M right singular vectors to be overwritten on A and */
2279 /* M left singular vectors to be computed in U */
2280 if (*lwork >= *m * *m + *m * 3)
2281 {
2282 /* Sufficient workspace for a fast algorithm */
2283 ir = 1;
2284 /* Computing MAX */
2285 i__3 = wrkbl;
2286 i__2 = *lda * *n; // , expr subst
2287 if (*lwork >= max(i__3,i__2) + *lda * *m)
2288 {
2289 /* WORK(IU) is LDA by N and WORK(IR) is LDA by M */
2290 ldwrku = *lda;
2291 chunk = *n;
2292 ldwrkr = *lda;
2293 }
2294 else /* if(complicated condition) */
2295 {
2296 /* Computing MAX */
2297 i__3 = wrkbl;
2298 i__2 = *lda * *n; // , expr subst
2299 if (*lwork >= max(i__3,i__2) + *m * *m)
2300 {
2301 /* WORK(IU) is LDA by N and WORK(IR) is M by M */
2302 ldwrku = *lda;
2303 chunk = *n;
2304 ldwrkr = *m;
2305 }
2306 else
2307 {
2308 /* WORK(IU) is M by CHUNK and WORK(IR) is M by M */
2309 ldwrku = *m;
2310 chunk = (*lwork - *m * *m) / *m;
2311 ldwrkr = *m;
2312 }
2313 }
2314 itau = ir + ldwrkr * *m;
2315 iwork = itau + *m;
2316 /* Compute A=L*Q */
2317 /* (CWorkspace: need M*M+2*M, prefer M*M+M+M*NB) */
2318 /* (RWorkspace: 0) */
2319 i__3 = *lwork - iwork + 1;
2320 zgelqf_(m, n, &a[a_offset], lda, &work[itau], &work[iwork] , &i__3, &ierr);
2321 /* Copy L to U, zeroing about above it */
2322 zlacpy_("L", m, m, &a[a_offset], lda, &u[u_offset], ldu);
2323 i__3 = *m - 1;
2324 i__2 = *m - 1;
2325 zlaset_("U", &i__3, &i__2, &c_b1, &c_b1, &u[(u_dim1 << 1) + 1], ldu);
2326 /* Generate Q in A */
2327 /* (CWorkspace: need M*M+2*M, prefer M*M+M+M*NB) */
2328 /* (RWorkspace: 0) */
2329 i__3 = *lwork - iwork + 1;
2330 zunglq_(m, n, m, &a[a_offset], lda, &work[itau], &work[ iwork], &i__3, &ierr);
2331 ie = 1;
2332 itauq = itau;
2333 itaup = itauq + *m;
2334 iwork = itaup + *m;
2335 /* Bidiagonalize L in U, copying result to WORK(IR) */
2336 /* (CWorkspace: need M*M+3*M, prefer M*M+2*M+2*M*NB) */
2337 /* (RWorkspace: need M) */
2338 i__3 = *lwork - iwork + 1;
2339 zgebrd_(m, m, &u[u_offset], ldu, &s[1], &rwork[ie], &work[ itauq], &work[itaup], &work[iwork], &i__3, &ierr);
2340 zlacpy_("U", m, m, &u[u_offset], ldu, &work[ir], &ldwrkr);
2341 /* Generate right vectors bidiagonalizing L in WORK(IR) */
2342 /* (CWorkspace: need M*M+3*M-1, prefer M*M+2*M+(M-1)*NB) */
2343 /* (RWorkspace: 0) */
2344 i__3 = *lwork - iwork + 1;
2345 zungbr_("P", m, m, m, &work[ir], &ldwrkr, &work[itaup], & work[iwork], &i__3, &ierr);
2346 /* Generate left vectors bidiagonalizing L in U */
2347 /* (CWorkspace: need M*M+3*M, prefer M*M+2*M+M*NB) */
2348 /* (RWorkspace: 0) */
2349 i__3 = *lwork - iwork + 1;
2350 zungbr_("Q", m, m, m, &u[u_offset], ldu, &work[itauq], & work[iwork], &i__3, &ierr);
2351 irwork = ie + *m;
2352 /* Perform bidiagonal QR iteration, computing left */
2353 /* singular vectors of L in U, and computing right */
2354 /* singular vectors of L in WORK(IR) */
2355 /* (CWorkspace: need M*M) */
2356 /* (RWorkspace: need BDSPAC) */
2357 zbdsqr_("U", m, m, m, &c__0, &s[1], &rwork[ie], &work[ir], &ldwrkr, &u[u_offset], ldu, cdum, &c__1, &rwork[ irwork], info);
2358 iu = itauq;
2359 /* Multiply right singular vectors of L in WORK(IR) by Q */
2360 /* in A, storing result in WORK(IU) and copying to A */
2361 /* (CWorkspace: need M*M+M, prefer M*M+M*N)) */
2362 /* (RWorkspace: 0) */
2363 i__3 = *n;
2364 i__2 = chunk;
2365 for (i__ = 1;
2366 i__2 < 0 ? i__ >= i__3 : i__ <= i__3;
2367 i__ += i__2)
2368 {
2369 /* Computing MIN */
2370 i__4 = *n - i__ + 1;
2371 blk = min(i__4,chunk);
2372 zgemm_("N", "N", m, &blk, m, &c_b2, &work[ir], & ldwrkr, &a[i__ * a_dim1 + 1], lda, &c_b1, & work[iu], &ldwrku);
2373 zlacpy_("F", m, &blk, &work[iu], &ldwrku, &a[i__ * a_dim1 + 1], lda);
2374 /* L40: */
2375 }
2376 }
2377 else
2378 {
2379 /* Insufficient workspace for a fast algorithm */
2380 itau = 1;
2381 iwork = itau + *m;
2382 /* Compute A=L*Q */
2383 /* (CWorkspace: need 2*M, prefer M+M*NB) */
2384 /* (RWorkspace: 0) */
2385 i__2 = *lwork - iwork + 1;
2386 zgelqf_(m, n, &a[a_offset], lda, &work[itau], &work[iwork] , &i__2, &ierr);
2387 /* Copy L to U, zeroing out above it */
2388 zlacpy_("L", m, m, &a[a_offset], lda, &u[u_offset], ldu);
2389 i__2 = *m - 1;
2390 i__3 = *m - 1;
2391 zlaset_("U", &i__2, &i__3, &c_b1, &c_b1, &u[(u_dim1 << 1) + 1], ldu);
2392 /* Generate Q in A */
2393 /* (CWorkspace: need 2*M, prefer M+M*NB) */
2394 /* (RWorkspace: 0) */
2395 i__2 = *lwork - iwork + 1;
2396 zunglq_(m, n, m, &a[a_offset], lda, &work[itau], &work[ iwork], &i__2, &ierr);
2397 ie = 1;
2398 itauq = itau;
2399 itaup = itauq + *m;
2400 iwork = itaup + *m;
2401 /* Bidiagonalize L in U */
2402 /* (CWorkspace: need 3*M, prefer 2*M+2*M*NB) */
2403 /* (RWorkspace: need M) */
2404 i__2 = *lwork - iwork + 1;
2405 zgebrd_(m, m, &u[u_offset], ldu, &s[1], &rwork[ie], &work[ itauq], &work[itaup], &work[iwork], &i__2, &ierr);
2406 /* Multiply right vectors bidiagonalizing L by Q in A */
2407 /* (CWorkspace: need 2*M+N, prefer 2*M+N*NB) */
2408 /* (RWorkspace: 0) */
2409 i__2 = *lwork - iwork + 1;
2410 zunmbr_("P", "L", "C", m, n, m, &u[u_offset], ldu, &work[ itaup], &a[a_offset], lda, &work[iwork], &i__2, & ierr);
2411 /* Generate left vectors bidiagonalizing L in U */
2412 /* (CWorkspace: need 3*M, prefer 2*M+M*NB) */
2413 /* (RWorkspace: 0) */
2414 i__2 = *lwork - iwork + 1;
2415 zungbr_("Q", m, m, m, &u[u_offset], ldu, &work[itauq], & work[iwork], &i__2, &ierr);
2416 irwork = ie + *m;
2417 /* Perform bidiagonal QR iteration, computing left */
2418 /* singular vectors of A in U and computing right */
2419 /* singular vectors of A in A */
2420 /* (CWorkspace: 0) */
2421 /* (RWorkspace: need BDSPAC) */
2422 zbdsqr_("U", m, n, m, &c__0, &s[1], &rwork[ie], &a[ a_offset], lda, &u[u_offset], ldu, cdum, &c__1, & rwork[irwork], info);
2423 }
2424 }
2425 else if (wntvs)
2426 {
2427 if (wntun)
2428 {
2429 /* Path 4t(N much larger than M, JOBU='N', JOBVT='S') */
2430 /* M right singular vectors to be computed in VT and */
2431 /* no left singular vectors to be computed */
2432 if (*lwork >= *m * *m + *m * 3)
2433 {
2434 /* Sufficient workspace for a fast algorithm */
2435 ir = 1;
2436 if (*lwork >= wrkbl + *lda * *m)
2437 {
2438 /* WORK(IR) is LDA by M */
2439 ldwrkr = *lda;
2440 }
2441 else
2442 {
2443 /* WORK(IR) is M by M */
2444 ldwrkr = *m;
2445 }
2446 itau = ir + ldwrkr * *m;
2447 iwork = itau + *m;
2448 /* Compute A=L*Q */
2449 /* (CWorkspace: need M*M+2*M, prefer M*M+M+M*NB) */
2450 /* (RWorkspace: 0) */
2451 i__2 = *lwork - iwork + 1;
2452 zgelqf_(m, n, &a[a_offset], lda, &work[itau], &work[ iwork], &i__2, &ierr);
2453 /* Copy L to WORK(IR), zeroing out above it */
2454 zlacpy_("L", m, m, &a[a_offset], lda, &work[ir], & ldwrkr);
2455 i__2 = *m - 1;
2456 i__3 = *m - 1;
2457 zlaset_("U", &i__2, &i__3, &c_b1, &c_b1, &work[ir + ldwrkr], &ldwrkr);
2458 /* Generate Q in A */
2459 /* (CWorkspace: need M*M+2*M, prefer M*M+M+M*NB) */
2460 /* (RWorkspace: 0) */
2461 i__2 = *lwork - iwork + 1;
2462 zunglq_(m, n, m, &a[a_offset], lda, &work[itau], & work[iwork], &i__2, &ierr);
2463 ie = 1;
2464 itauq = itau;
2465 itaup = itauq + *m;
2466 iwork = itaup + *m;
2467 /* Bidiagonalize L in WORK(IR) */
2468 /* (CWorkspace: need M*M+3*M, prefer M*M+2*M+2*M*NB) */
2469 /* (RWorkspace: need M) */
2470 i__2 = *lwork - iwork + 1;
2471 zgebrd_(m, m, &work[ir], &ldwrkr, &s[1], &rwork[ie], & work[itauq], &work[itaup], &work[iwork], & i__2, &ierr);
2472 /* Generate right vectors bidiagonalizing L in */
2473 /* WORK(IR) */
2474 /* (CWorkspace: need M*M+3*M, prefer M*M+2*M+(M-1)*NB) */
2475 /* (RWorkspace: 0) */
2476 i__2 = *lwork - iwork + 1;
2477 zungbr_("P", m, m, m, &work[ir], &ldwrkr, &work[itaup] , &work[iwork], &i__2, &ierr);
2478 irwork = ie + *m;
2479 /* Perform bidiagonal QR iteration, computing right */
2480 /* singular vectors of L in WORK(IR) */
2481 /* (CWorkspace: need M*M) */
2482 /* (RWorkspace: need BDSPAC) */
2483 zbdsqr_("U", m, m, &c__0, &c__0, &s[1], &rwork[ie], & work[ir], &ldwrkr, cdum, &c__1, cdum, &c__1, & rwork[irwork], info);
2484 /* Multiply right singular vectors of L in WORK(IR) by */
2485 /* Q in A, storing result in VT */
2486 /* (CWorkspace: need M*M) */
2487 /* (RWorkspace: 0) */
2488 zgemm_("N", "N", m, n, m, &c_b2, &work[ir], &ldwrkr, & a[a_offset], lda, &c_b1, &vt[vt_offset], ldvt);
2489 }
2490 else
2491 {
2492 /* Insufficient workspace for a fast algorithm */
2493 itau = 1;
2494 iwork = itau + *m;
2495 /* Compute A=L*Q */
2496 /* (CWorkspace: need 2*M, prefer M+M*NB) */
2497 /* (RWorkspace: 0) */
2498 i__2 = *lwork - iwork + 1;
2499 zgelqf_(m, n, &a[a_offset], lda, &work[itau], &work[ iwork], &i__2, &ierr);
2500 /* Copy result to VT */
2501 zlacpy_("U", m, n, &a[a_offset], lda, &vt[vt_offset], ldvt);
2502 /* Generate Q in VT */
2503 /* (CWorkspace: need 2*M, prefer M+M*NB) */
2504 /* (RWorkspace: 0) */
2505 i__2 = *lwork - iwork + 1;
2506 zunglq_(m, n, m, &vt[vt_offset], ldvt, &work[itau], & work[iwork], &i__2, &ierr);
2507 ie = 1;
2508 itauq = itau;
2509 itaup = itauq + *m;
2510 iwork = itaup + *m;
2511 /* Zero out above L in A */
2512 i__2 = *m - 1;
2513 i__3 = *m - 1;
2514 zlaset_("U", &i__2, &i__3, &c_b1, &c_b1, &a[(a_dim1 << 1) + 1], lda);
2515 /* Bidiagonalize L in A */
2516 /* (CWorkspace: need 3*M, prefer 2*M+2*M*NB) */
2517 /* (RWorkspace: need M) */
2518 i__2 = *lwork - iwork + 1;
2519 zgebrd_(m, m, &a[a_offset], lda, &s[1], &rwork[ie], & work[itauq], &work[itaup], &work[iwork], & i__2, &ierr);
2520 /* Multiply right vectors bidiagonalizing L by Q in VT */
2521 /* (CWorkspace: need 2*M+N, prefer 2*M+N*NB) */
2522 /* (RWorkspace: 0) */
2523 i__2 = *lwork - iwork + 1;
2524 zunmbr_("P", "L", "C", m, n, m, &a[a_offset], lda, & work[itaup], &vt[vt_offset], ldvt, &work[ iwork], &i__2, &ierr);
2525 irwork = ie + *m;
2526 /* Perform bidiagonal QR iteration, computing right */
2527 /* singular vectors of A in VT */
2528 /* (CWorkspace: 0) */
2529 /* (RWorkspace: need BDSPAC) */
2530 zbdsqr_("U", m, n, &c__0, &c__0, &s[1], &rwork[ie], & vt[vt_offset], ldvt, cdum, &c__1, cdum, &c__1, &rwork[irwork], info);
2531 }
2532 }
2533 else if (wntuo)
2534 {
2535 /* Path 5t(N much larger than M, JOBU='O', JOBVT='S') */
2536 /* M right singular vectors to be computed in VT and */
2537 /* M left singular vectors to be overwritten on A */
2538 if (*lwork >= (*m << 1) * *m + *m * 3)
2539 {
2540 /* Sufficient workspace for a fast algorithm */
2541 iu = 1;
2542 if (*lwork >= wrkbl + (*lda << 1) * *m)
2543 {
2544 /* WORK(IU) is LDA by M and WORK(IR) is LDA by M */
2545 ldwrku = *lda;
2546 ir = iu + ldwrku * *m;
2547 ldwrkr = *lda;
2548 }
2549 else if (*lwork >= wrkbl + (*lda + *m) * *m)
2550 {
2551 /* WORK(IU) is LDA by M and WORK(IR) is M by M */
2552 ldwrku = *lda;
2553 ir = iu + ldwrku * *m;
2554 ldwrkr = *m;
2555 }
2556 else
2557 {
2558 /* WORK(IU) is M by M and WORK(IR) is M by M */
2559 ldwrku = *m;
2560 ir = iu + ldwrku * *m;
2561 ldwrkr = *m;
2562 }
2563 itau = ir + ldwrkr * *m;
2564 iwork = itau + *m;
2565 /* Compute A=L*Q */
2566 /* (CWorkspace: need 2*M*M+2*M, prefer 2*M*M+M+M*NB) */
2567 /* (RWorkspace: 0) */
2568 i__2 = *lwork - iwork + 1;
2569 zgelqf_(m, n, &a[a_offset], lda, &work[itau], &work[ iwork], &i__2, &ierr);
2570 /* Copy L to WORK(IU), zeroing out below it */
2571 zlacpy_("L", m, m, &a[a_offset], lda, &work[iu], & ldwrku);
2572 i__2 = *m - 1;
2573 i__3 = *m - 1;
2574 zlaset_("U", &i__2, &i__3, &c_b1, &c_b1, &work[iu + ldwrku], &ldwrku);
2575 /* Generate Q in A */
2576 /* (CWorkspace: need 2*M*M+2*M, prefer 2*M*M+M+M*NB) */
2577 /* (RWorkspace: 0) */
2578 i__2 = *lwork - iwork + 1;
2579 zunglq_(m, n, m, &a[a_offset], lda, &work[itau], & work[iwork], &i__2, &ierr);
2580 ie = 1;
2581 itauq = itau;
2582 itaup = itauq + *m;
2583 iwork = itaup + *m;
2584 /* Bidiagonalize L in WORK(IU), copying result to */
2585 /* WORK(IR) */
2586 /* (CWorkspace: need 2*M*M+3*M, */
2587 /* prefer 2*M*M+2*M+2*M*NB) */
2588 /* (RWorkspace: need M) */
2589 i__2 = *lwork - iwork + 1;
2590 zgebrd_(m, m, &work[iu], &ldwrku, &s[1], &rwork[ie], & work[itauq], &work[itaup], &work[iwork], & i__2, &ierr);
2591 zlacpy_("L", m, m, &work[iu], &ldwrku, &work[ir], & ldwrkr);
2592 /* Generate right bidiagonalizing vectors in WORK(IU) */
2593 /* (CWorkspace: need 2*M*M+3*M-1, */
2594 /* prefer 2*M*M+2*M+(M-1)*NB) */
2595 /* (RWorkspace: 0) */
2596 i__2 = *lwork - iwork + 1;
2597 zungbr_("P", m, m, m, &work[iu], &ldwrku, &work[itaup] , &work[iwork], &i__2, &ierr);
2598 /* Generate left bidiagonalizing vectors in WORK(IR) */
2599 /* (CWorkspace: need 2*M*M+3*M, prefer 2*M*M+2*M+M*NB) */
2600 /* (RWorkspace: 0) */
2601 i__2 = *lwork - iwork + 1;
2602 zungbr_("Q", m, m, m, &work[ir], &ldwrkr, &work[itauq] , &work[iwork], &i__2, &ierr);
2603 irwork = ie + *m;
2604 /* Perform bidiagonal QR iteration, computing left */
2605 /* singular vectors of L in WORK(IR) and computing */
2606 /* right singular vectors of L in WORK(IU) */
2607 /* (CWorkspace: need 2*M*M) */
2608 /* (RWorkspace: need BDSPAC) */
2609 zbdsqr_("U", m, m, m, &c__0, &s[1], &rwork[ie], &work[ iu], &ldwrku, &work[ir], &ldwrkr, cdum, &c__1, &rwork[irwork], info);
2610 /* Multiply right singular vectors of L in WORK(IU) by */
2611 /* Q in A, storing result in VT */
2612 /* (CWorkspace: need M*M) */
2613 /* (RWorkspace: 0) */
2614 zgemm_("N", "N", m, n, m, &c_b2, &work[iu], &ldwrku, & a[a_offset], lda, &c_b1, &vt[vt_offset], ldvt);
2615 /* Copy left singular vectors of L to A */
2616 /* (CWorkspace: need M*M) */
2617 /* (RWorkspace: 0) */
2618 zlacpy_("F", m, m, &work[ir], &ldwrkr, &a[a_offset], lda);
2619 }
2620 else
2621 {
2622 /* Insufficient workspace for a fast algorithm */
2623 itau = 1;
2624 iwork = itau + *m;
2625 /* Compute A=L*Q, copying result to VT */
2626 /* (CWorkspace: need 2*M, prefer M+M*NB) */
2627 /* (RWorkspace: 0) */
2628 i__2 = *lwork - iwork + 1;
2629 zgelqf_(m, n, &a[a_offset], lda, &work[itau], &work[ iwork], &i__2, &ierr);
2630 zlacpy_("U", m, n, &a[a_offset], lda, &vt[vt_offset], ldvt);
2631 /* Generate Q in VT */
2632 /* (CWorkspace: need 2*M, prefer M+M*NB) */
2633 /* (RWorkspace: 0) */
2634 i__2 = *lwork - iwork + 1;
2635 zunglq_(m, n, m, &vt[vt_offset], ldvt, &work[itau], & work[iwork], &i__2, &ierr);
2636 ie = 1;
2637 itauq = itau;
2638 itaup = itauq + *m;
2639 iwork = itaup + *m;
2640 /* Zero out above L in A */
2641 i__2 = *m - 1;
2642 i__3 = *m - 1;
2643 zlaset_("U", &i__2, &i__3, &c_b1, &c_b1, &a[(a_dim1 << 1) + 1], lda);
2644 /* Bidiagonalize L in A */
2645 /* (CWorkspace: need 3*M, prefer 2*M+2*M*NB) */
2646 /* (RWorkspace: need M) */
2647 i__2 = *lwork - iwork + 1;
2648 zgebrd_(m, m, &a[a_offset], lda, &s[1], &rwork[ie], & work[itauq], &work[itaup], &work[iwork], & i__2, &ierr);
2649 /* Multiply right vectors bidiagonalizing L by Q in VT */
2650 /* (CWorkspace: need 2*M+N, prefer 2*M+N*NB) */
2651 /* (RWorkspace: 0) */
2652 i__2 = *lwork - iwork + 1;
2653 zunmbr_("P", "L", "C", m, n, m, &a[a_offset], lda, & work[itaup], &vt[vt_offset], ldvt, &work[ iwork], &i__2, &ierr);
2654 /* Generate left bidiagonalizing vectors of L in A */
2655 /* (CWorkspace: need 3*M, prefer 2*M+M*NB) */
2656 /* (RWorkspace: 0) */
2657 i__2 = *lwork - iwork + 1;
2658 zungbr_("Q", m, m, m, &a[a_offset], lda, &work[itauq], &work[iwork], &i__2, &ierr);
2659 irwork = ie + *m;
2660 /* Perform bidiagonal QR iteration, computing left */
2661 /* singular vectors of A in A and computing right */
2662 /* singular vectors of A in VT */
2663 /* (CWorkspace: 0) */
2664 /* (RWorkspace: need BDSPAC) */
2665 zbdsqr_("U", m, n, m, &c__0, &s[1], &rwork[ie], &vt[ vt_offset], ldvt, &a[a_offset], lda, cdum, & c__1, &rwork[irwork], info);
2666 }
2667 }
2668 else if (wntuas)
2669 {
2670 /* Path 6t(N much larger than M, JOBU='S' or 'A', */
2671 /* JOBVT='S') */
2672 /* M right singular vectors to be computed in VT and */
2673 /* M left singular vectors to be computed in U */
2674 if (*lwork >= *m * *m + *m * 3)
2675 {
2676 /* Sufficient workspace for a fast algorithm */
2677 iu = 1;
2678 if (*lwork >= wrkbl + *lda * *m)
2679 {
2680 /* WORK(IU) is LDA by N */
2681 ldwrku = *lda;
2682 }
2683 else
2684 {
2685 /* WORK(IU) is LDA by M */
2686 ldwrku = *m;
2687 }
2688 itau = iu + ldwrku * *m;
2689 iwork = itau + *m;
2690 /* Compute A=L*Q */
2691 /* (CWorkspace: need M*M+2*M, prefer M*M+M+M*NB) */
2692 /* (RWorkspace: 0) */
2693 i__2 = *lwork - iwork + 1;
2694 zgelqf_(m, n, &a[a_offset], lda, &work[itau], &work[ iwork], &i__2, &ierr);
2695 /* Copy L to WORK(IU), zeroing out above it */
2696 zlacpy_("L", m, m, &a[a_offset], lda, &work[iu], & ldwrku);
2697 i__2 = *m - 1;
2698 i__3 = *m - 1;
2699 zlaset_("U", &i__2, &i__3, &c_b1, &c_b1, &work[iu + ldwrku], &ldwrku);
2700 /* Generate Q in A */
2701 /* (CWorkspace: need M*M+2*M, prefer M*M+M+M*NB) */
2702 /* (RWorkspace: 0) */
2703 i__2 = *lwork - iwork + 1;
2704 zunglq_(m, n, m, &a[a_offset], lda, &work[itau], & work[iwork], &i__2, &ierr);
2705 ie = 1;
2706 itauq = itau;
2707 itaup = itauq + *m;
2708 iwork = itaup + *m;
2709 /* Bidiagonalize L in WORK(IU), copying result to U */
2710 /* (CWorkspace: need M*M+3*M, prefer M*M+2*M+2*M*NB) */
2711 /* (RWorkspace: need M) */
2712 i__2 = *lwork - iwork + 1;
2713 zgebrd_(m, m, &work[iu], &ldwrku, &s[1], &rwork[ie], & work[itauq], &work[itaup], &work[iwork], & i__2, &ierr);
2714 zlacpy_("L", m, m, &work[iu], &ldwrku, &u[u_offset], ldu);
2715 /* Generate right bidiagonalizing vectors in WORK(IU) */
2716 /* (CWorkspace: need M*M+3*M-1, */
2717 /* prefer M*M+2*M+(M-1)*NB) */
2718 /* (RWorkspace: 0) */
2719 i__2 = *lwork - iwork + 1;
2720 zungbr_("P", m, m, m, &work[iu], &ldwrku, &work[itaup] , &work[iwork], &i__2, &ierr);
2721 /* Generate left bidiagonalizing vectors in U */
2722 /* (CWorkspace: need M*M+3*M, prefer M*M+2*M+M*NB) */
2723 /* (RWorkspace: 0) */
2724 i__2 = *lwork - iwork + 1;
2725 zungbr_("Q", m, m, m, &u[u_offset], ldu, &work[itauq], &work[iwork], &i__2, &ierr);
2726 irwork = ie + *m;
2727 /* Perform bidiagonal QR iteration, computing left */
2728 /* singular vectors of L in U and computing right */
2729 /* singular vectors of L in WORK(IU) */
2730 /* (CWorkspace: need M*M) */
2731 /* (RWorkspace: need BDSPAC) */
2732 zbdsqr_("U", m, m, m, &c__0, &s[1], &rwork[ie], &work[ iu], &ldwrku, &u[u_offset], ldu, cdum, &c__1, &rwork[irwork], info);
2733 /* Multiply right singular vectors of L in WORK(IU) by */
2734 /* Q in A, storing result in VT */
2735 /* (CWorkspace: need M*M) */
2736 /* (RWorkspace: 0) */
2737 zgemm_("N", "N", m, n, m, &c_b2, &work[iu], &ldwrku, & a[a_offset], lda, &c_b1, &vt[vt_offset], ldvt);
2738 }
2739 else
2740 {
2741 /* Insufficient workspace for a fast algorithm */
2742 itau = 1;
2743 iwork = itau + *m;
2744 /* Compute A=L*Q, copying result to VT */
2745 /* (CWorkspace: need 2*M, prefer M+M*NB) */
2746 /* (RWorkspace: 0) */
2747 i__2 = *lwork - iwork + 1;
2748 zgelqf_(m, n, &a[a_offset], lda, &work[itau], &work[ iwork], &i__2, &ierr);
2749 zlacpy_("U", m, n, &a[a_offset], lda, &vt[vt_offset], ldvt);
2750 /* Generate Q in VT */
2751 /* (CWorkspace: need 2*M, prefer M+M*NB) */
2752 /* (RWorkspace: 0) */
2753 i__2 = *lwork - iwork + 1;
2754 zunglq_(m, n, m, &vt[vt_offset], ldvt, &work[itau], & work[iwork], &i__2, &ierr);
2755 /* Copy L to U, zeroing out above it */
2756 zlacpy_("L", m, m, &a[a_offset], lda, &u[u_offset], ldu);
2757 i__2 = *m - 1;
2758 i__3 = *m - 1;
2759 zlaset_("U", &i__2, &i__3, &c_b1, &c_b1, &u[(u_dim1 << 1) + 1], ldu);
2760 ie = 1;
2761 itauq = itau;
2762 itaup = itauq + *m;
2763 iwork = itaup + *m;
2764 /* Bidiagonalize L in U */
2765 /* (CWorkspace: need 3*M, prefer 2*M+2*M*NB) */
2766 /* (RWorkspace: need M) */
2767 i__2 = *lwork - iwork + 1;
2768 zgebrd_(m, m, &u[u_offset], ldu, &s[1], &rwork[ie], & work[itauq], &work[itaup], &work[iwork], & i__2, &ierr);
2769 /* Multiply right bidiagonalizing vectors in U by Q */
2770 /* in VT */
2771 /* (CWorkspace: need 2*M+N, prefer 2*M+N*NB) */
2772 /* (RWorkspace: 0) */
2773 i__2 = *lwork - iwork + 1;
2774 zunmbr_("P", "L", "C", m, n, m, &u[u_offset], ldu, & work[itaup], &vt[vt_offset], ldvt, &work[ iwork], &i__2, &ierr);
2775 /* Generate left bidiagonalizing vectors in U */
2776 /* (CWorkspace: need 3*M, prefer 2*M+M*NB) */
2777 /* (RWorkspace: 0) */
2778 i__2 = *lwork - iwork + 1;
2779 zungbr_("Q", m, m, m, &u[u_offset], ldu, &work[itauq], &work[iwork], &i__2, &ierr);
2780 irwork = ie + *m;
2781 /* Perform bidiagonal QR iteration, computing left */
2782 /* singular vectors of A in U and computing right */
2783 /* singular vectors of A in VT */
2784 /* (CWorkspace: 0) */
2785 /* (RWorkspace: need BDSPAC) */
2786 zbdsqr_("U", m, n, m, &c__0, &s[1], &rwork[ie], &vt[ vt_offset], ldvt, &u[u_offset], ldu, cdum, & c__1, &rwork[irwork], info);
2787 }
2788 }
2789 }
2790 else if (wntva)
2791 {
2792 if (wntun)
2793 {
2794 /* Path 7t(N much larger than M, JOBU='N', JOBVT='A') */
2795 /* N right singular vectors to be computed in VT and */
2796 /* no left singular vectors to be computed */
2797 /* Computing MAX */
2798 i__2 = *n + *m;
2799 i__3 = *m * 3; // , expr subst
2800 if (*lwork >= *m * *m + max(i__2,i__3))
2801 {
2802 /* Sufficient workspace for a fast algorithm */
2803 ir = 1;
2804 if (*lwork >= wrkbl + *lda * *m)
2805 {
2806 /* WORK(IR) is LDA by M */
2807 ldwrkr = *lda;
2808 }
2809 else
2810 {
2811 /* WORK(IR) is M by M */
2812 ldwrkr = *m;
2813 }
2814 itau = ir + ldwrkr * *m;
2815 iwork = itau + *m;
2816 /* Compute A=L*Q, copying result to VT */
2817 /* (CWorkspace: need M*M+2*M, prefer M*M+M+M*NB) */
2818 /* (RWorkspace: 0) */
2819 i__2 = *lwork - iwork + 1;
2820 zgelqf_(m, n, &a[a_offset], lda, &work[itau], &work[ iwork], &i__2, &ierr);
2821 zlacpy_("U", m, n, &a[a_offset], lda, &vt[vt_offset], ldvt);
2822 /* Copy L to WORK(IR), zeroing out above it */
2823 zlacpy_("L", m, m, &a[a_offset], lda, &work[ir], & ldwrkr);
2824 i__2 = *m - 1;
2825 i__3 = *m - 1;
2826 zlaset_("U", &i__2, &i__3, &c_b1, &c_b1, &work[ir + ldwrkr], &ldwrkr);
2827 /* Generate Q in VT */
2828 /* (CWorkspace: need M*M+M+N, prefer M*M+M+N*NB) */
2829 /* (RWorkspace: 0) */
2830 i__2 = *lwork - iwork + 1;
2831 zunglq_(n, n, m, &vt[vt_offset], ldvt, &work[itau], & work[iwork], &i__2, &ierr);
2832 ie = 1;
2833 itauq = itau;
2834 itaup = itauq + *m;
2835 iwork = itaup + *m;
2836 /* Bidiagonalize L in WORK(IR) */
2837 /* (CWorkspace: need M*M+3*M, prefer M*M+2*M+2*M*NB) */
2838 /* (RWorkspace: need M) */
2839 i__2 = *lwork - iwork + 1;
2840 zgebrd_(m, m, &work[ir], &ldwrkr, &s[1], &rwork[ie], & work[itauq], &work[itaup], &work[iwork], & i__2, &ierr);
2841 /* Generate right bidiagonalizing vectors in WORK(IR) */
2842 /* (CWorkspace: need M*M+3*M-1, */
2843 /* prefer M*M+2*M+(M-1)*NB) */
2844 /* (RWorkspace: 0) */
2845 i__2 = *lwork - iwork + 1;
2846 zungbr_("P", m, m, m, &work[ir], &ldwrkr, &work[itaup] , &work[iwork], &i__2, &ierr);
2847 irwork = ie + *m;
2848 /* Perform bidiagonal QR iteration, computing right */
2849 /* singular vectors of L in WORK(IR) */
2850 /* (CWorkspace: need M*M) */
2851 /* (RWorkspace: need BDSPAC) */
2852 zbdsqr_("U", m, m, &c__0, &c__0, &s[1], &rwork[ie], & work[ir], &ldwrkr, cdum, &c__1, cdum, &c__1, & rwork[irwork], info);
2853 /* Multiply right singular vectors of L in WORK(IR) by */
2854 /* Q in VT, storing result in A */
2855 /* (CWorkspace: need M*M) */
2856 /* (RWorkspace: 0) */
2857 zgemm_("N", "N", m, n, m, &c_b2, &work[ir], &ldwrkr, & vt[vt_offset], ldvt, &c_b1, &a[a_offset], lda);
2858 /* Copy right singular vectors of A from A to VT */
2859 zlacpy_("F", m, n, &a[a_offset], lda, &vt[vt_offset], ldvt);
2860 }
2861 else
2862 {
2863 /* Insufficient workspace for a fast algorithm */
2864 itau = 1;
2865 iwork = itau + *m;
2866 /* Compute A=L*Q, copying result to VT */
2867 /* (CWorkspace: need 2*M, prefer M+M*NB) */
2868 /* (RWorkspace: 0) */
2869 i__2 = *lwork - iwork + 1;
2870 zgelqf_(m, n, &a[a_offset], lda, &work[itau], &work[ iwork], &i__2, &ierr);
2871 zlacpy_("U", m, n, &a[a_offset], lda, &vt[vt_offset], ldvt);
2872 /* Generate Q in VT */
2873 /* (CWorkspace: need M+N, prefer M+N*NB) */
2874 /* (RWorkspace: 0) */
2875 i__2 = *lwork - iwork + 1;
2876 zunglq_(n, n, m, &vt[vt_offset], ldvt, &work[itau], & work[iwork], &i__2, &ierr);
2877 ie = 1;
2878 itauq = itau;
2879 itaup = itauq + *m;
2880 iwork = itaup + *m;
2881 /* Zero out above L in A */
2882 i__2 = *m - 1;
2883 i__3 = *m - 1;
2884 zlaset_("U", &i__2, &i__3, &c_b1, &c_b1, &a[(a_dim1 << 1) + 1], lda);
2885 /* Bidiagonalize L in A */
2886 /* (CWorkspace: need 3*M, prefer 2*M+2*M*NB) */
2887 /* (RWorkspace: need M) */
2888 i__2 = *lwork - iwork + 1;
2889 zgebrd_(m, m, &a[a_offset], lda, &s[1], &rwork[ie], & work[itauq], &work[itaup], &work[iwork], & i__2, &ierr);
2890 /* Multiply right bidiagonalizing vectors in A by Q */
2891 /* in VT */
2892 /* (CWorkspace: need 2*M+N, prefer 2*M+N*NB) */
2893 /* (RWorkspace: 0) */
2894 i__2 = *lwork - iwork + 1;
2895 zunmbr_("P", "L", "C", m, n, m, &a[a_offset], lda, & work[itaup], &vt[vt_offset], ldvt, &work[ iwork], &i__2, &ierr);
2896 irwork = ie + *m;
2897 /* Perform bidiagonal QR iteration, computing right */
2898 /* singular vectors of A in VT */
2899 /* (CWorkspace: 0) */
2900 /* (RWorkspace: need BDSPAC) */
2901 zbdsqr_("U", m, n, &c__0, &c__0, &s[1], &rwork[ie], & vt[vt_offset], ldvt, cdum, &c__1, cdum, &c__1, &rwork[irwork], info);
2902 }
2903 }
2904 else if (wntuo)
2905 {
2906 /* Path 8t(N much larger than M, JOBU='O', JOBVT='A') */
2907 /* N right singular vectors to be computed in VT and */
2908 /* M left singular vectors to be overwritten on A */
2909 /* Computing MAX */
2910 i__2 = *n + *m;
2911 i__3 = *m * 3; // , expr subst
2912 if (*lwork >= (*m << 1) * *m + max(i__2,i__3))
2913 {
2914 /* Sufficient workspace for a fast algorithm */
2915 iu = 1;
2916 if (*lwork >= wrkbl + (*lda << 1) * *m)
2917 {
2918 /* WORK(IU) is LDA by M and WORK(IR) is LDA by M */
2919 ldwrku = *lda;
2920 ir = iu + ldwrku * *m;
2921 ldwrkr = *lda;
2922 }
2923 else if (*lwork >= wrkbl + (*lda + *m) * *m)
2924 {
2925 /* WORK(IU) is LDA by M and WORK(IR) is M by M */
2926 ldwrku = *lda;
2927 ir = iu + ldwrku * *m;
2928 ldwrkr = *m;
2929 }
2930 else
2931 {
2932 /* WORK(IU) is M by M and WORK(IR) is M by M */
2933 ldwrku = *m;
2934 ir = iu + ldwrku * *m;
2935 ldwrkr = *m;
2936 }
2937 itau = ir + ldwrkr * *m;
2938 iwork = itau + *m;
2939 /* Compute A=L*Q, copying result to VT */
2940 /* (CWorkspace: need 2*M*M+2*M, prefer 2*M*M+M+M*NB) */
2941 /* (RWorkspace: 0) */
2942 i__2 = *lwork - iwork + 1;
2943 zgelqf_(m, n, &a[a_offset], lda, &work[itau], &work[ iwork], &i__2, &ierr);
2944 zlacpy_("U", m, n, &a[a_offset], lda, &vt[vt_offset], ldvt);
2945 /* Generate Q in VT */
2946 /* (CWorkspace: need 2*M*M+M+N, prefer 2*M*M+M+N*NB) */
2947 /* (RWorkspace: 0) */
2948 i__2 = *lwork - iwork + 1;
2949 zunglq_(n, n, m, &vt[vt_offset], ldvt, &work[itau], & work[iwork], &i__2, &ierr);
2950 /* Copy L to WORK(IU), zeroing out above it */
2951 zlacpy_("L", m, m, &a[a_offset], lda, &work[iu], & ldwrku);
2952 i__2 = *m - 1;
2953 i__3 = *m - 1;
2954 zlaset_("U", &i__2, &i__3, &c_b1, &c_b1, &work[iu + ldwrku], &ldwrku);
2955 ie = 1;
2956 itauq = itau;
2957 itaup = itauq + *m;
2958 iwork = itaup + *m;
2959 /* Bidiagonalize L in WORK(IU), copying result to */
2960 /* WORK(IR) */
2961 /* (CWorkspace: need 2*M*M+3*M, */
2962 /* prefer 2*M*M+2*M+2*M*NB) */
2963 /* (RWorkspace: need M) */
2964 i__2 = *lwork - iwork + 1;
2965 zgebrd_(m, m, &work[iu], &ldwrku, &s[1], &rwork[ie], & work[itauq], &work[itaup], &work[iwork], & i__2, &ierr);
2966 zlacpy_("L", m, m, &work[iu], &ldwrku, &work[ir], & ldwrkr);
2967 /* Generate right bidiagonalizing vectors in WORK(IU) */
2968 /* (CWorkspace: need 2*M*M+3*M-1, */
2969 /* prefer 2*M*M+2*M+(M-1)*NB) */
2970 /* (RWorkspace: 0) */
2971 i__2 = *lwork - iwork + 1;
2972 zungbr_("P", m, m, m, &work[iu], &ldwrku, &work[itaup] , &work[iwork], &i__2, &ierr);
2973 /* Generate left bidiagonalizing vectors in WORK(IR) */
2974 /* (CWorkspace: need 2*M*M+3*M, prefer 2*M*M+2*M+M*NB) */
2975 /* (RWorkspace: 0) */
2976 i__2 = *lwork - iwork + 1;
2977 zungbr_("Q", m, m, m, &work[ir], &ldwrkr, &work[itauq] , &work[iwork], &i__2, &ierr);
2978 irwork = ie + *m;
2979 /* Perform bidiagonal QR iteration, computing left */
2980 /* singular vectors of L in WORK(IR) and computing */
2981 /* right singular vectors of L in WORK(IU) */
2982 /* (CWorkspace: need 2*M*M) */
2983 /* (RWorkspace: need BDSPAC) */
2984 zbdsqr_("U", m, m, m, &c__0, &s[1], &rwork[ie], &work[ iu], &ldwrku, &work[ir], &ldwrkr, cdum, &c__1, &rwork[irwork], info);
2985 /* Multiply right singular vectors of L in WORK(IU) by */
2986 /* Q in VT, storing result in A */
2987 /* (CWorkspace: need M*M) */
2988 /* (RWorkspace: 0) */
2989 zgemm_("N", "N", m, n, m, &c_b2, &work[iu], &ldwrku, & vt[vt_offset], ldvt, &c_b1, &a[a_offset], lda);
2990 /* Copy right singular vectors of A from A to VT */
2991 zlacpy_("F", m, n, &a[a_offset], lda, &vt[vt_offset], ldvt);
2992 /* Copy left singular vectors of A from WORK(IR) to A */
2993 zlacpy_("F", m, m, &work[ir], &ldwrkr, &a[a_offset], lda);
2994 }
2995 else
2996 {
2997 /* Insufficient workspace for a fast algorithm */
2998 itau = 1;
2999 iwork = itau + *m;
3000 /* Compute A=L*Q, copying result to VT */
3001 /* (CWorkspace: need 2*M, prefer M+M*NB) */
3002 /* (RWorkspace: 0) */
3003 i__2 = *lwork - iwork + 1;
3004 zgelqf_(m, n, &a[a_offset], lda, &work[itau], &work[ iwork], &i__2, &ierr);
3005 zlacpy_("U", m, n, &a[a_offset], lda, &vt[vt_offset], ldvt);
3006 /* Generate Q in VT */
3007 /* (CWorkspace: need M+N, prefer M+N*NB) */
3008 /* (RWorkspace: 0) */
3009 i__2 = *lwork - iwork + 1;
3010 zunglq_(n, n, m, &vt[vt_offset], ldvt, &work[itau], & work[iwork], &i__2, &ierr);
3011 ie = 1;
3012 itauq = itau;
3013 itaup = itauq + *m;
3014 iwork = itaup + *m;
3015 /* Zero out above L in A */
3016 i__2 = *m - 1;
3017 i__3 = *m - 1;
3018 zlaset_("U", &i__2, &i__3, &c_b1, &c_b1, &a[(a_dim1 << 1) + 1], lda);
3019 /* Bidiagonalize L in A */
3020 /* (CWorkspace: need 3*M, prefer 2*M+2*M*NB) */
3021 /* (RWorkspace: need M) */
3022 i__2 = *lwork - iwork + 1;
3023 zgebrd_(m, m, &a[a_offset], lda, &s[1], &rwork[ie], & work[itauq], &work[itaup], &work[iwork], & i__2, &ierr);
3024 /* Multiply right bidiagonalizing vectors in A by Q */
3025 /* in VT */
3026 /* (CWorkspace: need 2*M+N, prefer 2*M+N*NB) */
3027 /* (RWorkspace: 0) */
3028 i__2 = *lwork - iwork + 1;
3029 zunmbr_("P", "L", "C", m, n, m, &a[a_offset], lda, & work[itaup], &vt[vt_offset], ldvt, &work[ iwork], &i__2, &ierr);
3030 /* Generate left bidiagonalizing vectors in A */
3031 /* (CWorkspace: need 3*M, prefer 2*M+M*NB) */
3032 /* (RWorkspace: 0) */
3033 i__2 = *lwork - iwork + 1;
3034 zungbr_("Q", m, m, m, &a[a_offset], lda, &work[itauq], &work[iwork], &i__2, &ierr);
3035 irwork = ie + *m;
3036 /* Perform bidiagonal QR iteration, computing left */
3037 /* singular vectors of A in A and computing right */
3038 /* singular vectors of A in VT */
3039 /* (CWorkspace: 0) */
3040 /* (RWorkspace: need BDSPAC) */
3041 zbdsqr_("U", m, n, m, &c__0, &s[1], &rwork[ie], &vt[ vt_offset], ldvt, &a[a_offset], lda, cdum, & c__1, &rwork[irwork], info);
3042 }
3043 }
3044 else if (wntuas)
3045 {
3046 /* Path 9t(N much larger than M, JOBU='S' or 'A', */
3047 /* JOBVT='A') */
3048 /* N right singular vectors to be computed in VT and */
3049 /* M left singular vectors to be computed in U */
3050 /* Computing MAX */
3051 i__2 = *n + *m;
3052 i__3 = *m * 3; // , expr subst
3053 if (*lwork >= *m * *m + max(i__2,i__3))
3054 {
3055 /* Sufficient workspace for a fast algorithm */
3056 iu = 1;
3057 if (*lwork >= wrkbl + *lda * *m)
3058 {
3059 /* WORK(IU) is LDA by M */
3060 ldwrku = *lda;
3061 }
3062 else
3063 {
3064 /* WORK(IU) is M by M */
3065 ldwrku = *m;
3066 }
3067 itau = iu + ldwrku * *m;
3068 iwork = itau + *m;
3069 /* Compute A=L*Q, copying result to VT */
3070 /* (CWorkspace: need M*M+2*M, prefer M*M+M+M*NB) */
3071 /* (RWorkspace: 0) */
3072 i__2 = *lwork - iwork + 1;
3073 zgelqf_(m, n, &a[a_offset], lda, &work[itau], &work[ iwork], &i__2, &ierr);
3074 zlacpy_("U", m, n, &a[a_offset], lda, &vt[vt_offset], ldvt);
3075 /* Generate Q in VT */
3076 /* (CWorkspace: need M*M+M+N, prefer M*M+M+N*NB) */
3077 /* (RWorkspace: 0) */
3078 i__2 = *lwork - iwork + 1;
3079 zunglq_(n, n, m, &vt[vt_offset], ldvt, &work[itau], & work[iwork], &i__2, &ierr);
3080 /* Copy L to WORK(IU), zeroing out above it */
3081 zlacpy_("L", m, m, &a[a_offset], lda, &work[iu], & ldwrku);
3082 i__2 = *m - 1;
3083 i__3 = *m - 1;
3084 zlaset_("U", &i__2, &i__3, &c_b1, &c_b1, &work[iu + ldwrku], &ldwrku);
3085 ie = 1;
3086 itauq = itau;
3087 itaup = itauq + *m;
3088 iwork = itaup + *m;
3089 /* Bidiagonalize L in WORK(IU), copying result to U */
3090 /* (CWorkspace: need M*M+3*M, prefer M*M+2*M+2*M*NB) */
3091 /* (RWorkspace: need M) */
3092 i__2 = *lwork - iwork + 1;
3093 zgebrd_(m, m, &work[iu], &ldwrku, &s[1], &rwork[ie], & work[itauq], &work[itaup], &work[iwork], & i__2, &ierr);
3094 zlacpy_("L", m, m, &work[iu], &ldwrku, &u[u_offset], ldu);
3095 /* Generate right bidiagonalizing vectors in WORK(IU) */
3096 /* (CWorkspace: need M*M+3*M, prefer M*M+2*M+(M-1)*NB) */
3097 /* (RWorkspace: 0) */
3098 i__2 = *lwork - iwork + 1;
3099 zungbr_("P", m, m, m, &work[iu], &ldwrku, &work[itaup] , &work[iwork], &i__2, &ierr);
3100 /* Generate left bidiagonalizing vectors in U */
3101 /* (CWorkspace: need M*M+3*M, prefer M*M+2*M+M*NB) */
3102 /* (RWorkspace: 0) */
3103 i__2 = *lwork - iwork + 1;
3104 zungbr_("Q", m, m, m, &u[u_offset], ldu, &work[itauq], &work[iwork], &i__2, &ierr);
3105 irwork = ie + *m;
3106 /* Perform bidiagonal QR iteration, computing left */
3107 /* singular vectors of L in U and computing right */
3108 /* singular vectors of L in WORK(IU) */
3109 /* (CWorkspace: need M*M) */
3110 /* (RWorkspace: need BDSPAC) */
3111 zbdsqr_("U", m, m, m, &c__0, &s[1], &rwork[ie], &work[ iu], &ldwrku, &u[u_offset], ldu, cdum, &c__1, &rwork[irwork], info);
3112 /* Multiply right singular vectors of L in WORK(IU) by */
3113 /* Q in VT, storing result in A */
3114 /* (CWorkspace: need M*M) */
3115 /* (RWorkspace: 0) */
3116 zgemm_("N", "N", m, n, m, &c_b2, &work[iu], &ldwrku, & vt[vt_offset], ldvt, &c_b1, &a[a_offset], lda);
3117 /* Copy right singular vectors of A from A to VT */
3118 zlacpy_("F", m, n, &a[a_offset], lda, &vt[vt_offset], ldvt);
3119 }
3120 else
3121 {
3122 /* Insufficient workspace for a fast algorithm */
3123 itau = 1;
3124 iwork = itau + *m;
3125 /* Compute A=L*Q, copying result to VT */
3126 /* (CWorkspace: need 2*M, prefer M+M*NB) */
3127 /* (RWorkspace: 0) */
3128 i__2 = *lwork - iwork + 1;
3129 zgelqf_(m, n, &a[a_offset], lda, &work[itau], &work[ iwork], &i__2, &ierr);
3130 zlacpy_("U", m, n, &a[a_offset], lda, &vt[vt_offset], ldvt);
3131 /* Generate Q in VT */
3132 /* (CWorkspace: need M+N, prefer M+N*NB) */
3133 /* (RWorkspace: 0) */
3134 i__2 = *lwork - iwork + 1;
3135 zunglq_(n, n, m, &vt[vt_offset], ldvt, &work[itau], & work[iwork], &i__2, &ierr);
3136 /* Copy L to U, zeroing out above it */
3137 zlacpy_("L", m, m, &a[a_offset], lda, &u[u_offset], ldu);
3138 i__2 = *m - 1;
3139 i__3 = *m - 1;
3140 zlaset_("U", &i__2, &i__3, &c_b1, &c_b1, &u[(u_dim1 << 1) + 1], ldu);
3141 ie = 1;
3142 itauq = itau;
3143 itaup = itauq + *m;
3144 iwork = itaup + *m;
3145 /* Bidiagonalize L in U */
3146 /* (CWorkspace: need 3*M, prefer 2*M+2*M*NB) */
3147 /* (RWorkspace: need M) */
3148 i__2 = *lwork - iwork + 1;
3149 zgebrd_(m, m, &u[u_offset], ldu, &s[1], &rwork[ie], & work[itauq], &work[itaup], &work[iwork], & i__2, &ierr);
3150 /* Multiply right bidiagonalizing vectors in U by Q */
3151 /* in VT */
3152 /* (CWorkspace: need 2*M+N, prefer 2*M+N*NB) */
3153 /* (RWorkspace: 0) */
3154 i__2 = *lwork - iwork + 1;
3155 zunmbr_("P", "L", "C", m, n, m, &u[u_offset], ldu, & work[itaup], &vt[vt_offset], ldvt, &work[ iwork], &i__2, &ierr);
3156 /* Generate left bidiagonalizing vectors in U */
3157 /* (CWorkspace: need 3*M, prefer 2*M+M*NB) */
3158 /* (RWorkspace: 0) */
3159 i__2 = *lwork - iwork + 1;
3160 zungbr_("Q", m, m, m, &u[u_offset], ldu, &work[itauq], &work[iwork], &i__2, &ierr);
3161 irwork = ie + *m;
3162 /* Perform bidiagonal QR iteration, computing left */
3163 /* singular vectors of A in U and computing right */
3164 /* singular vectors of A in VT */
3165 /* (CWorkspace: 0) */
3166 /* (RWorkspace: need BDSPAC) */
3167 zbdsqr_("U", m, n, m, &c__0, &s[1], &rwork[ie], &vt[ vt_offset], ldvt, &u[u_offset], ldu, cdum, & c__1, &rwork[irwork], info);
3168 }
3169 }
3170 }
3171 }
3172 else
3173 {
3174 /* N .LT. MNTHR */
3175 /* Path 10t(N greater than M, but not much larger) */
3176 /* Reduce to bidiagonal form without LQ decomposition */
3177 ie = 1;
3178 itauq = 1;
3179 itaup = itauq + *m;
3180 iwork = itaup + *m;
3181 /* Bidiagonalize A */
3182 /* (CWorkspace: need 2*M+N, prefer 2*M+(M+N)*NB) */
3183 /* (RWorkspace: M) */
3184 i__2 = *lwork - iwork + 1;
3185 zgebrd_(m, n, &a[a_offset], lda, &s[1], &rwork[ie], &work[itauq], &work[itaup], &work[iwork], &i__2, &ierr);
3186 if (wntuas)
3187 {
3188 /* If left singular vectors desired in U, copy result to U */
3189 /* and generate left bidiagonalizing vectors in U */
3190 /* (CWorkspace: need 3*M-1, prefer 2*M+(M-1)*NB) */
3191 /* (RWorkspace: 0) */
3192 zlacpy_("L", m, m, &a[a_offset], lda, &u[u_offset], ldu);
3193 i__2 = *lwork - iwork + 1;
3194 zungbr_("Q", m, m, n, &u[u_offset], ldu, &work[itauq], &work[ iwork], &i__2, &ierr);
3195 }
3196 if (wntvas)
3197 {
3198 /* If right singular vectors desired in VT, copy result to */
3199 /* VT and generate right bidiagonalizing vectors in VT */
3200 /* (CWorkspace: need 2*M+NRVT, prefer 2*M+NRVT*NB) */
3201 /* (RWorkspace: 0) */
3202 zlacpy_("U", m, n, &a[a_offset], lda, &vt[vt_offset], ldvt);
3203 if (wntva)
3204 {
3205 nrvt = *n;
3206 }
3207 if (wntvs)
3208 {
3209 nrvt = *m;
3210 }
3211 i__2 = *lwork - iwork + 1;
3212 zungbr_("P", &nrvt, n, m, &vt[vt_offset], ldvt, &work[itaup], &work[iwork], &i__2, &ierr);
3213 }
3214 if (wntuo)
3215 {
3216 /* If left singular vectors desired in A, generate left */
3217 /* bidiagonalizing vectors in A */
3218 /* (CWorkspace: need 3*M-1, prefer 2*M+(M-1)*NB) */
3219 /* (RWorkspace: 0) */
3220 i__2 = *lwork - iwork + 1;
3221 zungbr_("Q", m, m, n, &a[a_offset], lda, &work[itauq], &work[ iwork], &i__2, &ierr);
3222 }
3223 if (wntvo)
3224 {
3225 /* If right singular vectors desired in A, generate right */
3226 /* bidiagonalizing vectors in A */
3227 /* (CWorkspace: need 3*M, prefer 2*M+M*NB) */
3228 /* (RWorkspace: 0) */
3229 i__2 = *lwork - iwork + 1;
3230 zungbr_("P", m, n, m, &a[a_offset], lda, &work[itaup], &work[ iwork], &i__2, &ierr);
3231 }
3232 irwork = ie + *m;
3233 if (wntuas || wntuo)
3234 {
3235 nru = *m;
3236 }
3237 if (wntun)
3238 {
3239 nru = 0;
3240 }
3241 if (wntvas || wntvo)
3242 {
3243 ncvt = *n;
3244 }
3245 if (wntvn)
3246 {
3247 ncvt = 0;
3248 }
3249 if (! wntuo && ! wntvo)
3250 {
3251 /* Perform bidiagonal QR iteration, if desired, computing */
3252 /* left singular vectors in U and computing right singular */
3253 /* vectors in VT */
3254 /* (CWorkspace: 0) */
3255 /* (RWorkspace: need BDSPAC) */
3256 zbdsqr_("L", m, &ncvt, &nru, &c__0, &s[1], &rwork[ie], &vt[ vt_offset], ldvt, &u[u_offset], ldu, cdum, &c__1, & rwork[irwork], info);
3257 }
3258 else if (! wntuo && wntvo)
3259 {
3260 /* Perform bidiagonal QR iteration, if desired, computing */
3261 /* left singular vectors in U and computing right singular */
3262 /* vectors in A */
3263 /* (CWorkspace: 0) */
3264 /* (RWorkspace: need BDSPAC) */
3265 zbdsqr_("L", m, &ncvt, &nru, &c__0, &s[1], &rwork[ie], &a[ a_offset], lda, &u[u_offset], ldu, cdum, &c__1, & rwork[irwork], info);
3266 }
3267 else
3268 {
3269 /* Perform bidiagonal QR iteration, if desired, computing */
3270 /* left singular vectors in A and computing right singular */
3271 /* vectors in VT */
3272 /* (CWorkspace: 0) */
3273 /* (RWorkspace: need BDSPAC) */
3274 zbdsqr_("L", m, &ncvt, &nru, &c__0, &s[1], &rwork[ie], &vt[ vt_offset], ldvt, &a[a_offset], lda, cdum, &c__1, & rwork[irwork], info);
3275 }
3276 }
3277 }
3278 /* Undo scaling if necessary */
3279 if (iscl == 1)
3280 {
3281 if (anrm > bignum)
3282 {
3283 dlascl_("G", &c__0, &c__0, &bignum, &anrm, &minmn, &c__1, &s[1], & minmn, &ierr);
3284 }
3285 if (*info != 0 && anrm > bignum)
3286 {
3287 i__2 = minmn - 1;
3288 dlascl_("G", &c__0, &c__0, &bignum, &anrm, &i__2, &c__1, &rwork[ ie], &minmn, &ierr);
3289 }
3290 if (anrm < smlnum)
3291 {
3292 dlascl_("G", &c__0, &c__0, &smlnum, &anrm, &minmn, &c__1, &s[1], & minmn, &ierr);
3293 }
3294 if (*info != 0 && anrm < smlnum)
3295 {
3296 i__2 = minmn - 1;
3297 dlascl_("G", &c__0, &c__0, &smlnum, &anrm, &i__2, &c__1, &rwork[ ie], &minmn, &ierr);
3298 }
3299 }
3300 /* Return optimal workspace in WORK(1) */
3301 work[1].r = (doublereal) maxwrk;
3302 work[1].i = 0.; // , expr subst
3303 return 0;
3304 /* End of ZGESVD */
3305 }
3306 /* zgesvd_ */
3307