swiftnav: sgesvd.c Source File

Go to the documentation of this file.
00001 /* sgesvd.f -- translated by f2c (version 20061008).
00002    You must link the resulting object file with libf2c:
00003         on Microsoft Windows system, link with libf2c.lib;
00004         on Linux or Unix systems, link with .../path/to/libf2c.a -lm
00005         or, if you install libf2c.a in a standard place, with -lf2c -lm
00006         -- in that order, at the end of the command line, as in
00007                 cc *.o -lf2c -lm
00008         Source for libf2c is in /netlib/f2c/libf2c.zip, e.g.,
00009 
00010                 http://www.netlib.org/f2c/libf2c.zip
00011 */
00012 
00013 #include "f2c.h"
00014 #include "blaswrap.h"
00015 
00016 /* Table of constant values */
00017 
00018 static integer c__6 = 6;
00019 static integer c__0 = 0;
00020 static integer c__2 = 2;
00021 static integer c__1 = 1;
00022 static integer c_n1 = -1;
00023 static real c_b421 = 0.f;
00024 static real c_b443 = 1.f;
00025 
00026 /* Subroutine */ int sgesvd_(char *jobu, char *jobvt, integer *m, integer *n, 
00027         real *a, integer *lda, real *s, real *u, integer *ldu, real *vt, 
00028         integer *ldvt, real *work, integer *lwork, integer *info)
00029 {
00030     /* System generated locals */
00031     address a__1[2];
00032     integer a_dim1, a_offset, u_dim1, u_offset, vt_dim1, vt_offset, i__1[2], 
00033             i__2, i__3, i__4;
00034     char ch__1[2];
00035 
00036     /* Builtin functions */
00037     /* Subroutine */ int s_cat(char *, char **, integer *, integer *, ftnlen);
00038     double sqrt(doublereal);
00039 
00040     /* Local variables */
00041     integer i__, ie, ir, iu, blk, ncu;
00042     real dum[1], eps;
00043     integer nru, iscl;
00044     real anrm;
00045     integer ierr, itau, ncvt, nrvt;
00046     extern logical lsame_(char *, char *);
00047     integer chunk;
00048     extern /* Subroutine */ int sgemm_(char *, char *, integer *, integer *, 
00049             integer *, real *, real *, integer *, real *, integer *, real *, 
00050             real *, integer *);
00051     integer minmn, wrkbl, itaup, itauq, mnthr, iwork;
00052     logical wntua, wntva, wntun, wntuo, wntvn, wntvo, wntus, wntvs;
00053     integer bdspac;
00054     extern /* Subroutine */ int sgebrd_(integer *, integer *, real *, integer 
00055             *, real *, real *, real *, real *, real *, integer *, integer *);
00056     extern doublereal slamch_(char *), slange_(char *, integer *, 
00057             integer *, real *, integer *, real *);
00058     extern /* Subroutine */ int xerbla_(char *, integer *);
00059     extern integer ilaenv_(integer *, char *, char *, integer *, integer *, 
00060             integer *, integer *);
00061     real bignum;
00062     extern /* Subroutine */ int sgelqf_(integer *, integer *, real *, integer 
00063             *, real *, real *, integer *, integer *), slascl_(char *, integer 
00064             *, integer *, real *, real *, integer *, integer *, real *, 
00065             integer *, integer *), sgeqrf_(integer *, integer *, real 
00066             *, integer *, real *, real *, integer *, integer *), slacpy_(char 
00067             *, integer *, integer *, real *, integer *, real *, integer *), slaset_(char *, integer *, integer *, real *, real *, 
00068             real *, integer *), sbdsqr_(char *, integer *, integer *, 
00069             integer *, integer *, real *, real *, real *, integer *, real *, 
00070             integer *, real *, integer *, real *, integer *), sorgbr_(
00071             char *, integer *, integer *, integer *, real *, integer *, real *
00072 , real *, integer *, integer *), sormbr_(char *, char *, 
00073             char *, integer *, integer *, integer *, real *, integer *, real *
00074 , real *, integer *, real *, integer *, integer *);
00075     integer ldwrkr, minwrk, ldwrku, maxwrk;
00076     extern /* Subroutine */ int sorglq_(integer *, integer *, integer *, real 
00077             *, integer *, real *, real *, integer *, integer *);
00078     real smlnum;
00079     extern /* Subroutine */ int sorgqr_(integer *, integer *, integer *, real 
00080             *, integer *, real *, real *, integer *, integer *);
00081     logical lquery, wntuas, wntvas;
00082 
00083 
00084 /*  -- LAPACK driver routine (version 3.2) -- */
00085 /*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */
00086 /*     November 2006 */
00087 
00088 /*     .. Scalar Arguments .. */
00089 /*     .. */
00090 /*     .. Array Arguments .. */
00091 /*     .. */
00092 
00093 /*  Purpose */
00094 /*  ======= */
00095 
00096 /*  SGESVD computes the singular value decomposition (SVD) of a real */
00097 /*  M-by-N matrix A, optionally computing the left and/or right singular */
00098 /*  vectors. The SVD is written */
00099 
00100 /*       A = U * SIGMA * transpose(V) */
00101 
00102 /*  where SIGMA is an M-by-N matrix which is zero except for its */
00103 /*  min(m,n) diagonal elements, U is an M-by-M orthogonal matrix, and */
00104 /*  V is an N-by-N orthogonal matrix.  The diagonal elements of SIGMA */
00105 /*  are the singular values of A; they are real and non-negative, and */
00106 /*  are returned in descending order.  The first min(m,n) columns of */
00107 /*  U and V are the left and right singular vectors of A. */
00108 
00109 /*  Note that the routine returns V**T, not V. */
00110 
00111 /*  Arguments */
00112 /*  ========= */
00113 
00114 /*  JOBU    (input) CHARACTER*1 */
00115 /*          Specifies options for computing all or part of the matrix U: */
00116 /*          = 'A':  all M columns of U are returned in array U: */
00117 /*          = 'S':  the first min(m,n) columns of U (the left singular */
00118 /*                  vectors) are returned in the array U; */
00119 /*          = 'O':  the first min(m,n) columns of U (the left singular */
00120 /*                  vectors) are overwritten on the array A; */
00121 /*          = 'N':  no columns of U (no left singular vectors) are */
00122 /*                  computed. */
00123 
00124 /*  JOBVT   (input) CHARACTER*1 */
00125 /*          Specifies options for computing all or part of the matrix */
00126 /*          V**T: */
00127 /*          = 'A':  all N rows of V**T are returned in the array VT; */
00128 /*          = 'S':  the first min(m,n) rows of V**T (the right singular */
00129 /*                  vectors) are returned in the array VT; */
00130 /*          = 'O':  the first min(m,n) rows of V**T (the right singular */
00131 /*                  vectors) are overwritten on the array A; */
00132 /*          = 'N':  no rows of V**T (no right singular vectors) are */
00133 /*                  computed. */
00134 
00135 /*          JOBVT and JOBU cannot both be 'O'. */
00136 
00137 /*  M       (input) INTEGER */
00138 /*          The number of rows of the input matrix A.  M >= 0. */
00139 
00140 /*  N       (input) INTEGER */
00141 /*          The number of columns of the input matrix A.  N >= 0. */
00142 
00143 /*  A       (input/output) REAL array, dimension (LDA,N) */
00144 /*          On entry, the M-by-N matrix A. */
00145 /*          On exit, */
00146 /*          if JOBU = 'O',  A is overwritten with the first min(m,n) */
00147 /*                          columns of U (the left singular vectors, */
00148 /*                          stored columnwise); */
00149 /*          if JOBVT = 'O', A is overwritten with the first min(m,n) */
00150 /*                          rows of V**T (the right singular vectors, */
00151 /*                          stored rowwise); */
00152 /*          if JOBU .ne. 'O' and JOBVT .ne. 'O', the contents of A */
00153 /*                          are destroyed. */
00154 
00155 /*  LDA     (input) INTEGER */
00156 /*          The leading dimension of the array A.  LDA >= max(1,M). */
00157 
00158 /*  S       (output) REAL array, dimension (min(M,N)) */
00159 /*          The singular values of A, sorted so that S(i) >= S(i+1). */
00160 
00161 /*  U       (output) REAL array, dimension (LDU,UCOL) */
00162 /*          (LDU,M) if JOBU = 'A' or (LDU,min(M,N)) if JOBU = 'S'. */
00163 /*          If JOBU = 'A', U contains the M-by-M orthogonal matrix U; */
00164 /*          if JOBU = 'S', U contains the first min(m,n) columns of U */
00165 /*          (the left singular vectors, stored columnwise); */
00166 /*          if JOBU = 'N' or 'O', U is not referenced. */
00167 
00168 /*  LDU     (input) INTEGER */
00169 /*          The leading dimension of the array U.  LDU >= 1; if */
00170 /*          JOBU = 'S' or 'A', LDU >= M. */
00171 
00172 /*  VT      (output) REAL array, dimension (LDVT,N) */
00173 /*          If JOBVT = 'A', VT contains the N-by-N orthogonal matrix */
00174 /*          V**T; */
00175 /*          if JOBVT = 'S', VT contains the first min(m,n) rows of */
00176 /*          V**T (the right singular vectors, stored rowwise); */
00177 /*          if JOBVT = 'N' or 'O', VT is not referenced. */
00178 
00179 /*  LDVT    (input) INTEGER */
00180 /*          The leading dimension of the array VT.  LDVT >= 1; if */
00181 /*          JOBVT = 'A', LDVT >= N; if JOBVT = 'S', LDVT >= min(M,N). */
00182 
00183 /*  WORK    (workspace/output) REAL array, dimension (MAX(1,LWORK)) */
00184 /*          On exit, if INFO = 0, WORK(1) returns the optimal LWORK; */
00185 /*          if INFO > 0, WORK(2:MIN(M,N)) contains the unconverged */
00186 /*          superdiagonal elements of an upper bidiagonal matrix B */
00187 /*          whose diagonal is in S (not necessarily sorted). B */
00188 /*          satisfies A = U * B * VT, so it has the same singular values */
00189 /*          as A, and singular vectors related by U and VT. */
00190 
00191 /*  LWORK   (input) INTEGER */
00192 /*          The dimension of the array WORK. */
00193 /*          LWORK >= MAX(1,3*MIN(M,N)+MAX(M,N),5*MIN(M,N)). */
00194 /*          For good performance, LWORK should generally be larger. */
00195 
00196 /*          If LWORK = -1, then a workspace query is assumed; the routine */
00197 /*          only calculates the optimal size of the WORK array, returns */
00198 /*          this value as the first entry of the WORK array, and no error */
00199 /*          message related to LWORK is issued by XERBLA. */
00200 
00201 /*  INFO    (output) INTEGER */
00202 /*          = 0:  successful exit. */
00203 /*          < 0:  if INFO = -i, the i-th argument had an illegal value. */
00204 /*          > 0:  if SBDSQR did not converge, INFO specifies how many */
00205 /*                superdiagonals of an intermediate bidiagonal form B */
00206 /*                did not converge to zero. See the description of WORK */
00207 /*                above for details. */
00208 
00209 /*  ===================================================================== */
00210 
00211 /*     .. Parameters .. */
00212 /*     .. */
00213 /*     .. Local Scalars .. */
00214 /*     .. */
00215 /*     .. Local Arrays .. */
00216 /*     .. */
00217 /*     .. External Subroutines .. */
00218 /*     .. */
00219 /*     .. External Functions .. */
00220 /*     .. */
00221 /*     .. Intrinsic Functions .. */
00222 /*     .. */
00223 /*     .. Executable Statements .. */
00224 
00225 /*     Test the input arguments */
00226 
00227     /* Parameter adjustments */
00228     a_dim1 = *lda;
00229     a_offset = 1 + a_dim1;
00230     a -= a_offset;
00231     --s;
00232     u_dim1 = *ldu;
00233     u_offset = 1 + u_dim1;
00234     u -= u_offset;
00235     vt_dim1 = *ldvt;
00236     vt_offset = 1 + vt_dim1;
00237     vt -= vt_offset;
00238     --work;
00239 
00240     /* Function Body */
00241     *info = 0;
00242     minmn = min(*m,*n);
00243     wntua = lsame_(jobu, "A");
00244     wntus = lsame_(jobu, "S");
00245     wntuas = wntua || wntus;
00246     wntuo = lsame_(jobu, "O");
00247     wntun = lsame_(jobu, "N");
00248     wntva = lsame_(jobvt, "A");
00249     wntvs = lsame_(jobvt, "S");
00250     wntvas = wntva || wntvs;
00251     wntvo = lsame_(jobvt, "O");
00252     wntvn = lsame_(jobvt, "N");
00253     lquery = *lwork == -1;
00254 
00255     if (! (wntua || wntus || wntuo || wntun)) {
00256         *info = -1;
00257     } else if (! (wntva || wntvs || wntvo || wntvn) || wntvo && wntuo) {
00258         *info = -2;
00259     } else if (*m < 0) {
00260         *info = -3;
00261     } else if (*n < 0) {
00262         *info = -4;
00263     } else if (*lda < max(1,*m)) {
00264         *info = -6;
00265     } else if (*ldu < 1 || wntuas && *ldu < *m) {
00266         *info = -9;
00267     } else if (*ldvt < 1 || wntva && *ldvt < *n || wntvs && *ldvt < minmn) {
00268         *info = -11;
00269     }
00270 
00271 /*     Compute workspace */
00272 /*      (Note: Comments in the code beginning "Workspace:" describe the */
00273 /*       minimal amount of workspace needed at that point in the code, */
00274 /*       as well as the preferred amount for good performance. */
00275 /*       NB refers to the optimal block size for the immediately */
00276 /*       following subroutine, as returned by ILAENV.) */
00277 
00278     if (*info == 0) {
00279         minwrk = 1;
00280         maxwrk = 1;
00281         if (*m >= *n && minmn > 0) {
00282 
00283 /*           Compute space needed for SBDSQR */
00284 
00285 /* Writing concatenation */
00286             i__1[0] = 1, a__1[0] = jobu;
00287             i__1[1] = 1, a__1[1] = jobvt;
00288             s_cat(ch__1, a__1, i__1, &c__2, (ftnlen)2);
00289             mnthr = ilaenv_(&c__6, "SGESVD", ch__1, m, n, &c__0, &c__0);
00290             bdspac = *n * 5;
00291             if (*m >= mnthr) {
00292                 if (wntun) {
00293 
00294 /*                 Path 1 (M much larger than N, JOBU='N') */
00295 
00296                     maxwrk = *n + *n * ilaenv_(&c__1, "SGEQRF", " ", m, n, &
00297                             c_n1, &c_n1);
00298 /* Computing MAX */
00299                     i__2 = maxwrk, i__3 = *n * 3 + (*n << 1) * ilaenv_(&c__1, 
00300                             "SGEBRD", " ", n, n, &c_n1, &c_n1);
00301                     maxwrk = max(i__2,i__3);
00302                     if (wntvo || wntvas) {
00303 /* Computing MAX */
00304                         i__2 = maxwrk, i__3 = *n * 3 + (*n - 1) * ilaenv_(&
00305                                 c__1, "SORGBR", "P", n, n, n, &c_n1);
00306                         maxwrk = max(i__2,i__3);
00307                     }
00308                     maxwrk = max(maxwrk,bdspac);
00309 /* Computing MAX */
00310                     i__2 = *n << 2;
00311                     minwrk = max(i__2,bdspac);
00312                 } else if (wntuo && wntvn) {
00313 
00314 /*                 Path 2 (M much larger than N, JOBU='O', JOBVT='N') */
00315 
00316                     wrkbl = *n + *n * ilaenv_(&c__1, "SGEQRF", " ", m, n, &
00317                             c_n1, &c_n1);
00318 /* Computing MAX */
00319                     i__2 = wrkbl, i__3 = *n + *n * ilaenv_(&c__1, "SORGQR", 
00320                             " ", m, n, n, &c_n1);
00321                     wrkbl = max(i__2,i__3);
00322 /* Computing MAX */
00323                     i__2 = wrkbl, i__3 = *n * 3 + (*n << 1) * ilaenv_(&c__1, 
00324                             "SGEBRD", " ", n, n, &c_n1, &c_n1);
00325                     wrkbl = max(i__2,i__3);
00326 /* Computing MAX */
00327                     i__2 = wrkbl, i__3 = *n * 3 + *n * ilaenv_(&c__1, "SORGBR"
00328 , "Q", n, n, n, &c_n1);
00329                     wrkbl = max(i__2,i__3);
00330                     wrkbl = max(wrkbl,bdspac);
00331 /* Computing MAX */
00332                     i__2 = *n * *n + wrkbl, i__3 = *n * *n + *m * *n + *n;
00333                     maxwrk = max(i__2,i__3);
00334 /* Computing MAX */
00335                     i__2 = *n * 3 + *m;
00336                     minwrk = max(i__2,bdspac);
00337                 } else if (wntuo && wntvas) {
00338 
00339 /*                 Path 3 (M much larger than N, JOBU='O', JOBVT='S' or */
00340 /*                 'A') */
00341 
00342                     wrkbl = *n + *n * ilaenv_(&c__1, "SGEQRF", " ", m, n, &
00343                             c_n1, &c_n1);
00344 /* Computing MAX */
00345                     i__2 = wrkbl, i__3 = *n + *n * ilaenv_(&c__1, "SORGQR", 
00346                             " ", m, n, n, &c_n1);
00347                     wrkbl = max(i__2,i__3);
00348 /* Computing MAX */
00349                     i__2 = wrkbl, i__3 = *n * 3 + (*n << 1) * ilaenv_(&c__1, 
00350                             "SGEBRD", " ", n, n, &c_n1, &c_n1);
00351                     wrkbl = max(i__2,i__3);
00352 /* Computing MAX */
00353                     i__2 = wrkbl, i__3 = *n * 3 + *n * ilaenv_(&c__1, "SORGBR"
00354 , "Q", n, n, n, &c_n1);
00355                     wrkbl = max(i__2,i__3);
00356 /* Computing MAX */
00357                     i__2 = wrkbl, i__3 = *n * 3 + (*n - 1) * ilaenv_(&c__1, 
00358                             "SORGBR", "P", n, n, n, &c_n1);
00359                     wrkbl = max(i__2,i__3);
00360                     wrkbl = max(wrkbl,bdspac);
00361 /* Computing MAX */
00362                     i__2 = *n * *n + wrkbl, i__3 = *n * *n + *m * *n + *n;
00363                     maxwrk = max(i__2,i__3);
00364 /* Computing MAX */
00365                     i__2 = *n * 3 + *m;
00366                     minwrk = max(i__2,bdspac);
00367                 } else if (wntus && wntvn) {
00368 
00369 /*                 Path 4 (M much larger than N, JOBU='S', JOBVT='N') */
00370 
00371                     wrkbl = *n + *n * ilaenv_(&c__1, "SGEQRF", " ", m, n, &
00372                             c_n1, &c_n1);
00373 /* Computing MAX */
00374                     i__2 = wrkbl, i__3 = *n + *n * ilaenv_(&c__1, "SORGQR", 
00375                             " ", m, n, n, &c_n1);
00376                     wrkbl = max(i__2,i__3);
00377 /* Computing MAX */
00378                     i__2 = wrkbl, i__3 = *n * 3 + (*n << 1) * ilaenv_(&c__1, 
00379                             "SGEBRD", " ", n, n, &c_n1, &c_n1);
00380                     wrkbl = max(i__2,i__3);
00381 /* Computing MAX */
00382                     i__2 = wrkbl, i__3 = *n * 3 + *n * ilaenv_(&c__1, "SORGBR"
00383 , "Q", n, n, n, &c_n1);
00384                     wrkbl = max(i__2,i__3);
00385                     wrkbl = max(wrkbl,bdspac);
00386                     maxwrk = *n * *n + wrkbl;
00387 /* Computing MAX */
00388                     i__2 = *n * 3 + *m;
00389                     minwrk = max(i__2,bdspac);
00390                 } else if (wntus && wntvo) {
00391 
00392 /*                 Path 5 (M much larger than N, JOBU='S', JOBVT='O') */
00393 
00394                     wrkbl = *n + *n * ilaenv_(&c__1, "SGEQRF", " ", m, n, &
00395                             c_n1, &c_n1);
00396 /* Computing MAX */
00397                     i__2 = wrkbl, i__3 = *n + *n * ilaenv_(&c__1, "SORGQR", 
00398                             " ", m, n, n, &c_n1);
00399                     wrkbl = max(i__2,i__3);
00400 /* Computing MAX */
00401                     i__2 = wrkbl, i__3 = *n * 3 + (*n << 1) * ilaenv_(&c__1, 
00402                             "SGEBRD", " ", n, n, &c_n1, &c_n1);
00403                     wrkbl = max(i__2,i__3);
00404 /* Computing MAX */
00405                     i__2 = wrkbl, i__3 = *n * 3 + *n * ilaenv_(&c__1, "SORGBR"
00406 , "Q", n, n, n, &c_n1);
00407                     wrkbl = max(i__2,i__3);
00408 /* Computing MAX */
00409                     i__2 = wrkbl, i__3 = *n * 3 + (*n - 1) * ilaenv_(&c__1, 
00410                             "SORGBR", "P", n, n, n, &c_n1);
00411                     wrkbl = max(i__2,i__3);
00412                     wrkbl = max(wrkbl,bdspac);
00413                     maxwrk = (*n << 1) * *n + wrkbl;
00414 /* Computing MAX */
00415                     i__2 = *n * 3 + *m;
00416                     minwrk = max(i__2,bdspac);
00417                 } else if (wntus && wntvas) {
00418 
00419 /*                 Path 6 (M much larger than N, JOBU='S', JOBVT='S' or */
00420 /*                 'A') */
00421 
00422                     wrkbl = *n + *n * ilaenv_(&c__1, "SGEQRF", " ", m, n, &
00423                             c_n1, &c_n1);
00424 /* Computing MAX */
00425                     i__2 = wrkbl, i__3 = *n + *n * ilaenv_(&c__1, "SORGQR", 
00426                             " ", m, n, n, &c_n1);
00427                     wrkbl = max(i__2,i__3);
00428 /* Computing MAX */
00429                     i__2 = wrkbl, i__3 = *n * 3 + (*n << 1) * ilaenv_(&c__1, 
00430                             "SGEBRD", " ", n, n, &c_n1, &c_n1);
00431                     wrkbl = max(i__2,i__3);
00432 /* Computing MAX */
00433                     i__2 = wrkbl, i__3 = *n * 3 + *n * ilaenv_(&c__1, "SORGBR"
00434 , "Q", n, n, n, &c_n1);
00435                     wrkbl = max(i__2,i__3);
00436 /* Computing MAX */
00437                     i__2 = wrkbl, i__3 = *n * 3 + (*n - 1) * ilaenv_(&c__1, 
00438                             "SORGBR", "P", n, n, n, &c_n1);
00439                     wrkbl = max(i__2,i__3);
00440                     wrkbl = max(wrkbl,bdspac);
00441                     maxwrk = *n * *n + wrkbl;
00442 /* Computing MAX */
00443                     i__2 = *n * 3 + *m;
00444                     minwrk = max(i__2,bdspac);
00445                 } else if (wntua && wntvn) {
00446 
00447 /*                 Path 7 (M much larger than N, JOBU='A', JOBVT='N') */
00448 
00449                     wrkbl = *n + *n * ilaenv_(&c__1, "SGEQRF", " ", m, n, &
00450                             c_n1, &c_n1);
00451 /* Computing MAX */
00452                     i__2 = wrkbl, i__3 = *n + *m * ilaenv_(&c__1, "SORGQR", 
00453                             " ", m, m, n, &c_n1);
00454                     wrkbl = max(i__2,i__3);
00455 /* Computing MAX */
00456                     i__2 = wrkbl, i__3 = *n * 3 + (*n << 1) * ilaenv_(&c__1, 
00457                             "SGEBRD", " ", n, n, &c_n1, &c_n1);
00458                     wrkbl = max(i__2,i__3);
00459 /* Computing MAX */
00460                     i__2 = wrkbl, i__3 = *n * 3 + *n * ilaenv_(&c__1, "SORGBR"
00461 , "Q", n, n, n, &c_n1);
00462                     wrkbl = max(i__2,i__3);
00463                     wrkbl = max(wrkbl,bdspac);
00464                     maxwrk = *n * *n + wrkbl;
00465 /* Computing MAX */
00466                     i__2 = *n * 3 + *m;
00467                     minwrk = max(i__2,bdspac);
00468                 } else if (wntua && wntvo) {
00469 
00470 /*                 Path 8 (M much larger than N, JOBU='A', JOBVT='O') */
00471 
00472                     wrkbl = *n + *n * ilaenv_(&c__1, "SGEQRF", " ", m, n, &
00473                             c_n1, &c_n1);
00474 /* Computing MAX */
00475                     i__2 = wrkbl, i__3 = *n + *m * ilaenv_(&c__1, "SORGQR", 
00476                             " ", m, m, n, &c_n1);
00477                     wrkbl = max(i__2,i__3);
00478 /* Computing MAX */
00479                     i__2 = wrkbl, i__3 = *n * 3 + (*n << 1) * ilaenv_(&c__1, 
00480                             "SGEBRD", " ", n, n, &c_n1, &c_n1);
00481                     wrkbl = max(i__2,i__3);
00482 /* Computing MAX */
00483                     i__2 = wrkbl, i__3 = *n * 3 + *n * ilaenv_(&c__1, "SORGBR"
00484 , "Q", n, n, n, &c_n1);
00485                     wrkbl = max(i__2,i__3);
00486 /* Computing MAX */
00487                     i__2 = wrkbl, i__3 = *n * 3 + (*n - 1) * ilaenv_(&c__1, 
00488                             "SORGBR", "P", n, n, n, &c_n1);
00489                     wrkbl = max(i__2,i__3);
00490                     wrkbl = max(wrkbl,bdspac);
00491                     maxwrk = (*n << 1) * *n + wrkbl;
00492 /* Computing MAX */
00493                     i__2 = *n * 3 + *m;
00494                     minwrk = max(i__2,bdspac);
00495                 } else if (wntua && wntvas) {
00496 
00497 /*                 Path 9 (M much larger than N, JOBU='A', JOBVT='S' or */
00498 /*                 'A') */
00499 
00500                     wrkbl = *n + *n * ilaenv_(&c__1, "SGEQRF", " ", m, n, &
00501                             c_n1, &c_n1);
00502 /* Computing MAX */
00503                     i__2 = wrkbl, i__3 = *n + *m * ilaenv_(&c__1, "SORGQR", 
00504                             " ", m, m, n, &c_n1);
00505                     wrkbl = max(i__2,i__3);
00506 /* Computing MAX */
00507                     i__2 = wrkbl, i__3 = *n * 3 + (*n << 1) * ilaenv_(&c__1, 
00508                             "SGEBRD", " ", n, n, &c_n1, &c_n1);
00509                     wrkbl = max(i__2,i__3);
00510 /* Computing MAX */
00511                     i__2 = wrkbl, i__3 = *n * 3 + *n * ilaenv_(&c__1, "SORGBR"
00512 , "Q", n, n, n, &c_n1);
00513                     wrkbl = max(i__2,i__3);
00514 /* Computing MAX */
00515                     i__2 = wrkbl, i__3 = *n * 3 + (*n - 1) * ilaenv_(&c__1, 
00516                             "SORGBR", "P", n, n, n, &c_n1);
00517                     wrkbl = max(i__2,i__3);
00518                     wrkbl = max(wrkbl,bdspac);
00519                     maxwrk = *n * *n + wrkbl;
00520 /* Computing MAX */
00521                     i__2 = *n * 3 + *m;
00522                     minwrk = max(i__2,bdspac);
00523                 }
00524             } else {
00525 
00526 /*              Path 10 (M at least N, but not much larger) */
00527 
00528                 maxwrk = *n * 3 + (*m + *n) * ilaenv_(&c__1, "SGEBRD", " ", m, 
00529                          n, &c_n1, &c_n1);
00530                 if (wntus || wntuo) {
00531 /* Computing MAX */
00532                     i__2 = maxwrk, i__3 = *n * 3 + *n * ilaenv_(&c__1, "SORG"
00533                             "BR", "Q", m, n, n, &c_n1);
00534                     maxwrk = max(i__2,i__3);
00535                 }
00536                 if (wntua) {
00537 /* Computing MAX */
00538                     i__2 = maxwrk, i__3 = *n * 3 + *m * ilaenv_(&c__1, "SORG"
00539                             "BR", "Q", m, m, n, &c_n1);
00540                     maxwrk = max(i__2,i__3);
00541                 }
00542                 if (! wntvn) {
00543 /* Computing MAX */
00544                     i__2 = maxwrk, i__3 = *n * 3 + (*n - 1) * ilaenv_(&c__1, 
00545                             "SORGBR", "P", n, n, n, &c_n1);
00546                     maxwrk = max(i__2,i__3);
00547                 }
00548                 maxwrk = max(maxwrk,bdspac);
00549 /* Computing MAX */
00550                 i__2 = *n * 3 + *m;
00551                 minwrk = max(i__2,bdspac);
00552             }
00553         } else if (minmn > 0) {
00554 
00555 /*           Compute space needed for SBDSQR */
00556 
00557 /* Writing concatenation */
00558             i__1[0] = 1, a__1[0] = jobu;
00559             i__1[1] = 1, a__1[1] = jobvt;
00560             s_cat(ch__1, a__1, i__1, &c__2, (ftnlen)2);
00561             mnthr = ilaenv_(&c__6, "SGESVD", ch__1, m, n, &c__0, &c__0);
00562             bdspac = *m * 5;
00563             if (*n >= mnthr) {
00564                 if (wntvn) {
00565 
00566 /*                 Path 1t(N much larger than M, JOBVT='N') */
00567 
00568                     maxwrk = *m + *m * ilaenv_(&c__1, "SGELQF", " ", m, n, &
00569                             c_n1, &c_n1);
00570 /* Computing MAX */
00571                     i__2 = maxwrk, i__3 = *m * 3 + (*m << 1) * ilaenv_(&c__1, 
00572                             "SGEBRD", " ", m, m, &c_n1, &c_n1);
00573                     maxwrk = max(i__2,i__3);
00574                     if (wntuo || wntuas) {
00575 /* Computing MAX */
00576                         i__2 = maxwrk, i__3 = *m * 3 + *m * ilaenv_(&c__1, 
00577                                 "SORGBR", "Q", m, m, m, &c_n1);
00578                         maxwrk = max(i__2,i__3);
00579                     }
00580                     maxwrk = max(maxwrk,bdspac);
00581 /* Computing MAX */
00582                     i__2 = *m << 2;
00583                     minwrk = max(i__2,bdspac);
00584                 } else if (wntvo && wntun) {
00585 
00586 /*                 Path 2t(N much larger than M, JOBU='N', JOBVT='O') */
00587 
00588                     wrkbl = *m + *m * ilaenv_(&c__1, "SGELQF", " ", m, n, &
00589                             c_n1, &c_n1);
00590 /* Computing MAX */
00591                     i__2 = wrkbl, i__3 = *m + *m * ilaenv_(&c__1, "SORGLQ", 
00592                             " ", m, n, m, &c_n1);
00593                     wrkbl = max(i__2,i__3);
00594 /* Computing MAX */
00595                     i__2 = wrkbl, i__3 = *m * 3 + (*m << 1) * ilaenv_(&c__1, 
00596                             "SGEBRD", " ", m, m, &c_n1, &c_n1);
00597                     wrkbl = max(i__2,i__3);
00598 /* Computing MAX */
00599                     i__2 = wrkbl, i__3 = *m * 3 + (*m - 1) * ilaenv_(&c__1, 
00600                             "SORGBR", "P", m, m, m, &c_n1);
00601                     wrkbl = max(i__2,i__3);
00602                     wrkbl = max(wrkbl,bdspac);
00603 /* Computing MAX */
00604                     i__2 = *m * *m + wrkbl, i__3 = *m * *m + *m * *n + *m;
00605                     maxwrk = max(i__2,i__3);
00606 /* Computing MAX */
00607                     i__2 = *m * 3 + *n;
00608                     minwrk = max(i__2,bdspac);
00609                 } else if (wntvo && wntuas) {
00610 
00611 /*                 Path 3t(N much larger than M, JOBU='S' or 'A', */
00612 /*                 JOBVT='O') */
00613 
00614                     wrkbl = *m + *m * ilaenv_(&c__1, "SGELQF", " ", m, n, &
00615                             c_n1, &c_n1);
00616 /* Computing MAX */
00617                     i__2 = wrkbl, i__3 = *m + *m * ilaenv_(&c__1, "SORGLQ", 
00618                             " ", m, n, m, &c_n1);
00619                     wrkbl = max(i__2,i__3);
00620 /* Computing MAX */
00621                     i__2 = wrkbl, i__3 = *m * 3 + (*m << 1) * ilaenv_(&c__1, 
00622                             "SGEBRD", " ", m, m, &c_n1, &c_n1);
00623                     wrkbl = max(i__2,i__3);
00624 /* Computing MAX */
00625                     i__2 = wrkbl, i__3 = *m * 3 + (*m - 1) * ilaenv_(&c__1, 
00626                             "SORGBR", "P", m, m, m, &c_n1);
00627                     wrkbl = max(i__2,i__3);
00628 /* Computing MAX */
00629                     i__2 = wrkbl, i__3 = *m * 3 + *m * ilaenv_(&c__1, "SORGBR"
00630 , "Q", m, m, m, &c_n1);
00631                     wrkbl = max(i__2,i__3);
00632                     wrkbl = max(wrkbl,bdspac);
00633 /* Computing MAX */
00634                     i__2 = *m * *m + wrkbl, i__3 = *m * *m + *m * *n + *m;
00635                     maxwrk = max(i__2,i__3);
00636 /* Computing MAX */
00637                     i__2 = *m * 3 + *n;
00638                     minwrk = max(i__2,bdspac);
00639                 } else if (wntvs && wntun) {
00640 
00641 /*                 Path 4t(N much larger than M, JOBU='N', JOBVT='S') */
00642 
00643                     wrkbl = *m + *m * ilaenv_(&c__1, "SGELQF", " ", m, n, &
00644                             c_n1, &c_n1);
00645 /* Computing MAX */
00646                     i__2 = wrkbl, i__3 = *m + *m * ilaenv_(&c__1, "SORGLQ", 
00647                             " ", m, n, m, &c_n1);
00648                     wrkbl = max(i__2,i__3);
00649 /* Computing MAX */
00650                     i__2 = wrkbl, i__3 = *m * 3 + (*m << 1) * ilaenv_(&c__1, 
00651                             "SGEBRD", " ", m, m, &c_n1, &c_n1);
00652                     wrkbl = max(i__2,i__3);
00653 /* Computing MAX */
00654                     i__2 = wrkbl, i__3 = *m * 3 + (*m - 1) * ilaenv_(&c__1, 
00655                             "SORGBR", "P", m, m, m, &c_n1);
00656                     wrkbl = max(i__2,i__3);
00657                     wrkbl = max(wrkbl,bdspac);
00658                     maxwrk = *m * *m + wrkbl;
00659 /* Computing MAX */
00660                     i__2 = *m * 3 + *n;
00661                     minwrk = max(i__2,bdspac);
00662                 } else if (wntvs && wntuo) {
00663 
00664 /*                 Path 5t(N much larger than M, JOBU='O', JOBVT='S') */
00665 
00666                     wrkbl = *m + *m * ilaenv_(&c__1, "SGELQF", " ", m, n, &
00667                             c_n1, &c_n1);
00668 /* Computing MAX */
00669                     i__2 = wrkbl, i__3 = *m + *m * ilaenv_(&c__1, "SORGLQ", 
00670                             " ", m, n, m, &c_n1);
00671                     wrkbl = max(i__2,i__3);
00672 /* Computing MAX */
00673                     i__2 = wrkbl, i__3 = *m * 3 + (*m << 1) * ilaenv_(&c__1, 
00674                             "SGEBRD", " ", m, m, &c_n1, &c_n1);
00675                     wrkbl = max(i__2,i__3);
00676 /* Computing MAX */
00677                     i__2 = wrkbl, i__3 = *m * 3 + (*m - 1) * ilaenv_(&c__1, 
00678                             "SORGBR", "P", m, m, m, &c_n1);
00679                     wrkbl = max(i__2,i__3);
00680 /* Computing MAX */
00681                     i__2 = wrkbl, i__3 = *m * 3 + *m * ilaenv_(&c__1, "SORGBR"
00682 , "Q", m, m, m, &c_n1);
00683                     wrkbl = max(i__2,i__3);
00684                     wrkbl = max(wrkbl,bdspac);
00685                     maxwrk = (*m << 1) * *m + wrkbl;
00686 /* Computing MAX */
00687                     i__2 = *m * 3 + *n;
00688                     minwrk = max(i__2,bdspac);
00689                     maxwrk = max(maxwrk,minwrk);
00690                 } else if (wntvs && wntuas) {
00691 
00692 /*                 Path 6t(N much larger than M, JOBU='S' or 'A', */
00693 /*                 JOBVT='S') */
00694 
00695                     wrkbl = *m + *m * ilaenv_(&c__1, "SGELQF", " ", m, n, &
00696                             c_n1, &c_n1);
00697 /* Computing MAX */
00698                     i__2 = wrkbl, i__3 = *m + *m * ilaenv_(&c__1, "SORGLQ", 
00699                             " ", m, n, m, &c_n1);
00700                     wrkbl = max(i__2,i__3);
00701 /* Computing MAX */
00702                     i__2 = wrkbl, i__3 = *m * 3 + (*m << 1) * ilaenv_(&c__1, 
00703                             "SGEBRD", " ", m, m, &c_n1, &c_n1);
00704                     wrkbl = max(i__2,i__3);
00705 /* Computing MAX */
00706                     i__2 = wrkbl, i__3 = *m * 3 + (*m - 1) * ilaenv_(&c__1, 
00707                             "SORGBR", "P", m, m, m, &c_n1);
00708                     wrkbl = max(i__2,i__3);
00709 /* Computing MAX */
00710                     i__2 = wrkbl, i__3 = *m * 3 + *m * ilaenv_(&c__1, "SORGBR"
00711 , "Q", m, m, m, &c_n1);
00712                     wrkbl = max(i__2,i__3);
00713                     wrkbl = max(wrkbl,bdspac);
00714                     maxwrk = *m * *m + wrkbl;
00715 /* Computing MAX */
00716                     i__2 = *m * 3 + *n;
00717                     minwrk = max(i__2,bdspac);
00718                 } else if (wntva && wntun) {
00719 
00720 /*                 Path 7t(N much larger than M, JOBU='N', JOBVT='A') */
00721 
00722                     wrkbl = *m + *m * ilaenv_(&c__1, "SGELQF", " ", m, n, &
00723                             c_n1, &c_n1);
00724 /* Computing MAX */
00725                     i__2 = wrkbl, i__3 = *m + *n * ilaenv_(&c__1, "SORGLQ", 
00726                             " ", n, n, m, &c_n1);
00727                     wrkbl = max(i__2,i__3);
00728 /* Computing MAX */
00729                     i__2 = wrkbl, i__3 = *m * 3 + (*m << 1) * ilaenv_(&c__1, 
00730                             "SGEBRD", " ", m, m, &c_n1, &c_n1);
00731                     wrkbl = max(i__2,i__3);
00732 /* Computing MAX */
00733                     i__2 = wrkbl, i__3 = *m * 3 + (*m - 1) * ilaenv_(&c__1, 
00734                             "SORGBR", "P", m, m, m, &c_n1);
00735                     wrkbl = max(i__2,i__3);
00736                     wrkbl = max(wrkbl,bdspac);
00737                     maxwrk = *m * *m + wrkbl;
00738 /* Computing MAX */
00739                     i__2 = *m * 3 + *n;
00740                     minwrk = max(i__2,bdspac);
00741                 } else if (wntva && wntuo) {
00742 
00743 /*                 Path 8t(N much larger than M, JOBU='O', JOBVT='A') */
00744 
00745                     wrkbl = *m + *m * ilaenv_(&c__1, "SGELQF", " ", m, n, &
00746                             c_n1, &c_n1);
00747 /* Computing MAX */
00748                     i__2 = wrkbl, i__3 = *m + *n * ilaenv_(&c__1, "SORGLQ", 
00749                             " ", n, n, m, &c_n1);
00750                     wrkbl = max(i__2,i__3);
00751 /* Computing MAX */
00752                     i__2 = wrkbl, i__3 = *m * 3 + (*m << 1) * ilaenv_(&c__1, 
00753                             "SGEBRD", " ", m, m, &c_n1, &c_n1);
00754                     wrkbl = max(i__2,i__3);
00755 /* Computing MAX */
00756                     i__2 = wrkbl, i__3 = *m * 3 + (*m - 1) * ilaenv_(&c__1, 
00757                             "SORGBR", "P", m, m, m, &c_n1);
00758                     wrkbl = max(i__2,i__3);
00759 /* Computing MAX */
00760                     i__2 = wrkbl, i__3 = *m * 3 + *m * ilaenv_(&c__1, "SORGBR"
00761 , "Q", m, m, m, &c_n1);
00762                     wrkbl = max(i__2,i__3);
00763                     wrkbl = max(wrkbl,bdspac);
00764                     maxwrk = (*m << 1) * *m + wrkbl;
00765 /* Computing MAX */
00766                     i__2 = *m * 3 + *n;
00767                     minwrk = max(i__2,bdspac);
00768                 } else if (wntva && wntuas) {
00769 
00770 /*                 Path 9t(N much larger than M, JOBU='S' or 'A', */
00771 /*                 JOBVT='A') */
00772 
00773                     wrkbl = *m + *m * ilaenv_(&c__1, "SGELQF", " ", m, n, &
00774                             c_n1, &c_n1);
00775 /* Computing MAX */
00776                     i__2 = wrkbl, i__3 = *m + *n * ilaenv_(&c__1, "SORGLQ", 
00777                             " ", n, n, m, &c_n1);
00778                     wrkbl = max(i__2,i__3);
00779 /* Computing MAX */
00780                     i__2 = wrkbl, i__3 = *m * 3 + (*m << 1) * ilaenv_(&c__1, 
00781                             "SGEBRD", " ", m, m, &c_n1, &c_n1);
00782                     wrkbl = max(i__2,i__3);
00783 /* Computing MAX */
00784                     i__2 = wrkbl, i__3 = *m * 3 + (*m - 1) * ilaenv_(&c__1, 
00785                             "SORGBR", "P", m, m, m, &c_n1);
00786                     wrkbl = max(i__2,i__3);
00787 /* Computing MAX */
00788                     i__2 = wrkbl, i__3 = *m * 3 + *m * ilaenv_(&c__1, "SORGBR"
00789 , "Q", m, m, m, &c_n1);
00790                     wrkbl = max(i__2,i__3);
00791                     wrkbl = max(wrkbl,bdspac);
00792                     maxwrk = *m * *m + wrkbl;
00793 /* Computing MAX */
00794                     i__2 = *m * 3 + *n;
00795                     minwrk = max(i__2,bdspac);
00796                 }
00797             } else {
00798 
00799 /*              Path 10t(N greater than M, but not much larger) */
00800 
00801                 maxwrk = *m * 3 + (*m + *n) * ilaenv_(&c__1, "SGEBRD", " ", m, 
00802                          n, &c_n1, &c_n1);
00803                 if (wntvs || wntvo) {
00804 /* Computing MAX */
00805                     i__2 = maxwrk, i__3 = *m * 3 + *m * ilaenv_(&c__1, "SORG"
00806                             "BR", "P", m, n, m, &c_n1);
00807                     maxwrk = max(i__2,i__3);
00808                 }
00809                 if (wntva) {
00810 /* Computing MAX */
00811                     i__2 = maxwrk, i__3 = *m * 3 + *n * ilaenv_(&c__1, "SORG"
00812                             "BR", "P", n, n, m, &c_n1);
00813                     maxwrk = max(i__2,i__3);
00814                 }
00815                 if (! wntun) {
00816 /* Computing MAX */
00817                     i__2 = maxwrk, i__3 = *m * 3 + (*m - 1) * ilaenv_(&c__1, 
00818                             "SORGBR", "Q", m, m, m, &c_n1);
00819                     maxwrk = max(i__2,i__3);
00820                 }
00821                 maxwrk = max(maxwrk,bdspac);
00822 /* Computing MAX */
00823                 i__2 = *m * 3 + *n;
00824                 minwrk = max(i__2,bdspac);
00825             }
00826         }
00827         maxwrk = max(maxwrk,minwrk);
00828         work[1] = (real) maxwrk;
00829 
00830         if (*lwork < minwrk && ! lquery) {
00831             *info = -13;
00832         }
00833     }
00834 
00835     if (*info != 0) {
00836         i__2 = -(*info);
00837         xerbla_("SGESVD", &i__2);
00838         return 0;
00839     } else if (lquery) {
00840         return 0;
00841     }
00842 
00843 /*     Quick return if possible */
00844 
00845     if (*m == 0 || *n == 0) {
00846         return 0;
00847     }
00848 
00849 /*     Get machine constants */
00850 
00851     eps = slamch_("P");
00852     smlnum = sqrt(slamch_("S")) / eps;
00853     bignum = 1.f / smlnum;
00854 
00855 /*     Scale A if max element outside range [SMLNUM,BIGNUM] */
00856 
00857     anrm = slange_("M", m, n, &a[a_offset], lda, dum);
00858     iscl = 0;
00859     if (anrm > 0.f && anrm < smlnum) {
00860         iscl = 1;
00861         slascl_("G", &c__0, &c__0, &anrm, &smlnum, m, n, &a[a_offset], lda, &
00862                 ierr);
00863     } else if (anrm > bignum) {
00864         iscl = 1;
00865         slascl_("G", &c__0, &c__0, &anrm, &bignum, m, n, &a[a_offset], lda, &
00866                 ierr);
00867     }
00868 
00869     if (*m >= *n) {
00870 
00871 /*        A has at least as many rows as columns. If A has sufficiently */
00872 /*        more rows than columns, first reduce using the QR */
00873 /*        decomposition (if sufficient workspace available) */
00874 
00875         if (*m >= mnthr) {
00876 
00877             if (wntun) {
00878 
00879 /*              Path 1 (M much larger than N, JOBU='N') */
00880 /*              No left singular vectors to be computed */
00881 
00882                 itau = 1;
00883                 iwork = itau + *n;
00884 
00885 /*              Compute A=Q*R */
00886 /*              (Workspace: need 2*N, prefer N+N*NB) */
00887 
00888                 i__2 = *lwork - iwork + 1;
00889                 sgeqrf_(m, n, &a[a_offset], lda, &work[itau], &work[iwork], &
00890                         i__2, &ierr);
00891 
00892 /*              Zero out below R */
00893 
00894                 i__2 = *n - 1;
00895                 i__3 = *n - 1;
00896                 slaset_("L", &i__2, &i__3, &c_b421, &c_b421, &a[a_dim1 + 2], 
00897                         lda);
00898                 ie = 1;
00899                 itauq = ie + *n;
00900                 itaup = itauq + *n;
00901                 iwork = itaup + *n;
00902 
00903 /*              Bidiagonalize R in A */
00904 /*              (Workspace: need 4*N, prefer 3*N+2*N*NB) */
00905 
00906                 i__2 = *lwork - iwork + 1;
00907                 sgebrd_(n, n, &a[a_offset], lda, &s[1], &work[ie], &work[
00908                         itauq], &work[itaup], &work[iwork], &i__2, &ierr);
00909                 ncvt = 0;
00910                 if (wntvo || wntvas) {
00911 
00912 /*                 If right singular vectors desired, generate P'. */
00913 /*                 (Workspace: need 4*N-1, prefer 3*N+(N-1)*NB) */
00914 
00915                     i__2 = *lwork - iwork + 1;
00916                     sorgbr_("P", n, n, n, &a[a_offset], lda, &work[itaup], &
00917                             work[iwork], &i__2, &ierr);
00918                     ncvt = *n;
00919                 }
00920                 iwork = ie + *n;
00921 
00922 /*              Perform bidiagonal QR iteration, computing right */
00923 /*              singular vectors of A in A if desired */
00924 /*              (Workspace: need BDSPAC) */
00925 
00926                 sbdsqr_("U", n, &ncvt, &c__0, &c__0, &s[1], &work[ie], &a[
00927                         a_offset], lda, dum, &c__1, dum, &c__1, &work[iwork], 
00928                         info);
00929 
00930 /*              If right singular vectors desired in VT, copy them there */
00931 
00932                 if (wntvas) {
00933                     slacpy_("F", n, n, &a[a_offset], lda, &vt[vt_offset], 
00934                             ldvt);
00935                 }
00936 
00937             } else if (wntuo && wntvn) {
00938 
00939 /*              Path 2 (M much larger than N, JOBU='O', JOBVT='N') */
00940 /*              N left singular vectors to be overwritten on A and */
00941 /*              no right singular vectors to be computed */
00942 
00943 /* Computing MAX */
00944                 i__2 = *n << 2;
00945                 if (*lwork >= *n * *n + max(i__2,bdspac)) {
00946 
00947 /*                 Sufficient workspace for a fast algorithm */
00948 
00949                     ir = 1;
00950 /* Computing MAX */
00951                     i__2 = wrkbl, i__3 = *lda * *n + *n;
00952                     if (*lwork >= max(i__2,i__3) + *lda * *n) {
00953 
00954 /*                    WORK(IU) is LDA by N, WORK(IR) is LDA by N */
00955 
00956                         ldwrku = *lda;
00957                         ldwrkr = *lda;
00958                     } else /* if(complicated condition) */ {
00959 /* Computing MAX */
00960                         i__2 = wrkbl, i__3 = *lda * *n + *n;
00961                         if (*lwork >= max(i__2,i__3) + *n * *n) {
00962 
00963 /*                    WORK(IU) is LDA by N, WORK(IR) is N by N */
00964 
00965                             ldwrku = *lda;
00966                             ldwrkr = *n;
00967                         } else {
00968 
00969 /*                    WORK(IU) is LDWRKU by N, WORK(IR) is N by N */
00970 
00971                             ldwrku = (*lwork - *n * *n - *n) / *n;
00972                             ldwrkr = *n;
00973                         }
00974                     }
00975                     itau = ir + ldwrkr * *n;
00976                     iwork = itau + *n;
00977 
00978 /*                 Compute A=Q*R */
00979 /*                 (Workspace: need N*N+2*N, prefer N*N+N+N*NB) */
00980 
00981                     i__2 = *lwork - iwork + 1;
00982                     sgeqrf_(m, n, &a[a_offset], lda, &work[itau], &work[iwork]
00983 , &i__2, &ierr);
00984 
00985 /*                 Copy R to WORK(IR) and zero out below it */
00986 
00987                     slacpy_("U", n, n, &a[a_offset], lda, &work[ir], &ldwrkr);
00988                     i__2 = *n - 1;
00989                     i__3 = *n - 1;
00990                     slaset_("L", &i__2, &i__3, &c_b421, &c_b421, &work[ir + 1]
00991 , &ldwrkr);
00992 
00993 /*                 Generate Q in A */
00994 /*                 (Workspace: need N*N+2*N, prefer N*N+N+N*NB) */
00995 
00996                     i__2 = *lwork - iwork + 1;
00997                     sorgqr_(m, n, n, &a[a_offset], lda, &work[itau], &work[
00998                             iwork], &i__2, &ierr);
00999                     ie = itau;
01000                     itauq = ie + *n;
01001                     itaup = itauq + *n;
01002                     iwork = itaup + *n;
01003 
01004 /*                 Bidiagonalize R in WORK(IR) */
01005 /*                 (Workspace: need N*N+4*N, prefer N*N+3*N+2*N*NB) */
01006 
01007                     i__2 = *lwork - iwork + 1;
01008                     sgebrd_(n, n, &work[ir], &ldwrkr, &s[1], &work[ie], &work[
01009                             itauq], &work[itaup], &work[iwork], &i__2, &ierr);
01010 
01011 /*                 Generate left vectors bidiagonalizing R */
01012 /*                 (Workspace: need N*N+4*N, prefer N*N+3*N+N*NB) */
01013 
01014                     i__2 = *lwork - iwork + 1;
01015                     sorgbr_("Q", n, n, n, &work[ir], &ldwrkr, &work[itauq], &
01016                             work[iwork], &i__2, &ierr);
01017                     iwork = ie + *n;
01018 
01019 /*                 Perform bidiagonal QR iteration, computing left */
01020 /*                 singular vectors of R in WORK(IR) */
01021 /*                 (Workspace: need N*N+BDSPAC) */
01022 
01023                     sbdsqr_("U", n, &c__0, n, &c__0, &s[1], &work[ie], dum, &
01024                             c__1, &work[ir], &ldwrkr, dum, &c__1, &work[iwork]
01025 , info);
01026                     iu = ie + *n;
01027 
01028 /*                 Multiply Q in A by left singular vectors of R in */
01029 /*                 WORK(IR), storing result in WORK(IU) and copying to A */
01030 /*                 (Workspace: need N*N+2*N, prefer N*N+M*N+N) */
01031 
01032                     i__2 = *m;
01033                     i__3 = ldwrku;
01034                     for (i__ = 1; i__3 < 0 ? i__ >= i__2 : i__ <= i__2; i__ +=
01035                              i__3) {
01036 /* Computing MIN */
01037                         i__4 = *m - i__ + 1;
01038                         chunk = min(i__4,ldwrku);
01039                         sgemm_("N", "N", &chunk, n, n, &c_b443, &a[i__ + 
01040                                 a_dim1], lda, &work[ir], &ldwrkr, &c_b421, &
01041                                 work[iu], &ldwrku);
01042                         slacpy_("F", &chunk, n, &work[iu], &ldwrku, &a[i__ + 
01043                                 a_dim1], lda);
01044 /* L10: */
01045                     }
01046 
01047                 } else {
01048 
01049 /*                 Insufficient workspace for a fast algorithm */
01050 
01051                     ie = 1;
01052                     itauq = ie + *n;
01053                     itaup = itauq + *n;
01054                     iwork = itaup + *n;
01055 
01056 /*                 Bidiagonalize A */
01057 /*                 (Workspace: need 3*N+M, prefer 3*N+(M+N)*NB) */
01058 
01059                     i__3 = *lwork - iwork + 1;
01060                     sgebrd_(m, n, &a[a_offset], lda, &s[1], &work[ie], &work[
01061                             itauq], &work[itaup], &work[iwork], &i__3, &ierr);
01062 
01063 /*                 Generate left vectors bidiagonalizing A */
01064 /*                 (Workspace: need 4*N, prefer 3*N+N*NB) */
01065 
01066                     i__3 = *lwork - iwork + 1;
01067                     sorgbr_("Q", m, n, n, &a[a_offset], lda, &work[itauq], &
01068                             work[iwork], &i__3, &ierr);
01069                     iwork = ie + *n;
01070 
01071 /*                 Perform bidiagonal QR iteration, computing left */
01072 /*                 singular vectors of A in A */
01073 /*                 (Workspace: need BDSPAC) */
01074 
01075                     sbdsqr_("U", n, &c__0, m, &c__0, &s[1], &work[ie], dum, &
01076                             c__1, &a[a_offset], lda, dum, &c__1, &work[iwork], 
01077                              info);
01078 
01079                 }
01080 
01081             } else if (wntuo && wntvas) {
01082 
01083 /*              Path 3 (M much larger than N, JOBU='O', JOBVT='S' or 'A') */
01084 /*              N left singular vectors to be overwritten on A and */
01085 /*              N right singular vectors to be computed in VT */
01086 
01087 /* Computing MAX */
01088                 i__3 = *n << 2;
01089                 if (*lwork >= *n * *n + max(i__3,bdspac)) {
01090 
01091 /*                 Sufficient workspace for a fast algorithm */
01092 
01093                     ir = 1;
01094 /* Computing MAX */
01095                     i__3 = wrkbl, i__2 = *lda * *n + *n;
01096                     if (*lwork >= max(i__3,i__2) + *lda * *n) {
01097 
01098 /*                    WORK(IU) is LDA by N and WORK(IR) is LDA by N */
01099 
01100                         ldwrku = *lda;
01101                         ldwrkr = *lda;
01102                     } else /* if(complicated condition) */ {
01103 /* Computing MAX */
01104                         i__3 = wrkbl, i__2 = *lda * *n + *n;
01105                         if (*lwork >= max(i__3,i__2) + *n * *n) {
01106 
01107 /*                    WORK(IU) is LDA by N and WORK(IR) is N by N */
01108 
01109                             ldwrku = *lda;
01110                             ldwrkr = *n;
01111                         } else {
01112 
01113 /*                    WORK(IU) is LDWRKU by N and WORK(IR) is N by N */
01114 
01115                             ldwrku = (*lwork - *n * *n - *n) / *n;
01116                             ldwrkr = *n;
01117                         }
01118                     }
01119                     itau = ir + ldwrkr * *n;
01120                     iwork = itau + *n;
01121 
01122 /*                 Compute A=Q*R */
01123 /*                 (Workspace: need N*N+2*N, prefer N*N+N+N*NB) */
01124 
01125                     i__3 = *lwork - iwork + 1;
01126                     sgeqrf_(m, n, &a[a_offset], lda, &work[itau], &work[iwork]
01127 , &i__3, &ierr);
01128 
01129 /*                 Copy R to VT, zeroing out below it */
01130 
01131                     slacpy_("U", n, n, &a[a_offset], lda, &vt[vt_offset], 
01132                             ldvt);
01133                     if (*n > 1) {
01134                         i__3 = *n - 1;
01135                         i__2 = *n - 1;
01136                         slaset_("L", &i__3, &i__2, &c_b421, &c_b421, &vt[
01137                                 vt_dim1 + 2], ldvt);
01138                     }
01139 
01140 /*                 Generate Q in A */
01141 /*                 (Workspace: need N*N+2*N, prefer N*N+N+N*NB) */
01142 
01143                     i__3 = *lwork - iwork + 1;
01144                     sorgqr_(m, n, n, &a[a_offset], lda, &work[itau], &work[
01145                             iwork], &i__3, &ierr);
01146                     ie = itau;
01147                     itauq = ie + *n;
01148                     itaup = itauq + *n;
01149                     iwork = itaup + *n;
01150 
01151 /*                 Bidiagonalize R in VT, copying result to WORK(IR) */
01152 /*                 (Workspace: need N*N+4*N, prefer N*N+3*N+2*N*NB) */
01153 
01154                     i__3 = *lwork - iwork + 1;
01155                     sgebrd_(n, n, &vt[vt_offset], ldvt, &s[1], &work[ie], &
01156                             work[itauq], &work[itaup], &work[iwork], &i__3, &
01157                             ierr);
01158                     slacpy_("L", n, n, &vt[vt_offset], ldvt, &work[ir], &
01159                             ldwrkr);
01160 
01161 /*                 Generate left vectors bidiagonalizing R in WORK(IR) */
01162 /*                 (Workspace: need N*N+4*N, prefer N*N+3*N+N*NB) */
01163 
01164                     i__3 = *lwork - iwork + 1;
01165                     sorgbr_("Q", n, n, n, &work[ir], &ldwrkr, &work[itauq], &
01166                             work[iwork], &i__3, &ierr);
01167 
01168 /*                 Generate right vectors bidiagonalizing R in VT */
01169 /*                 (Workspace: need N*N+4*N-1, prefer N*N+3*N+(N-1)*NB) */
01170 
01171                     i__3 = *lwork - iwork + 1;
01172                     sorgbr_("P", n, n, n, &vt[vt_offset], ldvt, &work[itaup], 
01173                             &work[iwork], &i__3, &ierr);
01174                     iwork = ie + *n;
01175 
01176 /*                 Perform bidiagonal QR iteration, computing left */
01177 /*                 singular vectors of R in WORK(IR) and computing right */
01178 /*                 singular vectors of R in VT */
01179 /*                 (Workspace: need N*N+BDSPAC) */
01180 
01181                     sbdsqr_("U", n, n, n, &c__0, &s[1], &work[ie], &vt[
01182                             vt_offset], ldvt, &work[ir], &ldwrkr, dum, &c__1, 
01183                             &work[iwork], info);
01184                     iu = ie + *n;
01185 
01186 /*                 Multiply Q in A by left singular vectors of R in */
01187 /*                 WORK(IR), storing result in WORK(IU) and copying to A */
01188 /*                 (Workspace: need N*N+2*N, prefer N*N+M*N+N) */
01189 
01190                     i__3 = *m;
01191                     i__2 = ldwrku;
01192                     for (i__ = 1; i__2 < 0 ? i__ >= i__3 : i__ <= i__3; i__ +=
01193                              i__2) {
01194 /* Computing MIN */
01195                         i__4 = *m - i__ + 1;
01196                         chunk = min(i__4,ldwrku);
01197                         sgemm_("N", "N", &chunk, n, n, &c_b443, &a[i__ + 
01198                                 a_dim1], lda, &work[ir], &ldwrkr, &c_b421, &
01199                                 work[iu], &ldwrku);
01200                         slacpy_("F", &chunk, n, &work[iu], &ldwrku, &a[i__ + 
01201                                 a_dim1], lda);
01202 /* L20: */
01203                     }
01204 
01205                 } else {
01206 
01207 /*                 Insufficient workspace for a fast algorithm */
01208 
01209                     itau = 1;
01210                     iwork = itau + *n;
01211 
01212 /*                 Compute A=Q*R */
01213 /*                 (Workspace: need 2*N, prefer N+N*NB) */
01214 
01215                     i__2 = *lwork - iwork + 1;
01216                     sgeqrf_(m, n, &a[a_offset], lda, &work[itau], &work[iwork]
01217 , &i__2, &ierr);
01218 
01219 /*                 Copy R to VT, zeroing out below it */
01220 
01221                     slacpy_("U", n, n, &a[a_offset], lda, &vt[vt_offset], 
01222                             ldvt);
01223                     if (*n > 1) {
01224                         i__2 = *n - 1;
01225                         i__3 = *n - 1;
01226                         slaset_("L", &i__2, &i__3, &c_b421, &c_b421, &vt[
01227                                 vt_dim1 + 2], ldvt);
01228                     }
01229 
01230 /*                 Generate Q in A */
01231 /*                 (Workspace: need 2*N, prefer N+N*NB) */
01232 
01233                     i__2 = *lwork - iwork + 1;
01234                     sorgqr_(m, n, n, &a[a_offset], lda, &work[itau], &work[
01235                             iwork], &i__2, &ierr);
01236                     ie = itau;
01237                     itauq = ie + *n;
01238                     itaup = itauq + *n;
01239                     iwork = itaup + *n;
01240 
01241 /*                 Bidiagonalize R in VT */
01242 /*                 (Workspace: need 4*N, prefer 3*N+2*N*NB) */
01243 
01244                     i__2 = *lwork - iwork + 1;
01245                     sgebrd_(n, n, &vt[vt_offset], ldvt, &s[1], &work[ie], &
01246                             work[itauq], &work[itaup], &work[iwork], &i__2, &
01247                             ierr);
01248 
01249 /*                 Multiply Q in A by left vectors bidiagonalizing R */
01250 /*                 (Workspace: need 3*N+M, prefer 3*N+M*NB) */
01251 
01252                     i__2 = *lwork - iwork + 1;
01253                     sormbr_("Q", "R", "N", m, n, n, &vt[vt_offset], ldvt, &
01254                             work[itauq], &a[a_offset], lda, &work[iwork], &
01255                             i__2, &ierr);
01256 
01257 /*                 Generate right vectors bidiagonalizing R in VT */
01258 /*                 (Workspace: need 4*N-1, prefer 3*N+(N-1)*NB) */
01259 
01260                     i__2 = *lwork - iwork + 1;
01261                     sorgbr_("P", n, n, n, &vt[vt_offset], ldvt, &work[itaup], 
01262                             &work[iwork], &i__2, &ierr);
01263                     iwork = ie + *n;
01264 
01265 /*                 Perform bidiagonal QR iteration, computing left */
01266 /*                 singular vectors of A in A and computing right */
01267 /*                 singular vectors of A in VT */
01268 /*                 (Workspace: need BDSPAC) */
01269 
01270                     sbdsqr_("U", n, n, m, &c__0, &s[1], &work[ie], &vt[
01271                             vt_offset], ldvt, &a[a_offset], lda, dum, &c__1, &
01272                             work[iwork], info);
01273 
01274                 }
01275 
01276             } else if (wntus) {
01277 
01278                 if (wntvn) {
01279 
01280 /*                 Path 4 (M much larger than N, JOBU='S', JOBVT='N') */
01281 /*                 N left singular vectors to be computed in U and */
01282 /*                 no right singular vectors to be computed */
01283 
01284 /* Computing MAX */
01285                     i__2 = *n << 2;
01286                     if (*lwork >= *n * *n + max(i__2,bdspac)) {
01287 
01288 /*                    Sufficient workspace for a fast algorithm */
01289 
01290                         ir = 1;
01291                         if (*lwork >= wrkbl + *lda * *n) {
01292 
01293 /*                       WORK(IR) is LDA by N */
01294 
01295                             ldwrkr = *lda;
01296                         } else {
01297 
01298 /*                       WORK(IR) is N by N */
01299 
01300                             ldwrkr = *n;
01301                         }
01302                         itau = ir + ldwrkr * *n;
01303                         iwork = itau + *n;
01304 
01305 /*                    Compute A=Q*R */
01306 /*                    (Workspace: need N*N+2*N, prefer N*N+N+N*NB) */
01307 
01308                         i__2 = *lwork - iwork + 1;
01309                         sgeqrf_(m, n, &a[a_offset], lda, &work[itau], &work[
01310                                 iwork], &i__2, &ierr);
01311 
01312 /*                    Copy R to WORK(IR), zeroing out below it */
01313 
01314                         slacpy_("U", n, n, &a[a_offset], lda, &work[ir], &
01315                                 ldwrkr);
01316                         i__2 = *n - 1;
01317                         i__3 = *n - 1;
01318                         slaset_("L", &i__2, &i__3, &c_b421, &c_b421, &work[ir 
01319                                 + 1], &ldwrkr);
01320 
01321 /*                    Generate Q in A */
01322 /*                    (Workspace: need N*N+2*N, prefer N*N+N+N*NB) */
01323 
01324                         i__2 = *lwork - iwork + 1;
01325                         sorgqr_(m, n, n, &a[a_offset], lda, &work[itau], &
01326                                 work[iwork], &i__2, &ierr);
01327                         ie = itau;
01328                         itauq = ie + *n;
01329                         itaup = itauq + *n;
01330                         iwork = itaup + *n;
01331 
01332 /*                    Bidiagonalize R in WORK(IR) */
01333 /*                    (Workspace: need N*N+4*N, prefer N*N+3*N+2*N*NB) */
01334 
01335                         i__2 = *lwork - iwork + 1;
01336                         sgebrd_(n, n, &work[ir], &ldwrkr, &s[1], &work[ie], &
01337                                 work[itauq], &work[itaup], &work[iwork], &
01338                                 i__2, &ierr);
01339 
01340 /*                    Generate left vectors bidiagonalizing R in WORK(IR) */
01341 /*                    (Workspace: need N*N+4*N, prefer N*N+3*N+N*NB) */
01342 
01343                         i__2 = *lwork - iwork + 1;
01344                         sorgbr_("Q", n, n, n, &work[ir], &ldwrkr, &work[itauq]
01345 , &work[iwork], &i__2, &ierr);
01346                         iwork = ie + *n;
01347 
01348 /*                    Perform bidiagonal QR iteration, computing left */
01349 /*                    singular vectors of R in WORK(IR) */
01350 /*                    (Workspace: need N*N+BDSPAC) */
01351 
01352                         sbdsqr_("U", n, &c__0, n, &c__0, &s[1], &work[ie], 
01353                                 dum, &c__1, &work[ir], &ldwrkr, dum, &c__1, &
01354                                 work[iwork], info);
01355 
01356 /*                    Multiply Q in A by left singular vectors of R in */
01357 /*                    WORK(IR), storing result in U */
01358 /*                    (Workspace: need N*N) */
01359 
01360                         sgemm_("N", "N", m, n, n, &c_b443, &a[a_offset], lda, 
01361                                 &work[ir], &ldwrkr, &c_b421, &u[u_offset], 
01362                                 ldu);
01363 
01364                     } else {
01365 
01366 /*                    Insufficient workspace for a fast algorithm */
01367 
01368                         itau = 1;
01369                         iwork = itau + *n;
01370 
01371 /*                    Compute A=Q*R, copying result to U */
01372 /*                    (Workspace: need 2*N, prefer N+N*NB) */
01373 
01374                         i__2 = *lwork - iwork + 1;
01375                         sgeqrf_(m, n, &a[a_offset], lda, &work[itau], &work[
01376                                 iwork], &i__2, &ierr);
01377                         slacpy_("L", m, n, &a[a_offset], lda, &u[u_offset], 
01378                                 ldu);
01379 
01380 /*                    Generate Q in U */
01381 /*                    (Workspace: need 2*N, prefer N+N*NB) */
01382 
01383                         i__2 = *lwork - iwork + 1;
01384                         sorgqr_(m, n, n, &u[u_offset], ldu, &work[itau], &
01385                                 work[iwork], &i__2, &ierr);
01386                         ie = itau;
01387                         itauq = ie + *n;
01388                         itaup = itauq + *n;
01389                         iwork = itaup + *n;
01390 
01391 /*                    Zero out below R in A */
01392 
01393                         i__2 = *n - 1;
01394                         i__3 = *n - 1;
01395                         slaset_("L", &i__2, &i__3, &c_b421, &c_b421, &a[
01396                                 a_dim1 + 2], lda);
01397 
01398 /*                    Bidiagonalize R in A */
01399 /*                    (Workspace: need 4*N, prefer 3*N+2*N*NB) */
01400 
01401                         i__2 = *lwork - iwork + 1;
01402                         sgebrd_(n, n, &a[a_offset], lda, &s[1], &work[ie], &
01403                                 work[itauq], &work[itaup], &work[iwork], &
01404                                 i__2, &ierr);
01405 
01406 /*                    Multiply Q in U by left vectors bidiagonalizing R */
01407 /*                    (Workspace: need 3*N+M, prefer 3*N+M*NB) */
01408 
01409                         i__2 = *lwork - iwork + 1;
01410                         sormbr_("Q", "R", "N", m, n, n, &a[a_offset], lda, &
01411                                 work[itauq], &u[u_offset], ldu, &work[iwork], 
01412                                 &i__2, &ierr)
01413                                 ;
01414                         iwork = ie + *n;
01415 
01416 /*                    Perform bidiagonal QR iteration, computing left */
01417 /*                    singular vectors of A in U */
01418 /*                    (Workspace: need BDSPAC) */
01419 
01420                         sbdsqr_("U", n, &c__0, m, &c__0, &s[1], &work[ie], 
01421                                 dum, &c__1, &u[u_offset], ldu, dum, &c__1, &
01422                                 work[iwork], info);
01423 
01424                     }
01425 
01426                 } else if (wntvo) {
01427 
01428 /*                 Path 5 (M much larger than N, JOBU='S', JOBVT='O') */
01429 /*                 N left singular vectors to be computed in U and */
01430 /*                 N right singular vectors to be overwritten on A */
01431 
01432 /* Computing MAX */
01433                     i__2 = *n << 2;
01434                     if (*lwork >= (*n << 1) * *n + max(i__2,bdspac)) {
01435 
01436 /*                    Sufficient workspace for a fast algorithm */
01437 
01438                         iu = 1;
01439                         if (*lwork >= wrkbl + (*lda << 1) * *n) {
01440 
01441 /*                       WORK(IU) is LDA by N and WORK(IR) is LDA by N */
01442 
01443                             ldwrku = *lda;
01444                             ir = iu + ldwrku * *n;
01445                             ldwrkr = *lda;
01446                         } else if (*lwork >= wrkbl + (*lda + *n) * *n) {
01447 
01448 /*                       WORK(IU) is LDA by N and WORK(IR) is N by N */
01449 
01450                             ldwrku = *lda;
01451                             ir = iu + ldwrku * *n;
01452                             ldwrkr = *n;
01453                         } else {
01454 
01455 /*                       WORK(IU) is N by N and WORK(IR) is N by N */
01456 
01457                             ldwrku = *n;
01458                             ir = iu + ldwrku * *n;
01459                             ldwrkr = *n;
01460                         }
01461                         itau = ir + ldwrkr * *n;
01462                         iwork = itau + *n;
01463 
01464 /*                    Compute A=Q*R */
01465 /*                    (Workspace: need 2*N*N+2*N, prefer 2*N*N+N+N*NB) */
01466 
01467                         i__2 = *lwork - iwork + 1;
01468                         sgeqrf_(m, n, &a[a_offset], lda, &work[itau], &work[
01469                                 iwork], &i__2, &ierr);
01470 
01471 /*                    Copy R to WORK(IU), zeroing out below it */
01472 
01473                         slacpy_("U", n, n, &a[a_offset], lda, &work[iu], &
01474                                 ldwrku);
01475                         i__2 = *n - 1;
01476                         i__3 = *n - 1;
01477                         slaset_("L", &i__2, &i__3, &c_b421, &c_b421, &work[iu 
01478                                 + 1], &ldwrku);
01479 
01480 /*                    Generate Q in A */
01481 /*                    (Workspace: need 2*N*N+2*N, prefer 2*N*N+N+N*NB) */
01482 
01483                         i__2 = *lwork - iwork + 1;
01484                         sorgqr_(m, n, n, &a[a_offset], lda, &work[itau], &
01485                                 work[iwork], &i__2, &ierr);
01486                         ie = itau;
01487                         itauq = ie + *n;
01488                         itaup = itauq + *n;
01489                         iwork = itaup + *n;
01490 
01491 /*                    Bidiagonalize R in WORK(IU), copying result to */
01492 /*                    WORK(IR) */
01493 /*                    (Workspace: need 2*N*N+4*N, */
01494 /*                                prefer 2*N*N+3*N+2*N*NB) */
01495 
01496                         i__2 = *lwork - iwork + 1;
01497                         sgebrd_(n, n, &work[iu], &ldwrku, &s[1], &work[ie], &
01498                                 work[itauq], &work[itaup], &work[iwork], &
01499                                 i__2, &ierr);
01500                         slacpy_("U", n, n, &work[iu], &ldwrku, &work[ir], &
01501                                 ldwrkr);
01502 
01503 /*                    Generate left bidiagonalizing vectors in WORK(IU) */
01504 /*                    (Workspace: need 2*N*N+4*N, prefer 2*N*N+3*N+N*NB) */
01505 
01506                         i__2 = *lwork - iwork + 1;
01507                         sorgbr_("Q", n, n, n, &work[iu], &ldwrku, &work[itauq]
01508 , &work[iwork], &i__2, &ierr);
01509 
01510 /*                    Generate right bidiagonalizing vectors in WORK(IR) */
01511 /*                    (Workspace: need 2*N*N+4*N-1, */
01512 /*                                prefer 2*N*N+3*N+(N-1)*NB) */
01513 
01514                         i__2 = *lwork - iwork + 1;
01515                         sorgbr_("P", n, n, n, &work[ir], &ldwrkr, &work[itaup]
01516 , &work[iwork], &i__2, &ierr);
01517                         iwork = ie + *n;
01518 
01519 /*                    Perform bidiagonal QR iteration, computing left */
01520 /*                    singular vectors of R in WORK(IU) and computing */
01521 /*                    right singular vectors of R in WORK(IR) */
01522 /*                    (Workspace: need 2*N*N+BDSPAC) */
01523 
01524                         sbdsqr_("U", n, n, n, &c__0, &s[1], &work[ie], &work[
01525                                 ir], &ldwrkr, &work[iu], &ldwrku, dum, &c__1, 
01526                                 &work[iwork], info);
01527 
01528 /*                    Multiply Q in A by left singular vectors of R in */
01529 /*                    WORK(IU), storing result in U */
01530 /*                    (Workspace: need N*N) */
01531 
01532                         sgemm_("N", "N", m, n, n, &c_b443, &a[a_offset], lda, 
01533                                 &work[iu], &ldwrku, &c_b421, &u[u_offset], 
01534                                 ldu);
01535 
01536 /*                    Copy right singular vectors of R to A */
01537 /*                    (Workspace: need N*N) */
01538 
01539                         slacpy_("F", n, n, &work[ir], &ldwrkr, &a[a_offset], 
01540                                 lda);
01541 
01542                     } else {
01543 
01544 /*                    Insufficient workspace for a fast algorithm */
01545 
01546                         itau = 1;
01547                         iwork = itau + *n;
01548 
01549 /*                    Compute A=Q*R, copying result to U */
01550 /*                    (Workspace: need 2*N, prefer N+N*NB) */
01551 
01552                         i__2 = *lwork - iwork + 1;
01553                         sgeqrf_(m, n, &a[a_offset], lda, &work[itau], &work[
01554                                 iwork], &i__2, &ierr);
01555                         slacpy_("L", m, n, &a[a_offset], lda, &u[u_offset], 
01556                                 ldu);
01557 
01558 /*                    Generate Q in U */
01559 /*                    (Workspace: need 2*N, prefer N+N*NB) */
01560 
01561                         i__2 = *lwork - iwork + 1;
01562                         sorgqr_(m, n, n, &u[u_offset], ldu, &work[itau], &
01563                                 work[iwork], &i__2, &ierr);
01564                         ie = itau;
01565                         itauq = ie + *n;
01566                         itaup = itauq + *n;
01567                         iwork = itaup + *n;
01568 
01569 /*                    Zero out below R in A */
01570 
01571                         i__2 = *n - 1;
01572                         i__3 = *n - 1;
01573                         slaset_("L", &i__2, &i__3, &c_b421, &c_b421, &a[
01574                                 a_dim1 + 2], lda);
01575 
01576 /*                    Bidiagonalize R in A */
01577 /*                    (Workspace: need 4*N, prefer 3*N+2*N*NB) */
01578 
01579                         i__2 = *lwork - iwork + 1;
01580                         sgebrd_(n, n, &a[a_offset], lda, &s[1], &work[ie], &
01581                                 work[itauq], &work[itaup], &work[iwork], &
01582                                 i__2, &ierr);
01583 
01584 /*                    Multiply Q in U by left vectors bidiagonalizing R */
01585 /*                    (Workspace: need 3*N+M, prefer 3*N+M*NB) */
01586 
01587                         i__2 = *lwork - iwork + 1;
01588                         sormbr_("Q", "R", "N", m, n, n, &a[a_offset], lda, &
01589                                 work[itauq], &u[u_offset], ldu, &work[iwork], 
01590                                 &i__2, &ierr)
01591                                 ;
01592 
01593 /*                    Generate right vectors bidiagonalizing R in A */
01594 /*                    (Workspace: need 4*N-1, prefer 3*N+(N-1)*NB) */
01595 
01596                         i__2 = *lwork - iwork + 1;
01597                         sorgbr_("P", n, n, n, &a[a_offset], lda, &work[itaup], 
01598                                  &work[iwork], &i__2, &ierr);
01599                         iwork = ie + *n;
01600 
01601 /*                    Perform bidiagonal QR iteration, computing left */
01602 /*                    singular vectors of A in U and computing right */
01603 /*                    singular vectors of A in A */
01604 /*                    (Workspace: need BDSPAC) */
01605 
01606                         sbdsqr_("U", n, n, m, &c__0, &s[1], &work[ie], &a[
01607                                 a_offset], lda, &u[u_offset], ldu, dum, &c__1, 
01608                                  &work[iwork], info);
01609 
01610                     }
01611 
01612                 } else if (wntvas) {
01613 
01614 /*                 Path 6 (M much larger than N, JOBU='S', JOBVT='S' */
01615 /*                         or 'A') */
01616 /*                 N left singular vectors to be computed in U and */
01617 /*                 N right singular vectors to be computed in VT */
01618 
01619 /* Computing MAX */
01620                     i__2 = *n << 2;
01621                     if (*lwork >= *n * *n + max(i__2,bdspac)) {
01622 
01623 /*                    Sufficient workspace for a fast algorithm */
01624 
01625                         iu = 1;
01626                         if (*lwork >= wrkbl + *lda * *n) {
01627 
01628 /*                       WORK(IU) is LDA by N */
01629 
01630                             ldwrku = *lda;
01631                         } else {
01632 
01633 /*                       WORK(IU) is N by N */
01634 
01635                             ldwrku = *n;
01636                         }
01637                         itau = iu + ldwrku * *n;
01638                         iwork = itau + *n;
01639 
01640 /*                    Compute A=Q*R */
01641 /*                    (Workspace: need N*N+2*N, prefer N*N+N+N*NB) */
01642 
01643                         i__2 = *lwork - iwork + 1;
01644                         sgeqrf_(m, n, &a[a_offset], lda, &work[itau], &work[
01645                                 iwork], &i__2, &ierr);
01646 
01647 /*                    Copy R to WORK(IU), zeroing out below it */
01648 
01649                         slacpy_("U", n, n, &a[a_offset], lda, &work[iu], &
01650                                 ldwrku);
01651                         i__2 = *n - 1;
01652                         i__3 = *n - 1;
01653                         slaset_("L", &i__2, &i__3, &c_b421, &c_b421, &work[iu 
01654                                 + 1], &ldwrku);
01655 
01656 /*                    Generate Q in A */
01657 /*                    (Workspace: need N*N+2*N, prefer N*N+N+N*NB) */
01658 
01659                         i__2 = *lwork - iwork + 1;
01660                         sorgqr_(m, n, n, &a[a_offset], lda, &work[itau], &
01661                                 work[iwork], &i__2, &ierr);
01662                         ie = itau;
01663                         itauq = ie + *n;
01664                         itaup = itauq + *n;
01665                         iwork = itaup + *n;
01666 
01667 /*                    Bidiagonalize R in WORK(IU), copying result to VT */
01668 /*                    (Workspace: need N*N+4*N, prefer N*N+3*N+2*N*NB) */
01669 
01670                         i__2 = *lwork - iwork + 1;
01671                         sgebrd_(n, n, &work[iu], &ldwrku, &s[1], &work[ie], &
01672                                 work[itauq], &work[itaup], &work[iwork], &
01673                                 i__2, &ierr);
01674                         slacpy_("U", n, n, &work[iu], &ldwrku, &vt[vt_offset], 
01675                                  ldvt);
01676 
01677 /*                    Generate left bidiagonalizing vectors in WORK(IU) */
01678 /*                    (Workspace: need N*N+4*N, prefer N*N+3*N+N*NB) */
01679 
01680                         i__2 = *lwork - iwork + 1;
01681                         sorgbr_("Q", n, n, n, &work[iu], &ldwrku, &work[itauq]
01682 , &work[iwork], &i__2, &ierr);
01683 
01684 /*                    Generate right bidiagonalizing vectors in VT */
01685 /*                    (Workspace: need N*N+4*N-1, */
01686 /*                                prefer N*N+3*N+(N-1)*NB) */
01687 
01688                         i__2 = *lwork - iwork + 1;
01689                         sorgbr_("P", n, n, n, &vt[vt_offset], ldvt, &work[
01690                                 itaup], &work[iwork], &i__2, &ierr)
01691                                 ;
01692                         iwork = ie + *n;
01693 
01694 /*                    Perform bidiagonal QR iteration, computing left */
01695 /*                    singular vectors of R in WORK(IU) and computing */
01696 /*                    right singular vectors of R in VT */
01697 /*                    (Workspace: need N*N+BDSPAC) */
01698 
01699                         sbdsqr_("U", n, n, n, &c__0, &s[1], &work[ie], &vt[
01700                                 vt_offset], ldvt, &work[iu], &ldwrku, dum, &
01701                                 c__1, &work[iwork], info);
01702 
01703 /*                    Multiply Q in A by left singular vectors of R in */
01704 /*                    WORK(IU), storing result in U */
01705 /*                    (Workspace: need N*N) */
01706 
01707                         sgemm_("N", "N", m, n, n, &c_b443, &a[a_offset], lda, 
01708                                 &work[iu], &ldwrku, &c_b421, &u[u_offset], 
01709                                 ldu);
01710 
01711                     } else {
01712 
01713 /*                    Insufficient workspace for a fast algorithm */
01714 
01715                         itau = 1;
01716                         iwork = itau + *n;
01717 
01718 /*                    Compute A=Q*R, copying result to U */
01719 /*                    (Workspace: need 2*N, prefer N+N*NB) */
01720 
01721                         i__2 = *lwork - iwork + 1;
01722                         sgeqrf_(m, n, &a[a_offset], lda, &work[itau], &work[
01723                                 iwork], &i__2, &ierr);
01724                         slacpy_("L", m, n, &a[a_offset], lda, &u[u_offset], 
01725                                 ldu);
01726 
01727 /*                    Generate Q in U */
01728 /*                    (Workspace: need 2*N, prefer N+N*NB) */
01729 
01730                         i__2 = *lwork - iwork + 1;
01731                         sorgqr_(m, n, n, &u[u_offset], ldu, &work[itau], &
01732                                 work[iwork], &i__2, &ierr);
01733 
01734 /*                    Copy R to VT, zeroing out below it */
01735 
01736                         slacpy_("U", n, n, &a[a_offset], lda, &vt[vt_offset], 
01737                                 ldvt);
01738                         if (*n > 1) {
01739                             i__2 = *n - 1;
01740                             i__3 = *n - 1;
01741                             slaset_("L", &i__2, &i__3, &c_b421, &c_b421, &vt[
01742                                     vt_dim1 + 2], ldvt);
01743                         }
01744                         ie = itau;
01745                         itauq = ie + *n;
01746                         itaup = itauq + *n;
01747                         iwork = itaup + *n;
01748 
01749 /*                    Bidiagonalize R in VT */
01750 /*                    (Workspace: need 4*N, prefer 3*N+2*N*NB) */
01751 
01752                         i__2 = *lwork - iwork + 1;
01753                         sgebrd_(n, n, &vt[vt_offset], ldvt, &s[1], &work[ie], 
01754                                 &work[itauq], &work[itaup], &work[iwork], &
01755                                 i__2, &ierr);
01756 
01757 /*                    Multiply Q in U by left bidiagonalizing vectors */
01758 /*                    in VT */
01759 /*                    (Workspace: need 3*N+M, prefer 3*N+M*NB) */
01760 
01761                         i__2 = *lwork - iwork + 1;
01762                         sormbr_("Q", "R", "N", m, n, n, &vt[vt_offset], ldvt, 
01763                                 &work[itauq], &u[u_offset], ldu, &work[iwork], 
01764                                  &i__2, &ierr);
01765 
01766 /*                    Generate right bidiagonalizing vectors in VT */
01767 /*                    (Workspace: need 4*N-1, prefer 3*N+(N-1)*NB) */
01768 
01769                         i__2 = *lwork - iwork + 1;
01770                         sorgbr_("P", n, n, n, &vt[vt_offset], ldvt, &work[
01771                                 itaup], &work[iwork], &i__2, &ierr)
01772                                 ;
01773                         iwork = ie + *n;
01774 
01775 /*                    Perform bidiagonal QR iteration, computing left */
01776 /*                    singular vectors of A in U and computing right */
01777 /*                    singular vectors of A in VT */
01778 /*                    (Workspace: need BDSPAC) */
01779 
01780                         sbdsqr_("U", n, n, m, &c__0, &s[1], &work[ie], &vt[
01781                                 vt_offset], ldvt, &u[u_offset], ldu, dum, &
01782                                 c__1, &work[iwork], info);
01783 
01784                     }
01785 
01786                 }
01787 
01788             } else if (wntua) {
01789 
01790                 if (wntvn) {
01791 
01792 /*                 Path 7 (M much larger than N, JOBU='A', JOBVT='N') */
01793 /*                 M left singular vectors to be computed in U and */
01794 /*                 no right singular vectors to be computed */
01795 
01796 /* Computing MAX */
01797                     i__2 = *n + *m, i__3 = *n << 2, i__2 = max(i__2,i__3);
01798                     if (*lwork >= *n * *n + max(i__2,bdspac)) {
01799 
01800 /*                    Sufficient workspace for a fast algorithm */
01801 
01802                         ir = 1;
01803                         if (*lwork >= wrkbl + *lda * *n) {
01804 
01805 /*                       WORK(IR) is LDA by N */
01806 
01807                             ldwrkr = *lda;
01808                         } else {
01809 
01810 /*                       WORK(IR) is N by N */
01811 
01812                             ldwrkr = *n;
01813                         }
01814                         itau = ir + ldwrkr * *n;
01815                         iwork = itau + *n;
01816 
01817 /*                    Compute A=Q*R, copying result to U */
01818 /*                    (Workspace: need N*N+2*N, prefer N*N+N+N*NB) */
01819 
01820                         i__2 = *lwork - iwork + 1;
01821                         sgeqrf_(m, n, &a[a_offset], lda, &work[itau], &work[
01822                                 iwork], &i__2, &ierr);
01823                         slacpy_("L", m, n, &a[a_offset], lda, &u[u_offset], 
01824                                 ldu);
01825 
01826 /*                    Copy R to WORK(IR), zeroing out below it */
01827 
01828                         slacpy_("U", n, n, &a[a_offset], lda, &work[ir], &
01829                                 ldwrkr);
01830                         i__2 = *n - 1;
01831                         i__3 = *n - 1;
01832                         slaset_("L", &i__2, &i__3, &c_b421, &c_b421, &work[ir 
01833                                 + 1], &ldwrkr);
01834 
01835 /*                    Generate Q in U */
01836 /*                    (Workspace: need N*N+N+M, prefer N*N+N+M*NB) */
01837 
01838                         i__2 = *lwork - iwork + 1;
01839                         sorgqr_(m, m, n, &u[u_offset], ldu, &work[itau], &
01840                                 work[iwork], &i__2, &ierr);
01841                         ie = itau;
01842                         itauq = ie + *n;
01843                         itaup = itauq + *n;
01844                         iwork = itaup + *n;
01845 
01846 /*                    Bidiagonalize R in WORK(IR) */
01847 /*                    (Workspace: need N*N+4*N, prefer N*N+3*N+2*N*NB) */
01848 
01849                         i__2 = *lwork - iwork + 1;
01850                         sgebrd_(n, n, &work[ir], &ldwrkr, &s[1], &work[ie], &
01851                                 work[itauq], &work[itaup], &work[iwork], &
01852                                 i__2, &ierr);
01853 
01854 /*                    Generate left bidiagonalizing vectors in WORK(IR) */
01855 /*                    (Workspace: need N*N+4*N, prefer N*N+3*N+N*NB) */
01856 
01857                         i__2 = *lwork - iwork + 1;
01858                         sorgbr_("Q", n, n, n, &work[ir], &ldwrkr, &work[itauq]
01859 , &work[iwork], &i__2, &ierr);
01860                         iwork = ie + *n;
01861 
01862 /*                    Perform bidiagonal QR iteration, computing left */
01863 /*                    singular vectors of R in WORK(IR) */
01864 /*                    (Workspace: need N*N+BDSPAC) */
01865 
01866                         sbdsqr_("U", n, &c__0, n, &c__0, &s[1], &work[ie], 
01867                                 dum, &c__1, &work[ir], &ldwrkr, dum, &c__1, &
01868                                 work[iwork], info);
01869 
01870 /*                    Multiply Q in U by left singular vectors of R in */
01871 /*                    WORK(IR), storing result in A */
01872 /*                    (Workspace: need N*N) */
01873 
01874                         sgemm_("N", "N", m, n, n, &c_b443, &u[u_offset], ldu, 
01875                                 &work[ir], &ldwrkr, &c_b421, &a[a_offset], 
01876                                 lda);
01877 
01878 /*                    Copy left singular vectors of A from A to U */
01879 
01880                         slacpy_("F", m, n, &a[a_offset], lda, &u[u_offset], 
01881                                 ldu);
01882 
01883                     } else {
01884 
01885 /*                    Insufficient workspace for a fast algorithm */
01886 
01887                         itau = 1;
01888                         iwork = itau + *n;
01889 
01890 /*                    Compute A=Q*R, copying result to U */
01891 /*                    (Workspace: need 2*N, prefer N+N*NB) */
01892 
01893                         i__2 = *lwork - iwork + 1;
01894                         sgeqrf_(m, n, &a[a_offset], lda, &work[itau], &work[
01895                                 iwork], &i__2, &ierr);
01896                         slacpy_("L", m, n, &a[a_offset], lda, &u[u_offset], 
01897                                 ldu);
01898 
01899 /*                    Generate Q in U */
01900 /*                    (Workspace: need N+M, prefer N+M*NB) */
01901 
01902                         i__2 = *lwork - iwork + 1;
01903                         sorgqr_(m, m, n, &u[u_offset], ldu, &work[itau], &
01904                                 work[iwork], &i__2, &ierr);
01905                         ie = itau;
01906                         itauq = ie + *n;
01907                         itaup = itauq + *n;
01908                         iwork = itaup + *n;
01909 
01910 /*                    Zero out below R in A */
01911 
01912                         i__2 = *n - 1;
01913                         i__3 = *n - 1;
01914                         slaset_("L", &i__2, &i__3, &c_b421, &c_b421, &a[
01915                                 a_dim1 + 2], lda);
01916 
01917 /*                    Bidiagonalize R in A */
01918 /*                    (Workspace: need 4*N, prefer 3*N+2*N*NB) */
01919 
01920                         i__2 = *lwork - iwork + 1;
01921                         sgebrd_(n, n, &a[a_offset], lda, &s[1], &work[ie], &
01922                                 work[itauq], &work[itaup], &work[iwork], &
01923                                 i__2, &ierr);
01924 
01925 /*                    Multiply Q in U by left bidiagonalizing vectors */
01926 /*                    in A */
01927 /*                    (Workspace: need 3*N+M, prefer 3*N+M*NB) */
01928 
01929                         i__2 = *lwork - iwork + 1;
01930                         sormbr_("Q", "R", "N", m, n, n, &a[a_offset], lda, &
01931                                 work[itauq], &u[u_offset], ldu, &work[iwork], 
01932                                 &i__2, &ierr)
01933                                 ;
01934                         iwork = ie + *n;
01935 
01936 /*                    Perform bidiagonal QR iteration, computing left */
01937 /*                    singular vectors of A in U */
01938 /*                    (Workspace: need BDSPAC) */
01939 
01940                         sbdsqr_("U", n, &c__0, m, &c__0, &s[1], &work[ie], 
01941                                 dum, &c__1, &u[u_offset], ldu, dum, &c__1, &
01942                                 work[iwork], info);
01943 
01944                     }
01945 
01946                 } else if (wntvo) {
01947 
01948 /*                 Path 8 (M much larger than N, JOBU='A', JOBVT='O') */
01949 /*                 M left singular vectors to be computed in U and */
01950 /*                 N right singular vectors to be overwritten on A */
01951 
01952 /* Computing MAX */
01953                     i__2 = *n + *m, i__3 = *n << 2, i__2 = max(i__2,i__3);
01954                     if (*lwork >= (*n << 1) * *n + max(i__2,bdspac)) {
01955 
01956 /*                    Sufficient workspace for a fast algorithm */
01957 
01958                         iu = 1;
01959                         if (*lwork >= wrkbl + (*lda << 1) * *n) {
01960 
01961 /*                       WORK(IU) is LDA by N and WORK(IR) is LDA by N */
01962 
01963                             ldwrku = *lda;
01964                             ir = iu + ldwrku * *n;
01965                             ldwrkr = *lda;
01966                         } else if (*lwork >= wrkbl + (*lda + *n) * *n) {
01967 
01968 /*                       WORK(IU) is LDA by N and WORK(IR) is N by N */
01969 
01970                             ldwrku = *lda;
01971                             ir = iu + ldwrku * *n;
01972                             ldwrkr = *n;
01973                         } else {
01974 
01975 /*                       WORK(IU) is N by N and WORK(IR) is N by N */
01976 
01977                             ldwrku = *n;
01978                             ir = iu + ldwrku * *n;
01979                             ldwrkr = *n;
01980                         }
01981                         itau = ir + ldwrkr * *n;
01982                         iwork = itau + *n;
01983 
01984 /*                    Compute A=Q*R, copying result to U */
01985 /*                    (Workspace: need 2*N*N+2*N, prefer 2*N*N+N+N*NB) */
01986 
01987                         i__2 = *lwork - iwork + 1;
01988                         sgeqrf_(m, n, &a[a_offset], lda, &work[itau], &work[
01989                                 iwork], &i__2, &ierr);
01990                         slacpy_("L", m, n, &a[a_offset], lda, &u[u_offset], 
01991                                 ldu);
01992 
01993 /*                    Generate Q in U */
01994 /*                    (Workspace: need 2*N*N+N+M, prefer 2*N*N+N+M*NB) */
01995 
01996                         i__2 = *lwork - iwork + 1;
01997                         sorgqr_(m, m, n, &u[u_offset], ldu, &work[itau], &
01998                                 work[iwork], &i__2, &ierr);
01999 
02000 /*                    Copy R to WORK(IU), zeroing out below it */
02001 
02002                         slacpy_("U", n, n, &a[a_offset], lda, &work[iu], &
02003                                 ldwrku);
02004                         i__2 = *n - 1;
02005                         i__3 = *n - 1;
02006                         slaset_("L", &i__2, &i__3, &c_b421, &c_b421, &work[iu 
02007                                 + 1], &ldwrku);
02008                         ie = itau;
02009                         itauq = ie + *n;
02010                         itaup = itauq + *n;
02011                         iwork = itaup + *n;
02012 
02013 /*                    Bidiagonalize R in WORK(IU), copying result to */
02014 /*                    WORK(IR) */
02015 /*                    (Workspace: need 2*N*N+4*N, */
02016 /*                                prefer 2*N*N+3*N+2*N*NB) */
02017 
02018                         i__2 = *lwork - iwork + 1;
02019                         sgebrd_(n, n, &work[iu], &ldwrku, &s[1], &work[ie], &
02020                                 work[itauq], &work[itaup], &work[iwork], &
02021                                 i__2, &ierr);
02022                         slacpy_("U", n, n, &work[iu], &ldwrku, &work[ir], &
02023                                 ldwrkr);
02024 
02025 /*                    Generate left bidiagonalizing vectors in WORK(IU) */
02026 /*                    (Workspace: need 2*N*N+4*N, prefer 2*N*N+3*N+N*NB) */
02027 
02028                         i__2 = *lwork - iwork + 1;
02029                         sorgbr_("Q", n, n, n, &work[iu], &ldwrku, &work[itauq]
02030 , &work[iwork], &i__2, &ierr);
02031 
02032 /*                    Generate right bidiagonalizing vectors in WORK(IR) */
02033 /*                    (Workspace: need 2*N*N+4*N-1, */
02034 /*                                prefer 2*N*N+3*N+(N-1)*NB) */
02035 
02036                         i__2 = *lwork - iwork + 1;
02037                         sorgbr_("P", n, n, n, &work[ir], &ldwrkr, &work[itaup]
02038 , &work[iwork], &i__2, &ierr);
02039                         iwork = ie + *n;
02040 
02041 /*                    Perform bidiagonal QR iteration, computing left */
02042 /*                    singular vectors of R in WORK(IU) and computing */
02043 /*                    right singular vectors of R in WORK(IR) */
02044 /*                    (Workspace: need 2*N*N+BDSPAC) */
02045 
02046                         sbdsqr_("U", n, n, n, &c__0, &s[1], &work[ie], &work[
02047                                 ir], &ldwrkr, &work[iu], &ldwrku, dum, &c__1, 
02048                                 &work[iwork], info);
02049 
02050 /*                    Multiply Q in U by left singular vectors of R in */
02051 /*                    WORK(IU), storing result in A */
02052 /*                    (Workspace: need N*N) */
02053 
02054                         sgemm_("N", "N", m, n, n, &c_b443, &u[u_offset], ldu, 
02055                                 &work[iu], &ldwrku, &c_b421, &a[a_offset], 
02056                                 lda);
02057 
02058 /*                    Copy left singular vectors of A from A to U */
02059 
02060                         slacpy_("F", m, n, &a[a_offset], lda, &u[u_offset], 
02061                                 ldu);
02062 
02063 /*                    Copy right singular vectors of R from WORK(IR) to A */
02064 
02065                         slacpy_("F", n, n, &work[ir], &ldwrkr, &a[a_offset], 
02066                                 lda);
02067 
02068                     } else {
02069 
02070 /*                    Insufficient workspace for a fast algorithm */
02071 
02072                         itau = 1;
02073                         iwork = itau + *n;
02074 
02075 /*                    Compute A=Q*R, copying result to U */
02076 /*                    (Workspace: need 2*N, prefer N+N*NB) */
02077 
02078                         i__2 = *lwork - iwork + 1;
02079                         sgeqrf_(m, n, &a[a_offset], lda, &work[itau], &work[
02080                                 iwork], &i__2, &ierr);
02081                         slacpy_("L", m, n, &a[a_offset], lda, &u[u_offset], 
02082                                 ldu);
02083 
02084 /*                    Generate Q in U */
02085 /*                    (Workspace: need N+M, prefer N+M*NB) */
02086 
02087                         i__2 = *lwork - iwork + 1;
02088                         sorgqr_(m, m, n, &u[u_offset], ldu, &work[itau], &
02089                                 work[iwork], &i__2, &ierr);
02090                         ie = itau;
02091                         itauq = ie + *n;
02092                         itaup = itauq + *n;
02093                         iwork = itaup + *n;
02094 
02095 /*                    Zero out below R in A */
02096 
02097                         i__2 = *n - 1;
02098                         i__3 = *n - 1;
02099                         slaset_("L", &i__2, &i__3, &c_b421, &c_b421, &a[
02100                                 a_dim1 + 2], lda);
02101 
02102 /*                    Bidiagonalize R in A */
02103 /*                    (Workspace: need 4*N, prefer 3*N+2*N*NB) */
02104 
02105                         i__2 = *lwork - iwork + 1;
02106                         sgebrd_(n, n, &a[a_offset], lda, &s[1], &work[ie], &
02107                                 work[itauq], &work[itaup], &work[iwork], &
02108                                 i__2, &ierr);
02109 
02110 /*                    Multiply Q in U by left bidiagonalizing vectors */
02111 /*                    in A */
02112 /*                    (Workspace: need 3*N+M, prefer 3*N+M*NB) */
02113 
02114                         i__2 = *lwork - iwork + 1;
02115                         sormbr_("Q", "R", "N", m, n, n, &a[a_offset], lda, &
02116                                 work[itauq], &u[u_offset], ldu, &work[iwork], 
02117                                 &i__2, &ierr)
02118                                 ;
02119 
02120 /*                    Generate right bidiagonalizing vectors in A */
02121 /*                    (Workspace: need 4*N-1, prefer 3*N+(N-1)*NB) */
02122 
02123                         i__2 = *lwork - iwork + 1;
02124                         sorgbr_("P", n, n, n, &a[a_offset], lda, &work[itaup], 
02125                                  &work[iwork], &i__2, &ierr);
02126                         iwork = ie + *n;
02127 
02128 /*                    Perform bidiagonal QR iteration, computing left */
02129 /*                    singular vectors of A in U and computing right */
02130 /*                    singular vectors of A in A */
02131 /*                    (Workspace: need BDSPAC) */
02132 
02133                         sbdsqr_("U", n, n, m, &c__0, &s[1], &work[ie], &a[
02134                                 a_offset], lda, &u[u_offset], ldu, dum, &c__1, 
02135                                  &work[iwork], info);
02136 
02137                     }
02138 
02139                 } else if (wntvas) {
02140 
02141 /*                 Path 9 (M much larger than N, JOBU='A', JOBVT='S' */
02142 /*                         or 'A') */
02143 /*                 M left singular vectors to be computed in U and */
02144 /*                 N right singular vectors to be computed in VT */
02145 
02146 /* Computing MAX */
02147                     i__2 = *n + *m, i__3 = *n << 2, i__2 = max(i__2,i__3);
02148                     if (*lwork >= *n * *n + max(i__2,bdspac)) {
02149 
02150 /*                    Sufficient workspace for a fast algorithm */
02151 
02152                         iu = 1;
02153                         if (*lwork >= wrkbl + *lda * *n) {
02154 
02155 /*                       WORK(IU) is LDA by N */
02156 
02157                             ldwrku = *lda;
02158                         } else {
02159 
02160 /*                       WORK(IU) is N by N */
02161 
02162                             ldwrku = *n;
02163                         }
02164                         itau = iu + ldwrku * *n;
02165                         iwork = itau + *n;
02166 
02167 /*                    Compute A=Q*R, copying result to U */
02168 /*                    (Workspace: need N*N+2*N, prefer N*N+N+N*NB) */
02169 
02170                         i__2 = *lwork - iwork + 1;
02171                         sgeqrf_(m, n, &a[a_offset], lda, &work[itau], &work[
02172                                 iwork], &i__2, &ierr);
02173                         slacpy_("L", m, n, &a[a_offset], lda, &u[u_offset], 
02174                                 ldu);
02175 
02176 /*                    Generate Q in U */
02177 /*                    (Workspace: need N*N+N+M, prefer N*N+N+M*NB) */
02178 
02179                         i__2 = *lwork - iwork + 1;
02180                         sorgqr_(m, m, n, &u[u_offset], ldu, &work[itau], &
02181                                 work[iwork], &i__2, &ierr);
02182 
02183 /*                    Copy R to WORK(IU), zeroing out below it */
02184 
02185                         slacpy_("U", n, n, &a[a_offset], lda, &work[iu], &
02186                                 ldwrku);
02187                         i__2 = *n - 1;
02188                         i__3 = *n - 1;
02189                         slaset_("L", &i__2, &i__3, &c_b421, &c_b421, &work[iu 
02190                                 + 1], &ldwrku);
02191                         ie = itau;
02192                         itauq = ie + *n;
02193                         itaup = itauq + *n;
02194                         iwork = itaup + *n;
02195 
02196 /*                    Bidiagonalize R in WORK(IU), copying result to VT */
02197 /*                    (Workspace: need N*N+4*N, prefer N*N+3*N+2*N*NB) */
02198 
02199                         i__2 = *lwork - iwork + 1;
02200                         sgebrd_(n, n, &work[iu], &ldwrku, &s[1], &work[ie], &
02201                                 work[itauq], &work[itaup], &work[iwork], &
02202                                 i__2, &ierr);
02203                         slacpy_("U", n, n, &work[iu], &ldwrku, &vt[vt_offset], 
02204                                  ldvt);
02205 
02206 /*                    Generate left bidiagonalizing vectors in WORK(IU) */
02207 /*                    (Workspace: need N*N+4*N, prefer N*N+3*N+N*NB) */
02208 
02209                         i__2 = *lwork - iwork + 1;
02210                         sorgbr_("Q", n, n, n, &work[iu], &ldwrku, &work[itauq]
02211 , &work[iwork], &i__2, &ierr);
02212 
02213 /*                    Generate right bidiagonalizing vectors in VT */
02214 /*                    (Workspace: need N*N+4*N-1, */
02215 /*                                prefer N*N+3*N+(N-1)*NB) */
02216 
02217                         i__2 = *lwork - iwork + 1;
02218                         sorgbr_("P", n, n, n, &vt[vt_offset], ldvt, &work[
02219                                 itaup], &work[iwork], &i__2, &ierr)
02220                                 ;
02221                         iwork = ie + *n;
02222 
02223 /*                    Perform bidiagonal QR iteration, computing left */
02224 /*                    singular vectors of R in WORK(IU) and computing */
02225 /*                    right singular vectors of R in VT */
02226 /*                    (Workspace: need N*N+BDSPAC) */
02227 
02228                         sbdsqr_("U", n, n, n, &c__0, &s[1], &work[ie], &vt[
02229                                 vt_offset], ldvt, &work[iu], &ldwrku, dum, &
02230                                 c__1, &work[iwork], info);
02231 
02232 /*                    Multiply Q in U by left singular vectors of R in */
02233 /*                    WORK(IU), storing result in A */
02234 /*                    (Workspace: need N*N) */
02235 
02236                         sgemm_("N", "N", m, n, n, &c_b443, &u[u_offset], ldu, 
02237                                 &work[iu], &ldwrku, &c_b421, &a[a_offset], 
02238                                 lda);
02239 
02240 /*                    Copy left singular vectors of A from A to U */
02241 
02242                         slacpy_("F", m, n, &a[a_offset], lda, &u[u_offset], 
02243                                 ldu);
02244 
02245                     } else {
02246 
02247 /*                    Insufficient workspace for a fast algorithm */
02248 
02249                         itau = 1;
02250                         iwork = itau + *n;
02251 
02252 /*                    Compute A=Q*R, copying result to U */
02253 /*                    (Workspace: need 2*N, prefer N+N*NB) */
02254 
02255                         i__2 = *lwork - iwork + 1;
02256                         sgeqrf_(m, n, &a[a_offset], lda, &work[itau], &work[
02257                                 iwork], &i__2, &ierr);
02258                         slacpy_("L", m, n, &a[a_offset], lda, &u[u_offset], 
02259                                 ldu);
02260 
02261 /*                    Generate Q in U */
02262 /*                    (Workspace: need N+M, prefer N+M*NB) */
02263 
02264                         i__2 = *lwork - iwork + 1;
02265                         sorgqr_(m, m, n, &u[u_offset], ldu, &work[itau], &
02266                                 work[iwork], &i__2, &ierr);
02267 
02268 /*                    Copy R from A to VT, zeroing out below it */
02269 
02270                         slacpy_("U", n, n, &a[a_offset], lda, &vt[vt_offset], 
02271                                 ldvt);
02272                         if (*n > 1) {
02273                             i__2 = *n - 1;
02274                             i__3 = *n - 1;
02275                             slaset_("L", &i__2, &i__3, &c_b421, &c_b421, &vt[
02276                                     vt_dim1 + 2], ldvt);
02277                         }
02278                         ie = itau;
02279                         itauq = ie + *n;
02280                         itaup = itauq + *n;
02281                         iwork = itaup + *n;
02282 
02283 /*                    Bidiagonalize R in VT */
02284 /*                    (Workspace: need 4*N, prefer 3*N+2*N*NB) */
02285 
02286                         i__2 = *lwork - iwork + 1;
02287                         sgebrd_(n, n, &vt[vt_offset], ldvt, &s[1], &work[ie], 
02288                                 &work[itauq], &work[itaup], &work[iwork], &
02289                                 i__2, &ierr);
02290 
02291 /*                    Multiply Q in U by left bidiagonalizing vectors */
02292 /*                    in VT */
02293 /*                    (Workspace: need 3*N+M, prefer 3*N+M*NB) */
02294 
02295                         i__2 = *lwork - iwork + 1;
02296                         sormbr_("Q", "R", "N", m, n, n, &vt[vt_offset], ldvt, 
02297                                 &work[itauq], &u[u_offset], ldu, &work[iwork], 
02298                                  &i__2, &ierr);
02299 
02300 /*                    Generate right bidiagonalizing vectors in VT */
02301 /*                    (Workspace: need 4*N-1, prefer 3*N+(N-1)*NB) */
02302 
02303                         i__2 = *lwork - iwork + 1;
02304                         sorgbr_("P", n, n, n, &vt[vt_offset], ldvt, &work[
02305                                 itaup], &work[iwork], &i__2, &ierr)
02306                                 ;
02307                         iwork = ie + *n;
02308 
02309 /*                    Perform bidiagonal QR iteration, computing left */
02310 /*                    singular vectors of A in U and computing right */
02311 /*                    singular vectors of A in VT */
02312 /*                    (Workspace: need BDSPAC) */
02313 
02314                         sbdsqr_("U", n, n, m, &c__0, &s[1], &work[ie], &vt[
02315                                 vt_offset], ldvt, &u[u_offset], ldu, dum, &
02316                                 c__1, &work[iwork], info);
02317 
02318                     }
02319 
02320                 }
02321 
02322             }
02323 
02324         } else {
02325 
02326 /*           M .LT. MNTHR */
02327 
02328 /*           Path 10 (M at least N, but not much larger) */
02329 /*           Reduce to bidiagonal form without QR decomposition */
02330 
02331             ie = 1;
02332             itauq = ie + *n;
02333             itaup = itauq + *n;
02334             iwork = itaup + *n;
02335 
02336 /*           Bidiagonalize A */
02337 /*           (Workspace: need 3*N+M, prefer 3*N+(M+N)*NB) */
02338 
02339             i__2 = *lwork - iwork + 1;
02340             sgebrd_(m, n, &a[a_offset], lda, &s[1], &work[ie], &work[itauq], &
02341                     work[itaup], &work[iwork], &i__2, &ierr);
02342             if (wntuas) {
02343 
02344 /*              If left singular vectors desired in U, copy result to U */
02345 /*              and generate left bidiagonalizing vectors in U */
02346 /*              (Workspace: need 3*N+NCU, prefer 3*N+NCU*NB) */
02347 
02348                 slacpy_("L", m, n, &a[a_offset], lda, &u[u_offset], ldu);
02349                 if (wntus) {
02350                     ncu = *n;
02351                 }
02352                 if (wntua) {
02353                     ncu = *m;
02354                 }
02355                 i__2 = *lwork - iwork + 1;
02356                 sorgbr_("Q", m, &ncu, n, &u[u_offset], ldu, &work[itauq], &
02357                         work[iwork], &i__2, &ierr);
02358             }
02359             if (wntvas) {
02360 
02361 /*              If right singular vectors desired in VT, copy result to */
02362 /*              VT and generate right bidiagonalizing vectors in VT */
02363 /*              (Workspace: need 4*N-1, prefer 3*N+(N-1)*NB) */
02364 
02365                 slacpy_("U", n, n, &a[a_offset], lda, &vt[vt_offset], ldvt);
02366                 i__2 = *lwork - iwork + 1;
02367                 sorgbr_("P", n, n, n, &vt[vt_offset], ldvt, &work[itaup], &
02368                         work[iwork], &i__2, &ierr);
02369             }
02370             if (wntuo) {
02371 
02372 /*              If left singular vectors desired in A, generate left */
02373 /*              bidiagonalizing vectors in A */
02374 /*              (Workspace: need 4*N, prefer 3*N+N*NB) */
02375 
02376                 i__2 = *lwork - iwork + 1;
02377                 sorgbr_("Q", m, n, n, &a[a_offset], lda, &work[itauq], &work[
02378                         iwork], &i__2, &ierr);
02379             }
02380             if (wntvo) {
02381 
02382 /*              If right singular vectors desired in A, generate right */
02383 /*              bidiagonalizing vectors in A */
02384 /*              (Workspace: need 4*N-1, prefer 3*N+(N-1)*NB) */
02385 
02386                 i__2 = *lwork - iwork + 1;
02387                 sorgbr_("P", n, n, n, &a[a_offset], lda, &work[itaup], &work[
02388                         iwork], &i__2, &ierr);
02389             }
02390             iwork = ie + *n;
02391             if (wntuas || wntuo) {
02392                 nru = *m;
02393             }
02394             if (wntun) {
02395                 nru = 0;
02396             }
02397             if (wntvas || wntvo) {
02398                 ncvt = *n;
02399             }
02400             if (wntvn) {
02401                 ncvt = 0;
02402             }
02403             if (! wntuo && ! wntvo) {
02404 
02405 /*              Perform bidiagonal QR iteration, if desired, computing */
02406 /*              left singular vectors in U and computing right singular */
02407 /*              vectors in VT */
02408 /*              (Workspace: need BDSPAC) */
02409 
02410                 sbdsqr_("U", n, &ncvt, &nru, &c__0, &s[1], &work[ie], &vt[
02411                         vt_offset], ldvt, &u[u_offset], ldu, dum, &c__1, &
02412                         work[iwork], info);
02413             } else if (! wntuo && wntvo) {
02414 
02415 /*              Perform bidiagonal QR iteration, if desired, computing */
02416 /*              left singular vectors in U and computing right singular */
02417 /*              vectors in A */
02418 /*              (Workspace: need BDSPAC) */
02419 
02420                 sbdsqr_("U", n, &ncvt, &nru, &c__0, &s[1], &work[ie], &a[
02421                         a_offset], lda, &u[u_offset], ldu, dum, &c__1, &work[
02422                         iwork], info);
02423             } else {
02424 
02425 /*              Perform bidiagonal QR iteration, if desired, computing */
02426 /*              left singular vectors in A and computing right singular */
02427 /*              vectors in VT */
02428 /*              (Workspace: need BDSPAC) */
02429 
02430                 sbdsqr_("U", n, &ncvt, &nru, &c__0, &s[1], &work[ie], &vt[
02431                         vt_offset], ldvt, &a[a_offset], lda, dum, &c__1, &
02432                         work[iwork], info);
02433             }
02434 
02435         }
02436 
02437     } else {
02438 
02439 /*        A has more columns than rows. If A has sufficiently more */
02440 /*        columns than rows, first reduce using the LQ decomposition (if */
02441 /*        sufficient workspace available) */
02442 
02443         if (*n >= mnthr) {
02444 
02445             if (wntvn) {
02446 
02447 /*              Path 1t(N much larger than M, JOBVT='N') */
02448 /*              No right singular vectors to be computed */
02449 
02450                 itau = 1;
02451                 iwork = itau + *m;
02452 
02453 /*              Compute A=L*Q */
02454 /*              (Workspace: need 2*M, prefer M+M*NB) */
02455 
02456                 i__2 = *lwork - iwork + 1;
02457                 sgelqf_(m, n, &a[a_offset], lda, &work[itau], &work[iwork], &
02458                         i__2, &ierr);
02459 
02460 /*              Zero out above L */
02461 
02462                 i__2 = *m - 1;
02463                 i__3 = *m - 1;
02464                 slaset_("U", &i__2, &i__3, &c_b421, &c_b421, &a[(a_dim1 << 1) 
02465                         + 1], lda);
02466                 ie = 1;
02467                 itauq = ie + *m;
02468                 itaup = itauq + *m;
02469                 iwork = itaup + *m;
02470 
02471 /*              Bidiagonalize L in A */
02472 /*              (Workspace: need 4*M, prefer 3*M+2*M*NB) */
02473 
02474                 i__2 = *lwork - iwork + 1;
02475                 sgebrd_(m, m, &a[a_offset], lda, &s[1], &work[ie], &work[
02476                         itauq], &work[itaup], &work[iwork], &i__2, &ierr);
02477                 if (wntuo || wntuas) {
02478 
02479 /*                 If left singular vectors desired, generate Q */
02480 /*                 (Workspace: need 4*M, prefer 3*M+M*NB) */
02481 
02482                     i__2 = *lwork - iwork + 1;
02483                     sorgbr_("Q", m, m, m, &a[a_offset], lda, &work[itauq], &
02484                             work[iwork], &i__2, &ierr);
02485                 }
02486                 iwork = ie + *m;
02487                 nru = 0;
02488                 if (wntuo || wntuas) {
02489                     nru = *m;
02490                 }
02491 
02492 /*              Perform bidiagonal QR iteration, computing left singular */
02493 /*              vectors of A in A if desired */
02494 /*              (Workspace: need BDSPAC) */
02495 
02496                 sbdsqr_("U", m, &c__0, &nru, &c__0, &s[1], &work[ie], dum, &
02497                         c__1, &a[a_offset], lda, dum, &c__1, &work[iwork], 
02498                         info);
02499 
02500 /*              If left singular vectors desired in U, copy them there */
02501 
02502                 if (wntuas) {
02503                     slacpy_("F", m, m, &a[a_offset], lda, &u[u_offset], ldu);
02504                 }
02505 
02506             } else if (wntvo && wntun) {
02507 
02508 /*              Path 2t(N much larger than M, JOBU='N', JOBVT='O') */
02509 /*              M right singular vectors to be overwritten on A and */
02510 /*              no left singular vectors to be computed */
02511 
02512 /* Computing MAX */
02513                 i__2 = *m << 2;
02514                 if (*lwork >= *m * *m + max(i__2,bdspac)) {
02515 
02516 /*                 Sufficient workspace for a fast algorithm */
02517 
02518                     ir = 1;
02519 /* Computing MAX */
02520                     i__2 = wrkbl, i__3 = *lda * *n + *m;
02521                     if (*lwork >= max(i__2,i__3) + *lda * *m) {
02522 
02523 /*                    WORK(IU) is LDA by N and WORK(IR) is LDA by M */
02524 
02525                         ldwrku = *lda;
02526                         chunk = *n;
02527                         ldwrkr = *lda;
02528                     } else /* if(complicated condition) */ {
02529 /* Computing MAX */
02530                         i__2 = wrkbl, i__3 = *lda * *n + *m;
02531                         if (*lwork >= max(i__2,i__3) + *m * *m) {
02532 
02533 /*                    WORK(IU) is LDA by N and WORK(IR) is M by M */
02534 
02535                             ldwrku = *lda;
02536                             chunk = *n;
02537                             ldwrkr = *m;
02538                         } else {
02539 
02540 /*                    WORK(IU) is M by CHUNK and WORK(IR) is M by M */
02541 
02542                             ldwrku = *m;
02543                             chunk = (*lwork - *m * *m - *m) / *m;
02544                             ldwrkr = *m;
02545                         }
02546                     }
02547                     itau = ir + ldwrkr * *m;
02548                     iwork = itau + *m;
02549 
02550 /*                 Compute A=L*Q */
02551 /*                 (Workspace: need M*M+2*M, prefer M*M+M+M*NB) */
02552 
02553                     i__2 = *lwork - iwork + 1;
02554                     sgelqf_(m, n, &a[a_offset], lda, &work[itau], &work[iwork]
02555 , &i__2, &ierr);
02556 
02557 /*                 Copy L to WORK(IR) and zero out above it */
02558 
02559                     slacpy_("L", m, m, &a[a_offset], lda, &work[ir], &ldwrkr);
02560                     i__2 = *m - 1;
02561                     i__3 = *m - 1;
02562                     slaset_("U", &i__2, &i__3, &c_b421, &c_b421, &work[ir + 
02563                             ldwrkr], &ldwrkr);
02564 
02565 /*                 Generate Q in A */
02566 /*                 (Workspace: need M*M+2*M, prefer M*M+M+M*NB) */
02567 
02568                     i__2 = *lwork - iwork + 1;
02569                     sorglq_(m, n, m, &a[a_offset], lda, &work[itau], &work[
02570                             iwork], &i__2, &ierr);
02571                     ie = itau;
02572                     itauq = ie + *m;
02573                     itaup = itauq + *m;
02574                     iwork = itaup + *m;
02575 
02576 /*                 Bidiagonalize L in WORK(IR) */
02577 /*                 (Workspace: need M*M+4*M, prefer M*M+3*M+2*M*NB) */
02578 
02579                     i__2 = *lwork - iwork + 1;
02580                     sgebrd_(m, m, &work[ir], &ldwrkr, &s[1], &work[ie], &work[
02581                             itauq], &work[itaup], &work[iwork], &i__2, &ierr);
02582 
02583 /*                 Generate right vectors bidiagonalizing L */
02584 /*                 (Workspace: need M*M+4*M-1, prefer M*M+3*M+(M-1)*NB) */
02585 
02586                     i__2 = *lwork - iwork + 1;
02587                     sorgbr_("P", m, m, m, &work[ir], &ldwrkr, &work[itaup], &
02588                             work[iwork], &i__2, &ierr);
02589                     iwork = ie + *m;
02590 
02591 /*                 Perform bidiagonal QR iteration, computing right */
02592 /*                 singular vectors of L in WORK(IR) */
02593 /*                 (Workspace: need M*M+BDSPAC) */
02594 
02595                     sbdsqr_("U", m, m, &c__0, &c__0, &s[1], &work[ie], &work[
02596                             ir], &ldwrkr, dum, &c__1, dum, &c__1, &work[iwork]
02597 , info);
02598                     iu = ie + *m;
02599 
02600 /*                 Multiply right singular vectors of L in WORK(IR) by Q */
02601 /*                 in A, storing result in WORK(IU) and copying to A */
02602 /*                 (Workspace: need M*M+2*M, prefer M*M+M*N+M) */
02603 
02604                     i__2 = *n;
02605                     i__3 = chunk;
02606                     for (i__ = 1; i__3 < 0 ? i__ >= i__2 : i__ <= i__2; i__ +=
02607                              i__3) {
02608 /* Computing MIN */
02609                         i__4 = *n - i__ + 1;
02610                         blk = min(i__4,chunk);
02611                         sgemm_("N", "N", m, &blk, m, &c_b443, &work[ir], &
02612                                 ldwrkr, &a[i__ * a_dim1 + 1], lda, &c_b421, &
02613                                 work[iu], &ldwrku);
02614                         slacpy_("F", m, &blk, &work[iu], &ldwrku, &a[i__ * 
02615                                 a_dim1 + 1], lda);
02616 /* L30: */
02617                     }
02618 
02619                 } else {
02620 
02621 /*                 Insufficient workspace for a fast algorithm */
02622 
02623                     ie = 1;
02624                     itauq = ie + *m;
02625                     itaup = itauq + *m;
02626                     iwork = itaup + *m;
02627 
02628 /*                 Bidiagonalize A */
02629 /*                 (Workspace: need 3*M+N, prefer 3*M+(M+N)*NB) */
02630 
02631                     i__3 = *lwork - iwork + 1;
02632                     sgebrd_(m, n, &a[a_offset], lda, &s[1], &work[ie], &work[
02633                             itauq], &work[itaup], &work[iwork], &i__3, &ierr);
02634 
02635 /*                 Generate right vectors bidiagonalizing A */
02636 /*                 (Workspace: need 4*M, prefer 3*M+M*NB) */
02637 
02638                     i__3 = *lwork - iwork + 1;
02639                     sorgbr_("P", m, n, m, &a[a_offset], lda, &work[itaup], &
02640                             work[iwork], &i__3, &ierr);
02641                     iwork = ie + *m;
02642 
02643 /*                 Perform bidiagonal QR iteration, computing right */
02644 /*                 singular vectors of A in A */
02645 /*                 (Workspace: need BDSPAC) */
02646 
02647                     sbdsqr_("L", m, n, &c__0, &c__0, &s[1], &work[ie], &a[
02648                             a_offset], lda, dum, &c__1, dum, &c__1, &work[
02649                             iwork], info);
02650 
02651                 }
02652 
02653             } else if (wntvo && wntuas) {
02654 
02655 /*              Path 3t(N much larger than M, JOBU='S' or 'A', JOBVT='O') */
02656 /*              M right singular vectors to be overwritten on A and */
02657 /*              M left singular vectors to be computed in U */
02658 
02659 /* Computing MAX */
02660                 i__3 = *m << 2;
02661                 if (*lwork >= *m * *m + max(i__3,bdspac)) {
02662 
02663 /*                 Sufficient workspace for a fast algorithm */
02664 
02665                     ir = 1;
02666 /* Computing MAX */
02667                     i__3 = wrkbl, i__2 = *lda * *n + *m;
02668                     if (*lwork >= max(i__3,i__2) + *lda * *m) {
02669 
02670 /*                    WORK(IU) is LDA by N and WORK(IR) is LDA by M */
02671 
02672                         ldwrku = *lda;
02673                         chunk = *n;
02674                         ldwrkr = *lda;
02675                     } else /* if(complicated condition) */ {
02676 /* Computing MAX */
02677                         i__3 = wrkbl, i__2 = *lda * *n + *m;
02678                         if (*lwork >= max(i__3,i__2) + *m * *m) {
02679 
02680 /*                    WORK(IU) is LDA by N and WORK(IR) is M by M */
02681 
02682                             ldwrku = *lda;
02683                             chunk = *n;
02684                             ldwrkr = *m;
02685                         } else {
02686 
02687 /*                    WORK(IU) is M by CHUNK and WORK(IR) is M by M */
02688 
02689                             ldwrku = *m;
02690                             chunk = (*lwork - *m * *m - *m) / *m;
02691                             ldwrkr = *m;
02692                         }
02693                     }
02694                     itau = ir + ldwrkr * *m;
02695                     iwork = itau + *m;
02696 
02697 /*                 Compute A=L*Q */
02698 /*                 (Workspace: need M*M+2*M, prefer M*M+M+M*NB) */
02699 
02700                     i__3 = *lwork - iwork + 1;
02701                     sgelqf_(m, n, &a[a_offset], lda, &work[itau], &work[iwork]
02702 , &i__3, &ierr);
02703 
02704 /*                 Copy L to U, zeroing about above it */
02705 
02706                     slacpy_("L", m, m, &a[a_offset], lda, &u[u_offset], ldu);
02707                     i__3 = *m - 1;
02708                     i__2 = *m - 1;
02709                     slaset_("U", &i__3, &i__2, &c_b421, &c_b421, &u[(u_dim1 <<
02710                              1) + 1], ldu);
02711 
02712 /*                 Generate Q in A */
02713 /*                 (Workspace: need M*M+2*M, prefer M*M+M+M*NB) */
02714 
02715                     i__3 = *lwork - iwork + 1;
02716                     sorglq_(m, n, m, &a[a_offset], lda, &work[itau], &work[
02717                             iwork], &i__3, &ierr);
02718                     ie = itau;
02719                     itauq = ie + *m;
02720                     itaup = itauq + *m;
02721                     iwork = itaup + *m;
02722 
02723 /*                 Bidiagonalize L in U, copying result to WORK(IR) */
02724 /*                 (Workspace: need M*M+4*M, prefer M*M+3*M+2*M*NB) */
02725 
02726                     i__3 = *lwork - iwork + 1;
02727                     sgebrd_(m, m, &u[u_offset], ldu, &s[1], &work[ie], &work[
02728                             itauq], &work[itaup], &work[iwork], &i__3, &ierr);
02729                     slacpy_("U", m, m, &u[u_offset], ldu, &work[ir], &ldwrkr);
02730 
02731 /*                 Generate right vectors bidiagonalizing L in WORK(IR) */
02732 /*                 (Workspace: need M*M+4*M-1, prefer M*M+3*M+(M-1)*NB) */
02733 
02734                     i__3 = *lwork - iwork + 1;
02735                     sorgbr_("P", m, m, m, &work[ir], &ldwrkr, &work[itaup], &
02736                             work[iwork], &i__3, &ierr);
02737 
02738 /*                 Generate left vectors bidiagonalizing L in U */
02739 /*                 (Workspace: need M*M+4*M, prefer M*M+3*M+M*NB) */
02740 
02741                     i__3 = *lwork - iwork + 1;
02742                     sorgbr_("Q", m, m, m, &u[u_offset], ldu, &work[itauq], &
02743                             work[iwork], &i__3, &ierr);
02744                     iwork = ie + *m;
02745 
02746 /*                 Perform bidiagonal QR iteration, computing left */
02747 /*                 singular vectors of L in U, and computing right */
02748 /*                 singular vectors of L in WORK(IR) */
02749 /*                 (Workspace: need M*M+BDSPAC) */
02750 
02751                     sbdsqr_("U", m, m, m, &c__0, &s[1], &work[ie], &work[ir], 
02752                             &ldwrkr, &u[u_offset], ldu, dum, &c__1, &work[
02753                             iwork], info);
02754                     iu = ie + *m;
02755 
02756 /*                 Multiply right singular vectors of L in WORK(IR) by Q */
02757 /*                 in A, storing result in WORK(IU) and copying to A */
02758 /*                 (Workspace: need M*M+2*M, prefer M*M+M*N+M)) */
02759 
02760                     i__3 = *n;
02761                     i__2 = chunk;
02762                     for (i__ = 1; i__2 < 0 ? i__ >= i__3 : i__ <= i__3; i__ +=
02763                              i__2) {
02764 /* Computing MIN */
02765                         i__4 = *n - i__ + 1;
02766                         blk = min(i__4,chunk);
02767                         sgemm_("N", "N", m, &blk, m, &c_b443, &work[ir], &
02768                                 ldwrkr, &a[i__ * a_dim1 + 1], lda, &c_b421, &
02769                                 work[iu], &ldwrku);
02770                         slacpy_("F", m, &blk, &work[iu], &ldwrku, &a[i__ * 
02771                                 a_dim1 + 1], lda);
02772 /* L40: */
02773                     }
02774 
02775                 } else {
02776 
02777 /*                 Insufficient workspace for a fast algorithm */
02778 
02779                     itau = 1;
02780                     iwork = itau + *m;
02781 
02782 /*                 Compute A=L*Q */
02783 /*                 (Workspace: need 2*M, prefer M+M*NB) */
02784 
02785                     i__2 = *lwork - iwork + 1;
02786                     sgelqf_(m, n, &a[a_offset], lda, &work[itau], &work[iwork]
02787 , &i__2, &ierr);
02788 
02789 /*                 Copy L to U, zeroing out above it */
02790 
02791                     slacpy_("L", m, m, &a[a_offset], lda, &u[u_offset], ldu);
02792                     i__2 = *m - 1;
02793                     i__3 = *m - 1;
02794                     slaset_("U", &i__2, &i__3, &c_b421, &c_b421, &u[(u_dim1 <<
02795                              1) + 1], ldu);
02796 
02797 /*                 Generate Q in A */
02798 /*                 (Workspace: need 2*M, prefer M+M*NB) */
02799 
02800                     i__2 = *lwork - iwork + 1;
02801                     sorglq_(m, n, m, &a[a_offset], lda, &work[itau], &work[
02802                             iwork], &i__2, &ierr);
02803                     ie = itau;
02804                     itauq = ie + *m;
02805                     itaup = itauq + *m;
02806                     iwork = itaup + *m;
02807 
02808 /*                 Bidiagonalize L in U */
02809 /*                 (Workspace: need 4*M, prefer 3*M+2*M*NB) */
02810 
02811                     i__2 = *lwork - iwork + 1;
02812                     sgebrd_(m, m, &u[u_offset], ldu, &s[1], &work[ie], &work[
02813                             itauq], &work[itaup], &work[iwork], &i__2, &ierr);
02814 
02815 /*                 Multiply right vectors bidiagonalizing L by Q in A */
02816 /*                 (Workspace: need 3*M+N, prefer 3*M+N*NB) */
02817 
02818                     i__2 = *lwork - iwork + 1;
02819                     sormbr_("P", "L", "T", m, n, m, &u[u_offset], ldu, &work[
02820                             itaup], &a[a_offset], lda, &work[iwork], &i__2, &
02821                             ierr);
02822 
02823 /*                 Generate left vectors bidiagonalizing L in U */
02824 /*                 (Workspace: need 4*M, prefer 3*M+M*NB) */
02825 
02826                     i__2 = *lwork - iwork + 1;
02827                     sorgbr_("Q", m, m, m, &u[u_offset], ldu, &work[itauq], &
02828                             work[iwork], &i__2, &ierr);
02829                     iwork = ie + *m;
02830 
02831 /*                 Perform bidiagonal QR iteration, computing left */
02832 /*                 singular vectors of A in U and computing right */
02833 /*                 singular vectors of A in A */
02834 /*                 (Workspace: need BDSPAC) */
02835 
02836                     sbdsqr_("U", m, n, m, &c__0, &s[1], &work[ie], &a[
02837                             a_offset], lda, &u[u_offset], ldu, dum, &c__1, &
02838                             work[iwork], info);
02839 
02840                 }
02841 
02842             } else if (wntvs) {
02843 
02844                 if (wntun) {
02845 
02846 /*                 Path 4t(N much larger than M, JOBU='N', JOBVT='S') */
02847 /*                 M right singular vectors to be computed in VT and */
02848 /*                 no left singular vectors to be computed */
02849 
02850 /* Computing MAX */
02851                     i__2 = *m << 2;
02852                     if (*lwork >= *m * *m + max(i__2,bdspac)) {
02853 
02854 /*                    Sufficient workspace for a fast algorithm */
02855 
02856                         ir = 1;
02857                         if (*lwork >= wrkbl + *lda * *m) {
02858 
02859 /*                       WORK(IR) is LDA by M */
02860 
02861                             ldwrkr = *lda;
02862                         } else {
02863 
02864 /*                       WORK(IR) is M by M */
02865 
02866                             ldwrkr = *m;
02867                         }
02868                         itau = ir + ldwrkr * *m;
02869                         iwork = itau + *m;
02870 
02871 /*                    Compute A=L*Q */
02872 /*                    (Workspace: need M*M+2*M, prefer M*M+M+M*NB) */
02873 
02874                         i__2 = *lwork - iwork + 1;
02875                         sgelqf_(m, n, &a[a_offset], lda, &work[itau], &work[
02876                                 iwork], &i__2, &ierr);
02877 
02878 /*                    Copy L to WORK(IR), zeroing out above it */
02879 
02880                         slacpy_("L", m, m, &a[a_offset], lda, &work[ir], &
02881                                 ldwrkr);
02882                         i__2 = *m - 1;
02883                         i__3 = *m - 1;
02884                         slaset_("U", &i__2, &i__3, &c_b421, &c_b421, &work[ir 
02885                                 + ldwrkr], &ldwrkr);
02886 
02887 /*                    Generate Q in A */
02888 /*                    (Workspace: need M*M+2*M, prefer M*M+M+M*NB) */
02889 
02890                         i__2 = *lwork - iwork + 1;
02891                         sorglq_(m, n, m, &a[a_offset], lda, &work[itau], &
02892                                 work[iwork], &i__2, &ierr);
02893                         ie = itau;
02894                         itauq = ie + *m;
02895                         itaup = itauq + *m;
02896                         iwork = itaup + *m;
02897 
02898 /*                    Bidiagonalize L in WORK(IR) */
02899 /*                    (Workspace: need M*M+4*M, prefer M*M+3*M+2*M*NB) */
02900 
02901                         i__2 = *lwork - iwork + 1;
02902                         sgebrd_(m, m, &work[ir], &ldwrkr, &s[1], &work[ie], &
02903                                 work[itauq], &work[itaup], &work[iwork], &
02904                                 i__2, &ierr);
02905 
02906 /*                    Generate right vectors bidiagonalizing L in */
02907 /*                    WORK(IR) */
02908 /*                    (Workspace: need M*M+4*M, prefer M*M+3*M+(M-1)*NB) */
02909 
02910                         i__2 = *lwork - iwork + 1;
02911                         sorgbr_("P", m, m, m, &work[ir], &ldwrkr, &work[itaup]
02912 , &work[iwork], &i__2, &ierr);
02913                         iwork = ie + *m;
02914 
02915 /*                    Perform bidiagonal QR iteration, computing right */
02916 /*                    singular vectors of L in WORK(IR) */
02917 /*                    (Workspace: need M*M+BDSPAC) */
02918 
02919                         sbdsqr_("U", m, m, &c__0, &c__0, &s[1], &work[ie], &
02920                                 work[ir], &ldwrkr, dum, &c__1, dum, &c__1, &
02921                                 work[iwork], info);
02922 
02923 /*                    Multiply right singular vectors of L in WORK(IR) by */
02924 /*                    Q in A, storing result in VT */
02925 /*                    (Workspace: need M*M) */
02926 
02927                         sgemm_("N", "N", m, n, m, &c_b443, &work[ir], &ldwrkr, 
02928                                  &a[a_offset], lda, &c_b421, &vt[vt_offset], 
02929                                 ldvt);
02930 
02931                     } else {
02932 
02933 /*                    Insufficient workspace for a fast algorithm */
02934 
02935                         itau = 1;
02936                         iwork = itau + *m;
02937 
02938 /*                    Compute A=L*Q */
02939 /*                    (Workspace: need 2*M, prefer M+M*NB) */
02940 
02941                         i__2 = *lwork - iwork + 1;
02942                         sgelqf_(m, n, &a[a_offset], lda, &work[itau], &work[
02943                                 iwork], &i__2, &ierr);
02944 
02945 /*                    Copy result to VT */
02946 
02947                         slacpy_("U", m, n, &a[a_offset], lda, &vt[vt_offset], 
02948                                 ldvt);
02949 
02950 /*                    Generate Q in VT */
02951 /*                    (Workspace: need 2*M, prefer M+M*NB) */
02952 
02953                         i__2 = *lwork - iwork + 1;
02954                         sorglq_(m, n, m, &vt[vt_offset], ldvt, &work[itau], &
02955                                 work[iwork], &i__2, &ierr);
02956                         ie = itau;
02957                         itauq = ie + *m;
02958                         itaup = itauq + *m;
02959                         iwork = itaup + *m;
02960 
02961 /*                    Zero out above L in A */
02962 
02963                         i__2 = *m - 1;
02964                         i__3 = *m - 1;
02965                         slaset_("U", &i__2, &i__3, &c_b421, &c_b421, &a[(
02966                                 a_dim1 << 1) + 1], lda);
02967 
02968 /*                    Bidiagonalize L in A */
02969 /*                    (Workspace: need 4*M, prefer 3*M+2*M*NB) */
02970 
02971                         i__2 = *lwork - iwork + 1;
02972                         sgebrd_(m, m, &a[a_offset], lda, &s[1], &work[ie], &
02973                                 work[itauq], &work[itaup], &work[iwork], &
02974                                 i__2, &ierr);
02975 
02976 /*                    Multiply right vectors bidiagonalizing L by Q in VT */
02977 /*                    (Workspace: need 3*M+N, prefer 3*M+N*NB) */
02978 
02979                         i__2 = *lwork - iwork + 1;
02980                         sormbr_("P", "L", "T", m, n, m, &a[a_offset], lda, &
02981                                 work[itaup], &vt[vt_offset], ldvt, &work[
02982                                 iwork], &i__2, &ierr);
02983                         iwork = ie + *m;
02984 
02985 /*                    Perform bidiagonal QR iteration, computing right */
02986 /*                    singular vectors of A in VT */
02987 /*                    (Workspace: need BDSPAC) */
02988 
02989                         sbdsqr_("U", m, n, &c__0, &c__0, &s[1], &work[ie], &
02990                                 vt[vt_offset], ldvt, dum, &c__1, dum, &c__1, &
02991                                 work[iwork], info);
02992 
02993                     }
02994 
02995                 } else if (wntuo) {
02996 
02997 /*                 Path 5t(N much larger than M, JOBU='O', JOBVT='S') */
02998 /*                 M right singular vectors to be computed in VT and */
02999 /*                 M left singular vectors to be overwritten on A */
03000 
03001 /* Computing MAX */
03002                     i__2 = *m << 2;
03003                     if (*lwork >= (*m << 1) * *m + max(i__2,bdspac)) {
03004 
03005 /*                    Sufficient workspace for a fast algorithm */
03006 
03007                         iu = 1;
03008                         if (*lwork >= wrkbl + (*lda << 1) * *m) {
03009 
03010 /*                       WORK(IU) is LDA by M and WORK(IR) is LDA by M */
03011 
03012                             ldwrku = *lda;
03013                             ir = iu + ldwrku * *m;
03014                             ldwrkr = *lda;
03015                         } else if (*lwork >= wrkbl + (*lda + *m) * *m) {
03016 
03017 /*                       WORK(IU) is LDA by M and WORK(IR) is M by M */
03018 
03019                             ldwrku = *lda;
03020                             ir = iu + ldwrku * *m;
03021                             ldwrkr = *m;
03022                         } else {
03023 
03024 /*                       WORK(IU) is M by M and WORK(IR) is M by M */
03025 
03026                             ldwrku = *m;
03027                             ir = iu + ldwrku * *m;
03028                             ldwrkr = *m;
03029                         }
03030                         itau = ir + ldwrkr * *m;
03031                         iwork = itau + *m;
03032 
03033 /*                    Compute A=L*Q */
03034 /*                    (Workspace: need 2*M*M+2*M, prefer 2*M*M+M+M*NB) */
03035 
03036                         i__2 = *lwork - iwork + 1;
03037                         sgelqf_(m, n, &a[a_offset], lda, &work[itau], &work[
03038                                 iwork], &i__2, &ierr);
03039 
03040 /*                    Copy L to WORK(IU), zeroing out below it */
03041 
03042                         slacpy_("L", m, m, &a[a_offset], lda, &work[iu], &
03043                                 ldwrku);
03044                         i__2 = *m - 1;
03045                         i__3 = *m - 1;
03046                         slaset_("U", &i__2, &i__3, &c_b421, &c_b421, &work[iu 
03047                                 + ldwrku], &ldwrku);
03048 
03049 /*                    Generate Q in A */
03050 /*                    (Workspace: need 2*M*M+2*M, prefer 2*M*M+M+M*NB) */
03051 
03052                         i__2 = *lwork - iwork + 1;
03053                         sorglq_(m, n, m, &a[a_offset], lda, &work[itau], &
03054                                 work[iwork], &i__2, &ierr);
03055                         ie = itau;
03056                         itauq = ie + *m;
03057                         itaup = itauq + *m;
03058                         iwork = itaup + *m;
03059 
03060 /*                    Bidiagonalize L in WORK(IU), copying result to */
03061 /*                    WORK(IR) */
03062 /*                    (Workspace: need 2*M*M+4*M, */
03063 /*                                prefer 2*M*M+3*M+2*M*NB) */
03064 
03065                         i__2 = *lwork - iwork + 1;
03066                         sgebrd_(m, m, &work[iu], &ldwrku, &s[1], &work[ie], &
03067                                 work[itauq], &work[itaup], &work[iwork], &
03068                                 i__2, &ierr);
03069                         slacpy_("L", m, m, &work[iu], &ldwrku, &work[ir], &
03070                                 ldwrkr);
03071 
03072 /*                    Generate right bidiagonalizing vectors in WORK(IU) */
03073 /*                    (Workspace: need 2*M*M+4*M-1, */
03074 /*                                prefer 2*M*M+3*M+(M-1)*NB) */
03075 
03076                         i__2 = *lwork - iwork + 1;
03077                         sorgbr_("P", m, m, m, &work[iu], &ldwrku, &work[itaup]
03078 , &work[iwork], &i__2, &ierr);
03079 
03080 /*                    Generate left bidiagonalizing vectors in WORK(IR) */
03081 /*                    (Workspace: need 2*M*M+4*M, prefer 2*M*M+3*M+M*NB) */
03082 
03083                         i__2 = *lwork - iwork + 1;
03084                         sorgbr_("Q", m, m, m, &work[ir], &ldwrkr, &work[itauq]
03085 , &work[iwork], &i__2, &ierr);
03086                         iwork = ie + *m;
03087 
03088 /*                    Perform bidiagonal QR iteration, computing left */
03089 /*                    singular vectors of L in WORK(IR) and computing */
03090 /*                    right singular vectors of L in WORK(IU) */
03091 /*                    (Workspace: need 2*M*M+BDSPAC) */
03092 
03093                         sbdsqr_("U", m, m, m, &c__0, &s[1], &work[ie], &work[
03094                                 iu], &ldwrku, &work[ir], &ldwrkr, dum, &c__1, 
03095                                 &work[iwork], info);
03096 
03097 /*                    Multiply right singular vectors of L in WORK(IU) by */
03098 /*                    Q in A, storing result in VT */
03099 /*                    (Workspace: need M*M) */
03100 
03101                         sgemm_("N", "N", m, n, m, &c_b443, &work[iu], &ldwrku, 
03102                                  &a[a_offset], lda, &c_b421, &vt[vt_offset], 
03103                                 ldvt);
03104 
03105 /*                    Copy left singular vectors of L to A */
03106 /*                    (Workspace: need M*M) */
03107 
03108                         slacpy_("F", m, m, &work[ir], &ldwrkr, &a[a_offset], 
03109                                 lda);
03110 
03111                     } else {
03112 
03113 /*                    Insufficient workspace for a fast algorithm */
03114 
03115                         itau = 1;
03116                         iwork = itau + *m;
03117 
03118 /*                    Compute A=L*Q, copying result to VT */
03119 /*                    (Workspace: need 2*M, prefer M+M*NB) */
03120 
03121                         i__2 = *lwork - iwork + 1;
03122                         sgelqf_(m, n, &a[a_offset], lda, &work[itau], &work[
03123                                 iwork], &i__2, &ierr);
03124                         slacpy_("U", m, n, &a[a_offset], lda, &vt[vt_offset], 
03125                                 ldvt);
03126 
03127 /*                    Generate Q in VT */
03128 /*                    (Workspace: need 2*M, prefer M+M*NB) */
03129 
03130                         i__2 = *lwork - iwork + 1;
03131                         sorglq_(m, n, m, &vt[vt_offset], ldvt, &work[itau], &
03132                                 work[iwork], &i__2, &ierr);
03133                         ie = itau;
03134                         itauq = ie + *m;
03135                         itaup = itauq + *m;
03136                         iwork = itaup + *m;
03137 
03138 /*                    Zero out above L in A */
03139 
03140                         i__2 = *m - 1;
03141                         i__3 = *m - 1;
03142                         slaset_("U", &i__2, &i__3, &c_b421, &c_b421, &a[(
03143                                 a_dim1 << 1) + 1], lda);
03144 
03145 /*                    Bidiagonalize L in A */
03146 /*                    (Workspace: need 4*M, prefer 3*M+2*M*NB) */
03147 
03148                         i__2 = *lwork - iwork + 1;
03149                         sgebrd_(m, m, &a[a_offset], lda, &s[1], &work[ie], &
03150                                 work[itauq], &work[itaup], &work[iwork], &
03151                                 i__2, &ierr);
03152 
03153 /*                    Multiply right vectors bidiagonalizing L by Q in VT */
03154 /*                    (Workspace: need 3*M+N, prefer 3*M+N*NB) */
03155 
03156                         i__2 = *lwork - iwork + 1;
03157                         sormbr_("P", "L", "T", m, n, m, &a[a_offset], lda, &
03158                                 work[itaup], &vt[vt_offset], ldvt, &work[
03159                                 iwork], &i__2, &ierr);
03160 
03161 /*                    Generate left bidiagonalizing vectors of L in A */
03162 /*                    (Workspace: need 4*M, prefer 3*M+M*NB) */
03163 
03164                         i__2 = *lwork - iwork + 1;
03165                         sorgbr_("Q", m, m, m, &a[a_offset], lda, &work[itauq], 
03166                                  &work[iwork], &i__2, &ierr);
03167                         iwork = ie + *m;
03168 
03169 /*                    Perform bidiagonal QR iteration, compute left */
03170 /*                    singular vectors of A in A and compute right */
03171 /*                    singular vectors of A in VT */
03172 /*                    (Workspace: need BDSPAC) */
03173 
03174                         sbdsqr_("U", m, n, m, &c__0, &s[1], &work[ie], &vt[
03175                                 vt_offset], ldvt, &a[a_offset], lda, dum, &
03176                                 c__1, &work[iwork], info);
03177 
03178                     }
03179 
03180                 } else if (wntuas) {
03181 
03182 /*                 Path 6t(N much larger than M, JOBU='S' or 'A', */
03183 /*                         JOBVT='S') */
03184 /*                 M right singular vectors to be computed in VT and */
03185 /*                 M left singular vectors to be computed in U */
03186 
03187 /* Computing MAX */
03188                     i__2 = *m << 2;
03189                     if (*lwork >= *m * *m + max(i__2,bdspac)) {
03190 
03191 /*                    Sufficient workspace for a fast algorithm */
03192 
03193                         iu = 1;
03194                         if (*lwork >= wrkbl + *lda * *m) {
03195 
03196 /*                       WORK(IU) is LDA by N */
03197 
03198                             ldwrku = *lda;
03199                         } else {
03200 
03201 /*                       WORK(IU) is LDA by M */
03202 
03203                             ldwrku = *m;
03204                         }
03205                         itau = iu + ldwrku * *m;
03206                         iwork = itau + *m;
03207 
03208 /*                    Compute A=L*Q */
03209 /*                    (Workspace: need M*M+2*M, prefer M*M+M+M*NB) */
03210 
03211                         i__2 = *lwork - iwork + 1;
03212                         sgelqf_(m, n, &a[a_offset], lda, &work[itau], &work[
03213                                 iwork], &i__2, &ierr);
03214 
03215 /*                    Copy L to WORK(IU), zeroing out above it */
03216 
03217                         slacpy_("L", m, m, &a[a_offset], lda, &work[iu], &
03218                                 ldwrku);
03219                         i__2 = *m - 1;
03220                         i__3 = *m - 1;
03221                         slaset_("U", &i__2, &i__3, &c_b421, &c_b421, &work[iu 
03222                                 + ldwrku], &ldwrku);
03223 
03224 /*                    Generate Q in A */
03225 /*                    (Workspace: need M*M+2*M, prefer M*M+M+M*NB) */
03226 
03227                         i__2 = *lwork - iwork + 1;
03228                         sorglq_(m, n, m, &a[a_offset], lda, &work[itau], &
03229                                 work[iwork], &i__2, &ierr);
03230                         ie = itau;
03231                         itauq = ie + *m;
03232                         itaup = itauq + *m;
03233                         iwork = itaup + *m;
03234 
03235 /*                    Bidiagonalize L in WORK(IU), copying result to U */
03236 /*                    (Workspace: need M*M+4*M, prefer M*M+3*M+2*M*NB) */
03237 
03238                         i__2 = *lwork - iwork + 1;
03239                         sgebrd_(m, m, &work[iu], &ldwrku, &s[1], &work[ie], &
03240                                 work[itauq], &work[itaup], &work[iwork], &
03241                                 i__2, &ierr);
03242                         slacpy_("L", m, m, &work[iu], &ldwrku, &u[u_offset], 
03243                                 ldu);
03244 
03245 /*                    Generate right bidiagonalizing vectors in WORK(IU) */
03246 /*                    (Workspace: need M*M+4*M-1, */
03247 /*                                prefer M*M+3*M+(M-1)*NB) */
03248 
03249                         i__2 = *lwork - iwork + 1;
03250                         sorgbr_("P", m, m, m, &work[iu], &ldwrku, &work[itaup]
03251 , &work[iwork], &i__2, &ierr);
03252 
03253 /*                    Generate left bidiagonalizing vectors in U */
03254 /*                    (Workspace: need M*M+4*M, prefer M*M+3*M+M*NB) */
03255 
03256                         i__2 = *lwork - iwork + 1;
03257                         sorgbr_("Q", m, m, m, &u[u_offset], ldu, &work[itauq], 
03258                                  &work[iwork], &i__2, &ierr);
03259                         iwork = ie + *m;
03260 
03261 /*                    Perform bidiagonal QR iteration, computing left */
03262 /*                    singular vectors of L in U and computing right */
03263 /*                    singular vectors of L in WORK(IU) */
03264 /*                    (Workspace: need M*M+BDSPAC) */
03265 
03266                         sbdsqr_("U", m, m, m, &c__0, &s[1], &work[ie], &work[
03267                                 iu], &ldwrku, &u[u_offset], ldu, dum, &c__1, &
03268                                 work[iwork], info);
03269 
03270 /*                    Multiply right singular vectors of L in WORK(IU) by */
03271 /*                    Q in A, storing result in VT */
03272 /*                    (Workspace: need M*M) */
03273 
03274                         sgemm_("N", "N", m, n, m, &c_b443, &work[iu], &ldwrku, 
03275                                  &a[a_offset], lda, &c_b421, &vt[vt_offset], 
03276                                 ldvt);
03277 
03278                     } else {
03279 
03280 /*                    Insufficient workspace for a fast algorithm */
03281 
03282                         itau = 1;
03283                         iwork = itau + *m;
03284 
03285 /*                    Compute A=L*Q, copying result to VT */
03286 /*                    (Workspace: need 2*M, prefer M+M*NB) */
03287 
03288                         i__2 = *lwork - iwork + 1;
03289                         sgelqf_(m, n, &a[a_offset], lda, &work[itau], &work[
03290                                 iwork], &i__2, &ierr);
03291                         slacpy_("U", m, n, &a[a_offset], lda, &vt[vt_offset], 
03292                                 ldvt);
03293 
03294 /*                    Generate Q in VT */
03295 /*                    (Workspace: need 2*M, prefer M+M*NB) */
03296 
03297                         i__2 = *lwork - iwork + 1;
03298                         sorglq_(m, n, m, &vt[vt_offset], ldvt, &work[itau], &
03299                                 work[iwork], &i__2, &ierr);
03300 
03301 /*                    Copy L to U, zeroing out above it */
03302 
03303                         slacpy_("L", m, m, &a[a_offset], lda, &u[u_offset], 
03304                                 ldu);
03305                         i__2 = *m - 1;
03306                         i__3 = *m - 1;
03307                         slaset_("U", &i__2, &i__3, &c_b421, &c_b421, &u[(
03308                                 u_dim1 << 1) + 1], ldu);
03309                         ie = itau;
03310                         itauq = ie + *m;
03311                         itaup = itauq + *m;
03312                         iwork = itaup + *m;
03313 
03314 /*                    Bidiagonalize L in U */
03315 /*                    (Workspace: need 4*M, prefer 3*M+2*M*NB) */
03316 
03317                         i__2 = *lwork - iwork + 1;
03318                         sgebrd_(m, m, &u[u_offset], ldu, &s[1], &work[ie], &
03319                                 work[itauq], &work[itaup], &work[iwork], &
03320                                 i__2, &ierr);
03321 
03322 /*                    Multiply right bidiagonalizing vectors in U by Q */
03323 /*                    in VT */
03324 /*                    (Workspace: need 3*M+N, prefer 3*M+N*NB) */
03325 
03326                         i__2 = *lwork - iwork + 1;
03327                         sormbr_("P", "L", "T", m, n, m, &u[u_offset], ldu, &
03328                                 work[itaup], &vt[vt_offset], ldvt, &work[
03329                                 iwork], &i__2, &ierr);
03330 
03331 /*                    Generate left bidiagonalizing vectors in U */
03332 /*                    (Workspace: need 4*M, prefer 3*M+M*NB) */
03333 
03334                         i__2 = *lwork - iwork + 1;
03335                         sorgbr_("Q", m, m, m, &u[u_offset], ldu, &work[itauq], 
03336                                  &work[iwork], &i__2, &ierr);
03337                         iwork = ie + *m;
03338 
03339 /*                    Perform bidiagonal QR iteration, computing left */
03340 /*                    singular vectors of A in U and computing right */
03341 /*                    singular vectors of A in VT */
03342 /*                    (Workspace: need BDSPAC) */
03343 
03344                         sbdsqr_("U", m, n, m, &c__0, &s[1], &work[ie], &vt[
03345                                 vt_offset], ldvt, &u[u_offset], ldu, dum, &
03346                                 c__1, &work[iwork], info);
03347 
03348                     }
03349 
03350                 }
03351 
03352             } else if (wntva) {
03353 
03354                 if (wntun) {
03355 
03356 /*                 Path 7t(N much larger than M, JOBU='N', JOBVT='A') */
03357 /*                 N right singular vectors to be computed in VT and */
03358 /*                 no left singular vectors to be computed */
03359 
03360 /* Computing MAX */
03361                     i__2 = *n + *m, i__3 = *m << 2, i__2 = max(i__2,i__3);
03362                     if (*lwork >= *m * *m + max(i__2,bdspac)) {
03363 
03364 /*                    Sufficient workspace for a fast algorithm */
03365 
03366                         ir = 1;
03367                         if (*lwork >= wrkbl + *lda * *m) {
03368 
03369 /*                       WORK(IR) is LDA by M */
03370 
03371                             ldwrkr = *lda;
03372                         } else {
03373 
03374 /*                       WORK(IR) is M by M */
03375 
03376                             ldwrkr = *m;
03377                         }
03378                         itau = ir + ldwrkr * *m;
03379                         iwork = itau + *m;
03380 
03381 /*                    Compute A=L*Q, copying result to VT */
03382 /*                    (Workspace: need M*M+2*M, prefer M*M+M+M*NB) */
03383 
03384                         i__2 = *lwork - iwork + 1;
03385                         sgelqf_(m, n, &a[a_offset], lda, &work[itau], &work[
03386                                 iwork], &i__2, &ierr);
03387                         slacpy_("U", m, n, &a[a_offset], lda, &vt[vt_offset], 
03388                                 ldvt);
03389 
03390 /*                    Copy L to WORK(IR), zeroing out above it */
03391 
03392                         slacpy_("L", m, m, &a[a_offset], lda, &work[ir], &
03393                                 ldwrkr);
03394                         i__2 = *m - 1;
03395                         i__3 = *m - 1;
03396                         slaset_("U", &i__2, &i__3, &c_b421, &c_b421, &work[ir 
03397                                 + ldwrkr], &ldwrkr);
03398 
03399 /*                    Generate Q in VT */
03400 /*                    (Workspace: need M*M+M+N, prefer M*M+M+N*NB) */
03401 
03402                         i__2 = *lwork - iwork + 1;
03403                         sorglq_(n, n, m, &vt[vt_offset], ldvt, &work[itau], &
03404                                 work[iwork], &i__2, &ierr);
03405                         ie = itau;
03406                         itauq = ie + *m;
03407                         itaup = itauq + *m;
03408                         iwork = itaup + *m;
03409 
03410 /*                    Bidiagonalize L in WORK(IR) */
03411 /*                    (Workspace: need M*M+4*M, prefer M*M+3*M+2*M*NB) */
03412 
03413                         i__2 = *lwork - iwork + 1;
03414                         sgebrd_(m, m, &work[ir], &ldwrkr, &s[1], &work[ie], &
03415                                 work[itauq], &work[itaup], &work[iwork], &
03416                                 i__2, &ierr);
03417 
03418 /*                    Generate right bidiagonalizing vectors in WORK(IR) */
03419 /*                    (Workspace: need M*M+4*M-1, */
03420 /*                                prefer M*M+3*M+(M-1)*NB) */
03421 
03422                         i__2 = *lwork - iwork + 1;
03423                         sorgbr_("P", m, m, m, &work[ir], &ldwrkr, &work[itaup]
03424 , &work[iwork], &i__2, &ierr);
03425                         iwork = ie + *m;
03426 
03427 /*                    Perform bidiagonal QR iteration, computing right */
03428 /*                    singular vectors of L in WORK(IR) */
03429 /*                    (Workspace: need M*M+BDSPAC) */
03430 
03431                         sbdsqr_("U", m, m, &c__0, &c__0, &s[1], &work[ie], &
03432                                 work[ir], &ldwrkr, dum, &c__1, dum, &c__1, &
03433                                 work[iwork], info);
03434 
03435 /*                    Multiply right singular vectors of L in WORK(IR) by */
03436 /*                    Q in VT, storing result in A */
03437 /*                    (Workspace: need M*M) */
03438 
03439                         sgemm_("N", "N", m, n, m, &c_b443, &work[ir], &ldwrkr, 
03440                                  &vt[vt_offset], ldvt, &c_b421, &a[a_offset], 
03441                                 lda);
03442 
03443 /*                    Copy right singular vectors of A from A to VT */
03444 
03445                         slacpy_("F", m, n, &a[a_offset], lda, &vt[vt_offset], 
03446                                 ldvt);
03447 
03448                     } else {
03449 
03450 /*                    Insufficient workspace for a fast algorithm */
03451 
03452                         itau = 1;
03453                         iwork = itau + *m;
03454 
03455 /*                    Compute A=L*Q, copying result to VT */
03456 /*                    (Workspace: need 2*M, prefer M+M*NB) */
03457 
03458                         i__2 = *lwork - iwork + 1;
03459                         sgelqf_(m, n, &a[a_offset], lda, &work[itau], &work[
03460                                 iwork], &i__2, &ierr);
03461                         slacpy_("U", m, n, &a[a_offset], lda, &vt[vt_offset], 
03462                                 ldvt);
03463 
03464 /*                    Generate Q in VT */
03465 /*                    (Workspace: need M+N, prefer M+N*NB) */
03466 
03467                         i__2 = *lwork - iwork + 1;
03468                         sorglq_(n, n, m, &vt[vt_offset], ldvt, &work[itau], &
03469                                 work[iwork], &i__2, &ierr);
03470                         ie = itau;
03471                         itauq = ie + *m;
03472                         itaup = itauq + *m;
03473                         iwork = itaup + *m;
03474 
03475 /*                    Zero out above L in A */
03476 
03477                         i__2 = *m - 1;
03478                         i__3 = *m - 1;
03479                         slaset_("U", &i__2, &i__3, &c_b421, &c_b421, &a[(
03480                                 a_dim1 << 1) + 1], lda);
03481 
03482 /*                    Bidiagonalize L in A */
03483 /*                    (Workspace: need 4*M, prefer 3*M+2*M*NB) */
03484 
03485                         i__2 = *lwork - iwork + 1;
03486                         sgebrd_(m, m, &a[a_offset], lda, &s[1], &work[ie], &
03487                                 work[itauq], &work[itaup], &work[iwork], &
03488                                 i__2, &ierr);
03489 
03490 /*                    Multiply right bidiagonalizing vectors in A by Q */
03491 /*                    in VT */
03492 /*                    (Workspace: need 3*M+N, prefer 3*M+N*NB) */
03493 
03494                         i__2 = *lwork - iwork + 1;
03495                         sormbr_("P", "L", "T", m, n, m, &a[a_offset], lda, &
03496                                 work[itaup], &vt[vt_offset], ldvt, &work[
03497                                 iwork], &i__2, &ierr);
03498                         iwork = ie + *m;
03499 
03500 /*                    Perform bidiagonal QR iteration, computing right */
03501 /*                    singular vectors of A in VT */
03502 /*                    (Workspace: need BDSPAC) */
03503 
03504                         sbdsqr_("U", m, n, &c__0, &c__0, &s[1], &work[ie], &
03505                                 vt[vt_offset], ldvt, dum, &c__1, dum, &c__1, &
03506                                 work[iwork], info);
03507 
03508                     }
03509 
03510                 } else if (wntuo) {
03511 
03512 /*                 Path 8t(N much larger than M, JOBU='O', JOBVT='A') */
03513 /*                 N right singular vectors to be computed in VT and */
03514 /*                 M left singular vectors to be overwritten on A */
03515 
03516 /* Computing MAX */
03517                     i__2 = *n + *m, i__3 = *m << 2, i__2 = max(i__2,i__3);
03518                     if (*lwork >= (*m << 1) * *m + max(i__2,bdspac)) {
03519 
03520 /*                    Sufficient workspace for a fast algorithm */
03521 
03522                         iu = 1;
03523                         if (*lwork >= wrkbl + (*lda << 1) * *m) {
03524 
03525 /*                       WORK(IU) is LDA by M and WORK(IR) is LDA by M */
03526 
03527                             ldwrku = *lda;
03528                             ir = iu + ldwrku * *m;
03529                             ldwrkr = *lda;
03530                         } else if (*lwork >= wrkbl + (*lda + *m) * *m) {
03531 
03532 /*                       WORK(IU) is LDA by M and WORK(IR) is M by M */
03533 
03534                             ldwrku = *lda;
03535                             ir = iu + ldwrku * *m;
03536                             ldwrkr = *m;
03537                         } else {
03538 
03539 /*                       WORK(IU) is M by M and WORK(IR) is M by M */
03540 
03541                             ldwrku = *m;
03542                             ir = iu + ldwrku * *m;
03543                             ldwrkr = *m;
03544                         }
03545                         itau = ir + ldwrkr * *m;
03546                         iwork = itau + *m;
03547 
03548 /*                    Compute A=L*Q, copying result to VT */
03549 /*                    (Workspace: need 2*M*M+2*M, prefer 2*M*M+M+M*NB) */
03550 
03551                         i__2 = *lwork - iwork + 1;
03552                         sgelqf_(m, n, &a[a_offset], lda, &work[itau], &work[
03553                                 iwork], &i__2, &ierr);
03554                         slacpy_("U", m, n, &a[a_offset], lda, &vt[vt_offset], 
03555                                 ldvt);
03556 
03557 /*                    Generate Q in VT */
03558 /*                    (Workspace: need 2*M*M+M+N, prefer 2*M*M+M+N*NB) */
03559 
03560                         i__2 = *lwork - iwork + 1;
03561                         sorglq_(n, n, m, &vt[vt_offset], ldvt, &work[itau], &
03562                                 work[iwork], &i__2, &ierr);
03563 
03564 /*                    Copy L to WORK(IU), zeroing out above it */
03565 
03566                         slacpy_("L", m, m, &a[a_offset], lda, &work[iu], &
03567                                 ldwrku);
03568                         i__2 = *m - 1;
03569                         i__3 = *m - 1;
03570                         slaset_("U", &i__2, &i__3, &c_b421, &c_b421, &work[iu 
03571                                 + ldwrku], &ldwrku);
03572                         ie = itau;
03573                         itauq = ie + *m;
03574                         itaup = itauq + *m;
03575                         iwork = itaup + *m;
03576 
03577 /*                    Bidiagonalize L in WORK(IU), copying result to */
03578 /*                    WORK(IR) */
03579 /*                    (Workspace: need 2*M*M+4*M, */
03580 /*                                prefer 2*M*M+3*M+2*M*NB) */
03581 
03582                         i__2 = *lwork - iwork + 1;
03583                         sgebrd_(m, m, &work[iu], &ldwrku, &s[1], &work[ie], &
03584                                 work[itauq], &work[itaup], &work[iwork], &
03585                                 i__2, &ierr);
03586                         slacpy_("L", m, m, &work[iu], &ldwrku, &work[ir], &
03587                                 ldwrkr);
03588 
03589 /*                    Generate right bidiagonalizing vectors in WORK(IU) */
03590 /*                    (Workspace: need 2*M*M+4*M-1, */
03591 /*                                prefer 2*M*M+3*M+(M-1)*NB) */
03592 
03593                         i__2 = *lwork - iwork + 1;
03594                         sorgbr_("P", m, m, m, &work[iu], &ldwrku, &work[itaup]
03595 , &work[iwork], &i__2, &ierr);
03596 
03597 /*                    Generate left bidiagonalizing vectors in WORK(IR) */
03598 /*                    (Workspace: need 2*M*M+4*M, prefer 2*M*M+3*M+M*NB) */
03599 
03600                         i__2 = *lwork - iwork + 1;
03601                         sorgbr_("Q", m, m, m, &work[ir], &ldwrkr, &work[itauq]
03602 , &work[iwork], &i__2, &ierr);
03603                         iwork = ie + *m;
03604 
03605 /*                    Perform bidiagonal QR iteration, computing left */
03606 /*                    singular vectors of L in WORK(IR) and computing */
03607 /*                    right singular vectors of L in WORK(IU) */
03608 /*                    (Workspace: need 2*M*M+BDSPAC) */
03609 
03610                         sbdsqr_("U", m, m, m, &c__0, &s[1], &work[ie], &work[
03611                                 iu], &ldwrku, &work[ir], &ldwrkr, dum, &c__1, 
03612                                 &work[iwork], info);
03613 
03614 /*                    Multiply right singular vectors of L in WORK(IU) by */
03615 /*                    Q in VT, storing result in A */
03616 /*                    (Workspace: need M*M) */
03617 
03618                         sgemm_("N", "N", m, n, m, &c_b443, &work[iu], &ldwrku, 
03619                                  &vt[vt_offset], ldvt, &c_b421, &a[a_offset], 
03620                                 lda);
03621 
03622 /*                    Copy right singular vectors of A from A to VT */
03623 
03624                         slacpy_("F", m, n, &a[a_offset], lda, &vt[vt_offset], 
03625                                 ldvt);
03626 
03627 /*                    Copy left singular vectors of A from WORK(IR) to A */
03628 
03629                         slacpy_("F", m, m, &work[ir], &ldwrkr, &a[a_offset], 
03630                                 lda);
03631 
03632                     } else {
03633 
03634 /*                    Insufficient workspace for a fast algorithm */
03635 
03636                         itau = 1;
03637                         iwork = itau + *m;
03638 
03639 /*                    Compute A=L*Q, copying result to VT */
03640 /*                    (Workspace: need 2*M, prefer M+M*NB) */
03641 
03642                         i__2 = *lwork - iwork + 1;
03643                         sgelqf_(m, n, &a[a_offset], lda, &work[itau], &work[
03644                                 iwork], &i__2, &ierr);
03645                         slacpy_("U", m, n, &a[a_offset], lda, &vt[vt_offset], 
03646                                 ldvt);
03647 
03648 /*                    Generate Q in VT */
03649 /*                    (Workspace: need M+N, prefer M+N*NB) */
03650 
03651                         i__2 = *lwork - iwork + 1;
03652                         sorglq_(n, n, m, &vt[vt_offset], ldvt, &work[itau], &
03653                                 work[iwork], &i__2, &ierr);
03654                         ie = itau;
03655                         itauq = ie + *m;
03656                         itaup = itauq + *m;
03657                         iwork = itaup + *m;
03658 
03659 /*                    Zero out above L in A */
03660 
03661                         i__2 = *m - 1;
03662                         i__3 = *m - 1;
03663                         slaset_("U", &i__2, &i__3, &c_b421, &c_b421, &a[(
03664                                 a_dim1 << 1) + 1], lda);
03665 
03666 /*                    Bidiagonalize L in A */
03667 /*                    (Workspace: need 4*M, prefer 3*M+2*M*NB) */
03668 
03669                         i__2 = *lwork - iwork + 1;
03670                         sgebrd_(m, m, &a[a_offset], lda, &s[1], &work[ie], &
03671                                 work[itauq], &work[itaup], &work[iwork], &
03672                                 i__2, &ierr);
03673 
03674 /*                    Multiply right bidiagonalizing vectors in A by Q */
03675 /*                    in VT */
03676 /*                    (Workspace: need 3*M+N, prefer 3*M+N*NB) */
03677 
03678                         i__2 = *lwork - iwork + 1;
03679                         sormbr_("P", "L", "T", m, n, m, &a[a_offset], lda, &
03680                                 work[itaup], &vt[vt_offset], ldvt, &work[
03681                                 iwork], &i__2, &ierr);
03682 
03683 /*                    Generate left bidiagonalizing vectors in A */
03684 /*                    (Workspace: need 4*M, prefer 3*M+M*NB) */
03685 
03686                         i__2 = *lwork - iwork + 1;
03687                         sorgbr_("Q", m, m, m, &a[a_offset], lda, &work[itauq], 
03688                                  &work[iwork], &i__2, &ierr);
03689                         iwork = ie + *m;
03690 
03691 /*                    Perform bidiagonal QR iteration, computing left */
03692 /*                    singular vectors of A in A and computing right */
03693 /*                    singular vectors of A in VT */
03694 /*                    (Workspace: need BDSPAC) */
03695 
03696                         sbdsqr_("U", m, n, m, &c__0, &s[1], &work[ie], &vt[
03697                                 vt_offset], ldvt, &a[a_offset], lda, dum, &
03698                                 c__1, &work[iwork], info);
03699 
03700                     }
03701 
03702                 } else if (wntuas) {
03703 
03704 /*                 Path 9t(N much larger than M, JOBU='S' or 'A', */
03705 /*                         JOBVT='A') */
03706 /*                 N right singular vectors to be computed in VT and */
03707 /*                 M left singular vectors to be computed in U */
03708 
03709 /* Computing MAX */
03710                     i__2 = *n + *m, i__3 = *m << 2, i__2 = max(i__2,i__3);
03711                     if (*lwork >= *m * *m + max(i__2,bdspac)) {
03712 
03713 /*                    Sufficient workspace for a fast algorithm */
03714 
03715                         iu = 1;
03716                         if (*lwork >= wrkbl + *lda * *m) {
03717 
03718 /*                       WORK(IU) is LDA by M */
03719 
03720                             ldwrku = *lda;
03721                         } else {
03722 
03723 /*                       WORK(IU) is M by M */
03724 
03725                             ldwrku = *m;
03726                         }
03727                         itau = iu + ldwrku * *m;
03728                         iwork = itau + *m;
03729 
03730 /*                    Compute A=L*Q, copying result to VT */
03731 /*                    (Workspace: need M*M+2*M, prefer M*M+M+M*NB) */
03732 
03733                         i__2 = *lwork - iwork + 1;
03734                         sgelqf_(m, n, &a[a_offset], lda, &work[itau], &work[
03735                                 iwork], &i__2, &ierr);
03736                         slacpy_("U", m, n, &a[a_offset], lda, &vt[vt_offset], 
03737                                 ldvt);
03738 
03739 /*                    Generate Q in VT */
03740 /*                    (Workspace: need M*M+M+N, prefer M*M+M+N*NB) */
03741 
03742                         i__2 = *lwork - iwork + 1;
03743                         sorglq_(n, n, m, &vt[vt_offset], ldvt, &work[itau], &
03744                                 work[iwork], &i__2, &ierr);
03745 
03746 /*                    Copy L to WORK(IU), zeroing out above it */
03747 
03748                         slacpy_("L", m, m, &a[a_offset], lda, &work[iu], &
03749                                 ldwrku);
03750                         i__2 = *m - 1;
03751                         i__3 = *m - 1;
03752                         slaset_("U", &i__2, &i__3, &c_b421, &c_b421, &work[iu 
03753                                 + ldwrku], &ldwrku);
03754                         ie = itau;
03755                         itauq = ie + *m;
03756                         itaup = itauq + *m;
03757                         iwork = itaup + *m;
03758 
03759 /*                    Bidiagonalize L in WORK(IU), copying result to U */
03760 /*                    (Workspace: need M*M+4*M, prefer M*M+3*M+2*M*NB) */
03761 
03762                         i__2 = *lwork - iwork + 1;
03763                         sgebrd_(m, m, &work[iu], &ldwrku, &s[1], &work[ie], &
03764                                 work[itauq], &work[itaup], &work[iwork], &
03765                                 i__2, &ierr);
03766                         slacpy_("L", m, m, &work[iu], &ldwrku, &u[u_offset], 
03767                                 ldu);
03768 
03769 /*                    Generate right bidiagonalizing vectors in WORK(IU) */
03770 /*                    (Workspace: need M*M+4*M, prefer M*M+3*M+(M-1)*NB) */
03771 
03772                         i__2 = *lwork - iwork + 1;
03773                         sorgbr_("P", m, m, m, &work[iu], &ldwrku, &work[itaup]
03774 , &work[iwork], &i__2, &ierr);
03775 
03776 /*                    Generate left bidiagonalizing vectors in U */
03777 /*                    (Workspace: need M*M+4*M, prefer M*M+3*M+M*NB) */
03778 
03779                         i__2 = *lwork - iwork + 1;
03780                         sorgbr_("Q", m, m, m, &u[u_offset], ldu, &work[itauq], 
03781                                  &work[iwork], &i__2, &ierr);
03782                         iwork = ie + *m;
03783 
03784 /*                    Perform bidiagonal QR iteration, computing left */
03785 /*                    singular vectors of L in U and computing right */
03786 /*                    singular vectors of L in WORK(IU) */
03787 /*                    (Workspace: need M*M+BDSPAC) */
03788 
03789                         sbdsqr_("U", m, m, m, &c__0, &s[1], &work[ie], &work[
03790                                 iu], &ldwrku, &u[u_offset], ldu, dum, &c__1, &
03791                                 work[iwork], info);
03792 
03793 /*                    Multiply right singular vectors of L in WORK(IU) by */
03794 /*                    Q in VT, storing result in A */
03795 /*                    (Workspace: need M*M) */
03796 
03797                         sgemm_("N", "N", m, n, m, &c_b443, &work[iu], &ldwrku, 
03798                                  &vt[vt_offset], ldvt, &c_b421, &a[a_offset], 
03799                                 lda);
03800 
03801 /*                    Copy right singular vectors of A from A to VT */
03802 
03803                         slacpy_("F", m, n, &a[a_offset], lda, &vt[vt_offset], 
03804                                 ldvt);
03805 
03806                     } else {
03807 
03808 /*                    Insufficient workspace for a fast algorithm */
03809 
03810                         itau = 1;
03811                         iwork = itau + *m;
03812 
03813 /*                    Compute A=L*Q, copying result to VT */
03814 /*                    (Workspace: need 2*M, prefer M+M*NB) */
03815 
03816                         i__2 = *lwork - iwork + 1;
03817                         sgelqf_(m, n, &a[a_offset], lda, &work[itau], &work[
03818                                 iwork], &i__2, &ierr);
03819                         slacpy_("U", m, n, &a[a_offset], lda, &vt[vt_offset], 
03820                                 ldvt);
03821 
03822 /*                    Generate Q in VT */
03823 /*                    (Workspace: need M+N, prefer M+N*NB) */
03824 
03825                         i__2 = *lwork - iwork + 1;
03826                         sorglq_(n, n, m, &vt[vt_offset], ldvt, &work[itau], &
03827                                 work[iwork], &i__2, &ierr);
03828 
03829 /*                    Copy L to U, zeroing out above it */
03830 
03831                         slacpy_("L", m, m, &a[a_offset], lda, &u[u_offset], 
03832                                 ldu);
03833                         i__2 = *m - 1;
03834                         i__3 = *m - 1;
03835                         slaset_("U", &i__2, &i__3, &c_b421, &c_b421, &u[(
03836                                 u_dim1 << 1) + 1], ldu);
03837                         ie = itau;
03838                         itauq = ie + *m;
03839                         itaup = itauq + *m;
03840                         iwork = itaup + *m;
03841 
03842 /*                    Bidiagonalize L in U */
03843 /*                    (Workspace: need 4*M, prefer 3*M+2*M*NB) */
03844 
03845                         i__2 = *lwork - iwork + 1;
03846                         sgebrd_(m, m, &u[u_offset], ldu, &s[1], &work[ie], &
03847                                 work[itauq], &work[itaup], &work[iwork], &
03848                                 i__2, &ierr);
03849 
03850 /*                    Multiply right bidiagonalizing vectors in U by Q */
03851 /*                    in VT */
03852 /*                    (Workspace: need 3*M+N, prefer 3*M+N*NB) */
03853 
03854                         i__2 = *lwork - iwork + 1;
03855                         sormbr_("P", "L", "T", m, n, m, &u[u_offset], ldu, &
03856                                 work[itaup], &vt[vt_offset], ldvt, &work[
03857                                 iwork], &i__2, &ierr);
03858 
03859 /*                    Generate left bidiagonalizing vectors in U */
03860 /*                    (Workspace: need 4*M, prefer 3*M+M*NB) */
03861 
03862                         i__2 = *lwork - iwork + 1;
03863                         sorgbr_("Q", m, m, m, &u[u_offset], ldu, &work[itauq], 
03864                                  &work[iwork], &i__2, &ierr);
03865                         iwork = ie + *m;
03866 
03867 /*                    Perform bidiagonal QR iteration, computing left */
03868 /*                    singular vectors of A in U and computing right */
03869 /*                    singular vectors of A in VT */
03870 /*                    (Workspace: need BDSPAC) */
03871 
03872                         sbdsqr_("U", m, n, m, &c__0, &s[1], &work[ie], &vt[
03873                                 vt_offset], ldvt, &u[u_offset], ldu, dum, &
03874                                 c__1, &work[iwork], info);
03875 
03876                     }
03877 
03878                 }
03879 
03880             }
03881 
03882         } else {
03883 
03884 /*           N .LT. MNTHR */
03885 
03886 /*           Path 10t(N greater than M, but not much larger) */
03887 /*           Reduce to bidiagonal form without LQ decomposition */
03888 
03889             ie = 1;
03890             itauq = ie + *m;
03891             itaup = itauq + *m;
03892             iwork = itaup + *m;
03893 
03894 /*           Bidiagonalize A */
03895 /*           (Workspace: need 3*M+N, prefer 3*M+(M+N)*NB) */
03896 
03897             i__2 = *lwork - iwork + 1;
03898             sgebrd_(m, n, &a[a_offset], lda, &s[1], &work[ie], &work[itauq], &
03899                     work[itaup], &work[iwork], &i__2, &ierr);
03900             if (wntuas) {
03901 
03902 /*              If left singular vectors desired in U, copy result to U */
03903 /*              and generate left bidiagonalizing vectors in U */
03904 /*              (Workspace: need 4*M-1, prefer 3*M+(M-1)*NB) */
03905 
03906                 slacpy_("L", m, m, &a[a_offset], lda, &u[u_offset], ldu);
03907                 i__2 = *lwork - iwork + 1;
03908                 sorgbr_("Q", m, m, n, &u[u_offset], ldu, &work[itauq], &work[
03909                         iwork], &i__2, &ierr);
03910             }
03911             if (wntvas) {
03912 
03913 /*              If right singular vectors desired in VT, copy result to */
03914 /*              VT and generate right bidiagonalizing vectors in VT */
03915 /*              (Workspace: need 3*M+NRVT, prefer 3*M+NRVT*NB) */
03916 
03917                 slacpy_("U", m, n, &a[a_offset], lda, &vt[vt_offset], ldvt);
03918                 if (wntva) {
03919                     nrvt = *n;
03920                 }
03921                 if (wntvs) {
03922                     nrvt = *m;
03923                 }
03924                 i__2 = *lwork - iwork + 1;
03925                 sorgbr_("P", &nrvt, n, m, &vt[vt_offset], ldvt, &work[itaup], 
03926                         &work[iwork], &i__2, &ierr);
03927             }
03928             if (wntuo) {
03929 
03930 /*              If left singular vectors desired in A, generate left */
03931 /*              bidiagonalizing vectors in A */
03932 /*              (Workspace: need 4*M-1, prefer 3*M+(M-1)*NB) */
03933 
03934                 i__2 = *lwork - iwork + 1;
03935                 sorgbr_("Q", m, m, n, &a[a_offset], lda, &work[itauq], &work[
03936                         iwork], &i__2, &ierr);
03937             }
03938             if (wntvo) {
03939 
03940 /*              If right singular vectors desired in A, generate right */
03941 /*              bidiagonalizing vectors in A */
03942 /*              (Workspace: need 4*M, prefer 3*M+M*NB) */
03943 
03944                 i__2 = *lwork - iwork + 1;
03945                 sorgbr_("P", m, n, m, &a[a_offset], lda, &work[itaup], &work[
03946                         iwork], &i__2, &ierr);
03947             }
03948             iwork = ie + *m;
03949             if (wntuas || wntuo) {
03950                 nru = *m;
03951             }
03952             if (wntun) {
03953                 nru = 0;
03954             }
03955             if (wntvas || wntvo) {
03956                 ncvt = *n;
03957             }
03958             if (wntvn) {
03959                 ncvt = 0;
03960             }
03961             if (! wntuo && ! wntvo) {
03962 
03963 /*              Perform bidiagonal QR iteration, if desired, computing */
03964 /*              left singular vectors in U and computing right singular */
03965 /*              vectors in VT */
03966 /*              (Workspace: need BDSPAC) */
03967 
03968                 sbdsqr_("L", m, &ncvt, &nru, &c__0, &s[1], &work[ie], &vt[
03969                         vt_offset], ldvt, &u[u_offset], ldu, dum, &c__1, &
03970                         work[iwork], info);
03971             } else if (! wntuo && wntvo) {
03972 
03973 /*              Perform bidiagonal QR iteration, if desired, computing */
03974 /*              left singular vectors in U and computing right singular */
03975 /*              vectors in A */
03976 /*              (Workspace: need BDSPAC) */
03977 
03978                 sbdsqr_("L", m, &ncvt, &nru, &c__0, &s[1], &work[ie], &a[
03979                         a_offset], lda, &u[u_offset], ldu, dum, &c__1, &work[
03980                         iwork], info);
03981             } else {
03982 
03983 /*              Perform bidiagonal QR iteration, if desired, computing */
03984 /*              left singular vectors in A and computing right singular */
03985 /*              vectors in VT */
03986 /*              (Workspace: need BDSPAC) */
03987 
03988                 sbdsqr_("L", m, &ncvt, &nru, &c__0, &s[1], &work[ie], &vt[
03989                         vt_offset], ldvt, &a[a_offset], lda, dum, &c__1, &
03990                         work[iwork], info);
03991             }
03992 
03993         }
03994 
03995     }
03996 
03997 /*     If SBDSQR failed to converge, copy unconverged superdiagonals */
03998 /*     to WORK( 2:MINMN ) */
03999 
04000     if (*info != 0) {
04001         if (ie > 2) {
04002             i__2 = minmn - 1;
04003             for (i__ = 1; i__ <= i__2; ++i__) {
04004                 work[i__ + 1] = work[i__ + ie - 1];
04005 /* L50: */
04006             }
04007         }
04008         if (ie < 2) {
04009             for (i__ = minmn - 1; i__ >= 1; --i__) {
04010                 work[i__ + 1] = work[i__ + ie - 1];
04011 /* L60: */
04012             }
04013         }
04014     }
04015 
04016 /*     Undo scaling if necessary */
04017 
04018     if (iscl == 1) {
04019         if (anrm > bignum) {
04020             slascl_("G", &c__0, &c__0, &bignum, &anrm, &minmn, &c__1, &s[1], &
04021                     minmn, &ierr);
04022         }
04023         if (*info != 0 && anrm > bignum) {
04024             i__2 = minmn - 1;
04025             slascl_("G", &c__0, &c__0, &bignum, &anrm, &i__2, &c__1, &work[2], 
04026                      &minmn, &ierr);
04027         }
04028         if (anrm < smlnum) {
04029             slascl_("G", &c__0, &c__0, &smlnum, &anrm, &minmn, &c__1, &s[1], &
04030                     minmn, &ierr);
04031         }
04032         if (*info != 0 && anrm < smlnum) {
04033             i__2 = minmn - 1;
04034             slascl_("G", &c__0, &c__0, &smlnum, &anrm, &i__2, &c__1, &work[2], 
04035                      &minmn, &ierr);
04036         }
04037     }
04038 
04039 /*     Return optimal workspace in WORK(1) */
04040 
04041     work[1] = (real) maxwrk;
04042 
04043     return 0;
04044 
04045 /*     End of SGESVD */
04046 
04047 } /* sgesvd_ */