nsz Git - libfirm/blob - ir/be/test/zgesvd/zgesvd.c

   1 #include "blaswrap.h"
   2 /*  -- translated by f2c (version 19990503).
   3    You must link the resulting object file with the libraries:
   4         -lf2c -lm   (in that order)
   5 */
   6
   7 #include "f2c.h"
   8
   9 /* Table of constant values */
  10
  11 static doublecomplex c_b1 = {0.,0.};
  12 static doublecomplex c_b2 = {1.,0.};
  13 static integer c__6 = 6;
  14 static integer c__0 = 0;
  15 static integer c__2 = 2;
  16 static integer c__1 = 1;
  17 static integer c_n1 = -1;
  18
  19 /* Subroutine */ int zgesvd_(char *jobu, char *jobvt, integer *m, integer *n,
  20         doublecomplex *a, integer *lda, doublereal *s, doublecomplex *u,
  21         integer *ldu, doublecomplex *vt, integer *ldvt, doublecomplex *work,
  22         integer *lwork, doublereal *rwork, integer *info)
  23 {
  24     /* System generated locals */
  25     address a__1[2];
  26     integer a_dim1, a_offset, u_dim1, u_offset, vt_dim1, vt_offset, i__1[2],
  27             i__2, i__3, i__4;
  28     char ch__1[2];
  29
  30     /* Builtin functions
  31        Subroutine */ int s_cat(char *, char **, integer *, integer *, ftnlen);
  32     double sqrt(doublereal);
  33
  34     /* Local variables */
  35     static doublecomplex cdum[1];
  36     static integer iscl;
  37     static doublereal anrm;
  38     static integer ierr, itau, ncvt, nrvt, i__;
  39     extern logical lsame_(char *, char *);
  40     static integer chunk, minmn;
  41     extern /* Subroutine */ int zgemm_(char *, char *, integer *, integer *,
  42             integer *, doublecomplex *, doublecomplex *, integer *,
  43             doublecomplex *, integer *, doublecomplex *, doublecomplex *,
  44             integer *);
  45     static integer wrkbl, itaup, itauq, mnthr, iwork;
  46     static logical wntua, wntva, wntun, wntuo, wntvn, wntvo, wntus, wntvs;
  47     static integer ie;
  48     extern doublereal dlamch_(char *);
  49     static integer ir, iu;
  50     extern /* Subroutine */ int dlascl_(char *, integer *, integer *,
  51             doublereal *, doublereal *, integer *, integer *, doublereal *,
  52             integer *, integer *), xerbla_(char *, integer *),
  53              zgebrd_(integer *, integer *, doublecomplex *, integer *,
  54             doublereal *, doublereal *, doublecomplex *, doublecomplex *,
  55             doublecomplex *, integer *, integer *);
  56     extern integer ilaenv_(integer *, char *, char *, integer *, integer *,
  57             integer *, integer *, ftnlen, ftnlen);
  58     extern doublereal zlange_(char *, integer *, integer *, doublecomplex *,
  59             integer *, doublereal *);
  60     static doublereal bignum;
  61     extern /* Subroutine */ int zgelqf_(integer *, integer *, doublecomplex *,
  62              integer *, doublecomplex *, doublecomplex *, integer *, integer *
  63             ), zlascl_(char *, integer *, integer *, doublereal *, doublereal
  64             *, integer *, integer *, doublecomplex *, integer *, integer *), zgeqrf_(integer *, integer *, doublecomplex *, integer *,
  65              doublecomplex *, doublecomplex *, integer *, integer *), zlacpy_(
  66             char *, integer *, integer *, doublecomplex *, integer *,
  67             doublecomplex *, integer *), zlaset_(char *, integer *,
  68             integer *, doublecomplex *, doublecomplex *, doublecomplex *,
  69             integer *);
  70     static integer ldwrkr;
  71     extern /* Subroutine */ int zbdsqr_(char *, integer *, integer *, integer
  72             *, integer *, doublereal *, doublereal *, doublecomplex *,
  73             integer *, doublecomplex *, integer *, doublecomplex *, integer *,
  74              doublereal *, integer *);
  75     static integer minwrk, ldwrku, maxwrk;
  76     extern /* Subroutine */ int zungbr_(char *, integer *, integer *, integer
  77             *, doublecomplex *, integer *, doublecomplex *, doublecomplex *,
  78             integer *, integer *);
  79     static doublereal smlnum;
  80     static integer irwork;
  81     extern /* Subroutine */ int zunmbr_(char *, char *, char *, integer *,
  82             integer *, integer *, doublecomplex *, integer *, doublecomplex *,
  83              doublecomplex *, integer *, doublecomplex *, integer *, integer *
  84             ), zunglq_(integer *, integer *, integer *
  85             , doublecomplex *, integer *, doublecomplex *, doublecomplex *,
  86             integer *, integer *);
  87     static logical lquery, wntuas, wntvas;
  88     extern /* Subroutine */ int zungqr_(integer *, integer *, integer *,
  89             doublecomplex *, integer *, doublecomplex *, doublecomplex *,
  90             integer *, integer *);
  91     static integer blk, ncu;
  92     static doublereal dum[1], eps;
  93     static integer nru;
  94
  95
  96 #define a_subscr(a_1,a_2) (a_2)*a_dim1 + a_1
  97 #define a_ref(a_1,a_2) a[a_subscr(a_1,a_2)]
  98 #define u_subscr(a_1,a_2) (a_2)*u_dim1 + a_1
  99 #define u_ref(a_1,a_2) u[u_subscr(a_1,a_2)]
 100 #define vt_subscr(a_1,a_2) (a_2)*vt_dim1 + a_1
 101 #define vt_ref(a_1,a_2) vt[vt_subscr(a_1,a_2)]
 102
 103
 104 /*  -- LAPACK driver routine (version 3.0) --
 105        Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
 106        Courant Institute, Argonne National Lab, and Rice University
 107        October 31, 1999
 108
 109
 110     Purpose
 111     =======
 112
 113     ZGESVD computes the singular value decomposition (SVD) of a complex
 114     M-by-N matrix A, optionally computing the left and/or right singular
 115     vectors. The SVD is written
 116
 117          A = U * SIGMA * conjugate-transpose(V)
 118
 119     where SIGMA is an M-by-N matrix which is zero except for its
 120     min(m,n) diagonal elements, U is an M-by-M unitary matrix, and
 121     V is an N-by-N unitary matrix.  The diagonal elements of SIGMA
 122     are the singular values of A; they are real and non-negative, and
 123     are returned in descending order.  The first min(m,n) columns of
 124     U and V are the left and right singular vectors of A.
 125
 126     Note that the routine returns V**H, not V.
 127
 128     Arguments
 129     =========
 130
 131     JOBU    (input) CHARACTER*1
 132             Specifies options for computing all or part of the matrix U:
 133             = 'A':  all M columns of U are returned in array U:
 134             = 'S':  the first min(m,n) columns of U (the left singular
 135                     vectors) are returned in the array U;
 136             = 'O':  the first min(m,n) columns of U (the left singular
 137                     vectors) are overwritten on the array A;
 138             = 'N':  no columns of U (no left singular vectors) are
 139                     computed.
 140
 141     JOBVT   (input) CHARACTER*1
 142             Specifies options for computing all or part of the matrix
 143             V**H:
 144             = 'A':  all N rows of V**H are returned in the array VT;
 145             = 'S':  the first min(m,n) rows of V**H (the right singular
 146                     vectors) are returned in the array VT;
 147             = 'O':  the first min(m,n) rows of V**H (the right singular
 148                     vectors) are overwritten on the array A;
 149             = 'N':  no rows of V**H (no right singular vectors) are
 150                     computed.
 151
 152             JOBVT and JOBU cannot both be 'O'.
 153
 154     M       (input) INTEGER
 155             The number of rows of the input matrix A.  M >= 0.
 156
 157     N       (input) INTEGER
 158             The number of columns of the input matrix A.  N >= 0.
 159
 160     A       (input/output) COMPLEX*16 array, dimension (LDA,N)
 161             On entry, the M-by-N matrix A.
 162             On exit,
 163             if JOBU = 'O',  A is overwritten with the first min(m,n)
 164                             columns of U (the left singular vectors,
 165                             stored columnwise);
 166             if JOBVT = 'O', A is overwritten with the first min(m,n)
 167                             rows of V**H (the right singular vectors,
 168                             stored rowwise);
 169             if JOBU .ne. 'O' and JOBVT .ne. 'O', the contents of A
 170                             are destroyed.
 171
 172     LDA     (input) INTEGER
 173             The leading dimension of the array A.  LDA >= max(1,M).
 174
 175     S       (output) DOUBLE PRECISION array, dimension (min(M,N))
 176             The singular values of A, sorted so that S(i) >= S(i+1).
 177
 178     U       (output) COMPLEX*16 array, dimension (LDU,UCOL)
 179             (LDU,M) if JOBU = 'A' or (LDU,min(M,N)) if JOBU = 'S'.
 180             If JOBU = 'A', U contains the M-by-M unitary matrix U;
 181             if JOBU = 'S', U contains the first min(m,n) columns of U
 182             (the left singular vectors, stored columnwise);
 183             if JOBU = 'N' or 'O', U is not referenced.
 184
 185     LDU     (input) INTEGER
 186             The leading dimension of the array U.  LDU >= 1; if
 187             JOBU = 'S' or 'A', LDU >= M.
 188
 189     VT      (output) COMPLEX*16 array, dimension (LDVT,N)
 190             If JOBVT = 'A', VT contains the N-by-N unitary matrix
 191             V**H;
 192             if JOBVT = 'S', VT contains the first min(m,n) rows of
 193             V**H (the right singular vectors, stored rowwise);
 194             if JOBVT = 'N' or 'O', VT is not referenced.
 195
 196     LDVT    (input) INTEGER
 197             The leading dimension of the array VT.  LDVT >= 1; if
 198             JOBVT = 'A', LDVT >= N; if JOBVT = 'S', LDVT >= min(M,N).
 199
 200     WORK    (workspace/output) COMPLEX*16 array, dimension (LWORK)
 201             On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
 202
 203     LWORK   (input) INTEGER
 204             The dimension of the array WORK. LWORK >= 1.
 205             LWORK >=  2*MIN(M,N)+MAX(M,N).
 206             For good performance, LWORK should generally be larger.
 207
 208             If LWORK = -1, then a workspace query is assumed; the routine
 209             only calculates the optimal size of the WORK array, returns
 210             this value as the first entry of the WORK array, and no error
 211             message related to LWORK is issued by XERBLA.
 212
 213     RWORK   (workspace) DOUBLE PRECISION array, dimension (5*min(M,N))
 214             On exit, if INFO > 0, RWORK(1:MIN(M,N)-1) contains the
 215             unconverged superdiagonal elements of an upper bidiagonal
 216             matrix B whose diagonal is in S (not necessarily sorted).
 217             B satisfies A = U * B * VT, so it has the same singular
 218             values as A, and singular vectors related by U and VT.
 219
 220     INFO    (output) INTEGER
 221             = 0:  successful exit.
 222             < 0:  if INFO = -i, the i-th argument had an illegal value.
 223             > 0:  if ZBDSQR did not converge, INFO specifies how many
 224                   superdiagonals of an intermediate bidiagonal form B
 225                   did not converge to zero. See the description of RWORK
 226                   above for details.
 227
 228     =====================================================================
 229
 230
 231        Test the input arguments
 232
 233        Parameter adjustments */
 234     a_dim1 = *lda;
 235     a_offset = 1 + a_dim1 * 1;
 236     a -= a_offset;
 237     --s;
 238     u_dim1 = *ldu;
 239     u_offset = 1 + u_dim1 * 1;
 240     u -= u_offset;
 241     vt_dim1 = *ldvt;
 242     vt_offset = 1 + vt_dim1 * 1;
 243     vt -= vt_offset;
 244     --work;
 245     --rwork;
 246
 247     /* Function Body */
 248     *info = 0;
 249     minmn = min(*m,*n);
 250 /* Writing concatenation */
 251     i__1[0] = 1, a__1[0] = jobu;
 252     i__1[1] = 1, a__1[1] = jobvt;
 253     s_cat(ch__1, a__1, i__1, &c__2, (ftnlen)2);
 254     mnthr = ilaenv_(&c__6, "ZGESVD", ch__1, m, n, &c__0, &c__0, (ftnlen)6, (
 255             ftnlen)2);
 256     wntua = lsame_(jobu, "A");
 257     wntus = lsame_(jobu, "S");
 258     wntuas = wntua || wntus;
 259     wntuo = lsame_(jobu, "O");
 260     wntun = lsame_(jobu, "N");
 261     wntva = lsame_(jobvt, "A");
 262     wntvs = lsame_(jobvt, "S");
 263     wntvas = wntva || wntvs;
 264     wntvo = lsame_(jobvt, "O");
 265     wntvn = lsame_(jobvt, "N");
 266     minwrk = 1;
 267     lquery = *lwork == -1;
 268
 269     if (! (wntua || wntus || wntuo || wntun)) {
 270         *info = -1;
 271     } else if (! (wntva || wntvs || wntvo || wntvn) || wntvo && wntuo) {
 272         *info = -2;
 273     } else if (*m < 0) {
 274         *info = -3;
 275     } else if (*n < 0) {
 276         *info = -4;
 277     } else if (*lda < max(1,*m)) {
 278         *info = -6;
 279     } else if (*ldu < 1 || wntuas && *ldu < *m) {
 280         *info = -9;
 281     } else if (*ldvt < 1 || wntva && *ldvt < *n || wntvs && *ldvt < minmn) {
 282         *info = -11;
 283     }
 284
 285 /*     Compute workspace
 286         (Note: Comments in the code beginning "Workspace:" describe the
 287          minimal amount of workspace needed at that point in the code,
 288          as well as the preferred amount for good performance.
 289          CWorkspace refers to complex workspace, and RWorkspace to
 290          real workspace. NB refers to the optimal block size for the
 291          immediately following subroutine, as returned by ILAENV.) */
 292
 293     if (*info == 0 && (*lwork >= 1 || lquery) && *m > 0 && *n > 0) {
 294         if (*m >= *n) {
 295
 296 /*           Space needed for ZBDSQR is BDSPAC = 5*N */
 297
 298             if (*m >= mnthr) {
 299                 if (wntun) {
 300
 301 /*                 Path 1 (M much larger than N, JOBU='N') */
 302
 303                     maxwrk = *n + *n * ilaenv_(&c__1, "ZGEQRF", " ", m, n, &
 304                             c_n1, &c_n1, (ftnlen)6, (ftnlen)1);
 305 /* Computing MAX */
 306                     i__2 = maxwrk, i__3 = (*n << 1) + (*n << 1) * ilaenv_(&
 307                             c__1, "ZGEBRD", " ", n, n, &c_n1, &c_n1, (ftnlen)
 308                             6, (ftnlen)1);
 309                     maxwrk = max(i__2,i__3);
 310                     if (wntvo || wntvas) {
 311 /* Computing MAX */
 312                         i__2 = maxwrk, i__3 = (*n << 1) + (*n - 1) * ilaenv_(&
 313                                 c__1, "ZUNGBR", "P", n, n, n, &c_n1, (ftnlen)
 314                                 6, (ftnlen)1);
 315                         maxwrk = max(i__2,i__3);
 316                     }
 317                     minwrk = *n * 3;
 318                     maxwrk = max(minwrk,maxwrk);
 319                 } else if (wntuo && wntvn) {
 320
 321 /*                 Path 2 (M much larger than N, JOBU='O', JOBVT='N') */
 322
 323                     wrkbl = *n + *n * ilaenv_(&c__1, "ZGEQRF", " ", m, n, &
 324                             c_n1, &c_n1, (ftnlen)6, (ftnlen)1);
 325 /* Computing MAX */
 326                     i__2 = wrkbl, i__3 = *n + *n * ilaenv_(&c__1, "ZUNGQR",
 327                             " ", m, n, n, &c_n1, (ftnlen)6, (ftnlen)1);
 328                     wrkbl = max(i__2,i__3);
 329 /* Computing MAX */
 330                     i__2 = wrkbl, i__3 = (*n << 1) + (*n << 1) * ilaenv_(&
 331                             c__1, "ZGEBRD", " ", n, n, &c_n1, &c_n1, (ftnlen)
 332                             6, (ftnlen)1);
 333                     wrkbl = max(i__2,i__3);
 334 /* Computing MAX */
 335                     i__2 = wrkbl, i__3 = (*n << 1) + *n * ilaenv_(&c__1,
 336                             "ZUNGBR", "Q", n, n, n, &c_n1, (ftnlen)6, (ftnlen)
 337                             1);
 338                     wrkbl = max(i__2,i__3);
 339 /* Computing MAX */
 340                     i__2 = *n * *n + wrkbl, i__3 = *n * *n + *m * *n;
 341                     maxwrk = max(i__2,i__3);
 342                     minwrk = (*n << 1) + *m;
 343                     maxwrk = max(minwrk,maxwrk);
 344                 } else if (wntuo && wntvas) {
 345
 346 /*                 Path 3 (M much larger than N, JOBU='O', JOBVT='S' or
 347                    'A') */
 348
 349                     wrkbl = *n + *n * ilaenv_(&c__1, "ZGEQRF", " ", m, n, &
 350                             c_n1, &c_n1, (ftnlen)6, (ftnlen)1);
 351 /* Computing MAX */
 352                     i__2 = wrkbl, i__3 = *n + *n * ilaenv_(&c__1, "ZUNGQR",
 353                             " ", m, n, n, &c_n1, (ftnlen)6, (ftnlen)1);
 354                     wrkbl = max(i__2,i__3);
 355 /* Computing MAX */
 356                     i__2 = wrkbl, i__3 = (*n << 1) + (*n << 1) * ilaenv_(&
 357                             c__1, "ZGEBRD", " ", n, n, &c_n1, &c_n1, (ftnlen)
 358                             6, (ftnlen)1);
 359                     wrkbl = max(i__2,i__3);
 360 /* Computing MAX */
 361                     i__2 = wrkbl, i__3 = (*n << 1) + *n * ilaenv_(&c__1,
 362                             "ZUNGBR", "Q", n, n, n, &c_n1, (ftnlen)6, (ftnlen)
 363                             1);
 364                     wrkbl = max(i__2,i__3);
 365 /* Computing MAX */
 366                     i__2 = wrkbl, i__3 = (*n << 1) + (*n - 1) * ilaenv_(&c__1,
 367                              "ZUNGBR", "P", n, n, n, &c_n1, (ftnlen)6, (
 368                             ftnlen)1);
 369                     wrkbl = max(i__2,i__3);
 370 /* Computing MAX */
 371                     i__2 = *n * *n + wrkbl, i__3 = *n * *n + *m * *n;
 372                     maxwrk = max(i__2,i__3);
 373                     minwrk = (*n << 1) + *m;
 374                     maxwrk = max(minwrk,maxwrk);
 375                 } else if (wntus && wntvn) {
 376
 377 /*                 Path 4 (M much larger than N, JOBU='S', JOBVT='N') */
 378
 379                     wrkbl = *n + *n * ilaenv_(&c__1, "ZGEQRF", " ", m, n, &
 380                             c_n1, &c_n1, (ftnlen)6, (ftnlen)1);
 381 /* Computing MAX */
 382                     i__2 = wrkbl, i__3 = *n + *n * ilaenv_(&c__1, "ZUNGQR",
 383                             " ", m, n, n, &c_n1, (ftnlen)6, (ftnlen)1);
 384                     wrkbl = max(i__2,i__3);
 385 /* Computing MAX */
 386                     i__2 = wrkbl, i__3 = (*n << 1) + (*n << 1) * ilaenv_(&
 387                             c__1, "ZGEBRD", " ", n, n, &c_n1, &c_n1, (ftnlen)
 388                             6, (ftnlen)1);
 389                     wrkbl = max(i__2,i__3);
 390 /* Computing MAX */
 391                     i__2 = wrkbl, i__3 = (*n << 1) + *n * ilaenv_(&c__1,
 392                             "ZUNGBR", "Q", n, n, n, &c_n1, (ftnlen)6, (ftnlen)
 393                             1);
 394                     wrkbl = max(i__2,i__3);
 395                     maxwrk = *n * *n + wrkbl;
 396                     minwrk = (*n << 1) + *m;
 397                     maxwrk = max(minwrk,maxwrk);
 398                 } else if (wntus && wntvo) {
 399
 400 /*                 Path 5 (M much larger than N, JOBU='S', JOBVT='O') */
 401
 402                     wrkbl = *n + *n * ilaenv_(&c__1, "ZGEQRF", " ", m, n, &
 403                             c_n1, &c_n1, (ftnlen)6, (ftnlen)1);
 404 /* Computing MAX */
 405                     i__2 = wrkbl, i__3 = *n + *n * ilaenv_(&c__1, "ZUNGQR",
 406                             " ", m, n, n, &c_n1, (ftnlen)6, (ftnlen)1);
 407                     wrkbl = max(i__2,i__3);
 408 /* Computing MAX */
 409                     i__2 = wrkbl, i__3 = (*n << 1) + (*n << 1) * ilaenv_(&
 410                             c__1, "ZGEBRD", " ", n, n, &c_n1, &c_n1, (ftnlen)
 411                             6, (ftnlen)1);
 412                     wrkbl = max(i__2,i__3);
 413 /* Computing MAX */
 414                     i__2 = wrkbl, i__3 = (*n << 1) + *n * ilaenv_(&c__1,
 415                             "ZUNGBR", "Q", n, n, n, &c_n1, (ftnlen)6, (ftnlen)
 416                             1);
 417                     wrkbl = max(i__2,i__3);
 418 /* Computing MAX */
 419                     i__2 = wrkbl, i__3 = (*n << 1) + (*n - 1) * ilaenv_(&c__1,
 420                              "ZUNGBR", "P", n, n, n, &c_n1, (ftnlen)6, (
 421                             ftnlen)1);
 422                     wrkbl = max(i__2,i__3);
 423                     maxwrk = (*n << 1) * *n + wrkbl;
 424                     minwrk = (*n << 1) + *m;
 425                     maxwrk = max(minwrk,maxwrk);
 426                 } else if (wntus && wntvas) {
 427
 428 /*                 Path 6 (M much larger than N, JOBU='S', JOBVT='S' or
 429                    'A') */
 430
 431                     wrkbl = *n + *n * ilaenv_(&c__1, "ZGEQRF", " ", m, n, &
 432                             c_n1, &c_n1, (ftnlen)6, (ftnlen)1);
 433 /* Computing MAX */
 434                     i__2 = wrkbl, i__3 = *n + *n * ilaenv_(&c__1, "ZUNGQR",
 435                             " ", m, n, n, &c_n1, (ftnlen)6, (ftnlen)1);
 436                     wrkbl = max(i__2,i__3);
 437 /* Computing MAX */
 438                     i__2 = wrkbl, i__3 = (*n << 1) + (*n << 1) * ilaenv_(&
 439                             c__1, "ZGEBRD", " ", n, n, &c_n1, &c_n1, (ftnlen)
 440                             6, (ftnlen)1);
 441                     wrkbl = max(i__2,i__3);
 442 /* Computing MAX */
 443                     i__2 = wrkbl, i__3 = (*n << 1) + *n * ilaenv_(&c__1,
 444                             "ZUNGBR", "Q", n, n, n, &c_n1, (ftnlen)6, (ftnlen)
 445                             1);
 446                     wrkbl = max(i__2,i__3);
 447 /* Computing MAX */
 448                     i__2 = wrkbl, i__3 = (*n << 1) + (*n - 1) * ilaenv_(&c__1,
 449                              "ZUNGBR", "P", n, n, n, &c_n1, (ftnlen)6, (
 450                             ftnlen)1);
 451                     wrkbl = max(i__2,i__3);
 452                     maxwrk = *n * *n + wrkbl;
 453                     minwrk = (*n << 1) + *m;
 454                     maxwrk = max(minwrk,maxwrk);
 455                 } else if (wntua && wntvn) {
 456
 457 /*                 Path 7 (M much larger than N, JOBU='A', JOBVT='N') */
 458
 459                     wrkbl = *n + *n * ilaenv_(&c__1, "ZGEQRF", " ", m, n, &
 460                             c_n1, &c_n1, (ftnlen)6, (ftnlen)1);
 461 /* Computing MAX */
 462                     i__2 = wrkbl, i__3 = *n + *m * ilaenv_(&c__1, "ZUNGQR",
 463                             " ", m, m, n, &c_n1, (ftnlen)6, (ftnlen)1);
 464                     wrkbl = max(i__2,i__3);
 465 /* Computing MAX */
 466                     i__2 = wrkbl, i__3 = (*n << 1) + (*n << 1) * ilaenv_(&
 467                             c__1, "ZGEBRD", " ", n, n, &c_n1, &c_n1, (ftnlen)
 468                             6, (ftnlen)1);
 469                     wrkbl = max(i__2,i__3);
 470 /* Computing MAX */
 471                     i__2 = wrkbl, i__3 = (*n << 1) + *n * ilaenv_(&c__1,
 472                             "ZUNGBR", "Q", n, n, n, &c_n1, (ftnlen)6, (ftnlen)
 473                             1);
 474                     wrkbl = max(i__2,i__3);
 475                     maxwrk = *n * *n + wrkbl;
 476                     minwrk = (*n << 1) + *m;
 477                     maxwrk = max(minwrk,maxwrk);
 478                 } else if (wntua && wntvo) {
 479
 480 /*                 Path 8 (M much larger than N, JOBU='A', JOBVT='O') */
 481
 482                     wrkbl = *n + *n * ilaenv_(&c__1, "ZGEQRF", " ", m, n, &
 483                             c_n1, &c_n1, (ftnlen)6, (ftnlen)1);
 484 /* Computing MAX */
 485                     i__2 = wrkbl, i__3 = *n + *m * ilaenv_(&c__1, "ZUNGQR",
 486                             " ", m, m, n, &c_n1, (ftnlen)6, (ftnlen)1);
 487                     wrkbl = max(i__2,i__3);
 488 /* Computing MAX */
 489                     i__2 = wrkbl, i__3 = (*n << 1) + (*n << 1) * ilaenv_(&
 490                             c__1, "ZGEBRD", " ", n, n, &c_n1, &c_n1, (ftnlen)
 491                             6, (ftnlen)1);
 492                     wrkbl = max(i__2,i__3);
 493 /* Computing MAX */
 494                     i__2 = wrkbl, i__3 = (*n << 1) + *n * ilaenv_(&c__1,
 495                             "ZUNGBR", "Q", n, n, n, &c_n1, (ftnlen)6, (ftnlen)
 496                             1);
 497                     wrkbl = max(i__2,i__3);
 498 /* Computing MAX */
 499                     i__2 = wrkbl, i__3 = (*n << 1) + (*n - 1) * ilaenv_(&c__1,
 500                              "ZUNGBR", "P", n, n, n, &c_n1, (ftnlen)6, (
 501                             ftnlen)1);
 502                     wrkbl = max(i__2,i__3);
 503                     maxwrk = (*n << 1) * *n + wrkbl;
 504                     minwrk = (*n << 1) + *m;
 505                     maxwrk = max(minwrk,maxwrk);
 506                 } else if (wntua && wntvas) {
 507
 508 /*                 Path 9 (M much larger than N, JOBU='A', JOBVT='S' or
 509                    'A') */
 510
 511                     wrkbl = *n + *n * ilaenv_(&c__1, "ZGEQRF", " ", m, n, &
 512                             c_n1, &c_n1, (ftnlen)6, (ftnlen)1);
 513 /* Computing MAX */
 514                     i__2 = wrkbl, i__3 = *n + *m * ilaenv_(&c__1, "ZUNGQR",
 515                             " ", m, m, n, &c_n1, (ftnlen)6, (ftnlen)1);
 516                     wrkbl = max(i__2,i__3);
 517 /* Computing MAX */
 518                     i__2 = wrkbl, i__3 = (*n << 1) + (*n << 1) * ilaenv_(&
 519                             c__1, "ZGEBRD", " ", n, n, &c_n1, &c_n1, (ftnlen)
 520                             6, (ftnlen)1);
 521                     wrkbl = max(i__2,i__3);
 522 /* Computing MAX */
 523                     i__2 = wrkbl, i__3 = (*n << 1) + *n * ilaenv_(&c__1,
 524                             "ZUNGBR", "Q", n, n, n, &c_n1, (ftnlen)6, (ftnlen)
 525                             1);
 526                     wrkbl = max(i__2,i__3);
 527 /* Computing MAX */
 528                     i__2 = wrkbl, i__3 = (*n << 1) + (*n - 1) * ilaenv_(&c__1,
 529                              "ZUNGBR", "P", n, n, n, &c_n1, (ftnlen)6, (
 530                             ftnlen)1);
 531                     wrkbl = max(i__2,i__3);
 532                     maxwrk = *n * *n + wrkbl;
 533                     minwrk = (*n << 1) + *m;
 534                     maxwrk = max(minwrk,maxwrk);
 535                 }
 536             } else {
 537
 538 /*              Path 10 (M at least N, but not much larger) */
 539
 540                 maxwrk = (*n << 1) + (*m + *n) * ilaenv_(&c__1, "ZGEBRD",
 541                         " ", m, n, &c_n1, &c_n1, (ftnlen)6, (ftnlen)1);
 542                 if (wntus || wntuo) {
 543 /* Computing MAX */
 544                     i__2 = maxwrk, i__3 = (*n << 1) + *n * ilaenv_(&c__1,
 545                             "ZUNGBR", "Q", m, n, n, &c_n1, (ftnlen)6, (ftnlen)
 546                             1);
 547                     maxwrk = max(i__2,i__3);
 548                 }
 549                 if (wntua) {
 550 /* Computing MAX */
 551                     i__2 = maxwrk, i__3 = (*n << 1) + *m * ilaenv_(&c__1,
 552                             "ZUNGBR", "Q", m, m, n, &c_n1, (ftnlen)6, (ftnlen)
 553                             1);
 554                     maxwrk = max(i__2,i__3);
 555                 }
 556                 if (! wntvn) {
 557 /* Computing MAX */
 558                     i__2 = maxwrk, i__3 = (*n << 1) + (*n - 1) * ilaenv_(&
 559                             c__1, "ZUNGBR", "P", n, n, n, &c_n1, (ftnlen)6, (
 560                             ftnlen)1);
 561                     maxwrk = max(i__2,i__3);
 562                 }
 563                 minwrk = (*n << 1) + *m;
 564                 maxwrk = max(minwrk,maxwrk);
 565             }
 566         } else {
 567
 568 /*           Space needed for ZBDSQR is BDSPAC = 5*M */
 569
 570             if (*n >= mnthr) {
 571                 if (wntvn) {
 572
 573 /*                 Path 1t(N much larger than M, JOBVT='N') */
 574
 575                     maxwrk = *m + *m * ilaenv_(&c__1, "ZGELQF", " ", m, n, &
 576                             c_n1, &c_n1, (ftnlen)6, (ftnlen)1);
 577 /* Computing MAX */
 578                     i__2 = maxwrk, i__3 = (*m << 1) + (*m << 1) * ilaenv_(&
 579                             c__1, "ZGEBRD", " ", m, m, &c_n1, &c_n1, (ftnlen)
 580                             6, (ftnlen)1);
 581                     maxwrk = max(i__2,i__3);
 582                     if (wntuo || wntuas) {
 583 /* Computing MAX */
 584                         i__2 = maxwrk, i__3 = (*m << 1) + *m * ilaenv_(&c__1,
 585                                 "ZUNGBR", "Q", m, m, m, &c_n1, (ftnlen)6, (
 586                                 ftnlen)1);
 587                         maxwrk = max(i__2,i__3);
 588                     }
 589                     minwrk = *m * 3;
 590                     maxwrk = max(minwrk,maxwrk);
 591                 } else if (wntvo && wntun) {
 592
 593 /*                 Path 2t(N much larger than M, JOBU='N', JOBVT='O') */
 594
 595                     wrkbl = *m + *m * ilaenv_(&c__1, "ZGELQF", " ", m, n, &
 596                             c_n1, &c_n1, (ftnlen)6, (ftnlen)1);
 597 /* Computing MAX */
 598                     i__2 = wrkbl, i__3 = *m + *m * ilaenv_(&c__1, "ZUNGLQ",
 599                             " ", m, n, m, &c_n1, (ftnlen)6, (ftnlen)1);
 600                     wrkbl = max(i__2,i__3);
 601 /* Computing MAX */
 602                     i__2 = wrkbl, i__3 = (*m << 1) + (*m << 1) * ilaenv_(&
 603                             c__1, "ZGEBRD", " ", m, m, &c_n1, &c_n1, (ftnlen)
 604                             6, (ftnlen)1);
 605                     wrkbl = max(i__2,i__3);
 606 /* Computing MAX */
 607                     i__2 = wrkbl, i__3 = (*m << 1) + (*m - 1) * ilaenv_(&c__1,
 608                              "ZUNGBR", "P", m, m, m, &c_n1, (ftnlen)6, (
 609                             ftnlen)1);
 610                     wrkbl = max(i__2,i__3);
 611 /* Computing MAX */
 612                     i__2 = *m * *m + wrkbl, i__3 = *m * *m + *m * *n;
 613                     maxwrk = max(i__2,i__3);
 614                     minwrk = (*m << 1) + *n;
 615                     maxwrk = max(minwrk,maxwrk);
 616                 } else if (wntvo && wntuas) {
 617
 618 /*                 Path 3t(N much larger than M, JOBU='S' or 'A',
 619                    JOBVT='O') */
 620
 621                     wrkbl = *m + *m * ilaenv_(&c__1, "ZGELQF", " ", m, n, &
 622                             c_n1, &c_n1, (ftnlen)6, (ftnlen)1);
 623 /* Computing MAX */
 624                     i__2 = wrkbl, i__3 = *m + *m * ilaenv_(&c__1, "ZUNGLQ",
 625                             " ", m, n, m, &c_n1, (ftnlen)6, (ftnlen)1);
 626                     wrkbl = max(i__2,i__3);
 627 /* Computing MAX */
 628                     i__2 = wrkbl, i__3 = (*m << 1) + (*m << 1) * ilaenv_(&
 629                             c__1, "ZGEBRD", " ", m, m, &c_n1, &c_n1, (ftnlen)
 630                             6, (ftnlen)1);
 631                     wrkbl = max(i__2,i__3);
 632 /* Computing MAX */
 633                     i__2 = wrkbl, i__3 = (*m << 1) + (*m - 1) * ilaenv_(&c__1,
 634                              "ZUNGBR", "P", m, m, m, &c_n1, (ftnlen)6, (
 635                             ftnlen)1);
 636                     wrkbl = max(i__2,i__3);
 637 /* Computing MAX */
 638                     i__2 = wrkbl, i__3 = (*m << 1) + *m * ilaenv_(&c__1,
 639                             "ZUNGBR", "Q", m, m, m, &c_n1, (ftnlen)6, (ftnlen)
 640                             1);
 641                     wrkbl = max(i__2,i__3);
 642 /* Computing MAX */
 643                     i__2 = *m * *m + wrkbl, i__3 = *m * *m + *m * *n;
 644                     maxwrk = max(i__2,i__3);
 645                     minwrk = (*m << 1) + *n;
 646                     maxwrk = max(minwrk,maxwrk);
 647                 } else if (wntvs && wntun) {
 648
 649 /*                 Path 4t(N much larger than M, JOBU='N', JOBVT='S') */
 650
 651                     wrkbl = *m + *m * ilaenv_(&c__1, "ZGELQF", " ", m, n, &
 652                             c_n1, &c_n1, (ftnlen)6, (ftnlen)1);
 653 /* Computing MAX */
 654                     i__2 = wrkbl, i__3 = *m + *m * ilaenv_(&c__1, "ZUNGLQ",
 655                             " ", m, n, m, &c_n1, (ftnlen)6, (ftnlen)1);
 656                     wrkbl = max(i__2,i__3);
 657 /* Computing MAX */
 658                     i__2 = wrkbl, i__3 = (*m << 1) + (*m << 1) * ilaenv_(&
 659                             c__1, "ZGEBRD", " ", m, m, &c_n1, &c_n1, (ftnlen)
 660                             6, (ftnlen)1);
 661                     wrkbl = max(i__2,i__3);
 662 /* Computing MAX */
 663                     i__2 = wrkbl, i__3 = (*m << 1) + (*m - 1) * ilaenv_(&c__1,
 664                              "ZUNGBR", "P", m, m, m, &c_n1, (ftnlen)6, (
 665                             ftnlen)1);
 666                     wrkbl = max(i__2,i__3);
 667                     maxwrk = *m * *m + wrkbl;
 668                     minwrk = (*m << 1) + *n;
 669                     maxwrk = max(minwrk,maxwrk);
 670                 } else if (wntvs && wntuo) {
 671
 672 /*                 Path 5t(N much larger than M, JOBU='O', JOBVT='S') */
 673
 674                     wrkbl = *m + *m * ilaenv_(&c__1, "ZGELQF", " ", m, n, &
 675                             c_n1, &c_n1, (ftnlen)6, (ftnlen)1);
 676 /* Computing MAX */
 677                     i__2 = wrkbl, i__3 = *m + *m * ilaenv_(&c__1, "ZUNGLQ",
 678                             " ", m, n, m, &c_n1, (ftnlen)6, (ftnlen)1);
 679                     wrkbl = max(i__2,i__3);
 680 /* Computing MAX */
 681                     i__2 = wrkbl, i__3 = (*m << 1) + (*m << 1) * ilaenv_(&
 682                             c__1, "ZGEBRD", " ", m, m, &c_n1, &c_n1, (ftnlen)
 683                             6, (ftnlen)1);
 684                     wrkbl = max(i__2,i__3);
 685 /* Computing MAX */
 686                     i__2 = wrkbl, i__3 = (*m << 1) + (*m - 1) * ilaenv_(&c__1,
 687                              "ZUNGBR", "P", m, m, m, &c_n1, (ftnlen)6, (
 688                             ftnlen)1);
 689                     wrkbl = max(i__2,i__3);
 690 /* Computing MAX */
 691                     i__2 = wrkbl, i__3 = (*m << 1) + *m * ilaenv_(&c__1,
 692                             "ZUNGBR", "Q", m, m, m, &c_n1, (ftnlen)6, (ftnlen)
 693                             1);
 694                     wrkbl = max(i__2,i__3);
 695                     maxwrk = (*m << 1) * *m + wrkbl;
 696                     minwrk = (*m << 1) + *n;
 697                     maxwrk = max(minwrk,maxwrk);
 698                 } else if (wntvs && wntuas) {
 699
 700 /*                 Path 6t(N much larger than M, JOBU='S' or 'A',
 701                    JOBVT='S') */
 702
 703                     wrkbl = *m + *m * ilaenv_(&c__1, "ZGELQF", " ", m, n, &
 704                             c_n1, &c_n1, (ftnlen)6, (ftnlen)1);
 705 /* Computing MAX */
 706                     i__2 = wrkbl, i__3 = *m + *m * ilaenv_(&c__1, "ZUNGLQ",
 707                             " ", m, n, m, &c_n1, (ftnlen)6, (ftnlen)1);
 708                     wrkbl = max(i__2,i__3);
 709 /* Computing MAX */
 710                     i__2 = wrkbl, i__3 = (*m << 1) + (*m << 1) * ilaenv_(&
 711                             c__1, "ZGEBRD", " ", m, m, &c_n1, &c_n1, (ftnlen)
 712                             6, (ftnlen)1);
 713                     wrkbl = max(i__2,i__3);
 714 /* Computing MAX */
 715                     i__2 = wrkbl, i__3 = (*m << 1) + (*m - 1) * ilaenv_(&c__1,
 716                              "ZUNGBR", "P", m, m, m, &c_n1, (ftnlen)6, (
 717                             ftnlen)1);
 718                     wrkbl = max(i__2,i__3);
 719 /* Computing MAX */
 720                     i__2 = wrkbl, i__3 = (*m << 1) + *m * ilaenv_(&c__1,
 721                             "ZUNGBR", "Q", m, m, m, &c_n1, (ftnlen)6, (ftnlen)
 722                             1);
 723                     wrkbl = max(i__2,i__3);
 724                     maxwrk = *m * *m + wrkbl;
 725                     minwrk = (*m << 1) + *n;
 726                     maxwrk = max(minwrk,maxwrk);
 727                 } else if (wntva && wntun) {
 728
 729 /*                 Path 7t(N much larger than M, JOBU='N', JOBVT='A') */
 730
 731                     wrkbl = *m + *m * ilaenv_(&c__1, "ZGELQF", " ", m, n, &
 732                             c_n1, &c_n1, (ftnlen)6, (ftnlen)1);
 733 /* Computing MAX */
 734                     i__2 = wrkbl, i__3 = *m + *n * ilaenv_(&c__1, "ZUNGLQ",
 735                             " ", n, n, m, &c_n1, (ftnlen)6, (ftnlen)1);
 736                     wrkbl = max(i__2,i__3);
 737 /* Computing MAX */
 738                     i__2 = wrkbl, i__3 = (*m << 1) + (*m << 1) * ilaenv_(&
 739                             c__1, "ZGEBRD", " ", m, m, &c_n1, &c_n1, (ftnlen)
 740                             6, (ftnlen)1);
 741                     wrkbl = max(i__2,i__3);
 742 /* Computing MAX */
 743                     i__2 = wrkbl, i__3 = (*m << 1) + (*m - 1) * ilaenv_(&c__1,
 744                              "ZUNGBR", "P", m, m, m, &c_n1, (ftnlen)6, (
 745                             ftnlen)1);
 746                     wrkbl = max(i__2,i__3);
 747                     maxwrk = *m * *m + wrkbl;
 748                     minwrk = (*m << 1) + *n;
 749                     maxwrk = max(minwrk,maxwrk);
 750                 } else if (wntva && wntuo) {
 751
 752 /*                 Path 8t(N much larger than M, JOBU='O', JOBVT='A') */
 753
 754                     wrkbl = *m + *m * ilaenv_(&c__1, "ZGELQF", " ", m, n, &
 755                             c_n1, &c_n1, (ftnlen)6, (ftnlen)1);
 756 /* Computing MAX */
 757                     i__2 = wrkbl, i__3 = *m + *n * ilaenv_(&c__1, "ZUNGLQ",
 758                             " ", n, n, m, &c_n1, (ftnlen)6, (ftnlen)1);
 759                     wrkbl = max(i__2,i__3);
 760 /* Computing MAX */
 761                     i__2 = wrkbl, i__3 = (*m << 1) + (*m << 1) * ilaenv_(&
 762                             c__1, "ZGEBRD", " ", m, m, &c_n1, &c_n1, (ftnlen)
 763                             6, (ftnlen)1);
 764                     wrkbl = max(i__2,i__3);
 765 /* Computing MAX */
 766                     i__2 = wrkbl, i__3 = (*m << 1) + (*m - 1) * ilaenv_(&c__1,
 767                              "ZUNGBR", "P", m, m, m, &c_n1, (ftnlen)6, (
 768                             ftnlen)1);
 769                     wrkbl = max(i__2,i__3);
 770 /* Computing MAX */
 771                     i__2 = wrkbl, i__3 = (*m << 1) + *m * ilaenv_(&c__1,
 772                             "ZUNGBR", "Q", m, m, m, &c_n1, (ftnlen)6, (ftnlen)
 773                             1);
 774                     wrkbl = max(i__2,i__3);
 775                     maxwrk = (*m << 1) * *m + wrkbl;
 776                     minwrk = (*m << 1) + *n;
 777                     maxwrk = max(minwrk,maxwrk);
 778                 } else if (wntva && wntuas) {
 779
 780 /*                 Path 9t(N much larger than M, JOBU='S' or 'A',
 781                    JOBVT='A') */
 782
 783                     wrkbl = *m + *m * ilaenv_(&c__1, "ZGELQF", " ", m, n, &
 784                             c_n1, &c_n1, (ftnlen)6, (ftnlen)1);
 785 /* Computing MAX */
 786                     i__2 = wrkbl, i__3 = *m + *n * ilaenv_(&c__1, "ZUNGLQ",
 787                             " ", n, n, m, &c_n1, (ftnlen)6, (ftnlen)1);
 788                     wrkbl = max(i__2,i__3);
 789 /* Computing MAX */
 790                     i__2 = wrkbl, i__3 = (*m << 1) + (*m << 1) * ilaenv_(&
 791                             c__1, "ZGEBRD", " ", m, m, &c_n1, &c_n1, (ftnlen)
 792                             6, (ftnlen)1);
 793                     wrkbl = max(i__2,i__3);
 794 /* Computing MAX */
 795                     i__2 = wrkbl, i__3 = (*m << 1) + (*m - 1) * ilaenv_(&c__1,
 796                              "ZUNGBR", "P", m, m, m, &c_n1, (ftnlen)6, (
 797                             ftnlen)1);
 798                     wrkbl = max(i__2,i__3);
 799 /* Computing MAX */
 800                     i__2 = wrkbl, i__3 = (*m << 1) + *m * ilaenv_(&c__1,
 801                             "ZUNGBR", "Q", m, m, m, &c_n1, (ftnlen)6, (ftnlen)
 802                             1);
 803                     wrkbl = max(i__2,i__3);
 804                     maxwrk = *m * *m + wrkbl;
 805                     minwrk = (*m << 1) + *n;
 806                     maxwrk = max(minwrk,maxwrk);
 807                 }
 808             } else {
 809
 810 /*              Path 10t(N greater than M, but not much larger) */
 811
 812                 maxwrk = (*m << 1) + (*m + *n) * ilaenv_(&c__1, "ZGEBRD",
 813                         " ", m, n, &c_n1, &c_n1, (ftnlen)6, (ftnlen)1);
 814                 if (wntvs || wntvo) {
 815 /* Computing MAX */
 816                     i__2 = maxwrk, i__3 = (*m << 1) + *m * ilaenv_(&c__1,
 817                             "ZUNGBR", "P", m, n, m, &c_n1, (ftnlen)6, (ftnlen)
 818                             1);
 819                     maxwrk = max(i__2,i__3);
 820                 }
 821                 if (wntva) {
 822 /* Computing MAX */
 823                     i__2 = maxwrk, i__3 = (*m << 1) + *n * ilaenv_(&c__1,
 824                             "ZUNGBR", "P", n, n, m, &c_n1, (ftnlen)6, (ftnlen)
 825                             1);
 826                     maxwrk = max(i__2,i__3);
 827                 }
 828                 if (! wntun) {
 829 /* Computing MAX */
 830                     i__2 = maxwrk, i__3 = (*m << 1) + (*m - 1) * ilaenv_(&
 831                             c__1, "ZUNGBR", "Q", m, m, m, &c_n1, (ftnlen)6, (
 832                             ftnlen)1);
 833                     maxwrk = max(i__2,i__3);
 834                 }
 835                 minwrk = (*m << 1) + *n;
 836                 maxwrk = max(minwrk,maxwrk);
 837             }
 838         }
 839         work[1].r = (doublereal) maxwrk, work[1].i = 0.;
 840     }
 841
 842     if (*lwork < minwrk && ! lquery) {
 843         *info = -13;
 844     }
 845     if (*info != 0) {
 846         i__2 = -(*info);
 847         xerbla_("ZGESVD", &i__2);
 848         return 0;
 849     } else if (lquery) {
 850         return 0;
 851     }
 852
 853 /*     Quick return if possible */
 854
 855     if (*m == 0 || *n == 0) {
 856         if (*lwork >= 1) {
 857             work[1].r = 1., work[1].i = 0.;
 858         }
 859         return 0;
 860     }
 861
 862 /*     Get machine constants */
 863
 864     eps = dlamch_("P");
 865     smlnum = sqrt(dlamch_("S")) / eps;
 866     bignum = 1. / smlnum;
 867
 868 /*     Scale A if max element outside range [SMLNUM,BIGNUM] */
 869
 870     anrm = zlange_("M", m, n, &a[a_offset], lda, dum);
 871     iscl = 0;
 872     if (anrm > 0. && anrm < smlnum) {
 873         iscl = 1;
 874         zlascl_("G", &c__0, &c__0, &anrm, &smlnum, m, n, &a[a_offset], lda, &
 875                 ierr);
 876     } else if (anrm > bignum) {
 877         iscl = 1;
 878         zlascl_("G", &c__0, &c__0, &anrm, &bignum, m, n, &a[a_offset], lda, &
 879                 ierr);
 880     }
 881
 882     if (*m >= *n) {
 883
 884 /*        A has at least as many rows as columns. If A has sufficiently
 885           more rows than columns, first reduce using the QR
 886           decomposition (if sufficient workspace available) */
 887
 888         if (*m >= mnthr) {
 889
 890             if (wntun) {
 891
 892 /*              Path 1 (M much larger than N, JOBU='N')
 893                 No left singular vectors to be computed */
 894
 895                 itau = 1;
 896                 iwork = itau + *n;
 897
 898 /*              Compute A=Q*R
 899                 (CWorkspace: need 2*N, prefer N+N*NB)
 900                 (RWorkspace: need 0) */
 901
 902                 i__2 = *lwork - iwork + 1;
 903                 zgeqrf_(m, n, &a[a_offset], lda, &work[itau], &work[iwork], &
 904                         i__2, &ierr);
 905
 906 /*              Zero out below R */
 907
 908                 i__2 = *n - 1;
 909                 i__3 = *n - 1;
 910                 zlaset_("L", &i__2, &i__3, &c_b1, &c_b1, &a_ref(2, 1), lda);
 911                 ie = 1;
 912                 itauq = 1;
 913                 itaup = itauq + *n;
 914                 iwork = itaup + *n;
 915
 916 /*              Bidiagonalize R in A
 917                 (CWorkspace: need 3*N, prefer 2*N+2*N*NB)
 918                 (RWorkspace: need N) */
 919
 920                 i__2 = *lwork - iwork + 1;
 921                 zgebrd_(n, n, &a[a_offset], lda, &s[1], &rwork[ie], &work[
 922                         itauq], &work[itaup], &work[iwork], &i__2, &ierr);
 923                 ncvt = 0;
 924                 if (wntvo || wntvas) {
 925
 926 /*                 If right singular vectors desired, generate P'.
 927                    (CWorkspace: need 3*N-1, prefer 2*N+(N-1)*NB)
 928                    (RWorkspace: 0) */
 929
 930                     i__2 = *lwork - iwork + 1;
 931                     zungbr_("P", n, n, n, &a[a_offset], lda, &work[itaup], &
 932                             work[iwork], &i__2, &ierr);
 933                     ncvt = *n;
 934                 }
 935                 irwork = ie + *n;
 936
 937 /*              Perform bidiagonal QR iteration, computing right
 938                 singular vectors of A in A if desired
 939                 (CWorkspace: 0)
 940                 (RWorkspace: need BDSPAC) */
 941
 942                 zbdsqr_("U", n, &ncvt, &c__0, &c__0, &s[1], &rwork[ie], &a[
 943                         a_offset], lda, cdum, &c__1, cdum, &c__1, &rwork[
 944                         irwork], info);
 945
 946 /*              If right singular vectors desired in VT, copy them there */
 947
 948                 if (wntvas) {
 949                     zlacpy_("F", n, n, &a[a_offset], lda, &vt[vt_offset],
 950                             ldvt);
 951                 }
 952
 953             } else if (wntuo && wntvn) {
 954
 955 /*              Path 2 (M much larger than N, JOBU='O', JOBVT='N')
 956                 N left singular vectors to be overwritten on A and
 957                 no right singular vectors to be computed */
 958
 959                 if (*lwork >= *n * *n + *n * 3) {
 960
 961 /*                 Sufficient workspace for a fast algorithm */
 962
 963                     ir = 1;
 964 /* Computing MAX */
 965                     i__2 = wrkbl, i__3 = *lda * *n;
 966                     if (*lwork >= max(i__2,i__3) + *lda * *n) {
 967
 968 /*                    WORK(IU) is LDA by N, WORK(IR) is LDA by N */
 969
 970                         ldwrku = *lda;
 971                         ldwrkr = *lda;
 972                     } else /* if(complicated condition) */ {
 973 /* Computing MAX */
 974                         i__2 = wrkbl, i__3 = *lda * *n;
 975                         if (*lwork >= max(i__2,i__3) + *n * *n) {
 976
 977 /*                    WORK(IU) is LDA by N, WORK(IR) is N by N */
 978
 979                             ldwrku = *lda;
 980                             ldwrkr = *n;
 981                         } else {
 982
 983 /*                    WORK(IU) is LDWRKU by N, WORK(IR) is N by N */
 984
 985                             ldwrku = (*lwork - *n * *n) / *n;
 986                             ldwrkr = *n;
 987                         }
 988                     }
 989                     itau = ir + ldwrkr * *n;
 990                     iwork = itau + *n;
 991
 992 /*                 Compute A=Q*R
 993                    (CWorkspace: need N*N+2*N, prefer N*N+N+N*NB)
 994                    (RWorkspace: 0) */
 995
 996                     i__2 = *lwork - iwork + 1;
 997                     zgeqrf_(m, n, &a[a_offset], lda, &work[itau], &work[iwork]
 998                             , &i__2, &ierr);
 999
1000 /*                 Copy R to WORK(IR) and zero out below it */
1001
1002                     zlacpy_("U", n, n, &a[a_offset], lda, &work[ir], &ldwrkr);
1003                     i__2 = *n - 1;
1004                     i__3 = *n - 1;
1005                     zlaset_("L", &i__2, &i__3, &c_b1, &c_b1, &work[ir + 1], &
1006                             ldwrkr);
1007
1008 /*                 Generate Q in A
1009                    (CWorkspace: need N*N+2*N, prefer N*N+N+N*NB)
1010                    (RWorkspace: 0) */
1011
1012                     i__2 = *lwork - iwork + 1;
1013                     zungqr_(m, n, n, &a[a_offset], lda, &work[itau], &work[
1014                             iwork], &i__2, &ierr);
1015                     ie = 1;
1016                     itauq = itau;
1017                     itaup = itauq + *n;
1018                     iwork = itaup + *n;
1019
1020 /*                 Bidiagonalize R in WORK(IR)
1021                    (CWorkspace: need N*N+3*N, prefer N*N+2*N+2*N*NB)
1022                    (RWorkspace: need N) */
1023
1024                     i__2 = *lwork - iwork + 1;
1025                     zgebrd_(n, n, &work[ir], &ldwrkr, &s[1], &rwork[ie], &
1026                             work[itauq], &work[itaup], &work[iwork], &i__2, &
1027                             ierr);
1028
1029 /*                 Generate left vectors bidiagonalizing R
1030                    (CWorkspace: need N*N+3*N, prefer N*N+2*N+N*NB)
1031                    (RWorkspace: need 0) */
1032
1033                     i__2 = *lwork - iwork + 1;
1034                     zungbr_("Q", n, n, n, &work[ir], &ldwrkr, &work[itauq], &
1035                             work[iwork], &i__2, &ierr);
1036                     irwork = ie + *n;
1037
1038 /*                 Perform bidiagonal QR iteration, computing left
1039                    singular vectors of R in WORK(IR)
1040                    (CWorkspace: need N*N)
1041                    (RWorkspace: need BDSPAC) */
1042
1043                     zbdsqr_("U", n, &c__0, n, &c__0, &s[1], &rwork[ie], cdum,
1044                             &c__1, &work[ir], &ldwrkr, cdum, &c__1, &rwork[
1045                             irwork], info);
1046                     iu = itauq;
1047
1048 /*                 Multiply Q in A by left singular vectors of R in
1049                    WORK(IR), storing result in WORK(IU) and copying to A
1050                    (CWorkspace: need N*N+N, prefer N*N+M*N)
1051                    (RWorkspace: 0) */
1052
1053                     i__2 = *m;
1054                     i__3 = ldwrku;
1055                     for (i__ = 1; i__3 < 0 ? i__ >= i__2 : i__ <= i__2; i__ +=
1056                              i__3) {
1057 /* Computing MIN */
1058                         i__4 = *m - i__ + 1;
1059                         chunk = min(i__4,ldwrku);
1060                         zgemm_("N", "N", &chunk, n, n, &c_b2, &a_ref(i__, 1),
1061                                 lda, &work[ir], &ldwrkr, &c_b1, &work[iu], &
1062                                 ldwrku);
1063                         zlacpy_("F", &chunk, n, &work[iu], &ldwrku, &a_ref(
1064                                 i__, 1), lda);
1065 /* L10: */
1066                     }
1067
1068                 } else {
1069
1070 /*                 Insufficient workspace for a fast algorithm */
1071
1072                     ie = 1;
1073                     itauq = 1;
1074                     itaup = itauq + *n;
1075                     iwork = itaup + *n;
1076
1077 /*                 Bidiagonalize A
1078                    (CWorkspace: need 2*N+M, prefer 2*N+(M+N)*NB)
1079                    (RWorkspace: N) */
1080
1081                     i__3 = *lwork - iwork + 1;
1082                     zgebrd_(m, n, &a[a_offset], lda, &s[1], &rwork[ie], &work[
1083                             itauq], &work[itaup], &work[iwork], &i__3, &ierr);
1084
1085 /*                 Generate left vectors bidiagonalizing A
1086                    (CWorkspace: need 3*N, prefer 2*N+N*NB)
1087                    (RWorkspace: 0) */
1088
1089                     i__3 = *lwork - iwork + 1;
1090                     zungbr_("Q", m, n, n, &a[a_offset], lda, &work[itauq], &
1091                             work[iwork], &i__3, &ierr);
1092                     irwork = ie + *n;
1093
1094 /*                 Perform bidiagonal QR iteration, computing left
1095                    singular vectors of A in A
1096                    (CWorkspace: need 0)
1097                    (RWorkspace: need BDSPAC) */
1098
1099                     zbdsqr_("U", n, &c__0, m, &c__0, &s[1], &rwork[ie], cdum,
1100                             &c__1, &a[a_offset], lda, cdum, &c__1, &rwork[
1101                             irwork], info);
1102
1103                 }
1104
1105             } else if (wntuo && wntvas) {
1106
1107 /*              Path 3 (M much larger than N, JOBU='O', JOBVT='S' or 'A')
1108                 N left singular vectors to be overwritten on A and
1109                 N right singular vectors to be computed in VT */
1110
1111                 if (*lwork >= *n * *n + *n * 3) {
1112
1113 /*                 Sufficient workspace for a fast algorithm */
1114
1115                     ir = 1;
1116 /* Computing MAX */
1117                     i__3 = wrkbl, i__2 = *lda * *n;
1118                     if (*lwork >= max(i__3,i__2) + *lda * *n) {
1119
1120 /*                    WORK(IU) is LDA by N and WORK(IR) is LDA by N */
1121
1122                         ldwrku = *lda;
1123                         ldwrkr = *lda;
1124                     } else /* if(complicated condition) */ {
1125 /* Computing MAX */
1126                         i__3 = wrkbl, i__2 = *lda * *n;
1127                         if (*lwork >= max(i__3,i__2) + *n * *n) {
1128
1129 /*                    WORK(IU) is LDA by N and WORK(IR) is N by N */
1130
1131                             ldwrku = *lda;
1132                             ldwrkr = *n;
1133                         } else {
1134
1135 /*                    WORK(IU) is LDWRKU by N and WORK(IR) is N by N */
1136
1137                             ldwrku = (*lwork - *n * *n) / *n;
1138                             ldwrkr = *n;
1139                         }
1140                     }
1141                     itau = ir + ldwrkr * *n;
1142                     iwork = itau + *n;
1143
1144 /*                 Compute A=Q*R
1145                    (CWorkspace: need N*N+2*N, prefer N*N+N+N*NB)
1146                    (RWorkspace: 0) */
1147
1148                     i__3 = *lwork - iwork + 1;
1149                     zgeqrf_(m, n, &a[a_offset], lda, &work[itau], &work[iwork]
1150                             , &i__3, &ierr);
1151
1152 /*                 Copy R to VT, zeroing out below it */
1153
1154                     zlacpy_("U", n, n, &a[a_offset], lda, &vt[vt_offset],
1155                             ldvt);
1156                     i__3 = *n - 1;
1157                     i__2 = *n - 1;
1158                     zlaset_("L", &i__3, &i__2, &c_b1, &c_b1, &vt_ref(2, 1),
1159                             ldvt);
1160
1161 /*                 Generate Q in A
1162                    (CWorkspace: need N*N+2*N, prefer N*N+N+N*NB)
1163                    (RWorkspace: 0) */
1164
1165                     i__3 = *lwork - iwork + 1;
1166                     zungqr_(m, n, n, &a[a_offset], lda, &work[itau], &work[
1167                             iwork], &i__3, &ierr);
1168                     ie = 1;
1169                     itauq = itau;
1170                     itaup = itauq + *n;
1171                     iwork = itaup + *n;
1172
1173 /*                 Bidiagonalize R in VT, copying result to WORK(IR)
1174                    (CWorkspace: need N*N+3*N, prefer N*N+2*N+2*N*NB)
1175                    (RWorkspace: need N) */
1176
1177                     i__3 = *lwork - iwork + 1;
1178                     zgebrd_(n, n, &vt[vt_offset], ldvt, &s[1], &rwork[ie], &
1179                             work[itauq], &work[itaup], &work[iwork], &i__3, &
1180                             ierr);
1181                     zlacpy_("L", n, n, &vt[vt_offset], ldvt, &work[ir], &
1182                             ldwrkr);
1183
1184 /*                 Generate left vectors bidiagonalizing R in WORK(IR)
1185                    (CWorkspace: need N*N+3*N, prefer N*N+2*N+N*NB)
1186                    (RWorkspace: 0) */
1187
1188                     i__3 = *lwork - iwork + 1;
1189                     zungbr_("Q", n, n, n, &work[ir], &ldwrkr, &work[itauq], &
1190                             work[iwork], &i__3, &ierr);
1191
1192 /*                 Generate right vectors bidiagonalizing R in VT
1193                    (CWorkspace: need N*N+3*N-1, prefer N*N+2*N+(N-1)*NB)
1194                    (RWorkspace: 0) */
1195
1196                     i__3 = *lwork - iwork + 1;
1197                     zungbr_("P", n, n, n, &vt[vt_offset], ldvt, &work[itaup],
1198                             &work[iwork], &i__3, &ierr);
1199                     irwork = ie + *n;
1200
1201 /*                 Perform bidiagonal QR iteration, computing left
1202                    singular vectors of R in WORK(IR) and computing right
1203                    singular vectors of R in VT
1204                    (CWorkspace: need N*N)
1205                    (RWorkspace: need BDSPAC) */
1206
1207                     zbdsqr_("U", n, n, n, &c__0, &s[1], &rwork[ie], &vt[
1208                             vt_offset], ldvt, &work[ir], &ldwrkr, cdum, &c__1,
1209                              &rwork[irwork], info);
1210                     iu = itauq;
1211
1212 /*                 Multiply Q in A by left singular vectors of R in
1213                    WORK(IR), storing result in WORK(IU) and copying to A
1214                    (CWorkspace: need N*N+N, prefer N*N+M*N)
1215                    (RWorkspace: 0) */
1216
1217                     i__3 = *m;
1218                     i__2 = ldwrku;
1219                     for (i__ = 1; i__2 < 0 ? i__ >= i__3 : i__ <= i__3; i__ +=
1220                              i__2) {
1221 /* Computing MIN */
1222                         i__4 = *m - i__ + 1;
1223                         chunk = min(i__4,ldwrku);
1224                         zgemm_("N", "N", &chunk, n, n, &c_b2, &a_ref(i__, 1),
1225                                 lda, &work[ir], &ldwrkr, &c_b1, &work[iu], &
1226                                 ldwrku);
1227                         zlacpy_("F", &chunk, n, &work[iu], &ldwrku, &a_ref(
1228                                 i__, 1), lda);
1229 /* L20: */
1230                     }
1231
1232                 } else {
1233
1234 /*                 Insufficient workspace for a fast algorithm */
1235
1236                     itau = 1;
1237                     iwork = itau + *n;
1238
1239 /*                 Compute A=Q*R
1240                    (CWorkspace: need 2*N, prefer N+N*NB)
1241                    (RWorkspace: 0) */
1242
1243                     i__2 = *lwork - iwork + 1;
1244                     zgeqrf_(m, n, &a[a_offset], lda, &work[itau], &work[iwork]
1245                             , &i__2, &ierr);
1246
1247 /*                 Copy R to VT, zeroing out below it */
1248
1249                     zlacpy_("U", n, n, &a[a_offset], lda, &vt[vt_offset],
1250                             ldvt);
1251                     i__2 = *n - 1;
1252                     i__3 = *n - 1;
1253                     zlaset_("L", &i__2, &i__3, &c_b1, &c_b1, &vt_ref(2, 1),
1254                             ldvt);
1255
1256 /*                 Generate Q in A
1257                    (CWorkspace: need 2*N, prefer N+N*NB)
1258                    (RWorkspace: 0) */
1259
1260                     i__2 = *lwork - iwork + 1;
1261                     zungqr_(m, n, n, &a[a_offset], lda, &work[itau], &work[
1262                             iwork], &i__2, &ierr);
1263                     ie = 1;
1264                     itauq = itau;
1265                     itaup = itauq + *n;
1266                     iwork = itaup + *n;
1267
1268 /*                 Bidiagonalize R in VT
1269                    (CWorkspace: need 3*N, prefer 2*N+2*N*NB)
1270                    (RWorkspace: N) */
1271
1272                     i__2 = *lwork - iwork + 1;
1273                     zgebrd_(n, n, &vt[vt_offset], ldvt, &s[1], &rwork[ie], &
1274                             work[itauq], &work[itaup], &work[iwork], &i__2, &
1275                             ierr);
1276
1277 /*                 Multiply Q in A by left vectors bidiagonalizing R
1278                    (CWorkspace: need 2*N+M, prefer 2*N+M*NB)
1279                    (RWorkspace: 0) */
1280
1281                     i__2 = *lwork - iwork + 1;
1282                     zunmbr_("Q", "R", "N", m, n, n, &vt[vt_offset], ldvt, &
1283                             work[itauq], &a[a_offset], lda, &work[iwork], &
1284                             i__2, &ierr);
1285
1286 /*                 Generate right vectors bidiagonalizing R in VT
1287                    (CWorkspace: need 3*N-1, prefer 2*N+(N-1)*NB)
1288                    (RWorkspace: 0) */
1289
1290                     i__2 = *lwork - iwork + 1;
1291                     zungbr_("P", n, n, n, &vt[vt_offset], ldvt, &work[itaup],
1292                             &work[iwork], &i__2, &ierr);
1293                     irwork = ie + *n;
1294
1295 /*                 Perform bidiagonal QR iteration, computing left
1296                    singular vectors of A in A and computing right
1297                    singular vectors of A in VT
1298                    (CWorkspace: 0)
1299                    (RWorkspace: need BDSPAC) */
1300
1301                     zbdsqr_("U", n, n, m, &c__0, &s[1], &rwork[ie], &vt[
1302                             vt_offset], ldvt, &a[a_offset], lda, cdum, &c__1,
1303                             &rwork[irwork], info);
1304
1305                 }
1306
1307             } else if (wntus) {
1308
1309                 if (wntvn) {
1310
1311 /*                 Path 4 (M much larger than N, JOBU='S', JOBVT='N')
1312                    N left singular vectors to be computed in U and
1313                    no right singular vectors to be computed */
1314
1315                     if (*lwork >= *n * *n + *n * 3) {
1316
1317 /*                    Sufficient workspace for a fast algorithm */
1318
1319                         ir = 1;
1320                         if (*lwork >= wrkbl + *lda * *n) {
1321
1322 /*                       WORK(IR) is LDA by N */
1323
1324                             ldwrkr = *lda;
1325                         } else {
1326
1327 /*                       WORK(IR) is N by N */
1328
1329                             ldwrkr = *n;
1330                         }
1331                         itau = ir + ldwrkr * *n;
1332                         iwork = itau + *n;
1333
1334 /*                    Compute A=Q*R
1335                       (CWorkspace: need N*N+2*N, prefer N*N+N+N*NB)
1336                       (RWorkspace: 0) */
1337
1338                         i__2 = *lwork - iwork + 1;
1339                         zgeqrf_(m, n, &a[a_offset], lda, &work[itau], &work[
1340                                 iwork], &i__2, &ierr);
1341
1342 /*                    Copy R to WORK(IR), zeroing out below it */
1343
1344                         zlacpy_("U", n, n, &a[a_offset], lda, &work[ir], &
1345                                 ldwrkr);
1346                         i__2 = *n - 1;
1347                         i__3 = *n - 1;
1348                         zlaset_("L", &i__2, &i__3, &c_b1, &c_b1, &work[ir + 1]
1349                                 , &ldwrkr);
1350
1351 /*                    Generate Q in A
1352                       (CWorkspace: need N*N+2*N, prefer N*N+N+N*NB)
1353                       (RWorkspace: 0) */
1354
1355                         i__2 = *lwork - iwork + 1;
1356                         zungqr_(m, n, n, &a[a_offset], lda, &work[itau], &
1357                                 work[iwork], &i__2, &ierr);
1358                         ie = 1;
1359                         itauq = itau;
1360                         itaup = itauq + *n;
1361                         iwork = itaup + *n;
1362
1363 /*                    Bidiagonalize R in WORK(IR)
1364                       (CWorkspace: need N*N+3*N, prefer N*N+2*N+2*N*NB)
1365                       (RWorkspace: need N) */
1366
1367                         i__2 = *lwork - iwork + 1;
1368                         zgebrd_(n, n, &work[ir], &ldwrkr, &s[1], &rwork[ie], &
1369                                 work[itauq], &work[itaup], &work[iwork], &
1370                                 i__2, &ierr);
1371
1372 /*                    Generate left vectors bidiagonalizing R in WORK(IR)
1373                       (CWorkspace: need N*N+3*N, prefer N*N+2*N+N*NB)
1374                       (RWorkspace: 0) */
1375
1376                         i__2 = *lwork - iwork + 1;
1377                         zungbr_("Q", n, n, n, &work[ir], &ldwrkr, &work[itauq]
1378                                 , &work[iwork], &i__2, &ierr);
1379                         irwork = ie + *n;
1380
1381 /*                    Perform bidiagonal QR iteration, computing left
1382                       singular vectors of R in WORK(IR)
1383                       (CWorkspace: need N*N)
1384                       (RWorkspace: need BDSPAC) */
1385
1386                         zbdsqr_("U", n, &c__0, n, &c__0, &s[1], &rwork[ie],
1387                                 cdum, &c__1, &work[ir], &ldwrkr, cdum, &c__1,
1388                                 &rwork[irwork], info);
1389
1390 /*                    Multiply Q in A by left singular vectors of R in
1391                       WORK(IR), storing result in U
1392                       (CWorkspace: need N*N)
1393                       (RWorkspace: 0) */
1394
1395                         zgemm_("N", "N", m, n, n, &c_b2, &a[a_offset], lda, &
1396                                 work[ir], &ldwrkr, &c_b1, &u[u_offset], ldu);
1397
1398                     } else {
1399
1400 /*                    Insufficient workspace for a fast algorithm */
1401
1402                         itau = 1;
1403                         iwork = itau + *n;
1404
1405 /*                    Compute A=Q*R, copying result to U
1406                       (CWorkspace: need 2*N, prefer N+N*NB)
1407                       (RWorkspace: 0) */
1408
1409                         i__2 = *lwork - iwork + 1;
1410                         zgeqrf_(m, n, &a[a_offset], lda, &work[itau], &work[
1411                                 iwork], &i__2, &ierr);
1412                         zlacpy_("L", m, n, &a[a_offset], lda, &u[u_offset],
1413                                 ldu);
1414
1415 /*                    Generate Q in U
1416                       (CWorkspace: need 2*N, prefer N+N*NB)
1417                       (RWorkspace: 0) */
1418
1419                         i__2 = *lwork - iwork + 1;
1420                         zungqr_(m, n, n, &u[u_offset], ldu, &work[itau], &
1421                                 work[iwork], &i__2, &ierr);
1422                         ie = 1;
1423                         itauq = itau;
1424                         itaup = itauq + *n;
1425                         iwork = itaup + *n;
1426
1427 /*                    Zero out below R in A */
1428
1429                         i__2 = *n - 1;
1430                         i__3 = *n - 1;
1431                         zlaset_("L", &i__2, &i__3, &c_b1, &c_b1, &a_ref(2, 1),
1432                                  lda);
1433
1434 /*                    Bidiagonalize R in A
1435                       (CWorkspace: need 3*N, prefer 2*N+2*N*NB)
1436                       (RWorkspace: need N) */
1437
1438                         i__2 = *lwork - iwork + 1;
1439                         zgebrd_(n, n, &a[a_offset], lda, &s[1], &rwork[ie], &
1440                                 work[itauq], &work[itaup], &work[iwork], &
1441                                 i__2, &ierr);
1442
1443 /*                    Multiply Q in U by left vectors bidiagonalizing R
1444                       (CWorkspace: need 2*N+M, prefer 2*N+M*NB)
1445                       (RWorkspace: 0) */
1446
1447                         i__2 = *lwork - iwork + 1;
1448                         zunmbr_("Q", "R", "N", m, n, n, &a[a_offset], lda, &
1449                                 work[itauq], &u[u_offset], ldu, &work[iwork],
1450                                 &i__2, &ierr)
1451                                 ;
1452                         irwork = ie + *n;
1453
1454 /*                    Perform bidiagonal QR iteration, computing left
1455                       singular vectors of A in U
1456                       (CWorkspace: 0)
1457                       (RWorkspace: need BDSPAC) */
1458
1459                         zbdsqr_("U", n, &c__0, m, &c__0, &s[1], &rwork[ie],
1460                                 cdum, &c__1, &u[u_offset], ldu, cdum, &c__1, &
1461                                 rwork[irwork], info);
1462
1463                     }
1464
1465                 } else if (wntvo) {
1466
1467 /*                 Path 5 (M much larger than N, JOBU='S', JOBVT='O')
1468                    N left singular vectors to be computed in U and
1469                    N right singular vectors to be overwritten on A */
1470
1471                     if (*lwork >= (*n << 1) * *n + *n * 3) {
1472
1473 /*                    Sufficient workspace for a fast algorithm */
1474
1475                         iu = 1;
1476                         if (*lwork >= wrkbl + (*lda << 1) * *n) {
1477
1478 /*                       WORK(IU) is LDA by N and WORK(IR) is LDA by N */
1479
1480                             ldwrku = *lda;
1481                             ir = iu + ldwrku * *n;
1482                             ldwrkr = *lda;
1483                         } else if (*lwork >= wrkbl + (*lda + *n) * *n) {
1484
1485 /*                       WORK(IU) is LDA by N and WORK(IR) is N by N */
1486
1487                             ldwrku = *lda;
1488                             ir = iu + ldwrku * *n;
1489                             ldwrkr = *n;
1490                         } else {
1491
1492 /*                       WORK(IU) is N by N and WORK(IR) is N by N */
1493
1494                             ldwrku = *n;
1495                             ir = iu + ldwrku * *n;
1496                             ldwrkr = *n;
1497                         }
1498                         itau = ir + ldwrkr * *n;
1499                         iwork = itau + *n;
1500
1501 /*                    Compute A=Q*R
1502                       (CWorkspace: need 2*N*N+2*N, prefer 2*N*N+N+N*NB)
1503                       (RWorkspace: 0) */
1504
1505                         i__2 = *lwork - iwork + 1;
1506                         zgeqrf_(m, n, &a[a_offset], lda, &work[itau], &work[
1507                                 iwork], &i__2, &ierr);
1508
1509 /*                    Copy R to WORK(IU), zeroing out below it */
1510
1511                         zlacpy_("U", n, n, &a[a_offset], lda, &work[iu], &
1512                                 ldwrku);
1513                         i__2 = *n - 1;
1514                         i__3 = *n - 1;
1515                         zlaset_("L", &i__2, &i__3, &c_b1, &c_b1, &work[iu + 1]
1516                                 , &ldwrku);
1517
1518 /*                    Generate Q in A
1519                       (CWorkspace: need 2*N*N+2*N, prefer 2*N*N+N+N*NB)
1520                       (RWorkspace: 0) */
1521
1522                         i__2 = *lwork - iwork + 1;
1523                         zungqr_(m, n, n, &a[a_offset], lda, &work[itau], &
1524                                 work[iwork], &i__2, &ierr);
1525                         ie = 1;
1526                         itauq = itau;
1527                         itaup = itauq + *n;
1528                         iwork = itaup + *n;
1529
1530 /*                    Bidiagonalize R in WORK(IU), copying result to
1531                       WORK(IR)
1532                       (CWorkspace: need   2*N*N+3*N,
1533                                    prefer 2*N*N+2*N+2*N*NB)
1534                       (RWorkspace: need   N) */
1535
1536                         i__2 = *lwork - iwork + 1;
1537                         zgebrd_(n, n, &work[iu], &ldwrku, &s[1], &rwork[ie], &
1538                                 work[itauq], &work[itaup], &work[iwork], &
1539                                 i__2, &ierr);
1540                         zlacpy_("U", n, n, &work[iu], &ldwrku, &work[ir], &
1541                                 ldwrkr);
1542
1543 /*                    Generate left bidiagonalizing vectors in WORK(IU)
1544                       (CWorkspace: need 2*N*N+3*N, prefer 2*N*N+2*N+N*NB)
1545                       (RWorkspace: 0) */
1546
1547                         i__2 = *lwork - iwork + 1;
1548                         zungbr_("Q", n, n, n, &work[iu], &ldwrku, &work[itauq]
1549                                 , &work[iwork], &i__2, &ierr);
1550
1551 /*                    Generate right bidiagonalizing vectors in WORK(IR)
1552                       (CWorkspace: need   2*N*N+3*N-1,
1553                                    prefer 2*N*N+2*N+(N-1)*NB)
1554                       (RWorkspace: 0) */
1555
1556                         i__2 = *lwork - iwork + 1;
1557                         zungbr_("P", n, n, n, &work[ir], &ldwrkr, &work[itaup]
1558                                 , &work[iwork], &i__2, &ierr);
1559                         irwork = ie + *n;
1560
1561 /*                    Perform bidiagonal QR iteration, computing left
1562                       singular vectors of R in WORK(IU) and computing
1563                       right singular vectors of R in WORK(IR)
1564                       (CWorkspace: need 2*N*N)
1565                       (RWorkspace: need BDSPAC) */
1566
1567                         zbdsqr_("U", n, n, n, &c__0, &s[1], &rwork[ie], &work[
1568                                 ir], &ldwrkr, &work[iu], &ldwrku, cdum, &c__1,
1569                                  &rwork[irwork], info);
1570
1571 /*                    Multiply Q in A by left singular vectors of R in
1572                       WORK(IU), storing result in U
1573                       (CWorkspace: need N*N)
1574                       (RWorkspace: 0) */
1575
1576                         zgemm_("N", "N", m, n, n, &c_b2, &a[a_offset], lda, &
1577                                 work[iu], &ldwrku, &c_b1, &u[u_offset], ldu);
1578
1579 /*                    Copy right singular vectors of R to A
1580                       (CWorkspace: need N*N)
1581                       (RWorkspace: 0) */
1582
1583                         zlacpy_("F", n, n, &work[ir], &ldwrkr, &a[a_offset],
1584                                 lda);
1585
1586                     } else {
1587
1588 /*                    Insufficient workspace for a fast algorithm */
1589
1590                         itau = 1;
1591                         iwork = itau + *n;
1592
1593 /*                    Compute A=Q*R, copying result to U
1594                       (CWorkspace: need 2*N, prefer N+N*NB)
1595                       (RWorkspace: 0) */
1596
1597                         i__2 = *lwork - iwork + 1;
1598                         zgeqrf_(m, n, &a[a_offset], lda, &work[itau], &work[
1599                                 iwork], &i__2, &ierr);
1600                         zlacpy_("L", m, n, &a[a_offset], lda, &u[u_offset],
1601                                 ldu);
1602
1603 /*                    Generate Q in U
1604                       (CWorkspace: need 2*N, prefer N+N*NB)
1605                       (RWorkspace: 0) */
1606
1607                         i__2 = *lwork - iwork + 1;
1608                         zungqr_(m, n, n, &u[u_offset], ldu, &work[itau], &
1609                                 work[iwork], &i__2, &ierr);
1610                         ie = 1;
1611                         itauq = itau;
1612                         itaup = itauq + *n;
1613                         iwork = itaup + *n;
1614
1615 /*                    Zero out below R in A */
1616
1617                         i__2 = *n - 1;
1618                         i__3 = *n - 1;
1619                         zlaset_("L", &i__2, &i__3, &c_b1, &c_b1, &a_ref(2, 1),
1620                                  lda);
1621
1622 /*                    Bidiagonalize R in A
1623                       (CWorkspace: need 3*N, prefer 2*N+2*N*NB)
1624                       (RWorkspace: need N) */
1625
1626                         i__2 = *lwork - iwork + 1;
1627                         zgebrd_(n, n, &a[a_offset], lda, &s[1], &rwork[ie], &
1628                                 work[itauq], &work[itaup], &work[iwork], &
1629                                 i__2, &ierr);
1630
1631 /*                    Multiply Q in U by left vectors bidiagonalizing R
1632                       (CWorkspace: need 2*N+M, prefer 2*N+M*NB)
1633                       (RWorkspace: 0) */
1634
1635                         i__2 = *lwork - iwork + 1;
1636                         zunmbr_("Q", "R", "N", m, n, n, &a[a_offset], lda, &
1637                                 work[itauq], &u[u_offset], ldu, &work[iwork],
1638                                 &i__2, &ierr)
1639                                 ;
1640
1641 /*                    Generate right vectors bidiagonalizing R in A
1642                       (CWorkspace: need 3*N-1, prefer 2*N+(N-1)*NB)
1643                       (RWorkspace: 0) */
1644
1645                         i__2 = *lwork - iwork + 1;
1646                         zungbr_("P", n, n, n, &a[a_offset], lda, &work[itaup],
1647                                  &work[iwork], &i__2, &ierr);
1648                         irwork = ie + *n;
1649
1650 /*                    Perform bidiagonal QR iteration, computing left
1651                       singular vectors of A in U and computing right
1652                       singular vectors of A in A
1653                       (CWorkspace: 0)
1654                       (RWorkspace: need BDSPAC) */
1655
1656                         zbdsqr_("U", n, n, m, &c__0, &s[1], &rwork[ie], &a[
1657                                 a_offset], lda, &u[u_offset], ldu, cdum, &
1658                                 c__1, &rwork[irwork], info);
1659
1660                     }
1661
1662                 } else if (wntvas) {
1663
1664 /*                 Path 6 (M much larger than N, JOBU='S', JOBVT='S'
1665                            or 'A')
1666                    N left singular vectors to be computed in U and
1667                    N right singular vectors to be computed in VT */
1668
1669                     if (*lwork >= *n * *n + *n * 3) {
1670
1671 /*                    Sufficient workspace for a fast algorithm */
1672
1673                         iu = 1;
1674                         if (*lwork >= wrkbl + *lda * *n) {
1675
1676 /*                       WORK(IU) is LDA by N */
1677
1678                             ldwrku = *lda;
1679                         } else {
1680
1681 /*                       WORK(IU) is N by N */
1682
1683                             ldwrku = *n;
1684                         }
1685                         itau = iu + ldwrku * *n;
1686                         iwork = itau + *n;
1687
1688 /*                    Compute A=Q*R
1689                       (CWorkspace: need N*N+2*N, prefer N*N+N+N*NB)
1690                       (RWorkspace: 0) */
1691
1692                         i__2 = *lwork - iwork + 1;
1693                         zgeqrf_(m, n, &a[a_offset], lda, &work[itau], &work[
1694                                 iwork], &i__2, &ierr);
1695
1696 /*                    Copy R to WORK(IU), zeroing out below it */
1697
1698                         zlacpy_("U", n, n, &a[a_offset], lda, &work[iu], &
1699                                 ldwrku);
1700                         i__2 = *n - 1;
1701                         i__3 = *n - 1;
1702                         zlaset_("L", &i__2, &i__3, &c_b1, &c_b1, &work[iu + 1]
1703                                 , &ldwrku);
1704
1705 /*                    Generate Q in A
1706                       (CWorkspace: need N*N+2*N, prefer N*N+N+N*NB)
1707                       (RWorkspace: 0) */
1708
1709                         i__2 = *lwork - iwork + 1;
1710                         zungqr_(m, n, n, &a[a_offset], lda, &work[itau], &
1711                                 work[iwork], &i__2, &ierr);
1712                         ie = 1;
1713                         itauq = itau;
1714                         itaup = itauq + *n;
1715                         iwork = itaup + *n;
1716
1717 /*                    Bidiagonalize R in WORK(IU), copying result to VT
1718                       (CWorkspace: need N*N+3*N, prefer N*N+2*N+2*N*NB)
1719                       (RWorkspace: need N) */
1720
1721                         i__2 = *lwork - iwork + 1;
1722                         zgebrd_(n, n, &work[iu], &ldwrku, &s[1], &rwork[ie], &
1723                                 work[itauq], &work[itaup], &work[iwork], &
1724                                 i__2, &ierr);
1725                         zlacpy_("U", n, n, &work[iu], &ldwrku, &vt[vt_offset],
1726                                  ldvt);
1727
1728 /*                    Generate left bidiagonalizing vectors in WORK(IU)
1729                       (CWorkspace: need N*N+3*N, prefer N*N+2*N+N*NB)
1730                       (RWorkspace: 0) */
1731
1732                         i__2 = *lwork - iwork + 1;
1733                         zungbr_("Q", n, n, n, &work[iu], &ldwrku, &work[itauq]
1734                                 , &work[iwork], &i__2, &ierr);
1735
1736 /*                    Generate right bidiagonalizing vectors in VT
1737                       (CWorkspace: need   N*N+3*N-1,
1738                                    prefer N*N+2*N+(N-1)*NB)
1739                       (RWorkspace: 0) */
1740
1741                         i__2 = *lwork - iwork + 1;
1742                         zungbr_("P", n, n, n, &vt[vt_offset], ldvt, &work[
1743                                 itaup], &work[iwork], &i__2, &ierr)
1744                                 ;
1745                         irwork = ie + *n;
1746
1747 /*                    Perform bidiagonal QR iteration, computing left
1748                       singular vectors of R in WORK(IU) and computing
1749                       right singular vectors of R in VT
1750                       (CWorkspace: need N*N)
1751                       (RWorkspace: need BDSPAC) */
1752
1753                         zbdsqr_("U", n, n, n, &c__0, &s[1], &rwork[ie], &vt[
1754                                 vt_offset], ldvt, &work[iu], &ldwrku, cdum, &
1755                                 c__1, &rwork[irwork], info);
1756
1757 /*                    Multiply Q in A by left singular vectors of R in
1758                       WORK(IU), storing result in U
1759                       (CWorkspace: need N*N)
1760                       (RWorkspace: 0) */
1761
1762                         zgemm_("N", "N", m, n, n, &c_b2, &a[a_offset], lda, &
1763                                 work[iu], &ldwrku, &c_b1, &u[u_offset], ldu);
1764
1765                     } else {
1766
1767 /*                    Insufficient workspace for a fast algorithm */
1768
1769                         itau = 1;
1770                         iwork = itau + *n;
1771
1772 /*                    Compute A=Q*R, copying result to U
1773                       (CWorkspace: need 2*N, prefer N+N*NB)
1774                       (RWorkspace: 0) */
1775
1776                         i__2 = *lwork - iwork + 1;
1777                         zgeqrf_(m, n, &a[a_offset], lda, &work[itau], &work[
1778                                 iwork], &i__2, &ierr);
1779                         zlacpy_("L", m, n, &a[a_offset], lda, &u[u_offset],
1780                                 ldu);
1781
1782 /*                    Generate Q in U
1783                       (CWorkspace: need 2*N, prefer N+N*NB)
1784                       (RWorkspace: 0) */
1785
1786                         i__2 = *lwork - iwork + 1;
1787                         zungqr_(m, n, n, &u[u_offset], ldu, &work[itau], &
1788                                 work[iwork], &i__2, &ierr);
1789
1790 /*                    Copy R to VT, zeroing out below it */
1791
1792                         zlacpy_("U", n, n, &a[a_offset], lda, &vt[vt_offset],
1793                                 ldvt);
1794                         i__2 = *n - 1;
1795                         i__3 = *n - 1;
1796                         zlaset_("L", &i__2, &i__3, &c_b1, &c_b1, &vt_ref(2, 1)
1797                                 , ldvt);
1798                         ie = 1;
1799                         itauq = itau;
1800                         itaup = itauq + *n;
1801                         iwork = itaup + *n;
1802
1803 /*                    Bidiagonalize R in VT
1804                       (CWorkspace: need 3*N, prefer 2*N+2*N*NB)
1805                       (RWorkspace: need N) */
1806
1807                         i__2 = *lwork - iwork + 1;
1808                         zgebrd_(n, n, &vt[vt_offset], ldvt, &s[1], &rwork[ie],
1809                                  &work[itauq], &work[itaup], &work[iwork], &
1810                                 i__2, &ierr);
1811
1812 /*                    Multiply Q in U by left bidiagonalizing vectors
1813                       in VT
1814                       (CWorkspace: need 2*N+M, prefer 2*N+M*NB)
1815                       (RWorkspace: 0) */
1816
1817                         i__2 = *lwork - iwork + 1;
1818                         zunmbr_("Q", "R", "N", m, n, n, &vt[vt_offset], ldvt,
1819                                 &work[itauq], &u[u_offset], ldu, &work[iwork],
1820                                  &i__2, &ierr);
1821
1822 /*                    Generate right bidiagonalizing vectors in VT
1823                       (CWorkspace: need 3*N-1, prefer 2*N+(N-1)*NB)
1824                       (RWorkspace: 0) */
1825
1826                         i__2 = *lwork - iwork + 1;
1827                         zungbr_("P", n, n, n, &vt[vt_offset], ldvt, &work[
1828                                 itaup], &work[iwork], &i__2, &ierr)
1829                                 ;
1830                         irwork = ie + *n;
1831
1832 /*                    Perform bidiagonal QR iteration, computing left
1833                       singular vectors of A in U and computing right
1834                       singular vectors of A in VT
1835                       (CWorkspace: 0)
1836                       (RWorkspace: need BDSPAC) */
1837
1838                         zbdsqr_("U", n, n, m, &c__0, &s[1], &rwork[ie], &vt[
1839                                 vt_offset], ldvt, &u[u_offset], ldu, cdum, &
1840                                 c__1, &rwork[irwork], info);
1841
1842                     }
1843
1844                 }
1845
1846             } else if (wntua) {
1847
1848                 if (wntvn) {
1849
1850 /*                 Path 7 (M much larger than N, JOBU='A', JOBVT='N')
1851                    M left singular vectors to be computed in U and
1852                    no right singular vectors to be computed
1853
1854    Computing MAX */
1855                     i__2 = *n + *m, i__3 = *n * 3;
1856                     if (*lwork >= *n * *n + max(i__2,i__3)) {
1857
1858 /*                    Sufficient workspace for a fast algorithm */
1859
1860                         ir = 1;
1861                         if (*lwork >= wrkbl + *lda * *n) {
1862
1863 /*                       WORK(IR) is LDA by N */
1864
1865                             ldwrkr = *lda;
1866                         } else {
1867
1868 /*                       WORK(IR) is N by N */
1869
1870                             ldwrkr = *n;
1871                         }
1872                         itau = ir + ldwrkr * *n;
1873                         iwork = itau + *n;
1874
1875 /*                    Compute A=Q*R, copying result to U
1876                       (CWorkspace: need N*N+2*N, prefer N*N+N+N*NB)
1877                       (RWorkspace: 0) */
1878
1879                         i__2 = *lwork - iwork + 1;
1880                         zgeqrf_(m, n, &a[a_offset], lda, &work[itau], &work[
1881                                 iwork], &i__2, &ierr);
1882                         zlacpy_("L", m, n, &a[a_offset], lda, &u[u_offset],
1883                                 ldu);
1884
1885 /*                    Copy R to WORK(IR), zeroing out below it */
1886
1887                         zlacpy_("U", n, n, &a[a_offset], lda, &work[ir], &
1888                                 ldwrkr);
1889                         i__2 = *n - 1;
1890                         i__3 = *n - 1;
1891                         zlaset_("L", &i__2, &i__3, &c_b1, &c_b1, &work[ir + 1]
1892                                 , &ldwrkr);
1893
1894 /*                    Generate Q in U
1895                       (CWorkspace: need N*N+N+M, prefer N*N+N+M*NB)
1896                       (RWorkspace: 0) */
1897
1898                         i__2 = *lwork - iwork + 1;
1899                         zungqr_(m, m, n, &u[u_offset], ldu, &work[itau], &
1900                                 work[iwork], &i__2, &ierr);
1901                         ie = 1;
1902                         itauq = itau;
1903                         itaup = itauq + *n;
1904                         iwork = itaup + *n;
1905
1906 /*                    Bidiagonalize R in WORK(IR)
1907                       (CWorkspace: need N*N+3*N, prefer N*N+2*N+2*N*NB)
1908                       (RWorkspace: need N) */
1909
1910                         i__2 = *lwork - iwork + 1;
1911                         zgebrd_(n, n, &work[ir], &ldwrkr, &s[1], &rwork[ie], &
1912                                 work[itauq], &work[itaup], &work[iwork], &
1913                                 i__2, &ierr);
1914
1915 /*                    Generate left bidiagonalizing vectors in WORK(IR)
1916                       (CWorkspace: need N*N+3*N, prefer N*N+2*N+N*NB)
1917                       (RWorkspace: 0) */
1918
1919                         i__2 = *lwork - iwork + 1;
1920                         zungbr_("Q", n, n, n, &work[ir], &ldwrkr, &work[itauq]
1921                                 , &work[iwork], &i__2, &ierr);
1922                         irwork = ie + *n;
1923
1924 /*                    Perform bidiagonal QR iteration, computing left
1925                       singular vectors of R in WORK(IR)
1926                       (CWorkspace: need N*N)
1927                       (RWorkspace: need BDSPAC) */
1928
1929                         zbdsqr_("U", n, &c__0, n, &c__0, &s[1], &rwork[ie],
1930                                 cdum, &c__1, &work[ir], &ldwrkr, cdum, &c__1,
1931                                 &rwork[irwork], info);
1932
1933 /*                    Multiply Q in U by left singular vectors of R in
1934                       WORK(IR), storing result in A
1935                       (CWorkspace: need N*N)
1936                       (RWorkspace: 0) */
1937
1938                         zgemm_("N", "N", m, n, n, &c_b2, &u[u_offset], ldu, &
1939                                 work[ir], &ldwrkr, &c_b1, &a[a_offset], lda);
1940
1941 /*                    Copy left singular vectors of A from A to U */
1942
1943                         zlacpy_("F", m, n, &a[a_offset], lda, &u[u_offset],
1944                                 ldu);
1945
1946                     } else {
1947
1948 /*                    Insufficient workspace for a fast algorithm */
1949
1950                         itau = 1;
1951                         iwork = itau + *n;
1952
1953 /*                    Compute A=Q*R, copying result to U
1954                       (CWorkspace: need 2*N, prefer N+N*NB)
1955                       (RWorkspace: 0) */
1956
1957                         i__2 = *lwork - iwork + 1;
1958                         zgeqrf_(m, n, &a[a_offset], lda, &work[itau], &work[
1959                                 iwork], &i__2, &ierr);
1960                         zlacpy_("L", m, n, &a[a_offset], lda, &u[u_offset],
1961                                 ldu);
1962
1963 /*                    Generate Q in U
1964                       (CWorkspace: need N+M, prefer N+M*NB)
1965                       (RWorkspace: 0) */
1966
1967                         i__2 = *lwork - iwork + 1;
1968                         zungqr_(m, m, n, &u[u_offset], ldu, &work[itau], &
1969                                 work[iwork], &i__2, &ierr);
1970                         ie = 1;
1971                         itauq = itau;
1972                         itaup = itauq + *n;
1973                         iwork = itaup + *n;
1974
1975 /*                    Zero out below R in A */
1976
1977                         i__2 = *n - 1;
1978                         i__3 = *n - 1;
1979                         zlaset_("L", &i__2, &i__3, &c_b1, &c_b1, &a_ref(2, 1),
1980                                  lda);
1981
1982 /*                    Bidiagonalize R in A
1983                       (CWorkspace: need 3*N, prefer 2*N+2*N*NB)
1984                       (RWorkspace: need N) */
1985
1986                         i__2 = *lwork - iwork + 1;
1987                         zgebrd_(n, n, &a[a_offset], lda, &s[1], &rwork[ie], &
1988                                 work[itauq], &work[itaup], &work[iwork], &
1989                                 i__2, &ierr);
1990
1991 /*                    Multiply Q in U by left bidiagonalizing vectors
1992                       in A
1993                       (CWorkspace: need 2*N+M, prefer 2*N+M*NB)
1994                       (RWorkspace: 0) */
1995
1996                         i__2 = *lwork - iwork + 1;
1997                         zunmbr_("Q", "R", "N", m, n, n, &a[a_offset], lda, &
1998                                 work[itauq], &u[u_offset], ldu, &work[iwork],
1999                                 &i__2, &ierr)
2000                                 ;
2001                         irwork = ie + *n;
2002
2003 /*                    Perform bidiagonal QR iteration, computing left
2004                       singular vectors of A in U
2005                       (CWorkspace: 0)
2006                       (RWorkspace: need BDSPAC) */
2007
2008                         zbdsqr_("U", n, &c__0, m, &c__0, &s[1], &rwork[ie],
2009                                 cdum, &c__1, &u[u_offset], ldu, cdum, &c__1, &
2010                                 rwork[irwork], info);
2011
2012                     }
2013
2014                 } else if (wntvo) {
2015
2016 /*                 Path 8 (M much larger than N, JOBU='A', JOBVT='O')
2017                    M left singular vectors to be computed in U and
2018                    N right singular vectors to be overwritten on A
2019
2020    Computing MAX */
2021                     i__2 = *n + *m, i__3 = *n * 3;
2022                     if (*lwork >= (*n << 1) * *n + max(i__2,i__3)) {
2023
2024 /*                    Sufficient workspace for a fast algorithm */
2025
2026                         iu = 1;
2027                         if (*lwork >= wrkbl + (*lda << 1) * *n) {
2028
2029 /*                       WORK(IU) is LDA by N and WORK(IR) is LDA by N */
2030
2031                             ldwrku = *lda;
2032                             ir = iu + ldwrku * *n;
2033                             ldwrkr = *lda;
2034                         } else if (*lwork >= wrkbl + (*lda + *n) * *n) {
2035
2036 /*                       WORK(IU) is LDA by N and WORK(IR) is N by N */
2037
2038                             ldwrku = *lda;
2039                             ir = iu + ldwrku * *n;
2040                             ldwrkr = *n;
2041                         } else {
2042
2043 /*                       WORK(IU) is N by N and WORK(IR) is N by N */
2044
2045                             ldwrku = *n;
2046                             ir = iu + ldwrku * *n;
2047                             ldwrkr = *n;
2048                         }
2049                         itau = ir + ldwrkr * *n;
2050                         iwork = itau + *n;
2051
2052 /*                    Compute A=Q*R, copying result to U
2053                       (CWorkspace: need 2*N*N+2*N, prefer 2*N*N+N+N*NB)
2054                       (RWorkspace: 0) */
2055
2056                         i__2 = *lwork - iwork + 1;
2057                         zgeqrf_(m, n, &a[a_offset], lda, &work[itau], &work[
2058                                 iwork], &i__2, &ierr);
2059                         zlacpy_("L", m, n, &a[a_offset], lda, &u[u_offset],
2060                                 ldu);
2061
2062 /*                    Generate Q in U
2063                       (CWorkspace: need 2*N*N+N+M, prefer 2*N*N+N+M*NB)
2064                       (RWorkspace: 0) */
2065
2066                         i__2 = *lwork - iwork + 1;
2067                         zungqr_(m, m, n, &u[u_offset], ldu, &work[itau], &
2068                                 work[iwork], &i__2, &ierr);
2069
2070 /*                    Copy R to WORK(IU), zeroing out below it */
2071
2072                         zlacpy_("U", n, n, &a[a_offset], lda, &work[iu], &
2073                                 ldwrku);
2074                         i__2 = *n - 1;
2075                         i__3 = *n - 1;
2076                         zlaset_("L", &i__2, &i__3, &c_b1, &c_b1, &work[iu + 1]
2077                                 , &ldwrku);
2078                         ie = 1;
2079                         itauq = itau;
2080                         itaup = itauq + *n;
2081                         iwork = itaup + *n;
2082
2083 /*                    Bidiagonalize R in WORK(IU), copying result to
2084                       WORK(IR)
2085                       (CWorkspace: need   2*N*N+3*N,
2086                                    prefer 2*N*N+2*N+2*N*NB)
2087                       (RWorkspace: need   N) */
2088
2089                         i__2 = *lwork - iwork + 1;
2090                         zgebrd_(n, n, &work[iu], &ldwrku, &s[1], &rwork[ie], &
2091                                 work[itauq], &work[itaup], &work[iwork], &
2092                                 i__2, &ierr);
2093                         zlacpy_("U", n, n, &work[iu], &ldwrku, &work[ir], &
2094                                 ldwrkr);
2095
2096 /*                    Generate left bidiagonalizing vectors in WORK(IU)
2097                       (CWorkspace: need 2*N*N+3*N, prefer 2*N*N+2*N+N*NB)
2098                       (RWorkspace: 0) */
2099
2100                         i__2 = *lwork - iwork + 1;
2101                         zungbr_("Q", n, n, n, &work[iu], &ldwrku, &work[itauq]
2102                                 , &work[iwork], &i__2, &ierr);
2103
2104 /*                    Generate right bidiagonalizing vectors in WORK(IR)
2105                       (CWorkspace: need   2*N*N+3*N-1,
2106                                    prefer 2*N*N+2*N+(N-1)*NB)
2107                       (RWorkspace: 0) */
2108
2109                         i__2 = *lwork - iwork + 1;
2110                         zungbr_("P", n, n, n, &work[ir], &ldwrkr, &work[itaup]
2111                                 , &work[iwork], &i__2, &ierr);
2112                         irwork = ie + *n;
2113
2114 /*                    Perform bidiagonal QR iteration, computing left
2115                       singular vectors of R in WORK(IU) and computing
2116                       right singular vectors of R in WORK(IR)
2117                       (CWorkspace: need 2*N*N)
2118                       (RWorkspace: need BDSPAC) */
2119
2120                         zbdsqr_("U", n, n, n, &c__0, &s[1], &rwork[ie], &work[
2121                                 ir], &ldwrkr, &work[iu], &ldwrku, cdum, &c__1,
2122                                  &rwork[irwork], info);
2123
2124 /*                    Multiply Q in U by left singular vectors of R in
2125                       WORK(IU), storing result in A
2126                       (CWorkspace: need N*N)
2127                       (RWorkspace: 0) */
2128
2129                         zgemm_("N", "N", m, n, n, &c_b2, &u[u_offset], ldu, &
2130                                 work[iu], &ldwrku, &c_b1, &a[a_offset], lda);
2131
2132 /*                    Copy left singular vectors of A from A to U */
2133
2134                         zlacpy_("F", m, n, &a[a_offset], lda, &u[u_offset],
2135                                 ldu);
2136
2137 /*                    Copy right singular vectors of R from WORK(IR) to A */
2138
2139                         zlacpy_("F", n, n, &work[ir], &ldwrkr, &a[a_offset],
2140                                 lda);
2141
2142                     } else {
2143
2144 /*                    Insufficient workspace for a fast algorithm */
2145
2146                         itau = 1;
2147                         iwork = itau + *n;
2148
2149 /*                    Compute A=Q*R, copying result to U
2150                       (CWorkspace: need 2*N, prefer N+N*NB)
2151                       (RWorkspace: 0) */
2152
2153                         i__2 = *lwork - iwork + 1;
2154                         zgeqrf_(m, n, &a[a_offset], lda, &work[itau], &work[
2155                                 iwork], &i__2, &ierr);
2156                         zlacpy_("L", m, n, &a[a_offset], lda, &u[u_offset],
2157                                 ldu);
2158
2159 /*                    Generate Q in U
2160                       (CWorkspace: need N+M, prefer N+M*NB)
2161                       (RWorkspace: 0) */
2162
2163                         i__2 = *lwork - iwork + 1;
2164                         zungqr_(m, m, n, &u[u_offset], ldu, &work[itau], &
2165                                 work[iwork], &i__2, &ierr);
2166                         ie = 1;
2167                         itauq = itau;
2168                         itaup = itauq + *n;
2169                         iwork = itaup + *n;
2170
2171 /*                    Zero out below R in A */
2172
2173                         i__2 = *n - 1;
2174                         i__3 = *n - 1;
2175                         zlaset_("L", &i__2, &i__3, &c_b1, &c_b1, &a_ref(2, 1),
2176                                  lda);
2177
2178 /*                    Bidiagonalize R in A
2179                       (CWorkspace: need 3*N, prefer 2*N+2*N*NB)
2180                       (RWorkspace: need N) */
2181
2182                         i__2 = *lwork - iwork + 1;
2183                         zgebrd_(n, n, &a[a_offset], lda, &s[1], &rwork[ie], &
2184                                 work[itauq], &work[itaup], &work[iwork], &
2185                                 i__2, &ierr);
2186
2187 /*                    Multiply Q in U by left bidiagonalizing vectors
2188                       in A
2189                       (CWorkspace: need 2*N+M, prefer 2*N+M*NB)
2190                       (RWorkspace: 0) */
2191
2192                         i__2 = *lwork - iwork + 1;
2193                         zunmbr_("Q", "R", "N", m, n, n, &a[a_offset], lda, &
2194                                 work[itauq], &u[u_offset], ldu, &work[iwork],
2195                                 &i__2, &ierr)
2196                                 ;
2197
2198 /*                    Generate right bidiagonalizing vectors in A
2199                       (CWorkspace: need 3*N-1, prefer 2*N+(N-1)*NB)
2200                       (RWorkspace: 0) */
2201
2202                         i__2 = *lwork - iwork + 1;
2203                         zungbr_("P", n, n, n, &a[a_offset], lda, &work[itaup],
2204                                  &work[iwork], &i__2, &ierr);
2205                         irwork = ie + *n;
2206
2207 /*                    Perform bidiagonal QR iteration, computing left
2208                       singular vectors of A in U and computing right
2209                       singular vectors of A in A
2210                       (CWorkspace: 0)
2211                       (RWorkspace: need BDSPAC) */
2212
2213                         zbdsqr_("U", n, n, m, &c__0, &s[1], &rwork[ie], &a[
2214                                 a_offset], lda, &u[u_offset], ldu, cdum, &
2215                                 c__1, &rwork[irwork], info);
2216
2217                     }
2218
2219                 } else if (wntvas) {
2220
2221 /*                 Path 9 (M much larger than N, JOBU='A', JOBVT='S'
2222                            or 'A')
2223                    M left singular vectors to be computed in U and
2224                    N right singular vectors to be computed in VT
2225
2226    Computing MAX */
2227                     i__2 = *n + *m, i__3 = *n * 3;
2228                     if (*lwork >= *n * *n + max(i__2,i__3)) {
2229
2230 /*                    Sufficient workspace for a fast algorithm */
2231
2232                         iu = 1;
2233                         if (*lwork >= wrkbl + *lda * *n) {
2234
2235 /*                       WORK(IU) is LDA by N */
2236
2237                             ldwrku = *lda;
2238                         } else {
2239
2240 /*                       WORK(IU) is N by N */
2241
2242                             ldwrku = *n;
2243                         }
2244                         itau = iu + ldwrku * *n;
2245                         iwork = itau + *n;
2246
2247 /*                    Compute A=Q*R, copying result to U
2248                       (CWorkspace: need N*N+2*N, prefer N*N+N+N*NB)
2249                       (RWorkspace: 0) */
2250
2251                         i__2 = *lwork - iwork + 1;
2252                         zgeqrf_(m, n, &a[a_offset], lda, &work[itau], &work[
2253                                 iwork], &i__2, &ierr);
2254                         zlacpy_("L", m, n, &a[a_offset], lda, &u[u_offset],
2255                                 ldu);
2256
2257 /*                    Generate Q in U
2258                       (CWorkspace: need N*N+N+M, prefer N*N+N+M*NB)
2259                       (RWorkspace: 0) */
2260
2261                         i__2 = *lwork - iwork + 1;
2262                         zungqr_(m, m, n, &u[u_offset], ldu, &work[itau], &
2263                                 work[iwork], &i__2, &ierr);
2264
2265 /*                    Copy R to WORK(IU), zeroing out below it */
2266
2267                         zlacpy_("U", n, n, &a[a_offset], lda, &work[iu], &
2268                                 ldwrku);
2269                         i__2 = *n - 1;
2270                         i__3 = *n - 1;
2271                         zlaset_("L", &i__2, &i__3, &c_b1, &c_b1, &work[iu + 1]
2272                                 , &ldwrku);
2273                         ie = 1;
2274                         itauq = itau;
2275                         itaup = itauq + *n;
2276                         iwork = itaup + *n;
2277
2278 /*                    Bidiagonalize R in WORK(IU), copying result to VT
2279                       (CWorkspace: need N*N+3*N, prefer N*N+2*N+2*N*NB)
2280                       (RWorkspace: need N) */
2281
2282                         i__2 = *lwork - iwork + 1;
2283                         zgebrd_(n, n, &work[iu], &ldwrku, &s[1], &rwork[ie], &
2284                                 work[itauq], &work[itaup], &work[iwork], &
2285                                 i__2, &ierr);
2286                         zlacpy_("U", n, n, &work[iu], &ldwrku, &vt[vt_offset],
2287                                  ldvt);
2288
2289 /*                    Generate left bidiagonalizing vectors in WORK(IU)
2290                       (CWorkspace: need N*N+3*N, prefer N*N+2*N+N*NB)
2291                       (RWorkspace: 0) */
2292
2293                         i__2 = *lwork - iwork + 1;
2294                         zungbr_("Q", n, n, n, &work[iu], &ldwrku, &work[itauq]
2295                                 , &work[iwork], &i__2, &ierr);
2296
2297 /*                    Generate right bidiagonalizing vectors in VT
2298                       (CWorkspace: need   N*N+3*N-1,
2299                                    prefer N*N+2*N+(N-1)*NB)
2300                       (RWorkspace: need   0) */
2301
2302                         i__2 = *lwork - iwork + 1;
2303                         zungbr_("P", n, n, n, &vt[vt_offset], ldvt, &work[
2304                                 itaup], &work[iwork], &i__2, &ierr)
2305                                 ;
2306                         irwork = ie + *n;
2307
2308 /*                    Perform bidiagonal QR iteration, computing left
2309                       singular vectors of R in WORK(IU) and computing
2310                       right singular vectors of R in VT
2311                       (CWorkspace: need N*N)
2312                       (RWorkspace: need BDSPAC) */
2313
2314                         zbdsqr_("U", n, n, n, &c__0, &s[1], &rwork[ie], &vt[
2315                                 vt_offset], ldvt, &work[iu], &ldwrku, cdum, &
2316                                 c__1, &rwork[irwork], info);
2317
2318 /*                    Multiply Q in U by left singular vectors of R in
2319                       WORK(IU), storing result in A
2320                       (CWorkspace: need N*N)
2321                       (RWorkspace: 0) */
2322
2323                         zgemm_("N", "N", m, n, n, &c_b2, &u[u_offset], ldu, &
2324                                 work[iu], &ldwrku, &c_b1, &a[a_offset], lda);
2325
2326 /*                    Copy left singular vectors of A from A to U */
2327
2328                         zlacpy_("F", m, n, &a[a_offset], lda, &u[u_offset],
2329                                 ldu);
2330
2331                     } else {
2332
2333 /*                    Insufficient workspace for a fast algorithm */
2334
2335                         itau = 1;
2336                         iwork = itau + *n;
2337
2338 /*                    Compute A=Q*R, copying result to U
2339                       (CWorkspace: need 2*N, prefer N+N*NB)
2340                       (RWorkspace: 0) */
2341
2342                         i__2 = *lwork - iwork + 1;
2343                         zgeqrf_(m, n, &a[a_offset], lda, &work[itau], &work[
2344                                 iwork], &i__2, &ierr);
2345                         zlacpy_("L", m, n, &a[a_offset], lda, &u[u_offset],
2346                                 ldu);
2347
2348 /*                    Generate Q in U
2349                       (CWorkspace: need N+M, prefer N+M*NB)
2350                       (RWorkspace: 0) */
2351
2352                         i__2 = *lwork - iwork + 1;
2353                         zungqr_(m, m, n, &u[u_offset], ldu, &work[itau], &
2354                                 work[iwork], &i__2, &ierr);
2355
2356 /*                    Copy R from A to VT, zeroing out below it */
2357
2358                         zlacpy_("U", n, n, &a[a_offset], lda, &vt[vt_offset],
2359                                 ldvt);
2360                         i__2 = *n - 1;
2361                         i__3 = *n - 1;
2362                         zlaset_("L", &i__2, &i__3, &c_b1, &c_b1, &vt_ref(2, 1)
2363                                 , ldvt);
2364                         ie = 1;
2365                         itauq = itau;
2366                         itaup = itauq + *n;
2367                         iwork = itaup + *n;
2368
2369 /*                    Bidiagonalize R in VT
2370                       (CWorkspace: need 3*N, prefer 2*N+2*N*NB)
2371                       (RWorkspace: need N) */
2372
2373                         i__2 = *lwork - iwork + 1;
2374                         zgebrd_(n, n, &vt[vt_offset], ldvt, &s[1], &rwork[ie],
2375                                  &work[itauq], &work[itaup], &work[iwork], &
2376                                 i__2, &ierr);
2377
2378 /*                    Multiply Q in U by left bidiagonalizing vectors
2379                       in VT
2380                       (CWorkspace: need 2*N+M, prefer 2*N+M*NB)
2381                       (RWorkspace: 0) */
2382
2383                         i__2 = *lwork - iwork + 1;
2384                         zunmbr_("Q", "R", "N", m, n, n, &vt[vt_offset], ldvt,
2385                                 &work[itauq], &u[u_offset], ldu, &work[iwork],
2386                                  &i__2, &ierr);
2387
2388 /*                    Generate right bidiagonalizing vectors in VT
2389                       (CWorkspace: need 3*N-1, prefer 2*N+(N-1)*NB)
2390                       (RWorkspace: 0) */
2391
2392                         i__2 = *lwork - iwork + 1;
2393                         zungbr_("P", n, n, n, &vt[vt_offset], ldvt, &work[
2394                                 itaup], &work[iwork], &i__2, &ierr)
2395                                 ;
2396                         irwork = ie + *n;
2397
2398 /*                    Perform bidiagonal QR iteration, computing left
2399                       singular vectors of A in U and computing right
2400                       singular vectors of A in VT
2401                       (CWorkspace: 0)
2402                       (RWorkspace: need BDSPAC) */
2403
2404                         zbdsqr_("U", n, n, m, &c__0, &s[1], &rwork[ie], &vt[
2405                                 vt_offset], ldvt, &u[u_offset], ldu, cdum, &
2406                                 c__1, &rwork[irwork], info);
2407
2408                     }
2409
2410                 }
2411
2412             }
2413
2414         } else {
2415
2416 /*           M .LT. MNTHR
2417
2418              Path 10 (M at least N, but not much larger)
2419              Reduce to bidiagonal form without QR decomposition */
2420
2421             ie = 1;
2422             itauq = 1;
2423             itaup = itauq + *n;
2424             iwork = itaup + *n;
2425
2426 /*           Bidiagonalize A
2427              (CWorkspace: need 2*N+M, prefer 2*N+(M+N)*NB)
2428              (RWorkspace: need N) */
2429
2430             i__2 = *lwork - iwork + 1;
2431             zgebrd_(m, n, &a[a_offset], lda, &s[1], &rwork[ie], &work[itauq],
2432                     &work[itaup], &work[iwork], &i__2, &ierr);
2433             if (wntuas) {
2434
2435 /*              If left singular vectors desired in U, copy result to U
2436                 and generate left bidiagonalizing vectors in U
2437                 (CWorkspace: need 2*N+NCU, prefer 2*N+NCU*NB)
2438                 (RWorkspace: 0) */
2439
2440                 zlacpy_("L", m, n, &a[a_offset], lda, &u[u_offset], ldu);
2441                 if (wntus) {
2442                     ncu = *n;
2443                 }
2444                 if (wntua) {
2445                     ncu = *m;
2446                 }
2447                 i__2 = *lwork - iwork + 1;
2448                 zungbr_("Q", m, &ncu, n, &u[u_offset], ldu, &work[itauq], &
2449                         work[iwork], &i__2, &ierr);
2450             }
2451             if (wntvas) {
2452
2453 /*              If right singular vectors desired in VT, copy result to
2454                 VT and generate right bidiagonalizing vectors in VT
2455                 (CWorkspace: need 3*N-1, prefer 2*N+(N-1)*NB)
2456                 (RWorkspace: 0) */
2457
2458                 zlacpy_("U", n, n, &a[a_offset], lda, &vt[vt_offset], ldvt);
2459                 i__2 = *lwork - iwork + 1;
2460                 zungbr_("P", n, n, n, &vt[vt_offset], ldvt, &work[itaup], &
2461                         work[iwork], &i__2, &ierr);
2462             }
2463             if (wntuo) {
2464
2465 /*              If left singular vectors desired in A, generate left
2466                 bidiagonalizing vectors in A
2467                 (CWorkspace: need 3*N, prefer 2*N+N*NB)
2468                 (RWorkspace: 0) */
2469
2470                 i__2 = *lwork - iwork + 1;
2471                 zungbr_("Q", m, n, n, &a[a_offset], lda, &work[itauq], &work[
2472                         iwork], &i__2, &ierr);
2473             }
2474             if (wntvo) {
2475
2476 /*              If right singular vectors desired in A, generate right
2477                 bidiagonalizing vectors in A
2478                 (CWorkspace: need 3*N-1, prefer 2*N+(N-1)*NB)
2479                 (RWorkspace: 0) */
2480
2481                 i__2 = *lwork - iwork + 1;
2482                 zungbr_("P", n, n, n, &a[a_offset], lda, &work[itaup], &work[
2483                         iwork], &i__2, &ierr);
2484             }
2485             irwork = ie + *n;
2486             if (wntuas || wntuo) {
2487                 nru = *m;
2488             }
2489             if (wntun) {
2490                 nru = 0;
2491             }
2492             if (wntvas || wntvo) {
2493                 ncvt = *n;
2494             }
2495             if (wntvn) {
2496                 ncvt = 0;
2497             }
2498             if (! wntuo && ! wntvo) {
2499
2500 /*              Perform bidiagonal QR iteration, if desired, computing
2501                 left singular vectors in U and computing right singular
2502                 vectors in VT
2503                 (CWorkspace: 0)
2504                 (RWorkspace: need BDSPAC) */
2505
2506                 zbdsqr_("U", n, &ncvt, &nru, &c__0, &s[1], &rwork[ie], &vt[
2507                         vt_offset], ldvt, &u[u_offset], ldu, cdum, &c__1, &
2508                         rwork[irwork], info);
2509             } else if (! wntuo && wntvo) {
2510
2511 /*              Perform bidiagonal QR iteration, if desired, computing
2512                 left singular vectors in U and computing right singular
2513                 vectors in A
2514                 (CWorkspace: 0)
2515                 (RWorkspace: need BDSPAC) */
2516
2517                 zbdsqr_("U", n, &ncvt, &nru, &c__0, &s[1], &rwork[ie], &a[
2518                         a_offset], lda, &u[u_offset], ldu, cdum, &c__1, &
2519                         rwork[irwork], info);
2520             } else {
2521
2522 /*              Perform bidiagonal QR iteration, if desired, computing
2523                 left singular vectors in A and computing right singular
2524                 vectors in VT
2525                 (CWorkspace: 0)
2526                 (RWorkspace: need BDSPAC) */
2527
2528                 zbdsqr_("U", n, &ncvt, &nru, &c__0, &s[1], &rwork[ie], &vt[
2529                         vt_offset], ldvt, &a[a_offset], lda, cdum, &c__1, &
2530                         rwork[irwork], info);
2531             }
2532
2533         }
2534
2535     } else {
2536
2537 /*        A has more columns than rows. If A has sufficiently more
2538           columns than rows, first reduce using the LQ decomposition (if
2539           sufficient workspace available) */
2540
2541         if (*n >= mnthr) {
2542
2543             if (wntvn) {
2544
2545 /*              Path 1t(N much larger than M, JOBVT='N')
2546                 No right singular vectors to be computed */
2547
2548                 itau = 1;
2549                 iwork = itau + *m;
2550
2551 /*              Compute A=L*Q
2552                 (CWorkspace: need 2*M, prefer M+M*NB)
2553                 (RWorkspace: 0) */
2554
2555                 i__2 = *lwork - iwork + 1;
2556                 zgelqf_(m, n, &a[a_offset], lda, &work[itau], &work[iwork], &
2557                         i__2, &ierr);
2558
2559 /*              Zero out above L */
2560
2561                 i__2 = *m - 1;
2562                 i__3 = *m - 1;
2563                 zlaset_("U", &i__2, &i__3, &c_b1, &c_b1, &a_ref(1, 2), lda);
2564                 ie = 1;
2565                 itauq = 1;
2566                 itaup = itauq + *m;
2567                 iwork = itaup + *m;
2568
2569 /*              Bidiagonalize L in A
2570                 (CWorkspace: need 3*M, prefer 2*M+2*M*NB)
2571                 (RWorkspace: need M) */
2572
2573                 i__2 = *lwork - iwork + 1;
2574                 zgebrd_(m, m, &a[a_offset], lda, &s[1], &rwork[ie], &work[
2575                         itauq], &work[itaup], &work[iwork], &i__2, &ierr);
2576                 if (wntuo || wntuas) {
2577
2578 /*                 If left singular vectors desired, generate Q
2579                    (CWorkspace: need 3*M, prefer 2*M+M*NB)
2580                    (RWorkspace: 0) */
2581
2582                     i__2 = *lwork - iwork + 1;
2583                     zungbr_("Q", m, m, m, &a[a_offset], lda, &work[itauq], &
2584                             work[iwork], &i__2, &ierr);
2585                 }
2586                 irwork = ie + *m;
2587                 nru = 0;
2588                 if (wntuo || wntuas) {
2589                     nru = *m;
2590                 }
2591
2592 /*              Perform bidiagonal QR iteration, computing left singular
2593                 vectors of A in A if desired
2594                 (CWorkspace: 0)
2595                 (RWorkspace: need BDSPAC) */
2596
2597                 zbdsqr_("U", m, &c__0, &nru, &c__0, &s[1], &rwork[ie], cdum, &
2598                         c__1, &a[a_offset], lda, cdum, &c__1, &rwork[irwork],
2599                         info);
2600
2601 /*              If left singular vectors desired in U, copy them there */
2602
2603                 if (wntuas) {
2604                     zlacpy_("F", m, m, &a[a_offset], lda, &u[u_offset], ldu);
2605                 }
2606
2607             } else if (wntvo && wntun) {
2608
2609 /*              Path 2t(N much larger than M, JOBU='N', JOBVT='O')
2610                 M right singular vectors to be overwritten on A and
2611                 no left singular vectors to be computed */
2612
2613                 if (*lwork >= *m * *m + *m * 3) {
2614
2615 /*                 Sufficient workspace for a fast algorithm */
2616
2617                     ir = 1;
2618 /* Computing MAX */
2619                     i__2 = wrkbl, i__3 = *lda * *n;
2620                     if (*lwork >= max(i__2,i__3) + *lda * *m) {
2621
2622 /*                    WORK(IU) is LDA by N and WORK(IR) is LDA by M */
2623
2624                         ldwrku = *lda;
2625                         chunk = *n;
2626                         ldwrkr = *lda;
2627                     } else /* if(complicated condition) */ {
2628 /* Computing MAX */
2629                         i__2 = wrkbl, i__3 = *lda * *n;
2630                         if (*lwork >= max(i__2,i__3) + *m * *m) {
2631
2632 /*                    WORK(IU) is LDA by N and WORK(IR) is M by M */
2633
2634                             ldwrku = *lda;
2635                             chunk = *n;
2636                             ldwrkr = *m;
2637                         } else {
2638
2639 /*                    WORK(IU) is M by CHUNK and WORK(IR) is M by M */
2640
2641                             ldwrku = *m;
2642                             chunk = (*lwork - *m * *m) / *m;
2643                             ldwrkr = *m;
2644                         }
2645                     }
2646                     itau = ir + ldwrkr * *m;
2647                     iwork = itau + *m;
2648
2649 /*                 Compute A=L*Q
2650                    (CWorkspace: need M*M+2*M, prefer M*M+M+M*NB)
2651                    (RWorkspace: 0) */
2652
2653                     i__2 = *lwork - iwork + 1;
2654                     zgelqf_(m, n, &a[a_offset], lda, &work[itau], &work[iwork]
2655                             , &i__2, &ierr);
2656
2657 /*                 Copy L to WORK(IR) and zero out above it */
2658
2659                     zlacpy_("L", m, m, &a[a_offset], lda, &work[ir], &ldwrkr);
2660                     i__2 = *m - 1;
2661                     i__3 = *m - 1;
2662                     zlaset_("U", &i__2, &i__3, &c_b1, &c_b1, &work[ir +
2663                             ldwrkr], &ldwrkr);
2664
2665 /*                 Generate Q in A
2666                    (CWorkspace: need M*M+2*M, prefer M*M+M+M*NB)
2667                    (RWorkspace: 0) */
2668
2669                     i__2 = *lwork - iwork + 1;
2670                     zunglq_(m, n, m, &a[a_offset], lda, &work[itau], &work[
2671                             iwork], &i__2, &ierr);
2672                     ie = 1;
2673                     itauq = itau;
2674                     itaup = itauq + *m;
2675                     iwork = itaup + *m;
2676
2677 /*                 Bidiagonalize L in WORK(IR)
2678                    (CWorkspace: need M*M+3*M, prefer M*M+2*M+2*M*NB)
2679                    (RWorkspace: need M) */
2680
2681                     i__2 = *lwork - iwork + 1;
2682                     zgebrd_(m, m, &work[ir], &ldwrkr, &s[1], &rwork[ie], &
2683                             work[itauq], &work[itaup], &work[iwork], &i__2, &
2684                             ierr);
2685
2686 /*                 Generate right vectors bidiagonalizing L
2687                    (CWorkspace: need M*M+3*M-1, prefer M*M+2*M+(M-1)*NB)
2688                    (RWorkspace: 0) */
2689
2690                     i__2 = *lwork - iwork + 1;
2691                     zungbr_("P", m, m, m, &work[ir], &ldwrkr, &work[itaup], &
2692                             work[iwork], &i__2, &ierr);
2693                     irwork = ie + *m;
2694
2695 /*                 Perform bidiagonal QR iteration, computing right
2696                    singular vectors of L in WORK(IR)
2697                    (CWorkspace: need M*M)
2698                    (RWorkspace: need BDSPAC) */
2699
2700                     zbdsqr_("U", m, m, &c__0, &c__0, &s[1], &rwork[ie], &work[
2701                             ir], &ldwrkr, cdum, &c__1, cdum, &c__1, &rwork[
2702                             irwork], info);
2703                     iu = itauq;
2704
2705 /*                 Multiply right singular vectors of L in WORK(IR) by Q
2706                    in A, storing result in WORK(IU) and copying to A
2707                    (CWorkspace: need M*M+M, prefer M*M+M*N)
2708                    (RWorkspace: 0) */
2709
2710                     i__2 = *n;
2711                     i__3 = chunk;
2712                     for (i__ = 1; i__3 < 0 ? i__ >= i__2 : i__ <= i__2; i__ +=
2713                              i__3) {
2714 /* Computing MIN */
2715                         i__4 = *n - i__ + 1;
2716                         blk = min(i__4,chunk);
2717                         zgemm_("N", "N", m, &blk, m, &c_b2, &work[ir], &
2718                                 ldwrkr, &a_ref(1, i__), lda, &c_b1, &work[iu],
2719                                  &ldwrku);
2720                         zlacpy_("F", m, &blk, &work[iu], &ldwrku, &a_ref(1,
2721                                 i__), lda);
2722 /* L30: */
2723                     }
2724
2725                 } else {
2726
2727 /*                 Insufficient workspace for a fast algorithm */
2728
2729                     ie = 1;
2730                     itauq = 1;
2731                     itaup = itauq + *m;
2732                     iwork = itaup + *m;
2733
2734 /*                 Bidiagonalize A
2735                    (CWorkspace: need 2*M+N, prefer 2*M+(M+N)*NB)
2736                    (RWorkspace: need M) */
2737
2738                     i__3 = *lwork - iwork + 1;
2739                     zgebrd_(m, n, &a[a_offset], lda, &s[1], &rwork[ie], &work[
2740                             itauq], &work[itaup], &work[iwork], &i__3, &ierr);
2741
2742 /*                 Generate right vectors bidiagonalizing A
2743                    (CWorkspace: need 3*M, prefer 2*M+M*NB)
2744                    (RWorkspace: 0) */
2745
2746                     i__3 = *lwork - iwork + 1;
2747                     zungbr_("P", m, n, m, &a[a_offset], lda, &work[itaup], &
2748                             work[iwork], &i__3, &ierr);
2749                     irwork = ie + *m;
2750
2751 /*                 Perform bidiagonal QR iteration, computing right
2752                    singular vectors of A in A
2753                    (CWorkspace: 0)
2754                    (RWorkspace: need BDSPAC) */
2755
2756                     zbdsqr_("L", m, n, &c__0, &c__0, &s[1], &rwork[ie], &a[
2757                             a_offset], lda, cdum, &c__1, cdum, &c__1, &rwork[
2758                             irwork], info);
2759
2760                 }
2761
2762             } else if (wntvo && wntuas) {
2763
2764 /*              Path 3t(N much larger than M, JOBU='S' or 'A', JOBVT='O')
2765                 M right singular vectors to be overwritten on A and
2766                 M left singular vectors to be computed in U */
2767
2768                 if (*lwork >= *m * *m + *m * 3) {
2769
2770 /*                 Sufficient workspace for a fast algorithm */
2771
2772                     ir = 1;
2773 /* Computing MAX */
2774                     i__3 = wrkbl, i__2 = *lda * *n;
2775                     if (*lwork >= max(i__3,i__2) + *lda * *m) {
2776
2777 /*                    WORK(IU) is LDA by N and WORK(IR) is LDA by M */
2778
2779                         ldwrku = *lda;
2780                         chunk = *n;
2781                         ldwrkr = *lda;
2782                     } else /* if(complicated condition) */ {
2783 /* Computing MAX */
2784                         i__3 = wrkbl, i__2 = *lda * *n;
2785                         if (*lwork >= max(i__3,i__2) + *m * *m) {
2786
2787 /*                    WORK(IU) is LDA by N and WORK(IR) is M by M */
2788
2789                             ldwrku = *lda;
2790                             chunk = *n;
2791                             ldwrkr = *m;
2792                         } else {
2793
2794 /*                    WORK(IU) is M by CHUNK and WORK(IR) is M by M */
2795
2796                             ldwrku = *m;
2797                             chunk = (*lwork - *m * *m) / *m;
2798                             ldwrkr = *m;
2799                         }
2800                     }
2801                     itau = ir + ldwrkr * *m;
2802                     iwork = itau + *m;
2803
2804 /*                 Compute A=L*Q
2805                    (CWorkspace: need M*M+2*M, prefer M*M+M+M*NB)
2806                    (RWorkspace: 0) */
2807
2808                     i__3 = *lwork - iwork + 1;
2809                     zgelqf_(m, n, &a[a_offset], lda, &work[itau], &work[iwork]
2810                             , &i__3, &ierr);
2811
2812 /*                 Copy L to U, zeroing about above it */
2813
2814                     zlacpy_("L", m, m, &a[a_offset], lda, &u[u_offset], ldu);
2815                     i__3 = *m - 1;
2816                     i__2 = *m - 1;
2817                     zlaset_("U", &i__3, &i__2, &c_b1, &c_b1, &u_ref(1, 2),
2818                             ldu);
2819
2820 /*                 Generate Q in A
2821                    (CWorkspace: need M*M+2*M, prefer M*M+M+M*NB)
2822                    (RWorkspace: 0) */
2823
2824                     i__3 = *lwork - iwork + 1;
2825                     zunglq_(m, n, m, &a[a_offset], lda, &work[itau], &work[
2826                             iwork], &i__3, &ierr);
2827                     ie = 1;
2828                     itauq = itau;
2829                     itaup = itauq + *m;
2830                     iwork = itaup + *m;
2831
2832 /*                 Bidiagonalize L in U, copying result to WORK(IR)
2833                    (CWorkspace: need M*M+3*M, prefer M*M+2*M+2*M*NB)
2834                    (RWorkspace: need M) */
2835
2836                     i__3 = *lwork - iwork + 1;
2837                     zgebrd_(m, m, &u[u_offset], ldu, &s[1], &rwork[ie], &work[
2838                             itauq], &work[itaup], &work[iwork], &i__3, &ierr);
2839                     zlacpy_("U", m, m, &u[u_offset], ldu, &work[ir], &ldwrkr);
2840
2841 /*                 Generate right vectors bidiagonalizing L in WORK(IR)
2842                    (CWorkspace: need M*M+3*M-1, prefer M*M+2*M+(M-1)*NB)
2843                    (RWorkspace: 0) */
2844
2845                     i__3 = *lwork - iwork + 1;
2846                     zungbr_("P", m, m, m, &work[ir], &ldwrkr, &work[itaup], &
2847                             work[iwork], &i__3, &ierr);
2848
2849 /*                 Generate left vectors bidiagonalizing L in U
2850                    (CWorkspace: need M*M+3*M, prefer M*M+2*M+M*NB)
2851                    (RWorkspace: 0) */
2852
2853                     i__3 = *lwork - iwork + 1;
2854                     zungbr_("Q", m, m, m, &u[u_offset], ldu, &work[itauq], &
2855                             work[iwork], &i__3, &ierr);
2856                     irwork = ie + *m;
2857
2858 /*                 Perform bidiagonal QR iteration, computing left
2859                    singular vectors of L in U, and computing right
2860                    singular vectors of L in WORK(IR)
2861                    (CWorkspace: need M*M)
2862                    (RWorkspace: need BDSPAC) */
2863
2864                     zbdsqr_("U", m, m, m, &c__0, &s[1], &rwork[ie], &work[ir],
2865                              &ldwrkr, &u[u_offset], ldu, cdum, &c__1, &rwork[
2866                             irwork], info);
2867                     iu = itauq;
2868
2869 /*                 Multiply right singular vectors of L in WORK(IR) by Q
2870                    in A, storing result in WORK(IU) and copying to A
2871                    (CWorkspace: need M*M+M, prefer M*M+M*N))
2872                    (RWorkspace: 0) */
2873
2874                     i__3 = *n;
2875                     i__2 = chunk;
2876                     for (i__ = 1; i__2 < 0 ? i__ >= i__3 : i__ <= i__3; i__ +=
2877                              i__2) {
2878 /* Computing MIN */
2879                         i__4 = *n - i__ + 1;
2880                         blk = min(i__4,chunk);
2881                         zgemm_("N", "N", m, &blk, m, &c_b2, &work[ir], &
2882                                 ldwrkr, &a_ref(1, i__), lda, &c_b1, &work[iu],
2883                                  &ldwrku);
2884                         zlacpy_("F", m, &blk, &work[iu], &ldwrku, &a_ref(1,
2885                                 i__), lda);
2886 /* L40: */
2887                     }
2888
2889                 } else {
2890
2891 /*                 Insufficient workspace for a fast algorithm */
2892
2893                     itau = 1;
2894                     iwork = itau + *m;
2895
2896 /*                 Compute A=L*Q
2897                    (CWorkspace: need 2*M, prefer M+M*NB)
2898                    (RWorkspace: 0) */
2899
2900                     i__2 = *lwork - iwork + 1;
2901                     zgelqf_(m, n, &a[a_offset], lda, &work[itau], &work[iwork]
2902                             , &i__2, &ierr);
2903
2904 /*                 Copy L to U, zeroing out above it */
2905
2906                     zlacpy_("L", m, m, &a[a_offset], lda, &u[u_offset], ldu);
2907                     i__2 = *m - 1;
2908                     i__3 = *m - 1;
2909                     zlaset_("U", &i__2, &i__3, &c_b1, &c_b1, &u_ref(1, 2),
2910                             ldu);
2911
2912 /*                 Generate Q in A
2913                    (CWorkspace: need 2*M, prefer M+M*NB)
2914                    (RWorkspace: 0) */
2915
2916                     i__2 = *lwork - iwork + 1;
2917                     zunglq_(m, n, m, &a[a_offset], lda, &work[itau], &work[
2918                             iwork], &i__2, &ierr);
2919                     ie = 1;
2920                     itauq = itau;
2921                     itaup = itauq + *m;
2922                     iwork = itaup + *m;
2923
2924 /*                 Bidiagonalize L in U
2925                    (CWorkspace: need 3*M, prefer 2*M+2*M*NB)
2926                    (RWorkspace: need M) */
2927
2928                     i__2 = *lwork - iwork + 1;
2929                     zgebrd_(m, m, &u[u_offset], ldu, &s[1], &rwork[ie], &work[
2930                             itauq], &work[itaup], &work[iwork], &i__2, &ierr);
2931
2932 /*                 Multiply right vectors bidiagonalizing L by Q in A
2933                    (CWorkspace: need 2*M+N, prefer 2*M+N*NB)
2934                    (RWorkspace: 0) */
2935
2936                     i__2 = *lwork - iwork + 1;
2937                     zunmbr_("P", "L", "C", m, n, m, &u[u_offset], ldu, &work[
2938                             itaup], &a[a_offset], lda, &work[iwork], &i__2, &
2939                             ierr);
2940
2941 /*                 Generate left vectors bidiagonalizing L in U
2942                    (CWorkspace: need 3*M, prefer 2*M+M*NB)
2943                    (RWorkspace: 0) */
2944
2945                     i__2 = *lwork - iwork + 1;
2946                     zungbr_("Q", m, m, m, &u[u_offset], ldu, &work[itauq], &
2947                             work[iwork], &i__2, &ierr);
2948                     irwork = ie + *m;
2949
2950 /*                 Perform bidiagonal QR iteration, computing left
2951                    singular vectors of A in U and computing right
2952                    singular vectors of A in A
2953                    (CWorkspace: 0)
2954                    (RWorkspace: need BDSPAC) */
2955
2956                     zbdsqr_("U", m, n, m, &c__0, &s[1], &rwork[ie], &a[
2957                             a_offset], lda, &u[u_offset], ldu, cdum, &c__1, &
2958                             rwork[irwork], info);
2959
2960                 }
2961
2962             } else if (wntvs) {
2963
2964                 if (wntun) {
2965
2966 /*                 Path 4t(N much larger than M, JOBU='N', JOBVT='S')
2967                    M right singular vectors to be computed in VT and
2968                    no left singular vectors to be computed */
2969
2970                     if (*lwork >= *m * *m + *m * 3) {
2971
2972 /*                    Sufficient workspace for a fast algorithm */
2973
2974                         ir = 1;
2975                         if (*lwork >= wrkbl + *lda * *m) {
2976
2977 /*                       WORK(IR) is LDA by M */
2978
2979                             ldwrkr = *lda;
2980                         } else {
2981
2982 /*                       WORK(IR) is M by M */
2983
2984                             ldwrkr = *m;
2985                         }
2986                         itau = ir + ldwrkr * *m;
2987                         iwork = itau + *m;
2988
2989 /*                    Compute A=L*Q
2990                       (CWorkspace: need M*M+2*M, prefer M*M+M+M*NB)
2991                       (RWorkspace: 0) */
2992
2993                         i__2 = *lwork - iwork + 1;
2994                         zgelqf_(m, n, &a[a_offset], lda, &work[itau], &work[
2995                                 iwork], &i__2, &ierr);
2996
2997 /*                    Copy L to WORK(IR), zeroing out above it */
2998
2999                         zlacpy_("L", m, m, &a[a_offset], lda, &work[ir], &
3000                                 ldwrkr);
3001                         i__2 = *m - 1;
3002                         i__3 = *m - 1;
3003                         zlaset_("U", &i__2, &i__3, &c_b1, &c_b1, &work[ir +
3004                                 ldwrkr], &ldwrkr);
3005
3006 /*                    Generate Q in A
3007                       (CWorkspace: need M*M+2*M, prefer M*M+M+M*NB)
3008                       (RWorkspace: 0) */
3009
3010                         i__2 = *lwork - iwork + 1;
3011                         zunglq_(m, n, m, &a[a_offset], lda, &work[itau], &
3012                                 work[iwork], &i__2, &ierr);
3013                         ie = 1;
3014                         itauq = itau;
3015                         itaup = itauq + *m;
3016                         iwork = itaup + *m;
3017
3018 /*                    Bidiagonalize L in WORK(IR)
3019                       (CWorkspace: need M*M+3*M, prefer M*M+2*M+2*M*NB)
3020                       (RWorkspace: need M) */
3021
3022                         i__2 = *lwork - iwork + 1;
3023                         zgebrd_(m, m, &work[ir], &ldwrkr, &s[1], &rwork[ie], &
3024                                 work[itauq], &work[itaup], &work[iwork], &
3025                                 i__2, &ierr);
3026
3027 /*                    Generate right vectors bidiagonalizing L in
3028                       WORK(IR)
3029                       (CWorkspace: need M*M+3*M, prefer M*M+2*M+(M-1)*NB)
3030                       (RWorkspace: 0) */
3031
3032                         i__2 = *lwork - iwork + 1;
3033                         zungbr_("P", m, m, m, &work[ir], &ldwrkr, &work[itaup]
3034                                 , &work[iwork], &i__2, &ierr);
3035                         irwork = ie + *m;
3036
3037 /*                    Perform bidiagonal QR iteration, computing right
3038                       singular vectors of L in WORK(IR)
3039                       (CWorkspace: need M*M)
3040                       (RWorkspace: need BDSPAC) */
3041
3042                         zbdsqr_("U", m, m, &c__0, &c__0, &s[1], &rwork[ie], &
3043                                 work[ir], &ldwrkr, cdum, &c__1, cdum, &c__1, &
3044                                 rwork[irwork], info);
3045
3046 /*                    Multiply right singular vectors of L in WORK(IR) by
3047                       Q in A, storing result in VT
3048                       (CWorkspace: need M*M)
3049                       (RWorkspace: 0) */
3050
3051                         zgemm_("N", "N", m, n, m, &c_b2, &work[ir], &ldwrkr, &
3052                                 a[a_offset], lda, &c_b1, &vt[vt_offset], ldvt);
3053
3054                     } else {
3055
3056 /*                    Insufficient workspace for a fast algorithm */
3057
3058                         itau = 1;
3059                         iwork = itau + *m;
3060
3061 /*                    Compute A=L*Q
3062                       (CWorkspace: need 2*M, prefer M+M*NB)
3063                       (RWorkspace: 0) */
3064
3065                         i__2 = *lwork - iwork + 1;
3066                         zgelqf_(m, n, &a[a_offset], lda, &work[itau], &work[
3067                                 iwork], &i__2, &ierr);
3068
3069 /*                    Copy result to VT */
3070
3071                         zlacpy_("U", m, n, &a[a_offset], lda, &vt[vt_offset],
3072                                 ldvt);
3073
3074 /*                    Generate Q in VT
3075                       (CWorkspace: need 2*M, prefer M+M*NB)
3076                       (RWorkspace: 0) */
3077
3078                         i__2 = *lwork - iwork + 1;
3079                         zunglq_(m, n, m, &vt[vt_offset], ldvt, &work[itau], &
3080                                 work[iwork], &i__2, &ierr);
3081                         ie = 1;
3082                         itauq = itau;
3083                         itaup = itauq + *m;
3084                         iwork = itaup + *m;
3085
3086 /*                    Zero out above L in A */
3087
3088                         i__2 = *m - 1;
3089                         i__3 = *m - 1;
3090                         zlaset_("U", &i__2, &i__3, &c_b1, &c_b1, &a_ref(1, 2),
3091                                  lda);
3092
3093 /*                    Bidiagonalize L in A
3094                       (CWorkspace: need 3*M, prefer 2*M+2*M*NB)
3095                       (RWorkspace: need M) */
3096
3097                         i__2 = *lwork - iwork + 1;
3098                         zgebrd_(m, m, &a[a_offset], lda, &s[1], &rwork[ie], &
3099                                 work[itauq], &work[itaup], &work[iwork], &
3100                                 i__2, &ierr);
3101
3102 /*                    Multiply right vectors bidiagonalizing L by Q in VT
3103                       (CWorkspace: need 2*M+N, prefer 2*M+N*NB)
3104                       (RWorkspace: 0) */
3105
3106                         i__2 = *lwork - iwork + 1;
3107                         zunmbr_("P", "L", "C", m, n, m, &a[a_offset], lda, &
3108                                 work[itaup], &vt[vt_offset], ldvt, &work[
3109                                 iwork], &i__2, &ierr);
3110                         irwork = ie + *m;
3111
3112 /*                    Perform bidiagonal QR iteration, computing right
3113                       singular vectors of A in VT
3114                       (CWorkspace: 0)
3115                       (RWorkspace: need BDSPAC) */
3116
3117                         zbdsqr_("U", m, n, &c__0, &c__0, &s[1], &rwork[ie], &
3118                                 vt[vt_offset], ldvt, cdum, &c__1, cdum, &c__1,
3119                                  &rwork[irwork], info);
3120
3121                     }
3122
3123                 } else if (wntuo) {
3124
3125 /*                 Path 5t(N much larger than M, JOBU='O', JOBVT='S')
3126                    M right singular vectors to be computed in VT and
3127                    M left singular vectors to be overwritten on A */
3128
3129                     if (*lwork >= (*m << 1) * *m + *m * 3) {
3130
3131 /*                    Sufficient workspace for a fast algorithm */
3132
3133                         iu = 1;
3134                         if (*lwork >= wrkbl + (*lda << 1) * *m) {
3135
3136 /*                       WORK(IU) is LDA by M and WORK(IR) is LDA by M */
3137
3138                             ldwrku = *lda;
3139                             ir = iu + ldwrku * *m;
3140                             ldwrkr = *lda;
3141                         } else if (*lwork >= wrkbl + (*lda + *m) * *m) {
3142
3143 /*                       WORK(IU) is LDA by M and WORK(IR) is M by M */
3144
3145                             ldwrku = *lda;
3146                             ir = iu + ldwrku * *m;
3147                             ldwrkr = *m;
3148                         } else {
3149
3150 /*                       WORK(IU) is M by M and WORK(IR) is M by M */
3151
3152                             ldwrku = *m;
3153                             ir = iu + ldwrku * *m;
3154                             ldwrkr = *m;
3155                         }
3156                         itau = ir + ldwrkr * *m;
3157                         iwork = itau + *m;
3158
3159 /*                    Compute A=L*Q
3160                       (CWorkspace: need 2*M*M+2*M, prefer 2*M*M+M+M*NB)
3161                       (RWorkspace: 0) */
3162
3163                         i__2 = *lwork - iwork + 1;
3164                         zgelqf_(m, n, &a[a_offset], lda, &work[itau], &work[
3165                                 iwork], &i__2, &ierr);
3166
3167 /*                    Copy L to WORK(IU), zeroing out below it */
3168
3169                         zlacpy_("L", m, m, &a[a_offset], lda, &work[iu], &
3170                                 ldwrku);
3171                         i__2 = *m - 1;
3172                         i__3 = *m - 1;
3173                         zlaset_("U", &i__2, &i__3, &c_b1, &c_b1, &work[iu +
3174                                 ldwrku], &ldwrku);
3175
3176 /*                    Generate Q in A
3177                       (CWorkspace: need 2*M*M+2*M, prefer 2*M*M+M+M*NB)
3178                       (RWorkspace: 0) */
3179
3180                         i__2 = *lwork - iwork + 1;
3181                         zunglq_(m, n, m, &a[a_offset], lda, &work[itau], &
3182                                 work[iwork], &i__2, &ierr);
3183                         ie = 1;
3184                         itauq = itau;
3185                         itaup = itauq + *m;
3186                         iwork = itaup + *m;
3187
3188 /*                    Bidiagonalize L in WORK(IU), copying result to
3189                       WORK(IR)
3190                       (CWorkspace: need   2*M*M+3*M,
3191                                    prefer 2*M*M+2*M+2*M*NB)
3192                       (RWorkspace: need   M) */
3193
3194                         i__2 = *lwork - iwork + 1;
3195                         zgebrd_(m, m, &work[iu], &ldwrku, &s[1], &rwork[ie], &
3196                                 work[itauq], &work[itaup], &work[iwork], &
3197                                 i__2, &ierr);
3198                         zlacpy_("L", m, m, &work[iu], &ldwrku, &work[ir], &
3199                                 ldwrkr);
3200
3201 /*                    Generate right bidiagonalizing vectors in WORK(IU)
3202                       (CWorkspace: need   2*M*M+3*M-1,
3203                                    prefer 2*M*M+2*M+(M-1)*NB)
3204                       (RWorkspace: 0) */
3205
3206                         i__2 = *lwork - iwork + 1;
3207                         zungbr_("P", m, m, m, &work[iu], &ldwrku, &work[itaup]
3208                                 , &work[iwork], &i__2, &ierr);
3209
3210 /*                    Generate left bidiagonalizing vectors in WORK(IR)
3211                       (CWorkspace: need 2*M*M+3*M, prefer 2*M*M+2*M+M*NB)
3212                       (RWorkspace: 0) */
3213
3214                         i__2 = *lwork - iwork + 1;
3215                         zungbr_("Q", m, m, m, &work[ir], &ldwrkr, &work[itauq]
3216                                 , &work[iwork], &i__2, &ierr);
3217                         irwork = ie + *m;
3218
3219 /*                    Perform bidiagonal QR iteration, computing left
3220                       singular vectors of L in WORK(IR) and computing
3221                       right singular vectors of L in WORK(IU)
3222                       (CWorkspace: need 2*M*M)
3223                       (RWorkspace: need BDSPAC) */
3224
3225                         zbdsqr_("U", m, m, m, &c__0, &s[1], &rwork[ie], &work[
3226                                 iu], &ldwrku, &work[ir], &ldwrkr, cdum, &c__1,
3227                                  &rwork[irwork], info);
3228
3229 /*                    Multiply right singular vectors of L in WORK(IU) by
3230                       Q in A, storing result in VT
3231                       (CWorkspace: need M*M)
3232                       (RWorkspace: 0) */
3233
3234                         zgemm_("N", "N", m, n, m, &c_b2, &work[iu], &ldwrku, &
3235                                 a[a_offset], lda, &c_b1, &vt[vt_offset], ldvt);
3236
3237 /*                    Copy left singular vectors of L to A
3238                       (CWorkspace: need M*M)
3239                       (RWorkspace: 0) */
3240
3241                         zlacpy_("F", m, m, &work[ir], &ldwrkr, &a[a_offset],
3242                                 lda);
3243
3244                     } else {
3245
3246 /*                    Insufficient workspace for a fast algorithm */
3247
3248                         itau = 1;
3249                         iwork = itau + *m;
3250
3251 /*                    Compute A=L*Q, copying result to VT
3252                       (CWorkspace: need 2*M, prefer M+M*NB)
3253                       (RWorkspace: 0) */
3254
3255                         i__2 = *lwork - iwork + 1;
3256                         zgelqf_(m, n, &a[a_offset], lda, &work[itau], &work[
3257                                 iwork], &i__2, &ierr);
3258                         zlacpy_("U", m, n, &a[a_offset], lda, &vt[vt_offset],
3259                                 ldvt);
3260
3261 /*                    Generate Q in VT
3262                       (CWorkspace: need 2*M, prefer M+M*NB)
3263                       (RWorkspace: 0) */
3264
3265                         i__2 = *lwork - iwork + 1;
3266                         zunglq_(m, n, m, &vt[vt_offset], ldvt, &work[itau], &
3267                                 work[iwork], &i__2, &ierr);
3268                         ie = 1;
3269                         itauq = itau;
3270                         itaup = itauq + *m;
3271                         iwork = itaup + *m;
3272
3273 /*                    Zero out above L in A */
3274
3275                         i__2 = *m - 1;
3276                         i__3 = *m - 1;
3277                         zlaset_("U", &i__2, &i__3, &c_b1, &c_b1, &a_ref(1, 2),
3278                                  lda);
3279
3280 /*                    Bidiagonalize L in A
3281                       (CWorkspace: need 3*M, prefer 2*M+2*M*NB)
3282                       (RWorkspace: need M) */
3283
3284                         i__2 = *lwork - iwork + 1;
3285                         zgebrd_(m, m, &a[a_offset], lda, &s[1], &rwork[ie], &
3286                                 work[itauq], &work[itaup], &work[iwork], &
3287                                 i__2, &ierr);
3288
3289 /*                    Multiply right vectors bidiagonalizing L by Q in VT
3290                       (CWorkspace: need 2*M+N, prefer 2*M+N*NB)
3291                       (RWorkspace: 0) */
3292
3293                         i__2 = *lwork - iwork + 1;
3294                         zunmbr_("P", "L", "C", m, n, m, &a[a_offset], lda, &
3295                                 work[itaup], &vt[vt_offset], ldvt, &work[
3296                                 iwork], &i__2, &ierr);
3297
3298 /*                    Generate left bidiagonalizing vectors of L in A
3299                       (CWorkspace: need 3*M, prefer 2*M+M*NB)
3300                       (RWorkspace: 0) */
3301
3302                         i__2 = *lwork - iwork + 1;
3303                         zungbr_("Q", m, m, m, &a[a_offset], lda, &work[itauq],
3304                                  &work[iwork], &i__2, &ierr);
3305                         irwork = ie + *m;
3306
3307 /*                    Perform bidiagonal QR iteration, computing left
3308                       singular vectors of A in A and computing right
3309                       singular vectors of A in VT
3310                       (CWorkspace: 0)
3311                       (RWorkspace: need BDSPAC) */
3312
3313                         zbdsqr_("U", m, n, m, &c__0, &s[1], &rwork[ie], &vt[
3314                                 vt_offset], ldvt, &a[a_offset], lda, cdum, &
3315                                 c__1, &rwork[irwork], info);
3316
3317                     }
3318
3319                 } else if (wntuas) {
3320
3321 /*                 Path 6t(N much larger than M, JOBU='S' or 'A',
3322                            JOBVT='S')
3323                    M right singular vectors to be computed in VT and
3324                    M left singular vectors to be computed in U */
3325
3326                     if (*lwork >= *m * *m + *m * 3) {
3327
3328 /*                    Sufficient workspace for a fast algorithm */
3329
3330                         iu = 1;
3331                         if (*lwork >= wrkbl + *lda * *m) {
3332
3333 /*                       WORK(IU) is LDA by N */
3334
3335                             ldwrku = *lda;
3336                         } else {
3337
3338 /*                       WORK(IU) is LDA by M */
3339
3340                             ldwrku = *m;
3341                         }
3342                         itau = iu + ldwrku * *m;
3343                         iwork = itau + *m;
3344
3345 /*                    Compute A=L*Q
3346                       (CWorkspace: need M*M+2*M, prefer M*M+M+M*NB)
3347                       (RWorkspace: 0) */
3348
3349                         i__2 = *lwork - iwork + 1;
3350                         zgelqf_(m, n, &a[a_offset], lda, &work[itau], &work[
3351                                 iwork], &i__2, &ierr);
3352
3353 /*                    Copy L to WORK(IU), zeroing out above it */
3354
3355                         zlacpy_("L", m, m, &a[a_offset], lda, &work[iu], &
3356                                 ldwrku);
3357                         i__2 = *m - 1;
3358                         i__3 = *m - 1;
3359                         zlaset_("U", &i__2, &i__3, &c_b1, &c_b1, &work[iu +
3360                                 ldwrku], &ldwrku);
3361
3362 /*                    Generate Q in A
3363                       (CWorkspace: need M*M+2*M, prefer M*M+M+M*NB)
3364                       (RWorkspace: 0) */
3365
3366                         i__2 = *lwork - iwork + 1;
3367                         zunglq_(m, n, m, &a[a_offset], lda, &work[itau], &
3368                                 work[iwork], &i__2, &ierr);
3369                         ie = 1;
3370                         itauq = itau;
3371                         itaup = itauq + *m;
3372                         iwork = itaup + *m;
3373
3374 /*                    Bidiagonalize L in WORK(IU), copying result to U
3375                       (CWorkspace: need M*M+3*M, prefer M*M+2*M+2*M*NB)
3376                       (RWorkspace: need M) */
3377
3378                         i__2 = *lwork - iwork + 1;
3379                         zgebrd_(m, m, &work[iu], &ldwrku, &s[1], &rwork[ie], &
3380                                 work[itauq], &work[itaup], &work[iwork], &
3381                                 i__2, &ierr);
3382                         zlacpy_("L", m, m, &work[iu], &ldwrku, &u[u_offset],
3383                                 ldu);
3384
3385 /*                    Generate right bidiagonalizing vectors in WORK(IU)
3386                       (CWorkspace: need   M*M+3*M-1,
3387                                    prefer M*M+2*M+(M-1)*NB)
3388                       (RWorkspace: 0) */
3389
3390                         i__2 = *lwork - iwork + 1;
3391                         zungbr_("P", m, m, m, &work[iu], &ldwrku, &work[itaup]
3392                                 , &work[iwork], &i__2, &ierr);
3393
3394 /*                    Generate left bidiagonalizing vectors in U
3395                       (CWorkspace: need M*M+3*M, prefer M*M+2*M+M*NB)
3396                       (RWorkspace: 0) */
3397
3398                         i__2 = *lwork - iwork + 1;
3399                         zungbr_("Q", m, m, m, &u[u_offset], ldu, &work[itauq],
3400                                  &work[iwork], &i__2, &ierr);
3401                         irwork = ie + *m;
3402
3403 /*                    Perform bidiagonal QR iteration, computing left
3404                       singular vectors of L in U and computing right
3405                       singular vectors of L in WORK(IU)
3406                       (CWorkspace: need M*M)
3407                       (RWorkspace: need BDSPAC) */
3408
3409                         zbdsqr_("U", m, m, m, &c__0, &s[1], &rwork[ie], &work[
3410                                 iu], &ldwrku, &u[u_offset], ldu, cdum, &c__1,
3411                                 &rwork[irwork], info);
3412
3413 /*                    Multiply right singular vectors of L in WORK(IU) by
3414                       Q in A, storing result in VT
3415                       (CWorkspace: need M*M)
3416                       (RWorkspace: 0) */
3417
3418                         zgemm_("N", "N", m, n, m, &c_b2, &work[iu], &ldwrku, &
3419                                 a[a_offset], lda, &c_b1, &vt[vt_offset], ldvt);
3420
3421                     } else {
3422
3423 /*                    Insufficient workspace for a fast algorithm */
3424
3425                         itau = 1;
3426                         iwork = itau + *m;
3427
3428 /*                    Compute A=L*Q, copying result to VT
3429                       (CWorkspace: need 2*M, prefer M+M*NB)
3430                       (RWorkspace: 0) */
3431
3432                         i__2 = *lwork - iwork + 1;
3433                         zgelqf_(m, n, &a[a_offset], lda, &work[itau], &work[
3434                                 iwork], &i__2, &ierr);
3435                         zlacpy_("U", m, n, &a[a_offset], lda, &vt[vt_offset],
3436                                 ldvt);
3437
3438 /*                    Generate Q in VT
3439                       (CWorkspace: need 2*M, prefer M+M*NB)
3440                       (RWorkspace: 0) */
3441
3442                         i__2 = *lwork - iwork + 1;
3443                         zunglq_(m, n, m, &vt[vt_offset], ldvt, &work[itau], &
3444                                 work[iwork], &i__2, &ierr);
3445
3446 /*                    Copy L to U, zeroing out above it */
3447
3448                         zlacpy_("L", m, m, &a[a_offset], lda, &u[u_offset],
3449                                 ldu);
3450                         i__2 = *m - 1;
3451                         i__3 = *m - 1;
3452                         zlaset_("U", &i__2, &i__3, &c_b1, &c_b1, &u_ref(1, 2),
3453                                  ldu);
3454                         ie = 1;
3455                         itauq = itau;
3456                         itaup = itauq + *m;
3457                         iwork = itaup + *m;
3458
3459 /*                    Bidiagonalize L in U
3460                       (CWorkspace: need 3*M, prefer 2*M+2*M*NB)
3461                       (RWorkspace: need M) */
3462
3463                         i__2 = *lwork - iwork + 1;
3464                         zgebrd_(m, m, &u[u_offset], ldu, &s[1], &rwork[ie], &
3465                                 work[itauq], &work[itaup], &work[iwork], &
3466                                 i__2, &ierr);
3467
3468 /*                    Multiply right bidiagonalizing vectors in U by Q
3469                       in VT
3470                       (CWorkspace: need 2*M+N, prefer 2*M+N*NB)
3471                       (RWorkspace: 0) */
3472
3473                         i__2 = *lwork - iwork + 1;
3474                         zunmbr_("P", "L", "C", m, n, m, &u[u_offset], ldu, &
3475                                 work[itaup], &vt[vt_offset], ldvt, &work[
3476                                 iwork], &i__2, &ierr);
3477
3478 /*                    Generate left bidiagonalizing vectors in U
3479                       (CWorkspace: need 3*M, prefer 2*M+M*NB)
3480                       (RWorkspace: 0) */
3481
3482                         i__2 = *lwork - iwork + 1;
3483                         zungbr_("Q", m, m, m, &u[u_offset], ldu, &work[itauq],
3484                                  &work[iwork], &i__2, &ierr);
3485                         irwork = ie + *m;
3486
3487 /*                    Perform bidiagonal QR iteration, computing left
3488                       singular vectors of A in U and computing right
3489                       singular vectors of A in VT
3490                       (CWorkspace: 0)
3491                       (RWorkspace: need BDSPAC) */
3492
3493                         zbdsqr_("U", m, n, m, &c__0, &s[1], &rwork[ie], &vt[
3494                                 vt_offset], ldvt, &u[u_offset], ldu, cdum, &
3495                                 c__1, &rwork[irwork], info);
3496
3497                     }
3498
3499                 }
3500
3501             } else if (wntva) {
3502
3503                 if (wntun) {
3504
3505 /*                 Path 7t(N much larger than M, JOBU='N', JOBVT='A')
3506                    N right singular vectors to be computed in VT and
3507                    no left singular vectors to be computed
3508
3509    Computing MAX */
3510                     i__2 = *n + *m, i__3 = *m * 3;
3511                     if (*lwork >= *m * *m + max(i__2,i__3)) {
3512
3513 /*                    Sufficient workspace for a fast algorithm */
3514
3515                         ir = 1;
3516                         if (*lwork >= wrkbl + *lda * *m) {
3517
3518 /*                       WORK(IR) is LDA by M */
3519
3520                             ldwrkr = *lda;
3521                         } else {
3522
3523 /*                       WORK(IR) is M by M */
3524
3525                             ldwrkr = *m;
3526                         }
3527                         itau = ir + ldwrkr * *m;
3528                         iwork = itau + *m;
3529
3530 /*                    Compute A=L*Q, copying result to VT
3531                       (CWorkspace: need M*M+2*M, prefer M*M+M+M*NB)
3532                       (RWorkspace: 0) */
3533
3534                         i__2 = *lwork - iwork + 1;
3535                         zgelqf_(m, n, &a[a_offset], lda, &work[itau], &work[
3536                                 iwork], &i__2, &ierr);
3537                         zlacpy_("U", m, n, &a[a_offset], lda, &vt[vt_offset],
3538                                 ldvt);
3539
3540 /*                    Copy L to WORK(IR), zeroing out above it */
3541
3542                         zlacpy_("L", m, m, &a[a_offset], lda, &work[ir], &
3543                                 ldwrkr);
3544                         i__2 = *m - 1;
3545                         i__3 = *m - 1;
3546                         zlaset_("U", &i__2, &i__3, &c_b1, &c_b1, &work[ir +
3547                                 ldwrkr], &ldwrkr);
3548
3549 /*                    Generate Q in VT
3550                       (CWorkspace: need M*M+M+N, prefer M*M+M+N*NB)
3551                       (RWorkspace: 0) */
3552
3553                         i__2 = *lwork - iwork + 1;
3554                         zunglq_(n, n, m, &vt[vt_offset], ldvt, &work[itau], &
3555                                 work[iwork], &i__2, &ierr);
3556                         ie = 1;
3557                         itauq = itau;
3558                         itaup = itauq + *m;
3559                         iwork = itaup + *m;
3560
3561 /*                    Bidiagonalize L in WORK(IR)
3562                       (CWorkspace: need M*M+3*M, prefer M*M+2*M+2*M*NB)
3563                       (RWorkspace: need M) */
3564
3565                         i__2 = *lwork - iwork + 1;
3566                         zgebrd_(m, m, &work[ir], &ldwrkr, &s[1], &rwork[ie], &
3567                                 work[itauq], &work[itaup], &work[iwork], &
3568                                 i__2, &ierr);
3569
3570 /*                    Generate right bidiagonalizing vectors in WORK(IR)
3571                       (CWorkspace: need   M*M+3*M-1,
3572                                    prefer M*M+2*M+(M-1)*NB)
3573                       (RWorkspace: 0) */
3574
3575                         i__2 = *lwork - iwork + 1;
3576                         zungbr_("P", m, m, m, &work[ir], &ldwrkr, &work[itaup]
3577                                 , &work[iwork], &i__2, &ierr);
3578                         irwork = ie + *m;
3579
3580 /*                    Perform bidiagonal QR iteration, computing right
3581                       singular vectors of L in WORK(IR)
3582                       (CWorkspace: need M*M)
3583                       (RWorkspace: need BDSPAC) */
3584
3585                         zbdsqr_("U", m, m, &c__0, &c__0, &s[1], &rwork[ie], &
3586                                 work[ir], &ldwrkr, cdum, &c__1, cdum, &c__1, &
3587                                 rwork[irwork], info);
3588
3589 /*                    Multiply right singular vectors of L in WORK(IR) by
3590                       Q in VT, storing result in A
3591                       (CWorkspace: need M*M)
3592                       (RWorkspace: 0) */
3593
3594                         zgemm_("N", "N", m, n, m, &c_b2, &work[ir], &ldwrkr, &
3595                                 vt[vt_offset], ldvt, &c_b1, &a[a_offset], lda);
3596
3597 /*                    Copy right singular vectors of A from A to VT */
3598
3599                         zlacpy_("F", m, n, &a[a_offset], lda, &vt[vt_offset],
3600                                 ldvt);
3601
3602                     } else {
3603
3604 /*                    Insufficient workspace for a fast algorithm */
3605
3606                         itau = 1;
3607                         iwork = itau + *m;
3608
3609 /*                    Compute A=L*Q, copying result to VT
3610                       (CWorkspace: need 2*M, prefer M+M*NB)
3611                       (RWorkspace: 0) */
3612
3613                         i__2 = *lwork - iwork + 1;
3614                         zgelqf_(m, n, &a[a_offset], lda, &work[itau], &work[
3615                                 iwork], &i__2, &ierr);
3616                         zlacpy_("U", m, n, &a[a_offset], lda, &vt[vt_offset],
3617                                 ldvt);
3618
3619 /*                    Generate Q in VT
3620                       (CWorkspace: need M+N, prefer M+N*NB)
3621                       (RWorkspace: 0) */
3622
3623                         i__2 = *lwork - iwork + 1;
3624                         zunglq_(n, n, m, &vt[vt_offset], ldvt, &work[itau], &
3625                                 work[iwork], &i__2, &ierr);
3626                         ie = 1;
3627                         itauq = itau;
3628                         itaup = itauq + *m;
3629                         iwork = itaup + *m;
3630
3631 /*                    Zero out above L in A */
3632
3633                         i__2 = *m - 1;
3634                         i__3 = *m - 1;
3635                         zlaset_("U", &i__2, &i__3, &c_b1, &c_b1, &a_ref(1, 2),
3636                                  lda);
3637
3638 /*                    Bidiagonalize L in A
3639                       (CWorkspace: need 3*M, prefer 2*M+2*M*NB)
3640                       (RWorkspace: need M) */
3641
3642                         i__2 = *lwork - iwork + 1;
3643                         zgebrd_(m, m, &a[a_offset], lda, &s[1], &rwork[ie], &
3644                                 work[itauq], &work[itaup], &work[iwork], &
3645                                 i__2, &ierr);
3646
3647 /*                    Multiply right bidiagonalizing vectors in A by Q
3648                       in VT
3649                       (CWorkspace: need 2*M+N, prefer 2*M+N*NB)
3650                       (RWorkspace: 0) */
3651
3652                         i__2 = *lwork - iwork + 1;
3653                         zunmbr_("P", "L", "C", m, n, m, &a[a_offset], lda, &
3654                                 work[itaup], &vt[vt_offset], ldvt, &work[
3655                                 iwork], &i__2, &ierr);
3656                         irwork = ie + *m;
3657
3658 /*                    Perform bidiagonal QR iteration, computing right
3659                       singular vectors of A in VT
3660                       (CWorkspace: 0)
3661                       (RWorkspace: need BDSPAC) */
3662
3663                         zbdsqr_("U", m, n, &c__0, &c__0, &s[1], &rwork[ie], &
3664                                 vt[vt_offset], ldvt, cdum, &c__1, cdum, &c__1,
3665                                  &rwork[irwork], info);
3666
3667                     }
3668
3669                 } else if (wntuo) {
3670
3671 /*                 Path 8t(N much larger than M, JOBU='O', JOBVT='A')
3672                    N right singular vectors to be computed in VT and
3673                    M left singular vectors to be overwritten on A
3674
3675    Computing MAX */
3676                     i__2 = *n + *m, i__3 = *m * 3;
3677                     if (*lwork >= (*m << 1) * *m + max(i__2,i__3)) {
3678
3679 /*                    Sufficient workspace for a fast algorithm */
3680
3681                         iu = 1;
3682                         if (*lwork >= wrkbl + (*lda << 1) * *m) {
3683
3684 /*                       WORK(IU) is LDA by M and WORK(IR) is LDA by M */
3685
3686                             ldwrku = *lda;
3687                             ir = iu + ldwrku * *m;
3688                             ldwrkr = *lda;
3689                         } else if (*lwork >= wrkbl + (*lda + *m) * *m) {
3690
3691 /*                       WORK(IU) is LDA by M and WORK(IR) is M by M */
3692
3693                             ldwrku = *lda;
3694                             ir = iu + ldwrku * *m;
3695                             ldwrkr = *m;
3696                         } else {
3697
3698 /*                       WORK(IU) is M by M and WORK(IR) is M by M */
3699
3700                             ldwrku = *m;
3701                             ir = iu + ldwrku * *m;
3702                             ldwrkr = *m;
3703                         }
3704                         itau = ir + ldwrkr * *m;
3705                         iwork = itau + *m;
3706
3707 /*                    Compute A=L*Q, copying result to VT
3708                       (CWorkspace: need 2*M*M+2*M, prefer 2*M*M+M+M*NB)
3709                       (RWorkspace: 0) */
3710
3711                         i__2 = *lwork - iwork + 1;
3712                         zgelqf_(m, n, &a[a_offset], lda, &work[itau], &work[
3713                                 iwork], &i__2, &ierr);
3714                         zlacpy_("U", m, n, &a[a_offset], lda, &vt[vt_offset],
3715                                 ldvt);
3716
3717 /*                    Generate Q in VT
3718                       (CWorkspace: need 2*M*M+M+N, prefer 2*M*M+M+N*NB)
3719                       (RWorkspace: 0) */
3720
3721                         i__2 = *lwork - iwork + 1;
3722                         zunglq_(n, n, m, &vt[vt_offset], ldvt, &work[itau], &
3723                                 work[iwork], &i__2, &ierr);
3724
3725 /*                    Copy L to WORK(IU), zeroing out above it */
3726
3727                         zlacpy_("L", m, m, &a[a_offset], lda, &work[iu], &
3728                                 ldwrku);
3729                         i__2 = *m - 1;
3730                         i__3 = *m - 1;
3731                         zlaset_("U", &i__2, &i__3, &c_b1, &c_b1, &work[iu +
3732                                 ldwrku], &ldwrku);
3733                         ie = 1;
3734                         itauq = itau;
3735                         itaup = itauq + *m;
3736                         iwork = itaup + *m;
3737
3738 /*                    Bidiagonalize L in WORK(IU), copying result to
3739                       WORK(IR)
3740                       (CWorkspace: need   2*M*M+3*M,
3741                                    prefer 2*M*M+2*M+2*M*NB)
3742                       (RWorkspace: need   M) */
3743
3744                         i__2 = *lwork - iwork + 1;
3745                         zgebrd_(m, m, &work[iu], &ldwrku, &s[1], &rwork[ie], &
3746                                 work[itauq], &work[itaup], &work[iwork], &
3747                                 i__2, &ierr);
3748                         zlacpy_("L", m, m, &work[iu], &ldwrku, &work[ir], &
3749                                 ldwrkr);
3750
3751 /*                    Generate right bidiagonalizing vectors in WORK(IU)
3752                       (CWorkspace: need   2*M*M+3*M-1,
3753                                    prefer 2*M*M+2*M+(M-1)*NB)
3754                       (RWorkspace: 0) */
3755
3756                         i__2 = *lwork - iwork + 1;
3757                         zungbr_("P", m, m, m, &work[iu], &ldwrku, &work[itaup]
3758                                 , &work[iwork], &i__2, &ierr);
3759
3760 /*                    Generate left bidiagonalizing vectors in WORK(IR)
3761                       (CWorkspace: need 2*M*M+3*M, prefer 2*M*M+2*M+M*NB)
3762                       (RWorkspace: 0) */
3763
3764                         i__2 = *lwork - iwork + 1;
3765                         zungbr_("Q", m, m, m, &work[ir], &ldwrkr, &work[itauq]
3766                                 , &work[iwork], &i__2, &ierr);
3767                         irwork = ie + *m;
3768
3769 /*                    Perform bidiagonal QR iteration, computing left
3770                       singular vectors of L in WORK(IR) and computing
3771                       right singular vectors of L in WORK(IU)
3772                       (CWorkspace: need 2*M*M)
3773                       (RWorkspace: need BDSPAC) */
3774
3775                         zbdsqr_("U", m, m, m, &c__0, &s[1], &rwork[ie], &work[
3776                                 iu], &ldwrku, &work[ir], &ldwrkr, cdum, &c__1,
3777                                  &rwork[irwork], info);
3778
3779 /*                    Multiply right singular vectors of L in WORK(IU) by
3780                       Q in VT, storing result in A
3781                       (CWorkspace: need M*M)
3782                       (RWorkspace: 0) */
3783
3784                         zgemm_("N", "N", m, n, m, &c_b2, &work[iu], &ldwrku, &
3785                                 vt[vt_offset], ldvt, &c_b1, &a[a_offset], lda);
3786
3787 /*                    Copy right singular vectors of A from A to VT */
3788
3789                         zlacpy_("F", m, n, &a[a_offset], lda, &vt[vt_offset],
3790                                 ldvt);
3791
3792 /*                    Copy left singular vectors of A from WORK(IR) to A */
3793
3794                         zlacpy_("F", m, m, &work[ir], &ldwrkr, &a[a_offset],
3795                                 lda);
3796
3797                     } else {
3798
3799 /*                    Insufficient workspace for a fast algorithm */
3800
3801                         itau = 1;
3802                         iwork = itau + *m;
3803
3804 /*                    Compute A=L*Q, copying result to VT
3805                       (CWorkspace: need 2*M, prefer M+M*NB)
3806                       (RWorkspace: 0) */
3807
3808                         i__2 = *lwork - iwork + 1;
3809                         zgelqf_(m, n, &a[a_offset], lda, &work[itau], &work[
3810                                 iwork], &i__2, &ierr);
3811                         zlacpy_("U", m, n, &a[a_offset], lda, &vt[vt_offset],
3812                                 ldvt);
3813
3814 /*                    Generate Q in VT
3815                       (CWorkspace: need M+N, prefer M+N*NB)
3816                       (RWorkspace: 0) */
3817
3818                         i__2 = *lwork - iwork + 1;
3819                         zunglq_(n, n, m, &vt[vt_offset], ldvt, &work[itau], &
3820                                 work[iwork], &i__2, &ierr);
3821                         ie = 1;
3822                         itauq = itau;
3823                         itaup = itauq + *m;
3824                         iwork = itaup + *m;
3825
3826 /*                    Zero out above L in A */
3827
3828                         i__2 = *m - 1;
3829                         i__3 = *m - 1;
3830                         zlaset_("U", &i__2, &i__3, &c_b1, &c_b1, &a_ref(1, 2),
3831                                  lda);
3832
3833 /*                    Bidiagonalize L in A
3834                       (CWorkspace: need 3*M, prefer 2*M+2*M*NB)
3835                       (RWorkspace: need M) */
3836
3837                         i__2 = *lwork - iwork + 1;
3838                         zgebrd_(m, m, &a[a_offset], lda, &s[1], &rwork[ie], &
3839                                 work[itauq], &work[itaup], &work[iwork], &
3840                                 i__2, &ierr);
3841
3842 /*                    Multiply right bidiagonalizing vectors in A by Q
3843                       in VT
3844                       (CWorkspace: need 2*M+N, prefer 2*M+N*NB)
3845                       (RWorkspace: 0) */
3846
3847                         i__2 = *lwork - iwork + 1;
3848                         zunmbr_("P", "L", "C", m, n, m, &a[a_offset], lda, &
3849                                 work[itaup], &vt[vt_offset], ldvt, &work[
3850                                 iwork], &i__2, &ierr);
3851
3852 /*                    Generate left bidiagonalizing vectors in A
3853                       (CWorkspace: need 3*M, prefer 2*M+M*NB)
3854                       (RWorkspace: 0) */
3855
3856                         i__2 = *lwork - iwork + 1;
3857                         zungbr_("Q", m, m, m, &a[a_offset], lda, &work[itauq],
3858                                  &work[iwork], &i__2, &ierr);
3859                         irwork = ie + *m;
3860
3861 /*                    Perform bidiagonal QR iteration, computing left
3862                       singular vectors of A in A and computing right
3863                       singular vectors of A in VT
3864                       (CWorkspace: 0)
3865                       (RWorkspace: need BDSPAC) */
3866
3867                         zbdsqr_("U", m, n, m, &c__0, &s[1], &rwork[ie], &vt[
3868                                 vt_offset], ldvt, &a[a_offset], lda, cdum, &
3869                                 c__1, &rwork[irwork], info);
3870
3871                     }
3872
3873                 } else if (wntuas) {
3874
3875 /*                 Path 9t(N much larger than M, JOBU='S' or 'A',
3876                            JOBVT='A')
3877                    N right singular vectors to be computed in VT and
3878                    M left singular vectors to be computed in U
3879
3880    Computing MAX */
3881                     i__2 = *n + *m, i__3 = *m * 3;
3882                     if (*lwork >= *m * *m + max(i__2,i__3)) {
3883
3884 /*                    Sufficient workspace for a fast algorithm */
3885
3886                         iu = 1;
3887                         if (*lwork >= wrkbl + *lda * *m) {
3888
3889 /*                       WORK(IU) is LDA by M */
3890
3891                             ldwrku = *lda;
3892                         } else {
3893
3894 /*                       WORK(IU) is M by M */
3895
3896                             ldwrku = *m;
3897                         }
3898                         itau = iu + ldwrku * *m;
3899                         iwork = itau + *m;
3900
3901 /*                    Compute A=L*Q, copying result to VT
3902                       (CWorkspace: need M*M+2*M, prefer M*M+M+M*NB)
3903                       (RWorkspace: 0) */
3904
3905                         i__2 = *lwork - iwork + 1;
3906                         zgelqf_(m, n, &a[a_offset], lda, &work[itau], &work[
3907                                 iwork], &i__2, &ierr);
3908                         zlacpy_("U", m, n, &a[a_offset], lda, &vt[vt_offset],
3909                                 ldvt);
3910
3911 /*                    Generate Q in VT
3912                       (CWorkspace: need M*M+M+N, prefer M*M+M+N*NB)
3913                       (RWorkspace: 0) */
3914
3915                         i__2 = *lwork - iwork + 1;
3916                         zunglq_(n, n, m, &vt[vt_offset], ldvt, &work[itau], &
3917                                 work[iwork], &i__2, &ierr);
3918
3919 /*                    Copy L to WORK(IU), zeroing out above it */
3920
3921                         zlacpy_("L", m, m, &a[a_offset], lda, &work[iu], &
3922                                 ldwrku);
3923                         i__2 = *m - 1;
3924                         i__3 = *m - 1;
3925                         zlaset_("U", &i__2, &i__3, &c_b1, &c_b1, &work[iu +
3926                                 ldwrku], &ldwrku);
3927                         ie = 1;
3928                         itauq = itau;
3929                         itaup = itauq + *m;
3930                         iwork = itaup + *m;
3931
3932 /*                    Bidiagonalize L in WORK(IU), copying result to U
3933                       (CWorkspace: need M*M+3*M, prefer M*M+2*M+2*M*NB)
3934                       (RWorkspace: need M) */
3935
3936                         i__2 = *lwork - iwork + 1;
3937                         zgebrd_(m, m, &work[iu], &ldwrku, &s[1], &rwork[ie], &
3938                                 work[itauq], &work[itaup], &work[iwork], &
3939                                 i__2, &ierr);
3940                         zlacpy_("L", m, m, &work[iu], &ldwrku, &u[u_offset],
3941                                 ldu);
3942
3943 /*                    Generate right bidiagonalizing vectors in WORK(IU)
3944                       (CWorkspace: need M*M+3*M, prefer M*M+2*M+(M-1)*NB)
3945                       (RWorkspace: 0) */
3946
3947                         i__2 = *lwork - iwork + 1;
3948                         zungbr_("P", m, m, m, &work[iu], &ldwrku, &work[itaup]
3949                                 , &work[iwork], &i__2, &ierr);
3950
3951 /*                    Generate left bidiagonalizing vectors in U
3952                       (CWorkspace: need M*M+3*M, prefer M*M+2*M+M*NB)
3953                       (RWorkspace: 0) */
3954
3955                         i__2 = *lwork - iwork + 1;
3956                         zungbr_("Q", m, m, m, &u[u_offset], ldu, &work[itauq],
3957                                  &work[iwork], &i__2, &ierr);
3958                         irwork = ie + *m;
3959
3960 /*                    Perform bidiagonal QR iteration, computing left
3961                       singular vectors of L in U and computing right
3962                       singular vectors of L in WORK(IU)
3963                       (CWorkspace: need M*M)
3964                       (RWorkspace: need BDSPAC) */
3965
3966                         zbdsqr_("U", m, m, m, &c__0, &s[1], &rwork[ie], &work[
3967                                 iu], &ldwrku, &u[u_offset], ldu, cdum, &c__1,
3968                                 &rwork[irwork], info);
3969
3970 /*                    Multiply right singular vectors of L in WORK(IU) by
3971                       Q in VT, storing result in A
3972                       (CWorkspace: need M*M)
3973                       (RWorkspace: 0) */
3974
3975                         zgemm_("N", "N", m, n, m, &c_b2, &work[iu], &ldwrku, &
3976                                 vt[vt_offset], ldvt, &c_b1, &a[a_offset], lda);
3977
3978 /*                    Copy right singular vectors of A from A to VT */
3979
3980                         zlacpy_("F", m, n, &a[a_offset], lda, &vt[vt_offset],
3981                                 ldvt);
3982
3983                     } else {
3984
3985 /*                    Insufficient workspace for a fast algorithm */
3986
3987                         itau = 1;
3988                         iwork = itau + *m;
3989
3990 /*                    Compute A=L*Q, copying result to VT
3991                       (CWorkspace: need 2*M, prefer M+M*NB)
3992                       (RWorkspace: 0) */
3993
3994                         i__2 = *lwork - iwork + 1;
3995                         zgelqf_(m, n, &a[a_offset], lda, &work[itau], &work[
3996                                 iwork], &i__2, &ierr);
3997                         zlacpy_("U", m, n, &a[a_offset], lda, &vt[vt_offset],
3998                                 ldvt);
3999
4000 /*                    Generate Q in VT
4001                       (CWorkspace: need M+N, prefer M+N*NB)
4002                       (RWorkspace: 0) */
4003
4004                         i__2 = *lwork - iwork + 1;
4005                         zunglq_(n, n, m, &vt[vt_offset], ldvt, &work[itau], &
4006                                 work[iwork], &i__2, &ierr);
4007
4008 /*                    Copy L to U, zeroing out above it */
4009
4010                         zlacpy_("L", m, m, &a[a_offset], lda, &u[u_offset],
4011                                 ldu);
4012                         i__2 = *m - 1;
4013                         i__3 = *m - 1;
4014                         zlaset_("U", &i__2, &i__3, &c_b1, &c_b1, &u_ref(1, 2),
4015                                  ldu);
4016                         ie = 1;
4017                         itauq = itau;
4018                         itaup = itauq + *m;
4019                         iwork = itaup + *m;
4020
4021 /*                    Bidiagonalize L in U
4022                       (CWorkspace: need 3*M, prefer 2*M+2*M*NB)
4023                       (RWorkspace: need M) */
4024
4025                         i__2 = *lwork - iwork + 1;
4026                         zgebrd_(m, m, &u[u_offset], ldu, &s[1], &rwork[ie], &
4027                                 work[itauq], &work[itaup], &work[iwork], &
4028                                 i__2, &ierr);
4029
4030 /*                    Multiply right bidiagonalizing vectors in U by Q
4031                       in VT
4032                       (CWorkspace: need 2*M+N, prefer 2*M+N*NB)
4033                       (RWorkspace: 0) */
4034
4035                         i__2 = *lwork - iwork + 1;
4036                         zunmbr_("P", "L", "C", m, n, m, &u[u_offset], ldu, &
4037                                 work[itaup], &vt[vt_offset], ldvt, &work[
4038                                 iwork], &i__2, &ierr);
4039
4040 /*                    Generate left bidiagonalizing vectors in U
4041                       (CWorkspace: need 3*M, prefer 2*M+M*NB)
4042                       (RWorkspace: 0) */
4043
4044                         i__2 = *lwork - iwork + 1;
4045                         zungbr_("Q", m, m, m, &u[u_offset], ldu, &work[itauq],
4046                                  &work[iwork], &i__2, &ierr);
4047                         irwork = ie + *m;
4048
4049 /*                    Perform bidiagonal QR iteration, computing left
4050                       singular vectors of A in U and computing right
4051                       singular vectors of A in VT
4052                       (CWorkspace: 0)
4053                       (RWorkspace: need BDSPAC) */
4054
4055                         zbdsqr_("U", m, n, m, &c__0, &s[1], &rwork[ie], &vt[
4056                                 vt_offset], ldvt, &u[u_offset], ldu, cdum, &
4057                                 c__1, &rwork[irwork], info);
4058
4059                     }
4060
4061                 }
4062
4063             }
4064
4065         } else {
4066
4067 /*           N .LT. MNTHR
4068
4069              Path 10t(N greater than M, but not much larger)
4070              Reduce to bidiagonal form without LQ decomposition */
4071
4072             ie = 1;
4073             itauq = 1;
4074             itaup = itauq + *m;
4075             iwork = itaup + *m;
4076
4077 /*           Bidiagonalize A
4078              (CWorkspace: need 2*M+N, prefer 2*M+(M+N)*NB)
4079              (RWorkspace: M) */
4080
4081             i__2 = *lwork - iwork + 1;
4082             zgebrd_(m, n, &a[a_offset], lda, &s[1], &rwork[ie], &work[itauq],
4083                     &work[itaup], &work[iwork], &i__2, &ierr);
4084             if (wntuas) {
4085
4086 /*              If left singular vectors desired in U, copy result to U
4087                 and generate left bidiagonalizing vectors in U
4088                 (CWorkspace: need 3*M-1, prefer 2*M+(M-1)*NB)
4089                 (RWorkspace: 0) */
4090
4091                 zlacpy_("L", m, m, &a[a_offset], lda, &u[u_offset], ldu);
4092                 i__2 = *lwork - iwork + 1;
4093                 zungbr_("Q", m, m, n, &u[u_offset], ldu, &work[itauq], &work[
4094                         iwork], &i__2, &ierr);
4095             }
4096             if (wntvas) {
4097
4098 /*              If right singular vectors desired in VT, copy result to
4099                 VT and generate right bidiagonalizing vectors in VT
4100                 (CWorkspace: need 2*M+NRVT, prefer 2*M+NRVT*NB)
4101                 (RWorkspace: 0) */
4102
4103                 zlacpy_("U", m, n, &a[a_offset], lda, &vt[vt_offset], ldvt);
4104                 if (wntva) {
4105                     nrvt = *n;
4106                 }
4107                 if (wntvs) {
4108                     nrvt = *m;
4109                 }
4110                 i__2 = *lwork - iwork + 1;
4111                 zungbr_("P", &nrvt, n, m, &vt[vt_offset], ldvt, &work[itaup],
4112                         &work[iwork], &i__2, &ierr);
4113             }
4114             if (wntuo) {
4115
4116 /*              If left singular vectors desired in A, generate left
4117                 bidiagonalizing vectors in A
4118                 (CWorkspace: need 3*M-1, prefer 2*M+(M-1)*NB)
4119                 (RWorkspace: 0) */
4120
4121                 i__2 = *lwork - iwork + 1;
4122                 zungbr_("Q", m, m, n, &a[a_offset], lda, &work[itauq], &work[
4123                         iwork], &i__2, &ierr);
4124             }
4125             if (wntvo) {
4126
4127 /*              If right singular vectors desired in A, generate right
4128                 bidiagonalizing vectors in A
4129                 (CWorkspace: need 3*M, prefer 2*M+M*NB)
4130                 (RWorkspace: 0) */
4131
4132                 i__2 = *lwork - iwork + 1;
4133                 zungbr_("P", m, n, m, &a[a_offset], lda, &work[itaup], &work[
4134                         iwork], &i__2, &ierr);
4135             }
4136             irwork = ie + *m;
4137             if (wntuas || wntuo) {
4138                 nru = *m;
4139             }
4140             if (wntun) {
4141                 nru = 0;
4142             }
4143             if (wntvas || wntvo) {
4144                 ncvt = *n;
4145             }
4146             if (wntvn) {
4147                 ncvt = 0;
4148             }
4149             if (! wntuo && ! wntvo) {
4150
4151 /*              Perform bidiagonal QR iteration, if desired, computing
4152                 left singular vectors in U and computing right singular
4153                 vectors in VT
4154                 (CWorkspace: 0)
4155                 (RWorkspace: need BDSPAC) */
4156
4157                 zbdsqr_("L", m, &ncvt, &nru, &c__0, &s[1], &rwork[ie], &vt[
4158                         vt_offset], ldvt, &u[u_offset], ldu, cdum, &c__1, &
4159                         rwork[irwork], info);
4160             } else if (! wntuo && wntvo) {
4161
4162 /*              Perform bidiagonal QR iteration, if desired, computing
4163                 left singular vectors in U and computing right singular
4164                 vectors in A
4165                 (CWorkspace: 0)
4166                 (RWorkspace: need BDSPAC) */
4167
4168                 zbdsqr_("L", m, &ncvt, &nru, &c__0, &s[1], &rwork[ie], &a[
4169                         a_offset], lda, &u[u_offset], ldu, cdum, &c__1, &
4170                         rwork[irwork], info);
4171             } else {
4172
4173 /*              Perform bidiagonal QR iteration, if desired, computing
4174                 left singular vectors in A and computing right singular
4175                 vectors in VT
4176                 (CWorkspace: 0)
4177                 (RWorkspace: need BDSPAC) */
4178
4179                 zbdsqr_("L", m, &ncvt, &nru, &c__0, &s[1], &rwork[ie], &vt[
4180                         vt_offset], ldvt, &a[a_offset], lda, cdum, &c__1, &
4181                         rwork[irwork], info);
4182             }
4183
4184         }
4185
4186     }
4187
4188 /*     Undo scaling if necessary */
4189
4190     if (iscl == 1) {
4191         if (anrm > bignum) {
4192             dlascl_("G", &c__0, &c__0, &bignum, &anrm, &minmn, &c__1, &s[1], &
4193                     minmn, &ierr);
4194         }
4195         if (*info != 0 && anrm > bignum) {
4196             i__2 = minmn - 1;
4197             dlascl_("G", &c__0, &c__0, &bignum, &anrm, &i__2, &c__1, &rwork[
4198                     ie], &minmn, &ierr);
4199         }
4200         if (anrm < smlnum) {
4201             dlascl_("G", &c__0, &c__0, &smlnum, &anrm, &minmn, &c__1, &s[1], &
4202                     minmn, &ierr);
4203         }
4204         if (*info != 0 && anrm < smlnum) {
4205             i__2 = minmn - 1;
4206             dlascl_("G", &c__0, &c__0, &smlnum, &anrm, &i__2, &c__1, &rwork[
4207                     ie], &minmn, &ierr);
4208         }
4209     }
4210
4211 /*     Return optimal workspace in WORK(1) */
4212
4213     work[1].r = (doublereal) maxwrk, work[1].i = 0.;
4214
4215     return 0;
4216
4217 /*     End of ZGESVD */
4218
4219 } /* zgesvd_ */
4220
4221 #undef vt_ref
4222 #undef vt_subscr
4223 #undef u_ref
4224 #undef u_subscr
4225 #undef a_ref
4226 #undef a_subscr