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