cgesdd.c
Go to the documentation of this file.
00001 /* cgesdd.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 complex c_b1 = {0.f,0.f};
00019 static complex c_b2 = {1.f,0.f};
00020 static integer c__1 = 1;
00021 static integer c_n1 = -1;
00022 static integer c__0 = 0;
00023 
00024 /* Subroutine */ int cgesdd_(char *jobz, integer *m, integer *n, complex *a, 
00025         integer *lda, real *s, complex *u, integer *ldu, complex *vt, integer 
00026         *ldvt, complex *work, integer *lwork, real *rwork, integer *iwork, 
00027         integer *info)
00028 {
00029     /* System generated locals */
00030     integer a_dim1, a_offset, u_dim1, u_offset, vt_dim1, vt_offset, i__1, 
00031             i__2, i__3;
00032 
00033     /* Builtin functions */
00034     double sqrt(doublereal);
00035 
00036     /* Local variables */
00037     integer i__, ie, il, ir, iu, blk;
00038     real dum[1], eps;
00039     integer iru, ivt, iscl;
00040     real anrm;
00041     integer idum[1], ierr, itau, irvt;
00042     extern /* Subroutine */ int cgemm_(char *, char *, integer *, integer *, 
00043             integer *, complex *, complex *, integer *, complex *, integer *, 
00044             complex *, complex *, integer *);
00045     extern logical lsame_(char *, char *);
00046     integer chunk, minmn, wrkbl, itaup, itauq;
00047     logical wntqa;
00048     integer nwork;
00049     extern /* Subroutine */ int clacp2_(char *, integer *, integer *, real *, 
00050             integer *, complex *, integer *);
00051     logical wntqn, wntqo, wntqs;
00052     integer mnthr1, mnthr2;
00053     extern /* Subroutine */ int cgebrd_(integer *, integer *, complex *, 
00054             integer *, real *, real *, complex *, complex *, complex *, 
00055             integer *, integer *);
00056     extern doublereal clange_(char *, integer *, integer *, complex *, 
00057             integer *, real *);
00058     extern /* Subroutine */ int cgelqf_(integer *, integer *, complex *, 
00059             integer *, complex *, complex *, integer *, integer *), clacrm_(
00060             integer *, integer *, complex *, integer *, real *, integer *, 
00061             complex *, integer *, real *), clarcm_(integer *, integer *, real 
00062             *, integer *, complex *, integer *, complex *, integer *, real *),
00063              clascl_(char *, integer *, integer *, real *, real *, integer *, 
00064             integer *, complex *, integer *, integer *), sbdsdc_(char 
00065             *, char *, integer *, real *, real *, real *, integer *, real *, 
00066             integer *, real *, integer *, real *, integer *, integer *), cgeqrf_(integer *, integer *, complex *, integer 
00067             *, complex *, complex *, integer *, integer *);
00068     extern doublereal slamch_(char *);
00069     extern /* Subroutine */ int clacpy_(char *, integer *, integer *, complex 
00070             *, integer *, complex *, integer *), claset_(char *, 
00071             integer *, integer *, complex *, complex *, complex *, integer *), xerbla_(char *, integer *);
00072     extern integer ilaenv_(integer *, char *, char *, integer *, integer *, 
00073             integer *, integer *);
00074     extern /* Subroutine */ int cungbr_(char *, integer *, integer *, integer 
00075             *, complex *, integer *, complex *, complex *, integer *, integer 
00076             *);
00077     real bignum;
00078     extern /* Subroutine */ int slascl_(char *, integer *, integer *, real *, 
00079             real *, integer *, integer *, real *, integer *, integer *), cunmbr_(char *, char *, char *, integer *, integer *, 
00080             integer *, complex *, integer *, complex *, complex *, integer *, 
00081             complex *, integer *, integer *), cunglq_(
00082             integer *, integer *, integer *, complex *, integer *, complex *, 
00083             complex *, integer *, integer *);
00084     integer ldwrkl;
00085     extern /* Subroutine */ int cungqr_(integer *, integer *, integer *, 
00086             complex *, integer *, complex *, complex *, integer *, integer *);
00087     integer ldwrkr, minwrk, ldwrku, maxwrk, ldwkvt;
00088     real smlnum;
00089     logical wntqas;
00090     integer nrwork;
00091 
00092 
00093 /*  -- LAPACK driver routine (version 3.2) -- */
00094 /*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */
00095 /*     November 2006 */
00096 /*     8-15-00:  Improve consistency of WS calculations (eca) */
00097 
00098 /*     .. Scalar Arguments .. */
00099 /*     .. */
00100 /*     .. Array Arguments .. */
00101 /*     .. */
00102 
00103 /*  Purpose */
00104 /*  ======= */
00105 
00106 /*  CGESDD computes the singular value decomposition (SVD) of a complex */
00107 /*  M-by-N matrix A, optionally computing the left and/or right singular */
00108 /*  vectors, by using divide-and-conquer method. The SVD is written */
00109 
00110 /*       A = U * SIGMA * conjugate-transpose(V) */
00111 
00112 /*  where SIGMA is an M-by-N matrix which is zero except for its */
00113 /*  min(m,n) diagonal elements, U is an M-by-M unitary matrix, and */
00114 /*  V is an N-by-N unitary matrix.  The diagonal elements of SIGMA */
00115 /*  are the singular values of A; they are real and non-negative, and */
00116 /*  are returned in descending order.  The first min(m,n) columns of */
00117 /*  U and V are the left and right singular vectors of A. */
00118 
00119 /*  Note that the routine returns VT = V**H, not V. */
00120 
00121 /*  The divide and conquer algorithm makes very mild assumptions about */
00122 /*  floating point arithmetic. It will work on machines with a guard */
00123 /*  digit in add/subtract, or on those binary machines without guard */
00124 /*  digits which subtract like the Cray X-MP, Cray Y-MP, Cray C-90, or */
00125 /*  Cray-2. It could conceivably fail on hexadecimal or decimal machines */
00126 /*  without guard digits, but we know of none. */
00127 
00128 /*  Arguments */
00129 /*  ========= */
00130 
00131 /*  JOBZ    (input) CHARACTER*1 */
00132 /*          Specifies options for computing all or part of the matrix U: */
00133 /*          = 'A':  all M columns of U and all N rows of V**H are */
00134 /*                  returned in the arrays U and VT; */
00135 /*          = 'S':  the first min(M,N) columns of U and the first */
00136 /*                  min(M,N) rows of V**H are returned in the arrays U */
00137 /*                  and VT; */
00138 /*          = 'O':  If M >= N, the first N columns of U are overwritten */
00139 /*                  in the array A and all rows of V**H are returned in */
00140 /*                  the array VT; */
00141 /*                  otherwise, all columns of U are returned in the */
00142 /*                  array U and the first M rows of V**H are overwritten */
00143 /*                  in the array A; */
00144 /*          = 'N':  no columns of U or rows of V**H are computed. */
00145 
00146 /*  M       (input) INTEGER */
00147 /*          The number of rows of the input matrix A.  M >= 0. */
00148 
00149 /*  N       (input) INTEGER */
00150 /*          The number of columns of the input matrix A.  N >= 0. */
00151 
00152 /*  A       (input/output) COMPLEX array, dimension (LDA,N) */
00153 /*          On entry, the M-by-N matrix A. */
00154 /*          On exit, */
00155 /*          if JOBZ = 'O',  A is overwritten with the first N columns */
00156 /*                          of U (the left singular vectors, stored */
00157 /*                          columnwise) if M >= N; */
00158 /*                          A is overwritten with the first M rows */
00159 /*                          of V**H (the right singular vectors, stored */
00160 /*                          rowwise) otherwise. */
00161 /*          if JOBZ .ne. 'O', the contents of A are destroyed. */
00162 
00163 /*  LDA     (input) INTEGER */
00164 /*          The leading dimension of the array A.  LDA >= max(1,M). */
00165 
00166 /*  S       (output) REAL array, dimension (min(M,N)) */
00167 /*          The singular values of A, sorted so that S(i) >= S(i+1). */
00168 
00169 /*  U       (output) COMPLEX array, dimension (LDU,UCOL) */
00170 /*          UCOL = M if JOBZ = 'A' or JOBZ = 'O' and M < N; */
00171 /*          UCOL = min(M,N) if JOBZ = 'S'. */
00172 /*          If JOBZ = 'A' or JOBZ = 'O' and M < N, U contains the M-by-M */
00173 /*          unitary matrix U; */
00174 /*          if JOBZ = 'S', U contains the first min(M,N) columns of U */
00175 /*          (the left singular vectors, stored columnwise); */
00176 /*          if JOBZ = 'O' and M >= N, or JOBZ = 'N', U is not referenced. */
00177 
00178 /*  LDU     (input) INTEGER */
00179 /*          The leading dimension of the array U.  LDU >= 1; if */
00180 /*          JOBZ = 'S' or 'A' or JOBZ = 'O' and M < N, LDU >= M. */
00181 
00182 /*  VT      (output) COMPLEX array, dimension (LDVT,N) */
00183 /*          If JOBZ = 'A' or JOBZ = 'O' and M >= N, VT contains the */
00184 /*          N-by-N unitary matrix V**H; */
00185 /*          if JOBZ = 'S', VT contains the first min(M,N) rows of */
00186 /*          V**H (the right singular vectors, stored rowwise); */
00187 /*          if JOBZ = 'O' and M < N, or JOBZ = 'N', VT is not referenced. */
00188 
00189 /*  LDVT    (input) INTEGER */
00190 /*          The leading dimension of the array VT.  LDVT >= 1; if */
00191 /*          JOBZ = 'A' or JOBZ = 'O' and M >= N, LDVT >= N; */
00192 /*          if JOBZ = 'S', LDVT >= min(M,N). */
00193 
00194 /*  WORK    (workspace/output) COMPLEX array, dimension (MAX(1,LWORK)) */
00195 /*          On exit, if INFO = 0, WORK(1) returns the optimal LWORK. */
00196 
00197 /*  LWORK   (input) INTEGER */
00198 /*          The dimension of the array WORK. LWORK >= 1. */
00199 /*          if JOBZ = 'N', LWORK >= 2*min(M,N)+max(M,N). */
00200 /*          if JOBZ = 'O', */
00201 /*                LWORK >= 2*min(M,N)*min(M,N)+2*min(M,N)+max(M,N). */
00202 /*          if JOBZ = 'S' or 'A', */
00203 /*                LWORK >= min(M,N)*min(M,N)+2*min(M,N)+max(M,N). */
00204 /*          For good performance, LWORK should generally be larger. */
00205 
00206 /*          If LWORK = -1, a workspace query is assumed.  The optimal */
00207 /*          size for the WORK array is calculated and stored in WORK(1), */
00208 /*          and no other work except argument checking is performed. */
00209 
00210 /*  RWORK   (workspace) REAL array, dimension (MAX(1,LRWORK)) */
00211 /*          If JOBZ = 'N', LRWORK >= 5*min(M,N). */
00212 /*          Otherwise, LRWORK >= 5*min(M,N)*min(M,N) + 7*min(M,N) */
00213 
00214 /*  IWORK   (workspace) INTEGER array, dimension (8*min(M,N)) */
00215 
00216 /*  INFO    (output) INTEGER */
00217 /*          = 0:  successful exit. */
00218 /*          < 0:  if INFO = -i, the i-th argument had an illegal value. */
00219 /*          > 0:  The updating process of SBDSDC did not converge. */
00220 
00221 /*  Further Details */
00222 /*  =============== */
00223 
00224 /*  Based on contributions by */
00225 /*     Ming Gu and Huan Ren, Computer Science Division, University of */
00226 /*     California at Berkeley, USA */
00227 
00228 /*  ===================================================================== */
00229 
00230 /*     .. Parameters .. */
00231 /*     .. */
00232 /*     .. Local Scalars .. */
00233 /*     .. */
00234 /*     .. Local Arrays .. */
00235 /*     .. */
00236 /*     .. External Subroutines .. */
00237 /*     .. */
00238 /*     .. External Functions .. */
00239 /*     .. */
00240 /*     .. Intrinsic Functions .. */
00241 /*     .. */
00242 /*     .. Executable Statements .. */
00243 
00244 /*     Test the input arguments */
00245 
00246     /* Parameter adjustments */
00247     a_dim1 = *lda;
00248     a_offset = 1 + a_dim1;
00249     a -= a_offset;
00250     --s;
00251     u_dim1 = *ldu;
00252     u_offset = 1 + u_dim1;
00253     u -= u_offset;
00254     vt_dim1 = *ldvt;
00255     vt_offset = 1 + vt_dim1;
00256     vt -= vt_offset;
00257     --work;
00258     --rwork;
00259     --iwork;
00260 
00261     /* Function Body */
00262     *info = 0;
00263     minmn = min(*m,*n);
00264     mnthr1 = (integer) (minmn * 17.f / 9.f);
00265     mnthr2 = (integer) (minmn * 5.f / 3.f);
00266     wntqa = lsame_(jobz, "A");
00267     wntqs = lsame_(jobz, "S");
00268     wntqas = wntqa || wntqs;
00269     wntqo = lsame_(jobz, "O");
00270     wntqn = lsame_(jobz, "N");
00271     minwrk = 1;
00272     maxwrk = 1;
00273 
00274     if (! (wntqa || wntqs || wntqo || wntqn)) {
00275         *info = -1;
00276     } else if (*m < 0) {
00277         *info = -2;
00278     } else if (*n < 0) {
00279         *info = -3;
00280     } else if (*lda < max(1,*m)) {
00281         *info = -5;
00282     } else if (*ldu < 1 || wntqas && *ldu < *m || wntqo && *m < *n && *ldu < *
00283             m) {
00284         *info = -8;
00285     } else if (*ldvt < 1 || wntqa && *ldvt < *n || wntqs && *ldvt < minmn || 
00286             wntqo && *m >= *n && *ldvt < *n) {
00287         *info = -10;
00288     }
00289 
00290 /*     Compute workspace */
00291 /*      (Note: Comments in the code beginning "Workspace:" describe the */
00292 /*       minimal amount of workspace needed at that point in the code, */
00293 /*       as well as the preferred amount for good performance. */
00294 /*       CWorkspace refers to complex workspace, and RWorkspace to */
00295 /*       real workspace. NB refers to the optimal block size for the */
00296 /*       immediately following subroutine, as returned by ILAENV.) */
00297 
00298     if (*info == 0 && *m > 0 && *n > 0) {
00299         if (*m >= *n) {
00300 
00301 /*           There is no complex work space needed for bidiagonal SVD */
00302 /*           The real work space needed for bidiagonal SVD is BDSPAC */
00303 /*           for computing singular values and singular vectors; BDSPAN */
00304 /*           for computing singular values only. */
00305 /*           BDSPAC = 5*N*N + 7*N */
00306 /*           BDSPAN = MAX(7*N+4, 3*N+2+SMLSIZ*(SMLSIZ+8)) */
00307 
00308             if (*m >= mnthr1) {
00309                 if (wntqn) {
00310 
00311 /*                 Path 1 (M much larger than N, JOBZ='N') */
00312 
00313                     maxwrk = *n + *n * ilaenv_(&c__1, "CGEQRF", " ", m, n, &
00314                             c_n1, &c_n1);
00315 /* Computing MAX */
00316                     i__1 = maxwrk, i__2 = (*n << 1) + (*n << 1) * ilaenv_(&
00317                             c__1, "CGEBRD", " ", n, n, &c_n1, &c_n1);
00318                     maxwrk = max(i__1,i__2);
00319                     minwrk = *n * 3;
00320                 } else if (wntqo) {
00321 
00322 /*                 Path 2 (M much larger than N, JOBZ='O') */
00323 
00324                     wrkbl = *n + *n * ilaenv_(&c__1, "CGEQRF", " ", m, n, &
00325                             c_n1, &c_n1);
00326 /* Computing MAX */
00327                     i__1 = wrkbl, i__2 = *n + *n * ilaenv_(&c__1, "CUNGQR", 
00328                             " ", m, n, n, &c_n1);
00329                     wrkbl = max(i__1,i__2);
00330 /* Computing MAX */
00331                     i__1 = wrkbl, i__2 = (*n << 1) + (*n << 1) * ilaenv_(&
00332                             c__1, "CGEBRD", " ", n, n, &c_n1, &c_n1);
00333                     wrkbl = max(i__1,i__2);
00334 /* Computing MAX */
00335                     i__1 = wrkbl, i__2 = (*n << 1) + *n * ilaenv_(&c__1, 
00336                             "CUNMBR", "QLN", n, n, n, &c_n1);
00337                     wrkbl = max(i__1,i__2);
00338 /* Computing MAX */
00339                     i__1 = wrkbl, i__2 = (*n << 1) + *n * ilaenv_(&c__1, 
00340                             "CUNMBR", "PRC", n, n, n, &c_n1);
00341                     wrkbl = max(i__1,i__2);
00342                     maxwrk = *m * *n + *n * *n + wrkbl;
00343                     minwrk = (*n << 1) * *n + *n * 3;
00344                 } else if (wntqs) {
00345 
00346 /*                 Path 3 (M much larger than N, JOBZ='S') */
00347 
00348                     wrkbl = *n + *n * ilaenv_(&c__1, "CGEQRF", " ", m, n, &
00349                             c_n1, &c_n1);
00350 /* Computing MAX */
00351                     i__1 = wrkbl, i__2 = *n + *n * ilaenv_(&c__1, "CUNGQR", 
00352                             " ", m, n, n, &c_n1);
00353                     wrkbl = max(i__1,i__2);
00354 /* Computing MAX */
00355                     i__1 = wrkbl, i__2 = (*n << 1) + (*n << 1) * ilaenv_(&
00356                             c__1, "CGEBRD", " ", n, n, &c_n1, &c_n1);
00357                     wrkbl = max(i__1,i__2);
00358 /* Computing MAX */
00359                     i__1 = wrkbl, i__2 = (*n << 1) + *n * ilaenv_(&c__1, 
00360                             "CUNMBR", "QLN", n, n, n, &c_n1);
00361                     wrkbl = max(i__1,i__2);
00362 /* Computing MAX */
00363                     i__1 = wrkbl, i__2 = (*n << 1) + *n * ilaenv_(&c__1, 
00364                             "CUNMBR", "PRC", n, n, n, &c_n1);
00365                     wrkbl = max(i__1,i__2);
00366                     maxwrk = *n * *n + wrkbl;
00367                     minwrk = *n * *n + *n * 3;
00368                 } else if (wntqa) {
00369 
00370 /*                 Path 4 (M much larger than N, JOBZ='A') */
00371 
00372                     wrkbl = *n + *n * ilaenv_(&c__1, "CGEQRF", " ", m, n, &
00373                             c_n1, &c_n1);
00374 /* Computing MAX */
00375                     i__1 = wrkbl, i__2 = *n + *m * ilaenv_(&c__1, "CUNGQR", 
00376                             " ", m, m, n, &c_n1);
00377                     wrkbl = max(i__1,i__2);
00378 /* Computing MAX */
00379                     i__1 = wrkbl, i__2 = (*n << 1) + (*n << 1) * ilaenv_(&
00380                             c__1, "CGEBRD", " ", n, n, &c_n1, &c_n1);
00381                     wrkbl = max(i__1,i__2);
00382 /* Computing MAX */
00383                     i__1 = wrkbl, i__2 = (*n << 1) + *n * ilaenv_(&c__1, 
00384                             "CUNMBR", "QLN", n, n, n, &c_n1);
00385                     wrkbl = max(i__1,i__2);
00386 /* Computing MAX */
00387                     i__1 = wrkbl, i__2 = (*n << 1) + *n * ilaenv_(&c__1, 
00388                             "CUNMBR", "PRC", n, n, n, &c_n1);
00389                     wrkbl = max(i__1,i__2);
00390                     maxwrk = *n * *n + wrkbl;
00391                     minwrk = *n * *n + (*n << 1) + *m;
00392                 }
00393             } else if (*m >= mnthr2) {
00394 
00395 /*              Path 5 (M much larger than N, but not as much as MNTHR1) */
00396 
00397                 maxwrk = (*n << 1) + (*m + *n) * ilaenv_(&c__1, "CGEBRD", 
00398                         " ", m, n, &c_n1, &c_n1);
00399                 minwrk = (*n << 1) + *m;
00400                 if (wntqo) {
00401 /* Computing MAX */
00402                     i__1 = maxwrk, i__2 = (*n << 1) + *n * ilaenv_(&c__1, 
00403                             "CUNGBR", "P", n, n, n, &c_n1);
00404                     maxwrk = max(i__1,i__2);
00405 /* Computing MAX */
00406                     i__1 = maxwrk, i__2 = (*n << 1) + *n * ilaenv_(&c__1, 
00407                             "CUNGBR", "Q", m, n, n, &c_n1);
00408                     maxwrk = max(i__1,i__2);
00409                     maxwrk += *m * *n;
00410                     minwrk += *n * *n;
00411                 } else if (wntqs) {
00412 /* Computing MAX */
00413                     i__1 = maxwrk, i__2 = (*n << 1) + *n * ilaenv_(&c__1, 
00414                             "CUNGBR", "P", n, n, n, &c_n1);
00415                     maxwrk = max(i__1,i__2);
00416 /* Computing MAX */
00417                     i__1 = maxwrk, i__2 = (*n << 1) + *n * ilaenv_(&c__1, 
00418                             "CUNGBR", "Q", m, n, n, &c_n1);
00419                     maxwrk = max(i__1,i__2);
00420                 } else if (wntqa) {
00421 /* Computing MAX */
00422                     i__1 = maxwrk, i__2 = (*n << 1) + *n * ilaenv_(&c__1, 
00423                             "CUNGBR", "P", n, n, n, &c_n1);
00424                     maxwrk = max(i__1,i__2);
00425 /* Computing MAX */
00426                     i__1 = maxwrk, i__2 = (*n << 1) + *m * ilaenv_(&c__1, 
00427                             "CUNGBR", "Q", m, m, n, &c_n1);
00428                     maxwrk = max(i__1,i__2);
00429                 }
00430             } else {
00431 
00432 /*              Path 6 (M at least N, but not much larger) */
00433 
00434                 maxwrk = (*n << 1) + (*m + *n) * ilaenv_(&c__1, "CGEBRD", 
00435                         " ", m, n, &c_n1, &c_n1);
00436                 minwrk = (*n << 1) + *m;
00437                 if (wntqo) {
00438 /* Computing MAX */
00439                     i__1 = maxwrk, i__2 = (*n << 1) + *n * ilaenv_(&c__1, 
00440                             "CUNMBR", "PRC", n, n, n, &c_n1);
00441                     maxwrk = max(i__1,i__2);
00442 /* Computing MAX */
00443                     i__1 = maxwrk, i__2 = (*n << 1) + *n * ilaenv_(&c__1, 
00444                             "CUNMBR", "QLN", m, n, n, &c_n1);
00445                     maxwrk = max(i__1,i__2);
00446                     maxwrk += *m * *n;
00447                     minwrk += *n * *n;
00448                 } else if (wntqs) {
00449 /* Computing MAX */
00450                     i__1 = maxwrk, i__2 = (*n << 1) + *n * ilaenv_(&c__1, 
00451                             "CUNMBR", "PRC", n, n, n, &c_n1);
00452                     maxwrk = max(i__1,i__2);
00453 /* Computing MAX */
00454                     i__1 = maxwrk, i__2 = (*n << 1) + *n * ilaenv_(&c__1, 
00455                             "CUNMBR", "QLN", m, n, n, &c_n1);
00456                     maxwrk = max(i__1,i__2);
00457                 } else if (wntqa) {
00458 /* Computing MAX */
00459                     i__1 = maxwrk, i__2 = (*n << 1) + *n * ilaenv_(&c__1, 
00460                             "CUNGBR", "PRC", n, n, n, &c_n1);
00461                     maxwrk = max(i__1,i__2);
00462 /* Computing MAX */
00463                     i__1 = maxwrk, i__2 = (*n << 1) + *m * ilaenv_(&c__1, 
00464                             "CUNGBR", "QLN", m, m, n, &c_n1);
00465                     maxwrk = max(i__1,i__2);
00466                 }
00467             }
00468         } else {
00469 
00470 /*           There is no complex work space needed for bidiagonal SVD */
00471 /*           The real work space needed for bidiagonal SVD is BDSPAC */
00472 /*           for computing singular values and singular vectors; BDSPAN */
00473 /*           for computing singular values only. */
00474 /*           BDSPAC = 5*M*M + 7*M */
00475 /*           BDSPAN = MAX(7*M+4, 3*M+2+SMLSIZ*(SMLSIZ+8)) */
00476 
00477             if (*n >= mnthr1) {
00478                 if (wntqn) {
00479 
00480 /*                 Path 1t (N much larger than M, JOBZ='N') */
00481 
00482                     maxwrk = *m + *m * ilaenv_(&c__1, "CGELQF", " ", m, n, &
00483                             c_n1, &c_n1);
00484 /* Computing MAX */
00485                     i__1 = maxwrk, i__2 = (*m << 1) + (*m << 1) * ilaenv_(&
00486                             c__1, "CGEBRD", " ", m, m, &c_n1, &c_n1);
00487                     maxwrk = max(i__1,i__2);
00488                     minwrk = *m * 3;
00489                 } else if (wntqo) {
00490 
00491 /*                 Path 2t (N much larger than M, JOBZ='O') */
00492 
00493                     wrkbl = *m + *m * ilaenv_(&c__1, "CGELQF", " ", m, n, &
00494                             c_n1, &c_n1);
00495 /* Computing MAX */
00496                     i__1 = wrkbl, i__2 = *m + *m * ilaenv_(&c__1, "CUNGLQ", 
00497                             " ", m, n, m, &c_n1);
00498                     wrkbl = max(i__1,i__2);
00499 /* Computing MAX */
00500                     i__1 = wrkbl, i__2 = (*m << 1) + (*m << 1) * ilaenv_(&
00501                             c__1, "CGEBRD", " ", m, m, &c_n1, &c_n1);
00502                     wrkbl = max(i__1,i__2);
00503 /* Computing MAX */
00504                     i__1 = wrkbl, i__2 = (*m << 1) + *m * ilaenv_(&c__1, 
00505                             "CUNMBR", "PRC", m, m, m, &c_n1);
00506                     wrkbl = max(i__1,i__2);
00507 /* Computing MAX */
00508                     i__1 = wrkbl, i__2 = (*m << 1) + *m * ilaenv_(&c__1, 
00509                             "CUNMBR", "QLN", m, m, m, &c_n1);
00510                     wrkbl = max(i__1,i__2);
00511                     maxwrk = *m * *n + *m * *m + wrkbl;
00512                     minwrk = (*m << 1) * *m + *m * 3;
00513                 } else if (wntqs) {
00514 
00515 /*                 Path 3t (N much larger than M, JOBZ='S') */
00516 
00517                     wrkbl = *m + *m * ilaenv_(&c__1, "CGELQF", " ", m, n, &
00518                             c_n1, &c_n1);
00519 /* Computing MAX */
00520                     i__1 = wrkbl, i__2 = *m + *m * ilaenv_(&c__1, "CUNGLQ", 
00521                             " ", m, n, m, &c_n1);
00522                     wrkbl = max(i__1,i__2);
00523 /* Computing MAX */
00524                     i__1 = wrkbl, i__2 = (*m << 1) + (*m << 1) * ilaenv_(&
00525                             c__1, "CGEBRD", " ", m, m, &c_n1, &c_n1);
00526                     wrkbl = max(i__1,i__2);
00527 /* Computing MAX */
00528                     i__1 = wrkbl, i__2 = (*m << 1) + *m * ilaenv_(&c__1, 
00529                             "CUNMBR", "PRC", m, m, m, &c_n1);
00530                     wrkbl = max(i__1,i__2);
00531 /* Computing MAX */
00532                     i__1 = wrkbl, i__2 = (*m << 1) + *m * ilaenv_(&c__1, 
00533                             "CUNMBR", "QLN", m, m, m, &c_n1);
00534                     wrkbl = max(i__1,i__2);
00535                     maxwrk = *m * *m + wrkbl;
00536                     minwrk = *m * *m + *m * 3;
00537                 } else if (wntqa) {
00538 
00539 /*                 Path 4t (N much larger than M, JOBZ='A') */
00540 
00541                     wrkbl = *m + *m * ilaenv_(&c__1, "CGELQF", " ", m, n, &
00542                             c_n1, &c_n1);
00543 /* Computing MAX */
00544                     i__1 = wrkbl, i__2 = *m + *n * ilaenv_(&c__1, "CUNGLQ", 
00545                             " ", n, n, m, &c_n1);
00546                     wrkbl = max(i__1,i__2);
00547 /* Computing MAX */
00548                     i__1 = wrkbl, i__2 = (*m << 1) + (*m << 1) * ilaenv_(&
00549                             c__1, "CGEBRD", " ", m, m, &c_n1, &c_n1);
00550                     wrkbl = max(i__1,i__2);
00551 /* Computing MAX */
00552                     i__1 = wrkbl, i__2 = (*m << 1) + *m * ilaenv_(&c__1, 
00553                             "CUNMBR", "PRC", m, m, m, &c_n1);
00554                     wrkbl = max(i__1,i__2);
00555 /* Computing MAX */
00556                     i__1 = wrkbl, i__2 = (*m << 1) + *m * ilaenv_(&c__1, 
00557                             "CUNMBR", "QLN", m, m, m, &c_n1);
00558                     wrkbl = max(i__1,i__2);
00559                     maxwrk = *m * *m + wrkbl;
00560                     minwrk = *m * *m + (*m << 1) + *n;
00561                 }
00562             } else if (*n >= mnthr2) {
00563 
00564 /*              Path 5t (N much larger than M, but not as much as MNTHR1) */
00565 
00566                 maxwrk = (*m << 1) + (*m + *n) * ilaenv_(&c__1, "CGEBRD", 
00567                         " ", m, n, &c_n1, &c_n1);
00568                 minwrk = (*m << 1) + *n;
00569                 if (wntqo) {
00570 /* Computing MAX */
00571                     i__1 = maxwrk, i__2 = (*m << 1) + *m * ilaenv_(&c__1, 
00572                             "CUNGBR", "P", m, n, m, &c_n1);
00573                     maxwrk = max(i__1,i__2);
00574 /* Computing MAX */
00575                     i__1 = maxwrk, i__2 = (*m << 1) + *m * ilaenv_(&c__1, 
00576                             "CUNGBR", "Q", m, m, n, &c_n1);
00577                     maxwrk = max(i__1,i__2);
00578                     maxwrk += *m * *n;
00579                     minwrk += *m * *m;
00580                 } else if (wntqs) {
00581 /* Computing MAX */
00582                     i__1 = maxwrk, i__2 = (*m << 1) + *m * ilaenv_(&c__1, 
00583                             "CUNGBR", "P", m, n, m, &c_n1);
00584                     maxwrk = max(i__1,i__2);
00585 /* Computing MAX */
00586                     i__1 = maxwrk, i__2 = (*m << 1) + *m * ilaenv_(&c__1, 
00587                             "CUNGBR", "Q", m, m, n, &c_n1);
00588                     maxwrk = max(i__1,i__2);
00589                 } else if (wntqa) {
00590 /* Computing MAX */
00591                     i__1 = maxwrk, i__2 = (*m << 1) + *n * ilaenv_(&c__1, 
00592                             "CUNGBR", "P", n, n, m, &c_n1);
00593                     maxwrk = max(i__1,i__2);
00594 /* Computing MAX */
00595                     i__1 = maxwrk, i__2 = (*m << 1) + *m * ilaenv_(&c__1, 
00596                             "CUNGBR", "Q", m, m, n, &c_n1);
00597                     maxwrk = max(i__1,i__2);
00598                 }
00599             } else {
00600 
00601 /*              Path 6t (N greater than M, but not much larger) */
00602 
00603                 maxwrk = (*m << 1) + (*m + *n) * ilaenv_(&c__1, "CGEBRD", 
00604                         " ", m, n, &c_n1, &c_n1);
00605                 minwrk = (*m << 1) + *n;
00606                 if (wntqo) {
00607 /* Computing MAX */
00608                     i__1 = maxwrk, i__2 = (*m << 1) + *m * ilaenv_(&c__1, 
00609                             "CUNMBR", "PRC", m, n, m, &c_n1);
00610                     maxwrk = max(i__1,i__2);
00611 /* Computing MAX */
00612                     i__1 = maxwrk, i__2 = (*m << 1) + *m * ilaenv_(&c__1, 
00613                             "CUNMBR", "QLN", m, m, n, &c_n1);
00614                     maxwrk = max(i__1,i__2);
00615                     maxwrk += *m * *n;
00616                     minwrk += *m * *m;
00617                 } else if (wntqs) {
00618 /* Computing MAX */
00619                     i__1 = maxwrk, i__2 = (*m << 1) + *m * ilaenv_(&c__1, 
00620                             "CUNGBR", "PRC", m, n, m, &c_n1);
00621                     maxwrk = max(i__1,i__2);
00622 /* Computing MAX */
00623                     i__1 = maxwrk, i__2 = (*m << 1) + *m * ilaenv_(&c__1, 
00624                             "CUNGBR", "QLN", m, m, n, &c_n1);
00625                     maxwrk = max(i__1,i__2);
00626                 } else if (wntqa) {
00627 /* Computing MAX */
00628                     i__1 = maxwrk, i__2 = (*m << 1) + *n * ilaenv_(&c__1, 
00629                             "CUNGBR", "PRC", n, n, m, &c_n1);
00630                     maxwrk = max(i__1,i__2);
00631 /* Computing MAX */
00632                     i__1 = maxwrk, i__2 = (*m << 1) + *m * ilaenv_(&c__1, 
00633                             "CUNGBR", "QLN", m, m, n, &c_n1);
00634                     maxwrk = max(i__1,i__2);
00635                 }
00636             }
00637         }
00638         maxwrk = max(maxwrk,minwrk);
00639     }
00640     if (*info == 0) {
00641         work[1].r = (real) maxwrk, work[1].i = 0.f;
00642         if (*lwork < minwrk && *lwork != -1) {
00643             *info = -13;
00644         }
00645     }
00646 
00647 /*     Quick returns */
00648 
00649     if (*info != 0) {
00650         i__1 = -(*info);
00651         xerbla_("CGESDD", &i__1);
00652         return 0;
00653     }
00654     if (*lwork == -1) {
00655         return 0;
00656     }
00657     if (*m == 0 || *n == 0) {
00658         return 0;
00659     }
00660 
00661 /*     Get machine constants */
00662 
00663     eps = slamch_("P");
00664     smlnum = sqrt(slamch_("S")) / eps;
00665     bignum = 1.f / smlnum;
00666 
00667 /*     Scale A if max element outside range [SMLNUM,BIGNUM] */
00668 
00669     anrm = clange_("M", m, n, &a[a_offset], lda, dum);
00670     iscl = 0;
00671     if (anrm > 0.f && anrm < smlnum) {
00672         iscl = 1;
00673         clascl_("G", &c__0, &c__0, &anrm, &smlnum, m, n, &a[a_offset], lda, &
00674                 ierr);
00675     } else if (anrm > bignum) {
00676         iscl = 1;
00677         clascl_("G", &c__0, &c__0, &anrm, &bignum, m, n, &a[a_offset], lda, &
00678                 ierr);
00679     }
00680 
00681     if (*m >= *n) {
00682 
00683 /*        A has at least as many rows as columns. If A has sufficiently */
00684 /*        more rows than columns, first reduce using the QR */
00685 /*        decomposition (if sufficient workspace available) */
00686 
00687         if (*m >= mnthr1) {
00688 
00689             if (wntqn) {
00690 
00691 /*              Path 1 (M much larger than N, JOBZ='N') */
00692 /*              No singular vectors to be computed */
00693 
00694                 itau = 1;
00695                 nwork = itau + *n;
00696 
00697 /*              Compute A=Q*R */
00698 /*              (CWorkspace: need 2*N, prefer N+N*NB) */
00699 /*              (RWorkspace: need 0) */
00700 
00701                 i__1 = *lwork - nwork + 1;
00702                 cgeqrf_(m, n, &a[a_offset], lda, &work[itau], &work[nwork], &
00703                         i__1, &ierr);
00704 
00705 /*              Zero out below R */
00706 
00707                 i__1 = *n - 1;
00708                 i__2 = *n - 1;
00709                 claset_("L", &i__1, &i__2, &c_b1, &c_b1, &a[a_dim1 + 2], lda);
00710                 ie = 1;
00711                 itauq = 1;
00712                 itaup = itauq + *n;
00713                 nwork = itaup + *n;
00714 
00715 /*              Bidiagonalize R in A */
00716 /*              (CWorkspace: need 3*N, prefer 2*N+2*N*NB) */
00717 /*              (RWorkspace: need N) */
00718 
00719                 i__1 = *lwork - nwork + 1;
00720                 cgebrd_(n, n, &a[a_offset], lda, &s[1], &rwork[ie], &work[
00721                         itauq], &work[itaup], &work[nwork], &i__1, &ierr);
00722                 nrwork = ie + *n;
00723 
00724 /*              Perform bidiagonal SVD, compute singular values only */
00725 /*              (CWorkspace: 0) */
00726 /*              (RWorkspace: need BDSPAN) */
00727 
00728                 sbdsdc_("U", "N", n, &s[1], &rwork[ie], dum, &c__1, dum, &
00729                         c__1, dum, idum, &rwork[nrwork], &iwork[1], info);
00730 
00731             } else if (wntqo) {
00732 
00733 /*              Path 2 (M much larger than N, JOBZ='O') */
00734 /*              N left singular vectors to be overwritten on A and */
00735 /*              N right singular vectors to be computed in VT */
00736 
00737                 iu = 1;
00738 
00739 /*              WORK(IU) is N by N */
00740 
00741                 ldwrku = *n;
00742                 ir = iu + ldwrku * *n;
00743                 if (*lwork >= *m * *n + *n * *n + *n * 3) {
00744 
00745 /*                 WORK(IR) is M by N */
00746 
00747                     ldwrkr = *m;
00748                 } else {
00749                     ldwrkr = (*lwork - *n * *n - *n * 3) / *n;
00750                 }
00751                 itau = ir + ldwrkr * *n;
00752                 nwork = itau + *n;
00753 
00754 /*              Compute A=Q*R */
00755 /*              (CWorkspace: need N*N+2*N, prefer M*N+N+N*NB) */
00756 /*              (RWorkspace: 0) */
00757 
00758                 i__1 = *lwork - nwork + 1;
00759                 cgeqrf_(m, n, &a[a_offset], lda, &work[itau], &work[nwork], &
00760                         i__1, &ierr);
00761 
00762 /*              Copy R to WORK( IR ), zeroing out below it */
00763 
00764                 clacpy_("U", n, n, &a[a_offset], lda, &work[ir], &ldwrkr);
00765                 i__1 = *n - 1;
00766                 i__2 = *n - 1;
00767                 claset_("L", &i__1, &i__2, &c_b1, &c_b1, &work[ir + 1], &
00768                         ldwrkr);
00769 
00770 /*              Generate Q in A */
00771 /*              (CWorkspace: need 2*N, prefer N+N*NB) */
00772 /*              (RWorkspace: 0) */
00773 
00774                 i__1 = *lwork - nwork + 1;
00775                 cungqr_(m, n, n, &a[a_offset], lda, &work[itau], &work[nwork], 
00776                          &i__1, &ierr);
00777                 ie = 1;
00778                 itauq = itau;
00779                 itaup = itauq + *n;
00780                 nwork = itaup + *n;
00781 
00782 /*              Bidiagonalize R in WORK(IR) */
00783 /*              (CWorkspace: need N*N+3*N, prefer M*N+2*N+2*N*NB) */
00784 /*              (RWorkspace: need N) */
00785 
00786                 i__1 = *lwork - nwork + 1;
00787                 cgebrd_(n, n, &work[ir], &ldwrkr, &s[1], &rwork[ie], &work[
00788                         itauq], &work[itaup], &work[nwork], &i__1, &ierr);
00789 
00790 /*              Perform bidiagonal SVD, computing left singular vectors */
00791 /*              of R in WORK(IRU) and computing right singular vectors */
00792 /*              of R in WORK(IRVT) */
00793 /*              (CWorkspace: need 0) */
00794 /*              (RWorkspace: need BDSPAC) */
00795 
00796                 iru = ie + *n;
00797                 irvt = iru + *n * *n;
00798                 nrwork = irvt + *n * *n;
00799                 sbdsdc_("U", "I", n, &s[1], &rwork[ie], &rwork[iru], n, &
00800                         rwork[irvt], n, dum, idum, &rwork[nrwork], &iwork[1], 
00801                         info);
00802 
00803 /*              Copy real matrix RWORK(IRU) to complex matrix WORK(IU) */
00804 /*              Overwrite WORK(IU) by the left singular vectors of R */
00805 /*              (CWorkspace: need 2*N*N+3*N, prefer M*N+N*N+2*N+N*NB) */
00806 /*              (RWorkspace: 0) */
00807 
00808                 clacp2_("F", n, n, &rwork[iru], n, &work[iu], &ldwrku);
00809                 i__1 = *lwork - nwork + 1;
00810                 cunmbr_("Q", "L", "N", n, n, n, &work[ir], &ldwrkr, &work[
00811                         itauq], &work[iu], &ldwrku, &work[nwork], &i__1, &
00812                         ierr);
00813 
00814 /*              Copy real matrix RWORK(IRVT) to complex matrix VT */
00815 /*              Overwrite VT by the right singular vectors of R */
00816 /*              (CWorkspace: need N*N+3*N, prefer M*N+2*N+N*NB) */
00817 /*              (RWorkspace: 0) */
00818 
00819                 clacp2_("F", n, n, &rwork[irvt], n, &vt[vt_offset], ldvt);
00820                 i__1 = *lwork - nwork + 1;
00821                 cunmbr_("P", "R", "C", n, n, n, &work[ir], &ldwrkr, &work[
00822                         itaup], &vt[vt_offset], ldvt, &work[nwork], &i__1, &
00823                         ierr);
00824 
00825 /*              Multiply Q in A by left singular vectors of R in */
00826 /*              WORK(IU), storing result in WORK(IR) and copying to A */
00827 /*              (CWorkspace: need 2*N*N, prefer N*N+M*N) */
00828 /*              (RWorkspace: 0) */
00829 
00830                 i__1 = *m;
00831                 i__2 = ldwrkr;
00832                 for (i__ = 1; i__2 < 0 ? i__ >= i__1 : i__ <= i__1; i__ += 
00833                         i__2) {
00834 /* Computing MIN */
00835                     i__3 = *m - i__ + 1;
00836                     chunk = min(i__3,ldwrkr);
00837                     cgemm_("N", "N", &chunk, n, n, &c_b2, &a[i__ + a_dim1], 
00838                             lda, &work[iu], &ldwrku, &c_b1, &work[ir], &
00839                             ldwrkr);
00840                     clacpy_("F", &chunk, n, &work[ir], &ldwrkr, &a[i__ + 
00841                             a_dim1], lda);
00842 /* L10: */
00843                 }
00844 
00845             } else if (wntqs) {
00846 
00847 /*              Path 3 (M much larger than N, JOBZ='S') */
00848 /*              N left singular vectors to be computed in U and */
00849 /*              N right singular vectors to be computed in VT */
00850 
00851                 ir = 1;
00852 
00853 /*              WORK(IR) is N by N */
00854 
00855                 ldwrkr = *n;
00856                 itau = ir + ldwrkr * *n;
00857                 nwork = itau + *n;
00858 
00859 /*              Compute A=Q*R */
00860 /*              (CWorkspace: need N*N+2*N, prefer N*N+N+N*NB) */
00861 /*              (RWorkspace: 0) */
00862 
00863                 i__2 = *lwork - nwork + 1;
00864                 cgeqrf_(m, n, &a[a_offset], lda, &work[itau], &work[nwork], &
00865                         i__2, &ierr);
00866 
00867 /*              Copy R to WORK(IR), zeroing out below it */
00868 
00869                 clacpy_("U", n, n, &a[a_offset], lda, &work[ir], &ldwrkr);
00870                 i__2 = *n - 1;
00871                 i__1 = *n - 1;
00872                 claset_("L", &i__2, &i__1, &c_b1, &c_b1, &work[ir + 1], &
00873                         ldwrkr);
00874 
00875 /*              Generate Q in A */
00876 /*              (CWorkspace: need 2*N, prefer N+N*NB) */
00877 /*              (RWorkspace: 0) */
00878 
00879                 i__2 = *lwork - nwork + 1;
00880                 cungqr_(m, n, n, &a[a_offset], lda, &work[itau], &work[nwork], 
00881                          &i__2, &ierr);
00882                 ie = 1;
00883                 itauq = itau;
00884                 itaup = itauq + *n;
00885                 nwork = itaup + *n;
00886 
00887 /*              Bidiagonalize R in WORK(IR) */
00888 /*              (CWorkspace: need N*N+3*N, prefer N*N+2*N+2*N*NB) */
00889 /*              (RWorkspace: need N) */
00890 
00891                 i__2 = *lwork - nwork + 1;
00892                 cgebrd_(n, n, &work[ir], &ldwrkr, &s[1], &rwork[ie], &work[
00893                         itauq], &work[itaup], &work[nwork], &i__2, &ierr);
00894 
00895 /*              Perform bidiagonal SVD, computing left singular vectors */
00896 /*              of bidiagonal matrix in RWORK(IRU) and computing right */
00897 /*              singular vectors of bidiagonal matrix in RWORK(IRVT) */
00898 /*              (CWorkspace: need 0) */
00899 /*              (RWorkspace: need BDSPAC) */
00900 
00901                 iru = ie + *n;
00902                 irvt = iru + *n * *n;
00903                 nrwork = irvt + *n * *n;
00904                 sbdsdc_("U", "I", n, &s[1], &rwork[ie], &rwork[iru], n, &
00905                         rwork[irvt], n, dum, idum, &rwork[nrwork], &iwork[1], 
00906                         info);
00907 
00908 /*              Copy real matrix RWORK(IRU) to complex matrix U */
00909 /*              Overwrite U by left singular vectors of R */
00910 /*              (CWorkspace: need N*N+3*N, prefer N*N+2*N+N*NB) */
00911 /*              (RWorkspace: 0) */
00912 
00913                 clacp2_("F", n, n, &rwork[iru], n, &u[u_offset], ldu);
00914                 i__2 = *lwork - nwork + 1;
00915                 cunmbr_("Q", "L", "N", n, n, n, &work[ir], &ldwrkr, &work[
00916                         itauq], &u[u_offset], ldu, &work[nwork], &i__2, &ierr);
00917 
00918 /*              Copy real matrix RWORK(IRVT) to complex matrix VT */
00919 /*              Overwrite VT by right singular vectors of R */
00920 /*              (CWorkspace: need N*N+3*N, prefer N*N+2*N+N*NB) */
00921 /*              (RWorkspace: 0) */
00922 
00923                 clacp2_("F", n, n, &rwork[irvt], n, &vt[vt_offset], ldvt);
00924                 i__2 = *lwork - nwork + 1;
00925                 cunmbr_("P", "R", "C", n, n, n, &work[ir], &ldwrkr, &work[
00926                         itaup], &vt[vt_offset], ldvt, &work[nwork], &i__2, &
00927                         ierr);
00928 
00929 /*              Multiply Q in A by left singular vectors of R in */
00930 /*              WORK(IR), storing result in U */
00931 /*              (CWorkspace: need N*N) */
00932 /*              (RWorkspace: 0) */
00933 
00934                 clacpy_("F", n, n, &u[u_offset], ldu, &work[ir], &ldwrkr);
00935                 cgemm_("N", "N", m, n, n, &c_b2, &a[a_offset], lda, &work[ir], 
00936                          &ldwrkr, &c_b1, &u[u_offset], ldu);
00937 
00938             } else if (wntqa) {
00939 
00940 /*              Path 4 (M much larger than N, JOBZ='A') */
00941 /*              M left singular vectors to be computed in U and */
00942 /*              N right singular vectors to be computed in VT */
00943 
00944                 iu = 1;
00945 
00946 /*              WORK(IU) is N by N */
00947 
00948                 ldwrku = *n;
00949                 itau = iu + ldwrku * *n;
00950                 nwork = itau + *n;
00951 
00952 /*              Compute A=Q*R, copying result to U */
00953 /*              (CWorkspace: need 2*N, prefer N+N*NB) */
00954 /*              (RWorkspace: 0) */
00955 
00956                 i__2 = *lwork - nwork + 1;
00957                 cgeqrf_(m, n, &a[a_offset], lda, &work[itau], &work[nwork], &
00958                         i__2, &ierr);
00959                 clacpy_("L", m, n, &a[a_offset], lda, &u[u_offset], ldu);
00960 
00961 /*              Generate Q in U */
00962 /*              (CWorkspace: need N+M, prefer N+M*NB) */
00963 /*              (RWorkspace: 0) */
00964 
00965                 i__2 = *lwork - nwork + 1;
00966                 cungqr_(m, m, n, &u[u_offset], ldu, &work[itau], &work[nwork], 
00967                          &i__2, &ierr);
00968 
00969 /*              Produce R in A, zeroing out below it */
00970 
00971                 i__2 = *n - 1;
00972                 i__1 = *n - 1;
00973                 claset_("L", &i__2, &i__1, &c_b1, &c_b1, &a[a_dim1 + 2], lda);
00974                 ie = 1;
00975                 itauq = itau;
00976                 itaup = itauq + *n;
00977                 nwork = itaup + *n;
00978 
00979 /*              Bidiagonalize R in A */
00980 /*              (CWorkspace: need 3*N, prefer 2*N+2*N*NB) */
00981 /*              (RWorkspace: need N) */
00982 
00983                 i__2 = *lwork - nwork + 1;
00984                 cgebrd_(n, n, &a[a_offset], lda, &s[1], &rwork[ie], &work[
00985                         itauq], &work[itaup], &work[nwork], &i__2, &ierr);
00986                 iru = ie + *n;
00987                 irvt = iru + *n * *n;
00988                 nrwork = irvt + *n * *n;
00989 
00990 /*              Perform bidiagonal SVD, computing left singular vectors */
00991 /*              of bidiagonal matrix in RWORK(IRU) and computing right */
00992 /*              singular vectors of bidiagonal matrix in RWORK(IRVT) */
00993 /*              (CWorkspace: need 0) */
00994 /*              (RWorkspace: need BDSPAC) */
00995 
00996                 sbdsdc_("U", "I", n, &s[1], &rwork[ie], &rwork[iru], n, &
00997                         rwork[irvt], n, dum, idum, &rwork[nrwork], &iwork[1], 
00998                         info);
00999 
01000 /*              Copy real matrix RWORK(IRU) to complex matrix WORK(IU) */
01001 /*              Overwrite WORK(IU) by left singular vectors of R */
01002 /*              (CWorkspace: need N*N+3*N, prefer N*N+2*N+N*NB) */
01003 /*              (RWorkspace: 0) */
01004 
01005                 clacp2_("F", n, n, &rwork[iru], n, &work[iu], &ldwrku);
01006                 i__2 = *lwork - nwork + 1;
01007                 cunmbr_("Q", "L", "N", n, n, n, &a[a_offset], lda, &work[
01008                         itauq], &work[iu], &ldwrku, &work[nwork], &i__2, &
01009                         ierr);
01010 
01011 /*              Copy real matrix RWORK(IRVT) to complex matrix VT */
01012 /*              Overwrite VT by right singular vectors of R */
01013 /*              (CWorkspace: need 3*N, prefer 2*N+N*NB) */
01014 /*              (RWorkspace: 0) */
01015 
01016                 clacp2_("F", n, n, &rwork[irvt], n, &vt[vt_offset], ldvt);
01017                 i__2 = *lwork - nwork + 1;
01018                 cunmbr_("P", "R", "C", n, n, n, &a[a_offset], lda, &work[
01019                         itaup], &vt[vt_offset], ldvt, &work[nwork], &i__2, &
01020                         ierr);
01021 
01022 /*              Multiply Q in U by left singular vectors of R in */
01023 /*              WORK(IU), storing result in A */
01024 /*              (CWorkspace: need N*N) */
01025 /*              (RWorkspace: 0) */
01026 
01027                 cgemm_("N", "N", m, n, n, &c_b2, &u[u_offset], ldu, &work[iu], 
01028                          &ldwrku, &c_b1, &a[a_offset], lda);
01029 
01030 /*              Copy left singular vectors of A from A to U */
01031 
01032                 clacpy_("F", m, n, &a[a_offset], lda, &u[u_offset], ldu);
01033 
01034             }
01035 
01036         } else if (*m >= mnthr2) {
01037 
01038 /*           MNTHR2 <= M < MNTHR1 */
01039 
01040 /*           Path 5 (M much larger than N, but not as much as MNTHR1) */
01041 /*           Reduce to bidiagonal form without QR decomposition, use */
01042 /*           CUNGBR and matrix multiplication to compute singular vectors */
01043 
01044             ie = 1;
01045             nrwork = ie + *n;
01046             itauq = 1;
01047             itaup = itauq + *n;
01048             nwork = itaup + *n;
01049 
01050 /*           Bidiagonalize A */
01051 /*           (CWorkspace: need 2*N+M, prefer 2*N+(M+N)*NB) */
01052 /*           (RWorkspace: need N) */
01053 
01054             i__2 = *lwork - nwork + 1;
01055             cgebrd_(m, n, &a[a_offset], lda, &s[1], &rwork[ie], &work[itauq], 
01056                     &work[itaup], &work[nwork], &i__2, &ierr);
01057             if (wntqn) {
01058 
01059 /*              Compute singular values only */
01060 /*              (Cworkspace: 0) */
01061 /*              (Rworkspace: need BDSPAN) */
01062 
01063                 sbdsdc_("U", "N", n, &s[1], &rwork[ie], dum, &c__1, dum, &
01064                         c__1, dum, idum, &rwork[nrwork], &iwork[1], info);
01065             } else if (wntqo) {
01066                 iu = nwork;
01067                 iru = nrwork;
01068                 irvt = iru + *n * *n;
01069                 nrwork = irvt + *n * *n;
01070 
01071 /*              Copy A to VT, generate P**H */
01072 /*              (Cworkspace: need 2*N, prefer N+N*NB) */
01073 /*              (Rworkspace: 0) */
01074 
01075                 clacpy_("U", n, n, &a[a_offset], lda, &vt[vt_offset], ldvt);
01076                 i__2 = *lwork - nwork + 1;
01077                 cungbr_("P", n, n, n, &vt[vt_offset], ldvt, &work[itaup], &
01078                         work[nwork], &i__2, &ierr);
01079 
01080 /*              Generate Q in A */
01081 /*              (CWorkspace: need 2*N, prefer N+N*NB) */
01082 /*              (RWorkspace: 0) */
01083 
01084                 i__2 = *lwork - nwork + 1;
01085                 cungbr_("Q", m, n, n, &a[a_offset], lda, &work[itauq], &work[
01086                         nwork], &i__2, &ierr);
01087 
01088                 if (*lwork >= *m * *n + *n * 3) {
01089 
01090 /*                 WORK( IU ) is M by N */
01091 
01092                     ldwrku = *m;
01093                 } else {
01094 
01095 /*                 WORK(IU) is LDWRKU by N */
01096 
01097                     ldwrku = (*lwork - *n * 3) / *n;
01098                 }
01099                 nwork = iu + ldwrku * *n;
01100 
01101 /*              Perform bidiagonal SVD, computing left singular vectors */
01102 /*              of bidiagonal matrix in RWORK(IRU) and computing right */
01103 /*              singular vectors of bidiagonal matrix in RWORK(IRVT) */
01104 /*              (CWorkspace: need 0) */
01105 /*              (RWorkspace: need BDSPAC) */
01106 
01107                 sbdsdc_("U", "I", n, &s[1], &rwork[ie], &rwork[iru], n, &
01108                         rwork[irvt], n, dum, idum, &rwork[nrwork], &iwork[1], 
01109                         info);
01110 
01111 /*              Multiply real matrix RWORK(IRVT) by P**H in VT, */
01112 /*              storing the result in WORK(IU), copying to VT */
01113 /*              (Cworkspace: need 0) */
01114 /*              (Rworkspace: need 3*N*N) */
01115 
01116                 clarcm_(n, n, &rwork[irvt], n, &vt[vt_offset], ldvt, &work[iu]
01117 , &ldwrku, &rwork[nrwork]);
01118                 clacpy_("F", n, n, &work[iu], &ldwrku, &vt[vt_offset], ldvt);
01119 
01120 /*              Multiply Q in A by real matrix RWORK(IRU), storing the */
01121 /*              result in WORK(IU), copying to A */
01122 /*              (CWorkspace: need N*N, prefer M*N) */
01123 /*              (Rworkspace: need 3*N*N, prefer N*N+2*M*N) */
01124 
01125                 nrwork = irvt;
01126                 i__2 = *m;
01127                 i__1 = ldwrku;
01128                 for (i__ = 1; i__1 < 0 ? i__ >= i__2 : i__ <= i__2; i__ += 
01129                         i__1) {
01130 /* Computing MIN */
01131                     i__3 = *m - i__ + 1;
01132                     chunk = min(i__3,ldwrku);
01133                     clacrm_(&chunk, n, &a[i__ + a_dim1], lda, &rwork[iru], n, 
01134                             &work[iu], &ldwrku, &rwork[nrwork]);
01135                     clacpy_("F", &chunk, n, &work[iu], &ldwrku, &a[i__ + 
01136                             a_dim1], lda);
01137 /* L20: */
01138                 }
01139 
01140             } else if (wntqs) {
01141 
01142 /*              Copy A to VT, generate P**H */
01143 /*              (Cworkspace: need 2*N, prefer N+N*NB) */
01144 /*              (Rworkspace: 0) */
01145 
01146                 clacpy_("U", n, n, &a[a_offset], lda, &vt[vt_offset], ldvt);
01147                 i__1 = *lwork - nwork + 1;
01148                 cungbr_("P", n, n, n, &vt[vt_offset], ldvt, &work[itaup], &
01149                         work[nwork], &i__1, &ierr);
01150 
01151 /*              Copy A to U, generate Q */
01152 /*              (Cworkspace: need 2*N, prefer N+N*NB) */
01153 /*              (Rworkspace: 0) */
01154 
01155                 clacpy_("L", m, n, &a[a_offset], lda, &u[u_offset], ldu);
01156                 i__1 = *lwork - nwork + 1;
01157                 cungbr_("Q", m, n, n, &u[u_offset], ldu, &work[itauq], &work[
01158                         nwork], &i__1, &ierr);
01159 
01160 /*              Perform bidiagonal SVD, computing left singular vectors */
01161 /*              of bidiagonal matrix in RWORK(IRU) and computing right */
01162 /*              singular vectors of bidiagonal matrix in RWORK(IRVT) */
01163 /*              (CWorkspace: need 0) */
01164 /*              (RWorkspace: need BDSPAC) */
01165 
01166                 iru = nrwork;
01167                 irvt = iru + *n * *n;
01168                 nrwork = irvt + *n * *n;
01169                 sbdsdc_("U", "I", n, &s[1], &rwork[ie], &rwork[iru], n, &
01170                         rwork[irvt], n, dum, idum, &rwork[nrwork], &iwork[1], 
01171                         info);
01172 
01173 /*              Multiply real matrix RWORK(IRVT) by P**H in VT, */
01174 /*              storing the result in A, copying to VT */
01175 /*              (Cworkspace: need 0) */
01176 /*              (Rworkspace: need 3*N*N) */
01177 
01178                 clarcm_(n, n, &rwork[irvt], n, &vt[vt_offset], ldvt, &a[
01179                         a_offset], lda, &rwork[nrwork]);
01180                 clacpy_("F", n, n, &a[a_offset], lda, &vt[vt_offset], ldvt);
01181 
01182 /*              Multiply Q in U by real matrix RWORK(IRU), storing the */
01183 /*              result in A, copying to U */
01184 /*              (CWorkspace: need 0) */
01185 /*              (Rworkspace: need N*N+2*M*N) */
01186 
01187                 nrwork = irvt;
01188                 clacrm_(m, n, &u[u_offset], ldu, &rwork[iru], n, &a[a_offset], 
01189                          lda, &rwork[nrwork]);
01190                 clacpy_("F", m, n, &a[a_offset], lda, &u[u_offset], ldu);
01191             } else {
01192 
01193 /*              Copy A to VT, generate P**H */
01194 /*              (Cworkspace: need 2*N, prefer N+N*NB) */
01195 /*              (Rworkspace: 0) */
01196 
01197                 clacpy_("U", n, n, &a[a_offset], lda, &vt[vt_offset], ldvt);
01198                 i__1 = *lwork - nwork + 1;
01199                 cungbr_("P", n, n, n, &vt[vt_offset], ldvt, &work[itaup], &
01200                         work[nwork], &i__1, &ierr);
01201 
01202 /*              Copy A to U, generate Q */
01203 /*              (Cworkspace: need 2*N, prefer N+N*NB) */
01204 /*              (Rworkspace: 0) */
01205 
01206                 clacpy_("L", m, n, &a[a_offset], lda, &u[u_offset], ldu);
01207                 i__1 = *lwork - nwork + 1;
01208                 cungbr_("Q", m, m, n, &u[u_offset], ldu, &work[itauq], &work[
01209                         nwork], &i__1, &ierr);
01210 
01211 /*              Perform bidiagonal SVD, computing left singular vectors */
01212 /*              of bidiagonal matrix in RWORK(IRU) and computing right */
01213 /*              singular vectors of bidiagonal matrix in RWORK(IRVT) */
01214 /*              (CWorkspace: need 0) */
01215 /*              (RWorkspace: need BDSPAC) */
01216 
01217                 iru = nrwork;
01218                 irvt = iru + *n * *n;
01219                 nrwork = irvt + *n * *n;
01220                 sbdsdc_("U", "I", n, &s[1], &rwork[ie], &rwork[iru], n, &
01221                         rwork[irvt], n, dum, idum, &rwork[nrwork], &iwork[1], 
01222                         info);
01223 
01224 /*              Multiply real matrix RWORK(IRVT) by P**H in VT, */
01225 /*              storing the result in A, copying to VT */
01226 /*              (Cworkspace: need 0) */
01227 /*              (Rworkspace: need 3*N*N) */
01228 
01229                 clarcm_(n, n, &rwork[irvt], n, &vt[vt_offset], ldvt, &a[
01230                         a_offset], lda, &rwork[nrwork]);
01231                 clacpy_("F", n, n, &a[a_offset], lda, &vt[vt_offset], ldvt);
01232 
01233 /*              Multiply Q in U by real matrix RWORK(IRU), storing the */
01234 /*              result in A, copying to U */
01235 /*              (CWorkspace: 0) */
01236 /*              (Rworkspace: need 3*N*N) */
01237 
01238                 nrwork = irvt;
01239                 clacrm_(m, n, &u[u_offset], ldu, &rwork[iru], n, &a[a_offset], 
01240                          lda, &rwork[nrwork]);
01241                 clacpy_("F", m, n, &a[a_offset], lda, &u[u_offset], ldu);
01242             }
01243 
01244         } else {
01245 
01246 /*           M .LT. MNTHR2 */
01247 
01248 /*           Path 6 (M at least N, but not much larger) */
01249 /*           Reduce to bidiagonal form without QR decomposition */
01250 /*           Use CUNMBR to compute singular vectors */
01251 
01252             ie = 1;
01253             nrwork = ie + *n;
01254             itauq = 1;
01255             itaup = itauq + *n;
01256             nwork = itaup + *n;
01257 
01258 /*           Bidiagonalize A */
01259 /*           (CWorkspace: need 2*N+M, prefer 2*N+(M+N)*NB) */
01260 /*           (RWorkspace: need N) */
01261 
01262             i__1 = *lwork - nwork + 1;
01263             cgebrd_(m, n, &a[a_offset], lda, &s[1], &rwork[ie], &work[itauq], 
01264                     &work[itaup], &work[nwork], &i__1, &ierr);
01265             if (wntqn) {
01266 
01267 /*              Compute singular values only */
01268 /*              (Cworkspace: 0) */
01269 /*              (Rworkspace: need BDSPAN) */
01270 
01271                 sbdsdc_("U", "N", n, &s[1], &rwork[ie], dum, &c__1, dum, &
01272                         c__1, dum, idum, &rwork[nrwork], &iwork[1], info);
01273             } else if (wntqo) {
01274                 iu = nwork;
01275                 iru = nrwork;
01276                 irvt = iru + *n * *n;
01277                 nrwork = irvt + *n * *n;
01278                 if (*lwork >= *m * *n + *n * 3) {
01279 
01280 /*                 WORK( IU ) is M by N */
01281 
01282                     ldwrku = *m;
01283                 } else {
01284 
01285 /*                 WORK( IU ) is LDWRKU by N */
01286 
01287                     ldwrku = (*lwork - *n * 3) / *n;
01288                 }
01289                 nwork = iu + ldwrku * *n;
01290 
01291 /*              Perform bidiagonal SVD, computing left singular vectors */
01292 /*              of bidiagonal matrix in RWORK(IRU) and computing right */
01293 /*              singular vectors of bidiagonal matrix in RWORK(IRVT) */
01294 /*              (CWorkspace: need 0) */
01295 /*              (RWorkspace: need BDSPAC) */
01296 
01297                 sbdsdc_("U", "I", n, &s[1], &rwork[ie], &rwork[iru], n, &
01298                         rwork[irvt], n, dum, idum, &rwork[nrwork], &iwork[1], 
01299                         info);
01300 
01301 /*              Copy real matrix RWORK(IRVT) to complex matrix VT */
01302 /*              Overwrite VT by right singular vectors of A */
01303 /*              (Cworkspace: need 2*N, prefer N+N*NB) */
01304 /*              (Rworkspace: need 0) */
01305 
01306                 clacp2_("F", n, n, &rwork[irvt], n, &vt[vt_offset], ldvt);
01307                 i__1 = *lwork - nwork + 1;
01308                 cunmbr_("P", "R", "C", n, n, n, &a[a_offset], lda, &work[
01309                         itaup], &vt[vt_offset], ldvt, &work[nwork], &i__1, &
01310                         ierr);
01311 
01312                 if (*lwork >= *m * *n + *n * 3) {
01313 
01314 /*              Copy real matrix RWORK(IRU) to complex matrix WORK(IU) */
01315 /*              Overwrite WORK(IU) by left singular vectors of A, copying */
01316 /*              to A */
01317 /*              (Cworkspace: need M*N+2*N, prefer M*N+N+N*NB) */
01318 /*              (Rworkspace: need 0) */
01319 
01320                     claset_("F", m, n, &c_b1, &c_b1, &work[iu], &ldwrku);
01321                     clacp2_("F", n, n, &rwork[iru], n, &work[iu], &ldwrku);
01322                     i__1 = *lwork - nwork + 1;
01323                     cunmbr_("Q", "L", "N", m, n, n, &a[a_offset], lda, &work[
01324                             itauq], &work[iu], &ldwrku, &work[nwork], &i__1, &
01325                             ierr);
01326                     clacpy_("F", m, n, &work[iu], &ldwrku, &a[a_offset], lda);
01327                 } else {
01328 
01329 /*                 Generate Q in A */
01330 /*                 (Cworkspace: need 2*N, prefer N+N*NB) */
01331 /*                 (Rworkspace: need 0) */
01332 
01333                     i__1 = *lwork - nwork + 1;
01334                     cungbr_("Q", m, n, n, &a[a_offset], lda, &work[itauq], &
01335                             work[nwork], &i__1, &ierr);
01336 
01337 /*                 Multiply Q in A by real matrix RWORK(IRU), storing the */
01338 /*                 result in WORK(IU), copying to A */
01339 /*                 (CWorkspace: need N*N, prefer M*N) */
01340 /*                 (Rworkspace: need 3*N*N, prefer N*N+2*M*N) */
01341 
01342                     nrwork = irvt;
01343                     i__1 = *m;
01344                     i__2 = ldwrku;
01345                     for (i__ = 1; i__2 < 0 ? i__ >= i__1 : i__ <= i__1; i__ +=
01346                              i__2) {
01347 /* Computing MIN */
01348                         i__3 = *m - i__ + 1;
01349                         chunk = min(i__3,ldwrku);
01350                         clacrm_(&chunk, n, &a[i__ + a_dim1], lda, &rwork[iru], 
01351                                  n, &work[iu], &ldwrku, &rwork[nrwork]);
01352                         clacpy_("F", &chunk, n, &work[iu], &ldwrku, &a[i__ + 
01353                                 a_dim1], lda);
01354 /* L30: */
01355                     }
01356                 }
01357 
01358             } else if (wntqs) {
01359 
01360 /*              Perform bidiagonal SVD, computing left singular vectors */
01361 /*              of bidiagonal matrix in RWORK(IRU) and computing right */
01362 /*              singular vectors of bidiagonal matrix in RWORK(IRVT) */
01363 /*              (CWorkspace: need 0) */
01364 /*              (RWorkspace: need BDSPAC) */
01365 
01366                 iru = nrwork;
01367                 irvt = iru + *n * *n;
01368                 nrwork = irvt + *n * *n;
01369                 sbdsdc_("U", "I", n, &s[1], &rwork[ie], &rwork[iru], n, &
01370                         rwork[irvt], n, dum, idum, &rwork[nrwork], &iwork[1], 
01371                         info);
01372 
01373 /*              Copy real matrix RWORK(IRU) to complex matrix U */
01374 /*              Overwrite U by left singular vectors of A */
01375 /*              (CWorkspace: need 3*N, prefer 2*N+N*NB) */
01376 /*              (RWorkspace: 0) */
01377 
01378                 claset_("F", m, n, &c_b1, &c_b1, &u[u_offset], ldu)
01379                         ;
01380                 clacp2_("F", n, n, &rwork[iru], n, &u[u_offset], ldu);
01381                 i__2 = *lwork - nwork + 1;
01382                 cunmbr_("Q", "L", "N", m, n, n, &a[a_offset], lda, &work[
01383                         itauq], &u[u_offset], ldu, &work[nwork], &i__2, &ierr);
01384 
01385 /*              Copy real matrix RWORK(IRVT) to complex matrix VT */
01386 /*              Overwrite VT by right singular vectors of A */
01387 /*              (CWorkspace: need 3*N, prefer 2*N+N*NB) */
01388 /*              (RWorkspace: 0) */
01389 
01390                 clacp2_("F", n, n, &rwork[irvt], n, &vt[vt_offset], ldvt);
01391                 i__2 = *lwork - nwork + 1;
01392                 cunmbr_("P", "R", "C", n, n, n, &a[a_offset], lda, &work[
01393                         itaup], &vt[vt_offset], ldvt, &work[nwork], &i__2, &
01394                         ierr);
01395             } else {
01396 
01397 /*              Perform bidiagonal SVD, computing left singular vectors */
01398 /*              of bidiagonal matrix in RWORK(IRU) and computing right */
01399 /*              singular vectors of bidiagonal matrix in RWORK(IRVT) */
01400 /*              (CWorkspace: need 0) */
01401 /*              (RWorkspace: need BDSPAC) */
01402 
01403                 iru = nrwork;
01404                 irvt = iru + *n * *n;
01405                 nrwork = irvt + *n * *n;
01406                 sbdsdc_("U", "I", n, &s[1], &rwork[ie], &rwork[iru], n, &
01407                         rwork[irvt], n, dum, idum, &rwork[nrwork], &iwork[1], 
01408                         info);
01409 
01410 /*              Set the right corner of U to identity matrix */
01411 
01412                 claset_("F", m, m, &c_b1, &c_b1, &u[u_offset], ldu)
01413                         ;
01414                 if (*m > *n) {
01415                     i__2 = *m - *n;
01416                     i__1 = *m - *n;
01417                     claset_("F", &i__2, &i__1, &c_b1, &c_b2, &u[*n + 1 + (*n 
01418                             + 1) * u_dim1], ldu);
01419                 }
01420 
01421 /*              Copy real matrix RWORK(IRU) to complex matrix U */
01422 /*              Overwrite U by left singular vectors of A */
01423 /*              (CWorkspace: need 2*N+M, prefer 2*N+M*NB) */
01424 /*              (RWorkspace: 0) */
01425 
01426                 clacp2_("F", n, n, &rwork[iru], n, &u[u_offset], ldu);
01427                 i__2 = *lwork - nwork + 1;
01428                 cunmbr_("Q", "L", "N", m, m, n, &a[a_offset], lda, &work[
01429                         itauq], &u[u_offset], ldu, &work[nwork], &i__2, &ierr);
01430 
01431 /*              Copy real matrix RWORK(IRVT) to complex matrix VT */
01432 /*              Overwrite VT by right singular vectors of A */
01433 /*              (CWorkspace: need 3*N, prefer 2*N+N*NB) */
01434 /*              (RWorkspace: 0) */
01435 
01436                 clacp2_("F", n, n, &rwork[irvt], n, &vt[vt_offset], ldvt);
01437                 i__2 = *lwork - nwork + 1;
01438                 cunmbr_("P", "R", "C", n, n, n, &a[a_offset], lda, &work[
01439                         itaup], &vt[vt_offset], ldvt, &work[nwork], &i__2, &
01440                         ierr);
01441             }
01442 
01443         }
01444 
01445     } else {
01446 
01447 /*        A has more columns than rows. If A has sufficiently more */
01448 /*        columns than rows, first reduce using the LQ decomposition (if */
01449 /*        sufficient workspace available) */
01450 
01451         if (*n >= mnthr1) {
01452 
01453             if (wntqn) {
01454 
01455 /*              Path 1t (N much larger than M, JOBZ='N') */
01456 /*              No singular vectors to be computed */
01457 
01458                 itau = 1;
01459                 nwork = itau + *m;
01460 
01461 /*              Compute A=L*Q */
01462 /*              (CWorkspace: need 2*M, prefer M+M*NB) */
01463 /*              (RWorkspace: 0) */
01464 
01465                 i__2 = *lwork - nwork + 1;
01466                 cgelqf_(m, n, &a[a_offset], lda, &work[itau], &work[nwork], &
01467                         i__2, &ierr);
01468 
01469 /*              Zero out above L */
01470 
01471                 i__2 = *m - 1;
01472                 i__1 = *m - 1;
01473                 claset_("U", &i__2, &i__1, &c_b1, &c_b1, &a[(a_dim1 << 1) + 1]
01474 , lda);
01475                 ie = 1;
01476                 itauq = 1;
01477                 itaup = itauq + *m;
01478                 nwork = itaup + *m;
01479 
01480 /*              Bidiagonalize L in A */
01481 /*              (CWorkspace: need 3*M, prefer 2*M+2*M*NB) */
01482 /*              (RWorkspace: need M) */
01483 
01484                 i__2 = *lwork - nwork + 1;
01485                 cgebrd_(m, m, &a[a_offset], lda, &s[1], &rwork[ie], &work[
01486                         itauq], &work[itaup], &work[nwork], &i__2, &ierr);
01487                 nrwork = ie + *m;
01488 
01489 /*              Perform bidiagonal SVD, compute singular values only */
01490 /*              (CWorkspace: 0) */
01491 /*              (RWorkspace: need BDSPAN) */
01492 
01493                 sbdsdc_("U", "N", m, &s[1], &rwork[ie], dum, &c__1, dum, &
01494                         c__1, dum, idum, &rwork[nrwork], &iwork[1], info);
01495 
01496             } else if (wntqo) {
01497 
01498 /*              Path 2t (N much larger than M, JOBZ='O') */
01499 /*              M right singular vectors to be overwritten on A and */
01500 /*              M left singular vectors to be computed in U */
01501 
01502                 ivt = 1;
01503                 ldwkvt = *m;
01504 
01505 /*              WORK(IVT) is M by M */
01506 
01507                 il = ivt + ldwkvt * *m;
01508                 if (*lwork >= *m * *n + *m * *m + *m * 3) {
01509 
01510 /*                 WORK(IL) M by N */
01511 
01512                     ldwrkl = *m;
01513                     chunk = *n;
01514                 } else {
01515 
01516 /*                 WORK(IL) is M by CHUNK */
01517 
01518                     ldwrkl = *m;
01519                     chunk = (*lwork - *m * *m - *m * 3) / *m;
01520                 }
01521                 itau = il + ldwrkl * chunk;
01522                 nwork = itau + *m;
01523 
01524 /*              Compute A=L*Q */
01525 /*              (CWorkspace: need 2*M, prefer M+M*NB) */
01526 /*              (RWorkspace: 0) */
01527 
01528                 i__2 = *lwork - nwork + 1;
01529                 cgelqf_(m, n, &a[a_offset], lda, &work[itau], &work[nwork], &
01530                         i__2, &ierr);
01531 
01532 /*              Copy L to WORK(IL), zeroing about above it */
01533 
01534                 clacpy_("L", m, m, &a[a_offset], lda, &work[il], &ldwrkl);
01535                 i__2 = *m - 1;
01536                 i__1 = *m - 1;
01537                 claset_("U", &i__2, &i__1, &c_b1, &c_b1, &work[il + ldwrkl], &
01538                         ldwrkl);
01539 
01540 /*              Generate Q in A */
01541 /*              (CWorkspace: need M*M+2*M, prefer M*M+M+M*NB) */
01542 /*              (RWorkspace: 0) */
01543 
01544                 i__2 = *lwork - nwork + 1;
01545                 cunglq_(m, n, m, &a[a_offset], lda, &work[itau], &work[nwork], 
01546                          &i__2, &ierr);
01547                 ie = 1;
01548                 itauq = itau;
01549                 itaup = itauq + *m;
01550                 nwork = itaup + *m;
01551 
01552 /*              Bidiagonalize L in WORK(IL) */
01553 /*              (CWorkspace: need M*M+3*M, prefer M*M+2*M+2*M*NB) */
01554 /*              (RWorkspace: need M) */
01555 
01556                 i__2 = *lwork - nwork + 1;
01557                 cgebrd_(m, m, &work[il], &ldwrkl, &s[1], &rwork[ie], &work[
01558                         itauq], &work[itaup], &work[nwork], &i__2, &ierr);
01559 
01560 /*              Perform bidiagonal SVD, computing left singular vectors */
01561 /*              of bidiagonal matrix in RWORK(IRU) and computing right */
01562 /*              singular vectors of bidiagonal matrix in RWORK(IRVT) */
01563 /*              (CWorkspace: need 0) */
01564 /*              (RWorkspace: need BDSPAC) */
01565 
01566                 iru = ie + *m;
01567                 irvt = iru + *m * *m;
01568                 nrwork = irvt + *m * *m;
01569                 sbdsdc_("U", "I", m, &s[1], &rwork[ie], &rwork[iru], m, &
01570                         rwork[irvt], m, dum, idum, &rwork[nrwork], &iwork[1], 
01571                         info);
01572 
01573 /*              Copy real matrix RWORK(IRU) to complex matrix WORK(IU) */
01574 /*              Overwrite WORK(IU) by the left singular vectors of L */
01575 /*              (CWorkspace: need N*N+3*N, prefer M*N+2*N+N*NB) */
01576 /*              (RWorkspace: 0) */
01577 
01578                 clacp2_("F", m, m, &rwork[iru], m, &u[u_offset], ldu);
01579                 i__2 = *lwork - nwork + 1;
01580                 cunmbr_("Q", "L", "N", m, m, m, &work[il], &ldwrkl, &work[
01581                         itauq], &u[u_offset], ldu, &work[nwork], &i__2, &ierr);
01582 
01583 /*              Copy real matrix RWORK(IRVT) to complex matrix WORK(IVT) */
01584 /*              Overwrite WORK(IVT) by the right singular vectors of L */
01585 /*              (CWorkspace: need N*N+3*N, prefer M*N+2*N+N*NB) */
01586 /*              (RWorkspace: 0) */
01587 
01588                 clacp2_("F", m, m, &rwork[irvt], m, &work[ivt], &ldwkvt);
01589                 i__2 = *lwork - nwork + 1;
01590                 cunmbr_("P", "R", "C", m, m, m, &work[il], &ldwrkl, &work[
01591                         itaup], &work[ivt], &ldwkvt, &work[nwork], &i__2, &
01592                         ierr);
01593 
01594 /*              Multiply right singular vectors of L in WORK(IL) by Q */
01595 /*              in A, storing result in WORK(IL) and copying to A */
01596 /*              (CWorkspace: need 2*M*M, prefer M*M+M*N)) */
01597 /*              (RWorkspace: 0) */
01598 
01599                 i__2 = *n;
01600                 i__1 = chunk;
01601                 for (i__ = 1; i__1 < 0 ? i__ >= i__2 : i__ <= i__2; i__ += 
01602                         i__1) {
01603 /* Computing MIN */
01604                     i__3 = *n - i__ + 1;
01605                     blk = min(i__3,chunk);
01606                     cgemm_("N", "N", m, &blk, m, &c_b2, &work[ivt], m, &a[i__ 
01607                             * a_dim1 + 1], lda, &c_b1, &work[il], &ldwrkl);
01608                     clacpy_("F", m, &blk, &work[il], &ldwrkl, &a[i__ * a_dim1 
01609                             + 1], lda);
01610 /* L40: */
01611                 }
01612 
01613             } else if (wntqs) {
01614 
01615 /*             Path 3t (N much larger than M, JOBZ='S') */
01616 /*             M right singular vectors to be computed in VT and */
01617 /*             M left singular vectors to be computed in U */
01618 
01619                 il = 1;
01620 
01621 /*              WORK(IL) is M by M */
01622 
01623                 ldwrkl = *m;
01624                 itau = il + ldwrkl * *m;
01625                 nwork = itau + *m;
01626 
01627 /*              Compute A=L*Q */
01628 /*              (CWorkspace: need 2*M, prefer M+M*NB) */
01629 /*              (RWorkspace: 0) */
01630 
01631                 i__1 = *lwork - nwork + 1;
01632                 cgelqf_(m, n, &a[a_offset], lda, &work[itau], &work[nwork], &
01633                         i__1, &ierr);
01634 
01635 /*              Copy L to WORK(IL), zeroing out above it */
01636 
01637                 clacpy_("L", m, m, &a[a_offset], lda, &work[il], &ldwrkl);
01638                 i__1 = *m - 1;
01639                 i__2 = *m - 1;
01640                 claset_("U", &i__1, &i__2, &c_b1, &c_b1, &work[il + ldwrkl], &
01641                         ldwrkl);
01642 
01643 /*              Generate Q in A */
01644 /*              (CWorkspace: need M*M+2*M, prefer M*M+M+M*NB) */
01645 /*              (RWorkspace: 0) */
01646 
01647                 i__1 = *lwork - nwork + 1;
01648                 cunglq_(m, n, m, &a[a_offset], lda, &work[itau], &work[nwork], 
01649                          &i__1, &ierr);
01650                 ie = 1;
01651                 itauq = itau;
01652                 itaup = itauq + *m;
01653                 nwork = itaup + *m;
01654 
01655 /*              Bidiagonalize L in WORK(IL) */
01656 /*              (CWorkspace: need M*M+3*M, prefer M*M+2*M+2*M*NB) */
01657 /*              (RWorkspace: need M) */
01658 
01659                 i__1 = *lwork - nwork + 1;
01660                 cgebrd_(m, m, &work[il], &ldwrkl, &s[1], &rwork[ie], &work[
01661                         itauq], &work[itaup], &work[nwork], &i__1, &ierr);
01662 
01663 /*              Perform bidiagonal SVD, computing left singular vectors */
01664 /*              of bidiagonal matrix in RWORK(IRU) and computing right */
01665 /*              singular vectors of bidiagonal matrix in RWORK(IRVT) */
01666 /*              (CWorkspace: need 0) */
01667 /*              (RWorkspace: need BDSPAC) */
01668 
01669                 iru = ie + *m;
01670                 irvt = iru + *m * *m;
01671                 nrwork = irvt + *m * *m;
01672                 sbdsdc_("U", "I", m, &s[1], &rwork[ie], &rwork[iru], m, &
01673                         rwork[irvt], m, dum, idum, &rwork[nrwork], &iwork[1], 
01674                         info);
01675 
01676 /*              Copy real matrix RWORK(IRU) to complex matrix U */
01677 /*              Overwrite U by left singular vectors of L */
01678 /*              (CWorkspace: need M*M+3*M, prefer M*M+2*M+M*NB) */
01679 /*              (RWorkspace: 0) */
01680 
01681                 clacp2_("F", m, m, &rwork[iru], m, &u[u_offset], ldu);
01682                 i__1 = *lwork - nwork + 1;
01683                 cunmbr_("Q", "L", "N", m, m, m, &work[il], &ldwrkl, &work[
01684                         itauq], &u[u_offset], ldu, &work[nwork], &i__1, &ierr);
01685 
01686 /*              Copy real matrix RWORK(IRVT) to complex matrix VT */
01687 /*              Overwrite VT by left singular vectors of L */
01688 /*              (CWorkspace: need M*M+3*M, prefer M*M+2*M+M*NB) */
01689 /*              (RWorkspace: 0) */
01690 
01691                 clacp2_("F", m, m, &rwork[irvt], m, &vt[vt_offset], ldvt);
01692                 i__1 = *lwork - nwork + 1;
01693                 cunmbr_("P", "R", "C", m, m, m, &work[il], &ldwrkl, &work[
01694                         itaup], &vt[vt_offset], ldvt, &work[nwork], &i__1, &
01695                         ierr);
01696 
01697 /*              Copy VT to WORK(IL), multiply right singular vectors of L */
01698 /*              in WORK(IL) by Q in A, storing result in VT */
01699 /*              (CWorkspace: need M*M) */
01700 /*              (RWorkspace: 0) */
01701 
01702                 clacpy_("F", m, m, &vt[vt_offset], ldvt, &work[il], &ldwrkl);
01703                 cgemm_("N", "N", m, n, m, &c_b2, &work[il], &ldwrkl, &a[
01704                         a_offset], lda, &c_b1, &vt[vt_offset], ldvt);
01705 
01706             } else if (wntqa) {
01707 
01708 /*              Path 9t (N much larger than M, JOBZ='A') */
01709 /*              N right singular vectors to be computed in VT and */
01710 /*              M left singular vectors to be computed in U */
01711 
01712                 ivt = 1;
01713 
01714 /*              WORK(IVT) is M by M */
01715 
01716                 ldwkvt = *m;
01717                 itau = ivt + ldwkvt * *m;
01718                 nwork = itau + *m;
01719 
01720 /*              Compute A=L*Q, copying result to VT */
01721 /*              (CWorkspace: need 2*M, prefer M+M*NB) */
01722 /*              (RWorkspace: 0) */
01723 
01724                 i__1 = *lwork - nwork + 1;
01725                 cgelqf_(m, n, &a[a_offset], lda, &work[itau], &work[nwork], &
01726                         i__1, &ierr);
01727                 clacpy_("U", m, n, &a[a_offset], lda, &vt[vt_offset], ldvt);
01728 
01729 /*              Generate Q in VT */
01730 /*              (CWorkspace: need M+N, prefer M+N*NB) */
01731 /*              (RWorkspace: 0) */
01732 
01733                 i__1 = *lwork - nwork + 1;
01734                 cunglq_(n, n, m, &vt[vt_offset], ldvt, &work[itau], &work[
01735                         nwork], &i__1, &ierr);
01736 
01737 /*              Produce L in A, zeroing out above it */
01738 
01739                 i__1 = *m - 1;
01740                 i__2 = *m - 1;
01741                 claset_("U", &i__1, &i__2, &c_b1, &c_b1, &a[(a_dim1 << 1) + 1]
01742 , lda);
01743                 ie = 1;
01744                 itauq = itau;
01745                 itaup = itauq + *m;
01746                 nwork = itaup + *m;
01747 
01748 /*              Bidiagonalize L in A */
01749 /*              (CWorkspace: need M*M+3*M, prefer M*M+2*M+2*M*NB) */
01750 /*              (RWorkspace: need M) */
01751 
01752                 i__1 = *lwork - nwork + 1;
01753                 cgebrd_(m, m, &a[a_offset], lda, &s[1], &rwork[ie], &work[
01754                         itauq], &work[itaup], &work[nwork], &i__1, &ierr);
01755 
01756 /*              Perform bidiagonal SVD, computing left singular vectors */
01757 /*              of bidiagonal matrix in RWORK(IRU) and computing right */
01758 /*              singular vectors of bidiagonal matrix in RWORK(IRVT) */
01759 /*              (CWorkspace: need 0) */
01760 /*              (RWorkspace: need BDSPAC) */
01761 
01762                 iru = ie + *m;
01763                 irvt = iru + *m * *m;
01764                 nrwork = irvt + *m * *m;
01765                 sbdsdc_("U", "I", m, &s[1], &rwork[ie], &rwork[iru], m, &
01766                         rwork[irvt], m, dum, idum, &rwork[nrwork], &iwork[1], 
01767                         info);
01768 
01769 /*              Copy real matrix RWORK(IRU) to complex matrix U */
01770 /*              Overwrite U by left singular vectors of L */
01771 /*              (CWorkspace: need 3*M, prefer 2*M+M*NB) */
01772 /*              (RWorkspace: 0) */
01773 
01774                 clacp2_("F", m, m, &rwork[iru], m, &u[u_offset], ldu);
01775                 i__1 = *lwork - nwork + 1;
01776                 cunmbr_("Q", "L", "N", m, m, m, &a[a_offset], lda, &work[
01777                         itauq], &u[u_offset], ldu, &work[nwork], &i__1, &ierr);
01778 
01779 /*              Copy real matrix RWORK(IRVT) to complex matrix WORK(IVT) */
01780 /*              Overwrite WORK(IVT) by right singular vectors of L */
01781 /*              (CWorkspace: need M*M+3*M, prefer M*M+2*M+M*NB) */
01782 /*              (RWorkspace: 0) */
01783 
01784                 clacp2_("F", m, m, &rwork[irvt], m, &work[ivt], &ldwkvt);
01785                 i__1 = *lwork - nwork + 1;
01786                 cunmbr_("P", "R", "C", m, m, m, &a[a_offset], lda, &work[
01787                         itaup], &work[ivt], &ldwkvt, &work[nwork], &i__1, &
01788                         ierr);
01789 
01790 /*              Multiply right singular vectors of L in WORK(IVT) by */
01791 /*              Q in VT, storing result in A */
01792 /*              (CWorkspace: need M*M) */
01793 /*              (RWorkspace: 0) */
01794 
01795                 cgemm_("N", "N", m, n, m, &c_b2, &work[ivt], &ldwkvt, &vt[
01796                         vt_offset], ldvt, &c_b1, &a[a_offset], lda);
01797 
01798 /*              Copy right singular vectors of A from A to VT */
01799 
01800                 clacpy_("F", m, n, &a[a_offset], lda, &vt[vt_offset], ldvt);
01801 
01802             }
01803 
01804         } else if (*n >= mnthr2) {
01805 
01806 /*           MNTHR2 <= N < MNTHR1 */
01807 
01808 /*           Path 5t (N much larger than M, but not as much as MNTHR1) */
01809 /*           Reduce to bidiagonal form without QR decomposition, use */
01810 /*           CUNGBR and matrix multiplication to compute singular vectors */
01811 
01812 
01813             ie = 1;
01814             nrwork = ie + *m;
01815             itauq = 1;
01816             itaup = itauq + *m;
01817             nwork = itaup + *m;
01818 
01819 /*           Bidiagonalize A */
01820 /*           (CWorkspace: need 2*M+N, prefer 2*M+(M+N)*NB) */
01821 /*           (RWorkspace: M) */
01822 
01823             i__1 = *lwork - nwork + 1;
01824             cgebrd_(m, n, &a[a_offset], lda, &s[1], &rwork[ie], &work[itauq], 
01825                     &work[itaup], &work[nwork], &i__1, &ierr);
01826 
01827             if (wntqn) {
01828 
01829 /*              Compute singular values only */
01830 /*              (Cworkspace: 0) */
01831 /*              (Rworkspace: need BDSPAN) */
01832 
01833                 sbdsdc_("L", "N", m, &s[1], &rwork[ie], dum, &c__1, dum, &
01834                         c__1, dum, idum, &rwork[nrwork], &iwork[1], info);
01835             } else if (wntqo) {
01836                 irvt = nrwork;
01837                 iru = irvt + *m * *m;
01838                 nrwork = iru + *m * *m;
01839                 ivt = nwork;
01840 
01841 /*              Copy A to U, generate Q */
01842 /*              (Cworkspace: need 2*M, prefer M+M*NB) */
01843 /*              (Rworkspace: 0) */
01844 
01845                 clacpy_("L", m, m, &a[a_offset], lda, &u[u_offset], ldu);
01846                 i__1 = *lwork - nwork + 1;
01847                 cungbr_("Q", m, m, n, &u[u_offset], ldu, &work[itauq], &work[
01848                         nwork], &i__1, &ierr);
01849 
01850 /*              Generate P**H in A */
01851 /*              (Cworkspace: need 2*M, prefer M+M*NB) */
01852 /*              (Rworkspace: 0) */
01853 
01854                 i__1 = *lwork - nwork + 1;
01855                 cungbr_("P", m, n, m, &a[a_offset], lda, &work[itaup], &work[
01856                         nwork], &i__1, &ierr);
01857 
01858                 ldwkvt = *m;
01859                 if (*lwork >= *m * *n + *m * 3) {
01860 
01861 /*                 WORK( IVT ) is M by N */
01862 
01863                     nwork = ivt + ldwkvt * *n;
01864                     chunk = *n;
01865                 } else {
01866 
01867 /*                 WORK( IVT ) is M by CHUNK */
01868 
01869                     chunk = (*lwork - *m * 3) / *m;
01870                     nwork = ivt + ldwkvt * chunk;
01871                 }
01872 
01873 /*              Perform bidiagonal SVD, computing left singular vectors */
01874 /*              of bidiagonal matrix in RWORK(IRU) and computing right */
01875 /*              singular vectors of bidiagonal matrix in RWORK(IRVT) */
01876 /*              (CWorkspace: need 0) */
01877 /*              (RWorkspace: need BDSPAC) */
01878 
01879                 sbdsdc_("L", "I", m, &s[1], &rwork[ie], &rwork[iru], m, &
01880                         rwork[irvt], m, dum, idum, &rwork[nrwork], &iwork[1], 
01881                         info);
01882 
01883 /*              Multiply Q in U by real matrix RWORK(IRVT) */
01884 /*              storing the result in WORK(IVT), copying to U */
01885 /*              (Cworkspace: need 0) */
01886 /*              (Rworkspace: need 2*M*M) */
01887 
01888                 clacrm_(m, m, &u[u_offset], ldu, &rwork[iru], m, &work[ivt], &
01889                         ldwkvt, &rwork[nrwork]);
01890                 clacpy_("F", m, m, &work[ivt], &ldwkvt, &u[u_offset], ldu);
01891 
01892 /*              Multiply RWORK(IRVT) by P**H in A, storing the */
01893 /*              result in WORK(IVT), copying to A */
01894 /*              (CWorkspace: need M*M, prefer M*N) */
01895 /*              (Rworkspace: need 2*M*M, prefer 2*M*N) */
01896 
01897                 nrwork = iru;
01898                 i__1 = *n;
01899                 i__2 = chunk;
01900                 for (i__ = 1; i__2 < 0 ? i__ >= i__1 : i__ <= i__1; i__ += 
01901                         i__2) {
01902 /* Computing MIN */
01903                     i__3 = *n - i__ + 1;
01904                     blk = min(i__3,chunk);
01905                     clarcm_(m, &blk, &rwork[irvt], m, &a[i__ * a_dim1 + 1], 
01906                             lda, &work[ivt], &ldwkvt, &rwork[nrwork]);
01907                     clacpy_("F", m, &blk, &work[ivt], &ldwkvt, &a[i__ * 
01908                             a_dim1 + 1], lda);
01909 /* L50: */
01910                 }
01911             } else if (wntqs) {
01912 
01913 /*              Copy A to U, generate Q */
01914 /*              (Cworkspace: need 2*M, prefer M+M*NB) */
01915 /*              (Rworkspace: 0) */
01916 
01917                 clacpy_("L", m, m, &a[a_offset], lda, &u[u_offset], ldu);
01918                 i__2 = *lwork - nwork + 1;
01919                 cungbr_("Q", m, m, n, &u[u_offset], ldu, &work[itauq], &work[
01920                         nwork], &i__2, &ierr);
01921 
01922 /*              Copy A to VT, generate P**H */
01923 /*              (Cworkspace: need 2*M, prefer M+M*NB) */
01924 /*              (Rworkspace: 0) */
01925 
01926                 clacpy_("U", m, n, &a[a_offset], lda, &vt[vt_offset], ldvt);
01927                 i__2 = *lwork - nwork + 1;
01928                 cungbr_("P", m, n, m, &vt[vt_offset], ldvt, &work[itaup], &
01929                         work[nwork], &i__2, &ierr);
01930 
01931 /*              Perform bidiagonal SVD, computing left singular vectors */
01932 /*              of bidiagonal matrix in RWORK(IRU) and computing right */
01933 /*              singular vectors of bidiagonal matrix in RWORK(IRVT) */
01934 /*              (CWorkspace: need 0) */
01935 /*              (RWorkspace: need BDSPAC) */
01936 
01937                 irvt = nrwork;
01938                 iru = irvt + *m * *m;
01939                 nrwork = iru + *m * *m;
01940                 sbdsdc_("L", "I", m, &s[1], &rwork[ie], &rwork[iru], m, &
01941                         rwork[irvt], m, dum, idum, &rwork[nrwork], &iwork[1], 
01942                         info);
01943 
01944 /*              Multiply Q in U by real matrix RWORK(IRU), storing the */
01945 /*              result in A, copying to U */
01946 /*              (CWorkspace: need 0) */
01947 /*              (Rworkspace: need 3*M*M) */
01948 
01949                 clacrm_(m, m, &u[u_offset], ldu, &rwork[iru], m, &a[a_offset], 
01950                          lda, &rwork[nrwork]);
01951                 clacpy_("F", m, m, &a[a_offset], lda, &u[u_offset], ldu);
01952 
01953 /*              Multiply real matrix RWORK(IRVT) by P**H in VT, */
01954 /*              storing the result in A, copying to VT */
01955 /*              (Cworkspace: need 0) */
01956 /*              (Rworkspace: need M*M+2*M*N) */
01957 
01958                 nrwork = iru;
01959                 clarcm_(m, n, &rwork[irvt], m, &vt[vt_offset], ldvt, &a[
01960                         a_offset], lda, &rwork[nrwork]);
01961                 clacpy_("F", m, n, &a[a_offset], lda, &vt[vt_offset], ldvt);
01962             } else {
01963 
01964 /*              Copy A to U, generate Q */
01965 /*              (Cworkspace: need 2*M, prefer M+M*NB) */
01966 /*              (Rworkspace: 0) */
01967 
01968                 clacpy_("L", m, m, &a[a_offset], lda, &u[u_offset], ldu);
01969                 i__2 = *lwork - nwork + 1;
01970                 cungbr_("Q", m, m, n, &u[u_offset], ldu, &work[itauq], &work[
01971                         nwork], &i__2, &ierr);
01972 
01973 /*              Copy A to VT, generate P**H */
01974 /*              (Cworkspace: need 2*M, prefer M+M*NB) */
01975 /*              (Rworkspace: 0) */
01976 
01977                 clacpy_("U", m, n, &a[a_offset], lda, &vt[vt_offset], ldvt);
01978                 i__2 = *lwork - nwork + 1;
01979                 cungbr_("P", n, n, m, &vt[vt_offset], ldvt, &work[itaup], &
01980                         work[nwork], &i__2, &ierr);
01981 
01982 /*              Perform bidiagonal SVD, computing left singular vectors */
01983 /*              of bidiagonal matrix in RWORK(IRU) and computing right */
01984 /*              singular vectors of bidiagonal matrix in RWORK(IRVT) */
01985 /*              (CWorkspace: need 0) */
01986 /*              (RWorkspace: need BDSPAC) */
01987 
01988                 irvt = nrwork;
01989                 iru = irvt + *m * *m;
01990                 nrwork = iru + *m * *m;
01991                 sbdsdc_("L", "I", m, &s[1], &rwork[ie], &rwork[iru], m, &
01992                         rwork[irvt], m, dum, idum, &rwork[nrwork], &iwork[1], 
01993                         info);
01994 
01995 /*              Multiply Q in U by real matrix RWORK(IRU), storing the */
01996 /*              result in A, copying to U */
01997 /*              (CWorkspace: need 0) */
01998 /*              (Rworkspace: need 3*M*M) */
01999 
02000                 clacrm_(m, m, &u[u_offset], ldu, &rwork[iru], m, &a[a_offset], 
02001                          lda, &rwork[nrwork]);
02002                 clacpy_("F", m, m, &a[a_offset], lda, &u[u_offset], ldu);
02003 
02004 /*              Multiply real matrix RWORK(IRVT) by P**H in VT, */
02005 /*              storing the result in A, copying to VT */
02006 /*              (Cworkspace: need 0) */
02007 /*              (Rworkspace: need M*M+2*M*N) */
02008 
02009                 clarcm_(m, n, &rwork[irvt], m, &vt[vt_offset], ldvt, &a[
02010                         a_offset], lda, &rwork[nrwork]);
02011                 clacpy_("F", m, n, &a[a_offset], lda, &vt[vt_offset], ldvt);
02012             }
02013 
02014         } else {
02015 
02016 /*           N .LT. MNTHR2 */
02017 
02018 /*           Path 6t (N greater than M, but not much larger) */
02019 /*           Reduce to bidiagonal form without LQ decomposition */
02020 /*           Use CUNMBR to compute singular vectors */
02021 
02022             ie = 1;
02023             nrwork = ie + *m;
02024             itauq = 1;
02025             itaup = itauq + *m;
02026             nwork = itaup + *m;
02027 
02028 /*           Bidiagonalize A */
02029 /*           (CWorkspace: need 2*M+N, prefer 2*M+(M+N)*NB) */
02030 /*           (RWorkspace: M) */
02031 
02032             i__2 = *lwork - nwork + 1;
02033             cgebrd_(m, n, &a[a_offset], lda, &s[1], &rwork[ie], &work[itauq], 
02034                     &work[itaup], &work[nwork], &i__2, &ierr);
02035             if (wntqn) {
02036 
02037 /*              Compute singular values only */
02038 /*              (Cworkspace: 0) */
02039 /*              (Rworkspace: need BDSPAN) */
02040 
02041                 sbdsdc_("L", "N", m, &s[1], &rwork[ie], dum, &c__1, dum, &
02042                         c__1, dum, idum, &rwork[nrwork], &iwork[1], info);
02043             } else if (wntqo) {
02044                 ldwkvt = *m;
02045                 ivt = nwork;
02046                 if (*lwork >= *m * *n + *m * 3) {
02047 
02048 /*                 WORK( IVT ) is M by N */
02049 
02050                     claset_("F", m, n, &c_b1, &c_b1, &work[ivt], &ldwkvt);
02051                     nwork = ivt + ldwkvt * *n;
02052                 } else {
02053 
02054 /*                 WORK( IVT ) is M by CHUNK */
02055 
02056                     chunk = (*lwork - *m * 3) / *m;
02057                     nwork = ivt + ldwkvt * chunk;
02058                 }
02059 
02060 /*              Perform bidiagonal SVD, computing left singular vectors */
02061 /*              of bidiagonal matrix in RWORK(IRU) and computing right */
02062 /*              singular vectors of bidiagonal matrix in RWORK(IRVT) */
02063 /*              (CWorkspace: need 0) */
02064 /*              (RWorkspace: need BDSPAC) */
02065 
02066                 irvt = nrwork;
02067                 iru = irvt + *m * *m;
02068                 nrwork = iru + *m * *m;
02069                 sbdsdc_("L", "I", m, &s[1], &rwork[ie], &rwork[iru], m, &
02070                         rwork[irvt], m, dum, idum, &rwork[nrwork], &iwork[1], 
02071                         info);
02072 
02073 /*              Copy real matrix RWORK(IRU) to complex matrix U */
02074 /*              Overwrite U by left singular vectors of A */
02075 /*              (Cworkspace: need 2*M, prefer M+M*NB) */
02076 /*              (Rworkspace: need 0) */
02077 
02078                 clacp2_("F", m, m, &rwork[iru], m, &u[u_offset], ldu);
02079                 i__2 = *lwork - nwork + 1;
02080                 cunmbr_("Q", "L", "N", m, m, n, &a[a_offset], lda, &work[
02081                         itauq], &u[u_offset], ldu, &work[nwork], &i__2, &ierr);
02082 
02083                 if (*lwork >= *m * *n + *m * 3) {
02084 
02085 /*              Copy real matrix RWORK(IRVT) to complex matrix WORK(IVT) */
02086 /*              Overwrite WORK(IVT) by right singular vectors of A, */
02087 /*              copying to A */
02088 /*              (Cworkspace: need M*N+2*M, prefer M*N+M+M*NB) */
02089 /*              (Rworkspace: need 0) */
02090 
02091                     clacp2_("F", m, m, &rwork[irvt], m, &work[ivt], &ldwkvt);
02092                     i__2 = *lwork - nwork + 1;
02093                     cunmbr_("P", "R", "C", m, n, m, &a[a_offset], lda, &work[
02094                             itaup], &work[ivt], &ldwkvt, &work[nwork], &i__2, 
02095                             &ierr);
02096                     clacpy_("F", m, n, &work[ivt], &ldwkvt, &a[a_offset], lda);
02097                 } else {
02098 
02099 /*                 Generate P**H in A */
02100 /*                 (Cworkspace: need 2*M, prefer M+M*NB) */
02101 /*                 (Rworkspace: need 0) */
02102 
02103                     i__2 = *lwork - nwork + 1;
02104                     cungbr_("P", m, n, m, &a[a_offset], lda, &work[itaup], &
02105                             work[nwork], &i__2, &ierr);
02106 
02107 /*                 Multiply Q in A by real matrix RWORK(IRU), storing the */
02108 /*                 result in WORK(IU), copying to A */
02109 /*                 (CWorkspace: need M*M, prefer M*N) */
02110 /*                 (Rworkspace: need 3*M*M, prefer M*M+2*M*N) */
02111 
02112                     nrwork = iru;
02113                     i__2 = *n;
02114                     i__1 = chunk;
02115                     for (i__ = 1; i__1 < 0 ? i__ >= i__2 : i__ <= i__2; i__ +=
02116                              i__1) {
02117 /* Computing MIN */
02118                         i__3 = *n - i__ + 1;
02119                         blk = min(i__3,chunk);
02120                         clarcm_(m, &blk, &rwork[irvt], m, &a[i__ * a_dim1 + 1]
02121 , lda, &work[ivt], &ldwkvt, &rwork[nrwork]);
02122                         clacpy_("F", m, &blk, &work[ivt], &ldwkvt, &a[i__ * 
02123                                 a_dim1 + 1], lda);
02124 /* L60: */
02125                     }
02126                 }
02127             } else if (wntqs) {
02128 
02129 /*              Perform bidiagonal SVD, computing left singular vectors */
02130 /*              of bidiagonal matrix in RWORK(IRU) and computing right */
02131 /*              singular vectors of bidiagonal matrix in RWORK(IRVT) */
02132 /*              (CWorkspace: need 0) */
02133 /*              (RWorkspace: need BDSPAC) */
02134 
02135                 irvt = nrwork;
02136                 iru = irvt + *m * *m;
02137                 nrwork = iru + *m * *m;
02138                 sbdsdc_("L", "I", m, &s[1], &rwork[ie], &rwork[iru], m, &
02139                         rwork[irvt], m, dum, idum, &rwork[nrwork], &iwork[1], 
02140                         info);
02141 
02142 /*              Copy real matrix RWORK(IRU) to complex matrix U */
02143 /*              Overwrite U by left singular vectors of A */
02144 /*              (CWorkspace: need 3*M, prefer 2*M+M*NB) */
02145 /*              (RWorkspace: M*M) */
02146 
02147                 clacp2_("F", m, m, &rwork[iru], m, &u[u_offset], ldu);
02148                 i__1 = *lwork - nwork + 1;
02149                 cunmbr_("Q", "L", "N", m, m, n, &a[a_offset], lda, &work[
02150                         itauq], &u[u_offset], ldu, &work[nwork], &i__1, &ierr);
02151 
02152 /*              Copy real matrix RWORK(IRVT) to complex matrix VT */
02153 /*              Overwrite VT by right singular vectors of A */
02154 /*              (CWorkspace: need 3*M, prefer 2*M+M*NB) */
02155 /*              (RWorkspace: M*M) */
02156 
02157                 claset_("F", m, n, &c_b1, &c_b1, &vt[vt_offset], ldvt);
02158                 clacp2_("F", m, m, &rwork[irvt], m, &vt[vt_offset], ldvt);
02159                 i__1 = *lwork - nwork + 1;
02160                 cunmbr_("P", "R", "C", m, n, m, &a[a_offset], lda, &work[
02161                         itaup], &vt[vt_offset], ldvt, &work[nwork], &i__1, &
02162                         ierr);
02163             } else {
02164 
02165 /*              Perform bidiagonal SVD, computing left singular vectors */
02166 /*              of bidiagonal matrix in RWORK(IRU) and computing right */
02167 /*              singular vectors of bidiagonal matrix in RWORK(IRVT) */
02168 /*              (CWorkspace: need 0) */
02169 /*              (RWorkspace: need BDSPAC) */
02170 
02171                 irvt = nrwork;
02172                 iru = irvt + *m * *m;
02173                 nrwork = iru + *m * *m;
02174 
02175                 sbdsdc_("L", "I", m, &s[1], &rwork[ie], &rwork[iru], m, &
02176                         rwork[irvt], m, dum, idum, &rwork[nrwork], &iwork[1], 
02177                         info);
02178 
02179 /*              Copy real matrix RWORK(IRU) to complex matrix U */
02180 /*              Overwrite U by left singular vectors of A */
02181 /*              (CWorkspace: need 3*M, prefer 2*M+M*NB) */
02182 /*              (RWorkspace: M*M) */
02183 
02184                 clacp2_("F", m, m, &rwork[iru], m, &u[u_offset], ldu);
02185                 i__1 = *lwork - nwork + 1;
02186                 cunmbr_("Q", "L", "N", m, m, n, &a[a_offset], lda, &work[
02187                         itauq], &u[u_offset], ldu, &work[nwork], &i__1, &ierr);
02188 
02189 /*              Set all of VT to identity matrix */
02190 
02191                 claset_("F", n, n, &c_b1, &c_b2, &vt[vt_offset], ldvt);
02192 
02193 /*              Copy real matrix RWORK(IRVT) to complex matrix VT */
02194 /*              Overwrite VT by right singular vectors of A */
02195 /*              (CWorkspace: need 2*M+N, prefer 2*M+N*NB) */
02196 /*              (RWorkspace: M*M) */
02197 
02198                 clacp2_("F", m, m, &rwork[irvt], m, &vt[vt_offset], ldvt);
02199                 i__1 = *lwork - nwork + 1;
02200                 cunmbr_("P", "R", "C", n, n, m, &a[a_offset], lda, &work[
02201                         itaup], &vt[vt_offset], ldvt, &work[nwork], &i__1, &
02202                         ierr);
02203             }
02204 
02205         }
02206 
02207     }
02208 
02209 /*     Undo scaling if necessary */
02210 
02211     if (iscl == 1) {
02212         if (anrm > bignum) {
02213             slascl_("G", &c__0, &c__0, &bignum, &anrm, &minmn, &c__1, &s[1], &
02214                     minmn, &ierr);
02215         }
02216         if (*info != 0 && anrm > bignum) {
02217             i__1 = minmn - 1;
02218             slascl_("G", &c__0, &c__0, &bignum, &anrm, &i__1, &c__1, &rwork[
02219                     ie], &minmn, &ierr);
02220         }
02221         if (anrm < smlnum) {
02222             slascl_("G", &c__0, &c__0, &smlnum, &anrm, &minmn, &c__1, &s[1], &
02223                     minmn, &ierr);
02224         }
02225         if (*info != 0 && anrm < smlnum) {
02226             i__1 = minmn - 1;
02227             slascl_("G", &c__0, &c__0, &smlnum, &anrm, &i__1, &c__1, &rwork[
02228                     ie], &minmn, &ierr);
02229         }
02230     }
02231 
02232 /*     Return optimal workspace in WORK(1) */
02233 
02234     work[1].r = (real) maxwrk, work[1].i = 0.f;
02235 
02236     return 0;
02237 
02238 /*     End of CGESDD */
02239 
02240 } /* cgesdd_ */

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