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

swiftnav
Author(s):
autogenerated on Sat Jun 8 2019 18:55:44