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