git.maemo.org Git - opencv/blob - 3rdparty/lapack/dgesdd.c

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