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