0001:       SUBROUTINE CGESDD( JOBZ, M, N, A, LDA, S, U, LDU, VT, LDVT,
0002:      $                   WORK, LWORK, RWORK, IWORK, INFO )
0003: *
0004: *  -- LAPACK driver routine (version 3.2) --
0005: *     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
0006: *     November 2006
0007: *     8-15-00:  Improve consistency of WS calculations (eca)
0008: *
0009: *     .. Scalar Arguments ..
0010:       CHARACTER          JOBZ
0011:       INTEGER            INFO, LDA, LDU, LDVT, LWORK, M, N
0012: *     ..
0013: *     .. Array Arguments ..
0014:       INTEGER            IWORK( * )
0015:       REAL               RWORK( * ), S( * )
0016:       COMPLEX            A( LDA, * ), U( LDU, * ), VT( LDVT, * ),
0017:      $                   WORK( * )
0018: *     ..
0019: *
0020: *  Purpose
0021: *  =======
0022: *
0023: *  CGESDD computes the singular value decomposition (SVD) of a complex
0024: *  M-by-N matrix A, optionally computing the left and/or right singular
0025: *  vectors, by using divide-and-conquer method. The SVD is written
0026: *
0027: *       A = U * SIGMA * conjugate-transpose(V)
0028: *
0029: *  where SIGMA is an M-by-N matrix which is zero except for its
0030: *  min(m,n) diagonal elements, U is an M-by-M unitary matrix, and
0031: *  V is an N-by-N unitary matrix.  The diagonal elements of SIGMA
0032: *  are the singular values of A; they are real and non-negative, and
0033: *  are returned in descending order.  The first min(m,n) columns of
0034: *  U and V are the left and right singular vectors of A.
0035: *
0036: *  Note that the routine returns VT = V**H, not V.
0037: *
0038: *  The divide and conquer algorithm makes very mild assumptions about
0039: *  floating point arithmetic. It will work on machines with a guard
0040: *  digit in add/subtract, or on those binary machines without guard
0041: *  digits which subtract like the Cray X-MP, Cray Y-MP, Cray C-90, or
0042: *  Cray-2. It could conceivably fail on hexadecimal or decimal machines
0043: *  without guard digits, but we know of none.
0044: *
0045: *  Arguments
0046: *  =========
0047: *
0048: *  JOBZ    (input) CHARACTER*1
0049: *          Specifies options for computing all or part of the matrix U:
0050: *          = 'A':  all M columns of U and all N rows of V**H are
0051: *                  returned in the arrays U and VT;
0052: *          = 'S':  the first min(M,N) columns of U and the first
0053: *                  min(M,N) rows of V**H are returned in the arrays U
0054: *                  and VT;
0055: *          = 'O':  If M >= N, the first N columns of U are overwritten
0056: *                  in the array A and all rows of V**H are returned in
0057: *                  the array VT;
0058: *                  otherwise, all columns of U are returned in the
0059: *                  array U and the first M rows of V**H are overwritten
0060: *                  in the array A;
0061: *          = 'N':  no columns of U or rows of V**H are computed.
0062: *
0063: *  M       (input) INTEGER
0064: *          The number of rows of the input matrix A.  M >= 0.
0065: *
0066: *  N       (input) INTEGER
0067: *          The number of columns of the input matrix A.  N >= 0.
0068: *
0069: *  A       (input/output) COMPLEX array, dimension (LDA,N)
0070: *          On entry, the M-by-N matrix A.
0071: *          On exit,
0072: *          if JOBZ = 'O',  A is overwritten with the first N columns
0073: *                          of U (the left singular vectors, stored
0074: *                          columnwise) if M >= N;
0075: *                          A is overwritten with the first M rows
0076: *                          of V**H (the right singular vectors, stored
0077: *                          rowwise) otherwise.
0078: *          if JOBZ .ne. 'O', the contents of A 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: *          UCOL = M if JOBZ = 'A' or JOBZ = 'O' and M < N;
0088: *          UCOL = min(M,N) if JOBZ = 'S'.
0089: *          If JOBZ = 'A' or JOBZ = 'O' and M < N, U contains the M-by-M
0090: *          unitary matrix U;
0091: *          if JOBZ = 'S', U contains the first min(M,N) columns of U
0092: *          (the left singular vectors, stored columnwise);
0093: *          if JOBZ = 'O' and M >= N, or JOBZ = 'N', U is not referenced.
0094: *
0095: *  LDU     (input) INTEGER
0096: *          The leading dimension of the array U.  LDU >= 1; if
0097: *          JOBZ = 'S' or 'A' or JOBZ = 'O' and M < N, LDU >= M.
0098: *
0099: *  VT      (output) COMPLEX array, dimension (LDVT,N)
0100: *          If JOBZ = 'A' or JOBZ = 'O' and M >= N, VT contains the
0101: *          N-by-N unitary matrix V**H;
0102: *          if JOBZ = 'S', VT contains the first min(M,N) rows of
0103: *          V**H (the right singular vectors, stored rowwise);
0104: *          if JOBZ = 'O' and M < N, or JOBZ = 'N', VT is not referenced.
0105: *
0106: *  LDVT    (input) INTEGER
0107: *          The leading dimension of the array VT.  LDVT >= 1; if
0108: *          JOBZ = 'A' or JOBZ = 'O' and M >= N, LDVT >= N;
0109: *          if JOBZ = 'S', LDVT >= min(M,N).
0110: *
0111: *  WORK    (workspace/output) COMPLEX array, dimension (MAX(1,LWORK))
0112: *          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
0113: *
0114: *  LWORK   (input) INTEGER
0115: *          The dimension of the array WORK. LWORK >= 1.
0116: *          if JOBZ = 'N', LWORK >= 2*min(M,N)+max(M,N).
0117: *          if JOBZ = 'O',
0118: *                LWORK >= 2*min(M,N)*min(M,N)+2*min(M,N)+max(M,N).
0119: *          if JOBZ = 'S' or 'A',
0120: *                LWORK >= min(M,N)*min(M,N)+2*min(M,N)+max(M,N).
0121: *          For good performance, LWORK should generally be larger.
0122: *
0123: *          If LWORK = -1, a workspace query is assumed.  The optimal
0124: *          size for the WORK array is calculated and stored in WORK(1),
0125: *          and no other work except argument checking is performed.
0126: *
0127: *  RWORK   (workspace) REAL array, dimension (MAX(1,LRWORK))
0128: *          If JOBZ = 'N', LRWORK >= 5*min(M,N).
0129: *          Otherwise, LRWORK >= 5*min(M,N)*min(M,N) + 7*min(M,N)
0130: *
0131: *  IWORK   (workspace) INTEGER array, dimension (8*min(M,N))
0132: *
0133: *  INFO    (output) INTEGER
0134: *          = 0:  successful exit.
0135: *          < 0:  if INFO = -i, the i-th argument had an illegal value.
0136: *          > 0:  The updating process of SBDSDC did not converge.
0137: *
0138: *  Further Details
0139: *  ===============
0140: *
0141: *  Based on contributions by
0142: *     Ming Gu and Huan Ren, Computer Science Division, University of
0143: *     California at Berkeley, USA
0144: *
0145: *  =====================================================================
0146: *
0147: *     .. Parameters ..
0148:       INTEGER            LQUERV
0149:       PARAMETER          ( LQUERV = -1 )
0150:       COMPLEX            CZERO, CONE
0151:       PARAMETER          ( CZERO = ( 0.0E+0, 0.0E+0 ),
0152:      $                   CONE = ( 1.0E+0, 0.0E+0 ) )
0153:       REAL               ZERO, ONE
0154:       PARAMETER          ( ZERO = 0.0E+0, ONE = 1.0E+0 )
0155: *     ..
0156: *     .. Local Scalars ..
0157:       LOGICAL            WNTQA, WNTQAS, WNTQN, WNTQO, WNTQS
0158:       INTEGER            BLK, CHUNK, I, IE, IERR, IL, IR, IRU, IRVT,
0159:      $                   ISCL, ITAU, ITAUP, ITAUQ, IU, IVT, LDWKVT,
0160:      $                   LDWRKL, LDWRKR, LDWRKU, MAXWRK, MINMN, MINWRK,
0161:      $                   MNTHR1, MNTHR2, NRWORK, NWORK, WRKBL
0162:       REAL               ANRM, BIGNUM, EPS, SMLNUM
0163: *     ..
0164: *     .. Local Arrays ..
0165:       INTEGER            IDUM( 1 )
0166:       REAL               DUM( 1 )
0167: *     ..
0168: *     .. External Subroutines ..
0169:       EXTERNAL           CGEBRD, CGELQF, CGEMM, CGEQRF, CLACP2, CLACPY,
0170:      $                   CLACRM, CLARCM, CLASCL, CLASET, CUNGBR, CUNGLQ,
0171:      $                   CUNGQR, CUNMBR, SBDSDC, SLASCL, XERBLA
0172: *     ..
0173: *     .. External Functions ..
0174:       LOGICAL            LSAME
0175:       INTEGER            ILAENV
0176:       REAL               CLANGE, SLAMCH
0177:       EXTERNAL           CLANGE, SLAMCH, ILAENV, LSAME
0178: *     ..
0179: *     .. Intrinsic Functions ..
0180:       INTRINSIC          INT, MAX, MIN, SQRT
0181: *     ..
0182: *     .. Executable Statements ..
0183: *
0184: *     Test the input arguments
0185: *
0186:       INFO = 0
0187:       MINMN = MIN( M, N )
0188:       MNTHR1 = INT( MINMN*17.0 / 9.0 )
0189:       MNTHR2 = INT( MINMN*5.0 / 3.0 )
0190:       WNTQA = LSAME( JOBZ, 'A' )
0191:       WNTQS = LSAME( JOBZ, 'S' )
0192:       WNTQAS = WNTQA .OR. WNTQS
0193:       WNTQO = LSAME( JOBZ, 'O' )
0194:       WNTQN = LSAME( JOBZ, 'N' )
0195:       MINWRK = 1
0196:       MAXWRK = 1
0197: *
0198:       IF( .NOT.( WNTQA .OR. WNTQS .OR. WNTQO .OR. WNTQN ) ) THEN
0199:          INFO = -1
0200:       ELSE IF( M.LT.0 ) THEN
0201:          INFO = -2
0202:       ELSE IF( N.LT.0 ) THEN
0203:          INFO = -3
0204:       ELSE IF( LDA.LT.MAX( 1, M ) ) THEN
0205:          INFO = -5
0206:       ELSE IF( LDU.LT.1 .OR. ( WNTQAS .AND. LDU.LT.M ) .OR.
0207:      $         ( WNTQO .AND. M.LT.N .AND. LDU.LT.M ) ) THEN
0208:          INFO = -8
0209:       ELSE IF( LDVT.LT.1 .OR. ( WNTQA .AND. LDVT.LT.N ) .OR.
0210:      $         ( WNTQS .AND. LDVT.LT.MINMN ) .OR.
0211:      $         ( WNTQO .AND. M.GE.N .AND. LDVT.LT.N ) ) THEN
0212:          INFO = -10
0213:       END IF
0214: *
0215: *     Compute workspace
0216: *      (Note: Comments in the code beginning "Workspace:" describe the
0217: *       minimal amount of workspace needed at that point in the code,
0218: *       as well as the preferred amount for good performance.
0219: *       CWorkspace refers to complex workspace, and RWorkspace to
0220: *       real workspace. NB refers to the optimal block size for the
0221: *       immediately following subroutine, as returned by ILAENV.)
0222: *
0223:       IF( INFO.EQ.0 .AND. M.GT.0 .AND. N.GT.0 ) THEN
0224:          IF( M.GE.N ) THEN
0225: *
0226: *           There is no complex work space needed for bidiagonal SVD
0227: *           The real work space needed for bidiagonal SVD is BDSPAC
0228: *           for computing singular values and singular vectors; BDSPAN
0229: *           for computing singular values only.
0230: *           BDSPAC = 5*N*N + 7*N
0231: *           BDSPAN = MAX(7*N+4, 3*N+2+SMLSIZ*(SMLSIZ+8))
0232: *
0233:             IF( M.GE.MNTHR1 ) THEN
0234:                IF( WNTQN ) THEN
0235: *
0236: *                 Path 1 (M much larger than N, JOBZ='N')
0237: *
0238:                   MAXWRK = N + N*ILAENV( 1, 'CGEQRF', ' ', M, N, -1,
0239:      $                     -1 )
0240:                   MAXWRK = MAX( MAXWRK, 2*N+2*N*
0241:      $                     ILAENV( 1, 'CGEBRD', ' ', N, N, -1, -1 ) )
0242:                   MINWRK = 3*N
0243:                ELSE IF( WNTQO ) THEN
0244: *
0245: *                 Path 2 (M much larger than N, JOBZ='O')
0246: *
0247:                   WRKBL = N + N*ILAENV( 1, 'CGEQRF', ' ', M, N, -1, -1 )
0248:                   WRKBL = MAX( WRKBL, N+N*ILAENV( 1, 'CUNGQR', ' ', M,
0249:      $                    N, N, -1 ) )
0250:                   WRKBL = MAX( WRKBL, 2*N+2*N*
0251:      $                    ILAENV( 1, 'CGEBRD', ' ', N, N, -1, -1 ) )
0252:                   WRKBL = MAX( WRKBL, 2*N+N*
0253:      $                    ILAENV( 1, 'CUNMBR', 'QLN', N, N, N, -1 ) )
0254:                   WRKBL = MAX( WRKBL, 2*N+N*
0255:      $                    ILAENV( 1, 'CUNMBR', 'PRC', N, N, N, -1 ) )
0256:                   MAXWRK = M*N + N*N + WRKBL
0257:                   MINWRK = 2*N*N + 3*N
0258:                ELSE IF( WNTQS ) THEN
0259: *
0260: *                 Path 3 (M much larger than N, JOBZ='S')
0261: *
0262:                   WRKBL = N + N*ILAENV( 1, 'CGEQRF', ' ', M, N, -1, -1 )
0263:                   WRKBL = MAX( WRKBL, N+N*ILAENV( 1, 'CUNGQR', ' ', M,
0264:      $                    N, N, -1 ) )
0265:                   WRKBL = MAX( WRKBL, 2*N+2*N*
0266:      $                    ILAENV( 1, 'CGEBRD', ' ', N, N, -1, -1 ) )
0267:                   WRKBL = MAX( WRKBL, 2*N+N*
0268:      $                    ILAENV( 1, 'CUNMBR', 'QLN', N, N, N, -1 ) )
0269:                   WRKBL = MAX( WRKBL, 2*N+N*
0270:      $                    ILAENV( 1, 'CUNMBR', 'PRC', N, N, N, -1 ) )
0271:                   MAXWRK = N*N + WRKBL
0272:                   MINWRK = N*N + 3*N
0273:                ELSE IF( WNTQA ) THEN
0274: *
0275: *                 Path 4 (M much larger than N, JOBZ='A')
0276: *
0277:                   WRKBL = N + N*ILAENV( 1, 'CGEQRF', ' ', M, N, -1, -1 )
0278:                   WRKBL = MAX( WRKBL, N+M*ILAENV( 1, 'CUNGQR', ' ', M,
0279:      $                    M, N, -1 ) )
0280:                   WRKBL = MAX( WRKBL, 2*N+2*N*
0281:      $                    ILAENV( 1, 'CGEBRD', ' ', N, N, -1, -1 ) )
0282:                   WRKBL = MAX( WRKBL, 2*N+N*
0283:      $                    ILAENV( 1, 'CUNMBR', 'QLN', N, N, N, -1 ) )
0284:                   WRKBL = MAX( WRKBL, 2*N+N*
0285:      $                    ILAENV( 1, 'CUNMBR', 'PRC', N, N, N, -1 ) )
0286:                   MAXWRK = N*N + WRKBL
0287:                   MINWRK = N*N + 2*N + M
0288:                END IF
0289:             ELSE IF( M.GE.MNTHR2 ) THEN
0290: *
0291: *              Path 5 (M much larger than N, but not as much as MNTHR1)
0292: *
0293:                MAXWRK = 2*N + ( M+N )*ILAENV( 1, 'CGEBRD', ' ', M, N,
0294:      $                  -1, -1 )
0295:                MINWRK = 2*N + M
0296:                IF( WNTQO ) THEN
0297:                   MAXWRK = MAX( MAXWRK, 2*N+N*
0298:      $                     ILAENV( 1, 'CUNGBR', 'P', N, N, N, -1 ) )
0299:                   MAXWRK = MAX( MAXWRK, 2*N+N*
0300:      $                     ILAENV( 1, 'CUNGBR', 'Q', M, N, N, -1 ) )
0301:                   MAXWRK = MAXWRK + M*N
0302:                   MINWRK = MINWRK + N*N
0303:                ELSE IF( WNTQS ) THEN
0304:                   MAXWRK = MAX( MAXWRK, 2*N+N*
0305:      $                     ILAENV( 1, 'CUNGBR', 'P', N, N, N, -1 ) )
0306:                   MAXWRK = MAX( MAXWRK, 2*N+N*
0307:      $                     ILAENV( 1, 'CUNGBR', 'Q', M, N, N, -1 ) )
0308:                ELSE IF( WNTQA ) THEN
0309:                   MAXWRK = MAX( MAXWRK, 2*N+N*
0310:      $                     ILAENV( 1, 'CUNGBR', 'P', N, N, N, -1 ) )
0311:                   MAXWRK = MAX( MAXWRK, 2*N+M*
0312:      $                     ILAENV( 1, 'CUNGBR', 'Q', M, M, N, -1 ) )
0313:                END IF
0314:             ELSE
0315: *
0316: *              Path 6 (M at least N, but not much larger)
0317: *
0318:                MAXWRK = 2*N + ( M+N )*ILAENV( 1, 'CGEBRD', ' ', M, N,
0319:      $                  -1, -1 )
0320:                MINWRK = 2*N + M
0321:                IF( WNTQO ) THEN
0322:                   MAXWRK = MAX( MAXWRK, 2*N+N*
0323:      $                     ILAENV( 1, 'CUNMBR', 'PRC', N, N, N, -1 ) )
0324:                   MAXWRK = MAX( MAXWRK, 2*N+N*
0325:      $                     ILAENV( 1, 'CUNMBR', 'QLN', M, N, N, -1 ) )
0326:                   MAXWRK = MAXWRK + M*N
0327:                   MINWRK = MINWRK + N*N
0328:                ELSE IF( WNTQS ) THEN
0329:                   MAXWRK = MAX( MAXWRK, 2*N+N*
0330:      $                     ILAENV( 1, 'CUNMBR', 'PRC', N, N, N, -1 ) )
0331:                   MAXWRK = MAX( MAXWRK, 2*N+N*
0332:      $                     ILAENV( 1, 'CUNMBR', 'QLN', M, N, N, -1 ) )
0333:                ELSE IF( WNTQA ) THEN
0334:                   MAXWRK = MAX( MAXWRK, 2*N+N*
0335:      $                     ILAENV( 1, 'CUNGBR', 'PRC', N, N, N, -1 ) )
0336:                   MAXWRK = MAX( MAXWRK, 2*N+M*
0337:      $                     ILAENV( 1, 'CUNGBR', 'QLN', M, M, N, -1 ) )
0338:                END IF
0339:             END IF
0340:          ELSE
0341: *
0342: *           There is no complex work space needed for bidiagonal SVD
0343: *           The real work space needed for bidiagonal SVD is BDSPAC
0344: *           for computing singular values and singular vectors; BDSPAN
0345: *           for computing singular values only.
0346: *           BDSPAC = 5*M*M + 7*M
0347: *           BDSPAN = MAX(7*M+4, 3*M+2+SMLSIZ*(SMLSIZ+8))
0348: *
0349:             IF( N.GE.MNTHR1 ) THEN
0350:                IF( WNTQN ) THEN
0351: *
0352: *                 Path 1t (N much larger than M, JOBZ='N')
0353: *
0354:                   MAXWRK = M + M*ILAENV( 1, 'CGELQF', ' ', M, N, -1,
0355:      $                     -1 )
0356:                   MAXWRK = MAX( MAXWRK, 2*M+2*M*
0357:      $                     ILAENV( 1, 'CGEBRD', ' ', M, M, -1, -1 ) )
0358:                   MINWRK = 3*M
0359:                ELSE IF( WNTQO ) THEN
0360: *
0361: *                 Path 2t (N much larger than M, JOBZ='O')
0362: *
0363:                   WRKBL = M + M*ILAENV( 1, 'CGELQF', ' ', M, N, -1, -1 )
0364:                   WRKBL = MAX( WRKBL, M+M*ILAENV( 1, 'CUNGLQ', ' ', M,
0365:      $                    N, M, -1 ) )
0366:                   WRKBL = MAX( WRKBL, 2*M+2*M*
0367:      $                    ILAENV( 1, 'CGEBRD', ' ', M, M, -1, -1 ) )
0368:                   WRKBL = MAX( WRKBL, 2*M+M*
0369:      $                    ILAENV( 1, 'CUNMBR', 'PRC', M, M, M, -1 ) )
0370:                   WRKBL = MAX( WRKBL, 2*M+M*
0371:      $                    ILAENV( 1, 'CUNMBR', 'QLN', M, M, M, -1 ) )
0372:                   MAXWRK = M*N + M*M + WRKBL
0373:                   MINWRK = 2*M*M + 3*M
0374:                ELSE IF( WNTQS ) THEN
0375: *
0376: *                 Path 3t (N much larger than M, JOBZ='S')
0377: *
0378:                   WRKBL = M + M*ILAENV( 1, 'CGELQF', ' ', M, N, -1, -1 )
0379:                   WRKBL = MAX( WRKBL, M+M*ILAENV( 1, 'CUNGLQ', ' ', M,
0380:      $                    N, M, -1 ) )
0381:                   WRKBL = MAX( WRKBL, 2*M+2*M*
0382:      $                    ILAENV( 1, 'CGEBRD', ' ', M, M, -1, -1 ) )
0383:                   WRKBL = MAX( WRKBL, 2*M+M*
0384:      $                    ILAENV( 1, 'CUNMBR', 'PRC', M, M, M, -1 ) )
0385:                   WRKBL = MAX( WRKBL, 2*M+M*
0386:      $                    ILAENV( 1, 'CUNMBR', 'QLN', M, M, M, -1 ) )
0387:                   MAXWRK = M*M + WRKBL
0388:                   MINWRK = M*M + 3*M
0389:                ELSE IF( WNTQA ) THEN
0390: *
0391: *                 Path 4t (N much larger than M, JOBZ='A')
0392: *
0393:                   WRKBL = M + M*ILAENV( 1, 'CGELQF', ' ', M, N, -1, -1 )
0394:                   WRKBL = MAX( WRKBL, M+N*ILAENV( 1, 'CUNGLQ', ' ', N,
0395:      $                    N, M, -1 ) )
0396:                   WRKBL = MAX( WRKBL, 2*M+2*M*
0397:      $                    ILAENV( 1, 'CGEBRD', ' ', M, M, -1, -1 ) )
0398:                   WRKBL = MAX( WRKBL, 2*M+M*
0399:      $                    ILAENV( 1, 'CUNMBR', 'PRC', M, M, M, -1 ) )
0400:                   WRKBL = MAX( WRKBL, 2*M+M*
0401:      $                    ILAENV( 1, 'CUNMBR', 'QLN', M, M, M, -1 ) )
0402:                   MAXWRK = M*M + WRKBL
0403:                   MINWRK = M*M + 2*M + N
0404:                END IF
0405:             ELSE IF( N.GE.MNTHR2 ) THEN
0406: *
0407: *              Path 5t (N much larger than M, but not as much as MNTHR1)
0408: *
0409:                MAXWRK = 2*M + ( M+N )*ILAENV( 1, 'CGEBRD', ' ', M, N,
0410:      $                  -1, -1 )
0411:                MINWRK = 2*M + N
0412:                IF( WNTQO ) THEN
0413:                   MAXWRK = MAX( MAXWRK, 2*M+M*
0414:      $                     ILAENV( 1, 'CUNGBR', 'P', M, N, M, -1 ) )
0415:                   MAXWRK = MAX( MAXWRK, 2*M+M*
0416:      $                     ILAENV( 1, 'CUNGBR', 'Q', M, M, N, -1 ) )
0417:                   MAXWRK = MAXWRK + M*N
0418:                   MINWRK = MINWRK + M*M
0419:                ELSE IF( WNTQS ) THEN
0420:                   MAXWRK = MAX( MAXWRK, 2*M+M*
0421:      $                     ILAENV( 1, 'CUNGBR', 'P', M, N, M, -1 ) )
0422:                   MAXWRK = MAX( MAXWRK, 2*M+M*
0423:      $                     ILAENV( 1, 'CUNGBR', 'Q', M, M, N, -1 ) )
0424:                ELSE IF( WNTQA ) THEN
0425:                   MAXWRK = MAX( MAXWRK, 2*M+N*
0426:      $                     ILAENV( 1, 'CUNGBR', 'P', N, N, M, -1 ) )
0427:                   MAXWRK = MAX( MAXWRK, 2*M+M*
0428:      $                     ILAENV( 1, 'CUNGBR', 'Q', M, M, N, -1 ) )
0429:                END IF
0430:             ELSE
0431: *
0432: *              Path 6t (N greater than M, but not much larger)
0433: *
0434:                MAXWRK = 2*M + ( M+N )*ILAENV( 1, 'CGEBRD', ' ', M, N,
0435:      $                  -1, -1 )
0436:                MINWRK = 2*M + N
0437:                IF( WNTQO ) THEN
0438:                   MAXWRK = MAX( MAXWRK, 2*M+M*
0439:      $                     ILAENV( 1, 'CUNMBR', 'PRC', M, N, M, -1 ) )
0440:                   MAXWRK = MAX( MAXWRK, 2*M+M*
0441:      $                     ILAENV( 1, 'CUNMBR', 'QLN', M, M, N, -1 ) )
0442:                   MAXWRK = MAXWRK + M*N
0443:                   MINWRK = MINWRK + M*M
0444:                ELSE IF( WNTQS ) THEN
0445:                   MAXWRK = MAX( MAXWRK, 2*M+M*
0446:      $                     ILAENV( 1, 'CUNGBR', 'PRC', M, N, M, -1 ) )
0447:                   MAXWRK = MAX( MAXWRK, 2*M+M*
0448:      $                     ILAENV( 1, 'CUNGBR', 'QLN', M, M, N, -1 ) )
0449:                ELSE IF( WNTQA ) THEN
0450:                   MAXWRK = MAX( MAXWRK, 2*M+N*
0451:      $                     ILAENV( 1, 'CUNGBR', 'PRC', N, N, M, -1 ) )
0452:                   MAXWRK = MAX( MAXWRK, 2*M+M*
0453:      $                     ILAENV( 1, 'CUNGBR', 'QLN', M, M, N, -1 ) )
0454:                END IF
0455:             END IF
0456:          END IF
0457:          MAXWRK = MAX( MAXWRK, MINWRK )
0458:       END IF
0459:       IF( INFO.EQ.0 ) THEN
0460:          WORK( 1 ) = MAXWRK
0461:          IF( LWORK.LT.MINWRK .AND. LWORK.NE.LQUERV )
0462:      $      INFO = -13
0463:       END IF
0464: *
0465: *     Quick returns
0466: *
0467:       IF( INFO.NE.0 ) THEN
0468:          CALL XERBLA( 'CGESDD', -INFO )
0469:          RETURN
0470:       END IF
0471:       IF( LWORK.EQ.LQUERV )
0472:      $   RETURN
0473:       IF( M.EQ.0 .OR. N.EQ.0 ) THEN
0474:          RETURN
0475:       END IF
0476: *
0477: *     Get machine constants
0478: *
0479:       EPS = SLAMCH( 'P' )
0480:       SMLNUM = SQRT( SLAMCH( 'S' ) ) / EPS
0481:       BIGNUM = ONE / SMLNUM
0482: *
0483: *     Scale A if max element outside range [SMLNUM,BIGNUM]
0484: *
0485:       ANRM = CLANGE( 'M', M, N, A, LDA, DUM )
0486:       ISCL = 0
0487:       IF( ANRM.GT.ZERO .AND. ANRM.LT.SMLNUM ) THEN
0488:          ISCL = 1
0489:          CALL CLASCL( 'G', 0, 0, ANRM, SMLNUM, M, N, A, LDA, IERR )
0490:       ELSE IF( ANRM.GT.BIGNUM ) THEN
0491:          ISCL = 1
0492:          CALL CLASCL( 'G', 0, 0, ANRM, BIGNUM, M, N, A, LDA, IERR )
0493:       END IF
0494: *
0495:       IF( M.GE.N ) THEN
0496: *
0497: *        A has at least as many rows as columns. If A has sufficiently
0498: *        more rows than columns, first reduce using the QR
0499: *        decomposition (if sufficient workspace available)
0500: *
0501:          IF( M.GE.MNTHR1 ) THEN
0502: *
0503:             IF( WNTQN ) THEN
0504: *
0505: *              Path 1 (M much larger than N, JOBZ='N')
0506: *              No singular vectors to be computed
0507: *
0508:                ITAU = 1
0509:                NWORK = ITAU + N
0510: *
0511: *              Compute A=Q*R
0512: *              (CWorkspace: need 2*N, prefer N+N*NB)
0513: *              (RWorkspace: need 0)
0514: *
0515:                CALL CGEQRF( M, N, A, LDA, WORK( ITAU ), WORK( NWORK ),
0516:      $                      LWORK-NWORK+1, IERR )
0517: *
0518: *              Zero out below R
0519: *
0520:                CALL CLASET( 'L', N-1, N-1, CZERO, CZERO, A( 2, 1 ),
0521:      $                      LDA )
0522:                IE = 1
0523:                ITAUQ = 1
0524:                ITAUP = ITAUQ + N
0525:                NWORK = ITAUP + N
0526: *
0527: *              Bidiagonalize R in A
0528: *              (CWorkspace: need 3*N, prefer 2*N+2*N*NB)
0529: *              (RWorkspace: need N)
0530: *
0531:                CALL CGEBRD( N, N, A, LDA, S, RWORK( IE ), WORK( ITAUQ ),
0532:      $                      WORK( ITAUP ), WORK( NWORK ), LWORK-NWORK+1,
0533:      $                      IERR )
0534:                NRWORK = IE + N
0535: *
0536: *              Perform bidiagonal SVD, compute singular values only
0537: *              (CWorkspace: 0)
0538: *              (RWorkspace: need BDSPAN)
0539: *
0540:                CALL SBDSDC( 'U', 'N', N, S, RWORK( IE ), DUM, 1, DUM, 1,
0541:      $                      DUM, IDUM, RWORK( NRWORK ), IWORK, INFO )
0542: *
0543:             ELSE IF( WNTQO ) THEN
0544: *
0545: *              Path 2 (M much larger than N, JOBZ='O')
0546: *              N left singular vectors to be overwritten on A and
0547: *              N right singular vectors to be computed in VT
0548: *
0549:                IU = 1
0550: *
0551: *              WORK(IU) is N by N
0552: *
0553:                LDWRKU = N
0554:                IR = IU + LDWRKU*N
0555:                IF( LWORK.GE.M*N+N*N+3*N ) THEN
0556: *
0557: *                 WORK(IR) is M by N
0558: *
0559:                   LDWRKR = M
0560:                ELSE
0561:                   LDWRKR = ( LWORK-N*N-3*N ) / N
0562:                END IF
0563:                ITAU = IR + LDWRKR*N
0564:                NWORK = ITAU + N
0565: *
0566: *              Compute A=Q*R
0567: *              (CWorkspace: need N*N+2*N, prefer M*N+N+N*NB)
0568: *              (RWorkspace: 0)
0569: *
0570:                CALL CGEQRF( M, N, A, LDA, WORK( ITAU ), WORK( NWORK ),
0571:      $                      LWORK-NWORK+1, IERR )
0572: *
0573: *              Copy R to WORK( IR ), zeroing out below it
0574: *
0575:                CALL CLACPY( 'U', N, N, A, LDA, WORK( IR ), LDWRKR )
0576:                CALL CLASET( 'L', N-1, N-1, CZERO, CZERO, WORK( IR+1 ),
0577:      $                      LDWRKR )
0578: *
0579: *              Generate Q in A
0580: *              (CWorkspace: need 2*N, prefer N+N*NB)
0581: *              (RWorkspace: 0)
0582: *
0583:                CALL CUNGQR( M, N, N, A, LDA, WORK( ITAU ),
0584:      $                      WORK( NWORK ), LWORK-NWORK+1, IERR )
0585:                IE = 1
0586:                ITAUQ = ITAU
0587:                ITAUP = ITAUQ + N
0588:                NWORK = ITAUP + N
0589: *
0590: *              Bidiagonalize R in WORK(IR)
0591: *              (CWorkspace: need N*N+3*N, prefer M*N+2*N+2*N*NB)
0592: *              (RWorkspace: need N)
0593: *
0594:                CALL CGEBRD( N, N, WORK( IR ), LDWRKR, S, RWORK( IE ),
0595:      $                      WORK( ITAUQ ), WORK( ITAUP ), WORK( NWORK ),
0596:      $                      LWORK-NWORK+1, IERR )
0597: *
0598: *              Perform bidiagonal SVD, computing left singular vectors
0599: *              of R in WORK(IRU) and computing right singular vectors
0600: *              of R in WORK(IRVT)
0601: *              (CWorkspace: need 0)
0602: *              (RWorkspace: need BDSPAC)
0603: *
0604:                IRU = IE + N
0605:                IRVT = IRU + N*N
0606:                NRWORK = IRVT + N*N
0607:                CALL SBDSDC( 'U', 'I', N, S, RWORK( IE ), RWORK( IRU ),
0608:      $                      N, RWORK( IRVT ), N, DUM, IDUM,
0609:      $                      RWORK( NRWORK ), IWORK, INFO )
0610: *
0611: *              Copy real matrix RWORK(IRU) to complex matrix WORK(IU)
0612: *              Overwrite WORK(IU) by the left singular vectors of R
0613: *              (CWorkspace: need 2*N*N+3*N, prefer M*N+N*N+2*N+N*NB)
0614: *              (RWorkspace: 0)
0615: *
0616:                CALL CLACP2( 'F', N, N, RWORK( IRU ), N, WORK( IU ),
0617:      $                      LDWRKU )
0618:                CALL CUNMBR( 'Q', 'L', 'N', N, N, N, WORK( IR ), LDWRKR,
0619:      $                      WORK( ITAUQ ), WORK( IU ), LDWRKU,
0620:      $                      WORK( NWORK ), LWORK-NWORK+1, IERR )
0621: *
0622: *              Copy real matrix RWORK(IRVT) to complex matrix VT
0623: *              Overwrite VT by the right singular vectors of R
0624: *              (CWorkspace: need N*N+3*N, prefer M*N+2*N+N*NB)
0625: *              (RWorkspace: 0)
0626: *
0627:                CALL CLACP2( 'F', N, N, RWORK( IRVT ), N, VT, LDVT )
0628:                CALL CUNMBR( 'P', 'R', 'C', N, N, N, WORK( IR ), LDWRKR,
0629:      $                      WORK( ITAUP ), VT, LDVT, WORK( NWORK ),
0630:      $                      LWORK-NWORK+1, IERR )
0631: *
0632: *              Multiply Q in A by left singular vectors of R in
0633: *              WORK(IU), storing result in WORK(IR) and copying to A
0634: *              (CWorkspace: need 2*N*N, prefer N*N+M*N)
0635: *              (RWorkspace: 0)
0636: *
0637:                DO 10 I = 1, M, LDWRKR
0638:                   CHUNK = MIN( M-I+1, LDWRKR )
0639:                   CALL CGEMM( 'N', 'N', CHUNK, N, N, CONE, A( I, 1 ),
0640:      $                        LDA, WORK( IU ), LDWRKU, CZERO,
0641:      $                        WORK( IR ), LDWRKR )
0642:                   CALL CLACPY( 'F', CHUNK, N, WORK( IR ), LDWRKR,
0643:      $                         A( I, 1 ), LDA )
0644:    10          CONTINUE
0645: *
0646:             ELSE IF( WNTQS ) THEN
0647: *
0648: *              Path 3 (M much larger than N, JOBZ='S')
0649: *              N left singular vectors to be computed in U and
0650: *              N right singular vectors to be computed in VT
0651: *
0652:                IR = 1
0653: *
0654: *              WORK(IR) is N by N
0655: *
0656:                LDWRKR = N
0657:                ITAU = IR + LDWRKR*N
0658:                NWORK = ITAU + N
0659: *
0660: *              Compute A=Q*R
0661: *              (CWorkspace: need N*N+2*N, prefer N*N+N+N*NB)
0662: *              (RWorkspace: 0)
0663: *
0664:                CALL CGEQRF( M, N, A, LDA, WORK( ITAU ), WORK( NWORK ),
0665:      $                      LWORK-NWORK+1, IERR )
0666: *
0667: *              Copy R to WORK(IR), zeroing out below it
0668: *
0669:                CALL CLACPY( 'U', N, N, A, LDA, WORK( IR ), LDWRKR )
0670:                CALL CLASET( 'L', N-1, N-1, CZERO, CZERO, WORK( IR+1 ),
0671:      $                      LDWRKR )
0672: *
0673: *              Generate Q in A
0674: *              (CWorkspace: need 2*N, prefer N+N*NB)
0675: *              (RWorkspace: 0)
0676: *
0677:                CALL CUNGQR( M, N, N, A, LDA, WORK( ITAU ),
0678:      $                      WORK( NWORK ), LWORK-NWORK+1, IERR )
0679:                IE = 1
0680:                ITAUQ = ITAU
0681:                ITAUP = ITAUQ + N
0682:                NWORK = ITAUP + N
0683: *
0684: *              Bidiagonalize R in WORK(IR)
0685: *              (CWorkspace: need N*N+3*N, prefer N*N+2*N+2*N*NB)
0686: *              (RWorkspace: need N)
0687: *
0688:                CALL CGEBRD( N, N, WORK( IR ), LDWRKR, S, RWORK( IE ),
0689:      $                      WORK( ITAUQ ), WORK( ITAUP ), WORK( NWORK ),
0690:      $                      LWORK-NWORK+1, IERR )
0691: *
0692: *              Perform bidiagonal SVD, computing left singular vectors
0693: *              of bidiagonal matrix in RWORK(IRU) and computing right
0694: *              singular vectors of bidiagonal matrix in RWORK(IRVT)
0695: *              (CWorkspace: need 0)
0696: *              (RWorkspace: need BDSPAC)
0697: *
0698:                IRU = IE + N
0699:                IRVT = IRU + N*N
0700:                NRWORK = IRVT + N*N
0701:                CALL SBDSDC( 'U', 'I', N, S, RWORK( IE ), RWORK( IRU ),
0702:      $                      N, RWORK( IRVT ), N, DUM, IDUM,
0703:      $                      RWORK( NRWORK ), IWORK, INFO )
0704: *
0705: *              Copy real matrix RWORK(IRU) to complex matrix U
0706: *              Overwrite U by left singular vectors of R
0707: *              (CWorkspace: need N*N+3*N, prefer N*N+2*N+N*NB)
0708: *              (RWorkspace: 0)
0709: *
0710:                CALL CLACP2( 'F', N, N, RWORK( IRU ), N, U, LDU )
0711:                CALL CUNMBR( 'Q', 'L', 'N', N, N, N, WORK( IR ), LDWRKR,
0712:      $                      WORK( ITAUQ ), U, LDU, WORK( NWORK ),
0713:      $                      LWORK-NWORK+1, IERR )
0714: *
0715: *              Copy real matrix RWORK(IRVT) to complex matrix VT
0716: *              Overwrite VT by right singular vectors of R
0717: *              (CWorkspace: need N*N+3*N, prefer N*N+2*N+N*NB)
0718: *              (RWorkspace: 0)
0719: *
0720:                CALL CLACP2( 'F', N, N, RWORK( IRVT ), N, VT, LDVT )
0721:                CALL CUNMBR( 'P', 'R', 'C', N, N, N, WORK( IR ), LDWRKR,
0722:      $                      WORK( ITAUP ), VT, LDVT, WORK( NWORK ),
0723:      $                      LWORK-NWORK+1, IERR )
0724: *
0725: *              Multiply Q in A by left singular vectors of R in
0726: *              WORK(IR), storing result in U
0727: *              (CWorkspace: need N*N)
0728: *              (RWorkspace: 0)
0729: *
0730:                CALL CLACPY( 'F', N, N, U, LDU, WORK( IR ), LDWRKR )
0731:                CALL CGEMM( 'N', 'N', M, N, N, CONE, A, LDA, WORK( IR ),
0732:      $                     LDWRKR, CZERO, U, LDU )
0733: *
0734:             ELSE IF( WNTQA ) THEN
0735: *
0736: *              Path 4 (M much larger than N, JOBZ='A')
0737: *              M left singular vectors to be computed in U and
0738: *              N right singular vectors to be computed in VT
0739: *
0740:                IU = 1
0741: *
0742: *              WORK(IU) is N by N
0743: *
0744:                LDWRKU = N
0745:                ITAU = IU + LDWRKU*N
0746:                NWORK = ITAU + N
0747: *
0748: *              Compute A=Q*R, copying result to U
0749: *              (CWorkspace: need 2*N, prefer N+N*NB)
0750: *              (RWorkspace: 0)
0751: *
0752:                CALL CGEQRF( M, N, A, LDA, WORK( ITAU ), WORK( NWORK ),
0753:      $                      LWORK-NWORK+1, IERR )
0754:                CALL CLACPY( 'L', M, N, A, LDA, U, LDU )
0755: *
0756: *              Generate Q in U
0757: *              (CWorkspace: need N+M, prefer N+M*NB)
0758: *              (RWorkspace: 0)
0759: *
0760:                CALL CUNGQR( M, M, N, U, LDU, WORK( ITAU ),
0761:      $                      WORK( NWORK ), LWORK-NWORK+1, IERR )
0762: *
0763: *              Produce R in A, zeroing out below it
0764: *
0765:                CALL CLASET( 'L', N-1, N-1, CZERO, CZERO, A( 2, 1 ),
0766:      $                      LDA )
0767:                IE = 1
0768:                ITAUQ = ITAU
0769:                ITAUP = ITAUQ + N
0770:                NWORK = ITAUP + N
0771: *
0772: *              Bidiagonalize R in A
0773: *              (CWorkspace: need 3*N, prefer 2*N+2*N*NB)
0774: *              (RWorkspace: need N)
0775: *
0776:                CALL CGEBRD( N, N, A, LDA, S, RWORK( IE ), WORK( ITAUQ ),
0777:      $                      WORK( ITAUP ), WORK( NWORK ), LWORK-NWORK+1,
0778:      $                      IERR )
0779:                IRU = IE + N
0780:                IRVT = IRU + N*N
0781:                NRWORK = IRVT + N*N
0782: *
0783: *              Perform bidiagonal SVD, computing left singular vectors
0784: *              of bidiagonal matrix in RWORK(IRU) and computing right
0785: *              singular vectors of bidiagonal matrix in RWORK(IRVT)
0786: *              (CWorkspace: need 0)
0787: *              (RWorkspace: need BDSPAC)
0788: *
0789:                CALL SBDSDC( 'U', 'I', N, S, RWORK( IE ), RWORK( IRU ),
0790:      $                      N, RWORK( IRVT ), N, DUM, IDUM,
0791:      $                      RWORK( NRWORK ), IWORK, INFO )
0792: *
0793: *              Copy real matrix RWORK(IRU) to complex matrix WORK(IU)
0794: *              Overwrite WORK(IU) by left singular vectors of R
0795: *              (CWorkspace: need N*N+3*N, prefer N*N+2*N+N*NB)
0796: *              (RWorkspace: 0)
0797: *
0798:                CALL CLACP2( 'F', N, N, RWORK( IRU ), N, WORK( IU ),
0799:      $                      LDWRKU )
0800:                CALL CUNMBR( 'Q', 'L', 'N', N, N, N, A, LDA,
0801:      $                      WORK( ITAUQ ), WORK( IU ), LDWRKU,
0802:      $                      WORK( NWORK ), LWORK-NWORK+1, IERR )
0803: *
0804: *              Copy real matrix RWORK(IRVT) to complex matrix VT
0805: *              Overwrite VT by right singular vectors of R
0806: *              (CWorkspace: need 3*N, prefer 2*N+N*NB)
0807: *              (RWorkspace: 0)
0808: *
0809:                CALL CLACP2( 'F', N, N, RWORK( IRVT ), N, VT, LDVT )
0810:                CALL CUNMBR( 'P', 'R', 'C', N, N, N, A, LDA,
0811:      $                      WORK( ITAUP ), VT, LDVT, WORK( NWORK ),
0812:      $                      LWORK-NWORK+1, IERR )
0813: *
0814: *              Multiply Q in U by left singular vectors of R in
0815: *              WORK(IU), storing result in A
0816: *              (CWorkspace: need N*N)
0817: *              (RWorkspace: 0)
0818: *
0819:                CALL CGEMM( 'N', 'N', M, N, N, CONE, U, LDU, WORK( IU ),
0820:      $                     LDWRKU, CZERO, A, LDA )
0821: *
0822: *              Copy left singular vectors of A from A to U
0823: *
0824:                CALL CLACPY( 'F', M, N, A, LDA, U, LDU )
0825: *
0826:             END IF
0827: *
0828:          ELSE IF( M.GE.MNTHR2 ) THEN
0829: *
0830: *           MNTHR2 <= M < MNTHR1
0831: *
0832: *           Path 5 (M much larger than N, but not as much as MNTHR1)
0833: *           Reduce to bidiagonal form without QR decomposition, use
0834: *           CUNGBR and matrix multiplication to compute singular vectors
0835: *
0836:             IE = 1
0837:             NRWORK = IE + N
0838:             ITAUQ = 1
0839:             ITAUP = ITAUQ + N
0840:             NWORK = ITAUP + N
0841: *
0842: *           Bidiagonalize A
0843: *           (CWorkspace: need 2*N+M, prefer 2*N+(M+N)*NB)
0844: *           (RWorkspace: need N)
0845: *
0846:             CALL CGEBRD( M, N, A, LDA, S, RWORK( IE ), WORK( ITAUQ ),
0847:      $                   WORK( ITAUP ), WORK( NWORK ), LWORK-NWORK+1,
0848:      $                   IERR )
0849:             IF( WNTQN ) THEN
0850: *
0851: *              Compute singular values only
0852: *              (Cworkspace: 0)
0853: *              (Rworkspace: need BDSPAN)
0854: *
0855:                CALL SBDSDC( 'U', 'N', N, S, RWORK( IE ), DUM, 1, DUM, 1,
0856:      $                      DUM, IDUM, RWORK( NRWORK ), IWORK, INFO )
0857:             ELSE IF( WNTQO ) THEN
0858:                IU = NWORK
0859:                IRU = NRWORK
0860:                IRVT = IRU + N*N
0861:                NRWORK = IRVT + N*N
0862: *
0863: *              Copy A to VT, generate P**H
0864: *              (Cworkspace: need 2*N, prefer N+N*NB)
0865: *              (Rworkspace: 0)
0866: *
0867:                CALL CLACPY( 'U', N, N, A, LDA, VT, LDVT )
0868:                CALL CUNGBR( 'P', N, N, N, VT, LDVT, WORK( ITAUP ),
0869:      $                      WORK( NWORK ), LWORK-NWORK+1, IERR )
0870: *
0871: *              Generate Q in A
0872: *              (CWorkspace: need 2*N, prefer N+N*NB)
0873: *              (RWorkspace: 0)
0874: *
0875:                CALL CUNGBR( 'Q', M, N, N, A, LDA, WORK( ITAUQ ),
0876:      $                      WORK( NWORK ), LWORK-NWORK+1, IERR )
0877: *
0878:                IF( LWORK.GE.M*N+3*N ) THEN
0879: *
0880: *                 WORK( IU ) is M by N
0881: *
0882:                   LDWRKU = M
0883:                ELSE
0884: *
0885: *                 WORK(IU) is LDWRKU by N
0886: *
0887:                   LDWRKU = ( LWORK-3*N ) / N
0888:                END IF
0889:                NWORK = IU + LDWRKU*N
0890: *
0891: *              Perform bidiagonal SVD, computing left singular vectors
0892: *              of bidiagonal matrix in RWORK(IRU) and computing right
0893: *              singular vectors of bidiagonal matrix in RWORK(IRVT)
0894: *              (CWorkspace: need 0)
0895: *              (RWorkspace: need BDSPAC)
0896: *
0897:                CALL SBDSDC( 'U', 'I', N, S, RWORK( IE ), RWORK( IRU ),
0898:      $                      N, RWORK( IRVT ), N, DUM, IDUM,
0899:      $                      RWORK( NRWORK ), IWORK, INFO )
0900: *
0901: *              Multiply real matrix RWORK(IRVT) by P**H in VT,
0902: *              storing the result in WORK(IU), copying to VT
0903: *              (Cworkspace: need 0)
0904: *              (Rworkspace: need 3*N*N)
0905: *
0906:                CALL CLARCM( N, N, RWORK( IRVT ), N, VT, LDVT,
0907:      $                      WORK( IU ), LDWRKU, RWORK( NRWORK ) )
0908:                CALL CLACPY( 'F', N, N, WORK( IU ), LDWRKU, VT, LDVT )
0909: *
0910: *              Multiply Q in A by real matrix RWORK(IRU), storing the
0911: *              result in WORK(IU), copying to A
0912: *              (CWorkspace: need N*N, prefer M*N)
0913: *              (Rworkspace: need 3*N*N, prefer N*N+2*M*N)
0914: *
0915:                NRWORK = IRVT
0916:                DO 20 I = 1, M, LDWRKU
0917:                   CHUNK = MIN( M-I+1, LDWRKU )
0918:                   CALL CLACRM( CHUNK, N, A( I, 1 ), LDA, RWORK( IRU ),
0919:      $                         N, WORK( IU ), LDWRKU, RWORK( NRWORK ) )
0920:                   CALL CLACPY( 'F', CHUNK, N, WORK( IU ), LDWRKU,
0921:      $                         A( I, 1 ), LDA )
0922:    20          CONTINUE
0923: *
0924:             ELSE IF( WNTQS ) THEN
0925: *
0926: *              Copy A to VT, generate P**H
0927: *              (Cworkspace: need 2*N, prefer N+N*NB)
0928: *              (Rworkspace: 0)
0929: *
0930:                CALL CLACPY( 'U', N, N, A, LDA, VT, LDVT )
0931:                CALL CUNGBR( 'P', N, N, N, VT, LDVT, WORK( ITAUP ),
0932:      $                      WORK( NWORK ), LWORK-NWORK+1, IERR )
0933: *
0934: *              Copy A to U, generate Q
0935: *              (Cworkspace: need 2*N, prefer N+N*NB)
0936: *              (Rworkspace: 0)
0937: *
0938:                CALL CLACPY( 'L', M, N, A, LDA, U, LDU )
0939:                CALL CUNGBR( 'Q', M, N, N, U, LDU, WORK( ITAUQ ),
0940:      $                      WORK( NWORK ), LWORK-NWORK+1, IERR )
0941: *
0942: *              Perform bidiagonal SVD, computing left singular vectors
0943: *              of bidiagonal matrix in RWORK(IRU) and computing right
0944: *              singular vectors of bidiagonal matrix in RWORK(IRVT)
0945: *              (CWorkspace: need 0)
0946: *              (RWorkspace: need BDSPAC)
0947: *
0948:                IRU = NRWORK
0949:                IRVT = IRU + N*N
0950:                NRWORK = IRVT + N*N
0951:                CALL SBDSDC( 'U', 'I', N, S, RWORK( IE ), RWORK( IRU ),
0952:      $                      N, RWORK( IRVT ), N, DUM, IDUM,
0953:      $                      RWORK( NRWORK ), IWORK, INFO )
0954: *
0955: *              Multiply real matrix RWORK(IRVT) by P**H in VT,
0956: *              storing the result in A, copying to VT
0957: *              (Cworkspace: need 0)
0958: *              (Rworkspace: need 3*N*N)
0959: *
0960:                CALL CLARCM( N, N, RWORK( IRVT ), N, VT, LDVT, A, LDA,
0961:      $                      RWORK( NRWORK ) )
0962:                CALL CLACPY( 'F', N, N, A, LDA, VT, LDVT )
0963: *
0964: *              Multiply Q in U by real matrix RWORK(IRU), storing the
0965: *              result in A, copying to U
0966: *              (CWorkspace: need 0)
0967: *              (Rworkspace: need N*N+2*M*N)
0968: *
0969:                NRWORK = IRVT
0970:                CALL CLACRM( M, N, U, LDU, RWORK( IRU ), N, A, LDA,
0971:      $                      RWORK( NRWORK ) )
0972:                CALL CLACPY( 'F', M, N, A, LDA, U, LDU )
0973:             ELSE
0974: *
0975: *              Copy A to VT, generate P**H
0976: *              (Cworkspace: need 2*N, prefer N+N*NB)
0977: *              (Rworkspace: 0)
0978: *
0979:                CALL CLACPY( 'U', N, N, A, LDA, VT, LDVT )
0980:                CALL CUNGBR( 'P', N, N, N, VT, LDVT, WORK( ITAUP ),
0981:      $                      WORK( NWORK ), LWORK-NWORK+1, IERR )
0982: *
0983: *              Copy A to U, generate Q
0984: *              (Cworkspace: need 2*N, prefer N+N*NB)
0985: *              (Rworkspace: 0)
0986: *
0987:                CALL CLACPY( 'L', M, N, A, LDA, U, LDU )
0988:                CALL CUNGBR( 'Q', M, M, N, U, LDU, WORK( ITAUQ ),
0989:      $                      WORK( NWORK ), LWORK-NWORK+1, IERR )
0990: *
0991: *              Perform bidiagonal SVD, computing left singular vectors
0992: *              of bidiagonal matrix in RWORK(IRU) and computing right
0993: *              singular vectors of bidiagonal matrix in RWORK(IRVT)
0994: *              (CWorkspace: need 0)
0995: *              (RWorkspace: need BDSPAC)
0996: *
0997:                IRU = NRWORK
0998:                IRVT = IRU + N*N
0999:                NRWORK = IRVT + N*N
1000:                CALL SBDSDC( 'U', 'I', N, S, RWORK( IE ), RWORK( IRU ),
1001:      $                      N, RWORK( IRVT ), N, DUM, IDUM,
1002:      $                      RWORK( NRWORK ), IWORK, INFO )
1003: *
1004: *              Multiply real matrix RWORK(IRVT) by P**H in VT,
1005: *              storing the result in A, copying to VT
1006: *              (Cworkspace: need 0)
1007: *              (Rworkspace: need 3*N*N)
1008: *
1009:                CALL CLARCM( N, N, RWORK( IRVT ), N, VT, LDVT, A, LDA,
1010:      $                      RWORK( NRWORK ) )
1011:                CALL CLACPY( 'F', N, N, A, LDA, VT, LDVT )
1012: *
1013: *              Multiply Q in U by real matrix RWORK(IRU), storing the
1014: *              result in A, copying to U
1015: *              (CWorkspace: 0)
1016: *              (Rworkspace: need 3*N*N)
1017: *
1018:                NRWORK = IRVT
1019:                CALL CLACRM( M, N, U, LDU, RWORK( IRU ), N, A, LDA,
1020:      $                      RWORK( NRWORK ) )
1021:                CALL CLACPY( 'F', M, N, A, LDA, U, LDU )
1022:             END IF
1023: *
1024:          ELSE
1025: *
1026: *           M .LT. MNTHR2
1027: *
1028: *           Path 6 (M at least N, but not much larger)
1029: *           Reduce to bidiagonal form without QR decomposition
1030: *           Use CUNMBR to compute singular vectors
1031: *
1032:             IE = 1
1033:             NRWORK = IE + N
1034:             ITAUQ = 1
1035:             ITAUP = ITAUQ + N
1036:             NWORK = ITAUP + N
1037: *
1038: *           Bidiagonalize A
1039: *           (CWorkspace: need 2*N+M, prefer 2*N+(M+N)*NB)
1040: *           (RWorkspace: need N)
1041: *
1042:             CALL CGEBRD( M, N, A, LDA, S, RWORK( IE ), WORK( ITAUQ ),
1043:      $                   WORK( ITAUP ), WORK( NWORK ), LWORK-NWORK+1,
1044:      $                   IERR )
1045:             IF( WNTQN ) THEN
1046: *
1047: *              Compute singular values only
1048: *              (Cworkspace: 0)
1049: *              (Rworkspace: need BDSPAN)
1050: *
1051:                CALL SBDSDC( 'U', 'N', N, S, RWORK( IE ), DUM, 1, DUM, 1,
1052:      $                      DUM, IDUM, RWORK( NRWORK ), IWORK, INFO )
1053:             ELSE IF( WNTQO ) THEN
1054:                IU = NWORK
1055:                IRU = NRWORK
1056:                IRVT = IRU + N*N
1057:                NRWORK = IRVT + N*N
1058:                IF( LWORK.GE.M*N+3*N ) THEN
1059: *
1060: *                 WORK( IU ) is M by N
1061: *
1062:                   LDWRKU = M
1063:                ELSE
1064: *
1065: *                 WORK( IU ) is LDWRKU by N
1066: *
1067:                   LDWRKU = ( LWORK-3*N ) / N
1068:                END IF
1069:                NWORK = IU + LDWRKU*N
1070: *
1071: *              Perform bidiagonal SVD, computing left singular vectors
1072: *              of bidiagonal matrix in RWORK(IRU) and computing right
1073: *              singular vectors of bidiagonal matrix in RWORK(IRVT)
1074: *              (CWorkspace: need 0)
1075: *              (RWorkspace: need BDSPAC)
1076: *
1077:                CALL SBDSDC( 'U', 'I', N, S, RWORK( IE ), RWORK( IRU ),
1078:      $                      N, RWORK( IRVT ), N, DUM, IDUM,
1079:      $                      RWORK( NRWORK ), IWORK, INFO )
1080: *
1081: *              Copy real matrix RWORK(IRVT) to complex matrix VT
1082: *              Overwrite VT by right singular vectors of A
1083: *              (Cworkspace: need 2*N, prefer N+N*NB)
1084: *              (Rworkspace: need 0)
1085: *
1086:                CALL CLACP2( 'F', N, N, RWORK( IRVT ), N, VT, LDVT )
1087:                CALL CUNMBR( 'P', 'R', 'C', N, N, N, A, LDA,
1088:      $                      WORK( ITAUP ), VT, LDVT, WORK( NWORK ),
1089:      $                      LWORK-NWORK+1, IERR )
1090: *
1091:                IF( LWORK.GE.M*N+3*N ) THEN
1092: *
1093: *              Copy real matrix RWORK(IRU) to complex matrix WORK(IU)
1094: *              Overwrite WORK(IU) by left singular vectors of A, copying
1095: *              to A
1096: *              (Cworkspace: need M*N+2*N, prefer M*N+N+N*NB)
1097: *              (Rworkspace: need 0)
1098: *
1099:                   CALL CLASET( 'F', M, N, CZERO, CZERO, WORK( IU ),
1100:      $                         LDWRKU )
1101:                   CALL CLACP2( 'F', N, N, RWORK( IRU ), N, WORK( IU ),
1102:      $                         LDWRKU )
1103:                   CALL CUNMBR( 'Q', 'L', 'N', M, N, N, A, LDA,
1104:      $                         WORK( ITAUQ ), WORK( IU ), LDWRKU,
1105:      $                         WORK( NWORK ), LWORK-NWORK+1, IERR )
1106:                   CALL CLACPY( 'F', M, N, WORK( IU ), LDWRKU, A, LDA )
1107:                ELSE
1108: *
1109: *                 Generate Q in A
1110: *                 (Cworkspace: need 2*N, prefer N+N*NB)
1111: *                 (Rworkspace: need 0)
1112: *
1113:                   CALL CUNGBR( 'Q', M, N, N, A, LDA, WORK( ITAUQ ),
1114:      $                         WORK( NWORK ), LWORK-NWORK+1, IERR )
1115: *
1116: *                 Multiply Q in A by real matrix RWORK(IRU), storing the
1117: *                 result in WORK(IU), copying to A
1118: *                 (CWorkspace: need N*N, prefer M*N)
1119: *                 (Rworkspace: need 3*N*N, prefer N*N+2*M*N)
1120: *
1121:                   NRWORK = IRVT
1122:                   DO 30 I = 1, M, LDWRKU
1123:                      CHUNK = MIN( M-I+1, LDWRKU )
1124:                      CALL CLACRM( CHUNK, N, A( I, 1 ), LDA,
1125:      $                            RWORK( IRU ), N, WORK( IU ), LDWRKU,
1126:      $                            RWORK( NRWORK ) )
1127:                      CALL CLACPY( 'F', CHUNK, N, WORK( IU ), LDWRKU,
1128:      $                            A( I, 1 ), LDA )
1129:    30             CONTINUE
1130:                END IF
1131: *
1132:             ELSE IF( WNTQS ) THEN
1133: *
1134: *              Perform bidiagonal SVD, computing left singular vectors
1135: *              of bidiagonal matrix in RWORK(IRU) and computing right
1136: *              singular vectors of bidiagonal matrix in RWORK(IRVT)
1137: *              (CWorkspace: need 0)
1138: *              (RWorkspace: need BDSPAC)
1139: *
1140:                IRU = NRWORK
1141:                IRVT = IRU + N*N
1142:                NRWORK = IRVT + N*N
1143:                CALL SBDSDC( 'U', 'I', N, S, RWORK( IE ), RWORK( IRU ),
1144:      $                      N, RWORK( IRVT ), N, DUM, IDUM,
1145:      $                      RWORK( NRWORK ), IWORK, INFO )
1146: *
1147: *              Copy real matrix RWORK(IRU) to complex matrix U
1148: *              Overwrite U by left singular vectors of A
1149: *              (CWorkspace: need 3*N, prefer 2*N+N*NB)
1150: *              (RWorkspace: 0)
1151: *
1152:                CALL CLASET( 'F', M, N, CZERO, CZERO, U, LDU )
1153:                CALL CLACP2( 'F', N, N, RWORK( IRU ), N, U, LDU )
1154:                CALL CUNMBR( 'Q', 'L', 'N', M, N, N, A, LDA,
1155:      $                      WORK( ITAUQ ), U, LDU, WORK( NWORK ),
1156:      $                      LWORK-NWORK+1, IERR )
1157: *
1158: *              Copy real matrix RWORK(IRVT) to complex matrix VT
1159: *              Overwrite VT by right singular vectors of A
1160: *              (CWorkspace: need 3*N, prefer 2*N+N*NB)
1161: *              (RWorkspace: 0)
1162: *
1163:                CALL CLACP2( 'F', N, N, RWORK( IRVT ), N, VT, LDVT )
1164:                CALL CUNMBR( 'P', 'R', 'C', N, N, N, A, LDA,
1165:      $                      WORK( ITAUP ), VT, LDVT, WORK( NWORK ),
1166:      $                      LWORK-NWORK+1, IERR )
1167:             ELSE
1168: *
1169: *              Perform bidiagonal SVD, computing left singular vectors
1170: *              of bidiagonal matrix in RWORK(IRU) and computing right
1171: *              singular vectors of bidiagonal matrix in RWORK(IRVT)
1172: *              (CWorkspace: need 0)
1173: *              (RWorkspace: need BDSPAC)
1174: *
1175:                IRU = NRWORK
1176:                IRVT = IRU + N*N
1177:                NRWORK = IRVT + N*N
1178:                CALL SBDSDC( 'U', 'I', N, S, RWORK( IE ), RWORK( IRU ),
1179:      $                      N, RWORK( IRVT ), N, DUM, IDUM,
1180:      $                      RWORK( NRWORK ), IWORK, INFO )
1181: *
1182: *              Set the right corner of U to identity matrix
1183: *
1184:                CALL CLASET( 'F', M, M, CZERO, CZERO, U, LDU )
1185:                IF( M.GT.N ) THEN
1186:                   CALL CLASET( 'F', M-N, M-N, CZERO, CONE,
1187:      $                         U( N+1, N+1 ), LDU )
1188:                END IF
1189: *
1190: *              Copy real matrix RWORK(IRU) to complex matrix U
1191: *              Overwrite U by left singular vectors of A
1192: *              (CWorkspace: need 2*N+M, prefer 2*N+M*NB)
1193: *              (RWorkspace: 0)
1194: *
1195:                CALL CLACP2( 'F', N, N, RWORK( IRU ), N, U, LDU )
1196:                CALL CUNMBR( 'Q', 'L', 'N', M, M, N, A, LDA,
1197:      $                      WORK( ITAUQ ), U, LDU, WORK( NWORK ),
1198:      $                      LWORK-NWORK+1, IERR )
1199: *
1200: *              Copy real matrix RWORK(IRVT) to complex matrix VT
1201: *              Overwrite VT by right singular vectors of A
1202: *              (CWorkspace: need 3*N, prefer 2*N+N*NB)
1203: *              (RWorkspace: 0)
1204: *
1205:                CALL CLACP2( 'F', N, N, RWORK( IRVT ), N, VT, LDVT )
1206:                CALL CUNMBR( 'P', 'R', 'C', N, N, N, A, LDA,
1207:      $                      WORK( ITAUP ), VT, LDVT, WORK( NWORK ),
1208:      $                      LWORK-NWORK+1, IERR )
1209:             END IF
1210: *
1211:          END IF
1212: *
1213:       ELSE
1214: *
1215: *        A has more columns than rows. If A has sufficiently more
1216: *        columns than rows, first reduce using the LQ decomposition (if
1217: *        sufficient workspace available)
1218: *
1219:          IF( N.GE.MNTHR1 ) THEN
1220: *
1221:             IF( WNTQN ) THEN
1222: *
1223: *              Path 1t (N much larger than M, JOBZ='N')
1224: *              No singular vectors to be computed
1225: *
1226:                ITAU = 1
1227:                NWORK = ITAU + M
1228: *
1229: *              Compute A=L*Q
1230: *              (CWorkspace: need 2*M, prefer M+M*NB)
1231: *              (RWorkspace: 0)
1232: *
1233:                CALL CGELQF( M, N, A, LDA, WORK( ITAU ), WORK( NWORK ),
1234:      $                      LWORK-NWORK+1, IERR )
1235: *
1236: *              Zero out above L
1237: *
1238:                CALL CLASET( 'U', M-1, M-1, CZERO, CZERO, A( 1, 2 ),
1239:      $                      LDA )
1240:                IE = 1
1241:                ITAUQ = 1
1242:                ITAUP = ITAUQ + M
1243:                NWORK = ITAUP + M
1244: *
1245: *              Bidiagonalize L in A
1246: *              (CWorkspace: need 3*M, prefer 2*M+2*M*NB)
1247: *              (RWorkspace: need M)
1248: *
1249:                CALL CGEBRD( M, M, A, LDA, S, RWORK( IE ), WORK( ITAUQ ),
1250:      $                      WORK( ITAUP ), WORK( NWORK ), LWORK-NWORK+1,
1251:      $                      IERR )
1252:                NRWORK = IE + M
1253: *
1254: *              Perform bidiagonal SVD, compute singular values only
1255: *              (CWorkspace: 0)
1256: *              (RWorkspace: need BDSPAN)
1257: *
1258:                CALL SBDSDC( 'U', 'N', M, S, RWORK( IE ), DUM, 1, DUM, 1,
1259:      $                      DUM, IDUM, RWORK( NRWORK ), IWORK, INFO )
1260: *
1261:             ELSE IF( WNTQO ) THEN
1262: *
1263: *              Path 2t (N much larger than M, JOBZ='O')
1264: *              M right singular vectors to be overwritten on A and
1265: *              M left singular vectors to be computed in U
1266: *
1267:                IVT = 1
1268:                LDWKVT = M
1269: *
1270: *              WORK(IVT) is M by M
1271: *
1272:                IL = IVT + LDWKVT*M
1273:                IF( LWORK.GE.M*N+M*M+3*M ) THEN
1274: *
1275: *                 WORK(IL) M by N
1276: *
1277:                   LDWRKL = M
1278:                   CHUNK = N
1279:                ELSE
1280: *
1281: *                 WORK(IL) is M by CHUNK
1282: *
1283:                   LDWRKL = M
1284:                   CHUNK = ( LWORK-M*M-3*M ) / M
1285:                END IF
1286:                ITAU = IL + LDWRKL*CHUNK
1287:                NWORK = ITAU + M
1288: *
1289: *              Compute A=L*Q
1290: *              (CWorkspace: need 2*M, prefer M+M*NB)
1291: *              (RWorkspace: 0)
1292: *
1293:                CALL CGELQF( M, N, A, LDA, WORK( ITAU ), WORK( NWORK ),
1294:      $                      LWORK-NWORK+1, IERR )
1295: *
1296: *              Copy L to WORK(IL), zeroing about above it
1297: *
1298:                CALL CLACPY( 'L', M, M, A, LDA, WORK( IL ), LDWRKL )
1299:                CALL CLASET( 'U', M-1, M-1, CZERO, CZERO,
1300:      $                      WORK( IL+LDWRKL ), LDWRKL )
1301: *
1302: *              Generate Q in A
1303: *              (CWorkspace: need M*M+2*M, prefer M*M+M+M*NB)
1304: *              (RWorkspace: 0)
1305: *
1306:                CALL CUNGLQ( M, N, M, A, LDA, WORK( ITAU ),
1307:      $                      WORK( NWORK ), LWORK-NWORK+1, IERR )
1308:                IE = 1
1309:                ITAUQ = ITAU
1310:                ITAUP = ITAUQ + M
1311:                NWORK = ITAUP + M
1312: *
1313: *              Bidiagonalize L in WORK(IL)
1314: *              (CWorkspace: need M*M+3*M, prefer M*M+2*M+2*M*NB)
1315: *              (RWorkspace: need M)
1316: *
1317:                CALL CGEBRD( M, M, WORK( IL ), LDWRKL, S, RWORK( IE ),
1318:      $                      WORK( ITAUQ ), WORK( ITAUP ), WORK( NWORK ),
1319:      $                      LWORK-NWORK+1, IERR )
1320: *
1321: *              Perform bidiagonal SVD, computing left singular vectors
1322: *              of bidiagonal matrix in RWORK(IRU) and computing right
1323: *              singular vectors of bidiagonal matrix in RWORK(IRVT)
1324: *              (CWorkspace: need 0)
1325: *              (RWorkspace: need BDSPAC)
1326: *
1327:                IRU = IE + M
1328:                IRVT = IRU + M*M
1329:                NRWORK = IRVT + M*M
1330:                CALL SBDSDC( 'U', 'I', M, S, RWORK( IE ), RWORK( IRU ),
1331:      $                      M, RWORK( IRVT ), M, DUM, IDUM,
1332:      $                      RWORK( NRWORK ), IWORK, INFO )
1333: *
1334: *              Copy real matrix RWORK(IRU) to complex matrix WORK(IU)
1335: *              Overwrite WORK(IU) by the left singular vectors of L
1336: *              (CWorkspace: need N*N+3*N, prefer M*N+2*N+N*NB)
1337: *              (RWorkspace: 0)
1338: *
1339:                CALL CLACP2( 'F', M, M, RWORK( IRU ), M, U, LDU )
1340:                CALL CUNMBR( 'Q', 'L', 'N', M, M, M, WORK( IL ), LDWRKL,
1341:      $                      WORK( ITAUQ ), U, LDU, WORK( NWORK ),
1342:      $                      LWORK-NWORK+1, IERR )
1343: *
1344: *              Copy real matrix RWORK(IRVT) to complex matrix WORK(IVT)
1345: *              Overwrite WORK(IVT) by the right singular vectors of L
1346: *              (CWorkspace: need N*N+3*N, prefer M*N+2*N+N*NB)
1347: *              (RWorkspace: 0)
1348: *
1349:                CALL CLACP2( 'F', M, M, RWORK( IRVT ), M, WORK( IVT ),
1350:      $                      LDWKVT )
1351:                CALL CUNMBR( 'P', 'R', 'C', M, M, M, WORK( IL ), LDWRKL,
1352:      $                      WORK( ITAUP ), WORK( IVT ), LDWKVT,
1353:      $                      WORK( NWORK ), LWORK-NWORK+1, IERR )
1354: *
1355: *              Multiply right singular vectors of L in WORK(IL) by Q
1356: *              in A, storing result in WORK(IL) and copying to A
1357: *              (CWorkspace: need 2*M*M, prefer M*M+M*N))
1358: *              (RWorkspace: 0)
1359: *
1360:                DO 40 I = 1, N, CHUNK
1361:                   BLK = MIN( N-I+1, CHUNK )
1362:                   CALL CGEMM( 'N', 'N', M, BLK, M, CONE, WORK( IVT ), M,
1363:      $                        A( 1, I ), LDA, CZERO, WORK( IL ),
1364:      $                        LDWRKL )
1365:                   CALL CLACPY( 'F', M, BLK, WORK( IL ), LDWRKL,
1366:      $                         A( 1, I ), LDA )
1367:    40          CONTINUE
1368: *
1369:             ELSE IF( WNTQS ) THEN
1370: *
1371: *             Path 3t (N much larger than M, JOBZ='S')
1372: *             M right singular vectors to be computed in VT and
1373: *             M left singular vectors to be computed in U
1374: *
1375:                IL = 1
1376: *
1377: *              WORK(IL) is M by M
1378: *
1379:                LDWRKL = M
1380:                ITAU = IL + LDWRKL*M
1381:                NWORK = ITAU + M
1382: *
1383: *              Compute A=L*Q
1384: *              (CWorkspace: need 2*M, prefer M+M*NB)
1385: *              (RWorkspace: 0)
1386: *
1387:                CALL CGELQF( M, N, A, LDA, WORK( ITAU ), WORK( NWORK ),
1388:      $                      LWORK-NWORK+1, IERR )
1389: *
1390: *              Copy L to WORK(IL), zeroing out above it
1391: *
1392:                CALL CLACPY( 'L', M, M, A, LDA, WORK( IL ), LDWRKL )
1393:                CALL CLASET( 'U', M-1, M-1, CZERO, CZERO,
1394:      $                      WORK( IL+LDWRKL ), LDWRKL )
1395: *
1396: *              Generate Q in A
1397: *              (CWorkspace: need M*M+2*M, prefer M*M+M+M*NB)
1398: *              (RWorkspace: 0)
1399: *
1400:                CALL CUNGLQ( M, N, M, A, LDA, WORK( ITAU ),
1401:      $                      WORK( NWORK ), LWORK-NWORK+1, IERR )
1402:                IE = 1
1403:                ITAUQ = ITAU
1404:                ITAUP = ITAUQ + M
1405:                NWORK = ITAUP + M
1406: *
1407: *              Bidiagonalize L in WORK(IL)
1408: *              (CWorkspace: need M*M+3*M, prefer M*M+2*M+2*M*NB)
1409: *              (RWorkspace: need M)
1410: *
1411:                CALL CGEBRD( M, M, WORK( IL ), LDWRKL, S, RWORK( IE ),
1412:      $                      WORK( ITAUQ ), WORK( ITAUP ), WORK( NWORK ),
1413:      $                      LWORK-NWORK+1, IERR )
1414: *
1415: *              Perform bidiagonal SVD, computing left singular vectors
1416: *              of bidiagonal matrix in RWORK(IRU) and computing right
1417: *              singular vectors of bidiagonal matrix in RWORK(IRVT)
1418: *              (CWorkspace: need 0)
1419: *              (RWorkspace: need BDSPAC)
1420: *
1421:                IRU = IE + M
1422:                IRVT = IRU + M*M
1423:                NRWORK = IRVT + M*M
1424:                CALL SBDSDC( 'U', 'I', M, S, RWORK( IE ), RWORK( IRU ),
1425:      $                      M, RWORK( IRVT ), M, DUM, IDUM,
1426:      $                      RWORK( NRWORK ), IWORK, INFO )
1427: *
1428: *              Copy real matrix RWORK(IRU) to complex matrix U
1429: *              Overwrite U by left singular vectors of L
1430: *              (CWorkspace: need M*M+3*M, prefer M*M+2*M+M*NB)
1431: *              (RWorkspace: 0)
1432: *
1433:                CALL CLACP2( 'F', M, M, RWORK( IRU ), M, U, LDU )
1434:                CALL CUNMBR( 'Q', 'L', 'N', M, M, M, WORK( IL ), LDWRKL,
1435:      $                      WORK( ITAUQ ), U, LDU, WORK( NWORK ),
1436:      $                      LWORK-NWORK+1, IERR )
1437: *
1438: *              Copy real matrix RWORK(IRVT) to complex matrix VT
1439: *              Overwrite VT by left singular vectors of L
1440: *              (CWorkspace: need M*M+3*M, prefer M*M+2*M+M*NB)
1441: *              (RWorkspace: 0)
1442: *
1443:                CALL CLACP2( 'F', M, M, RWORK( IRVT ), M, VT, LDVT )
1444:                CALL CUNMBR( 'P', 'R', 'C', M, M, M, WORK( IL ), LDWRKL,
1445:      $                      WORK( ITAUP ), VT, LDVT, WORK( NWORK ),
1446:      $                      LWORK-NWORK+1, IERR )
1447: *
1448: *              Copy VT to WORK(IL), multiply right singular vectors of L
1449: *              in WORK(IL) by Q in A, storing result in VT
1450: *              (CWorkspace: need M*M)
1451: *              (RWorkspace: 0)
1452: *
1453:                CALL CLACPY( 'F', M, M, VT, LDVT, WORK( IL ), LDWRKL )
1454:                CALL CGEMM( 'N', 'N', M, N, M, CONE, WORK( IL ), LDWRKL,
1455:      $                     A, LDA, CZERO, VT, LDVT )
1456: *
1457:             ELSE IF( WNTQA ) THEN
1458: *
1459: *              Path 9t (N much larger than M, JOBZ='A')
1460: *              N right singular vectors to be computed in VT and
1461: *              M left singular vectors to be computed in U
1462: *
1463:                IVT = 1
1464: *
1465: *              WORK(IVT) is M by M
1466: *
1467:                LDWKVT = M
1468:                ITAU = IVT + LDWKVT*M
1469:                NWORK = ITAU + M
1470: *
1471: *              Compute A=L*Q, copying result to VT
1472: *              (CWorkspace: need 2*M, prefer M+M*NB)
1473: *              (RWorkspace: 0)
1474: *
1475:                CALL CGELQF( M, N, A, LDA, WORK( ITAU ), WORK( NWORK ),
1476:      $                      LWORK-NWORK+1, IERR )
1477:                CALL CLACPY( 'U', M, N, A, LDA, VT, LDVT )
1478: *
1479: *              Generate Q in VT
1480: *              (CWorkspace: need M+N, prefer M+N*NB)
1481: *              (RWorkspace: 0)
1482: *
1483:                CALL CUNGLQ( N, N, M, VT, LDVT, WORK( ITAU ),
1484:      $                      WORK( NWORK ), LWORK-NWORK+1, IERR )
1485: *
1486: *              Produce L in A, zeroing out above it
1487: *
1488:                CALL CLASET( 'U', M-1, M-1, CZERO, CZERO, A( 1, 2 ),
1489:      $                      LDA )
1490:                IE = 1
1491:                ITAUQ = ITAU
1492:                ITAUP = ITAUQ + M
1493:                NWORK = ITAUP + M
1494: *
1495: *              Bidiagonalize L in A
1496: *              (CWorkspace: need M*M+3*M, prefer M*M+2*M+2*M*NB)
1497: *              (RWorkspace: need M)
1498: *
1499:                CALL CGEBRD( M, M, A, LDA, S, RWORK( IE ), WORK( ITAUQ ),
1500:      $                      WORK( ITAUP ), WORK( NWORK ), LWORK-NWORK+1,
1501:      $                      IERR )
1502: *
1503: *              Perform bidiagonal SVD, computing left singular vectors
1504: *              of bidiagonal matrix in RWORK(IRU) and computing right
1505: *              singular vectors of bidiagonal matrix in RWORK(IRVT)
1506: *              (CWorkspace: need 0)
1507: *              (RWorkspace: need BDSPAC)
1508: *
1509:                IRU = IE + M
1510:                IRVT = IRU + M*M
1511:                NRWORK = IRVT + M*M
1512:                CALL SBDSDC( 'U', 'I', M, S, RWORK( IE ), RWORK( IRU ),
1513:      $                      M, RWORK( IRVT ), M, DUM, IDUM,
1514:      $                      RWORK( NRWORK ), IWORK, INFO )
1515: *
1516: *              Copy real matrix RWORK(IRU) to complex matrix U
1517: *              Overwrite U by left singular vectors of L
1518: *              (CWorkspace: need 3*M, prefer 2*M+M*NB)
1519: *              (RWorkspace: 0)
1520: *
1521:                CALL CLACP2( 'F', M, M, RWORK( IRU ), M, U, LDU )
1522:                CALL CUNMBR( 'Q', 'L', 'N', M, M, M, A, LDA,
1523:      $                      WORK( ITAUQ ), U, LDU, WORK( NWORK ),
1524:      $                      LWORK-NWORK+1, IERR )
1525: *
1526: *              Copy real matrix RWORK(IRVT) to complex matrix WORK(IVT)
1527: *              Overwrite WORK(IVT) by right singular vectors of L
1528: *              (CWorkspace: need M*M+3*M, prefer M*M+2*M+M*NB)
1529: *              (RWorkspace: 0)
1530: *
1531:                CALL CLACP2( 'F', M, M, RWORK( IRVT ), M, WORK( IVT ),
1532:      $                      LDWKVT )
1533:                CALL CUNMBR( 'P', 'R', 'C', M, M, M, A, LDA,
1534:      $                      WORK( ITAUP ), WORK( IVT ), LDWKVT,
1535:      $                      WORK( NWORK ), LWORK-NWORK+1, IERR )
1536: *
1537: *              Multiply right singular vectors of L in WORK(IVT) by
1538: *              Q in VT, storing result in A
1539: *              (CWorkspace: need M*M)
1540: *              (RWorkspace: 0)
1541: *
1542:                CALL CGEMM( 'N', 'N', M, N, M, CONE, WORK( IVT ),
1543:      $                     LDWKVT, VT, LDVT, CZERO, A, LDA )
1544: *
1545: *              Copy right singular vectors of A from A to VT
1546: *
1547:                CALL CLACPY( 'F', M, N, A, LDA, VT, LDVT )
1548: *
1549:             END IF
1550: *
1551:          ELSE IF( N.GE.MNTHR2 ) THEN
1552: *
1553: *           MNTHR2 <= N < MNTHR1
1554: *
1555: *           Path 5t (N much larger than M, but not as much as MNTHR1)
1556: *           Reduce to bidiagonal form without QR decomposition, use
1557: *           CUNGBR and matrix multiplication to compute singular vectors
1558: *
1559: *
1560:             IE = 1
1561:             NRWORK = IE + M
1562:             ITAUQ = 1
1563:             ITAUP = ITAUQ + M
1564:             NWORK = ITAUP + M
1565: *
1566: *           Bidiagonalize A
1567: *           (CWorkspace: need 2*M+N, prefer 2*M+(M+N)*NB)
1568: *           (RWorkspace: M)
1569: *
1570:             CALL CGEBRD( M, N, A, LDA, S, RWORK( IE ), WORK( ITAUQ ),
1571:      $                   WORK( ITAUP ), WORK( NWORK ), LWORK-NWORK+1,
1572:      $                   IERR )
1573: *
1574:             IF( WNTQN ) THEN
1575: *
1576: *              Compute singular values only
1577: *              (Cworkspace: 0)
1578: *              (Rworkspace: need BDSPAN)
1579: *
1580:                CALL SBDSDC( 'L', 'N', M, S, RWORK( IE ), DUM, 1, DUM, 1,
1581:      $                      DUM, IDUM, RWORK( NRWORK ), IWORK, INFO )
1582:             ELSE IF( WNTQO ) THEN
1583:                IRVT = NRWORK
1584:                IRU = IRVT + M*M
1585:                NRWORK = IRU + M*M
1586:                IVT = NWORK
1587: *
1588: *              Copy A to U, generate Q
1589: *              (Cworkspace: need 2*M, prefer M+M*NB)
1590: *              (Rworkspace: 0)
1591: *
1592:                CALL CLACPY( 'L', M, M, A, LDA, U, LDU )
1593:                CALL CUNGBR( 'Q', M, M, N, U, LDU, WORK( ITAUQ ),
1594:      $                      WORK( NWORK ), LWORK-NWORK+1, IERR )
1595: *
1596: *              Generate P**H in A
1597: *              (Cworkspace: need 2*M, prefer M+M*NB)
1598: *              (Rworkspace: 0)
1599: *
1600:                CALL CUNGBR( 'P', M, N, M, A, LDA, WORK( ITAUP ),
1601:      $                      WORK( NWORK ), LWORK-NWORK+1, IERR )
1602: *
1603:                LDWKVT = M
1604:                IF( LWORK.GE.M*N+3*M ) THEN
1605: *
1606: *                 WORK( IVT ) is M by N
1607: *
1608:                   NWORK = IVT + LDWKVT*N
1609:                   CHUNK = N
1610:                ELSE
1611: *
1612: *                 WORK( IVT ) is M by CHUNK
1613: *
1614:                   CHUNK = ( LWORK-3*M ) / M
1615:                   NWORK = IVT + LDWKVT*CHUNK
1616:                END IF
1617: *
1618: *              Perform bidiagonal SVD, computing left singular vectors
1619: *              of bidiagonal matrix in RWORK(IRU) and computing right
1620: *              singular vectors of bidiagonal matrix in RWORK(IRVT)
1621: *              (CWorkspace: need 0)
1622: *              (RWorkspace: need BDSPAC)
1623: *
1624:                CALL SBDSDC( 'L', 'I', M, S, RWORK( IE ), RWORK( IRU ),
1625:      $                      M, RWORK( IRVT ), M, DUM, IDUM,
1626:      $                      RWORK( NRWORK ), IWORK, INFO )
1627: *
1628: *              Multiply Q in U by real matrix RWORK(IRVT)
1629: *              storing the result in WORK(IVT), copying to U
1630: *              (Cworkspace: need 0)
1631: *              (Rworkspace: need 2*M*M)
1632: *
1633:                CALL CLACRM( M, M, U, LDU, RWORK( IRU ), M, WORK( IVT ),
1634:      $                      LDWKVT, RWORK( NRWORK ) )
1635:                CALL CLACPY( 'F', M, M, WORK( IVT ), LDWKVT, U, LDU )
1636: *
1637: *              Multiply RWORK(IRVT) by P**H in A, storing the
1638: *              result in WORK(IVT), copying to A
1639: *              (CWorkspace: need M*M, prefer M*N)
1640: *              (Rworkspace: need 2*M*M, prefer 2*M*N)
1641: *
1642:                NRWORK = IRU
1643:                DO 50 I = 1, N, CHUNK
1644:                   BLK = MIN( N-I+1, CHUNK )
1645:                   CALL CLARCM( M, BLK, RWORK( IRVT ), M, A( 1, I ), LDA,
1646:      $                         WORK( IVT ), LDWKVT, RWORK( NRWORK ) )
1647:                   CALL CLACPY( 'F', M, BLK, WORK( IVT ), LDWKVT,
1648:      $                         A( 1, I ), LDA )
1649:    50          CONTINUE
1650:             ELSE IF( WNTQS ) THEN
1651: *
1652: *              Copy A to U, generate Q
1653: *              (Cworkspace: need 2*M, prefer M+M*NB)
1654: *              (Rworkspace: 0)
1655: *
1656:                CALL CLACPY( 'L', M, M, A, LDA, U, LDU )
1657:                CALL CUNGBR( 'Q', M, M, N, U, LDU, WORK( ITAUQ ),
1658:      $                      WORK( NWORK ), LWORK-NWORK+1, IERR )
1659: *
1660: *              Copy A to VT, generate P**H
1661: *              (Cworkspace: need 2*M, prefer M+M*NB)
1662: *              (Rworkspace: 0)
1663: *
1664:                CALL CLACPY( 'U', M, N, A, LDA, VT, LDVT )
1665:                CALL CUNGBR( 'P', M, N, M, VT, LDVT, WORK( ITAUP ),
1666:      $                      WORK( NWORK ), LWORK-NWORK+1, IERR )
1667: *
1668: *              Perform bidiagonal SVD, computing left singular vectors
1669: *              of bidiagonal matrix in RWORK(IRU) and computing right
1670: *              singular vectors of bidiagonal matrix in RWORK(IRVT)
1671: *              (CWorkspace: need 0)
1672: *              (RWorkspace: need BDSPAC)
1673: *
1674:                IRVT = NRWORK
1675:                IRU = IRVT + M*M
1676:                NRWORK = IRU + M*M
1677:                CALL SBDSDC( 'L', 'I', M, S, RWORK( IE ), RWORK( IRU ),
1678:      $                      M, RWORK( IRVT ), M, DUM, IDUM,
1679:      $                      RWORK( NRWORK ), IWORK, INFO )
1680: *
1681: *              Multiply Q in U by real matrix RWORK(IRU), storing the
1682: *              result in A, copying to U
1683: *              (CWorkspace: need 0)
1684: *              (Rworkspace: need 3*M*M)
1685: *
1686:                CALL CLACRM( M, M, U, LDU, RWORK( IRU ), M, A, LDA,
1687:      $                      RWORK( NRWORK ) )
1688:                CALL CLACPY( 'F', M, M, A, LDA, U, LDU )
1689: *
1690: *              Multiply real matrix RWORK(IRVT) by P**H in VT,
1691: *              storing the result in A, copying to VT
1692: *              (Cworkspace: need 0)
1693: *              (Rworkspace: need M*M+2*M*N)
1694: *
1695:                NRWORK = IRU
1696:                CALL CLARCM( M, N, RWORK( IRVT ), M, VT, LDVT, A, LDA,
1697:      $                      RWORK( NRWORK ) )
1698:                CALL CLACPY( 'F', M, N, A, LDA, VT, LDVT )
1699:             ELSE
1700: *
1701: *              Copy A to U, generate Q
1702: *              (Cworkspace: need 2*M, prefer M+M*NB)
1703: *              (Rworkspace: 0)
1704: *
1705:                CALL CLACPY( 'L', M, M, A, LDA, U, LDU )
1706:                CALL CUNGBR( 'Q', M, M, N, U, LDU, WORK( ITAUQ ),
1707:      $                      WORK( NWORK ), LWORK-NWORK+1, IERR )
1708: *
1709: *              Copy A to VT, generate P**H
1710: *              (Cworkspace: need 2*M, prefer M+M*NB)
1711: *              (Rworkspace: 0)
1712: *
1713:                CALL CLACPY( 'U', M, N, A, LDA, VT, LDVT )
1714:                CALL CUNGBR( 'P', N, N, M, VT, LDVT, WORK( ITAUP ),
1715:      $                      WORK( NWORK ), LWORK-NWORK+1, IERR )
1716: *
1717: *              Perform bidiagonal SVD, computing left singular vectors
1718: *              of bidiagonal matrix in RWORK(IRU) and computing right
1719: *              singular vectors of bidiagonal matrix in RWORK(IRVT)
1720: *              (CWorkspace: need 0)
1721: *              (RWorkspace: need BDSPAC)
1722: *
1723:                IRVT = NRWORK
1724:                IRU = IRVT + M*M
1725:                NRWORK = IRU + M*M
1726:                CALL SBDSDC( 'L', 'I', M, S, RWORK( IE ), RWORK( IRU ),
1727:      $                      M, RWORK( IRVT ), M, DUM, IDUM,
1728:      $                      RWORK( NRWORK ), IWORK, INFO )
1729: *
1730: *              Multiply Q in U by real matrix RWORK(IRU), storing the
1731: *              result in A, copying to U
1732: *              (CWorkspace: need 0)
1733: *              (Rworkspace: need 3*M*M)
1734: *
1735:                CALL CLACRM( M, M, U, LDU, RWORK( IRU ), M, A, LDA,
1736:      $                      RWORK( NRWORK ) )
1737:                CALL CLACPY( 'F', M, M, A, LDA, U, LDU )
1738: *
1739: *              Multiply real matrix RWORK(IRVT) by P**H in VT,
1740: *              storing the result in A, copying to VT
1741: *              (Cworkspace: need 0)
1742: *              (Rworkspace: need M*M+2*M*N)
1743: *
1744:                CALL CLARCM( M, N, RWORK( IRVT ), M, VT, LDVT, A, LDA,
1745:      $                      RWORK( NRWORK ) )
1746:                CALL CLACPY( 'F', M, N, A, LDA, VT, LDVT )
1747:             END IF
1748: *
1749:          ELSE
1750: *
1751: *           N .LT. MNTHR2
1752: *
1753: *           Path 6t (N greater than M, but not much larger)
1754: *           Reduce to bidiagonal form without LQ decomposition
1755: *           Use CUNMBR to compute singular vectors
1756: *
1757:             IE = 1
1758:             NRWORK = IE + M
1759:             ITAUQ = 1
1760:             ITAUP = ITAUQ + M
1761:             NWORK = ITAUP + M
1762: *
1763: *           Bidiagonalize A
1764: *           (CWorkspace: need 2*M+N, prefer 2*M+(M+N)*NB)
1765: *           (RWorkspace: M)
1766: *
1767:             CALL CGEBRD( M, N, A, LDA, S, RWORK( IE ), WORK( ITAUQ ),
1768:      $                   WORK( ITAUP ), WORK( NWORK ), LWORK-NWORK+1,
1769:      $                   IERR )
1770:             IF( WNTQN ) THEN
1771: *
1772: *              Compute singular values only
1773: *              (Cworkspace: 0)
1774: *              (Rworkspace: need BDSPAN)
1775: *
1776:                CALL SBDSDC( 'L', 'N', M, S, RWORK( IE ), DUM, 1, DUM, 1,
1777:      $                      DUM, IDUM, RWORK( NRWORK ), IWORK, INFO )
1778:             ELSE IF( WNTQO ) THEN
1779:                LDWKVT = M
1780:                IVT = NWORK
1781:                IF( LWORK.GE.M*N+3*M ) THEN
1782: *
1783: *                 WORK( IVT ) is M by N
1784: *
1785:                   CALL CLASET( 'F', M, N, CZERO, CZERO, WORK( IVT ),
1786:      $                         LDWKVT )
1787:                   NWORK = IVT + LDWKVT*N
1788:                ELSE
1789: *
1790: *                 WORK( IVT ) is M by CHUNK
1791: *
1792:                   CHUNK = ( LWORK-3*M ) / M
1793:                   NWORK = IVT + LDWKVT*CHUNK
1794:                END IF
1795: *
1796: *              Perform bidiagonal SVD, computing left singular vectors
1797: *              of bidiagonal matrix in RWORK(IRU) and computing right
1798: *              singular vectors of bidiagonal matrix in RWORK(IRVT)
1799: *              (CWorkspace: need 0)
1800: *              (RWorkspace: need BDSPAC)
1801: *
1802:                IRVT = NRWORK
1803:                IRU = IRVT + M*M
1804:                NRWORK = IRU + M*M
1805:                CALL SBDSDC( 'L', 'I', M, S, RWORK( IE ), RWORK( IRU ),
1806:      $                      M, RWORK( IRVT ), M, DUM, IDUM,
1807:      $                      RWORK( NRWORK ), IWORK, INFO )
1808: *
1809: *              Copy real matrix RWORK(IRU) to complex matrix U
1810: *              Overwrite U by left singular vectors of A
1811: *              (Cworkspace: need 2*M, prefer M+M*NB)
1812: *              (Rworkspace: need 0)
1813: *
1814:                CALL CLACP2( 'F', M, M, RWORK( IRU ), M, U, LDU )
1815:                CALL CUNMBR( 'Q', 'L', 'N', M, M, N, A, LDA,
1816:      $                      WORK( ITAUQ ), U, LDU, WORK( NWORK ),
1817:      $                      LWORK-NWORK+1, IERR )
1818: *
1819:                IF( LWORK.GE.M*N+3*M ) THEN
1820: *
1821: *              Copy real matrix RWORK(IRVT) to complex matrix WORK(IVT)
1822: *              Overwrite WORK(IVT) by right singular vectors of A,
1823: *              copying to A
1824: *              (Cworkspace: need M*N+2*M, prefer M*N+M+M*NB)
1825: *              (Rworkspace: need 0)
1826: *
1827:                   CALL CLACP2( 'F', M, M, RWORK( IRVT ), M, WORK( IVT ),
1828:      $                         LDWKVT )
1829:                   CALL CUNMBR( 'P', 'R', 'C', M, N, M, A, LDA,
1830:      $                         WORK( ITAUP ), WORK( IVT ), LDWKVT,
1831:      $                         WORK( NWORK ), LWORK-NWORK+1, IERR )
1832:                   CALL CLACPY( 'F', M, N, WORK( IVT ), LDWKVT, A, LDA )
1833:                ELSE
1834: *
1835: *                 Generate P**H in A
1836: *                 (Cworkspace: need 2*M, prefer M+M*NB)
1837: *                 (Rworkspace: need 0)
1838: *
1839:                   CALL CUNGBR( 'P', M, N, M, A, LDA, WORK( ITAUP ),
1840:      $                         WORK( NWORK ), LWORK-NWORK+1, IERR )
1841: *
1842: *                 Multiply Q in A by real matrix RWORK(IRU), storing the
1843: *                 result in WORK(IU), copying to A
1844: *                 (CWorkspace: need M*M, prefer M*N)
1845: *                 (Rworkspace: need 3*M*M, prefer M*M+2*M*N)
1846: *
1847:                   NRWORK = IRU
1848:                   DO 60 I = 1, N, CHUNK
1849:                      BLK = MIN( N-I+1, CHUNK )
1850:                      CALL CLARCM( M, BLK, RWORK( IRVT ), M, A( 1, I ),
1851:      $                            LDA, WORK( IVT ), LDWKVT,
1852:      $                            RWORK( NRWORK ) )
1853:                      CALL CLACPY( 'F', M, BLK, WORK( IVT ), LDWKVT,
1854:      $                            A( 1, I ), LDA )
1855:    60             CONTINUE
1856:                END IF
1857:             ELSE IF( WNTQS ) THEN
1858: *
1859: *              Perform bidiagonal SVD, computing left singular vectors
1860: *              of bidiagonal matrix in RWORK(IRU) and computing right
1861: *              singular vectors of bidiagonal matrix in RWORK(IRVT)
1862: *              (CWorkspace: need 0)
1863: *              (RWorkspace: need BDSPAC)
1864: *
1865:                IRVT = NRWORK
1866:                IRU = IRVT + M*M
1867:                NRWORK = IRU + M*M
1868:                CALL SBDSDC( 'L', 'I', M, S, RWORK( IE ), RWORK( IRU ),
1869:      $                      M, RWORK( IRVT ), M, DUM, IDUM,
1870:      $                      RWORK( NRWORK ), IWORK, INFO )
1871: *
1872: *              Copy real matrix RWORK(IRU) to complex matrix U
1873: *              Overwrite U by left singular vectors of A
1874: *              (CWorkspace: need 3*M, prefer 2*M+M*NB)
1875: *              (RWorkspace: M*M)
1876: *
1877:                CALL CLACP2( 'F', M, M, RWORK( IRU ), M, U, LDU )
1878:                CALL CUNMBR( 'Q', 'L', 'N', M, M, N, A, LDA,
1879:      $                      WORK( ITAUQ ), U, LDU, WORK( NWORK ),
1880:      $                      LWORK-NWORK+1, IERR )
1881: *
1882: *              Copy real matrix RWORK(IRVT) to complex matrix VT
1883: *              Overwrite VT by right singular vectors of A
1884: *              (CWorkspace: need 3*M, prefer 2*M+M*NB)
1885: *              (RWorkspace: M*M)
1886: *
1887:                CALL CLASET( 'F', M, N, CZERO, CZERO, VT, LDVT )
1888:                CALL CLACP2( 'F', M, M, RWORK( IRVT ), M, VT, LDVT )
1889:                CALL CUNMBR( 'P', 'R', 'C', M, N, M, A, LDA,
1890:      $                      WORK( ITAUP ), VT, LDVT, WORK( NWORK ),
1891:      $                      LWORK-NWORK+1, IERR )
1892:             ELSE
1893: *
1894: *              Perform bidiagonal SVD, computing left singular vectors
1895: *              of bidiagonal matrix in RWORK(IRU) and computing right
1896: *              singular vectors of bidiagonal matrix in RWORK(IRVT)
1897: *              (CWorkspace: need 0)
1898: *              (RWorkspace: need BDSPAC)
1899: *
1900:                IRVT = NRWORK
1901:                IRU = IRVT + M*M
1902:                NRWORK = IRU + M*M
1903: *
1904:                CALL SBDSDC( 'L', 'I', M, S, RWORK( IE ), RWORK( IRU ),
1905:      $                      M, RWORK( IRVT ), M, DUM, IDUM,
1906:      $                      RWORK( NRWORK ), IWORK, INFO )
1907: *
1908: *              Copy real matrix RWORK(IRU) to complex matrix U
1909: *              Overwrite U by left singular vectors of A
1910: *              (CWorkspace: need 3*M, prefer 2*M+M*NB)
1911: *              (RWorkspace: M*M)
1912: *
1913:                CALL CLACP2( 'F', M, M, RWORK( IRU ), M, U, LDU )
1914:                CALL CUNMBR( 'Q', 'L', 'N', M, M, N, A, LDA,
1915:      $                      WORK( ITAUQ ), U, LDU, WORK( NWORK ),
1916:      $                      LWORK-NWORK+1, IERR )
1917: *
1918: *              Set all of VT to identity matrix
1919: *
1920:                CALL CLASET( 'F', N, N, CZERO, CONE, VT, LDVT )
1921: *
1922: *              Copy real matrix RWORK(IRVT) to complex matrix VT
1923: *              Overwrite VT by right singular vectors of A
1924: *              (CWorkspace: need 2*M+N, prefer 2*M+N*NB)
1925: *              (RWorkspace: M*M)
1926: *
1927:                CALL CLACP2( 'F', M, M, RWORK( IRVT ), M, VT, LDVT )
1928:                CALL CUNMBR( 'P', 'R', 'C', N, N, M, A, LDA,
1929:      $                      WORK( ITAUP ), VT, LDVT, WORK( NWORK ),
1930:      $                      LWORK-NWORK+1, IERR )
1931:             END IF
1932: *
1933:          END IF
1934: *
1935:       END IF
1936: *
1937: *     Undo scaling if necessary
1938: *
1939:       IF( ISCL.EQ.1 ) THEN
1940:          IF( ANRM.GT.BIGNUM )
1941:      $      CALL SLASCL( 'G', 0, 0, BIGNUM, ANRM, MINMN, 1, S, MINMN,
1942:      $                   IERR )
1943:          IF( INFO.NE.0 .AND. ANRM.GT.BIGNUM )
1944:      $      CALL SLASCL( 'G', 0, 0, BIGNUM, ANRM, MINMN-1, 1,
1945:      $                   RWORK( IE ), MINMN, IERR )
1946:          IF( ANRM.LT.SMLNUM )
1947:      $      CALL SLASCL( 'G', 0, 0, SMLNUM, ANRM, MINMN, 1, S, MINMN,
1948:      $                   IERR )
1949:          IF( INFO.NE.0 .AND. ANRM.LT.SMLNUM )
1950:      $      CALL SLASCL( 'G', 0, 0, SMLNUM, ANRM, MINMN-1, 1,
1951:      $                   RWORK( IE ), MINMN, IERR )
1952:       END IF
1953: *
1954: *     Return optimal workspace in WORK(1)
1955: *
1956:       WORK( 1 ) = MAXWRK
1957: *
1958:       RETURN
1959: *
1960: *     End of CGESDD
1961: *
1962:       END
1963: