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