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