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