0001:       SUBROUTINE CGESVD( JOBU, JOBVT, M, N, A, LDA, S, U, LDU, VT, LDVT,
0002:      $                   WORK, LWORK, RWORK, INFO )
0003: *
0004: *  -- LAPACK driver routine (version 3.2) --
0005: *     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
0006: *     November 2006
0007: *
0008: *     .. Scalar Arguments ..
0009:       CHARACTER          JOBU, JOBVT
0010:       INTEGER            INFO, LDA, LDU, LDVT, LWORK, M, N
0011: *     ..
0012: *     .. Array Arguments ..
0013:       REAL               RWORK( * ), S( * )
0014:       COMPLEX            A( LDA, * ), U( LDU, * ), VT( LDVT, * ),
0015:      $                   WORK( * )
0016: *     ..
0017: *
0018: *  Purpose
0019: *  =======
0020: *
0021: *  CGESVD computes the singular value decomposition (SVD) of a complex
0022: *  M-by-N matrix A, optionally computing the left and/or right singular
0023: *  vectors. The SVD is written
0024: *
0025: *       A = U * SIGMA * conjugate-transpose(V)
0026: *
0027: *  where SIGMA is an M-by-N matrix which is zero except for its
0028: *  min(m,n) diagonal elements, U is an M-by-M unitary matrix, and
0029: *  V is an N-by-N unitary matrix.  The diagonal elements of SIGMA
0030: *  are the singular values of A; they are real and non-negative, and
0031: *  are returned in descending order.  The first min(m,n) columns of
0032: *  U and V are the left and right singular vectors of A.
0033: *
0034: *  Note that the routine returns V**H, not V.
0035: *
0036: *  Arguments
0037: *  =========
0038: *
0039: *  JOBU    (input) CHARACTER*1
0040: *          Specifies options for computing all or part of the matrix U:
0041: *          = 'A':  all M columns of U are returned in array U:
0042: *          = 'S':  the first min(m,n) columns of U (the left singular
0043: *                  vectors) are returned in the array U;
0044: *          = 'O':  the first min(m,n) columns of U (the left singular
0045: *                  vectors) are overwritten on the array A;
0046: *          = 'N':  no columns of U (no left singular vectors) are
0047: *                  computed.
0048: *
0049: *  JOBVT   (input) CHARACTER*1
0050: *          Specifies options for computing all or part of the matrix
0051: *          V**H:
0052: *          = 'A':  all N rows of V**H are returned in the array VT;
0053: *          = 'S':  the first min(m,n) rows of V**H (the right singular
0054: *                  vectors) are returned in the array VT;
0055: *          = 'O':  the first min(m,n) rows of V**H (the right singular
0056: *                  vectors) are overwritten on the array A;
0057: *          = 'N':  no rows of V**H (no right singular vectors) are
0058: *                  computed.
0059: *
0060: *          JOBVT and JOBU cannot both be 'O'.
0061: *
0062: *  M       (input) INTEGER
0063: *          The number of rows of the input matrix A.  M >= 0.
0064: *
0065: *  N       (input) INTEGER
0066: *          The number of columns of the input matrix A.  N >= 0.
0067: *
0068: *  A       (input/output) COMPLEX array, dimension (LDA,N)
0069: *          On entry, the M-by-N matrix A.
0070: *          On exit,
0071: *          if JOBU = 'O',  A is overwritten with the first min(m,n)
0072: *                          columns of U (the left singular vectors,
0073: *                          stored columnwise);
0074: *          if JOBVT = 'O', A is overwritten with the first min(m,n)
0075: *                          rows of V**H (the right singular vectors,
0076: *                          stored rowwise);
0077: *          if JOBU .ne. 'O' and JOBVT .ne. 'O', the contents of A
0078: *                          are destroyed.
0079: *
0080: *  LDA     (input) INTEGER
0081: *          The leading dimension of the array A.  LDA >= max(1,M).
0082: *
0083: *  S       (output) REAL array, dimension (min(M,N))
0084: *          The singular values of A, sorted so that S(i) >= S(i+1).
0085: *
0086: *  U       (output) COMPLEX array, dimension (LDU,UCOL)
0087: *          (LDU,M) if JOBU = 'A' or (LDU,min(M,N)) if JOBU = 'S'.
0088: *          If JOBU = 'A', U contains the M-by-M unitary matrix U;
0089: *          if JOBU = 'S', U contains the first min(m,n) columns of U
0090: *          (the left singular vectors, stored columnwise);
0091: *          if JOBU = 'N' or 'O', U is not referenced.
0092: *
0093: *  LDU     (input) INTEGER
0094: *          The leading dimension of the array U.  LDU >= 1; if
0095: *          JOBU = 'S' or 'A', LDU >= M.
0096: *
0097: *  VT      (output) COMPLEX array, dimension (LDVT,N)
0098: *          If JOBVT = 'A', VT contains the N-by-N unitary matrix
0099: *          V**H;
0100: *          if JOBVT = 'S', VT contains the first min(m,n) rows of
0101: *          V**H (the right singular vectors, stored rowwise);
0102: *          if JOBVT = 'N' or 'O', VT is not referenced.
0103: *
0104: *  LDVT    (input) INTEGER
0105: *          The leading dimension of the array VT.  LDVT >= 1; if
0106: *          JOBVT = 'A', LDVT >= N; if JOBVT = 'S', LDVT >= min(M,N).
0107: *
0108: *  WORK    (workspace/output) COMPLEX array, dimension (MAX(1,LWORK))
0109: *          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
0110: *
0111: *  LWORK   (input) INTEGER
0112: *          The dimension of the array WORK.
0113: *          LWORK >=  MAX(1,2*MIN(M,N)+MAX(M,N)).
0114: *          For good performance, LWORK should generally be larger.
0115: *
0116: *          If LWORK = -1, then a workspace query is assumed; the routine
0117: *          only calculates the optimal size of the WORK array, returns
0118: *          this value as the first entry of the WORK array, and no error
0119: *          message related to LWORK is issued by XERBLA.
0120: *
0121: *  RWORK   (workspace) REAL array, dimension (5*min(M,N))
0122: *          On exit, if INFO > 0, RWORK(1:MIN(M,N)-1) contains the
0123: *          unconverged superdiagonal elements of an upper bidiagonal
0124: *          matrix B whose diagonal is in S (not necessarily sorted).
0125: *          B satisfies A = U * B * VT, so it has the same singular
0126: *          values as A, and singular vectors related by U and VT.
0127: *
0128: *  INFO    (output) INTEGER
0129: *          = 0:  successful exit.
0130: *          < 0:  if INFO = -i, the i-th argument had an illegal value.
0131: *          > 0:  if CBDSQR did not converge, INFO specifies how many
0132: *                superdiagonals of an intermediate bidiagonal form B
0133: *                did not converge to zero. See the description of RWORK
0134: *                above for details.
0135: *
0136: *  =====================================================================
0137: *
0138: *     .. Parameters ..
0139:       COMPLEX            CZERO, CONE
0140:       PARAMETER          ( CZERO = ( 0.0E0, 0.0E0 ),
0141:      $                   CONE = ( 1.0E0, 0.0E0 ) )
0142:       REAL               ZERO, ONE
0143:       PARAMETER          ( ZERO = 0.0E0, ONE = 1.0E0 )
0144: *     ..
0145: *     .. Local Scalars ..
0146:       LOGICAL            LQUERY, WNTUA, WNTUAS, WNTUN, WNTUO, WNTUS,
0147:      $                   WNTVA, WNTVAS, WNTVN, WNTVO, WNTVS
0148:       INTEGER            BLK, CHUNK, I, IE, IERR, IR, IRWORK, ISCL,
0149:      $                   ITAU, ITAUP, ITAUQ, IU, IWORK, LDWRKR, LDWRKU,
0150:      $                   MAXWRK, MINMN, MINWRK, MNTHR, NCU, NCVT, NRU,
0151:      $                   NRVT, WRKBL
0152:       REAL               ANRM, BIGNUM, EPS, SMLNUM
0153: *     ..
0154: *     .. Local Arrays ..
0155:       REAL               DUM( 1 )
0156:       COMPLEX            CDUM( 1 )
0157: *     ..
0158: *     .. External Subroutines ..
0159:       EXTERNAL           CBDSQR, CGEBRD, CGELQF, CGEMM, CGEQRF, CLACPY,
0160:      $                   CLASCL, CLASET, CUNGBR, CUNGLQ, CUNGQR, CUNMBR,
0161:      $                   SLASCL, XERBLA
0162: *     ..
0163: *     .. External Functions ..
0164:       LOGICAL            LSAME
0165:       INTEGER            ILAENV
0166:       REAL               CLANGE, SLAMCH
0167:       EXTERNAL           LSAME, ILAENV, CLANGE, SLAMCH
0168: *     ..
0169: *     .. Intrinsic Functions ..
0170:       INTRINSIC          MAX, MIN, SQRT
0171: *     ..
0172: *     .. Executable Statements ..
0173: *
0174: *     Test the input arguments
0175: *
0176:       INFO = 0
0177:       MINMN = MIN( M, N )
0178:       WNTUA = LSAME( JOBU, 'A' )
0179:       WNTUS = LSAME( JOBU, 'S' )
0180:       WNTUAS = WNTUA .OR. WNTUS
0181:       WNTUO = LSAME( JOBU, 'O' )
0182:       WNTUN = LSAME( JOBU, 'N' )
0183:       WNTVA = LSAME( JOBVT, 'A' )
0184:       WNTVS = LSAME( JOBVT, 'S' )
0185:       WNTVAS = WNTVA .OR. WNTVS
0186:       WNTVO = LSAME( JOBVT, 'O' )
0187:       WNTVN = LSAME( JOBVT, 'N' )
0188:       LQUERY = ( LWORK.EQ.-1 )
0189: *
0190:       IF( .NOT.( WNTUA .OR. WNTUS .OR. WNTUO .OR. WNTUN ) ) THEN
0191:          INFO = -1
0192:       ELSE IF( .NOT.( WNTVA .OR. WNTVS .OR. WNTVO .OR. WNTVN ) .OR.
0193:      $         ( WNTVO .AND. WNTUO ) ) THEN
0194:          INFO = -2
0195:       ELSE IF( M.LT.0 ) THEN
0196:          INFO = -3
0197:       ELSE IF( N.LT.0 ) THEN
0198:          INFO = -4
0199:       ELSE IF( LDA.LT.MAX( 1, M ) ) THEN
0200:          INFO = -6
0201:       ELSE IF( LDU.LT.1 .OR. ( WNTUAS .AND. LDU.LT.M ) ) THEN
0202:          INFO = -9
0203:       ELSE IF( LDVT.LT.1 .OR. ( WNTVA .AND. LDVT.LT.N ) .OR.
0204:      $         ( WNTVS .AND. LDVT.LT.MINMN ) ) THEN
0205:          INFO = -11
0206:       END IF
0207: *
0208: *     Compute workspace
0209: *      (Note: Comments in the code beginning "Workspace:" describe the
0210: *       minimal amount of workspace needed at that point in the code,
0211: *       as well as the preferred amount for good performance.
0212: *       CWorkspace refers to complex workspace, and RWorkspace to
0213: *       real workspace. NB refers to the optimal block size for the
0214: *       immediately following subroutine, as returned by ILAENV.)
0215: *
0216:       IF( INFO.EQ.0 ) THEN
0217:          MINWRK = 1
0218:          MAXWRK = 1
0219:          IF( M.GE.N .AND. MINMN.GT.0 ) THEN
0220: *
0221: *           Space needed for CBDSQR is BDSPAC = 5*N
0222: *
0223:             MNTHR = ILAENV( 6, 'CGESVD', JOBU // JOBVT, M, N, 0, 0 )
0224:             IF( M.GE.MNTHR ) THEN
0225:                IF( WNTUN ) THEN
0226: *
0227: *                 Path 1 (M much larger than N, JOBU='N')
0228: *
0229:                   MAXWRK = N + N*ILAENV( 1, 'CGEQRF', ' ', M, N, -1,
0230:      $                     -1 )
0231:                   MAXWRK = MAX( MAXWRK, 2*N+2*N*
0232:      $                     ILAENV( 1, 'CGEBRD', ' ', N, N, -1, -1 ) )
0233:                   IF( WNTVO .OR. WNTVAS )
0234:      $               MAXWRK = MAX( MAXWRK, 2*N+( N-1 )*
0235:      $                        ILAENV( 1, 'CUNGBR', 'P', N, N, N, -1 ) )
0236:                   MINWRK = 3*N
0237:                ELSE IF( WNTUO .AND. WNTVN ) THEN
0238: *
0239: *                 Path 2 (M much larger than N, JOBU='O', JOBVT='N')
0240: *
0241:                   WRKBL = N + N*ILAENV( 1, 'CGEQRF', ' ', M, N, -1, -1 )
0242:                   WRKBL = MAX( WRKBL, N+N*ILAENV( 1, 'CUNGQR', ' ', M,
0243:      $                    N, N, -1 ) )
0244:                   WRKBL = MAX( WRKBL, 2*N+2*N*
0245:      $                    ILAENV( 1, 'CGEBRD', ' ', N, N, -1, -1 ) )
0246:                   WRKBL = MAX( WRKBL, 2*N+N*
0247:      $                    ILAENV( 1, 'CUNGBR', 'Q', N, N, N, -1 ) )
0248:                   MAXWRK = MAX( N*N+WRKBL, N*N+M*N )
0249:                   MINWRK = 2*N + M
0250:                ELSE IF( WNTUO .AND. WNTVAS ) THEN
0251: *
0252: *                 Path 3 (M much larger than N, JOBU='O', JOBVT='S' or
0253: *                 'A')
0254: *
0255:                   WRKBL = N + N*ILAENV( 1, 'CGEQRF', ' ', M, N, -1, -1 )
0256:                   WRKBL = MAX( WRKBL, N+N*ILAENV( 1, 'CUNGQR', ' ', M,
0257:      $                    N, N, -1 ) )
0258:                   WRKBL = MAX( WRKBL, 2*N+2*N*
0259:      $                    ILAENV( 1, 'CGEBRD', ' ', N, N, -1, -1 ) )
0260:                   WRKBL = MAX( WRKBL, 2*N+N*
0261:      $                    ILAENV( 1, 'CUNGBR', 'Q', N, N, N, -1 ) )
0262:                   WRKBL = MAX( WRKBL, 2*N+( N-1 )*
0263:      $                    ILAENV( 1, 'CUNGBR', 'P', N, N, N, -1 ) )
0264:                   MAXWRK = MAX( N*N+WRKBL, N*N+M*N )
0265:                   MINWRK = 2*N + M
0266:                ELSE IF( WNTUS .AND. WNTVN ) THEN
0267: *
0268: *                 Path 4 (M much larger than N, JOBU='S', JOBVT='N')
0269: *
0270:                   WRKBL = N + N*ILAENV( 1, 'CGEQRF', ' ', M, N, -1, -1 )
0271:                   WRKBL = MAX( WRKBL, N+N*ILAENV( 1, 'CUNGQR', ' ', M,
0272:      $                    N, N, -1 ) )
0273:                   WRKBL = MAX( WRKBL, 2*N+2*N*
0274:      $                    ILAENV( 1, 'CGEBRD', ' ', N, N, -1, -1 ) )
0275:                   WRKBL = MAX( WRKBL, 2*N+N*
0276:      $                    ILAENV( 1, 'CUNGBR', 'Q', N, N, N, -1 ) )
0277:                   MAXWRK = N*N + WRKBL
0278:                   MINWRK = 2*N + M
0279:                ELSE IF( WNTUS .AND. WNTVO ) THEN
0280: *
0281: *                 Path 5 (M much larger than N, JOBU='S', JOBVT='O')
0282: *
0283:                   WRKBL = N + N*ILAENV( 1, 'CGEQRF', ' ', M, N, -1, -1 )
0284:                   WRKBL = MAX( WRKBL, N+N*ILAENV( 1, 'CUNGQR', ' ', M,
0285:      $                    N, N, -1 ) )
0286:                   WRKBL = MAX( WRKBL, 2*N+2*N*
0287:      $                    ILAENV( 1, 'CGEBRD', ' ', N, N, -1, -1 ) )
0288:                   WRKBL = MAX( WRKBL, 2*N+N*
0289:      $                    ILAENV( 1, 'CUNGBR', 'Q', N, N, N, -1 ) )
0290:                   WRKBL = MAX( WRKBL, 2*N+( N-1 )*
0291:      $                    ILAENV( 1, 'CUNGBR', 'P', N, N, N, -1 ) )
0292:                   MAXWRK = 2*N*N + WRKBL
0293:                   MINWRK = 2*N + M
0294:                ELSE IF( WNTUS .AND. WNTVAS ) THEN
0295: *
0296: *                 Path 6 (M much larger than N, JOBU='S', JOBVT='S' or
0297: *                 'A')
0298: *
0299:                   WRKBL = N + N*ILAENV( 1, 'CGEQRF', ' ', M, N, -1, -1 )
0300:                   WRKBL = MAX( WRKBL, N+N*ILAENV( 1, 'CUNGQR', ' ', M,
0301:      $                    N, N, -1 ) )
0302:                   WRKBL = MAX( WRKBL, 2*N+2*N*
0303:      $                    ILAENV( 1, 'CGEBRD', ' ', N, N, -1, -1 ) )
0304:                   WRKBL = MAX( WRKBL, 2*N+N*
0305:      $                    ILAENV( 1, 'CUNGBR', 'Q', N, N, N, -1 ) )
0306:                   WRKBL = MAX( WRKBL, 2*N+( N-1 )*
0307:      $                    ILAENV( 1, 'CUNGBR', 'P', N, N, N, -1 ) )
0308:                   MAXWRK = N*N + WRKBL
0309:                   MINWRK = 2*N + M
0310:                ELSE IF( WNTUA .AND. WNTVN ) THEN
0311: *
0312: *                 Path 7 (M much larger than N, JOBU='A', JOBVT='N')
0313: *
0314:                   WRKBL = N + N*ILAENV( 1, 'CGEQRF', ' ', M, N, -1, -1 )
0315:                   WRKBL = MAX( WRKBL, N+M*ILAENV( 1, 'CUNGQR', ' ', M,
0316:      $                    M, N, -1 ) )
0317:                   WRKBL = MAX( WRKBL, 2*N+2*N*
0318:      $                    ILAENV( 1, 'CGEBRD', ' ', N, N, -1, -1 ) )
0319:                   WRKBL = MAX( WRKBL, 2*N+N*
0320:      $                    ILAENV( 1, 'CUNGBR', 'Q', N, N, N, -1 ) )
0321:                   MAXWRK = N*N + WRKBL
0322:                   MINWRK = 2*N + M
0323:                ELSE IF( WNTUA .AND. WNTVO ) THEN
0324: *
0325: *                 Path 8 (M much larger than N, JOBU='A', JOBVT='O')
0326: *
0327:                   WRKBL = N + N*ILAENV( 1, 'CGEQRF', ' ', M, N, -1, -1 )
0328:                   WRKBL = MAX( WRKBL, N+M*ILAENV( 1, 'CUNGQR', ' ', M,
0329:      $                    M, N, -1 ) )
0330:                   WRKBL = MAX( WRKBL, 2*N+2*N*
0331:      $                    ILAENV( 1, 'CGEBRD', ' ', N, N, -1, -1 ) )
0332:                   WRKBL = MAX( WRKBL, 2*N+N*
0333:      $                    ILAENV( 1, 'CUNGBR', 'Q', N, N, N, -1 ) )
0334:                   WRKBL = MAX( WRKBL, 2*N+( N-1 )*
0335:      $                    ILAENV( 1, 'CUNGBR', 'P', N, N, N, -1 ) )
0336:                   MAXWRK = 2*N*N + WRKBL
0337:                   MINWRK = 2*N + M
0338:                ELSE IF( WNTUA .AND. WNTVAS ) THEN
0339: *
0340: *                 Path 9 (M much larger than N, JOBU='A', JOBVT='S' or
0341: *                 'A')
0342: *
0343:                   WRKBL = N + N*ILAENV( 1, 'CGEQRF', ' ', M, N, -1, -1 )
0344:                   WRKBL = MAX( WRKBL, N+M*ILAENV( 1, 'CUNGQR', ' ', M,
0345:      $                    M, N, -1 ) )
0346:                   WRKBL = MAX( WRKBL, 2*N+2*N*
0347:      $                    ILAENV( 1, 'CGEBRD', ' ', N, N, -1, -1 ) )
0348:                   WRKBL = MAX( WRKBL, 2*N+N*
0349:      $                    ILAENV( 1, 'CUNGBR', 'Q', N, N, N, -1 ) )
0350:                   WRKBL = MAX( WRKBL, 2*N+( N-1 )*
0351:      $                    ILAENV( 1, 'CUNGBR', 'P', N, N, N, -1 ) )
0352:                   MAXWRK = N*N + WRKBL
0353:                   MINWRK = 2*N + M
0354:                END IF
0355:             ELSE
0356: *
0357: *              Path 10 (M at least N, but not much larger)
0358: *
0359:                MAXWRK = 2*N + ( M+N )*ILAENV( 1, 'CGEBRD', ' ', M, N,
0360:      $                  -1, -1 )
0361:                IF( WNTUS .OR. WNTUO )
0362:      $            MAXWRK = MAX( MAXWRK, 2*N+N*
0363:      $                     ILAENV( 1, 'CUNGBR', 'Q', M, N, N, -1 ) )
0364:                IF( WNTUA )
0365:      $            MAXWRK = MAX( MAXWRK, 2*N+M*
0366:      $                     ILAENV( 1, 'CUNGBR', 'Q', M, M, N, -1 ) )
0367:                IF( .NOT.WNTVN )
0368:      $            MAXWRK = MAX( MAXWRK, 2*N+( N-1 )*
0369:      $                     ILAENV( 1, 'CUNGBR', 'P', N, N, N, -1 ) )
0370:                MINWRK = 2*N + M
0371:             END IF
0372:          ELSE IF( MINMN.GT.0 ) THEN
0373: *
0374: *           Space needed for CBDSQR is BDSPAC = 5*M
0375: *
0376:             MNTHR = ILAENV( 6, 'CGESVD', JOBU // JOBVT, M, N, 0, 0 )
0377:             IF( N.GE.MNTHR ) THEN
0378:                IF( WNTVN ) THEN
0379: *
0380: *                 Path 1t(N much larger than M, JOBVT='N')
0381: *
0382:                   MAXWRK = M + M*ILAENV( 1, 'CGELQF', ' ', M, N, -1,
0383:      $                     -1 )
0384:                   MAXWRK = MAX( MAXWRK, 2*M+2*M*
0385:      $                     ILAENV( 1, 'CGEBRD', ' ', M, M, -1, -1 ) )
0386:                   IF( WNTUO .OR. WNTUAS )
0387:      $               MAXWRK = MAX( MAXWRK, 2*M+M*
0388:      $                        ILAENV( 1, 'CUNGBR', 'Q', M, M, M, -1 ) )
0389:                   MINWRK = 3*M
0390:                ELSE IF( WNTVO .AND. WNTUN ) THEN
0391: *
0392: *                 Path 2t(N much larger than M, JOBU='N', JOBVT='O')
0393: *
0394:                   WRKBL = M + M*ILAENV( 1, 'CGELQF', ' ', M, N, -1, -1 )
0395:                   WRKBL = MAX( WRKBL, M+M*ILAENV( 1, 'CUNGLQ', ' ', M,
0396:      $                    N, M, -1 ) )
0397:                   WRKBL = MAX( WRKBL, 2*M+2*M*
0398:      $                    ILAENV( 1, 'CGEBRD', ' ', M, M, -1, -1 ) )
0399:                   WRKBL = MAX( WRKBL, 2*M+( M-1 )*
0400:      $                    ILAENV( 1, 'CUNGBR', 'P', M, M, M, -1 ) )
0401:                   MAXWRK = MAX( M*M+WRKBL, M*M+M*N )
0402:                   MINWRK = 2*M + N
0403:                ELSE IF( WNTVO .AND. WNTUAS ) THEN
0404: *
0405: *                 Path 3t(N much larger than M, JOBU='S' or 'A',
0406: *                 JOBVT='O')
0407: *
0408:                   WRKBL = M + M*ILAENV( 1, 'CGELQF', ' ', M, N, -1, -1 )
0409:                   WRKBL = MAX( WRKBL, M+M*ILAENV( 1, 'CUNGLQ', ' ', M,
0410:      $                    N, M, -1 ) )
0411:                   WRKBL = MAX( WRKBL, 2*M+2*M*
0412:      $                    ILAENV( 1, 'CGEBRD', ' ', M, M, -1, -1 ) )
0413:                   WRKBL = MAX( WRKBL, 2*M+( M-1 )*
0414:      $                    ILAENV( 1, 'CUNGBR', 'P', M, M, M, -1 ) )
0415:                   WRKBL = MAX( WRKBL, 2*M+M*
0416:      $                    ILAENV( 1, 'CUNGBR', 'Q', M, M, M, -1 ) )
0417:                   MAXWRK = MAX( M*M+WRKBL, M*M+M*N )
0418:                   MINWRK = 2*M + N
0419:                ELSE IF( WNTVS .AND. WNTUN ) THEN
0420: *
0421: *                 Path 4t(N much larger than M, JOBU='N', JOBVT='S')
0422: *
0423:                   WRKBL = M + M*ILAENV( 1, 'CGELQF', ' ', M, N, -1, -1 )
0424:                   WRKBL = MAX( WRKBL, M+M*ILAENV( 1, 'CUNGLQ', ' ', M,
0425:      $                    N, M, -1 ) )
0426:                   WRKBL = MAX( WRKBL, 2*M+2*M*
0427:      $                    ILAENV( 1, 'CGEBRD', ' ', M, M, -1, -1 ) )
0428:                   WRKBL = MAX( WRKBL, 2*M+( M-1 )*
0429:      $                    ILAENV( 1, 'CUNGBR', 'P', M, M, M, -1 ) )
0430:                   MAXWRK = M*M + WRKBL
0431:                   MINWRK = 2*M + N
0432:                ELSE IF( WNTVS .AND. WNTUO ) THEN
0433: *
0434: *                 Path 5t(N much larger than M, JOBU='O', JOBVT='S')
0435: *
0436:                   WRKBL = M + M*ILAENV( 1, 'CGELQF', ' ', M, N, -1, -1 )
0437:                   WRKBL = MAX( WRKBL, M+M*ILAENV( 1, 'CUNGLQ', ' ', M,
0438:      $                    N, M, -1 ) )
0439:                   WRKBL = MAX( WRKBL, 2*M+2*M*
0440:      $                    ILAENV( 1, 'CGEBRD', ' ', M, M, -1, -1 ) )
0441:                   WRKBL = MAX( WRKBL, 2*M+( M-1 )*
0442:      $                    ILAENV( 1, 'CUNGBR', 'P', M, M, M, -1 ) )
0443:                   WRKBL = MAX( WRKBL, 2*M+M*
0444:      $                    ILAENV( 1, 'CUNGBR', 'Q', M, M, M, -1 ) )
0445:                   MAXWRK = 2*M*M + WRKBL
0446:                   MINWRK = 2*M + N
0447:                ELSE IF( WNTVS .AND. WNTUAS ) THEN
0448: *
0449: *                 Path 6t(N much larger than M, JOBU='S' or 'A',
0450: *                 JOBVT='S')
0451: *
0452:                   WRKBL = M + M*ILAENV( 1, 'CGELQF', ' ', M, N, -1, -1 )
0453:                   WRKBL = MAX( WRKBL, M+M*ILAENV( 1, 'CUNGLQ', ' ', M,
0454:      $                    N, M, -1 ) )
0455:                   WRKBL = MAX( WRKBL, 2*M+2*M*
0456:      $                    ILAENV( 1, 'CGEBRD', ' ', M, M, -1, -1 ) )
0457:                   WRKBL = MAX( WRKBL, 2*M+( M-1 )*
0458:      $                    ILAENV( 1, 'CUNGBR', 'P', M, M, M, -1 ) )
0459:                   WRKBL = MAX( WRKBL, 2*M+M*
0460:      $                    ILAENV( 1, 'CUNGBR', 'Q', M, M, M, -1 ) )
0461:                   MAXWRK = M*M + WRKBL
0462:                   MINWRK = 2*M + N
0463:                ELSE IF( WNTVA .AND. WNTUN ) THEN
0464: *
0465: *                 Path 7t(N much larger than M, JOBU='N', JOBVT='A')
0466: *
0467:                   WRKBL = M + M*ILAENV( 1, 'CGELQF', ' ', M, N, -1, -1 )
0468:                   WRKBL = MAX( WRKBL, M+N*ILAENV( 1, 'CUNGLQ', ' ', N,
0469:      $                    N, M, -1 ) )
0470:                   WRKBL = MAX( WRKBL, 2*M+2*M*
0471:      $                    ILAENV( 1, 'CGEBRD', ' ', M, M, -1, -1 ) )
0472:                   WRKBL = MAX( WRKBL, 2*M+( M-1 )*
0473:      $                    ILAENV( 1, 'CUNGBR', 'P', M, M, M, -1 ) )
0474:                   MAXWRK = M*M + WRKBL
0475:                   MINWRK = 2*M + N
0476:                ELSE IF( WNTVA .AND. WNTUO ) THEN
0477: *
0478: *                 Path 8t(N much larger than M, JOBU='O', JOBVT='A')
0479: *
0480:                   WRKBL = M + M*ILAENV( 1, 'CGELQF', ' ', M, N, -1, -1 )
0481:                   WRKBL = MAX( WRKBL, M+N*ILAENV( 1, 'CUNGLQ', ' ', N,
0482:      $                    N, M, -1 ) )
0483:                   WRKBL = MAX( WRKBL, 2*M+2*M*
0484:      $                    ILAENV( 1, 'CGEBRD', ' ', M, M, -1, -1 ) )
0485:                   WRKBL = MAX( WRKBL, 2*M+( M-1 )*
0486:      $                    ILAENV( 1, 'CUNGBR', 'P', M, M, M, -1 ) )
0487:                   WRKBL = MAX( WRKBL, 2*M+M*
0488:      $                    ILAENV( 1, 'CUNGBR', 'Q', M, M, M, -1 ) )
0489:                   MAXWRK = 2*M*M + WRKBL
0490:                   MINWRK = 2*M + N
0491:                ELSE IF( WNTVA .AND. WNTUAS ) THEN
0492: *
0493: *                 Path 9t(N much larger than M, JOBU='S' or 'A',
0494: *                 JOBVT='A')
0495: *
0496:                   WRKBL = M + M*ILAENV( 1, 'CGELQF', ' ', M, N, -1, -1 )
0497:                   WRKBL = MAX( WRKBL, M+N*ILAENV( 1, 'CUNGLQ', ' ', N,
0498:      $                    N, M, -1 ) )
0499:                   WRKBL = MAX( WRKBL, 2*M+2*M*
0500:      $                    ILAENV( 1, 'CGEBRD', ' ', M, M, -1, -1 ) )
0501:                   WRKBL = MAX( WRKBL, 2*M+( M-1 )*
0502:      $                    ILAENV( 1, 'CUNGBR', 'P', M, M, M, -1 ) )
0503:                   WRKBL = MAX( WRKBL, 2*M+M*
0504:      $                    ILAENV( 1, 'CUNGBR', 'Q', M, M, M, -1 ) )
0505:                   MAXWRK = M*M + WRKBL
0506:                   MINWRK = 2*M + N
0507:                END IF
0508:             ELSE
0509: *
0510: *              Path 10t(N greater than M, but not much larger)
0511: *
0512:                MAXWRK = 2*M + ( M+N )*ILAENV( 1, 'CGEBRD', ' ', M, N,
0513:      $                  -1, -1 )
0514:                IF( WNTVS .OR. WNTVO )
0515:      $            MAXWRK = MAX( MAXWRK, 2*M+M*
0516:      $                     ILAENV( 1, 'CUNGBR', 'P', M, N, M, -1 ) )
0517:                IF( WNTVA )
0518:      $            MAXWRK = MAX( MAXWRK, 2*M+N*
0519:      $                     ILAENV( 1, 'CUNGBR', 'P', N, N, M, -1 ) )
0520:                IF( .NOT.WNTUN )
0521:      $            MAXWRK = MAX( MAXWRK, 2*M+( M-1 )*
0522:      $                     ILAENV( 1, 'CUNGBR', 'Q', M, M, M, -1 ) )
0523:                MINWRK = 2*M + N
0524:             END IF
0525:          END IF
0526:          MAXWRK = MAX( MINWRK, MAXWRK )
0527:          WORK( 1 ) = MAXWRK
0528: *
0529:          IF( LWORK.LT.MINWRK .AND. .NOT.LQUERY ) THEN
0530:             INFO = -13
0531:          END IF
0532:       END IF
0533: *
0534:       IF( INFO.NE.0 ) THEN
0535:          CALL XERBLA( 'CGESVD', -INFO )
0536:          RETURN
0537:       ELSE IF( LQUERY ) THEN
0538:          RETURN
0539:       END IF
0540: *
0541: *     Quick return if possible
0542: *
0543:       IF( M.EQ.0 .OR. N.EQ.0 ) THEN
0544:          RETURN
0545:       END IF
0546: *
0547: *     Get machine constants
0548: *
0549:       EPS = SLAMCH( 'P' )
0550:       SMLNUM = SQRT( SLAMCH( 'S' ) ) / EPS
0551:       BIGNUM = ONE / SMLNUM
0552: *
0553: *     Scale A if max element outside range [SMLNUM,BIGNUM]
0554: *
0555:       ANRM = CLANGE( 'M', M, N, A, LDA, DUM )
0556:       ISCL = 0
0557:       IF( ANRM.GT.ZERO .AND. ANRM.LT.SMLNUM ) THEN
0558:          ISCL = 1
0559:          CALL CLASCL( 'G', 0, 0, ANRM, SMLNUM, M, N, A, LDA, IERR )
0560:       ELSE IF( ANRM.GT.BIGNUM ) THEN
0561:          ISCL = 1
0562:          CALL CLASCL( 'G', 0, 0, ANRM, BIGNUM, M, N, A, LDA, IERR )
0563:       END IF
0564: *
0565:       IF( M.GE.N ) THEN
0566: *
0567: *        A has at least as many rows as columns. If A has sufficiently
0568: *        more rows than columns, first reduce using the QR
0569: *        decomposition (if sufficient workspace available)
0570: *
0571:          IF( M.GE.MNTHR ) THEN
0572: *
0573:             IF( WNTUN ) THEN
0574: *
0575: *              Path 1 (M much larger than N, JOBU='N')
0576: *              No left singular vectors to be computed
0577: *
0578:                ITAU = 1
0579:                IWORK = ITAU + N
0580: *
0581: *              Compute A=Q*R
0582: *              (CWorkspace: need 2*N, prefer N+N*NB)
0583: *              (RWorkspace: need 0)
0584: *
0585:                CALL CGEQRF( M, N, A, LDA, WORK( ITAU ), WORK( IWORK ),
0586:      $                      LWORK-IWORK+1, IERR )
0587: *
0588: *              Zero out below R
0589: *
0590:                CALL CLASET( 'L', N-1, N-1, CZERO, CZERO, A( 2, 1 ),
0591:      $                      LDA )
0592:                IE = 1
0593:                ITAUQ = 1
0594:                ITAUP = ITAUQ + N
0595:                IWORK = ITAUP + N
0596: *
0597: *              Bidiagonalize R in A
0598: *              (CWorkspace: need 3*N, prefer 2*N+2*N*NB)
0599: *              (RWorkspace: need N)
0600: *
0601:                CALL CGEBRD( N, N, A, LDA, S, RWORK( IE ), WORK( ITAUQ ),
0602:      $                      WORK( ITAUP ), WORK( IWORK ), LWORK-IWORK+1,
0603:      $                      IERR )
0604:                NCVT = 0
0605:                IF( WNTVO .OR. WNTVAS ) THEN
0606: *
0607: *                 If right singular vectors desired, generate P'.
0608: *                 (CWorkspace: need 3*N-1, prefer 2*N+(N-1)*NB)
0609: *                 (RWorkspace: 0)
0610: *
0611:                   CALL CUNGBR( 'P', N, N, N, A, LDA, WORK( ITAUP ),
0612:      $                         WORK( IWORK ), LWORK-IWORK+1, IERR )
0613:                   NCVT = N
0614:                END IF
0615:                IRWORK = IE + N
0616: *
0617: *              Perform bidiagonal QR iteration, computing right
0618: *              singular vectors of A in A if desired
0619: *              (CWorkspace: 0)
0620: *              (RWorkspace: need BDSPAC)
0621: *
0622:                CALL CBDSQR( 'U', N, NCVT, 0, 0, S, RWORK( IE ), A, LDA,
0623:      $                      CDUM, 1, CDUM, 1, RWORK( IRWORK ), INFO )
0624: *
0625: *              If right singular vectors desired in VT, copy them there
0626: *
0627:                IF( WNTVAS )
0628:      $            CALL CLACPY( 'F', N, N, A, LDA, VT, LDVT )
0629: *
0630:             ELSE IF( WNTUO .AND. WNTVN ) THEN
0631: *
0632: *              Path 2 (M much larger than N, JOBU='O', JOBVT='N')
0633: *              N left singular vectors to be overwritten on A and
0634: *              no right singular vectors to be computed
0635: *
0636:                IF( LWORK.GE.N*N+3*N ) THEN
0637: *
0638: *                 Sufficient workspace for a fast algorithm
0639: *
0640:                   IR = 1
0641:                   IF( LWORK.GE.MAX( WRKBL, LDA*N )+LDA*N ) THEN
0642: *
0643: *                    WORK(IU) is LDA by N, WORK(IR) is LDA by N
0644: *
0645:                      LDWRKU = LDA
0646:                      LDWRKR = LDA
0647:                   ELSE IF( LWORK.GE.MAX( WRKBL, LDA*N )+N*N ) THEN
0648: *
0649: *                    WORK(IU) is LDA by N, WORK(IR) is N by N
0650: *
0651:                      LDWRKU = LDA
0652:                      LDWRKR = N
0653:                   ELSE
0654: *
0655: *                    WORK(IU) is LDWRKU by N, WORK(IR) is N by N
0656: *
0657:                      LDWRKU = ( LWORK-N*N ) / N
0658:                      LDWRKR = N
0659:                   END IF
0660:                   ITAU = IR + LDWRKR*N
0661:                   IWORK = ITAU + N
0662: *
0663: *                 Compute A=Q*R
0664: *                 (CWorkspace: need N*N+2*N, prefer N*N+N+N*NB)
0665: *                 (RWorkspace: 0)
0666: *
0667:                   CALL CGEQRF( M, N, A, LDA, WORK( ITAU ),
0668:      $                         WORK( IWORK ), LWORK-IWORK+1, IERR )
0669: *
0670: *                 Copy R to WORK(IR) and zero out below it
0671: *
0672:                   CALL CLACPY( 'U', N, N, A, LDA, WORK( IR ), LDWRKR )
0673:                   CALL CLASET( 'L', N-1, N-1, CZERO, CZERO,
0674:      $                         WORK( IR+1 ), LDWRKR )
0675: *
0676: *                 Generate Q in A
0677: *                 (CWorkspace: need N*N+2*N, prefer N*N+N+N*NB)
0678: *                 (RWorkspace: 0)
0679: *
0680:                   CALL CUNGQR( M, N, N, A, LDA, WORK( ITAU ),
0681:      $                         WORK( IWORK ), LWORK-IWORK+1, IERR )
0682:                   IE = 1
0683:                   ITAUQ = ITAU
0684:                   ITAUP = ITAUQ + N
0685:                   IWORK = ITAUP + N
0686: *
0687: *                 Bidiagonalize R in WORK(IR)
0688: *                 (CWorkspace: need N*N+3*N, prefer N*N+2*N+2*N*NB)
0689: *                 (RWorkspace: need N)
0690: *
0691:                   CALL CGEBRD( N, N, WORK( IR ), LDWRKR, S, RWORK( IE ),
0692:      $                         WORK( ITAUQ ), WORK( ITAUP ),
0693:      $                         WORK( IWORK ), LWORK-IWORK+1, IERR )
0694: *
0695: *                 Generate left vectors bidiagonalizing R
0696: *                 (CWorkspace: need N*N+3*N, prefer N*N+2*N+N*NB)
0697: *                 (RWorkspace: need 0)
0698: *
0699:                   CALL CUNGBR( 'Q', N, N, N, WORK( IR ), LDWRKR,
0700:      $                         WORK( ITAUQ ), WORK( IWORK ),
0701:      $                         LWORK-IWORK+1, IERR )
0702:                   IRWORK = IE + N
0703: *
0704: *                 Perform bidiagonal QR iteration, computing left
0705: *                 singular vectors of R in WORK(IR)
0706: *                 (CWorkspace: need N*N)
0707: *                 (RWorkspace: need BDSPAC)
0708: *
0709:                   CALL CBDSQR( 'U', N, 0, N, 0, S, RWORK( IE ), CDUM, 1,
0710:      $                         WORK( IR ), LDWRKR, CDUM, 1,
0711:      $                         RWORK( IRWORK ), INFO )
0712:                   IU = ITAUQ
0713: *
0714: *                 Multiply Q in A by left singular vectors of R in
0715: *                 WORK(IR), storing result in WORK(IU) and copying to A
0716: *                 (CWorkspace: need N*N+N, prefer N*N+M*N)
0717: *                 (RWorkspace: 0)
0718: *
0719:                   DO 10 I = 1, M, LDWRKU
0720:                      CHUNK = MIN( M-I+1, LDWRKU )
0721:                      CALL CGEMM( 'N', 'N', CHUNK, N, N, CONE, A( I, 1 ),
0722:      $                           LDA, WORK( IR ), LDWRKR, CZERO,
0723:      $                           WORK( IU ), LDWRKU )
0724:                      CALL CLACPY( 'F', CHUNK, N, WORK( IU ), LDWRKU,
0725:      $                            A( I, 1 ), LDA )
0726:    10             CONTINUE
0727: *
0728:                ELSE
0729: *
0730: *                 Insufficient workspace for a fast algorithm
0731: *
0732:                   IE = 1
0733:                   ITAUQ = 1
0734:                   ITAUP = ITAUQ + N
0735:                   IWORK = ITAUP + N
0736: *
0737: *                 Bidiagonalize A
0738: *                 (CWorkspace: need 2*N+M, prefer 2*N+(M+N)*NB)
0739: *                 (RWorkspace: N)
0740: *
0741:                   CALL CGEBRD( M, N, A, LDA, S, RWORK( IE ),
0742:      $                         WORK( ITAUQ ), WORK( ITAUP ),
0743:      $                         WORK( IWORK ), LWORK-IWORK+1, IERR )
0744: *
0745: *                 Generate left vectors bidiagonalizing A
0746: *                 (CWorkspace: need 3*N, prefer 2*N+N*NB)
0747: *                 (RWorkspace: 0)
0748: *
0749:                   CALL CUNGBR( 'Q', M, N, N, A, LDA, WORK( ITAUQ ),
0750:      $                         WORK( IWORK ), LWORK-IWORK+1, IERR )
0751:                   IRWORK = IE + N
0752: *
0753: *                 Perform bidiagonal QR iteration, computing left
0754: *                 singular vectors of A in A
0755: *                 (CWorkspace: need 0)
0756: *                 (RWorkspace: need BDSPAC)
0757: *
0758:                   CALL CBDSQR( 'U', N, 0, M, 0, S, RWORK( IE ), CDUM, 1,
0759:      $                         A, LDA, CDUM, 1, RWORK( IRWORK ), INFO )
0760: *
0761:                END IF
0762: *
0763:             ELSE IF( WNTUO .AND. WNTVAS ) THEN
0764: *
0765: *              Path 3 (M much larger than N, JOBU='O', JOBVT='S' or 'A')
0766: *              N left singular vectors to be overwritten on A and
0767: *              N right singular vectors to be computed in VT
0768: *
0769:                IF( LWORK.GE.N*N+3*N ) THEN
0770: *
0771: *                 Sufficient workspace for a fast algorithm
0772: *
0773:                   IR = 1
0774:                   IF( LWORK.GE.MAX( WRKBL, LDA*N )+LDA*N ) THEN
0775: *
0776: *                    WORK(IU) is LDA by N and WORK(IR) is LDA by N
0777: *
0778:                      LDWRKU = LDA
0779:                      LDWRKR = LDA
0780:                   ELSE IF( LWORK.GE.MAX( WRKBL, LDA*N )+N*N ) THEN
0781: *
0782: *                    WORK(IU) is LDA by N and WORK(IR) is N by N
0783: *
0784:                      LDWRKU = LDA
0785:                      LDWRKR = N
0786:                   ELSE
0787: *
0788: *                    WORK(IU) is LDWRKU by N and WORK(IR) is N by N
0789: *
0790:                      LDWRKU = ( LWORK-N*N ) / N
0791:                      LDWRKR = N
0792:                   END IF
0793:                   ITAU = IR + LDWRKR*N
0794:                   IWORK = ITAU + N
0795: *
0796: *                 Compute A=Q*R
0797: *                 (CWorkspace: need N*N+2*N, prefer N*N+N+N*NB)
0798: *                 (RWorkspace: 0)
0799: *
0800:                   CALL CGEQRF( M, N, A, LDA, WORK( ITAU ),
0801:      $                         WORK( IWORK ), LWORK-IWORK+1, IERR )
0802: *
0803: *                 Copy R to VT, zeroing out below it
0804: *
0805:                   CALL CLACPY( 'U', N, N, A, LDA, VT, LDVT )
0806:                   IF( N.GT.1 )
0807:      $               CALL CLASET( 'L', N-1, N-1, CZERO, CZERO,
0808:      $                            VT( 2, 1 ), LDVT )
0809: *
0810: *                 Generate Q in A
0811: *                 (CWorkspace: need N*N+2*N, prefer N*N+N+N*NB)
0812: *                 (RWorkspace: 0)
0813: *
0814:                   CALL CUNGQR( M, N, N, A, LDA, WORK( ITAU ),
0815:      $                         WORK( IWORK ), LWORK-IWORK+1, IERR )
0816:                   IE = 1
0817:                   ITAUQ = ITAU
0818:                   ITAUP = ITAUQ + N
0819:                   IWORK = ITAUP + N
0820: *
0821: *                 Bidiagonalize R in VT, copying result to WORK(IR)
0822: *                 (CWorkspace: need N*N+3*N, prefer N*N+2*N+2*N*NB)
0823: *                 (RWorkspace: need N)
0824: *
0825:                   CALL CGEBRD( N, N, VT, LDVT, S, RWORK( IE ),
0826:      $                         WORK( ITAUQ ), WORK( ITAUP ),
0827:      $                         WORK( IWORK ), LWORK-IWORK+1, IERR )
0828:                   CALL CLACPY( 'L', N, N, VT, LDVT, WORK( IR ), LDWRKR )
0829: *
0830: *                 Generate left vectors bidiagonalizing R in WORK(IR)
0831: *                 (CWorkspace: need N*N+3*N, prefer N*N+2*N+N*NB)
0832: *                 (RWorkspace: 0)
0833: *
0834:                   CALL CUNGBR( 'Q', N, N, N, WORK( IR ), LDWRKR,
0835:      $                         WORK( ITAUQ ), WORK( IWORK ),
0836:      $                         LWORK-IWORK+1, IERR )
0837: *
0838: *                 Generate right vectors bidiagonalizing R in VT
0839: *                 (CWorkspace: need N*N+3*N-1, prefer N*N+2*N+(N-1)*NB)
0840: *                 (RWorkspace: 0)
0841: *
0842:                   CALL CUNGBR( 'P', N, N, N, VT, LDVT, WORK( ITAUP ),
0843:      $                         WORK( IWORK ), LWORK-IWORK+1, IERR )
0844:                   IRWORK = IE + N
0845: *
0846: *                 Perform bidiagonal QR iteration, computing left
0847: *                 singular vectors of R in WORK(IR) and computing right
0848: *                 singular vectors of R in VT
0849: *                 (CWorkspace: need N*N)
0850: *                 (RWorkspace: need BDSPAC)
0851: *
0852:                   CALL CBDSQR( 'U', N, N, N, 0, S, RWORK( IE ), VT,
0853:      $                         LDVT, WORK( IR ), LDWRKR, CDUM, 1,
0854:      $                         RWORK( IRWORK ), INFO )
0855:                   IU = ITAUQ
0856: *
0857: *                 Multiply Q in A by left singular vectors of R in
0858: *                 WORK(IR), storing result in WORK(IU) and copying to A
0859: *                 (CWorkspace: need N*N+N, prefer N*N+M*N)
0860: *                 (RWorkspace: 0)
0861: *
0862:                   DO 20 I = 1, M, LDWRKU
0863:                      CHUNK = MIN( M-I+1, LDWRKU )
0864:                      CALL CGEMM( 'N', 'N', CHUNK, N, N, CONE, A( I, 1 ),
0865:      $                           LDA, WORK( IR ), LDWRKR, CZERO,
0866:      $                           WORK( IU ), LDWRKU )
0867:                      CALL CLACPY( 'F', CHUNK, N, WORK( IU ), LDWRKU,
0868:      $                            A( I, 1 ), LDA )
0869:    20             CONTINUE
0870: *
0871:                ELSE
0872: *
0873: *                 Insufficient workspace for a fast algorithm
0874: *
0875:                   ITAU = 1
0876:                   IWORK = ITAU + N
0877: *
0878: *                 Compute A=Q*R
0879: *                 (CWorkspace: need 2*N, prefer N+N*NB)
0880: *                 (RWorkspace: 0)
0881: *
0882:                   CALL CGEQRF( M, N, A, LDA, WORK( ITAU ),
0883:      $                         WORK( IWORK ), LWORK-IWORK+1, IERR )
0884: *
0885: *                 Copy R to VT, zeroing out below it
0886: *
0887:                   CALL CLACPY( 'U', N, N, A, LDA, VT, LDVT )
0888:                   IF( N.GT.1 )
0889:      $               CALL CLASET( 'L', N-1, N-1, CZERO, CZERO,
0890:      $                            VT( 2, 1 ), LDVT )
0891: *
0892: *                 Generate Q in A
0893: *                 (CWorkspace: need 2*N, prefer N+N*NB)
0894: *                 (RWorkspace: 0)
0895: *
0896:                   CALL CUNGQR( M, N, N, A, LDA, WORK( ITAU ),
0897:      $                         WORK( IWORK ), LWORK-IWORK+1, IERR )
0898:                   IE = 1
0899:                   ITAUQ = ITAU
0900:                   ITAUP = ITAUQ + N
0901:                   IWORK = ITAUP + N
0902: *
0903: *                 Bidiagonalize R in VT
0904: *                 (CWorkspace: need 3*N, prefer 2*N+2*N*NB)
0905: *                 (RWorkspace: N)
0906: *
0907:                   CALL CGEBRD( N, N, VT, LDVT, S, RWORK( IE ),
0908:      $                         WORK( ITAUQ ), WORK( ITAUP ),
0909:      $                         WORK( IWORK ), LWORK-IWORK+1, IERR )
0910: *
0911: *                 Multiply Q in A by left vectors bidiagonalizing R
0912: *                 (CWorkspace: need 2*N+M, prefer 2*N+M*NB)
0913: *                 (RWorkspace: 0)
0914: *
0915:                   CALL CUNMBR( 'Q', 'R', 'N', M, N, N, VT, LDVT,
0916:      $                         WORK( ITAUQ ), A, LDA, WORK( IWORK ),
0917:      $                         LWORK-IWORK+1, IERR )
0918: *
0919: *                 Generate right vectors bidiagonalizing R in VT
0920: *                 (CWorkspace: need 3*N-1, prefer 2*N+(N-1)*NB)
0921: *                 (RWorkspace: 0)
0922: *
0923:                   CALL CUNGBR( 'P', N, N, N, VT, LDVT, WORK( ITAUP ),
0924:      $                         WORK( IWORK ), LWORK-IWORK+1, IERR )
0925:                   IRWORK = IE + N
0926: *
0927: *                 Perform bidiagonal QR iteration, computing left
0928: *                 singular vectors of A in A and computing right
0929: *                 singular vectors of A in VT
0930: *                 (CWorkspace: 0)
0931: *                 (RWorkspace: need BDSPAC)
0932: *
0933:                   CALL CBDSQR( 'U', N, N, M, 0, S, RWORK( IE ), VT,
0934:      $                         LDVT, A, LDA, CDUM, 1, RWORK( IRWORK ),
0935:      $                         INFO )
0936: *
0937:                END IF
0938: *
0939:             ELSE IF( WNTUS ) THEN
0940: *
0941:                IF( WNTVN ) THEN
0942: *
0943: *                 Path 4 (M much larger than N, JOBU='S', JOBVT='N')
0944: *                 N left singular vectors to be computed in U and
0945: *                 no right singular vectors to be computed
0946: *
0947:                   IF( LWORK.GE.N*N+3*N ) THEN
0948: *
0949: *                    Sufficient workspace for a fast algorithm
0950: *
0951:                      IR = 1
0952:                      IF( LWORK.GE.WRKBL+LDA*N ) THEN
0953: *
0954: *                       WORK(IR) is LDA by N
0955: *
0956:                         LDWRKR = LDA
0957:                      ELSE
0958: *
0959: *                       WORK(IR) is N by N
0960: *
0961:                         LDWRKR = N
0962:                      END IF
0963:                      ITAU = IR + LDWRKR*N
0964:                      IWORK = ITAU + N
0965: *
0966: *                    Compute A=Q*R
0967: *                    (CWorkspace: need N*N+2*N, prefer N*N+N+N*NB)
0968: *                    (RWorkspace: 0)
0969: *
0970:                      CALL CGEQRF( M, N, A, LDA, WORK( ITAU ),
0971:      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
0972: *
0973: *                    Copy R to WORK(IR), zeroing out below it
0974: *
0975:                      CALL CLACPY( 'U', N, N, A, LDA, WORK( IR ),
0976:      $                            LDWRKR )
0977:                      CALL CLASET( 'L', N-1, N-1, CZERO, CZERO,
0978:      $                            WORK( IR+1 ), LDWRKR )
0979: *
0980: *                    Generate Q in A
0981: *                    (CWorkspace: need N*N+2*N, prefer N*N+N+N*NB)
0982: *                    (RWorkspace: 0)
0983: *
0984:                      CALL CUNGQR( M, N, N, A, LDA, WORK( ITAU ),
0985:      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
0986:                      IE = 1
0987:                      ITAUQ = ITAU
0988:                      ITAUP = ITAUQ + N
0989:                      IWORK = ITAUP + N
0990: *
0991: *                    Bidiagonalize R in WORK(IR)
0992: *                    (CWorkspace: need N*N+3*N, prefer N*N+2*N+2*N*NB)
0993: *                    (RWorkspace: need N)
0994: *
0995:                      CALL CGEBRD( N, N, WORK( IR ), LDWRKR, S,
0996:      $                            RWORK( IE ), WORK( ITAUQ ),
0997:      $                            WORK( ITAUP ), WORK( IWORK ),
0998:      $                            LWORK-IWORK+1, IERR )
0999: *
1000: *                    Generate left vectors bidiagonalizing R in WORK(IR)
1001: *                    (CWorkspace: need N*N+3*N, prefer N*N+2*N+N*NB)
1002: *                    (RWorkspace: 0)
1003: *
1004:                      CALL CUNGBR( 'Q', N, N, N, WORK( IR ), LDWRKR,
1005:      $                            WORK( ITAUQ ), WORK( IWORK ),
1006:      $                            LWORK-IWORK+1, IERR )
1007:                      IRWORK = IE + N
1008: *
1009: *                    Perform bidiagonal QR iteration, computing left
1010: *                    singular vectors of R in WORK(IR)
1011: *                    (CWorkspace: need N*N)
1012: *                    (RWorkspace: need BDSPAC)
1013: *
1014:                      CALL CBDSQR( 'U', N, 0, N, 0, S, RWORK( IE ), CDUM,
1015:      $                            1, WORK( IR ), LDWRKR, CDUM, 1,
1016:      $                            RWORK( IRWORK ), INFO )
1017: *
1018: *                    Multiply Q in A by left singular vectors of R in
1019: *                    WORK(IR), storing result in U
1020: *                    (CWorkspace: need N*N)
1021: *                    (RWorkspace: 0)
1022: *
1023:                      CALL CGEMM( 'N', 'N', M, N, N, CONE, A, LDA,
1024:      $                           WORK( IR ), LDWRKR, CZERO, U, LDU )
1025: *
1026:                   ELSE
1027: *
1028: *                    Insufficient workspace for a fast algorithm
1029: *
1030:                      ITAU = 1
1031:                      IWORK = ITAU + N
1032: *
1033: *                    Compute A=Q*R, copying result to U
1034: *                    (CWorkspace: need 2*N, prefer N+N*NB)
1035: *                    (RWorkspace: 0)
1036: *
1037:                      CALL CGEQRF( M, N, A, LDA, WORK( ITAU ),
1038:      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
1039:                      CALL CLACPY( 'L', M, N, A, LDA, U, LDU )
1040: *
1041: *                    Generate Q in U
1042: *                    (CWorkspace: need 2*N, prefer N+N*NB)
1043: *                    (RWorkspace: 0)
1044: *
1045:                      CALL CUNGQR( M, N, N, U, LDU, WORK( ITAU ),
1046:      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
1047:                      IE = 1
1048:                      ITAUQ = ITAU
1049:                      ITAUP = ITAUQ + N
1050:                      IWORK = ITAUP + N
1051: *
1052: *                    Zero out below R in A
1053: *
1054:                      CALL CLASET( 'L', N-1, N-1, CZERO, CZERO,
1055:      $                            A( 2, 1 ), LDA )
1056: *
1057: *                    Bidiagonalize R in A
1058: *                    (CWorkspace: need 3*N, prefer 2*N+2*N*NB)
1059: *                    (RWorkspace: need N)
1060: *
1061:                      CALL CGEBRD( N, N, A, LDA, S, RWORK( IE ),
1062:      $                            WORK( ITAUQ ), WORK( ITAUP ),
1063:      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
1064: *
1065: *                    Multiply Q in U by left vectors bidiagonalizing R
1066: *                    (CWorkspace: need 2*N+M, prefer 2*N+M*NB)
1067: *                    (RWorkspace: 0)
1068: *
1069:                      CALL CUNMBR( 'Q', 'R', 'N', M, N, N, A, LDA,
1070:      $                            WORK( ITAUQ ), U, LDU, WORK( IWORK ),
1071:      $                            LWORK-IWORK+1, IERR )
1072:                      IRWORK = IE + N
1073: *
1074: *                    Perform bidiagonal QR iteration, computing left
1075: *                    singular vectors of A in U
1076: *                    (CWorkspace: 0)
1077: *                    (RWorkspace: need BDSPAC)
1078: *
1079:                      CALL CBDSQR( 'U', N, 0, M, 0, S, RWORK( IE ), CDUM,
1080:      $                            1, U, LDU, CDUM, 1, RWORK( IRWORK ),
1081:      $                            INFO )
1082: *
1083:                   END IF
1084: *
1085:                ELSE IF( WNTVO ) THEN
1086: *
1087: *                 Path 5 (M much larger than N, JOBU='S', JOBVT='O')
1088: *                 N left singular vectors to be computed in U and
1089: *                 N right singular vectors to be overwritten on A
1090: *
1091:                   IF( LWORK.GE.2*N*N+3*N ) THEN
1092: *
1093: *                    Sufficient workspace for a fast algorithm
1094: *
1095:                      IU = 1
1096:                      IF( LWORK.GE.WRKBL+2*LDA*N ) THEN
1097: *
1098: *                       WORK(IU) is LDA by N and WORK(IR) is LDA by N
1099: *
1100:                         LDWRKU = LDA
1101:                         IR = IU + LDWRKU*N
1102:                         LDWRKR = LDA
1103:                      ELSE IF( LWORK.GE.WRKBL+( LDA+N )*N ) THEN
1104: *
1105: *                       WORK(IU) is LDA by N and WORK(IR) is N by N
1106: *
1107:                         LDWRKU = LDA
1108:                         IR = IU + LDWRKU*N
1109:                         LDWRKR = N
1110:                      ELSE
1111: *
1112: *                       WORK(IU) is N by N and WORK(IR) is N by N
1113: *
1114:                         LDWRKU = N
1115:                         IR = IU + LDWRKU*N
1116:                         LDWRKR = N
1117:                      END IF
1118:                      ITAU = IR + LDWRKR*N
1119:                      IWORK = ITAU + N
1120: *
1121: *                    Compute A=Q*R
1122: *                    (CWorkspace: need 2*N*N+2*N, prefer 2*N*N+N+N*NB)
1123: *                    (RWorkspace: 0)
1124: *
1125:                      CALL CGEQRF( M, N, A, LDA, WORK( ITAU ),
1126:      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
1127: *
1128: *                    Copy R to WORK(IU), zeroing out below it
1129: *
1130:                      CALL CLACPY( 'U', N, N, A, LDA, WORK( IU ),
1131:      $                            LDWRKU )
1132:                      CALL CLASET( 'L', N-1, N-1, CZERO, CZERO,
1133:      $                            WORK( IU+1 ), LDWRKU )
1134: *
1135: *                    Generate Q in A
1136: *                    (CWorkspace: need 2*N*N+2*N, prefer 2*N*N+N+N*NB)
1137: *                    (RWorkspace: 0)
1138: *
1139:                      CALL CUNGQR( M, N, N, A, LDA, WORK( ITAU ),
1140:      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
1141:                      IE = 1
1142:                      ITAUQ = ITAU
1143:                      ITAUP = ITAUQ + N
1144:                      IWORK = ITAUP + N
1145: *
1146: *                    Bidiagonalize R in WORK(IU), copying result to
1147: *                    WORK(IR)
1148: *                    (CWorkspace: need   2*N*N+3*N,
1149: *                                 prefer 2*N*N+2*N+2*N*NB)
1150: *                    (RWorkspace: need   N)
1151: *
1152:                      CALL CGEBRD( N, N, WORK( IU ), LDWRKU, S,
1153:      $                            RWORK( IE ), WORK( ITAUQ ),
1154:      $                            WORK( ITAUP ), WORK( IWORK ),
1155:      $                            LWORK-IWORK+1, IERR )
1156:                      CALL CLACPY( 'U', N, N, WORK( IU ), LDWRKU,
1157:      $                            WORK( IR ), LDWRKR )
1158: *
1159: *                    Generate left bidiagonalizing vectors in WORK(IU)
1160: *                    (CWorkspace: need 2*N*N+3*N, prefer 2*N*N+2*N+N*NB)
1161: *                    (RWorkspace: 0)
1162: *
1163:                      CALL CUNGBR( 'Q', N, N, N, WORK( IU ), LDWRKU,
1164:      $                            WORK( ITAUQ ), WORK( IWORK ),
1165:      $                            LWORK-IWORK+1, IERR )
1166: *
1167: *                    Generate right bidiagonalizing vectors in WORK(IR)
1168: *                    (CWorkspace: need   2*N*N+3*N-1,
1169: *                                 prefer 2*N*N+2*N+(N-1)*NB)
1170: *                    (RWorkspace: 0)
1171: *
1172:                      CALL CUNGBR( 'P', N, N, N, WORK( IR ), LDWRKR,
1173:      $                            WORK( ITAUP ), WORK( IWORK ),
1174:      $                            LWORK-IWORK+1, IERR )
1175:                      IRWORK = IE + N
1176: *
1177: *                    Perform bidiagonal QR iteration, computing left
1178: *                    singular vectors of R in WORK(IU) and computing
1179: *                    right singular vectors of R in WORK(IR)
1180: *                    (CWorkspace: need 2*N*N)
1181: *                    (RWorkspace: need BDSPAC)
1182: *
1183:                      CALL CBDSQR( 'U', N, N, N, 0, S, RWORK( IE ),
1184:      $                            WORK( IR ), LDWRKR, WORK( IU ),
1185:      $                            LDWRKU, CDUM, 1, RWORK( IRWORK ),
1186:      $                            INFO )
1187: *
1188: *                    Multiply Q in A by left singular vectors of R in
1189: *                    WORK(IU), storing result in U
1190: *                    (CWorkspace: need N*N)
1191: *                    (RWorkspace: 0)
1192: *
1193:                      CALL CGEMM( 'N', 'N', M, N, N, CONE, A, LDA,
1194:      $                           WORK( IU ), LDWRKU, CZERO, U, LDU )
1195: *
1196: *                    Copy right singular vectors of R to A
1197: *                    (CWorkspace: need N*N)
1198: *                    (RWorkspace: 0)
1199: *
1200:                      CALL CLACPY( 'F', N, N, WORK( IR ), LDWRKR, A,
1201:      $                            LDA )
1202: *
1203:                   ELSE
1204: *
1205: *                    Insufficient workspace for a fast algorithm
1206: *
1207:                      ITAU = 1
1208:                      IWORK = ITAU + N
1209: *
1210: *                    Compute A=Q*R, copying result to U
1211: *                    (CWorkspace: need 2*N, prefer N+N*NB)
1212: *                    (RWorkspace: 0)
1213: *
1214:                      CALL CGEQRF( M, N, A, LDA, WORK( ITAU ),
1215:      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
1216:                      CALL CLACPY( 'L', M, N, A, LDA, U, LDU )
1217: *
1218: *                    Generate Q in U
1219: *                    (CWorkspace: need 2*N, prefer N+N*NB)
1220: *                    (RWorkspace: 0)
1221: *
1222:                      CALL CUNGQR( M, N, N, U, LDU, WORK( ITAU ),
1223:      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
1224:                      IE = 1
1225:                      ITAUQ = ITAU
1226:                      ITAUP = ITAUQ + N
1227:                      IWORK = ITAUP + N
1228: *
1229: *                    Zero out below R in A
1230: *
1231:                      CALL CLASET( 'L', N-1, N-1, CZERO, CZERO,
1232:      $                            A( 2, 1 ), LDA )
1233: *
1234: *                    Bidiagonalize R in A
1235: *                    (CWorkspace: need 3*N, prefer 2*N+2*N*NB)
1236: *                    (RWorkspace: need N)
1237: *
1238:                      CALL CGEBRD( N, N, A, LDA, S, RWORK( IE ),
1239:      $                            WORK( ITAUQ ), WORK( ITAUP ),
1240:      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
1241: *
1242: *                    Multiply Q in U by left vectors bidiagonalizing R
1243: *                    (CWorkspace: need 2*N+M, prefer 2*N+M*NB)
1244: *                    (RWorkspace: 0)
1245: *
1246:                      CALL CUNMBR( 'Q', 'R', 'N', M, N, N, A, LDA,
1247:      $                            WORK( ITAUQ ), U, LDU, WORK( IWORK ),
1248:      $                            LWORK-IWORK+1, IERR )
1249: *
1250: *                    Generate right vectors bidiagonalizing R in A
1251: *                    (CWorkspace: need 3*N-1, prefer 2*N+(N-1)*NB)
1252: *                    (RWorkspace: 0)
1253: *
1254:                      CALL CUNGBR( 'P', N, N, N, A, LDA, WORK( ITAUP ),
1255:      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
1256:                      IRWORK = IE + N
1257: *
1258: *                    Perform bidiagonal QR iteration, computing left
1259: *                    singular vectors of A in U and computing right
1260: *                    singular vectors of A in A
1261: *                    (CWorkspace: 0)
1262: *                    (RWorkspace: need BDSPAC)
1263: *
1264:                      CALL CBDSQR( 'U', N, N, M, 0, S, RWORK( IE ), A,
1265:      $                            LDA, U, LDU, CDUM, 1, RWORK( IRWORK ),
1266:      $                            INFO )
1267: *
1268:                   END IF
1269: *
1270:                ELSE IF( WNTVAS ) THEN
1271: *
1272: *                 Path 6 (M much larger than N, JOBU='S', JOBVT='S'
1273: *                         or 'A')
1274: *                 N left singular vectors to be computed in U and
1275: *                 N right singular vectors to be computed in VT
1276: *
1277:                   IF( LWORK.GE.N*N+3*N ) THEN
1278: *
1279: *                    Sufficient workspace for a fast algorithm
1280: *
1281:                      IU = 1
1282:                      IF( LWORK.GE.WRKBL+LDA*N ) THEN
1283: *
1284: *                       WORK(IU) is LDA by N
1285: *
1286:                         LDWRKU = LDA
1287:                      ELSE
1288: *
1289: *                       WORK(IU) is N by N
1290: *
1291:                         LDWRKU = N
1292:                      END IF
1293:                      ITAU = IU + LDWRKU*N
1294:                      IWORK = ITAU + N
1295: *
1296: *                    Compute A=Q*R
1297: *                    (CWorkspace: need N*N+2*N, prefer N*N+N+N*NB)
1298: *                    (RWorkspace: 0)
1299: *
1300:                      CALL CGEQRF( M, N, A, LDA, WORK( ITAU ),
1301:      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
1302: *
1303: *                    Copy R to WORK(IU), zeroing out below it
1304: *
1305:                      CALL CLACPY( 'U', N, N, A, LDA, WORK( IU ),
1306:      $                            LDWRKU )
1307:                      CALL CLASET( 'L', N-1, N-1, CZERO, CZERO,
1308:      $                            WORK( IU+1 ), LDWRKU )
1309: *
1310: *                    Generate Q in A
1311: *                    (CWorkspace: need N*N+2*N, prefer N*N+N+N*NB)
1312: *                    (RWorkspace: 0)
1313: *
1314:                      CALL CUNGQR( M, N, N, A, LDA, WORK( ITAU ),
1315:      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
1316:                      IE = 1
1317:                      ITAUQ = ITAU
1318:                      ITAUP = ITAUQ + N
1319:                      IWORK = ITAUP + N
1320: *
1321: *                    Bidiagonalize R in WORK(IU), copying result to VT
1322: *                    (CWorkspace: need N*N+3*N, prefer N*N+2*N+2*N*NB)
1323: *                    (RWorkspace: need N)
1324: *
1325:                      CALL CGEBRD( N, N, WORK( IU ), LDWRKU, S,
1326:      $                            RWORK( IE ), WORK( ITAUQ ),
1327:      $                            WORK( ITAUP ), WORK( IWORK ),
1328:      $                            LWORK-IWORK+1, IERR )
1329:                      CALL CLACPY( 'U', N, N, WORK( IU ), LDWRKU, VT,
1330:      $                            LDVT )
1331: *
1332: *                    Generate left bidiagonalizing vectors in WORK(IU)
1333: *                    (CWorkspace: need N*N+3*N, prefer N*N+2*N+N*NB)
1334: *                    (RWorkspace: 0)
1335: *
1336:                      CALL CUNGBR( 'Q', N, N, N, WORK( IU ), LDWRKU,
1337:      $                            WORK( ITAUQ ), WORK( IWORK ),
1338:      $                            LWORK-IWORK+1, IERR )
1339: *
1340: *                    Generate right bidiagonalizing vectors in VT
1341: *                    (CWorkspace: need   N*N+3*N-1,
1342: *                                 prefer N*N+2*N+(N-1)*NB)
1343: *                    (RWorkspace: 0)
1344: *
1345:                      CALL CUNGBR( 'P', N, N, N, VT, LDVT, WORK( ITAUP ),
1346:      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
1347:                      IRWORK = IE + N
1348: *
1349: *                    Perform bidiagonal QR iteration, computing left
1350: *                    singular vectors of R in WORK(IU) and computing
1351: *                    right singular vectors of R in VT
1352: *                    (CWorkspace: need N*N)
1353: *                    (RWorkspace: need BDSPAC)
1354: *
1355:                      CALL CBDSQR( 'U', N, N, N, 0, S, RWORK( IE ), VT,
1356:      $                            LDVT, WORK( IU ), LDWRKU, CDUM, 1,
1357:      $                            RWORK( IRWORK ), INFO )
1358: *
1359: *                    Multiply Q in A by left singular vectors of R in
1360: *                    WORK(IU), storing result in U
1361: *                    (CWorkspace: need N*N)
1362: *                    (RWorkspace: 0)
1363: *
1364:                      CALL CGEMM( 'N', 'N', M, N, N, CONE, A, LDA,
1365:      $                           WORK( IU ), LDWRKU, CZERO, U, LDU )
1366: *
1367:                   ELSE
1368: *
1369: *                    Insufficient workspace for a fast algorithm
1370: *
1371:                      ITAU = 1
1372:                      IWORK = ITAU + N
1373: *
1374: *                    Compute A=Q*R, copying result to U
1375: *                    (CWorkspace: need 2*N, prefer N+N*NB)
1376: *                    (RWorkspace: 0)
1377: *
1378:                      CALL CGEQRF( M, N, A, LDA, WORK( ITAU ),
1379:      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
1380:                      CALL CLACPY( 'L', M, N, A, LDA, U, LDU )
1381: *
1382: *                    Generate Q in U
1383: *                    (CWorkspace: need 2*N, prefer N+N*NB)
1384: *                    (RWorkspace: 0)
1385: *
1386:                      CALL CUNGQR( M, N, N, U, LDU, WORK( ITAU ),
1387:      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
1388: *
1389: *                    Copy R to VT, zeroing out below it
1390: *
1391:                      CALL CLACPY( 'U', N, N, A, LDA, VT, LDVT )
1392:                      IF( N.GT.1 )
1393:      $                  CALL CLASET( 'L', N-1, N-1, CZERO, CZERO,
1394:      $                               VT( 2, 1 ), LDVT )
1395:                      IE = 1
1396:                      ITAUQ = ITAU
1397:                      ITAUP = ITAUQ + N
1398:                      IWORK = ITAUP + N
1399: *
1400: *                    Bidiagonalize R in VT
1401: *                    (CWorkspace: need 3*N, prefer 2*N+2*N*NB)
1402: *                    (RWorkspace: need N)
1403: *
1404:                      CALL CGEBRD( N, N, VT, LDVT, S, RWORK( IE ),
1405:      $                            WORK( ITAUQ ), WORK( ITAUP ),
1406:      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
1407: *
1408: *                    Multiply Q in U by left bidiagonalizing vectors
1409: *                    in VT
1410: *                    (CWorkspace: need 2*N+M, prefer 2*N+M*NB)
1411: *                    (RWorkspace: 0)
1412: *
1413:                      CALL CUNMBR( 'Q', 'R', 'N', M, N, N, VT, LDVT,
1414:      $                            WORK( ITAUQ ), U, LDU, WORK( IWORK ),
1415:      $                            LWORK-IWORK+1, IERR )
1416: *
1417: *                    Generate right bidiagonalizing vectors in VT
1418: *                    (CWorkspace: need 3*N-1, prefer 2*N+(N-1)*NB)
1419: *                    (RWorkspace: 0)
1420: *
1421:                      CALL CUNGBR( 'P', N, N, N, VT, LDVT, WORK( ITAUP ),
1422:      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
1423:                      IRWORK = IE + N
1424: *
1425: *                    Perform bidiagonal QR iteration, computing left
1426: *                    singular vectors of A in U and computing right
1427: *                    singular vectors of A in VT
1428: *                    (CWorkspace: 0)
1429: *                    (RWorkspace: need BDSPAC)
1430: *
1431:                      CALL CBDSQR( 'U', N, N, M, 0, S, RWORK( IE ), VT,
1432:      $                            LDVT, U, LDU, CDUM, 1,
1433:      $                            RWORK( IRWORK ), INFO )
1434: *
1435:                   END IF
1436: *
1437:                END IF
1438: *
1439:             ELSE IF( WNTUA ) THEN
1440: *
1441:                IF( WNTVN ) THEN
1442: *
1443: *                 Path 7 (M much larger than N, JOBU='A', JOBVT='N')
1444: *                 M left singular vectors to be computed in U and
1445: *                 no right singular vectors to be computed
1446: *
1447:                   IF( LWORK.GE.N*N+MAX( N+M, 3*N ) ) THEN
1448: *
1449: *                    Sufficient workspace for a fast algorithm
1450: *
1451:                      IR = 1
1452:                      IF( LWORK.GE.WRKBL+LDA*N ) THEN
1453: *
1454: *                       WORK(IR) is LDA by N
1455: *
1456:                         LDWRKR = LDA
1457:                      ELSE
1458: *
1459: *                       WORK(IR) is N by N
1460: *
1461:                         LDWRKR = N
1462:                      END IF
1463:                      ITAU = IR + LDWRKR*N
1464:                      IWORK = ITAU + N
1465: *
1466: *                    Compute A=Q*R, copying result to U
1467: *                    (CWorkspace: need N*N+2*N, prefer N*N+N+N*NB)
1468: *                    (RWorkspace: 0)
1469: *
1470:                      CALL CGEQRF( M, N, A, LDA, WORK( ITAU ),
1471:      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
1472:                      CALL CLACPY( 'L', M, N, A, LDA, U, LDU )
1473: *
1474: *                    Copy R to WORK(IR), zeroing out below it
1475: *
1476:                      CALL CLACPY( 'U', N, N, A, LDA, WORK( IR ),
1477:      $                            LDWRKR )
1478:                      CALL CLASET( 'L', N-1, N-1, CZERO, CZERO,
1479:      $                            WORK( IR+1 ), LDWRKR )
1480: *
1481: *                    Generate Q in U
1482: *                    (CWorkspace: need N*N+N+M, prefer N*N+N+M*NB)
1483: *                    (RWorkspace: 0)
1484: *
1485:                      CALL CUNGQR( M, M, N, U, LDU, WORK( ITAU ),
1486:      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
1487:                      IE = 1
1488:                      ITAUQ = ITAU
1489:                      ITAUP = ITAUQ + N
1490:                      IWORK = ITAUP + N
1491: *
1492: *                    Bidiagonalize R in WORK(IR)
1493: *                    (CWorkspace: need N*N+3*N, prefer N*N+2*N+2*N*NB)
1494: *                    (RWorkspace: need N)
1495: *
1496:                      CALL CGEBRD( N, N, WORK( IR ), LDWRKR, S,
1497:      $                            RWORK( IE ), WORK( ITAUQ ),
1498:      $                            WORK( ITAUP ), WORK( IWORK ),
1499:      $                            LWORK-IWORK+1, IERR )
1500: *
1501: *                    Generate left bidiagonalizing vectors in WORK(IR)
1502: *                    (CWorkspace: need N*N+3*N, prefer N*N+2*N+N*NB)
1503: *                    (RWorkspace: 0)
1504: *
1505:                      CALL CUNGBR( 'Q', N, N, N, WORK( IR ), LDWRKR,
1506:      $                            WORK( ITAUQ ), WORK( IWORK ),
1507:      $                            LWORK-IWORK+1, IERR )
1508:                      IRWORK = IE + N
1509: *
1510: *                    Perform bidiagonal QR iteration, computing left
1511: *                    singular vectors of R in WORK(IR)
1512: *                    (CWorkspace: need N*N)
1513: *                    (RWorkspace: need BDSPAC)
1514: *
1515:                      CALL CBDSQR( 'U', N, 0, N, 0, S, RWORK( IE ), CDUM,
1516:      $                            1, WORK( IR ), LDWRKR, CDUM, 1,
1517:      $                            RWORK( IRWORK ), INFO )
1518: *
1519: *                    Multiply Q in U by left singular vectors of R in
1520: *                    WORK(IR), storing result in A
1521: *                    (CWorkspace: need N*N)
1522: *                    (RWorkspace: 0)
1523: *
1524:                      CALL CGEMM( 'N', 'N', M, N, N, CONE, U, LDU,
1525:      $                           WORK( IR ), LDWRKR, CZERO, A, LDA )
1526: *
1527: *                    Copy left singular vectors of A from A to U
1528: *
1529:                      CALL CLACPY( 'F', M, N, A, LDA, U, LDU )
1530: *
1531:                   ELSE
1532: *
1533: *                    Insufficient workspace for a fast algorithm
1534: *
1535:                      ITAU = 1
1536:                      IWORK = ITAU + N
1537: *
1538: *                    Compute A=Q*R, copying result to U
1539: *                    (CWorkspace: need 2*N, prefer N+N*NB)
1540: *                    (RWorkspace: 0)
1541: *
1542:                      CALL CGEQRF( M, N, A, LDA, WORK( ITAU ),
1543:      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
1544:                      CALL CLACPY( 'L', M, N, A, LDA, U, LDU )
1545: *
1546: *                    Generate Q in U
1547: *                    (CWorkspace: need N+M, prefer N+M*NB)
1548: *                    (RWorkspace: 0)
1549: *
1550:                      CALL CUNGQR( M, M, N, U, LDU, WORK( ITAU ),
1551:      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
1552:                      IE = 1
1553:                      ITAUQ = ITAU
1554:                      ITAUP = ITAUQ + N
1555:                      IWORK = ITAUP + N
1556: *
1557: *                    Zero out below R in A
1558: *
1559:                      CALL CLASET( 'L', N-1, N-1, CZERO, CZERO,
1560:      $                            A( 2, 1 ), LDA )
1561: *
1562: *                    Bidiagonalize R in A
1563: *                    (CWorkspace: need 3*N, prefer 2*N+2*N*NB)
1564: *                    (RWorkspace: need N)
1565: *
1566:                      CALL CGEBRD( N, N, A, LDA, S, RWORK( IE ),
1567:      $                            WORK( ITAUQ ), WORK( ITAUP ),
1568:      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
1569: *
1570: *                    Multiply Q in U by left bidiagonalizing vectors
1571: *                    in A
1572: *                    (CWorkspace: need 2*N+M, prefer 2*N+M*NB)
1573: *                    (RWorkspace: 0)
1574: *
1575:                      CALL CUNMBR( 'Q', 'R', 'N', M, N, N, A, LDA,
1576:      $                            WORK( ITAUQ ), U, LDU, WORK( IWORK ),
1577:      $                            LWORK-IWORK+1, IERR )
1578:                      IRWORK = IE + N
1579: *
1580: *                    Perform bidiagonal QR iteration, computing left
1581: *                    singular vectors of A in U
1582: *                    (CWorkspace: 0)
1583: *                    (RWorkspace: need BDSPAC)
1584: *
1585:                      CALL CBDSQR( 'U', N, 0, M, 0, S, RWORK( IE ), CDUM,
1586:      $                            1, U, LDU, CDUM, 1, RWORK( IRWORK ),
1587:      $                            INFO )
1588: *
1589:                   END IF
1590: *
1591:                ELSE IF( WNTVO ) THEN
1592: *
1593: *                 Path 8 (M much larger than N, JOBU='A', JOBVT='O')
1594: *                 M left singular vectors to be computed in U and
1595: *                 N right singular vectors to be overwritten on A
1596: *
1597:                   IF( LWORK.GE.2*N*N+MAX( N+M, 3*N ) ) THEN
1598: *
1599: *                    Sufficient workspace for a fast algorithm
1600: *
1601:                      IU = 1
1602:                      IF( LWORK.GE.WRKBL+2*LDA*N ) THEN
1603: *
1604: *                       WORK(IU) is LDA by N and WORK(IR) is LDA by N
1605: *
1606:                         LDWRKU = LDA
1607:                         IR = IU + LDWRKU*N
1608:                         LDWRKR = LDA
1609:                      ELSE IF( LWORK.GE.WRKBL+( LDA+N )*N ) THEN
1610: *
1611: *                       WORK(IU) is LDA by N and WORK(IR) is N by N
1612: *
1613:                         LDWRKU = LDA
1614:                         IR = IU + LDWRKU*N
1615:                         LDWRKR = N
1616:                      ELSE
1617: *
1618: *                       WORK(IU) is N by N and WORK(IR) is N by N
1619: *
1620:                         LDWRKU = N
1621:                         IR = IU + LDWRKU*N
1622:                         LDWRKR = N
1623:                      END IF
1624:                      ITAU = IR + LDWRKR*N
1625:                      IWORK = ITAU + N
1626: *
1627: *                    Compute A=Q*R, copying result to U
1628: *                    (CWorkspace: need 2*N*N+2*N, prefer 2*N*N+N+N*NB)
1629: *                    (RWorkspace: 0)
1630: *
1631:                      CALL CGEQRF( M, N, A, LDA, WORK( ITAU ),
1632:      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
1633:                      CALL CLACPY( 'L', M, N, A, LDA, U, LDU )
1634: *
1635: *                    Generate Q in U
1636: *                    (CWorkspace: need 2*N*N+N+M, prefer 2*N*N+N+M*NB)
1637: *                    (RWorkspace: 0)
1638: *
1639:                      CALL CUNGQR( M, M, N, U, LDU, WORK( ITAU ),
1640:      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
1641: *
1642: *                    Copy R to WORK(IU), zeroing out below it
1643: *
1644:                      CALL CLACPY( 'U', N, N, A, LDA, WORK( IU ),
1645:      $                            LDWRKU )
1646:                      CALL CLASET( 'L', N-1, N-1, CZERO, CZERO,
1647:      $                            WORK( IU+1 ), LDWRKU )
1648:                      IE = 1
1649:                      ITAUQ = ITAU
1650:                      ITAUP = ITAUQ + N
1651:                      IWORK = ITAUP + N
1652: *
1653: *                    Bidiagonalize R in WORK(IU), copying result to
1654: *                    WORK(IR)
1655: *                    (CWorkspace: need   2*N*N+3*N,
1656: *                                 prefer 2*N*N+2*N+2*N*NB)
1657: *                    (RWorkspace: need   N)
1658: *
1659:                      CALL CGEBRD( N, N, WORK( IU ), LDWRKU, S,
1660:      $                            RWORK( IE ), WORK( ITAUQ ),
1661:      $                            WORK( ITAUP ), WORK( IWORK ),
1662:      $                            LWORK-IWORK+1, IERR )
1663:                      CALL CLACPY( 'U', N, N, WORK( IU ), LDWRKU,
1664:      $                            WORK( IR ), LDWRKR )
1665: *
1666: *                    Generate left bidiagonalizing vectors in WORK(IU)
1667: *                    (CWorkspace: need 2*N*N+3*N, prefer 2*N*N+2*N+N*NB)
1668: *                    (RWorkspace: 0)
1669: *
1670:                      CALL CUNGBR( 'Q', N, N, N, WORK( IU ), LDWRKU,
1671:      $                            WORK( ITAUQ ), WORK( IWORK ),
1672:      $                            LWORK-IWORK+1, IERR )
1673: *
1674: *                    Generate right bidiagonalizing vectors in WORK(IR)
1675: *                    (CWorkspace: need   2*N*N+3*N-1,
1676: *                                 prefer 2*N*N+2*N+(N-1)*NB)
1677: *                    (RWorkspace: 0)
1678: *
1679:                      CALL CUNGBR( 'P', N, N, N, WORK( IR ), LDWRKR,
1680:      $                            WORK( ITAUP ), WORK( IWORK ),
1681:      $                            LWORK-IWORK+1, IERR )
1682:                      IRWORK = IE + N
1683: *
1684: *                    Perform bidiagonal QR iteration, computing left
1685: *                    singular vectors of R in WORK(IU) and computing
1686: *                    right singular vectors of R in WORK(IR)
1687: *                    (CWorkspace: need 2*N*N)
1688: *                    (RWorkspace: need BDSPAC)
1689: *
1690:                      CALL CBDSQR( 'U', N, N, N, 0, S, RWORK( IE ),
1691:      $                            WORK( IR ), LDWRKR, WORK( IU ),
1692:      $                            LDWRKU, CDUM, 1, RWORK( IRWORK ),
1693:      $                            INFO )
1694: *
1695: *                    Multiply Q in U by left singular vectors of R in
1696: *                    WORK(IU), storing result in A
1697: *                    (CWorkspace: need N*N)
1698: *                    (RWorkspace: 0)
1699: *
1700:                      CALL CGEMM( 'N', 'N', M, N, N, CONE, U, LDU,
1701:      $                           WORK( IU ), LDWRKU, CZERO, A, LDA )
1702: *
1703: *                    Copy left singular vectors of A from A to U
1704: *
1705:                      CALL CLACPY( 'F', M, N, A, LDA, U, LDU )
1706: *
1707: *                    Copy right singular vectors of R from WORK(IR) to A
1708: *
1709:                      CALL CLACPY( 'F', N, N, WORK( IR ), LDWRKR, A,
1710:      $                            LDA )
1711: *
1712:                   ELSE
1713: *
1714: *                    Insufficient workspace for a fast algorithm
1715: *
1716:                      ITAU = 1
1717:                      IWORK = ITAU + N
1718: *
1719: *                    Compute A=Q*R, copying result to U
1720: *                    (CWorkspace: need 2*N, prefer N+N*NB)
1721: *                    (RWorkspace: 0)
1722: *
1723:                      CALL CGEQRF( M, N, A, LDA, WORK( ITAU ),
1724:      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
1725:                      CALL CLACPY( 'L', M, N, A, LDA, U, LDU )
1726: *
1727: *                    Generate Q in U
1728: *                    (CWorkspace: need N+M, prefer N+M*NB)
1729: *                    (RWorkspace: 0)
1730: *
1731:                      CALL CUNGQR( M, M, N, U, LDU, WORK( ITAU ),
1732:      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
1733:                      IE = 1
1734:                      ITAUQ = ITAU
1735:                      ITAUP = ITAUQ + N
1736:                      IWORK = ITAUP + N
1737: *
1738: *                    Zero out below R in A
1739: *
1740:                      CALL CLASET( 'L', N-1, N-1, CZERO, CZERO,
1741:      $                            A( 2, 1 ), LDA )
1742: *
1743: *                    Bidiagonalize R in A
1744: *                    (CWorkspace: need 3*N, prefer 2*N+2*N*NB)
1745: *                    (RWorkspace: need N)
1746: *
1747:                      CALL CGEBRD( N, N, A, LDA, S, RWORK( IE ),
1748:      $                            WORK( ITAUQ ), WORK( ITAUP ),
1749:      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
1750: *
1751: *                    Multiply Q in U by left bidiagonalizing vectors
1752: *                    in A
1753: *                    (CWorkspace: need 2*N+M, prefer 2*N+M*NB)
1754: *                    (RWorkspace: 0)
1755: *
1756:                      CALL CUNMBR( 'Q', 'R', 'N', M, N, N, A, LDA,
1757:      $                            WORK( ITAUQ ), U, LDU, WORK( IWORK ),
1758:      $                            LWORK-IWORK+1, IERR )
1759: *
1760: *                    Generate right bidiagonalizing vectors in A
1761: *                    (CWorkspace: need 3*N-1, prefer 2*N+(N-1)*NB)
1762: *                    (RWorkspace: 0)
1763: *
1764:                      CALL CUNGBR( 'P', N, N, N, A, LDA, WORK( ITAUP ),
1765:      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
1766:                      IRWORK = IE + N
1767: *
1768: *                    Perform bidiagonal QR iteration, computing left
1769: *                    singular vectors of A in U and computing right
1770: *                    singular vectors of A in A
1771: *                    (CWorkspace: 0)
1772: *                    (RWorkspace: need BDSPAC)
1773: *
1774:                      CALL CBDSQR( 'U', N, N, M, 0, S, RWORK( IE ), A,
1775:      $                            LDA, U, LDU, CDUM, 1, RWORK( IRWORK ),
1776:      $                            INFO )
1777: *
1778:                   END IF
1779: *
1780:                ELSE IF( WNTVAS ) THEN
1781: *
1782: *                 Path 9 (M much larger than N, JOBU='A', JOBVT='S'
1783: *                         or 'A')
1784: *                 M left singular vectors to be computed in U and
1785: *                 N right singular vectors to be computed in VT
1786: *
1787:                   IF( LWORK.GE.N*N+MAX( N+M, 3*N ) ) THEN
1788: *
1789: *                    Sufficient workspace for a fast algorithm
1790: *
1791:                      IU = 1
1792:                      IF( LWORK.GE.WRKBL+LDA*N ) THEN
1793: *
1794: *                       WORK(IU) is LDA by N
1795: *
1796:                         LDWRKU = LDA
1797:                      ELSE
1798: *
1799: *                       WORK(IU) is N by N
1800: *
1801:                         LDWRKU = N
1802:                      END IF
1803:                      ITAU = IU + LDWRKU*N
1804:                      IWORK = ITAU + N
1805: *
1806: *                    Compute A=Q*R, copying result to U
1807: *                    (CWorkspace: need N*N+2*N, prefer N*N+N+N*NB)
1808: *                    (RWorkspace: 0)
1809: *
1810:                      CALL CGEQRF( M, N, A, LDA, WORK( ITAU ),
1811:      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
1812:                      CALL CLACPY( 'L', M, N, A, LDA, U, LDU )
1813: *
1814: *                    Generate Q in U
1815: *                    (CWorkspace: need N*N+N+M, prefer N*N+N+M*NB)
1816: *                    (RWorkspace: 0)
1817: *
1818:                      CALL CUNGQR( M, M, N, U, LDU, WORK( ITAU ),
1819:      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
1820: *
1821: *                    Copy R to WORK(IU), zeroing out below it
1822: *
1823:                      CALL CLACPY( 'U', N, N, A, LDA, WORK( IU ),
1824:      $                            LDWRKU )
1825:                      CALL CLASET( 'L', N-1, N-1, CZERO, CZERO,
1826:      $                            WORK( IU+1 ), LDWRKU )
1827:                      IE = 1
1828:                      ITAUQ = ITAU
1829:                      ITAUP = ITAUQ + N
1830:                      IWORK = ITAUP + N
1831: *
1832: *                    Bidiagonalize R in WORK(IU), copying result to VT
1833: *                    (CWorkspace: need N*N+3*N, prefer N*N+2*N+2*N*NB)
1834: *                    (RWorkspace: need N)
1835: *
1836:                      CALL CGEBRD( N, N, WORK( IU ), LDWRKU, S,
1837:      $                            RWORK( IE ), WORK( ITAUQ ),
1838:      $                            WORK( ITAUP ), WORK( IWORK ),
1839:      $                            LWORK-IWORK+1, IERR )
1840:                      CALL CLACPY( 'U', N, N, WORK( IU ), LDWRKU, VT,
1841:      $                            LDVT )
1842: *
1843: *                    Generate left bidiagonalizing vectors in WORK(IU)
1844: *                    (CWorkspace: need N*N+3*N, prefer N*N+2*N+N*NB)
1845: *                    (RWorkspace: 0)
1846: *
1847:                      CALL CUNGBR( 'Q', N, N, N, WORK( IU ), LDWRKU,
1848:      $                            WORK( ITAUQ ), WORK( IWORK ),
1849:      $                            LWORK-IWORK+1, IERR )
1850: *
1851: *                    Generate right bidiagonalizing vectors in VT
1852: *                    (CWorkspace: need   N*N+3*N-1,
1853: *                                 prefer N*N+2*N+(N-1)*NB)
1854: *                    (RWorkspace: need   0)
1855: *
1856:                      CALL CUNGBR( 'P', N, N, N, VT, LDVT, WORK( ITAUP ),
1857:      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
1858:                      IRWORK = IE + N
1859: *
1860: *                    Perform bidiagonal QR iteration, computing left
1861: *                    singular vectors of R in WORK(IU) and computing
1862: *                    right singular vectors of R in VT
1863: *                    (CWorkspace: need N*N)
1864: *                    (RWorkspace: need BDSPAC)
1865: *
1866:                      CALL CBDSQR( 'U', N, N, N, 0, S, RWORK( IE ), VT,
1867:      $                            LDVT, WORK( IU ), LDWRKU, CDUM, 1,
1868:      $                            RWORK( IRWORK ), INFO )
1869: *
1870: *                    Multiply Q in U by left singular vectors of R in
1871: *                    WORK(IU), storing result in A
1872: *                    (CWorkspace: need N*N)
1873: *                    (RWorkspace: 0)
1874: *
1875:                      CALL CGEMM( 'N', 'N', M, N, N, CONE, U, LDU,
1876:      $                           WORK( IU ), LDWRKU, CZERO, A, LDA )
1877: *
1878: *                    Copy left singular vectors of A from A to U
1879: *
1880:                      CALL CLACPY( 'F', M, N, A, LDA, U, LDU )
1881: *
1882:                   ELSE
1883: *
1884: *                    Insufficient workspace for a fast algorithm
1885: *
1886:                      ITAU = 1
1887:                      IWORK = ITAU + N
1888: *
1889: *                    Compute A=Q*R, copying result to U
1890: *                    (CWorkspace: need 2*N, prefer N+N*NB)
1891: *                    (RWorkspace: 0)
1892: *
1893:                      CALL CGEQRF( M, N, A, LDA, WORK( ITAU ),
1894:      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
1895:                      CALL CLACPY( 'L', M, N, A, LDA, U, LDU )
1896: *
1897: *                    Generate Q in U
1898: *                    (CWorkspace: need N+M, prefer N+M*NB)
1899: *                    (RWorkspace: 0)
1900: *
1901:                      CALL CUNGQR( M, M, N, U, LDU, WORK( ITAU ),
1902:      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
1903: *
1904: *                    Copy R from A to VT, zeroing out below it
1905: *
1906:                      CALL CLACPY( 'U', N, N, A, LDA, VT, LDVT )
1907:                      IF( N.GT.1 )
1908:      $                  CALL CLASET( 'L', N-1, N-1, CZERO, CZERO,
1909:      $                               VT( 2, 1 ), LDVT )
1910:                      IE = 1
1911:                      ITAUQ = ITAU
1912:                      ITAUP = ITAUQ + N
1913:                      IWORK = ITAUP + N
1914: *
1915: *                    Bidiagonalize R in VT
1916: *                    (CWorkspace: need 3*N, prefer 2*N+2*N*NB)
1917: *                    (RWorkspace: need N)
1918: *
1919:                      CALL CGEBRD( N, N, VT, LDVT, S, RWORK( IE ),
1920:      $                            WORK( ITAUQ ), WORK( ITAUP ),
1921:      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
1922: *
1923: *                    Multiply Q in U by left bidiagonalizing vectors
1924: *                    in VT
1925: *                    (CWorkspace: need 2*N+M, prefer 2*N+M*NB)
1926: *                    (RWorkspace: 0)
1927: *
1928:                      CALL CUNMBR( 'Q', 'R', 'N', M, N, N, VT, LDVT,
1929:      $                            WORK( ITAUQ ), U, LDU, WORK( IWORK ),
1930:      $                            LWORK-IWORK+1, IERR )
1931: *
1932: *                    Generate right bidiagonalizing vectors in VT
1933: *                    (CWorkspace: need 3*N-1, prefer 2*N+(N-1)*NB)
1934: *                    (RWorkspace: 0)
1935: *
1936:                      CALL CUNGBR( 'P', N, N, N, VT, LDVT, WORK( ITAUP ),
1937:      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
1938:                      IRWORK = IE + N
1939: *
1940: *                    Perform bidiagonal QR iteration, computing left
1941: *                    singular vectors of A in U and computing right
1942: *                    singular vectors of A in VT
1943: *                    (CWorkspace: 0)
1944: *                    (RWorkspace: need BDSPAC)
1945: *
1946:                      CALL CBDSQR( 'U', N, N, M, 0, S, RWORK( IE ), VT,
1947:      $                            LDVT, U, LDU, CDUM, 1,
1948:      $                            RWORK( IRWORK ), INFO )
1949: *
1950:                   END IF
1951: *
1952:                END IF
1953: *
1954:             END IF
1955: *
1956:          ELSE
1957: *
1958: *           M .LT. MNTHR
1959: *
1960: *           Path 10 (M at least N, but not much larger)
1961: *           Reduce to bidiagonal form without QR decomposition
1962: *
1963:             IE = 1
1964:             ITAUQ = 1
1965:             ITAUP = ITAUQ + N
1966:             IWORK = ITAUP + N
1967: *
1968: *           Bidiagonalize A
1969: *           (CWorkspace: need 2*N+M, prefer 2*N+(M+N)*NB)
1970: *           (RWorkspace: need N)
1971: *
1972:             CALL CGEBRD( M, N, A, LDA, S, RWORK( IE ), WORK( ITAUQ ),
1973:      $                   WORK( ITAUP ), WORK( IWORK ), LWORK-IWORK+1,
1974:      $                   IERR )
1975:             IF( WNTUAS ) THEN
1976: *
1977: *              If left singular vectors desired in U, copy result to U
1978: *              and generate left bidiagonalizing vectors in U
1979: *              (CWorkspace: need 2*N+NCU, prefer 2*N+NCU*NB)
1980: *              (RWorkspace: 0)
1981: *
1982:                CALL CLACPY( 'L', M, N, A, LDA, U, LDU )
1983:                IF( WNTUS )
1984:      $            NCU = N
1985:                IF( WNTUA )
1986:      $            NCU = M
1987:                CALL CUNGBR( 'Q', M, NCU, N, U, LDU, WORK( ITAUQ ),
1988:      $                      WORK( IWORK ), LWORK-IWORK+1, IERR )
1989:             END IF
1990:             IF( WNTVAS ) THEN
1991: *
1992: *              If right singular vectors desired in VT, copy result to
1993: *              VT and generate right bidiagonalizing vectors in VT
1994: *              (CWorkspace: need 3*N-1, prefer 2*N+(N-1)*NB)
1995: *              (RWorkspace: 0)
1996: *
1997:                CALL CLACPY( 'U', N, N, A, LDA, VT, LDVT )
1998:                CALL CUNGBR( 'P', N, N, N, VT, LDVT, WORK( ITAUP ),
1999:      $                      WORK( IWORK ), LWORK-IWORK+1, IERR )
2000:             END IF
2001:             IF( WNTUO ) THEN
2002: *
2003: *              If left singular vectors desired in A, generate left
2004: *              bidiagonalizing vectors in A
2005: *              (CWorkspace: need 3*N, prefer 2*N+N*NB)
2006: *              (RWorkspace: 0)
2007: *
2008:                CALL CUNGBR( 'Q', M, N, N, A, LDA, WORK( ITAUQ ),
2009:      $                      WORK( IWORK ), LWORK-IWORK+1, IERR )
2010:             END IF
2011:             IF( WNTVO ) THEN
2012: *
2013: *              If right singular vectors desired in A, generate right
2014: *              bidiagonalizing vectors in A
2015: *              (CWorkspace: need 3*N-1, prefer 2*N+(N-1)*NB)
2016: *              (RWorkspace: 0)
2017: *
2018:                CALL CUNGBR( 'P', N, N, N, A, LDA, WORK( ITAUP ),
2019:      $                      WORK( IWORK ), LWORK-IWORK+1, IERR )
2020:             END IF
2021:             IRWORK = IE + N
2022:             IF( WNTUAS .OR. WNTUO )
2023:      $         NRU = M
2024:             IF( WNTUN )
2025:      $         NRU = 0
2026:             IF( WNTVAS .OR. WNTVO )
2027:      $         NCVT = N
2028:             IF( WNTVN )
2029:      $         NCVT = 0
2030:             IF( ( .NOT.WNTUO ) .AND. ( .NOT.WNTVO ) ) THEN
2031: *
2032: *              Perform bidiagonal QR iteration, if desired, computing
2033: *              left singular vectors in U and computing right singular
2034: *              vectors in VT
2035: *              (CWorkspace: 0)
2036: *              (RWorkspace: need BDSPAC)
2037: *
2038:                CALL CBDSQR( 'U', N, NCVT, NRU, 0, S, RWORK( IE ), VT,
2039:      $                      LDVT, U, LDU, CDUM, 1, RWORK( IRWORK ),
2040:      $                      INFO )
2041:             ELSE IF( ( .NOT.WNTUO ) .AND. WNTVO ) THEN
2042: *
2043: *              Perform bidiagonal QR iteration, if desired, computing
2044: *              left singular vectors in U and computing right singular
2045: *              vectors in A
2046: *              (CWorkspace: 0)
2047: *              (RWorkspace: need BDSPAC)
2048: *
2049:                CALL CBDSQR( 'U', N, NCVT, NRU, 0, S, RWORK( IE ), A,
2050:      $                      LDA, U, LDU, CDUM, 1, RWORK( IRWORK ),
2051:      $                      INFO )
2052:             ELSE
2053: *
2054: *              Perform bidiagonal QR iteration, if desired, computing
2055: *              left singular vectors in A and computing right singular
2056: *              vectors in VT
2057: *              (CWorkspace: 0)
2058: *              (RWorkspace: need BDSPAC)
2059: *
2060:                CALL CBDSQR( 'U', N, NCVT, NRU, 0, S, RWORK( IE ), VT,
2061:      $                      LDVT, A, LDA, CDUM, 1, RWORK( IRWORK ),
2062:      $                      INFO )
2063:             END IF
2064: *
2065:          END IF
2066: *
2067:       ELSE
2068: *
2069: *        A has more columns than rows. If A has sufficiently more
2070: *        columns than rows, first reduce using the LQ decomposition (if
2071: *        sufficient workspace available)
2072: *
2073:          IF( N.GE.MNTHR ) THEN
2074: *
2075:             IF( WNTVN ) THEN
2076: *
2077: *              Path 1t(N much larger than M, JOBVT='N')
2078: *              No right singular vectors to be computed
2079: *
2080:                ITAU = 1
2081:                IWORK = ITAU + M
2082: *
2083: *              Compute A=L*Q
2084: *              (CWorkspace: need 2*M, prefer M+M*NB)
2085: *              (RWorkspace: 0)
2086: *
2087:                CALL CGELQF( M, N, A, LDA, WORK( ITAU ), WORK( IWORK ),
2088:      $                      LWORK-IWORK+1, IERR )
2089: *
2090: *              Zero out above L
2091: *
2092:                CALL CLASET( 'U', M-1, M-1, CZERO, CZERO, A( 1, 2 ),
2093:      $                      LDA )
2094:                IE = 1
2095:                ITAUQ = 1
2096:                ITAUP = ITAUQ + M
2097:                IWORK = ITAUP + M
2098: *
2099: *              Bidiagonalize L in A
2100: *              (CWorkspace: need 3*M, prefer 2*M+2*M*NB)
2101: *              (RWorkspace: need M)
2102: *
2103:                CALL CGEBRD( M, M, A, LDA, S, RWORK( IE ), WORK( ITAUQ ),
2104:      $                      WORK( ITAUP ), WORK( IWORK ), LWORK-IWORK+1,
2105:      $                      IERR )
2106:                IF( WNTUO .OR. WNTUAS ) THEN
2107: *
2108: *                 If left singular vectors desired, generate Q
2109: *                 (CWorkspace: need 3*M, prefer 2*M+M*NB)
2110: *                 (RWorkspace: 0)
2111: *
2112:                   CALL CUNGBR( 'Q', M, M, M, A, LDA, WORK( ITAUQ ),
2113:      $                         WORK( IWORK ), LWORK-IWORK+1, IERR )
2114:                END IF
2115:                IRWORK = IE + M
2116:                NRU = 0
2117:                IF( WNTUO .OR. WNTUAS )
2118:      $            NRU = M
2119: *
2120: *              Perform bidiagonal QR iteration, computing left singular
2121: *              vectors of A in A if desired
2122: *              (CWorkspace: 0)
2123: *              (RWorkspace: need BDSPAC)
2124: *
2125:                CALL CBDSQR( 'U', M, 0, NRU, 0, S, RWORK( IE ), CDUM, 1,
2126:      $                      A, LDA, CDUM, 1, RWORK( IRWORK ), INFO )
2127: *
2128: *              If left singular vectors desired in U, copy them there
2129: *
2130:                IF( WNTUAS )
2131:      $            CALL CLACPY( 'F', M, M, A, LDA, U, LDU )
2132: *
2133:             ELSE IF( WNTVO .AND. WNTUN ) THEN
2134: *
2135: *              Path 2t(N much larger than M, JOBU='N', JOBVT='O')
2136: *              M right singular vectors to be overwritten on A and
2137: *              no left singular vectors to be computed
2138: *
2139:                IF( LWORK.GE.M*M+3*M ) THEN
2140: *
2141: *                 Sufficient workspace for a fast algorithm
2142: *
2143:                   IR = 1
2144:                   IF( LWORK.GE.MAX( WRKBL, LDA*N )+LDA*M ) THEN
2145: *
2146: *                    WORK(IU) is LDA by N and WORK(IR) is LDA by M
2147: *
2148:                      LDWRKU = LDA
2149:                      CHUNK = N
2150:                      LDWRKR = LDA
2151:                   ELSE IF( LWORK.GE.MAX( WRKBL, LDA*N )+M*M ) THEN
2152: *
2153: *                    WORK(IU) is LDA by N and WORK(IR) is M by M
2154: *
2155:                      LDWRKU = LDA
2156:                      CHUNK = N
2157:                      LDWRKR = M
2158:                   ELSE
2159: *
2160: *                    WORK(IU) is M by CHUNK and WORK(IR) is M by M
2161: *
2162:                      LDWRKU = M
2163:                      CHUNK = ( LWORK-M*M ) / M
2164:                      LDWRKR = M
2165:                   END IF
2166:                   ITAU = IR + LDWRKR*M
2167:                   IWORK = ITAU + M
2168: *
2169: *                 Compute A=L*Q
2170: *                 (CWorkspace: need M*M+2*M, prefer M*M+M+M*NB)
2171: *                 (RWorkspace: 0)
2172: *
2173:                   CALL CGELQF( M, N, A, LDA, WORK( ITAU ),
2174:      $                         WORK( IWORK ), LWORK-IWORK+1, IERR )
2175: *
2176: *                 Copy L to WORK(IR) and zero out above it
2177: *
2178:                   CALL CLACPY( 'L', M, M, A, LDA, WORK( IR ), LDWRKR )
2179:                   CALL CLASET( 'U', M-1, M-1, CZERO, CZERO,
2180:      $                         WORK( IR+LDWRKR ), LDWRKR )
2181: *
2182: *                 Generate Q in A
2183: *                 (CWorkspace: need M*M+2*M, prefer M*M+M+M*NB)
2184: *                 (RWorkspace: 0)
2185: *
2186:                   CALL CUNGLQ( M, N, M, A, LDA, WORK( ITAU ),
2187:      $                         WORK( IWORK ), LWORK-IWORK+1, IERR )
2188:                   IE = 1
2189:                   ITAUQ = ITAU
2190:                   ITAUP = ITAUQ + M
2191:                   IWORK = ITAUP + M
2192: *
2193: *                 Bidiagonalize L in WORK(IR)
2194: *                 (CWorkspace: need M*M+3*M, prefer M*M+2*M+2*M*NB)
2195: *                 (RWorkspace: need M)
2196: *
2197:                   CALL CGEBRD( M, M, WORK( IR ), LDWRKR, S, RWORK( IE ),
2198:      $                         WORK( ITAUQ ), WORK( ITAUP ),
2199:      $                         WORK( IWORK ), LWORK-IWORK+1, IERR )
2200: *
2201: *                 Generate right vectors bidiagonalizing L
2202: *                 (CWorkspace: need M*M+3*M-1, prefer M*M+2*M+(M-1)*NB)
2203: *                 (RWorkspace: 0)
2204: *
2205:                   CALL CUNGBR( 'P', M, M, M, WORK( IR ), LDWRKR,
2206:      $                         WORK( ITAUP ), WORK( IWORK ),
2207:      $                         LWORK-IWORK+1, IERR )
2208:                   IRWORK = IE + M
2209: *
2210: *                 Perform bidiagonal QR iteration, computing right
2211: *                 singular vectors of L in WORK(IR)
2212: *                 (CWorkspace: need M*M)
2213: *                 (RWorkspace: need BDSPAC)
2214: *
2215:                   CALL CBDSQR( 'U', M, M, 0, 0, S, RWORK( IE ),
2216:      $                         WORK( IR ), LDWRKR, CDUM, 1, CDUM, 1,
2217:      $                         RWORK( IRWORK ), INFO )
2218:                   IU = ITAUQ
2219: *
2220: *                 Multiply right singular vectors of L in WORK(IR) by Q
2221: *                 in A, storing result in WORK(IU) and copying to A
2222: *                 (CWorkspace: need M*M+M, prefer M*M+M*N)
2223: *                 (RWorkspace: 0)
2224: *
2225:                   DO 30 I = 1, N, CHUNK
2226:                      BLK = MIN( N-I+1, CHUNK )
2227:                      CALL CGEMM( 'N', 'N', M, BLK, M, CONE, WORK( IR ),
2228:      $                           LDWRKR, A( 1, I ), LDA, CZERO,
2229:      $                           WORK( IU ), LDWRKU )
2230:                      CALL CLACPY( 'F', M, BLK, WORK( IU ), LDWRKU,
2231:      $                            A( 1, I ), LDA )
2232:    30             CONTINUE
2233: *
2234:                ELSE
2235: *
2236: *                 Insufficient workspace for a fast algorithm
2237: *
2238:                   IE = 1
2239:                   ITAUQ = 1
2240:                   ITAUP = ITAUQ + M
2241:                   IWORK = ITAUP + M
2242: *
2243: *                 Bidiagonalize A
2244: *                 (CWorkspace: need 2*M+N, prefer 2*M+(M+N)*NB)
2245: *                 (RWorkspace: need M)
2246: *
2247:                   CALL CGEBRD( M, N, A, LDA, S, RWORK( IE ),
2248:      $                         WORK( ITAUQ ), WORK( ITAUP ),
2249:      $                         WORK( IWORK ), LWORK-IWORK+1, IERR )
2250: *
2251: *                 Generate right vectors bidiagonalizing A
2252: *                 (CWorkspace: need 3*M, prefer 2*M+M*NB)
2253: *                 (RWorkspace: 0)
2254: *
2255:                   CALL CUNGBR( 'P', M, N, M, A, LDA, WORK( ITAUP ),
2256:      $                         WORK( IWORK ), LWORK-IWORK+1, IERR )
2257:                   IRWORK = IE + M
2258: *
2259: *                 Perform bidiagonal QR iteration, computing right
2260: *                 singular vectors of A in A
2261: *                 (CWorkspace: 0)
2262: *                 (RWorkspace: need BDSPAC)
2263: *
2264:                   CALL CBDSQR( 'L', M, N, 0, 0, S, RWORK( IE ), A, LDA,
2265:      $                         CDUM, 1, CDUM, 1, RWORK( IRWORK ), INFO )
2266: *
2267:                END IF
2268: *
2269:             ELSE IF( WNTVO .AND. WNTUAS ) THEN
2270: *
2271: *              Path 3t(N much larger than M, JOBU='S' or 'A', JOBVT='O')
2272: *              M right singular vectors to be overwritten on A and
2273: *              M left singular vectors to be computed in U
2274: *
2275:                IF( LWORK.GE.M*M+3*M ) THEN
2276: *
2277: *                 Sufficient workspace for a fast algorithm
2278: *
2279:                   IR = 1
2280:                   IF( LWORK.GE.MAX( WRKBL, LDA*N )+LDA*M ) THEN
2281: *
2282: *                    WORK(IU) is LDA by N and WORK(IR) is LDA by M
2283: *
2284:                      LDWRKU = LDA
2285:                      CHUNK = N
2286:                      LDWRKR = LDA
2287:                   ELSE IF( LWORK.GE.MAX( WRKBL, LDA*N )+M*M ) THEN
2288: *
2289: *                    WORK(IU) is LDA by N and WORK(IR) is M by M
2290: *
2291:                      LDWRKU = LDA
2292:                      CHUNK = N
2293:                      LDWRKR = M
2294:                   ELSE
2295: *
2296: *                    WORK(IU) is M by CHUNK and WORK(IR) is M by M
2297: *
2298:                      LDWRKU = M
2299:                      CHUNK = ( LWORK-M*M ) / M
2300:                      LDWRKR = M
2301:                   END IF
2302:                   ITAU = IR + LDWRKR*M
2303:                   IWORK = ITAU + M
2304: *
2305: *                 Compute A=L*Q
2306: *                 (CWorkspace: need M*M+2*M, prefer M*M+M+M*NB)
2307: *                 (RWorkspace: 0)
2308: *
2309:                   CALL CGELQF( M, N, A, LDA, WORK( ITAU ),
2310:      $                         WORK( IWORK ), LWORK-IWORK+1, IERR )
2311: *
2312: *                 Copy L to U, zeroing about above it
2313: *
2314:                   CALL CLACPY( 'L', M, M, A, LDA, U, LDU )
2315:                   CALL CLASET( 'U', M-1, M-1, CZERO, CZERO, U( 1, 2 ),
2316:      $                         LDU )
2317: *
2318: *                 Generate Q in A
2319: *                 (CWorkspace: need M*M+2*M, prefer M*M+M+M*NB)
2320: *                 (RWorkspace: 0)
2321: *
2322:                   CALL CUNGLQ( M, N, M, A, LDA, WORK( ITAU ),
2323:      $                         WORK( IWORK ), LWORK-IWORK+1, IERR )
2324:                   IE = 1
2325:                   ITAUQ = ITAU
2326:                   ITAUP = ITAUQ + M
2327:                   IWORK = ITAUP + M
2328: *
2329: *                 Bidiagonalize L in U, copying result to WORK(IR)
2330: *                 (CWorkspace: need M*M+3*M, prefer M*M+2*M+2*M*NB)
2331: *                 (RWorkspace: need M)
2332: *
2333:                   CALL CGEBRD( M, M, U, LDU, S, RWORK( IE ),
2334:      $                         WORK( ITAUQ ), WORK( ITAUP ),
2335:      $                         WORK( IWORK ), LWORK-IWORK+1, IERR )
2336:                   CALL CLACPY( 'U', M, M, U, LDU, WORK( IR ), LDWRKR )
2337: *
2338: *                 Generate right vectors bidiagonalizing L in WORK(IR)
2339: *                 (CWorkspace: need M*M+3*M-1, prefer M*M+2*M+(M-1)*NB)
2340: *                 (RWorkspace: 0)
2341: *
2342:                   CALL CUNGBR( 'P', M, M, M, WORK( IR ), LDWRKR,
2343:      $                         WORK( ITAUP ), WORK( IWORK ),
2344:      $                         LWORK-IWORK+1, IERR )
2345: *
2346: *                 Generate left vectors bidiagonalizing L in U
2347: *                 (CWorkspace: need M*M+3*M, prefer M*M+2*M+M*NB)
2348: *                 (RWorkspace: 0)
2349: *
2350:                   CALL CUNGBR( 'Q', M, M, M, U, LDU, WORK( ITAUQ ),
2351:      $                         WORK( IWORK ), LWORK-IWORK+1, IERR )
2352:                   IRWORK = IE + M
2353: *
2354: *                 Perform bidiagonal QR iteration, computing left
2355: *                 singular vectors of L in U, and computing right
2356: *                 singular vectors of L in WORK(IR)
2357: *                 (CWorkspace: need M*M)
2358: *                 (RWorkspace: need BDSPAC)
2359: *
2360:                   CALL CBDSQR( 'U', M, M, M, 0, S, RWORK( IE ),
2361:      $                         WORK( IR ), LDWRKR, U, LDU, CDUM, 1,
2362:      $                         RWORK( IRWORK ), INFO )
2363:                   IU = ITAUQ
2364: *
2365: *                 Multiply right singular vectors of L in WORK(IR) by Q
2366: *                 in A, storing result in WORK(IU) and copying to A
2367: *                 (CWorkspace: need M*M+M, prefer M*M+M*N))
2368: *                 (RWorkspace: 0)
2369: *
2370:                   DO 40 I = 1, N, CHUNK
2371:                      BLK = MIN( N-I+1, CHUNK )
2372:                      CALL CGEMM( 'N', 'N', M, BLK, M, CONE, WORK( IR ),
2373:      $                           LDWRKR, A( 1, I ), LDA, CZERO,
2374:      $                           WORK( IU ), LDWRKU )
2375:                      CALL CLACPY( 'F', M, BLK, WORK( IU ), LDWRKU,
2376:      $                            A( 1, I ), LDA )
2377:    40             CONTINUE
2378: *
2379:                ELSE
2380: *
2381: *                 Insufficient workspace for a fast algorithm
2382: *
2383:                   ITAU = 1
2384:                   IWORK = ITAU + M
2385: *
2386: *                 Compute A=L*Q
2387: *                 (CWorkspace: need 2*M, prefer M+M*NB)
2388: *                 (RWorkspace: 0)
2389: *
2390:                   CALL CGELQF( M, N, A, LDA, WORK( ITAU ),
2391:      $                         WORK( IWORK ), LWORK-IWORK+1, IERR )
2392: *
2393: *                 Copy L to U, zeroing out above it
2394: *
2395:                   CALL CLACPY( 'L', M, M, A, LDA, U, LDU )
2396:                   CALL CLASET( 'U', M-1, M-1, CZERO, CZERO, U( 1, 2 ),
2397:      $                         LDU )
2398: *
2399: *                 Generate Q in A
2400: *                 (CWorkspace: need 2*M, prefer M+M*NB)
2401: *                 (RWorkspace: 0)
2402: *
2403:                   CALL CUNGLQ( M, N, M, A, LDA, WORK( ITAU ),
2404:      $                         WORK( IWORK ), LWORK-IWORK+1, IERR )
2405:                   IE = 1
2406:                   ITAUQ = ITAU
2407:                   ITAUP = ITAUQ + M
2408:                   IWORK = ITAUP + M
2409: *
2410: *                 Bidiagonalize L in U
2411: *                 (CWorkspace: need 3*M, prefer 2*M+2*M*NB)
2412: *                 (RWorkspace: need M)
2413: *
2414:                   CALL CGEBRD( M, M, U, LDU, S, RWORK( IE ),
2415:      $                         WORK( ITAUQ ), WORK( ITAUP ),
2416:      $                         WORK( IWORK ), LWORK-IWORK+1, IERR )
2417: *
2418: *                 Multiply right vectors bidiagonalizing L by Q in A
2419: *                 (CWorkspace: need 2*M+N, prefer 2*M+N*NB)
2420: *                 (RWorkspace: 0)
2421: *
2422:                   CALL CUNMBR( 'P', 'L', 'C', M, N, M, U, LDU,
2423:      $                         WORK( ITAUP ), A, LDA, WORK( IWORK ),
2424:      $                         LWORK-IWORK+1, IERR )
2425: *
2426: *                 Generate left vectors bidiagonalizing L in U
2427: *                 (CWorkspace: need 3*M, prefer 2*M+M*NB)
2428: *                 (RWorkspace: 0)
2429: *
2430:                   CALL CUNGBR( 'Q', M, M, M, U, LDU, WORK( ITAUQ ),
2431:      $                         WORK( IWORK ), LWORK-IWORK+1, IERR )
2432:                   IRWORK = IE + M
2433: *
2434: *                 Perform bidiagonal QR iteration, computing left
2435: *                 singular vectors of A in U and computing right
2436: *                 singular vectors of A in A
2437: *                 (CWorkspace: 0)
2438: *                 (RWorkspace: need BDSPAC)
2439: *
2440:                   CALL CBDSQR( 'U', M, N, M, 0, S, RWORK( IE ), A, LDA,
2441:      $                         U, LDU, CDUM, 1, RWORK( IRWORK ), INFO )
2442: *
2443:                END IF
2444: *
2445:             ELSE IF( WNTVS ) THEN
2446: *
2447:                IF( WNTUN ) THEN
2448: *
2449: *                 Path 4t(N much larger than M, JOBU='N', JOBVT='S')
2450: *                 M right singular vectors to be computed in VT and
2451: *                 no left singular vectors to be computed
2452: *
2453:                   IF( LWORK.GE.M*M+3*M ) THEN
2454: *
2455: *                    Sufficient workspace for a fast algorithm
2456: *
2457:                      IR = 1
2458:                      IF( LWORK.GE.WRKBL+LDA*M ) THEN
2459: *
2460: *                       WORK(IR) is LDA by M
2461: *
2462:                         LDWRKR = LDA
2463:                      ELSE
2464: *
2465: *                       WORK(IR) is M by M
2466: *
2467:                         LDWRKR = M
2468:                      END IF
2469:                      ITAU = IR + LDWRKR*M
2470:                      IWORK = ITAU + M
2471: *
2472: *                    Compute A=L*Q
2473: *                    (CWorkspace: need M*M+2*M, prefer M*M+M+M*NB)
2474: *                    (RWorkspace: 0)
2475: *
2476:                      CALL CGELQF( M, N, A, LDA, WORK( ITAU ),
2477:      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
2478: *
2479: *                    Copy L to WORK(IR), zeroing out above it
2480: *
2481:                      CALL CLACPY( 'L', M, M, A, LDA, WORK( IR ),
2482:      $                            LDWRKR )
2483:                      CALL CLASET( 'U', M-1, M-1, CZERO, CZERO,
2484:      $                            WORK( IR+LDWRKR ), LDWRKR )
2485: *
2486: *                    Generate Q in A
2487: *                    (CWorkspace: need M*M+2*M, prefer M*M+M+M*NB)
2488: *                    (RWorkspace: 0)
2489: *
2490:                      CALL CUNGLQ( M, N, M, A, LDA, WORK( ITAU ),
2491:      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
2492:                      IE = 1
2493:                      ITAUQ = ITAU
2494:                      ITAUP = ITAUQ + M
2495:                      IWORK = ITAUP + M
2496: *
2497: *                    Bidiagonalize L in WORK(IR)
2498: *                    (CWorkspace: need M*M+3*M, prefer M*M+2*M+2*M*NB)
2499: *                    (RWorkspace: need M)
2500: *
2501:                      CALL CGEBRD( M, M, WORK( IR ), LDWRKR, S,
2502:      $                            RWORK( IE ), WORK( ITAUQ ),
2503:      $                            WORK( ITAUP ), WORK( IWORK ),
2504:      $                            LWORK-IWORK+1, IERR )
2505: *
2506: *                    Generate right vectors bidiagonalizing L in
2507: *                    WORK(IR)
2508: *                    (CWorkspace: need M*M+3*M, prefer M*M+2*M+(M-1)*NB)
2509: *                    (RWorkspace: 0)
2510: *
2511:                      CALL CUNGBR( 'P', M, M, M, WORK( IR ), LDWRKR,
2512:      $                            WORK( ITAUP ), WORK( IWORK ),
2513:      $                            LWORK-IWORK+1, IERR )
2514:                      IRWORK = IE + M
2515: *
2516: *                    Perform bidiagonal QR iteration, computing right
2517: *                    singular vectors of L in WORK(IR)
2518: *                    (CWorkspace: need M*M)
2519: *                    (RWorkspace: need BDSPAC)
2520: *
2521:                      CALL CBDSQR( 'U', M, M, 0, 0, S, RWORK( IE ),
2522:      $                            WORK( IR ), LDWRKR, CDUM, 1, CDUM, 1,
2523:      $                            RWORK( IRWORK ), INFO )
2524: *
2525: *                    Multiply right singular vectors of L in WORK(IR) by
2526: *                    Q in A, storing result in VT
2527: *                    (CWorkspace: need M*M)
2528: *                    (RWorkspace: 0)
2529: *
2530:                      CALL CGEMM( 'N', 'N', M, N, M, CONE, WORK( IR ),
2531:      $                           LDWRKR, A, LDA, CZERO, VT, LDVT )
2532: *
2533:                   ELSE
2534: *
2535: *                    Insufficient workspace for a fast algorithm
2536: *
2537:                      ITAU = 1
2538:                      IWORK = ITAU + M
2539: *
2540: *                    Compute A=L*Q
2541: *                    (CWorkspace: need 2*M, prefer M+M*NB)
2542: *                    (RWorkspace: 0)
2543: *
2544:                      CALL CGELQF( M, N, A, LDA, WORK( ITAU ),
2545:      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
2546: *
2547: *                    Copy result to VT
2548: *
2549:                      CALL CLACPY( 'U', M, N, A, LDA, VT, LDVT )
2550: *
2551: *                    Generate Q in VT
2552: *                    (CWorkspace: need 2*M, prefer M+M*NB)
2553: *                    (RWorkspace: 0)
2554: *
2555:                      CALL CUNGLQ( M, N, M, VT, LDVT, WORK( ITAU ),
2556:      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
2557:                      IE = 1
2558:                      ITAUQ = ITAU
2559:                      ITAUP = ITAUQ + M
2560:                      IWORK = ITAUP + M
2561: *
2562: *                    Zero out above L in A
2563: *
2564:                      CALL CLASET( 'U', M-1, M-1, CZERO, CZERO,
2565:      $                            A( 1, 2 ), LDA )
2566: *
2567: *                    Bidiagonalize L in A
2568: *                    (CWorkspace: need 3*M, prefer 2*M+2*M*NB)
2569: *                    (RWorkspace: need M)
2570: *
2571:                      CALL CGEBRD( M, M, A, LDA, S, RWORK( IE ),
2572:      $                            WORK( ITAUQ ), WORK( ITAUP ),
2573:      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
2574: *
2575: *                    Multiply right vectors bidiagonalizing L by Q in VT
2576: *                    (CWorkspace: need 2*M+N, prefer 2*M+N*NB)
2577: *                    (RWorkspace: 0)
2578: *
2579:                      CALL CUNMBR( 'P', 'L', 'C', M, N, M, A, LDA,
2580:      $                            WORK( ITAUP ), VT, LDVT,
2581:      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
2582:                      IRWORK = IE + M
2583: *
2584: *                    Perform bidiagonal QR iteration, computing right
2585: *                    singular vectors of A in VT
2586: *                    (CWorkspace: 0)
2587: *                    (RWorkspace: need BDSPAC)
2588: *
2589:                      CALL CBDSQR( 'U', M, N, 0, 0, S, RWORK( IE ), VT,
2590:      $                            LDVT, CDUM, 1, CDUM, 1,
2591:      $                            RWORK( IRWORK ), INFO )
2592: *
2593:                   END IF
2594: *
2595:                ELSE IF( WNTUO ) THEN
2596: *
2597: *                 Path 5t(N much larger than M, JOBU='O', JOBVT='S')
2598: *                 M right singular vectors to be computed in VT and
2599: *                 M left singular vectors to be overwritten on A
2600: *
2601:                   IF( LWORK.GE.2*M*M+3*M ) THEN
2602: *
2603: *                    Sufficient workspace for a fast algorithm
2604: *
2605:                      IU = 1
2606:                      IF( LWORK.GE.WRKBL+2*LDA*M ) THEN
2607: *
2608: *                       WORK(IU) is LDA by M and WORK(IR) is LDA by M
2609: *
2610:                         LDWRKU = LDA
2611:                         IR = IU + LDWRKU*M
2612:                         LDWRKR = LDA
2613:                      ELSE IF( LWORK.GE.WRKBL+( LDA+M )*M ) THEN
2614: *
2615: *                       WORK(IU) is LDA by M and WORK(IR) is M by M
2616: *
2617:                         LDWRKU = LDA
2618:                         IR = IU + LDWRKU*M
2619:                         LDWRKR = M
2620:                      ELSE
2621: *
2622: *                       WORK(IU) is M by M and WORK(IR) is M by M
2623: *
2624:                         LDWRKU = M
2625:                         IR = IU + LDWRKU*M
2626:                         LDWRKR = M
2627:                      END IF
2628:                      ITAU = IR + LDWRKR*M
2629:                      IWORK = ITAU + M
2630: *
2631: *                    Compute A=L*Q
2632: *                    (CWorkspace: need 2*M*M+2*M, prefer 2*M*M+M+M*NB)
2633: *                    (RWorkspace: 0)
2634: *
2635:                      CALL CGELQF( M, N, A, LDA, WORK( ITAU ),
2636:      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
2637: *
2638: *                    Copy L to WORK(IU), zeroing out below it
2639: *
2640:                      CALL CLACPY( 'L', M, M, A, LDA, WORK( IU ),
2641:      $                            LDWRKU )
2642:                      CALL CLASET( 'U', M-1, M-1, CZERO, CZERO,
2643:      $                            WORK( IU+LDWRKU ), LDWRKU )
2644: *
2645: *                    Generate Q in A
2646: *                    (CWorkspace: need 2*M*M+2*M, prefer 2*M*M+M+M*NB)
2647: *                    (RWorkspace: 0)
2648: *
2649:                      CALL CUNGLQ( M, N, M, A, LDA, WORK( ITAU ),
2650:      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
2651:                      IE = 1
2652:                      ITAUQ = ITAU
2653:                      ITAUP = ITAUQ + M
2654:                      IWORK = ITAUP + M
2655: *
2656: *                    Bidiagonalize L in WORK(IU), copying result to
2657: *                    WORK(IR)
2658: *                    (CWorkspace: need   2*M*M+3*M,
2659: *                                 prefer 2*M*M+2*M+2*M*NB)
2660: *                    (RWorkspace: need   M)
2661: *
2662:                      CALL CGEBRD( M, M, WORK( IU ), LDWRKU, S,
2663:      $                            RWORK( IE ), WORK( ITAUQ ),
2664:      $                            WORK( ITAUP ), WORK( IWORK ),
2665:      $                            LWORK-IWORK+1, IERR )
2666:                      CALL CLACPY( 'L', M, M, WORK( IU ), LDWRKU,
2667:      $                            WORK( IR ), LDWRKR )
2668: *
2669: *                    Generate right bidiagonalizing vectors in WORK(IU)
2670: *                    (CWorkspace: need   2*M*M+3*M-1,
2671: *                                 prefer 2*M*M+2*M+(M-1)*NB)
2672: *                    (RWorkspace: 0)
2673: *
2674:                      CALL CUNGBR( 'P', M, M, M, WORK( IU ), LDWRKU,
2675:      $                            WORK( ITAUP ), WORK( IWORK ),
2676:      $                            LWORK-IWORK+1, IERR )
2677: *
2678: *                    Generate left bidiagonalizing vectors in WORK(IR)
2679: *                    (CWorkspace: need 2*M*M+3*M, prefer 2*M*M+2*M+M*NB)
2680: *                    (RWorkspace: 0)
2681: *
2682:                      CALL CUNGBR( 'Q', M, M, M, WORK( IR ), LDWRKR,
2683:      $                            WORK( ITAUQ ), WORK( IWORK ),
2684:      $                            LWORK-IWORK+1, IERR )
2685:                      IRWORK = IE + M
2686: *
2687: *                    Perform bidiagonal QR iteration, computing left
2688: *                    singular vectors of L in WORK(IR) and computing
2689: *                    right singular vectors of L in WORK(IU)
2690: *                    (CWorkspace: need 2*M*M)
2691: *                    (RWorkspace: need BDSPAC)
2692: *
2693:                      CALL CBDSQR( 'U', M, M, M, 0, S, RWORK( IE ),
2694:      $                            WORK( IU ), LDWRKU, WORK( IR ),
2695:      $                            LDWRKR, CDUM, 1, RWORK( IRWORK ),
2696:      $                            INFO )
2697: *
2698: *                    Multiply right singular vectors of L in WORK(IU) by
2699: *                    Q in A, storing result in VT
2700: *                    (CWorkspace: need M*M)
2701: *                    (RWorkspace: 0)
2702: *
2703:                      CALL CGEMM( 'N', 'N', M, N, M, CONE, WORK( IU ),
2704:      $                           LDWRKU, A, LDA, CZERO, VT, LDVT )
2705: *
2706: *                    Copy left singular vectors of L to A
2707: *                    (CWorkspace: need M*M)
2708: *                    (RWorkspace: 0)
2709: *
2710:                      CALL CLACPY( 'F', M, M, WORK( IR ), LDWRKR, A,
2711:      $                            LDA )
2712: *
2713:                   ELSE
2714: *
2715: *                    Insufficient workspace for a fast algorithm
2716: *
2717:                      ITAU = 1
2718:                      IWORK = ITAU + M
2719: *
2720: *                    Compute A=L*Q, copying result to VT
2721: *                    (CWorkspace: need 2*M, prefer M+M*NB)
2722: *                    (RWorkspace: 0)
2723: *
2724:                      CALL CGELQF( M, N, A, LDA, WORK( ITAU ),
2725:      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
2726:                      CALL CLACPY( 'U', M, N, A, LDA, VT, LDVT )
2727: *
2728: *                    Generate Q in VT
2729: *                    (CWorkspace: need 2*M, prefer M+M*NB)
2730: *                    (RWorkspace: 0)
2731: *
2732:                      CALL CUNGLQ( M, N, M, VT, LDVT, WORK( ITAU ),
2733:      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
2734:                      IE = 1
2735:                      ITAUQ = ITAU
2736:                      ITAUP = ITAUQ + M
2737:                      IWORK = ITAUP + M
2738: *
2739: *                    Zero out above L in A
2740: *
2741:                      CALL CLASET( 'U', M-1, M-1, CZERO, CZERO,
2742:      $                            A( 1, 2 ), LDA )
2743: *
2744: *                    Bidiagonalize L in A
2745: *                    (CWorkspace: need 3*M, prefer 2*M+2*M*NB)
2746: *                    (RWorkspace: need M)
2747: *
2748:                      CALL CGEBRD( M, M, A, LDA, S, RWORK( IE ),
2749:      $                            WORK( ITAUQ ), WORK( ITAUP ),
2750:      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
2751: *
2752: *                    Multiply right vectors bidiagonalizing L by Q in VT
2753: *                    (CWorkspace: need 2*M+N, prefer 2*M+N*NB)
2754: *                    (RWorkspace: 0)
2755: *
2756:                      CALL CUNMBR( 'P', 'L', 'C', M, N, M, A, LDA,
2757:      $                            WORK( ITAUP ), VT, LDVT,
2758:      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
2759: *
2760: *                    Generate left bidiagonalizing vectors of L in A
2761: *                    (CWorkspace: need 3*M, prefer 2*M+M*NB)
2762: *                    (RWorkspace: 0)
2763: *
2764:                      CALL CUNGBR( 'Q', M, M, M, A, LDA, WORK( ITAUQ ),
2765:      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
2766:                      IRWORK = IE + M
2767: *
2768: *                    Perform bidiagonal QR iteration, computing left
2769: *                    singular vectors of A in A and computing right
2770: *                    singular vectors of A in VT
2771: *                    (CWorkspace: 0)
2772: *                    (RWorkspace: need BDSPAC)
2773: *
2774:                      CALL CBDSQR( 'U', M, N, M, 0, S, RWORK( IE ), VT,
2775:      $                            LDVT, A, LDA, CDUM, 1,
2776:      $                            RWORK( IRWORK ), INFO )
2777: *
2778:                   END IF
2779: *
2780:                ELSE IF( WNTUAS ) THEN
2781: *
2782: *                 Path 6t(N much larger than M, JOBU='S' or 'A',
2783: *                         JOBVT='S')
2784: *                 M right singular vectors to be computed in VT and
2785: *                 M left singular vectors to be computed in U
2786: *
2787:                   IF( LWORK.GE.M*M+3*M ) THEN
2788: *
2789: *                    Sufficient workspace for a fast algorithm
2790: *
2791:                      IU = 1
2792:                      IF( LWORK.GE.WRKBL+LDA*M ) THEN
2793: *
2794: *                       WORK(IU) is LDA by N
2795: *
2796:                         LDWRKU = LDA
2797:                      ELSE
2798: *
2799: *                       WORK(IU) is LDA by M
2800: *
2801:                         LDWRKU = M
2802:                      END IF
2803:                      ITAU = IU + LDWRKU*M
2804:                      IWORK = ITAU + M
2805: *
2806: *                    Compute A=L*Q
2807: *                    (CWorkspace: need M*M+2*M, prefer M*M+M+M*NB)
2808: *                    (RWorkspace: 0)
2809: *
2810:                      CALL CGELQF( M, N, A, LDA, WORK( ITAU ),
2811:      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
2812: *
2813: *                    Copy L to WORK(IU), zeroing out above it
2814: *
2815:                      CALL CLACPY( 'L', M, M, A, LDA, WORK( IU ),
2816:      $                            LDWRKU )
2817:                      CALL CLASET( 'U', M-1, M-1, CZERO, CZERO,
2818:      $                            WORK( IU+LDWRKU ), LDWRKU )
2819: *
2820: *                    Generate Q in A
2821: *                    (CWorkspace: need M*M+2*M, prefer M*M+M+M*NB)
2822: *                    (RWorkspace: 0)
2823: *
2824:                      CALL CUNGLQ( M, N, M, A, LDA, WORK( ITAU ),
2825:      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
2826:                      IE = 1
2827:                      ITAUQ = ITAU
2828:                      ITAUP = ITAUQ + M
2829:                      IWORK = ITAUP + M
2830: *
2831: *                    Bidiagonalize L in WORK(IU), copying result to U
2832: *                    (CWorkspace: need M*M+3*M, prefer M*M+2*M+2*M*NB)
2833: *                    (RWorkspace: need M)
2834: *
2835:                      CALL CGEBRD( M, M, WORK( IU ), LDWRKU, S,
2836:      $                            RWORK( IE ), WORK( ITAUQ ),
2837:      $                            WORK( ITAUP ), WORK( IWORK ),
2838:      $                            LWORK-IWORK+1, IERR )
2839:                      CALL CLACPY( 'L', M, M, WORK( IU ), LDWRKU, U,
2840:      $                            LDU )
2841: *
2842: *                    Generate right bidiagonalizing vectors in WORK(IU)
2843: *                    (CWorkspace: need   M*M+3*M-1,
2844: *                                 prefer M*M+2*M+(M-1)*NB)
2845: *                    (RWorkspace: 0)
2846: *
2847:                      CALL CUNGBR( 'P', M, M, M, WORK( IU ), LDWRKU,
2848:      $                            WORK( ITAUP ), WORK( IWORK ),
2849:      $                            LWORK-IWORK+1, IERR )
2850: *
2851: *                    Generate left bidiagonalizing vectors in U
2852: *                    (CWorkspace: need M*M+3*M, prefer M*M+2*M+M*NB)
2853: *                    (RWorkspace: 0)
2854: *
2855:                      CALL CUNGBR( 'Q', M, M, M, U, LDU, WORK( ITAUQ ),
2856:      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
2857:                      IRWORK = IE + M
2858: *
2859: *                    Perform bidiagonal QR iteration, computing left
2860: *                    singular vectors of L in U and computing right
2861: *                    singular vectors of L in WORK(IU)
2862: *                    (CWorkspace: need M*M)
2863: *                    (RWorkspace: need BDSPAC)
2864: *
2865:                      CALL CBDSQR( 'U', M, M, M, 0, S, RWORK( IE ),
2866:      $                            WORK( IU ), LDWRKU, U, LDU, CDUM, 1,
2867:      $                            RWORK( IRWORK ), INFO )
2868: *
2869: *                    Multiply right singular vectors of L in WORK(IU) by
2870: *                    Q in A, storing result in VT
2871: *                    (CWorkspace: need M*M)
2872: *                    (RWorkspace: 0)
2873: *
2874:                      CALL CGEMM( 'N', 'N', M, N, M, CONE, WORK( IU ),
2875:      $                           LDWRKU, A, LDA, CZERO, VT, LDVT )
2876: *
2877:                   ELSE
2878: *
2879: *                    Insufficient workspace for a fast algorithm
2880: *
2881:                      ITAU = 1
2882:                      IWORK = ITAU + M
2883: *
2884: *                    Compute A=L*Q, copying result to VT
2885: *                    (CWorkspace: need 2*M, prefer M+M*NB)
2886: *                    (RWorkspace: 0)
2887: *
2888:                      CALL CGELQF( M, N, A, LDA, WORK( ITAU ),
2889:      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
2890:                      CALL CLACPY( 'U', M, N, A, LDA, VT, LDVT )
2891: *
2892: *                    Generate Q in VT
2893: *                    (CWorkspace: need 2*M, prefer M+M*NB)
2894: *                    (RWorkspace: 0)
2895: *
2896:                      CALL CUNGLQ( M, N, M, VT, LDVT, WORK( ITAU ),
2897:      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
2898: *
2899: *                    Copy L to U, zeroing out above it
2900: *
2901:                      CALL CLACPY( 'L', M, M, A, LDA, U, LDU )
2902:                      CALL CLASET( 'U', M-1, M-1, CZERO, CZERO,
2903:      $                            U( 1, 2 ), LDU )
2904:                      IE = 1
2905:                      ITAUQ = ITAU
2906:                      ITAUP = ITAUQ + M
2907:                      IWORK = ITAUP + M
2908: *
2909: *                    Bidiagonalize L in U
2910: *                    (CWorkspace: need 3*M, prefer 2*M+2*M*NB)
2911: *                    (RWorkspace: need M)
2912: *
2913:                      CALL CGEBRD( M, M, U, LDU, S, RWORK( IE ),
2914:      $                            WORK( ITAUQ ), WORK( ITAUP ),
2915:      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
2916: *
2917: *                    Multiply right bidiagonalizing vectors in U by Q
2918: *                    in VT
2919: *                    (CWorkspace: need 2*M+N, prefer 2*M+N*NB)
2920: *                    (RWorkspace: 0)
2921: *
2922:                      CALL CUNMBR( 'P', 'L', 'C', M, N, M, U, LDU,
2923:      $                            WORK( ITAUP ), VT, LDVT,
2924:      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
2925: *
2926: *                    Generate left bidiagonalizing vectors in U
2927: *                    (CWorkspace: need 3*M, prefer 2*M+M*NB)
2928: *                    (RWorkspace: 0)
2929: *
2930:                      CALL CUNGBR( 'Q', M, M, M, U, LDU, WORK( ITAUQ ),
2931:      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
2932:                      IRWORK = IE + M
2933: *
2934: *                    Perform bidiagonal QR iteration, computing left
2935: *                    singular vectors of A in U and computing right
2936: *                    singular vectors of A in VT
2937: *                    (CWorkspace: 0)
2938: *                    (RWorkspace: need BDSPAC)
2939: *
2940:                      CALL CBDSQR( 'U', M, N, M, 0, S, RWORK( IE ), VT,
2941:      $                            LDVT, U, LDU, CDUM, 1,
2942:      $                            RWORK( IRWORK ), INFO )
2943: *
2944:                   END IF
2945: *
2946:                END IF
2947: *
2948:             ELSE IF( WNTVA ) THEN
2949: *
2950:                IF( WNTUN ) THEN
2951: *
2952: *                 Path 7t(N much larger than M, JOBU='N', JOBVT='A')
2953: *                 N right singular vectors to be computed in VT and
2954: *                 no left singular vectors to be computed
2955: *
2956:                   IF( LWORK.GE.M*M+MAX( N+M, 3*M ) ) THEN
2957: *
2958: *                    Sufficient workspace for a fast algorithm
2959: *
2960:                      IR = 1
2961:                      IF( LWORK.GE.WRKBL+LDA*M ) THEN
2962: *
2963: *                       WORK(IR) is LDA by M
2964: *
2965:                         LDWRKR = LDA
2966:                      ELSE
2967: *
2968: *                       WORK(IR) is M by M
2969: *
2970:                         LDWRKR = M
2971:                      END IF
2972:                      ITAU = IR + LDWRKR*M
2973:                      IWORK = ITAU + M
2974: *
2975: *                    Compute A=L*Q, copying result to VT
2976: *                    (CWorkspace: need M*M+2*M, prefer M*M+M+M*NB)
2977: *                    (RWorkspace: 0)
2978: *
2979:                      CALL CGELQF( M, N, A, LDA, WORK( ITAU ),
2980:      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
2981:                      CALL CLACPY( 'U', M, N, A, LDA, VT, LDVT )
2982: *
2983: *                    Copy L to WORK(IR), zeroing out above it
2984: *
2985:                      CALL CLACPY( 'L', M, M, A, LDA, WORK( IR ),
2986:      $                            LDWRKR )
2987:                      CALL CLASET( 'U', M-1, M-1, CZERO, CZERO,
2988:      $                            WORK( IR+LDWRKR ), LDWRKR )
2989: *
2990: *                    Generate Q in VT
2991: *                    (CWorkspace: need M*M+M+N, prefer M*M+M+N*NB)
2992: *                    (RWorkspace: 0)
2993: *
2994:                      CALL CUNGLQ( N, N, M, VT, LDVT, WORK( ITAU ),
2995:      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
2996:                      IE = 1
2997:                      ITAUQ = ITAU
2998:                      ITAUP = ITAUQ + M
2999:                      IWORK = ITAUP + M
3000: *
3001: *                    Bidiagonalize L in WORK(IR)
3002: *                    (CWorkspace: need M*M+3*M, prefer M*M+2*M+2*M*NB)
3003: *                    (RWorkspace: need M)
3004: *
3005:                      CALL CGEBRD( M, M, WORK( IR ), LDWRKR, S,
3006:      $                            RWORK( IE ), WORK( ITAUQ ),
3007:      $                            WORK( ITAUP ), WORK( IWORK ),
3008:      $                            LWORK-IWORK+1, IERR )
3009: *
3010: *                    Generate right bidiagonalizing vectors in WORK(IR)
3011: *                    (CWorkspace: need   M*M+3*M-1,
3012: *                                 prefer M*M+2*M+(M-1)*NB)
3013: *                    (RWorkspace: 0)
3014: *
3015:                      CALL CUNGBR( 'P', M, M, M, WORK( IR ), LDWRKR,
3016:      $                            WORK( ITAUP ), WORK( IWORK ),
3017:      $                            LWORK-IWORK+1, IERR )
3018:                      IRWORK = IE + M
3019: *
3020: *                    Perform bidiagonal QR iteration, computing right
3021: *                    singular vectors of L in WORK(IR)
3022: *                    (CWorkspace: need M*M)
3023: *                    (RWorkspace: need BDSPAC)
3024: *
3025:                      CALL CBDSQR( 'U', M, M, 0, 0, S, RWORK( IE ),
3026:      $                            WORK( IR ), LDWRKR, CDUM, 1, CDUM, 1,
3027:      $                            RWORK( IRWORK ), INFO )
3028: *
3029: *                    Multiply right singular vectors of L in WORK(IR) by
3030: *                    Q in VT, storing result in A
3031: *                    (CWorkspace: need M*M)
3032: *                    (RWorkspace: 0)
3033: *
3034:                      CALL CGEMM( 'N', 'N', M, N, M, CONE, WORK( IR ),
3035:      $                           LDWRKR, VT, LDVT, CZERO, A, LDA )
3036: *
3037: *                    Copy right singular vectors of A from A to VT
3038: *
3039:                      CALL CLACPY( 'F', M, N, A, LDA, VT, LDVT )
3040: *
3041:                   ELSE
3042: *
3043: *                    Insufficient workspace for a fast algorithm
3044: *
3045:                      ITAU = 1
3046:                      IWORK = ITAU + M
3047: *
3048: *                    Compute A=L*Q, copying result to VT
3049: *                    (CWorkspace: need 2*M, prefer M+M*NB)
3050: *                    (RWorkspace: 0)
3051: *
3052:                      CALL CGELQF( M, N, A, LDA, WORK( ITAU ),
3053:      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
3054:                      CALL CLACPY( 'U', M, N, A, LDA, VT, LDVT )
3055: *
3056: *                    Generate Q in VT
3057: *                    (CWorkspace: need M+N, prefer M+N*NB)
3058: *                    (RWorkspace: 0)
3059: *
3060:                      CALL CUNGLQ( N, N, M, VT, LDVT, WORK( ITAU ),
3061:      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
3062:                      IE = 1
3063:                      ITAUQ = ITAU
3064:                      ITAUP = ITAUQ + M
3065:                      IWORK = ITAUP + M
3066: *
3067: *                    Zero out above L in A
3068: *
3069:                      CALL CLASET( 'U', M-1, M-1, CZERO, CZERO,
3070:      $                            A( 1, 2 ), LDA )
3071: *
3072: *                    Bidiagonalize L in A
3073: *                    (CWorkspace: need 3*M, prefer 2*M+2*M*NB)
3074: *                    (RWorkspace: need M)
3075: *
3076:                      CALL CGEBRD( M, M, A, LDA, S, RWORK( IE ),
3077:      $                            WORK( ITAUQ ), WORK( ITAUP ),
3078:      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
3079: *
3080: *                    Multiply right bidiagonalizing vectors in A by Q
3081: *                    in VT
3082: *                    (CWorkspace: need 2*M+N, prefer 2*M+N*NB)
3083: *                    (RWorkspace: 0)
3084: *
3085:                      CALL CUNMBR( 'P', 'L', 'C', M, N, M, A, LDA,
3086:      $                            WORK( ITAUP ), VT, LDVT,
3087:      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
3088:                      IRWORK = IE + M
3089: *
3090: *                    Perform bidiagonal QR iteration, computing right
3091: *                    singular vectors of A in VT
3092: *                    (CWorkspace: 0)
3093: *                    (RWorkspace: need BDSPAC)
3094: *
3095:                      CALL CBDSQR( 'U', M, N, 0, 0, S, RWORK( IE ), VT,
3096:      $                            LDVT, CDUM, 1, CDUM, 1,
3097:      $                            RWORK( IRWORK ), INFO )
3098: *
3099:                   END IF
3100: *
3101:                ELSE IF( WNTUO ) THEN
3102: *
3103: *                 Path 8t(N much larger than M, JOBU='O', JOBVT='A')
3104: *                 N right singular vectors to be computed in VT and
3105: *                 M left singular vectors to be overwritten on A
3106: *
3107:                   IF( LWORK.GE.2*M*M+MAX( N+M, 3*M ) ) THEN
3108: *
3109: *                    Sufficient workspace for a fast algorithm
3110: *
3111:                      IU = 1
3112:                      IF( LWORK.GE.WRKBL+2*LDA*M ) THEN
3113: *
3114: *                       WORK(IU) is LDA by M and WORK(IR) is LDA by M
3115: *
3116:                         LDWRKU = LDA
3117:                         IR = IU + LDWRKU*M
3118:                         LDWRKR = LDA
3119:                      ELSE IF( LWORK.GE.WRKBL+( LDA+M )*M ) THEN
3120: *
3121: *                       WORK(IU) is LDA by M and WORK(IR) is M by M
3122: *
3123:                         LDWRKU = LDA
3124:                         IR = IU + LDWRKU*M
3125:                         LDWRKR = M
3126:                      ELSE
3127: *
3128: *                       WORK(IU) is M by M and WORK(IR) is M by M
3129: *
3130:                         LDWRKU = M
3131:                         IR = IU + LDWRKU*M
3132:                         LDWRKR = M
3133:                      END IF
3134:                      ITAU = IR + LDWRKR*M
3135:                      IWORK = ITAU + M
3136: *
3137: *                    Compute A=L*Q, copying result to VT
3138: *                    (CWorkspace: need 2*M*M+2*M, prefer 2*M*M+M+M*NB)
3139: *                    (RWorkspace: 0)
3140: *
3141:                      CALL CGELQF( M, N, A, LDA, WORK( ITAU ),
3142:      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
3143:                      CALL CLACPY( 'U', M, N, A, LDA, VT, LDVT )
3144: *
3145: *                    Generate Q in VT
3146: *                    (CWorkspace: need 2*M*M+M+N, prefer 2*M*M+M+N*NB)
3147: *                    (RWorkspace: 0)
3148: *
3149:                      CALL CUNGLQ( N, N, M, VT, LDVT, WORK( ITAU ),
3150:      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
3151: *
3152: *                    Copy L to WORK(IU), zeroing out above it
3153: *
3154:                      CALL CLACPY( 'L', M, M, A, LDA, WORK( IU ),
3155:      $                            LDWRKU )
3156:                      CALL CLASET( 'U', M-1, M-1, CZERO, CZERO,
3157:      $                            WORK( IU+LDWRKU ), LDWRKU )
3158:                      IE = 1
3159:                      ITAUQ = ITAU
3160:                      ITAUP = ITAUQ + M
3161:                      IWORK = ITAUP + M
3162: *
3163: *                    Bidiagonalize L in WORK(IU), copying result to
3164: *                    WORK(IR)
3165: *                    (CWorkspace: need   2*M*M+3*M,
3166: *                                 prefer 2*M*M+2*M+2*M*NB)
3167: *                    (RWorkspace: need   M)
3168: *
3169:                      CALL CGEBRD( M, M, WORK( IU ), LDWRKU, S,
3170:      $                            RWORK( IE ), WORK( ITAUQ ),
3171:      $                            WORK( ITAUP ), WORK( IWORK ),
3172:      $                            LWORK-IWORK+1, IERR )
3173:                      CALL CLACPY( 'L', M, M, WORK( IU ), LDWRKU,
3174:      $                            WORK( IR ), LDWRKR )
3175: *
3176: *                    Generate right bidiagonalizing vectors in WORK(IU)
3177: *                    (CWorkspace: need   2*M*M+3*M-1,
3178: *                                 prefer 2*M*M+2*M+(M-1)*NB)
3179: *                    (RWorkspace: 0)
3180: *
3181:                      CALL CUNGBR( 'P', M, M, M, WORK( IU ), LDWRKU,
3182:      $                            WORK( ITAUP ), WORK( IWORK ),
3183:      $                            LWORK-IWORK+1, IERR )
3184: *
3185: *                    Generate left bidiagonalizing vectors in WORK(IR)
3186: *                    (CWorkspace: need 2*M*M+3*M, prefer 2*M*M+2*M+M*NB)
3187: *                    (RWorkspace: 0)
3188: *
3189:                      CALL CUNGBR( 'Q', M, M, M, WORK( IR ), LDWRKR,
3190:      $                            WORK( ITAUQ ), WORK( IWORK ),
3191:      $                            LWORK-IWORK+1, IERR )
3192:                      IRWORK = IE + M
3193: *
3194: *                    Perform bidiagonal QR iteration, computing left
3195: *                    singular vectors of L in WORK(IR) and computing
3196: *                    right singular vectors of L in WORK(IU)
3197: *                    (CWorkspace: need 2*M*M)
3198: *                    (RWorkspace: need BDSPAC)
3199: *
3200:                      CALL CBDSQR( 'U', M, M, M, 0, S, RWORK( IE ),
3201:      $                            WORK( IU ), LDWRKU, WORK( IR ),
3202:      $                            LDWRKR, CDUM, 1, RWORK( IRWORK ),
3203:      $                            INFO )
3204: *
3205: *                    Multiply right singular vectors of L in WORK(IU) by
3206: *                    Q in VT, storing result in A
3207: *                    (CWorkspace: need M*M)
3208: *                    (RWorkspace: 0)
3209: *
3210:                      CALL CGEMM( 'N', 'N', M, N, M, CONE, WORK( IU ),
3211:      $                           LDWRKU, VT, LDVT, CZERO, A, LDA )
3212: *
3213: *                    Copy right singular vectors of A from A to VT
3214: *
3215:                      CALL CLACPY( 'F', M, N, A, LDA, VT, LDVT )
3216: *
3217: *                    Copy left singular vectors of A from WORK(IR) to A
3218: *
3219:                      CALL CLACPY( 'F', M, M, WORK( IR ), LDWRKR, A,
3220:      $                            LDA )
3221: *
3222:                   ELSE
3223: *
3224: *                    Insufficient workspace for a fast algorithm
3225: *
3226:                      ITAU = 1
3227:                      IWORK = ITAU + M
3228: *
3229: *                    Compute A=L*Q, copying result to VT
3230: *                    (CWorkspace: need 2*M, prefer M+M*NB)
3231: *                    (RWorkspace: 0)
3232: *
3233:                      CALL CGELQF( M, N, A, LDA, WORK( ITAU ),
3234:      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
3235:                      CALL CLACPY( 'U', M, N, A, LDA, VT, LDVT )
3236: *
3237: *                    Generate Q in VT
3238: *                    (CWorkspace: need M+N, prefer M+N*NB)
3239: *                    (RWorkspace: 0)
3240: *
3241:                      CALL CUNGLQ( N, N, M, VT, LDVT, WORK( ITAU ),
3242:      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
3243:                      IE = 1
3244:                      ITAUQ = ITAU
3245:                      ITAUP = ITAUQ + M
3246:                      IWORK = ITAUP + M
3247: *
3248: *                    Zero out above L in A
3249: *
3250:                      CALL CLASET( 'U', M-1, M-1, CZERO, CZERO,
3251:      $                            A( 1, 2 ), LDA )
3252: *
3253: *                    Bidiagonalize L in A
3254: *                    (CWorkspace: need 3*M, prefer 2*M+2*M*NB)
3255: *                    (RWorkspace: need M)
3256: *
3257:                      CALL CGEBRD( M, M, A, LDA, S, RWORK( IE ),
3258:      $                            WORK( ITAUQ ), WORK( ITAUP ),
3259:      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
3260: *
3261: *                    Multiply right bidiagonalizing vectors in A by Q
3262: *                    in VT
3263: *                    (CWorkspace: need 2*M+N, prefer 2*M+N*NB)
3264: *                    (RWorkspace: 0)
3265: *
3266:                      CALL CUNMBR( 'P', 'L', 'C', M, N, M, A, LDA,
3267:      $                            WORK( ITAUP ), VT, LDVT,
3268:      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
3269: *
3270: *                    Generate left bidiagonalizing vectors in A
3271: *                    (CWorkspace: need 3*M, prefer 2*M+M*NB)
3272: *                    (RWorkspace: 0)
3273: *
3274:                      CALL CUNGBR( 'Q', M, M, M, A, LDA, WORK( ITAUQ ),
3275:      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
3276:                      IRWORK = IE + M
3277: *
3278: *                    Perform bidiagonal QR iteration, computing left
3279: *                    singular vectors of A in A and computing right
3280: *                    singular vectors of A in VT
3281: *                    (CWorkspace: 0)
3282: *                    (RWorkspace: need BDSPAC)
3283: *
3284:                      CALL CBDSQR( 'U', M, N, M, 0, S, RWORK( IE ), VT,
3285:      $                            LDVT, A, LDA, CDUM, 1,
3286:      $                            RWORK( IRWORK ), INFO )
3287: *
3288:                   END IF
3289: *
3290:                ELSE IF( WNTUAS ) THEN
3291: *
3292: *                 Path 9t(N much larger than M, JOBU='S' or 'A',
3293: *                         JOBVT='A')
3294: *                 N right singular vectors to be computed in VT and
3295: *                 M left singular vectors to be computed in U
3296: *
3297:                   IF( LWORK.GE.M*M+MAX( N+M, 3*M ) ) THEN
3298: *
3299: *                    Sufficient workspace for a fast algorithm
3300: *
3301:                      IU = 1
3302:                      IF( LWORK.GE.WRKBL+LDA*M ) THEN
3303: *
3304: *                       WORK(IU) is LDA by M
3305: *
3306:                         LDWRKU = LDA
3307:                      ELSE
3308: *
3309: *                       WORK(IU) is M by M
3310: *
3311:                         LDWRKU = M
3312:                      END IF
3313:                      ITAU = IU + LDWRKU*M
3314:                      IWORK = ITAU + M
3315: *
3316: *                    Compute A=L*Q, copying result to VT
3317: *                    (CWorkspace: need M*M+2*M, prefer M*M+M+M*NB)
3318: *                    (RWorkspace: 0)
3319: *
3320:                      CALL CGELQF( M, N, A, LDA, WORK( ITAU ),
3321:      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
3322:                      CALL CLACPY( 'U', M, N, A, LDA, VT, LDVT )
3323: *
3324: *                    Generate Q in VT
3325: *                    (CWorkspace: need M*M+M+N, prefer M*M+M+N*NB)
3326: *                    (RWorkspace: 0)
3327: *
3328:                      CALL CUNGLQ( N, N, M, VT, LDVT, WORK( ITAU ),
3329:      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
3330: *
3331: *                    Copy L to WORK(IU), zeroing out above it
3332: *
3333:                      CALL CLACPY( 'L', M, M, A, LDA, WORK( IU ),
3334:      $                            LDWRKU )
3335:                      CALL CLASET( 'U', M-1, M-1, CZERO, CZERO,
3336:      $                            WORK( IU+LDWRKU ), LDWRKU )
3337:                      IE = 1
3338:                      ITAUQ = ITAU
3339:                      ITAUP = ITAUQ + M
3340:                      IWORK = ITAUP + M
3341: *
3342: *                    Bidiagonalize L in WORK(IU), copying result to U
3343: *                    (CWorkspace: need M*M+3*M, prefer M*M+2*M+2*M*NB)
3344: *                    (RWorkspace: need M)
3345: *
3346:                      CALL CGEBRD( M, M, WORK( IU ), LDWRKU, S,
3347:      $                            RWORK( IE ), WORK( ITAUQ ),
3348:      $                            WORK( ITAUP ), WORK( IWORK ),
3349:      $                            LWORK-IWORK+1, IERR )
3350:                      CALL CLACPY( 'L', M, M, WORK( IU ), LDWRKU, U,
3351:      $                            LDU )
3352: *
3353: *                    Generate right bidiagonalizing vectors in WORK(IU)
3354: *                    (CWorkspace: need M*M+3*M, prefer M*M+2*M+(M-1)*NB)
3355: *                    (RWorkspace: 0)
3356: *
3357:                      CALL CUNGBR( 'P', M, M, M, WORK( IU ), LDWRKU,
3358:      $                            WORK( ITAUP ), WORK( IWORK ),
3359:      $                            LWORK-IWORK+1, IERR )
3360: *
3361: *                    Generate left bidiagonalizing vectors in U
3362: *                    (CWorkspace: need M*M+3*M, prefer M*M+2*M+M*NB)
3363: *                    (RWorkspace: 0)
3364: *
3365:                      CALL CUNGBR( 'Q', M, M, M, U, LDU, WORK( ITAUQ ),
3366:      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
3367:                      IRWORK = IE + M
3368: *
3369: *                    Perform bidiagonal QR iteration, computing left
3370: *                    singular vectors of L in U and computing right
3371: *                    singular vectors of L in WORK(IU)
3372: *                    (CWorkspace: need M*M)
3373: *                    (RWorkspace: need BDSPAC)
3374: *
3375:                      CALL CBDSQR( 'U', M, M, M, 0, S, RWORK( IE ),
3376:      $                            WORK( IU ), LDWRKU, U, LDU, CDUM, 1,
3377:      $                            RWORK( IRWORK ), INFO )
3378: *
3379: *                    Multiply right singular vectors of L in WORK(IU) by
3380: *                    Q in VT, storing result in A
3381: *                    (CWorkspace: need M*M)
3382: *                    (RWorkspace: 0)
3383: *
3384:                      CALL CGEMM( 'N', 'N', M, N, M, CONE, WORK( IU ),
3385:      $                           LDWRKU, VT, LDVT, CZERO, A, LDA )
3386: *
3387: *                    Copy right singular vectors of A from A to VT
3388: *
3389:                      CALL CLACPY( 'F', M, N, A, LDA, VT, LDVT )
3390: *
3391:                   ELSE
3392: *
3393: *                    Insufficient workspace for a fast algorithm
3394: *
3395:                      ITAU = 1
3396:                      IWORK = ITAU + M
3397: *
3398: *                    Compute A=L*Q, copying result to VT
3399: *                    (CWorkspace: need 2*M, prefer M+M*NB)
3400: *                    (RWorkspace: 0)
3401: *
3402:                      CALL CGELQF( M, N, A, LDA, WORK( ITAU ),
3403:      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
3404:                      CALL CLACPY( 'U', M, N, A, LDA, VT, LDVT )
3405: *
3406: *                    Generate Q in VT
3407: *                    (CWorkspace: need M+N, prefer M+N*NB)
3408: *                    (RWorkspace: 0)
3409: *
3410:                      CALL CUNGLQ( N, N, M, VT, LDVT, WORK( ITAU ),
3411:      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
3412: *
3413: *                    Copy L to U, zeroing out above it
3414: *
3415:                      CALL CLACPY( 'L', M, M, A, LDA, U, LDU )
3416:                      CALL CLASET( 'U', M-1, M-1, CZERO, CZERO,
3417:      $                            U( 1, 2 ), LDU )
3418:                      IE = 1
3419:                      ITAUQ = ITAU
3420:                      ITAUP = ITAUQ + M
3421:                      IWORK = ITAUP + M
3422: *
3423: *                    Bidiagonalize L in U
3424: *                    (CWorkspace: need 3*M, prefer 2*M+2*M*NB)
3425: *                    (RWorkspace: need M)
3426: *
3427:                      CALL CGEBRD( M, M, U, LDU, S, RWORK( IE ),
3428:      $                            WORK( ITAUQ ), WORK( ITAUP ),
3429:      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
3430: *
3431: *                    Multiply right bidiagonalizing vectors in U by Q
3432: *                    in VT
3433: *                    (CWorkspace: need 2*M+N, prefer 2*M+N*NB)
3434: *                    (RWorkspace: 0)
3435: *
3436:                      CALL CUNMBR( 'P', 'L', 'C', M, N, M, U, LDU,
3437:      $                            WORK( ITAUP ), VT, LDVT,
3438:      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
3439: *
3440: *                    Generate left bidiagonalizing vectors in U
3441: *                    (CWorkspace: need 3*M, prefer 2*M+M*NB)
3442: *                    (RWorkspace: 0)
3443: *
3444:                      CALL CUNGBR( 'Q', M, M, M, U, LDU, WORK( ITAUQ ),
3445:      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
3446:                      IRWORK = IE + M
3447: *
3448: *                    Perform bidiagonal QR iteration, computing left
3449: *                    singular vectors of A in U and computing right
3450: *                    singular vectors of A in VT
3451: *                    (CWorkspace: 0)
3452: *                    (RWorkspace: need BDSPAC)
3453: *
3454:                      CALL CBDSQR( 'U', M, N, M, 0, S, RWORK( IE ), VT,
3455:      $                            LDVT, U, LDU, CDUM, 1,
3456:      $                            RWORK( IRWORK ), INFO )
3457: *
3458:                   END IF
3459: *
3460:                END IF
3461: *
3462:             END IF
3463: *
3464:          ELSE
3465: *
3466: *           N .LT. MNTHR
3467: *
3468: *           Path 10t(N greater than M, but not much larger)
3469: *           Reduce to bidiagonal form without LQ decomposition
3470: *
3471:             IE = 1
3472:             ITAUQ = 1
3473:             ITAUP = ITAUQ + M
3474:             IWORK = ITAUP + M
3475: *
3476: *           Bidiagonalize A
3477: *           (CWorkspace: need 2*M+N, prefer 2*M+(M+N)*NB)
3478: *           (RWorkspace: M)
3479: *
3480:             CALL CGEBRD( M, N, A, LDA, S, RWORK( IE ), WORK( ITAUQ ),
3481:      $                   WORK( ITAUP ), WORK( IWORK ), LWORK-IWORK+1,
3482:      $                   IERR )
3483:             IF( WNTUAS ) THEN
3484: *
3485: *              If left singular vectors desired in U, copy result to U
3486: *              and generate left bidiagonalizing vectors in U
3487: *              (CWorkspace: need 3*M-1, prefer 2*M+(M-1)*NB)
3488: *              (RWorkspace: 0)
3489: *
3490:                CALL CLACPY( 'L', M, M, A, LDA, U, LDU )
3491:                CALL CUNGBR( 'Q', M, M, N, U, LDU, WORK( ITAUQ ),
3492:      $                      WORK( IWORK ), LWORK-IWORK+1, IERR )
3493:             END IF
3494:             IF( WNTVAS ) THEN
3495: *
3496: *              If right singular vectors desired in VT, copy result to
3497: *              VT and generate right bidiagonalizing vectors in VT
3498: *              (CWorkspace: need 2*M+NRVT, prefer 2*M+NRVT*NB)
3499: *              (RWorkspace: 0)
3500: *
3501:                CALL CLACPY( 'U', M, N, A, LDA, VT, LDVT )
3502:                IF( WNTVA )
3503:      $            NRVT = N
3504:                IF( WNTVS )
3505:      $            NRVT = M
3506:                CALL CUNGBR( 'P', NRVT, N, M, VT, LDVT, WORK( ITAUP ),
3507:      $                      WORK( IWORK ), LWORK-IWORK+1, IERR )
3508:             END IF
3509:             IF( WNTUO ) THEN
3510: *
3511: *              If left singular vectors desired in A, generate left
3512: *              bidiagonalizing vectors in A
3513: *              (CWorkspace: need 3*M-1, prefer 2*M+(M-1)*NB)
3514: *              (RWorkspace: 0)
3515: *
3516:                CALL CUNGBR( 'Q', M, M, N, A, LDA, WORK( ITAUQ ),
3517:      $                      WORK( IWORK ), LWORK-IWORK+1, IERR )
3518:             END IF
3519:             IF( WNTVO ) THEN
3520: *
3521: *              If right singular vectors desired in A, generate right
3522: *              bidiagonalizing vectors in A
3523: *              (CWorkspace: need 3*M, prefer 2*M+M*NB)
3524: *              (RWorkspace: 0)
3525: *
3526:                CALL CUNGBR( 'P', M, N, M, A, LDA, WORK( ITAUP ),
3527:      $                      WORK( IWORK ), LWORK-IWORK+1, IERR )
3528:             END IF
3529:             IRWORK = IE + M
3530:             IF( WNTUAS .OR. WNTUO )
3531:      $         NRU = M
3532:             IF( WNTUN )
3533:      $         NRU = 0
3534:             IF( WNTVAS .OR. WNTVO )
3535:      $         NCVT = N
3536:             IF( WNTVN )
3537:      $         NCVT = 0
3538:             IF( ( .NOT.WNTUO ) .AND. ( .NOT.WNTVO ) ) THEN
3539: *
3540: *              Perform bidiagonal QR iteration, if desired, computing
3541: *              left singular vectors in U and computing right singular
3542: *              vectors in VT
3543: *              (CWorkspace: 0)
3544: *              (RWorkspace: need BDSPAC)
3545: *
3546:                CALL CBDSQR( 'L', M, NCVT, NRU, 0, S, RWORK( IE ), VT,
3547:      $                      LDVT, U, LDU, CDUM, 1, RWORK( IRWORK ),
3548:      $                      INFO )
3549:             ELSE IF( ( .NOT.WNTUO ) .AND. WNTVO ) THEN
3550: *
3551: *              Perform bidiagonal QR iteration, if desired, computing
3552: *              left singular vectors in U and computing right singular
3553: *              vectors in A
3554: *              (CWorkspace: 0)
3555: *              (RWorkspace: need BDSPAC)
3556: *
3557:                CALL CBDSQR( 'L', M, NCVT, NRU, 0, S, RWORK( IE ), A,
3558:      $                      LDA, U, LDU, CDUM, 1, RWORK( IRWORK ),
3559:      $                      INFO )
3560:             ELSE
3561: *
3562: *              Perform bidiagonal QR iteration, if desired, computing
3563: *              left singular vectors in A and computing right singular
3564: *              vectors in VT
3565: *              (CWorkspace: 0)
3566: *              (RWorkspace: need BDSPAC)
3567: *
3568:                CALL CBDSQR( 'L', M, NCVT, NRU, 0, S, RWORK( IE ), VT,
3569:      $                      LDVT, A, LDA, CDUM, 1, RWORK( IRWORK ),
3570:      $                      INFO )
3571:             END IF
3572: *
3573:          END IF
3574: *
3575:       END IF
3576: *
3577: *     Undo scaling if necessary
3578: *
3579:       IF( ISCL.EQ.1 ) THEN
3580:          IF( ANRM.GT.BIGNUM )
3581:      $      CALL SLASCL( 'G', 0, 0, BIGNUM, ANRM, MINMN, 1, S, MINMN,
3582:      $                   IERR )
3583:          IF( INFO.NE.0 .AND. ANRM.GT.BIGNUM )
3584:      $      CALL SLASCL( 'G', 0, 0, BIGNUM, ANRM, MINMN-1, 1,
3585:      $                   RWORK( IE ), MINMN, IERR )
3586:          IF( ANRM.LT.SMLNUM )
3587:      $      CALL SLASCL( 'G', 0, 0, SMLNUM, ANRM, MINMN, 1, S, MINMN,
3588:      $                   IERR )
3589:          IF( INFO.NE.0 .AND. ANRM.LT.SMLNUM )
3590:      $      CALL SLASCL( 'G', 0, 0, SMLNUM, ANRM, MINMN-1, 1,
3591:      $                   RWORK( IE ), MINMN, IERR )
3592:       END IF
3593: *
3594: *     Return optimal workspace in WORK(1)
3595: *
3596:       WORK( 1 ) = MAXWRK
3597: *
3598:       RETURN
3599: *
3600: *     End of CGESVD
3601: *
3602:       END
3603: