LAPACK: SRC/cgesvd.f Source File

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