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

swiftnav
Author(s):
autogenerated on Sat Jun 8 2019 18:56:33