0001:       SUBROUTINE SGESDD( JOBZ, M, N, A, LDA, S, U, LDU, VT, LDVT, WORK,
0002:      $                   LWORK, IWORK, INFO )
0003: *
0004: *  -- LAPACK driver routine (version 3.2.1)                                  --
0005: *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
0006: *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
0007: *     March 2009
0008: *
0009: *     .. Scalar Arguments ..
0010:       CHARACTER          JOBZ
0011:       INTEGER            INFO, LDA, LDU, LDVT, LWORK, M, N
0012: *     ..
0013: *     .. Array Arguments ..
0014:       INTEGER            IWORK( * )
0015:       REAL               A( LDA, * ), S( * ), U( LDU, * ),
0016:      $                   VT( LDVT, * ), WORK( * )
0017: *     ..
0018: *
0019: *  Purpose
0020: *  =======
0021: *
0022: *  SGESDD computes the singular value decomposition (SVD) of a real
0023: *  M-by-N matrix A, optionally computing the left and right singular
0024: *  vectors.  If singular vectors are desired, it uses a
0025: *  divide-and-conquer algorithm.
0026: *
0027: *  The SVD is written
0028: *
0029: *       A = U * SIGMA * transpose(V)
0030: *
0031: *  where SIGMA is an M-by-N matrix which is zero except for its
0032: *  min(m,n) diagonal elements, U is an M-by-M orthogonal matrix, and
0033: *  V is an N-by-N orthogonal matrix.  The diagonal elements of SIGMA
0034: *  are the singular values of A; they are real and non-negative, and
0035: *  are returned in descending order.  The first min(m,n) columns of
0036: *  U and V are the left and right singular vectors of A.
0037: *
0038: *  Note that the routine returns VT = V**T, not V.
0039: *
0040: *  The divide and conquer algorithm makes very mild assumptions about
0041: *  floating point arithmetic. It will work on machines with a guard
0042: *  digit in add/subtract, or on those binary machines without guard
0043: *  digits which subtract like the Cray X-MP, Cray Y-MP, Cray C-90, or
0044: *  Cray-2. It could conceivably fail on hexadecimal or decimal machines
0045: *  without guard digits, but we know of none.
0046: *
0047: *  Arguments
0048: *  =========
0049: *
0050: *  JOBZ    (input) CHARACTER*1
0051: *          Specifies options for computing all or part of the matrix U:
0052: *          = 'A':  all M columns of U and all N rows of V**T are
0053: *                  returned in the arrays U and VT;
0054: *          = 'S':  the first min(M,N) columns of U and the first
0055: *                  min(M,N) rows of V**T are returned in the arrays U
0056: *                  and VT;
0057: *          = 'O':  If M >= N, the first N columns of U are overwritten
0058: *                  on the array A and all rows of V**T are returned in
0059: *                  the array VT;
0060: *                  otherwise, all columns of U are returned in the
0061: *                  array U and the first M rows of V**T are overwritten
0062: *                  in the array A;
0063: *          = 'N':  no columns of U or rows of V**T are computed.
0064: *
0065: *  M       (input) INTEGER
0066: *          The number of rows of the input matrix A.  M >= 0.
0067: *
0068: *  N       (input) INTEGER
0069: *          The number of columns of the input matrix A.  N >= 0.
0070: *
0071: *  A       (input/output) REAL array, dimension (LDA,N)
0072: *          On entry, the M-by-N matrix A.
0073: *          On exit,
0074: *          if JOBZ = 'O',  A is overwritten with the first N columns
0075: *                          of U (the left singular vectors, stored
0076: *                          columnwise) if M >= N;
0077: *                          A is overwritten with the first M rows
0078: *                          of V**T (the right singular vectors, stored
0079: *                          rowwise) otherwise.
0080: *          if JOBZ .ne. 'O', the contents of A are destroyed.
0081: *
0082: *  LDA     (input) INTEGER
0083: *          The leading dimension of the array A.  LDA >= max(1,M).
0084: *
0085: *  S       (output) REAL array, dimension (min(M,N))
0086: *          The singular values of A, sorted so that S(i) >= S(i+1).
0087: *
0088: *  U       (output) REAL array, dimension (LDU,UCOL)
0089: *          UCOL = M if JOBZ = 'A' or JOBZ = 'O' and M < N;
0090: *          UCOL = min(M,N) if JOBZ = 'S'.
0091: *          If JOBZ = 'A' or JOBZ = 'O' and M < N, U contains the M-by-M
0092: *          orthogonal matrix U;
0093: *          if JOBZ = 'S', U contains the first min(M,N) columns of U
0094: *          (the left singular vectors, stored columnwise);
0095: *          if JOBZ = 'O' and M >= N, or JOBZ = 'N', U is not referenced.
0096: *
0097: *  LDU     (input) INTEGER
0098: *          The leading dimension of the array U.  LDU >= 1; if
0099: *          JOBZ = 'S' or 'A' or JOBZ = 'O' and M < N, LDU >= M.
0100: *
0101: *  VT      (output) REAL array, dimension (LDVT,N)
0102: *          If JOBZ = 'A' or JOBZ = 'O' and M >= N, VT contains the
0103: *          N-by-N orthogonal matrix V**T;
0104: *          if JOBZ = 'S', VT contains the first min(M,N) rows of
0105: *          V**T (the right singular vectors, stored rowwise);
0106: *          if JOBZ = 'O' and M < N, or JOBZ = 'N', VT is not referenced.
0107: *
0108: *  LDVT    (input) INTEGER
0109: *          The leading dimension of the array VT.  LDVT >= 1; if
0110: *          JOBZ = 'A' or JOBZ = 'O' and M >= N, LDVT >= N;
0111: *          if JOBZ = 'S', LDVT >= min(M,N).
0112: *
0113: *  WORK    (workspace/output) REAL array, dimension (MAX(1,LWORK))
0114: *          On exit, if INFO = 0, WORK(1) returns the optimal LWORK;
0115: *
0116: *  LWORK   (input) INTEGER
0117: *          The dimension of the array WORK. LWORK >= 1.
0118: *          If JOBZ = 'N',
0119: *            LWORK >= 3*min(M,N) + max(max(M,N),6*min(M,N)).
0120: *          If JOBZ = 'O',
0121: *            LWORK >= 3*min(M,N) + 
0122: *                     max(max(M,N),5*min(M,N)*min(M,N)+4*min(M,N)).
0123: *          If JOBZ = 'S' or 'A'
0124: *            LWORK >= 3*min(M,N) +
0125: *                     max(max(M,N),4*min(M,N)*min(M,N)+4*min(M,N)).
0126: *          For good performance, LWORK should generally be larger.
0127: *          If LWORK = -1 but other input arguments are legal, WORK(1)
0128: *          returns the optimal LWORK.
0129: *
0130: *  IWORK   (workspace) INTEGER array, dimension (8*min(M,N))
0131: *
0132: *  INFO    (output) INTEGER
0133: *          = 0:  successful exit.
0134: *          < 0:  if INFO = -i, the i-th argument had an illegal value.
0135: *          > 0:  SBDSDC did not converge, updating process failed.
0136: *
0137: *  Further Details
0138: *  ===============
0139: *
0140: *  Based on contributions by
0141: *     Ming Gu and Huan Ren, Computer Science Division, University of
0142: *     California at Berkeley, USA
0143: *
0144: *  =====================================================================
0145: *
0146: *     .. Parameters ..
0147:       REAL               ZERO, ONE
0148:       PARAMETER          ( ZERO = 0.0E0, ONE = 1.0E0 )
0149: *     ..
0150: *     .. Local Scalars ..
0151:       LOGICAL            LQUERY, WNTQA, WNTQAS, WNTQN, WNTQO, WNTQS
0152:       INTEGER            BDSPAC, BLK, CHUNK, I, IE, IERR, IL,
0153:      $                   IR, ISCL, ITAU, ITAUP, ITAUQ, IU, IVT, LDWKVT,
0154:      $                   LDWRKL, LDWRKR, LDWRKU, MAXWRK, MINMN, MINWRK,
0155:      $                   MNTHR, NWORK, WRKBL
0156:       REAL               ANRM, BIGNUM, EPS, SMLNUM
0157: *     ..
0158: *     .. Local Arrays ..
0159:       INTEGER            IDUM( 1 )
0160:       REAL               DUM( 1 )
0161: *     ..
0162: *     .. External Subroutines ..
0163:       EXTERNAL           SBDSDC, SGEBRD, SGELQF, SGEMM, SGEQRF, SLACPY,
0164:      $                   SLASCL, SLASET, SORGBR, SORGLQ, SORGQR, SORMBR,
0165:      $                   XERBLA
0166: *     ..
0167: *     .. External Functions ..
0168:       LOGICAL            LSAME
0169:       INTEGER            ILAENV
0170:       REAL               SLAMCH, SLANGE
0171:       EXTERNAL           ILAENV, LSAME, SLAMCH, SLANGE
0172: *     ..
0173: *     .. Intrinsic Functions ..
0174:       INTRINSIC          INT, MAX, MIN, SQRT
0175: *     ..
0176: *     .. Executable Statements ..
0177: *
0178: *     Test the input arguments
0179: *
0180:       INFO = 0
0181:       MINMN = MIN( M, N )
0182:       WNTQA = LSAME( JOBZ, 'A' )
0183:       WNTQS = LSAME( JOBZ, 'S' )
0184:       WNTQAS = WNTQA .OR. WNTQS
0185:       WNTQO = LSAME( JOBZ, 'O' )
0186:       WNTQN = LSAME( JOBZ, 'N' )
0187:       LQUERY = ( LWORK.EQ.-1 )
0188: *
0189:       IF( .NOT.( WNTQA .OR. WNTQS .OR. WNTQO .OR. WNTQN ) ) THEN
0190:          INFO = -1
0191:       ELSE IF( M.LT.0 ) THEN
0192:          INFO = -2
0193:       ELSE IF( N.LT.0 ) THEN
0194:          INFO = -3
0195:       ELSE IF( LDA.LT.MAX( 1, M ) ) THEN
0196:          INFO = -5
0197:       ELSE IF( LDU.LT.1 .OR. ( WNTQAS .AND. LDU.LT.M ) .OR.
0198:      $         ( WNTQO .AND. M.LT.N .AND. LDU.LT.M ) ) THEN
0199:          INFO = -8
0200:       ELSE IF( LDVT.LT.1 .OR. ( WNTQA .AND. LDVT.LT.N ) .OR.
0201:      $         ( WNTQS .AND. LDVT.LT.MINMN ) .OR.
0202:      $         ( WNTQO .AND. M.GE.N .AND. LDVT.LT.N ) ) THEN
0203:          INFO = -10
0204:       END IF
0205: *
0206: *     Compute workspace
0207: *      (Note: Comments in the code beginning "Workspace:" describe the
0208: *       minimal amount of workspace needed at that point in the code,
0209: *       as well as the preferred amount for good performance.
0210: *       NB refers to the optimal block size for the immediately
0211: *       following subroutine, as returned by ILAENV.)
0212: *
0213:       IF( INFO.EQ.0 ) THEN
0214:          MINWRK = 1
0215:          MAXWRK = 1
0216:          IF( M.GE.N .AND. MINMN.GT.0 ) THEN
0217: *
0218: *           Compute space needed for SBDSDC
0219: *
0220:             MNTHR = INT( MINMN*11.0E0 / 6.0E0 )
0221:             IF( WNTQN ) THEN
0222:                BDSPAC = 7*N
0223:             ELSE
0224:                BDSPAC = 3*N*N + 4*N
0225:             END IF
0226:             IF( M.GE.MNTHR ) THEN
0227:                IF( WNTQN ) THEN
0228: *
0229: *                 Path 1 (M much larger than N, JOBZ='N')
0230: *
0231:                   WRKBL = N + N*ILAENV( 1, 'SGEQRF', ' ', M, N, -1,
0232:      $                    -1 )
0233:                   WRKBL = MAX( WRKBL, 3*N+2*N*
0234:      $                    ILAENV( 1, 'SGEBRD', ' ', N, N, -1, -1 ) )
0235:                   MAXWRK = MAX( WRKBL, BDSPAC+N )
0236:                   MINWRK = BDSPAC + N
0237:                ELSE IF( WNTQO ) THEN
0238: *
0239: *                 Path 2 (M much larger than N, JOBZ='O')
0240: *
0241:                   WRKBL = N + N*ILAENV( 1, 'SGEQRF', ' ', M, N, -1, -1 )
0242:                   WRKBL = MAX( WRKBL, N+N*ILAENV( 1, 'SORGQR', ' ', M,
0243:      $                    N, N, -1 ) )
0244:                   WRKBL = MAX( WRKBL, 3*N+2*N*
0245:      $                    ILAENV( 1, 'SGEBRD', ' ', N, N, -1, -1 ) )
0246:                   WRKBL = MAX( WRKBL, 3*N+N*
0247:      $                    ILAENV( 1, 'SORMBR', 'QLN', N, N, N, -1 ) )
0248:                   WRKBL = MAX( WRKBL, 3*N+N*
0249:      $                    ILAENV( 1, 'SORMBR', 'PRT', N, N, N, -1 ) )
0250:                   WRKBL = MAX( WRKBL, BDSPAC+3*N )
0251:                   MAXWRK = WRKBL + 2*N*N
0252:                   MINWRK = BDSPAC + 2*N*N + 3*N
0253:                ELSE IF( WNTQS ) THEN
0254: *
0255: *                 Path 3 (M much larger than N, JOBZ='S')
0256: *
0257:                   WRKBL = N + N*ILAENV( 1, 'SGEQRF', ' ', M, N, -1, -1 )
0258:                   WRKBL = MAX( WRKBL, N+N*ILAENV( 1, 'SORGQR', ' ', M,
0259:      $                    N, N, -1 ) )
0260:                   WRKBL = MAX( WRKBL, 3*N+2*N*
0261:      $                    ILAENV( 1, 'SGEBRD', ' ', N, N, -1, -1 ) )
0262:                   WRKBL = MAX( WRKBL, 3*N+N*
0263:      $                    ILAENV( 1, 'SORMBR', 'QLN', N, N, N, -1 ) )
0264:                   WRKBL = MAX( WRKBL, 3*N+N*
0265:      $                    ILAENV( 1, 'SORMBR', 'PRT', N, N, N, -1 ) )
0266:                   WRKBL = MAX( WRKBL, BDSPAC+3*N )
0267:                   MAXWRK = WRKBL + N*N
0268:                   MINWRK = BDSPAC + N*N + 3*N
0269:                ELSE IF( WNTQA ) THEN
0270: *
0271: *                 Path 4 (M much larger than N, JOBZ='A')
0272: *
0273:                   WRKBL = N + N*ILAENV( 1, 'SGEQRF', ' ', M, N, -1, -1 )
0274:                   WRKBL = MAX( WRKBL, N+M*ILAENV( 1, 'SORGQR', ' ', M,
0275:      $                    M, N, -1 ) )
0276:                   WRKBL = MAX( WRKBL, 3*N+2*N*
0277:      $                    ILAENV( 1, 'SGEBRD', ' ', N, N, -1, -1 ) )
0278:                   WRKBL = MAX( WRKBL, 3*N+N*
0279:      $                    ILAENV( 1, 'SORMBR', 'QLN', N, N, N, -1 ) )
0280:                   WRKBL = MAX( WRKBL, 3*N+N*
0281:      $                    ILAENV( 1, 'SORMBR', 'PRT', N, N, N, -1 ) )
0282:                   WRKBL = MAX( WRKBL, BDSPAC+3*N )
0283:                   MAXWRK = WRKBL + N*N
0284:                   MINWRK = BDSPAC + N*N + 3*N
0285:                END IF
0286:             ELSE
0287: *
0288: *              Path 5 (M at least N, but not much larger)
0289: *
0290:                WRKBL = 3*N + ( M+N )*ILAENV( 1, 'SGEBRD', ' ', M, N, -1,
0291:      $                 -1 )
0292:                IF( WNTQN ) THEN
0293:                   MAXWRK = MAX( WRKBL, BDSPAC+3*N )
0294:                   MINWRK = 3*N + MAX( M, BDSPAC )
0295:                ELSE IF( WNTQO ) THEN
0296:                   WRKBL = MAX( WRKBL, 3*N+N*
0297:      $                    ILAENV( 1, 'SORMBR', 'QLN', M, N, N, -1 ) )
0298:                   WRKBL = MAX( WRKBL, 3*N+N*
0299:      $                    ILAENV( 1, 'SORMBR', 'PRT', N, N, N, -1 ) )
0300:                   WRKBL = MAX( WRKBL, BDSPAC+3*N )
0301:                   MAXWRK = WRKBL + M*N
0302:                   MINWRK = 3*N + MAX( M, N*N+BDSPAC )
0303:                ELSE IF( WNTQS ) THEN
0304:                   WRKBL = MAX( WRKBL, 3*N+N*
0305:      $                    ILAENV( 1, 'SORMBR', 'QLN', M, N, N, -1 ) )
0306:                   WRKBL = MAX( WRKBL, 3*N+N*
0307:      $                    ILAENV( 1, 'SORMBR', 'PRT', N, N, N, -1 ) )
0308:                   MAXWRK = MAX( WRKBL, BDSPAC+3*N )
0309:                   MINWRK = 3*N + MAX( M, BDSPAC )
0310:                ELSE IF( WNTQA ) THEN
0311:                   WRKBL = MAX( WRKBL, 3*N+M*
0312:      $                    ILAENV( 1, 'SORMBR', 'QLN', M, M, N, -1 ) )
0313:                   WRKBL = MAX( WRKBL, 3*N+N*
0314:      $                    ILAENV( 1, 'SORMBR', 'PRT', N, N, N, -1 ) )
0315:                   MAXWRK = MAX( MAXWRK, BDSPAC+3*N )
0316:                   MINWRK = 3*N + MAX( M, BDSPAC )
0317:                END IF
0318:             END IF
0319:          ELSE IF ( MINMN.GT.0 ) THEN
0320: *
0321: *           Compute space needed for SBDSDC
0322: *
0323:             MNTHR = INT( MINMN*11.0E0 / 6.0E0 )
0324:             IF( WNTQN ) THEN
0325:                BDSPAC = 7*M
0326:             ELSE
0327:                BDSPAC = 3*M*M + 4*M
0328:             END IF
0329:             IF( N.GE.MNTHR ) THEN
0330:                IF( WNTQN ) THEN
0331: *
0332: *                 Path 1t (N much larger than M, JOBZ='N')
0333: *
0334:                   WRKBL = M + M*ILAENV( 1, 'SGELQF', ' ', M, N, -1,
0335:      $                    -1 )
0336:                   WRKBL = MAX( WRKBL, 3*M+2*M*
0337:      $                    ILAENV( 1, 'SGEBRD', ' ', M, M, -1, -1 ) )
0338:                   MAXWRK = MAX( WRKBL, BDSPAC+M )
0339:                   MINWRK = BDSPAC + M
0340:                ELSE IF( WNTQO ) THEN
0341: *
0342: *                 Path 2t (N much larger than M, JOBZ='O')
0343: *
0344:                   WRKBL = M + M*ILAENV( 1, 'SGELQF', ' ', M, N, -1, -1 )
0345:                   WRKBL = MAX( WRKBL, M+M*ILAENV( 1, 'SORGLQ', ' ', M,
0346:      $                    N, M, -1 ) )
0347:                   WRKBL = MAX( WRKBL, 3*M+2*M*
0348:      $                    ILAENV( 1, 'SGEBRD', ' ', M, M, -1, -1 ) )
0349:                   WRKBL = MAX( WRKBL, 3*M+M*
0350:      $                    ILAENV( 1, 'SORMBR', 'QLN', M, M, M, -1 ) )
0351:                   WRKBL = MAX( WRKBL, 3*M+M*
0352:      $                    ILAENV( 1, 'SORMBR', 'PRT', M, M, M, -1 ) )
0353:                   WRKBL = MAX( WRKBL, BDSPAC+3*M )
0354:                   MAXWRK = WRKBL + 2*M*M
0355:                   MINWRK = BDSPAC + 2*M*M + 3*M
0356:                ELSE IF( WNTQS ) THEN
0357: *
0358: *                 Path 3t (N much larger than M, JOBZ='S')
0359: *
0360:                   WRKBL = M + M*ILAENV( 1, 'SGELQF', ' ', M, N, -1, -1 )
0361:                   WRKBL = MAX( WRKBL, M+M*ILAENV( 1, 'SORGLQ', ' ', M,
0362:      $                    N, M, -1 ) )
0363:                   WRKBL = MAX( WRKBL, 3*M+2*M*
0364:      $                    ILAENV( 1, 'SGEBRD', ' ', M, M, -1, -1 ) )
0365:                   WRKBL = MAX( WRKBL, 3*M+M*
0366:      $                    ILAENV( 1, 'SORMBR', 'QLN', M, M, M, -1 ) )
0367:                   WRKBL = MAX( WRKBL, 3*M+M*
0368:      $                    ILAENV( 1, 'SORMBR', 'PRT', M, M, M, -1 ) )
0369:                   WRKBL = MAX( WRKBL, BDSPAC+3*M )
0370:                   MAXWRK = WRKBL + M*M
0371:                   MINWRK = BDSPAC + M*M + 3*M
0372:                ELSE IF( WNTQA ) THEN
0373: *
0374: *                 Path 4t (N much larger than M, JOBZ='A')
0375: *
0376:                   WRKBL = M + M*ILAENV( 1, 'SGELQF', ' ', M, N, -1, -1 )
0377:                   WRKBL = MAX( WRKBL, M+N*ILAENV( 1, 'SORGLQ', ' ', N,
0378:      $                    N, M, -1 ) )
0379:                   WRKBL = MAX( WRKBL, 3*M+2*M*
0380:      $                    ILAENV( 1, 'SGEBRD', ' ', M, M, -1, -1 ) )
0381:                   WRKBL = MAX( WRKBL, 3*M+M*
0382:      $                    ILAENV( 1, 'SORMBR', 'QLN', M, M, M, -1 ) )
0383:                   WRKBL = MAX( WRKBL, 3*M+M*
0384:      $                    ILAENV( 1, 'SORMBR', 'PRT', M, M, M, -1 ) )
0385:                   WRKBL = MAX( WRKBL, BDSPAC+3*M )
0386:                   MAXWRK = WRKBL + M*M
0387:                   MINWRK = BDSPAC + M*M + 3*M
0388:                END IF
0389:             ELSE
0390: *
0391: *              Path 5t (N greater than M, but not much larger)
0392: *
0393:                WRKBL = 3*M + ( M+N )*ILAENV( 1, 'SGEBRD', ' ', M, N, -1,
0394:      $                 -1 )
0395:                IF( WNTQN ) THEN
0396:                   MAXWRK = MAX( WRKBL, BDSPAC+3*M )
0397:                   MINWRK = 3*M + MAX( N, BDSPAC )
0398:                ELSE IF( WNTQO ) THEN
0399:                   WRKBL = MAX( WRKBL, 3*M+M*
0400:      $                    ILAENV( 1, 'SORMBR', 'QLN', M, M, N, -1 ) )
0401:                   WRKBL = MAX( WRKBL, 3*M+M*
0402:      $                    ILAENV( 1, 'SORMBR', 'PRT', M, N, M, -1 ) )
0403:                   WRKBL = MAX( WRKBL, BDSPAC+3*M )
0404:                   MAXWRK = WRKBL + M*N
0405:                   MINWRK = 3*M + MAX( N, M*M+BDSPAC )
0406:                ELSE IF( WNTQS ) THEN
0407:                   WRKBL = MAX( WRKBL, 3*M+M*
0408:      $                    ILAENV( 1, 'SORMBR', 'QLN', M, M, N, -1 ) )
0409:                   WRKBL = MAX( WRKBL, 3*M+M*
0410:      $                    ILAENV( 1, 'SORMBR', 'PRT', M, N, M, -1 ) )
0411:                   MAXWRK = MAX( WRKBL, BDSPAC+3*M )
0412:                   MINWRK = 3*M + MAX( N, BDSPAC )
0413:                ELSE IF( WNTQA ) THEN
0414:                   WRKBL = MAX( WRKBL, 3*M+M*
0415:      $                    ILAENV( 1, 'SORMBR', 'QLN', M, M, N, -1 ) )
0416:                   WRKBL = MAX( WRKBL, 3*M+M*
0417:      $                    ILAENV( 1, 'SORMBR', 'PRT', N, N, M, -1 ) )
0418:                   MAXWRK = MAX( WRKBL, BDSPAC+3*M )
0419:                   MINWRK = 3*M + MAX( N, BDSPAC )
0420:                END IF
0421:             END IF
0422:          END IF
0423:          MAXWRK = MAX( MAXWRK, MINWRK )
0424:          WORK( 1 ) = MAXWRK
0425: *
0426:          IF( LWORK.LT.MINWRK .AND. .NOT.LQUERY ) THEN
0427:             INFO = -12
0428:          END IF
0429:       END IF
0430: *
0431:       IF( INFO.NE.0 ) THEN
0432:          CALL XERBLA( 'SGESDD', -INFO )
0433:          RETURN
0434:       ELSE IF( LQUERY ) THEN
0435:          RETURN
0436:       END IF
0437: *
0438: *     Quick return if possible
0439: *
0440:       IF( M.EQ.0 .OR. N.EQ.0 ) THEN
0441:          RETURN
0442:       END IF
0443: *
0444: *     Get machine constants
0445: *
0446:       EPS = SLAMCH( 'P' )
0447:       SMLNUM = SQRT( SLAMCH( 'S' ) ) / EPS
0448:       BIGNUM = ONE / SMLNUM
0449: *
0450: *     Scale A if max element outside range [SMLNUM,BIGNUM]
0451: *
0452:       ANRM = SLANGE( 'M', M, N, A, LDA, DUM )
0453:       ISCL = 0
0454:       IF( ANRM.GT.ZERO .AND. ANRM.LT.SMLNUM ) THEN
0455:          ISCL = 1
0456:          CALL SLASCL( 'G', 0, 0, ANRM, SMLNUM, M, N, A, LDA, IERR )
0457:       ELSE IF( ANRM.GT.BIGNUM ) THEN
0458:          ISCL = 1
0459:          CALL SLASCL( 'G', 0, 0, ANRM, BIGNUM, M, N, A, LDA, IERR )
0460:       END IF
0461: *
0462:       IF( M.GE.N ) THEN
0463: *
0464: *        A has at least as many rows as columns. If A has sufficiently
0465: *        more rows than columns, first reduce using the QR
0466: *        decomposition (if sufficient workspace available)
0467: *
0468:          IF( M.GE.MNTHR ) THEN
0469: *
0470:             IF( WNTQN ) THEN
0471: *
0472: *              Path 1 (M much larger than N, JOBZ='N')
0473: *              No singular vectors to be computed
0474: *
0475:                ITAU = 1
0476:                NWORK = ITAU + N
0477: *
0478: *              Compute A=Q*R
0479: *              (Workspace: need 2*N, prefer N+N*NB)
0480: *
0481:                CALL SGEQRF( M, N, A, LDA, WORK( ITAU ), WORK( NWORK ),
0482:      $                      LWORK-NWORK+1, IERR )
0483: *
0484: *              Zero out below R
0485: *
0486:                CALL SLASET( 'L', N-1, N-1, ZERO, ZERO, A( 2, 1 ), LDA )
0487:                IE = 1
0488:                ITAUQ = IE + N
0489:                ITAUP = ITAUQ + N
0490:                NWORK = ITAUP + N
0491: *
0492: *              Bidiagonalize R in A
0493: *              (Workspace: need 4*N, prefer 3*N+2*N*NB)
0494: *
0495:                CALL SGEBRD( N, N, A, LDA, S, WORK( IE ), WORK( ITAUQ ),
0496:      $                      WORK( ITAUP ), WORK( NWORK ), LWORK-NWORK+1,
0497:      $                      IERR )
0498:                NWORK = IE + N
0499: *
0500: *              Perform bidiagonal SVD, computing singular values only
0501: *              (Workspace: need N+BDSPAC)
0502: *
0503:                CALL SBDSDC( 'U', 'N', N, S, WORK( IE ), DUM, 1, DUM, 1,
0504:      $                      DUM, IDUM, WORK( NWORK ), IWORK, INFO )
0505: *
0506:             ELSE IF( WNTQO ) THEN
0507: *
0508: *              Path 2 (M much larger than N, JOBZ = 'O')
0509: *              N left singular vectors to be overwritten on A and
0510: *              N right singular vectors to be computed in VT
0511: *
0512:                IR = 1
0513: *
0514: *              WORK(IR) is LDWRKR by N
0515: *
0516:                IF( LWORK.GE.LDA*N+N*N+3*N+BDSPAC ) THEN
0517:                   LDWRKR = LDA
0518:                ELSE
0519:                   LDWRKR = ( LWORK-N*N-3*N-BDSPAC ) / N
0520:                END IF
0521:                ITAU = IR + LDWRKR*N
0522:                NWORK = ITAU + N
0523: *
0524: *              Compute A=Q*R
0525: *              (Workspace: need N*N+2*N, prefer N*N+N+N*NB)
0526: *
0527:                CALL SGEQRF( M, N, A, LDA, WORK( ITAU ), WORK( NWORK ),
0528:      $                      LWORK-NWORK+1, IERR )
0529: *
0530: *              Copy R to WORK(IR), zeroing out below it
0531: *
0532:                CALL SLACPY( 'U', N, N, A, LDA, WORK( IR ), LDWRKR )
0533:                CALL SLASET( 'L', N-1, N-1, ZERO, ZERO, WORK( IR+1 ),
0534:      $                      LDWRKR )
0535: *
0536: *              Generate Q in A
0537: *              (Workspace: need N*N+2*N, prefer N*N+N+N*NB)
0538: *
0539:                CALL SORGQR( M, N, N, A, LDA, WORK( ITAU ),
0540:      $                      WORK( NWORK ), LWORK-NWORK+1, IERR )
0541:                IE = ITAU
0542:                ITAUQ = IE + N
0543:                ITAUP = ITAUQ + N
0544:                NWORK = ITAUP + N
0545: *
0546: *              Bidiagonalize R in VT, copying result to WORK(IR)
0547: *              (Workspace: need N*N+4*N, prefer N*N+3*N+2*N*NB)
0548: *
0549:                CALL SGEBRD( N, N, WORK( IR ), LDWRKR, S, WORK( IE ),
0550:      $                      WORK( ITAUQ ), WORK( ITAUP ), WORK( NWORK ),
0551:      $                      LWORK-NWORK+1, IERR )
0552: *
0553: *              WORK(IU) is N by N
0554: *
0555:                IU = NWORK
0556:                NWORK = IU + N*N
0557: *
0558: *              Perform bidiagonal SVD, computing left singular vectors
0559: *              of bidiagonal matrix in WORK(IU) and computing right
0560: *              singular vectors of bidiagonal matrix in VT
0561: *              (Workspace: need N+N*N+BDSPAC)
0562: *
0563:                CALL SBDSDC( 'U', 'I', N, S, WORK( IE ), WORK( IU ), N,
0564:      $                      VT, LDVT, DUM, IDUM, WORK( NWORK ), IWORK,
0565:      $                      INFO )
0566: *
0567: *              Overwrite WORK(IU) by left singular vectors of R
0568: *              and VT by right singular vectors of R
0569: *              (Workspace: need 2*N*N+3*N, prefer 2*N*N+2*N+N*NB)
0570: *
0571:                CALL SORMBR( 'Q', 'L', 'N', N, N, N, WORK( IR ), LDWRKR,
0572:      $                      WORK( ITAUQ ), WORK( IU ), N, WORK( NWORK ),
0573:      $                      LWORK-NWORK+1, IERR )
0574:                CALL SORMBR( 'P', 'R', 'T', N, N, N, WORK( IR ), LDWRKR,
0575:      $                      WORK( ITAUP ), VT, LDVT, WORK( NWORK ),
0576:      $                      LWORK-NWORK+1, IERR )
0577: *
0578: *              Multiply Q in A by left singular vectors of R in
0579: *              WORK(IU), storing result in WORK(IR) and copying to A
0580: *              (Workspace: need 2*N*N, prefer N*N+M*N)
0581: *
0582:                DO 10 I = 1, M, LDWRKR
0583:                   CHUNK = MIN( M-I+1, LDWRKR )
0584:                   CALL SGEMM( 'N', 'N', CHUNK, N, N, ONE, A( I, 1 ),
0585:      $                        LDA, WORK( IU ), N, ZERO, WORK( IR ),
0586:      $                        LDWRKR )
0587:                   CALL SLACPY( 'F', CHUNK, N, WORK( IR ), LDWRKR,
0588:      $                         A( I, 1 ), LDA )
0589:    10          CONTINUE
0590: *
0591:             ELSE IF( WNTQS ) THEN
0592: *
0593: *              Path 3 (M much larger than N, JOBZ='S')
0594: *              N left singular vectors to be computed in U and
0595: *              N right singular vectors to be computed in VT
0596: *
0597:                IR = 1
0598: *
0599: *              WORK(IR) is N by N
0600: *
0601:                LDWRKR = N
0602:                ITAU = IR + LDWRKR*N
0603:                NWORK = ITAU + N
0604: *
0605: *              Compute A=Q*R
0606: *              (Workspace: need N*N+2*N, prefer N*N+N+N*NB)
0607: *
0608:                CALL SGEQRF( M, N, A, LDA, WORK( ITAU ), WORK( NWORK ),
0609:      $                      LWORK-NWORK+1, IERR )
0610: *
0611: *              Copy R to WORK(IR), zeroing out below it
0612: *
0613:                CALL SLACPY( 'U', N, N, A, LDA, WORK( IR ), LDWRKR )
0614:                CALL SLASET( 'L', N-1, N-1, ZERO, ZERO, WORK( IR+1 ),
0615:      $                      LDWRKR )
0616: *
0617: *              Generate Q in A
0618: *              (Workspace: need N*N+2*N, prefer N*N+N+N*NB)
0619: *
0620:                CALL SORGQR( M, N, N, A, LDA, WORK( ITAU ),
0621:      $                      WORK( NWORK ), LWORK-NWORK+1, IERR )
0622:                IE = ITAU
0623:                ITAUQ = IE + N
0624:                ITAUP = ITAUQ + N
0625:                NWORK = ITAUP + N
0626: *
0627: *              Bidiagonalize R in WORK(IR)
0628: *              (Workspace: need N*N+4*N, prefer N*N+3*N+2*N*NB)
0629: *
0630:                CALL SGEBRD( N, N, WORK( IR ), LDWRKR, S, WORK( IE ),
0631:      $                      WORK( ITAUQ ), WORK( ITAUP ), WORK( NWORK ),
0632:      $                      LWORK-NWORK+1, IERR )
0633: *
0634: *              Perform bidiagonal SVD, computing left singular vectors
0635: *              of bidiagoal matrix in U and computing right singular
0636: *              vectors of bidiagonal matrix in VT
0637: *              (Workspace: need N+BDSPAC)
0638: *
0639:                CALL SBDSDC( 'U', 'I', N, S, WORK( IE ), U, LDU, VT,
0640:      $                      LDVT, DUM, IDUM, WORK( NWORK ), IWORK,
0641:      $                      INFO )
0642: *
0643: *              Overwrite U by left singular vectors of R and VT
0644: *              by right singular vectors of R
0645: *              (Workspace: need N*N+3*N, prefer N*N+2*N+N*NB)
0646: *
0647:                CALL SORMBR( 'Q', 'L', 'N', N, N, N, WORK( IR ), LDWRKR,
0648:      $                      WORK( ITAUQ ), U, LDU, WORK( NWORK ),
0649:      $                      LWORK-NWORK+1, IERR )
0650: *
0651:                CALL SORMBR( 'P', 'R', 'T', N, N, N, WORK( IR ), LDWRKR,
0652:      $                      WORK( ITAUP ), VT, LDVT, WORK( NWORK ),
0653:      $                      LWORK-NWORK+1, IERR )
0654: *
0655: *              Multiply Q in A by left singular vectors of R in
0656: *              WORK(IR), storing result in U
0657: *              (Workspace: need N*N)
0658: *
0659:                CALL SLACPY( 'F', N, N, U, LDU, WORK( IR ), LDWRKR )
0660:                CALL SGEMM( 'N', 'N', M, N, N, ONE, A, LDA, WORK( IR ),
0661:      $                     LDWRKR, ZERO, U, LDU )
0662: *
0663:             ELSE IF( WNTQA ) THEN
0664: *
0665: *              Path 4 (M much larger than N, JOBZ='A')
0666: *              M left singular vectors to be computed in U and
0667: *              N right singular vectors to be computed in VT
0668: *
0669:                IU = 1
0670: *
0671: *              WORK(IU) is N by N
0672: *
0673:                LDWRKU = N
0674:                ITAU = IU + LDWRKU*N
0675:                NWORK = ITAU + N
0676: *
0677: *              Compute A=Q*R, copying result to U
0678: *              (Workspace: need N*N+2*N, prefer N*N+N+N*NB)
0679: *
0680:                CALL SGEQRF( M, N, A, LDA, WORK( ITAU ), WORK( NWORK ),
0681:      $                      LWORK-NWORK+1, IERR )
0682:                CALL SLACPY( 'L', M, N, A, LDA, U, LDU )
0683: *
0684: *              Generate Q in U
0685: *              (Workspace: need N*N+2*N, prefer N*N+N+N*NB)
0686:                CALL SORGQR( M, M, N, U, LDU, WORK( ITAU ),
0687:      $                      WORK( NWORK ), LWORK-NWORK+1, IERR )
0688: *
0689: *              Produce R in A, zeroing out other entries
0690: *
0691:                CALL SLASET( 'L', N-1, N-1, ZERO, ZERO, A( 2, 1 ), LDA )
0692:                IE = ITAU
0693:                ITAUQ = IE + N
0694:                ITAUP = ITAUQ + N
0695:                NWORK = ITAUP + N
0696: *
0697: *              Bidiagonalize R in A
0698: *              (Workspace: need N*N+4*N, prefer N*N+3*N+2*N*NB)
0699: *
0700:                CALL SGEBRD( N, N, A, LDA, S, WORK( IE ), WORK( ITAUQ ),
0701:      $                      WORK( ITAUP ), WORK( NWORK ), LWORK-NWORK+1,
0702:      $                      IERR )
0703: *
0704: *              Perform bidiagonal SVD, computing left singular vectors
0705: *              of bidiagonal matrix in WORK(IU) and computing right
0706: *              singular vectors of bidiagonal matrix in VT
0707: *              (Workspace: need N+N*N+BDSPAC)
0708: *
0709:                CALL SBDSDC( 'U', 'I', N, S, WORK( IE ), WORK( IU ), N,
0710:      $                      VT, LDVT, DUM, IDUM, WORK( NWORK ), IWORK,
0711:      $                      INFO )
0712: *
0713: *              Overwrite WORK(IU) by left singular vectors of R and VT
0714: *              by right singular vectors of R
0715: *              (Workspace: need N*N+3*N, prefer N*N+2*N+N*NB)
0716: *
0717:                CALL SORMBR( 'Q', 'L', 'N', N, N, N, A, LDA,
0718:      $                      WORK( ITAUQ ), WORK( IU ), LDWRKU,
0719:      $                      WORK( NWORK ), LWORK-NWORK+1, IERR )
0720:                CALL SORMBR( 'P', 'R', 'T', N, N, N, A, LDA,
0721:      $                      WORK( ITAUP ), VT, LDVT, WORK( NWORK ),
0722:      $                      LWORK-NWORK+1, IERR )
0723: *
0724: *              Multiply Q in U by left singular vectors of R in
0725: *              WORK(IU), storing result in A
0726: *              (Workspace: need N*N)
0727: *
0728:                CALL SGEMM( 'N', 'N', M, N, N, ONE, U, LDU, WORK( IU ),
0729:      $                     LDWRKU, ZERO, A, LDA )
0730: *
0731: *              Copy left singular vectors of A from A to U
0732: *
0733:                CALL SLACPY( 'F', M, N, A, LDA, U, LDU )
0734: *
0735:             END IF
0736: *
0737:          ELSE
0738: *
0739: *           M .LT. MNTHR
0740: *
0741: *           Path 5 (M at least N, but not much larger)
0742: *           Reduce to bidiagonal form without QR decomposition
0743: *
0744:             IE = 1
0745:             ITAUQ = IE + N
0746:             ITAUP = ITAUQ + N
0747:             NWORK = ITAUP + N
0748: *
0749: *           Bidiagonalize A
0750: *           (Workspace: need 3*N+M, prefer 3*N+(M+N)*NB)
0751: *
0752:             CALL SGEBRD( M, N, A, LDA, S, WORK( IE ), WORK( ITAUQ ),
0753:      $                   WORK( ITAUP ), WORK( NWORK ), LWORK-NWORK+1,
0754:      $                   IERR )
0755:             IF( WNTQN ) THEN
0756: *
0757: *              Perform bidiagonal SVD, only computing singular values
0758: *              (Workspace: need N+BDSPAC)
0759: *
0760:                CALL SBDSDC( 'U', 'N', N, S, WORK( IE ), DUM, 1, DUM, 1,
0761:      $                      DUM, IDUM, WORK( NWORK ), IWORK, INFO )
0762:             ELSE IF( WNTQO ) THEN
0763:                IU = NWORK
0764:                IF( LWORK.GE.M*N+3*N+BDSPAC ) THEN
0765: *
0766: *                 WORK( IU ) is M by N
0767: *
0768:                   LDWRKU = M
0769:                   NWORK = IU + LDWRKU*N
0770:                   CALL SLASET( 'F', M, N, ZERO, ZERO, WORK( IU ),
0771:      $                         LDWRKU )
0772:                ELSE
0773: *
0774: *                 WORK( IU ) is N by N
0775: *
0776:                   LDWRKU = N
0777:                   NWORK = IU + LDWRKU*N
0778: *
0779: *                 WORK(IR) is LDWRKR by N
0780: *
0781:                   IR = NWORK
0782:                   LDWRKR = ( LWORK-N*N-3*N ) / N
0783:                END IF
0784:                NWORK = IU + LDWRKU*N
0785: *
0786: *              Perform bidiagonal SVD, computing left singular vectors
0787: *              of bidiagonal matrix in WORK(IU) and computing right
0788: *              singular vectors of bidiagonal matrix in VT
0789: *              (Workspace: need N+N*N+BDSPAC)
0790: *
0791:                CALL SBDSDC( 'U', 'I', N, S, WORK( IE ), WORK( IU ),
0792:      $                      LDWRKU, VT, LDVT, DUM, IDUM, WORK( NWORK ),
0793:      $                      IWORK, INFO )
0794: *
0795: *              Overwrite VT by right singular vectors of A
0796: *              (Workspace: need N*N+2*N, prefer N*N+N+N*NB)
0797: *
0798:                CALL SORMBR( 'P', 'R', 'T', N, N, N, A, LDA,
0799:      $                      WORK( ITAUP ), VT, LDVT, WORK( NWORK ),
0800:      $                      LWORK-NWORK+1, IERR )
0801: *
0802:                IF( LWORK.GE.M*N+3*N+BDSPAC ) THEN
0803: *
0804: *                 Overwrite WORK(IU) by left singular vectors of A
0805: *                 (Workspace: need N*N+2*N, prefer N*N+N+N*NB)
0806: *
0807:                   CALL SORMBR( 'Q', 'L', 'N', M, N, N, A, LDA,
0808:      $                         WORK( ITAUQ ), WORK( IU ), LDWRKU,
0809:      $                         WORK( NWORK ), LWORK-NWORK+1, IERR )
0810: *
0811: *                 Copy left singular vectors of A from WORK(IU) to A
0812: *
0813:                   CALL SLACPY( 'F', M, N, WORK( IU ), LDWRKU, A, LDA )
0814:                ELSE
0815: *
0816: *                 Generate Q in A
0817: *                 (Workspace: need N*N+2*N, prefer N*N+N+N*NB)
0818: *
0819:                   CALL SORGBR( 'Q', M, N, N, A, LDA, WORK( ITAUQ ),
0820:      $                         WORK( NWORK ), LWORK-NWORK+1, IERR )
0821: *
0822: *                 Multiply Q in A by left singular vectors of
0823: *                 bidiagonal matrix in WORK(IU), storing result in
0824: *                 WORK(IR) and copying to A
0825: *                 (Workspace: need 2*N*N, prefer N*N+M*N)
0826: *
0827:                   DO 20 I = 1, M, LDWRKR
0828:                      CHUNK = MIN( M-I+1, LDWRKR )
0829:                      CALL SGEMM( 'N', 'N', CHUNK, N, N, ONE, A( I, 1 ),
0830:      $                           LDA, WORK( IU ), LDWRKU, ZERO,
0831:      $                           WORK( IR ), LDWRKR )
0832:                      CALL SLACPY( 'F', CHUNK, N, WORK( IR ), LDWRKR,
0833:      $                            A( I, 1 ), LDA )
0834:    20             CONTINUE
0835:                END IF
0836: *
0837:             ELSE IF( WNTQS ) THEN
0838: *
0839: *              Perform bidiagonal SVD, computing left singular vectors
0840: *              of bidiagonal matrix in U and computing right singular
0841: *              vectors of bidiagonal matrix in VT
0842: *              (Workspace: need N+BDSPAC)
0843: *
0844:                CALL SLASET( 'F', M, N, ZERO, ZERO, U, LDU )
0845:                CALL SBDSDC( 'U', 'I', N, S, WORK( IE ), U, LDU, VT,
0846:      $                      LDVT, DUM, IDUM, WORK( NWORK ), IWORK,
0847:      $                      INFO )
0848: *
0849: *              Overwrite U by left singular vectors of A and VT
0850: *              by right singular vectors of A
0851: *              (Workspace: need 3*N, prefer 2*N+N*NB)
0852: *
0853:                CALL SORMBR( 'Q', 'L', 'N', M, N, N, A, LDA,
0854:      $                      WORK( ITAUQ ), U, LDU, WORK( NWORK ),
0855:      $                      LWORK-NWORK+1, IERR )
0856:                CALL SORMBR( 'P', 'R', 'T', N, N, N, A, LDA,
0857:      $                      WORK( ITAUP ), VT, LDVT, WORK( NWORK ),
0858:      $                      LWORK-NWORK+1, IERR )
0859:             ELSE IF( WNTQA ) THEN
0860: *
0861: *              Perform bidiagonal SVD, computing left singular vectors
0862: *              of bidiagonal matrix in U and computing right singular
0863: *              vectors of bidiagonal matrix in VT
0864: *              (Workspace: need N+BDSPAC)
0865: *
0866:                CALL SLASET( 'F', M, M, ZERO, ZERO, U, LDU )
0867:                CALL SBDSDC( 'U', 'I', N, S, WORK( IE ), U, LDU, VT,
0868:      $                      LDVT, DUM, IDUM, WORK( NWORK ), IWORK,
0869:      $                      INFO )
0870: *
0871: *              Set the right corner of U to identity matrix
0872: *
0873:                IF( M.GT.N ) THEN
0874:                   CALL SLASET( 'F', M-N, M-N, ZERO, ONE, U( N+1, N+1 ),
0875:      $                         LDU )
0876:                END IF
0877: *
0878: *              Overwrite U by left singular vectors of A and VT
0879: *              by right singular vectors of A
0880: *              (Workspace: need N*N+2*N+M, prefer N*N+2*N+M*NB)
0881: *
0882:                CALL SORMBR( 'Q', 'L', 'N', M, M, N, A, LDA,
0883:      $                      WORK( ITAUQ ), U, LDU, WORK( NWORK ),
0884:      $                      LWORK-NWORK+1, IERR )
0885:                CALL SORMBR( 'P', 'R', 'T', N, N, M, A, LDA,
0886:      $                      WORK( ITAUP ), VT, LDVT, WORK( NWORK ),
0887:      $                      LWORK-NWORK+1, IERR )
0888:             END IF
0889: *
0890:          END IF
0891: *
0892:       ELSE
0893: *
0894: *        A has more columns than rows. If A has sufficiently more
0895: *        columns than rows, first reduce using the LQ decomposition (if
0896: *        sufficient workspace available)
0897: *
0898:          IF( N.GE.MNTHR ) THEN
0899: *
0900:             IF( WNTQN ) THEN
0901: *
0902: *              Path 1t (N much larger than M, JOBZ='N')
0903: *              No singular vectors to be computed
0904: *
0905:                ITAU = 1
0906:                NWORK = ITAU + M
0907: *
0908: *              Compute A=L*Q
0909: *              (Workspace: need 2*M, prefer M+M*NB)
0910: *
0911:                CALL SGELQF( M, N, A, LDA, WORK( ITAU ), WORK( NWORK ),
0912:      $                      LWORK-NWORK+1, IERR )
0913: *
0914: *              Zero out above L
0915: *
0916:                CALL SLASET( 'U', M-1, M-1, ZERO, ZERO, A( 1, 2 ), LDA )
0917:                IE = 1
0918:                ITAUQ = IE + M
0919:                ITAUP = ITAUQ + M
0920:                NWORK = ITAUP + M
0921: *
0922: *              Bidiagonalize L in A
0923: *              (Workspace: need 4*M, prefer 3*M+2*M*NB)
0924: *
0925:                CALL SGEBRD( M, M, A, LDA, S, WORK( IE ), WORK( ITAUQ ),
0926:      $                      WORK( ITAUP ), WORK( NWORK ), LWORK-NWORK+1,
0927:      $                      IERR )
0928:                NWORK = IE + M
0929: *
0930: *              Perform bidiagonal SVD, computing singular values only
0931: *              (Workspace: need M+BDSPAC)
0932: *
0933:                CALL SBDSDC( 'U', 'N', M, S, WORK( IE ), DUM, 1, DUM, 1,
0934:      $                      DUM, IDUM, WORK( NWORK ), IWORK, INFO )
0935: *
0936:             ELSE IF( WNTQO ) THEN
0937: *
0938: *              Path 2t (N much larger than M, JOBZ='O')
0939: *              M right singular vectors to be overwritten on A and
0940: *              M left singular vectors to be computed in U
0941: *
0942:                IVT = 1
0943: *
0944: *              IVT is M by M
0945: *
0946:                IL = IVT + M*M
0947:                IF( LWORK.GE.M*N+M*M+3*M+BDSPAC ) THEN
0948: *
0949: *                 WORK(IL) is M by N
0950: *
0951:                   LDWRKL = M
0952:                   CHUNK = N
0953:                ELSE
0954:                   LDWRKL = M
0955:                   CHUNK = ( LWORK-M*M ) / M
0956:                END IF
0957:                ITAU = IL + LDWRKL*M
0958:                NWORK = ITAU + M
0959: *
0960: *              Compute A=L*Q
0961: *              (Workspace: need M*M+2*M, prefer M*M+M+M*NB)
0962: *
0963:                CALL SGELQF( M, N, A, LDA, WORK( ITAU ), WORK( NWORK ),
0964:      $                      LWORK-NWORK+1, IERR )
0965: *
0966: *              Copy L to WORK(IL), zeroing about above it
0967: *
0968:                CALL SLACPY( 'L', M, M, A, LDA, WORK( IL ), LDWRKL )
0969:                CALL SLASET( 'U', M-1, M-1, ZERO, ZERO,
0970:      $                      WORK( IL+LDWRKL ), LDWRKL )
0971: *
0972: *              Generate Q in A
0973: *              (Workspace: need M*M+2*M, prefer M*M+M+M*NB)
0974: *
0975:                CALL SORGLQ( M, N, M, A, LDA, WORK( ITAU ),
0976:      $                      WORK( NWORK ), LWORK-NWORK+1, IERR )
0977:                IE = ITAU
0978:                ITAUQ = IE + M
0979:                ITAUP = ITAUQ + M
0980:                NWORK = ITAUP + M
0981: *
0982: *              Bidiagonalize L in WORK(IL)
0983: *              (Workspace: need M*M+4*M, prefer M*M+3*M+2*M*NB)
0984: *
0985:                CALL SGEBRD( M, M, WORK( IL ), LDWRKL, S, WORK( IE ),
0986:      $                      WORK( ITAUQ ), WORK( ITAUP ), WORK( NWORK ),
0987:      $                      LWORK-NWORK+1, IERR )
0988: *
0989: *              Perform bidiagonal SVD, computing left singular vectors
0990: *              of bidiagonal matrix in U, and computing right singular
0991: *              vectors of bidiagonal matrix in WORK(IVT)
0992: *              (Workspace: need M+M*M+BDSPAC)
0993: *
0994:                CALL SBDSDC( 'U', 'I', M, S, WORK( IE ), U, LDU,
0995:      $                      WORK( IVT ), M, DUM, IDUM, WORK( NWORK ),
0996:      $                      IWORK, INFO )
0997: *
0998: *              Overwrite U by left singular vectors of L and WORK(IVT)
0999: *              by right singular vectors of L
1000: *              (Workspace: need 2*M*M+3*M, prefer 2*M*M+2*M+M*NB)
1001: *
1002:                CALL SORMBR( 'Q', 'L', 'N', M, M, M, WORK( IL ), LDWRKL,
1003:      $                      WORK( ITAUQ ), U, LDU, WORK( NWORK ),
1004:      $                      LWORK-NWORK+1, IERR )
1005:                CALL SORMBR( 'P', 'R', 'T', M, M, M, WORK( IL ), LDWRKL,
1006:      $                      WORK( ITAUP ), WORK( IVT ), M,
1007:      $                      WORK( NWORK ), LWORK-NWORK+1, IERR )
1008: *
1009: *              Multiply right singular vectors of L in WORK(IVT) by Q
1010: *              in A, storing result in WORK(IL) and copying to A
1011: *              (Workspace: need 2*M*M, prefer M*M+M*N)
1012: *
1013:                DO 30 I = 1, N, CHUNK
1014:                   BLK = MIN( N-I+1, CHUNK )
1015:                   CALL SGEMM( 'N', 'N', M, BLK, M, ONE, WORK( IVT ), M,
1016:      $                        A( 1, I ), LDA, ZERO, WORK( IL ), LDWRKL )
1017:                   CALL SLACPY( 'F', M, BLK, WORK( IL ), LDWRKL,
1018:      $                         A( 1, I ), LDA )
1019:    30          CONTINUE
1020: *
1021:             ELSE IF( WNTQS ) THEN
1022: *
1023: *              Path 3t (N much larger than M, JOBZ='S')
1024: *              M right singular vectors to be computed in VT and
1025: *              M left singular vectors to be computed in U
1026: *
1027:                IL = 1
1028: *
1029: *              WORK(IL) is M by M
1030: *
1031:                LDWRKL = M
1032:                ITAU = IL + LDWRKL*M
1033:                NWORK = ITAU + M
1034: *
1035: *              Compute A=L*Q
1036: *              (Workspace: need M*M+2*M, prefer M*M+M+M*NB)
1037: *
1038:                CALL SGELQF( M, N, A, LDA, WORK( ITAU ), WORK( NWORK ),
1039:      $                      LWORK-NWORK+1, IERR )
1040: *
1041: *              Copy L to WORK(IL), zeroing out above it
1042: *
1043:                CALL SLACPY( 'L', M, M, A, LDA, WORK( IL ), LDWRKL )
1044:                CALL SLASET( 'U', M-1, M-1, ZERO, ZERO,
1045:      $                      WORK( IL+LDWRKL ), LDWRKL )
1046: *
1047: *              Generate Q in A
1048: *              (Workspace: need M*M+2*M, prefer M*M+M+M*NB)
1049: *
1050:                CALL SORGLQ( M, N, M, A, LDA, WORK( ITAU ),
1051:      $                      WORK( NWORK ), LWORK-NWORK+1, IERR )
1052:                IE = ITAU
1053:                ITAUQ = IE + M
1054:                ITAUP = ITAUQ + M
1055:                NWORK = ITAUP + M
1056: *
1057: *              Bidiagonalize L in WORK(IU), copying result to U
1058: *              (Workspace: need M*M+4*M, prefer M*M+3*M+2*M*NB)
1059: *
1060:                CALL SGEBRD( M, M, WORK( IL ), LDWRKL, S, WORK( IE ),
1061:      $                      WORK( ITAUQ ), WORK( ITAUP ), WORK( NWORK ),
1062:      $                      LWORK-NWORK+1, IERR )
1063: *
1064: *              Perform bidiagonal SVD, computing left singular vectors
1065: *              of bidiagonal matrix in U and computing right singular
1066: *              vectors of bidiagonal matrix in VT
1067: *              (Workspace: need M+BDSPAC)
1068: *
1069:                CALL SBDSDC( 'U', 'I', M, S, WORK( IE ), U, LDU, VT,
1070:      $                      LDVT, DUM, IDUM, WORK( NWORK ), IWORK,
1071:      $                      INFO )
1072: *
1073: *              Overwrite U by left singular vectors of L and VT
1074: *              by right singular vectors of L
1075: *              (Workspace: need M*M+3*M, prefer M*M+2*M+M*NB)
1076: *
1077:                CALL SORMBR( 'Q', 'L', 'N', M, M, M, WORK( IL ), LDWRKL,
1078:      $                      WORK( ITAUQ ), U, LDU, WORK( NWORK ),
1079:      $                      LWORK-NWORK+1, IERR )
1080:                CALL SORMBR( 'P', 'R', 'T', M, M, M, WORK( IL ), LDWRKL,
1081:      $                      WORK( ITAUP ), VT, LDVT, WORK( NWORK ),
1082:      $                      LWORK-NWORK+1, IERR )
1083: *
1084: *              Multiply right singular vectors of L in WORK(IL) by
1085: *              Q in A, storing result in VT
1086: *              (Workspace: need M*M)
1087: *
1088:                CALL SLACPY( 'F', M, M, VT, LDVT, WORK( IL ), LDWRKL )
1089:                CALL SGEMM( 'N', 'N', M, N, M, ONE, WORK( IL ), LDWRKL,
1090:      $                     A, LDA, ZERO, VT, LDVT )
1091: *
1092:             ELSE IF( WNTQA ) THEN
1093: *
1094: *              Path 4t (N much larger than M, JOBZ='A')
1095: *              N right singular vectors to be computed in VT and
1096: *              M left singular vectors to be computed in U
1097: *
1098:                IVT = 1
1099: *
1100: *              WORK(IVT) is M by M
1101: *
1102:                LDWKVT = M
1103:                ITAU = IVT + LDWKVT*M
1104:                NWORK = ITAU + M
1105: *
1106: *              Compute A=L*Q, copying result to VT
1107: *              (Workspace: need M*M+2*M, prefer M*M+M+M*NB)
1108: *
1109:                CALL SGELQF( M, N, A, LDA, WORK( ITAU ), WORK( NWORK ),
1110:      $                      LWORK-NWORK+1, IERR )
1111:                CALL SLACPY( 'U', M, N, A, LDA, VT, LDVT )
1112: *
1113: *              Generate Q in VT
1114: *              (Workspace: need M*M+2*M, prefer M*M+M+M*NB)
1115: *
1116:                CALL SORGLQ( N, N, M, VT, LDVT, WORK( ITAU ),
1117:      $                      WORK( NWORK ), LWORK-NWORK+1, IERR )
1118: *
1119: *              Produce L in A, zeroing out other entries
1120: *
1121:                CALL SLASET( 'U', M-1, M-1, ZERO, ZERO, A( 1, 2 ), LDA )
1122:                IE = ITAU
1123:                ITAUQ = IE + M
1124:                ITAUP = ITAUQ + M
1125:                NWORK = ITAUP + M
1126: *
1127: *              Bidiagonalize L in A
1128: *              (Workspace: need M*M+4*M, prefer M*M+3*M+2*M*NB)
1129: *
1130:                CALL SGEBRD( M, M, A, LDA, S, WORK( IE ), WORK( ITAUQ ),
1131:      $                      WORK( ITAUP ), WORK( NWORK ), LWORK-NWORK+1,
1132:      $                      IERR )
1133: *
1134: *              Perform bidiagonal SVD, computing left singular vectors
1135: *              of bidiagonal matrix in U and computing right singular
1136: *              vectors of bidiagonal matrix in WORK(IVT)
1137: *              (Workspace: need M+M*M+BDSPAC)
1138: *
1139:                CALL SBDSDC( 'U', 'I', M, S, WORK( IE ), U, LDU,
1140:      $                      WORK( IVT ), LDWKVT, DUM, IDUM,
1141:      $                      WORK( NWORK ), IWORK, INFO )
1142: *
1143: *              Overwrite U by left singular vectors of L and WORK(IVT)
1144: *              by right singular vectors of L
1145: *              (Workspace: need M*M+3*M, prefer M*M+2*M+M*NB)
1146: *
1147:                CALL SORMBR( 'Q', 'L', 'N', M, M, M, A, LDA,
1148:      $                      WORK( ITAUQ ), U, LDU, WORK( NWORK ),
1149:      $                      LWORK-NWORK+1, IERR )
1150:                CALL SORMBR( 'P', 'R', 'T', M, M, M, A, LDA,
1151:      $                      WORK( ITAUP ), WORK( IVT ), LDWKVT,
1152:      $                      WORK( NWORK ), LWORK-NWORK+1, IERR )
1153: *
1154: *              Multiply right singular vectors of L in WORK(IVT) by
1155: *              Q in VT, storing result in A
1156: *              (Workspace: need M*M)
1157: *
1158:                CALL SGEMM( 'N', 'N', M, N, M, ONE, WORK( IVT ), LDWKVT,
1159:      $                     VT, LDVT, ZERO, A, LDA )
1160: *
1161: *              Copy right singular vectors of A from A to VT
1162: *
1163:                CALL SLACPY( 'F', M, N, A, LDA, VT, LDVT )
1164: *
1165:             END IF
1166: *
1167:          ELSE
1168: *
1169: *           N .LT. MNTHR
1170: *
1171: *           Path 5t (N greater than M, but not much larger)
1172: *           Reduce to bidiagonal form without LQ decomposition
1173: *
1174:             IE = 1
1175:             ITAUQ = IE + M
1176:             ITAUP = ITAUQ + M
1177:             NWORK = ITAUP + M
1178: *
1179: *           Bidiagonalize A
1180: *           (Workspace: need 3*M+N, prefer 3*M+(M+N)*NB)
1181: *
1182:             CALL SGEBRD( M, N, A, LDA, S, WORK( IE ), WORK( ITAUQ ),
1183:      $                   WORK( ITAUP ), WORK( NWORK ), LWORK-NWORK+1,
1184:      $                   IERR )
1185:             IF( WNTQN ) THEN
1186: *
1187: *              Perform bidiagonal SVD, only computing singular values
1188: *              (Workspace: need M+BDSPAC)
1189: *
1190:                CALL SBDSDC( 'L', 'N', M, S, WORK( IE ), DUM, 1, DUM, 1,
1191:      $                      DUM, IDUM, WORK( NWORK ), IWORK, INFO )
1192:             ELSE IF( WNTQO ) THEN
1193:                LDWKVT = M
1194:                IVT = NWORK
1195:                IF( LWORK.GE.M*N+3*M+BDSPAC ) THEN
1196: *
1197: *                 WORK( IVT ) is M by N
1198: *
1199:                   CALL SLASET( 'F', M, N, ZERO, ZERO, WORK( IVT ),
1200:      $                         LDWKVT )
1201:                   NWORK = IVT + LDWKVT*N
1202:                ELSE
1203: *
1204: *                 WORK( IVT ) is M by M
1205: *
1206:                   NWORK = IVT + LDWKVT*M
1207:                   IL = NWORK
1208: *
1209: *                 WORK(IL) is M by CHUNK
1210: *
1211:                   CHUNK = ( LWORK-M*M-3*M ) / M
1212:                END IF
1213: *
1214: *              Perform bidiagonal SVD, computing left singular vectors
1215: *              of bidiagonal matrix in U and computing right singular
1216: *              vectors of bidiagonal matrix in WORK(IVT)
1217: *              (Workspace: need M*M+BDSPAC)
1218: *
1219:                CALL SBDSDC( 'L', 'I', M, S, WORK( IE ), U, LDU,
1220:      $                      WORK( IVT ), LDWKVT, DUM, IDUM,
1221:      $                      WORK( NWORK ), IWORK, INFO )
1222: *
1223: *              Overwrite U by left singular vectors of A
1224: *              (Workspace: need M*M+2*M, prefer M*M+M+M*NB)
1225: *
1226:                CALL SORMBR( 'Q', 'L', 'N', M, M, N, A, LDA,
1227:      $                      WORK( ITAUQ ), U, LDU, WORK( NWORK ),
1228:      $                      LWORK-NWORK+1, IERR )
1229: *
1230:                IF( LWORK.GE.M*N+3*M+BDSPAC ) THEN
1231: *
1232: *                 Overwrite WORK(IVT) by left singular vectors of A
1233: *                 (Workspace: need M*M+2*M, prefer M*M+M+M*NB)
1234: *
1235:                   CALL SORMBR( 'P', 'R', 'T', M, N, M, A, LDA,
1236:      $                         WORK( ITAUP ), WORK( IVT ), LDWKVT,
1237:      $                         WORK( NWORK ), LWORK-NWORK+1, IERR )
1238: *
1239: *                 Copy right singular vectors of A from WORK(IVT) to A
1240: *
1241:                   CALL SLACPY( 'F', M, N, WORK( IVT ), LDWKVT, A, LDA )
1242:                ELSE
1243: *
1244: *                 Generate P**T in A
1245: *                 (Workspace: need M*M+2*M, prefer M*M+M+M*NB)
1246: *
1247:                   CALL SORGBR( 'P', M, N, M, A, LDA, WORK( ITAUP ),
1248:      $                         WORK( NWORK ), LWORK-NWORK+1, IERR )
1249: *
1250: *                 Multiply Q in A by right singular vectors of
1251: *                 bidiagonal matrix in WORK(IVT), storing result in
1252: *                 WORK(IL) and copying to A
1253: *                 (Workspace: need 2*M*M, prefer M*M+M*N)
1254: *
1255:                   DO 40 I = 1, N, CHUNK
1256:                      BLK = MIN( N-I+1, CHUNK )
1257:                      CALL SGEMM( 'N', 'N', M, BLK, M, ONE, WORK( IVT ),
1258:      $                           LDWKVT, A( 1, I ), LDA, ZERO,
1259:      $                           WORK( IL ), M )
1260:                      CALL SLACPY( 'F', M, BLK, WORK( IL ), M, A( 1, I ),
1261:      $                            LDA )
1262:    40             CONTINUE
1263:                END IF
1264:             ELSE IF( WNTQS ) THEN
1265: *
1266: *              Perform bidiagonal SVD, computing left singular vectors
1267: *              of bidiagonal matrix in U and computing right singular
1268: *              vectors of bidiagonal matrix in VT
1269: *              (Workspace: need M+BDSPAC)
1270: *
1271:                CALL SLASET( 'F', M, N, ZERO, ZERO, VT, LDVT )
1272:                CALL SBDSDC( 'L', 'I', M, S, WORK( IE ), U, LDU, VT,
1273:      $                      LDVT, DUM, IDUM, WORK( NWORK ), IWORK,
1274:      $                      INFO )
1275: *
1276: *              Overwrite U by left singular vectors of A and VT
1277: *              by right singular vectors of A
1278: *              (Workspace: need 3*M, prefer 2*M+M*NB)
1279: *
1280:                CALL SORMBR( 'Q', 'L', 'N', M, M, N, A, LDA,
1281:      $                      WORK( ITAUQ ), U, LDU, WORK( NWORK ),
1282:      $                      LWORK-NWORK+1, IERR )
1283:                CALL SORMBR( 'P', 'R', 'T', M, N, M, A, LDA,
1284:      $                      WORK( ITAUP ), VT, LDVT, WORK( NWORK ),
1285:      $                      LWORK-NWORK+1, IERR )
1286:             ELSE IF( WNTQA ) THEN
1287: *
1288: *              Perform bidiagonal SVD, computing left singular vectors
1289: *              of bidiagonal matrix in U and computing right singular
1290: *              vectors of bidiagonal matrix in VT
1291: *              (Workspace: need M+BDSPAC)
1292: *
1293:                CALL SLASET( 'F', N, N, ZERO, ZERO, VT, LDVT )
1294:                CALL SBDSDC( 'L', 'I', M, S, WORK( IE ), U, LDU, VT,
1295:      $                      LDVT, DUM, IDUM, WORK( NWORK ), IWORK,
1296:      $                      INFO )
1297: *
1298: *              Set the right corner of VT to identity matrix
1299: *
1300:                IF( N.GT.M ) THEN
1301:                   CALL SLASET( 'F', N-M, N-M, ZERO, ONE, VT( M+1, M+1 ),
1302:      $                         LDVT )
1303:                END IF
1304: *
1305: *              Overwrite U by left singular vectors of A and VT
1306: *              by right singular vectors of A
1307: *              (Workspace: need 2*M+N, prefer 2*M+N*NB)
1308: *
1309:                CALL SORMBR( 'Q', 'L', 'N', M, M, N, A, LDA,
1310:      $                      WORK( ITAUQ ), U, LDU, WORK( NWORK ),
1311:      $                      LWORK-NWORK+1, IERR )
1312:                CALL SORMBR( 'P', 'R', 'T', N, N, M, A, LDA,
1313:      $                      WORK( ITAUP ), VT, LDVT, WORK( NWORK ),
1314:      $                      LWORK-NWORK+1, IERR )
1315:             END IF
1316: *
1317:          END IF
1318: *
1319:       END IF
1320: *
1321: *     Undo scaling if necessary
1322: *
1323:       IF( ISCL.EQ.1 ) THEN
1324:          IF( ANRM.GT.BIGNUM )
1325:      $      CALL SLASCL( 'G', 0, 0, BIGNUM, ANRM, MINMN, 1, S, MINMN,
1326:      $                   IERR )
1327:          IF( ANRM.LT.SMLNUM )
1328:      $      CALL SLASCL( 'G', 0, 0, SMLNUM, ANRM, MINMN, 1, S, MINMN,
1329:      $                   IERR )
1330:       END IF
1331: *
1332: *     Return optimal workspace in WORK(1)
1333: *
1334:       WORK( 1 ) = MAXWRK
1335: *
1336:       RETURN
1337: *
1338: *     End of SGESDD
1339: *
1340:       END
1341: