LAPACK: SRC/sgesvd.f Source File

Go to the documentation of this file.
00001       SUBROUTINE SGESVD( JOBU, JOBVT, M, N, A, LDA, S, U, LDU, VT, LDVT,
00002      $                   WORK, LWORK, INFO )
00003 *
00004 *  -- LAPACK driver routine (version 3.2) --
00005 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
00006 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
00007 *     November 2006
00008 *
00009 *     .. Scalar Arguments ..
00010       CHARACTER          JOBU, JOBVT
00011       INTEGER            INFO, LDA, LDU, LDVT, LWORK, M, N
00012 *     ..
00013 *     .. Array Arguments ..
00014       REAL               A( LDA, * ), S( * ), U( LDU, * ),
00015      $                   VT( LDVT, * ), WORK( * )
00016 *     ..
00017 *
00018 *  Purpose
00019 *  =======
00020 *
00021 *  SGESVD computes the singular value decomposition (SVD) of a real
00022 *  M-by-N matrix A, optionally computing the left and/or right singular
00023 *  vectors. The SVD is written
00024 *
00025 *       A = U * SIGMA * transpose(V)
00026 *
00027 *  where SIGMA is an M-by-N matrix which is zero except for its
00028 *  min(m,n) diagonal elements, U is an M-by-M orthogonal matrix, and
00029 *  V is an N-by-N orthogonal matrix.  The diagonal elements of SIGMA
00030 *  are the singular values of A; they are real and non-negative, and
00031 *  are returned in descending order.  The first min(m,n) columns of
00032 *  U and V are the left and right singular vectors of A.
00033 *
00034 *  Note that the routine returns V**T, not V.
00035 *
00036 *  Arguments
00037 *  =========
00038 *
00039 *  JOBU    (input) CHARACTER*1
00040 *          Specifies options for computing all or part of the matrix U:
00041 *          = 'A':  all M columns of U are returned in array U:
00042 *          = 'S':  the first min(m,n) columns of U (the left singular
00043 *                  vectors) are returned in the array U;
00044 *          = 'O':  the first min(m,n) columns of U (the left singular
00045 *                  vectors) are overwritten on the array A;
00046 *          = 'N':  no columns of U (no left singular vectors) are
00047 *                  computed.
00048 *
00049 *  JOBVT   (input) CHARACTER*1
00050 *          Specifies options for computing all or part of the matrix
00051 *          V**T:
00052 *          = 'A':  all N rows of V**T are returned in the array VT;
00053 *          = 'S':  the first min(m,n) rows of V**T (the right singular
00054 *                  vectors) are returned in the array VT;
00055 *          = 'O':  the first min(m,n) rows of V**T (the right singular
00056 *                  vectors) are overwritten on the array A;
00057 *          = 'N':  no rows of V**T (no right singular vectors) are
00058 *                  computed.
00059 *
00060 *          JOBVT and JOBU cannot both be 'O'.
00061 *
00062 *  M       (input) INTEGER
00063 *          The number of rows of the input matrix A.  M >= 0.
00064 *
00065 *  N       (input) INTEGER
00066 *          The number of columns of the input matrix A.  N >= 0.
00067 *
00068 *  A       (input/output) REAL array, dimension (LDA,N)
00069 *          On entry, the M-by-N matrix A.
00070 *          On exit,
00071 *          if JOBU = 'O',  A is overwritten with the first min(m,n)
00072 *                          columns of U (the left singular vectors,
00073 *                          stored columnwise);
00074 *          if JOBVT = 'O', A is overwritten with the first min(m,n)
00075 *                          rows of V**T (the right singular vectors,
00076 *                          stored rowwise);
00077 *          if JOBU .ne. 'O' and JOBVT .ne. 'O', the contents of A
00078 *                          are destroyed.
00079 *
00080 *  LDA     (input) INTEGER
00081 *          The leading dimension of the array A.  LDA >= max(1,M).
00082 *
00083 *  S       (output) REAL array, dimension (min(M,N))
00084 *          The singular values of A, sorted so that S(i) >= S(i+1).
00085 *
00086 *  U       (output) REAL array, dimension (LDU,UCOL)
00087 *          (LDU,M) if JOBU = 'A' or (LDU,min(M,N)) if JOBU = 'S'.
00088 *          If JOBU = 'A', U contains the M-by-M orthogonal matrix U;
00089 *          if JOBU = 'S', U contains the first min(m,n) columns of U
00090 *          (the left singular vectors, stored columnwise);
00091 *          if JOBU = 'N' or 'O', U is not referenced.
00092 *
00093 *  LDU     (input) INTEGER
00094 *          The leading dimension of the array U.  LDU >= 1; if
00095 *          JOBU = 'S' or 'A', LDU >= M.
00096 *
00097 *  VT      (output) REAL array, dimension (LDVT,N)
00098 *          If JOBVT = 'A', VT contains the N-by-N orthogonal matrix
00099 *          V**T;
00100 *          if JOBVT = 'S', VT contains the first min(m,n) rows of
00101 *          V**T (the right singular vectors, stored rowwise);
00102 *          if JOBVT = 'N' or 'O', VT is not referenced.
00103 *
00104 *  LDVT    (input) INTEGER
00105 *          The leading dimension of the array VT.  LDVT >= 1; if
00106 *          JOBVT = 'A', LDVT >= N; if JOBVT = 'S', LDVT >= min(M,N).
00107 *
00108 *  WORK    (workspace/output) REAL array, dimension (MAX(1,LWORK))
00109 *          On exit, if INFO = 0, WORK(1) returns the optimal LWORK;
00110 *          if INFO > 0, WORK(2:MIN(M,N)) contains the unconverged
00111 *          superdiagonal elements of an upper bidiagonal matrix B
00112 *          whose diagonal is in S (not necessarily sorted). B
00113 *          satisfies A = U * B * VT, so it has the same singular values
00114 *          as A, and singular vectors related by U and VT.
00115 *
00116 *  LWORK   (input) INTEGER
00117 *          The dimension of the array WORK.
00118 *          LWORK >= MAX(1,3*MIN(M,N)+MAX(M,N),5*MIN(M,N)).
00119 *          For good performance, LWORK should generally be larger.
00120 *
00121 *          If LWORK = -1, then a workspace query is assumed; the routine
00122 *          only calculates the optimal size of the WORK array, returns
00123 *          this value as the first entry of the WORK array, and no error
00124 *          message related to LWORK is issued by XERBLA.
00125 *
00126 *  INFO    (output) INTEGER
00127 *          = 0:  successful exit.
00128 *          < 0:  if INFO = -i, the i-th argument had an illegal value.
00129 *          > 0:  if SBDSQR did not converge, INFO specifies how many
00130 *                superdiagonals of an intermediate bidiagonal form B
00131 *                did not converge to zero. See the description of WORK
00132 *                above for details.
00133 *
00134 *  =====================================================================
00135 *
00136 *     .. Parameters ..
00137       REAL               ZERO, ONE
00138       PARAMETER          ( ZERO = 0.0E0, ONE = 1.0E0 )
00139 *     ..
00140 *     .. Local Scalars ..
00141       LOGICAL            LQUERY, WNTUA, WNTUAS, WNTUN, WNTUO, WNTUS,
00142      $                   WNTVA, WNTVAS, WNTVN, WNTVO, WNTVS
00143       INTEGER            BDSPAC, BLK, CHUNK, I, IE, IERR, IR, ISCL,
00144      $                   ITAU, ITAUP, ITAUQ, IU, IWORK, LDWRKR, LDWRKU,
00145      $                   MAXWRK, MINMN, MINWRK, MNTHR, NCU, NCVT, NRU,
00146      $                   NRVT, WRKBL
00147       REAL               ANRM, BIGNUM, EPS, SMLNUM
00148 *     ..
00149 *     .. Local Arrays ..
00150       REAL               DUM( 1 )
00151 *     ..
00152 *     .. External Subroutines ..
00153       EXTERNAL           SBDSQR, SGEBRD, SGELQF, SGEMM, SGEQRF, SLACPY,
00154      $                   SLASCL, SLASET, SORGBR, SORGLQ, SORGQR, SORMBR,
00155      $                   XERBLA
00156 *     ..
00157 *     .. External Functions ..
00158       LOGICAL            LSAME
00159       INTEGER            ILAENV
00160       REAL               SLAMCH, SLANGE
00161       EXTERNAL           LSAME, ILAENV, SLAMCH, SLANGE
00162 *     ..
00163 *     .. Intrinsic Functions ..
00164       INTRINSIC          MAX, MIN, SQRT
00165 *     ..
00166 *     .. Executable Statements ..
00167 *
00168 *     Test the input arguments
00169 *
00170       INFO = 0
00171       MINMN = MIN( M, N )
00172       WNTUA = LSAME( JOBU, 'A' )
00173       WNTUS = LSAME( JOBU, 'S' )
00174       WNTUAS = WNTUA .OR. WNTUS
00175       WNTUO = LSAME( JOBU, 'O' )
00176       WNTUN = LSAME( JOBU, 'N' )
00177       WNTVA = LSAME( JOBVT, 'A' )
00178       WNTVS = LSAME( JOBVT, 'S' )
00179       WNTVAS = WNTVA .OR. WNTVS
00180       WNTVO = LSAME( JOBVT, 'O' )
00181       WNTVN = LSAME( JOBVT, 'N' )
00182       LQUERY = ( LWORK.EQ.-1 )
00183 *
00184       IF( .NOT.( WNTUA .OR. WNTUS .OR. WNTUO .OR. WNTUN ) ) THEN
00185          INFO = -1
00186       ELSE IF( .NOT.( WNTVA .OR. WNTVS .OR. WNTVO .OR. WNTVN ) .OR.
00187      $         ( WNTVO .AND. WNTUO ) ) THEN
00188          INFO = -2
00189       ELSE IF( M.LT.0 ) THEN
00190          INFO = -3
00191       ELSE IF( N.LT.0 ) THEN
00192          INFO = -4
00193       ELSE IF( LDA.LT.MAX( 1, M ) ) THEN
00194          INFO = -6
00195       ELSE IF( LDU.LT.1 .OR. ( WNTUAS .AND. LDU.LT.M ) ) THEN
00196          INFO = -9
00197       ELSE IF( LDVT.LT.1 .OR. ( WNTVA .AND. LDVT.LT.N ) .OR.
00198      $         ( WNTVS .AND. LDVT.LT.MINMN ) ) THEN
00199          INFO = -11
00200       END IF
00201 *
00202 *     Compute workspace
00203 *      (Note: Comments in the code beginning "Workspace:" describe the
00204 *       minimal amount of workspace needed at that point in the code,
00205 *       as well as the preferred amount for good performance.
00206 *       NB refers to the optimal block size for the immediately
00207 *       following subroutine, as returned by ILAENV.)
00208 *
00209       IF( INFO.EQ.0 ) THEN
00210          MINWRK = 1
00211          MAXWRK = 1
00212          IF( M.GE.N .AND. MINMN.GT.0 ) THEN
00213 *
00214 *           Compute space needed for SBDSQR
00215 *
00216             MNTHR = ILAENV( 6, 'SGESVD', JOBU // JOBVT, M, N, 0, 0 )
00217             BDSPAC = 5*N
00218             IF( M.GE.MNTHR ) THEN
00219                IF( WNTUN ) THEN
00220 *
00221 *                 Path 1 (M much larger than N, JOBU='N')
00222 *
00223                   MAXWRK = N + N*ILAENV( 1, 'SGEQRF', ' ', M, N, -1,
00224      $                     -1 )
00225                   MAXWRK = MAX( MAXWRK, 3*N+2*N*
00226      $                     ILAENV( 1, 'SGEBRD', ' ', N, N, -1, -1 ) )
00227                   IF( WNTVO .OR. WNTVAS )
00228      $               MAXWRK = MAX( MAXWRK, 3*N+( N-1 )*
00229      $                        ILAENV( 1, 'SORGBR', 'P', N, N, N, -1 ) )
00230                   MAXWRK = MAX( MAXWRK, BDSPAC )
00231                   MINWRK = MAX( 4*N, BDSPAC )
00232                ELSE IF( WNTUO .AND. WNTVN ) THEN
00233 *
00234 *                 Path 2 (M much larger than N, JOBU='O', JOBVT='N')
00235 *
00236                   WRKBL = N + N*ILAENV( 1, 'SGEQRF', ' ', M, N, -1, -1 )
00237                   WRKBL = MAX( WRKBL, N+N*ILAENV( 1, 'SORGQR', ' ', M,
00238      $                    N, N, -1 ) )
00239                   WRKBL = MAX( WRKBL, 3*N+2*N*
00240      $                    ILAENV( 1, 'SGEBRD', ' ', N, N, -1, -1 ) )
00241                   WRKBL = MAX( WRKBL, 3*N+N*
00242      $                    ILAENV( 1, 'SORGBR', 'Q', N, N, N, -1 ) )
00243                   WRKBL = MAX( WRKBL, BDSPAC )
00244                   MAXWRK = MAX( N*N+WRKBL, N*N+M*N+N )
00245                   MINWRK = MAX( 3*N+M, BDSPAC )
00246                ELSE IF( WNTUO .AND. WNTVAS ) THEN
00247 *
00248 *                 Path 3 (M much larger than N, JOBU='O', JOBVT='S' or
00249 *                 'A')
00250 *
00251                   WRKBL = N + N*ILAENV( 1, 'SGEQRF', ' ', M, N, -1, -1 )
00252                   WRKBL = MAX( WRKBL, N+N*ILAENV( 1, 'SORGQR', ' ', M,
00253      $                    N, N, -1 ) )
00254                   WRKBL = MAX( WRKBL, 3*N+2*N*
00255      $                    ILAENV( 1, 'SGEBRD', ' ', N, N, -1, -1 ) )
00256                   WRKBL = MAX( WRKBL, 3*N+N*
00257      $                    ILAENV( 1, 'SORGBR', 'Q', N, N, N, -1 ) )
00258                   WRKBL = MAX( WRKBL, 3*N+( N-1 )*
00259      $                    ILAENV( 1, 'SORGBR', 'P', N, N, N, -1 ) )
00260                   WRKBL = MAX( WRKBL, BDSPAC )
00261                   MAXWRK = MAX( N*N+WRKBL, N*N+M*N+N )
00262                   MINWRK = MAX( 3*N+M, BDSPAC )
00263                ELSE IF( WNTUS .AND. WNTVN ) THEN
00264 *
00265 *                 Path 4 (M much larger than N, JOBU='S', JOBVT='N')
00266 *
00267                   WRKBL = N + N*ILAENV( 1, 'SGEQRF', ' ', M, N, -1, -1 )
00268                   WRKBL = MAX( WRKBL, N+N*ILAENV( 1, 'SORGQR', ' ', M,
00269      $                    N, N, -1 ) )
00270                   WRKBL = MAX( WRKBL, 3*N+2*N*
00271      $                    ILAENV( 1, 'SGEBRD', ' ', N, N, -1, -1 ) )
00272                   WRKBL = MAX( WRKBL, 3*N+N*
00273      $                    ILAENV( 1, 'SORGBR', 'Q', N, N, N, -1 ) )
00274                   WRKBL = MAX( WRKBL, BDSPAC )
00275                   MAXWRK = N*N + WRKBL
00276                   MINWRK = MAX( 3*N+M, BDSPAC )
00277                ELSE IF( WNTUS .AND. WNTVO ) THEN
00278 *
00279 *                 Path 5 (M much larger than N, JOBU='S', JOBVT='O')
00280 *
00281                   WRKBL = N + N*ILAENV( 1, 'SGEQRF', ' ', M, N, -1, -1 )
00282                   WRKBL = MAX( WRKBL, N+N*ILAENV( 1, 'SORGQR', ' ', M,
00283      $                    N, N, -1 ) )
00284                   WRKBL = MAX( WRKBL, 3*N+2*N*
00285      $                    ILAENV( 1, 'SGEBRD', ' ', N, N, -1, -1 ) )
00286                   WRKBL = MAX( WRKBL, 3*N+N*
00287      $                    ILAENV( 1, 'SORGBR', 'Q', N, N, N, -1 ) )
00288                   WRKBL = MAX( WRKBL, 3*N+( N-1 )*
00289      $                    ILAENV( 1, 'SORGBR', 'P', N, N, N, -1 ) )
00290                   WRKBL = MAX( WRKBL, BDSPAC )
00291                   MAXWRK = 2*N*N + WRKBL
00292                   MINWRK = MAX( 3*N+M, BDSPAC )
00293                ELSE IF( WNTUS .AND. WNTVAS ) THEN
00294 *
00295 *                 Path 6 (M much larger than N, JOBU='S', JOBVT='S' or
00296 *                 'A')
00297 *
00298                   WRKBL = N + N*ILAENV( 1, 'SGEQRF', ' ', M, N, -1, -1 )
00299                   WRKBL = MAX( WRKBL, N+N*ILAENV( 1, 'SORGQR', ' ', M,
00300      $                    N, N, -1 ) )
00301                   WRKBL = MAX( WRKBL, 3*N+2*N*
00302      $                    ILAENV( 1, 'SGEBRD', ' ', N, N, -1, -1 ) )
00303                   WRKBL = MAX( WRKBL, 3*N+N*
00304      $                    ILAENV( 1, 'SORGBR', 'Q', N, N, N, -1 ) )
00305                   WRKBL = MAX( WRKBL, 3*N+( N-1 )*
00306      $                    ILAENV( 1, 'SORGBR', 'P', N, N, N, -1 ) )
00307                   WRKBL = MAX( WRKBL, BDSPAC )
00308                   MAXWRK = N*N + WRKBL
00309                   MINWRK = MAX( 3*N+M, BDSPAC )
00310                ELSE IF( WNTUA .AND. WNTVN ) THEN
00311 *
00312 *                 Path 7 (M much larger than N, JOBU='A', JOBVT='N')
00313 *
00314                   WRKBL = N + N*ILAENV( 1, 'SGEQRF', ' ', M, N, -1, -1 )
00315                   WRKBL = MAX( WRKBL, N+M*ILAENV( 1, 'SORGQR', ' ', M,
00316      $                    M, N, -1 ) )
00317                   WRKBL = MAX( WRKBL, 3*N+2*N*
00318      $                    ILAENV( 1, 'SGEBRD', ' ', N, N, -1, -1 ) )
00319                   WRKBL = MAX( WRKBL, 3*N+N*
00320      $                    ILAENV( 1, 'SORGBR', 'Q', N, N, N, -1 ) )
00321                   WRKBL = MAX( WRKBL, BDSPAC )
00322                   MAXWRK = N*N + WRKBL
00323                   MINWRK = MAX( 3*N+M, BDSPAC )
00324                ELSE IF( WNTUA .AND. WNTVO ) THEN
00325 *
00326 *                 Path 8 (M much larger than N, JOBU='A', JOBVT='O')
00327 *
00328                   WRKBL = N + N*ILAENV( 1, 'SGEQRF', ' ', M, N, -1, -1 )
00329                   WRKBL = MAX( WRKBL, N+M*ILAENV( 1, 'SORGQR', ' ', M,
00330      $                    M, N, -1 ) )
00331                   WRKBL = MAX( WRKBL, 3*N+2*N*
00332      $                    ILAENV( 1, 'SGEBRD', ' ', N, N, -1, -1 ) )
00333                   WRKBL = MAX( WRKBL, 3*N+N*
00334      $                    ILAENV( 1, 'SORGBR', 'Q', N, N, N, -1 ) )
00335                   WRKBL = MAX( WRKBL, 3*N+( N-1 )*
00336      $                    ILAENV( 1, 'SORGBR', 'P', N, N, N, -1 ) )
00337                   WRKBL = MAX( WRKBL, BDSPAC )
00338                   MAXWRK = 2*N*N + WRKBL
00339                   MINWRK = MAX( 3*N+M, BDSPAC )
00340                ELSE IF( WNTUA .AND. WNTVAS ) THEN
00341 *
00342 *                 Path 9 (M much larger than N, JOBU='A', JOBVT='S' or
00343 *                 'A')
00344 *
00345                   WRKBL = N + N*ILAENV( 1, 'SGEQRF', ' ', M, N, -1, -1 )
00346                   WRKBL = MAX( WRKBL, N+M*ILAENV( 1, 'SORGQR', ' ', M,
00347      $                    M, N, -1 ) )
00348                   WRKBL = MAX( WRKBL, 3*N+2*N*
00349      $                    ILAENV( 1, 'SGEBRD', ' ', N, N, -1, -1 ) )
00350                   WRKBL = MAX( WRKBL, 3*N+N*
00351      $                    ILAENV( 1, 'SORGBR', 'Q', N, N, N, -1 ) )
00352                   WRKBL = MAX( WRKBL, 3*N+( N-1 )*
00353      $                    ILAENV( 1, 'SORGBR', 'P', N, N, N, -1 ) )
00354                   WRKBL = MAX( WRKBL, BDSPAC )
00355                   MAXWRK = N*N + WRKBL
00356                   MINWRK = MAX( 3*N+M, BDSPAC )
00357                END IF
00358             ELSE
00359 *
00360 *              Path 10 (M at least N, but not much larger)
00361 *
00362                MAXWRK = 3*N + ( M+N )*ILAENV( 1, 'SGEBRD', ' ', M, N,
00363      $                  -1, -1 )
00364                IF( WNTUS .OR. WNTUO )
00365      $            MAXWRK = MAX( MAXWRK, 3*N+N*
00366      $                     ILAENV( 1, 'SORGBR', 'Q', M, N, N, -1 ) )
00367                IF( WNTUA )
00368      $            MAXWRK = MAX( MAXWRK, 3*N+M*
00369      $                     ILAENV( 1, 'SORGBR', 'Q', M, M, N, -1 ) )
00370                IF( .NOT.WNTVN )
00371      $            MAXWRK = MAX( MAXWRK, 3*N+( N-1 )*
00372      $                     ILAENV( 1, 'SORGBR', 'P', N, N, N, -1 ) )
00373                MAXWRK = MAX( MAXWRK, BDSPAC )
00374                MINWRK = MAX( 3*N+M, BDSPAC )
00375             END IF
00376          ELSE IF( MINMN.GT.0 ) THEN
00377 *
00378 *           Compute space needed for SBDSQR
00379 *
00380             MNTHR = ILAENV( 6, 'SGESVD', JOBU // JOBVT, M, N, 0, 0 )
00381             BDSPAC = 5*M
00382             IF( N.GE.MNTHR ) THEN
00383                IF( WNTVN ) THEN
00384 *
00385 *                 Path 1t(N much larger than M, JOBVT='N')
00386 *
00387                   MAXWRK = M + M*ILAENV( 1, 'SGELQF', ' ', M, N, -1,
00388      $                     -1 )
00389                   MAXWRK = MAX( MAXWRK, 3*M+2*M*
00390      $                     ILAENV( 1, 'SGEBRD', ' ', M, M, -1, -1 ) )
00391                   IF( WNTUO .OR. WNTUAS )
00392      $               MAXWRK = MAX( MAXWRK, 3*M+M*
00393      $                        ILAENV( 1, 'SORGBR', 'Q', M, M, M, -1 ) )
00394                   MAXWRK = MAX( MAXWRK, BDSPAC )
00395                   MINWRK = MAX( 4*M, BDSPAC )
00396                ELSE IF( WNTVO .AND. WNTUN ) THEN
00397 *
00398 *                 Path 2t(N much larger than M, JOBU='N', JOBVT='O')
00399 *
00400                   WRKBL = M + M*ILAENV( 1, 'SGELQF', ' ', M, N, -1, -1 )
00401                   WRKBL = MAX( WRKBL, M+M*ILAENV( 1, 'SORGLQ', ' ', M,
00402      $                    N, M, -1 ) )
00403                   WRKBL = MAX( WRKBL, 3*M+2*M*
00404      $                    ILAENV( 1, 'SGEBRD', ' ', M, M, -1, -1 ) )
00405                   WRKBL = MAX( WRKBL, 3*M+( M-1 )*
00406      $                    ILAENV( 1, 'SORGBR', 'P', M, M, M, -1 ) )
00407                   WRKBL = MAX( WRKBL, BDSPAC )
00408                   MAXWRK = MAX( M*M+WRKBL, M*M+M*N+M )
00409                   MINWRK = MAX( 3*M+N, BDSPAC )
00410                ELSE IF( WNTVO .AND. WNTUAS ) THEN
00411 *
00412 *                 Path 3t(N much larger than M, JOBU='S' or 'A',
00413 *                 JOBVT='O')
00414 *
00415                   WRKBL = M + M*ILAENV( 1, 'SGELQF', ' ', M, N, -1, -1 )
00416                   WRKBL = MAX( WRKBL, M+M*ILAENV( 1, 'SORGLQ', ' ', M,
00417      $                    N, M, -1 ) )
00418                   WRKBL = MAX( WRKBL, 3*M+2*M*
00419      $                    ILAENV( 1, 'SGEBRD', ' ', M, M, -1, -1 ) )
00420                   WRKBL = MAX( WRKBL, 3*M+( M-1 )*
00421      $                    ILAENV( 1, 'SORGBR', 'P', M, M, M, -1 ) )
00422                   WRKBL = MAX( WRKBL, 3*M+M*
00423      $                    ILAENV( 1, 'SORGBR', 'Q', M, M, M, -1 ) )
00424                   WRKBL = MAX( WRKBL, BDSPAC )
00425                   MAXWRK = MAX( M*M+WRKBL, M*M+M*N+M )
00426                   MINWRK = MAX( 3*M+N, BDSPAC )
00427                ELSE IF( WNTVS .AND. WNTUN ) THEN
00428 *
00429 *                 Path 4t(N much larger than M, JOBU='N', JOBVT='S')
00430 *
00431                   WRKBL = M + M*ILAENV( 1, 'SGELQF', ' ', M, N, -1, -1 )
00432                   WRKBL = MAX( WRKBL, M+M*ILAENV( 1, 'SORGLQ', ' ', M,
00433      $                    N, M, -1 ) )
00434                   WRKBL = MAX( WRKBL, 3*M+2*M*
00435      $                    ILAENV( 1, 'SGEBRD', ' ', M, M, -1, -1 ) )
00436                   WRKBL = MAX( WRKBL, 3*M+( M-1 )*
00437      $                    ILAENV( 1, 'SORGBR', 'P', M, M, M, -1 ) )
00438                   WRKBL = MAX( WRKBL, BDSPAC )
00439                   MAXWRK = M*M + WRKBL
00440                   MINWRK = MAX( 3*M+N, BDSPAC )
00441                ELSE IF( WNTVS .AND. WNTUO ) THEN
00442 *
00443 *                 Path 5t(N much larger than M, JOBU='O', JOBVT='S')
00444 *
00445                   WRKBL = M + M*ILAENV( 1, 'SGELQF', ' ', M, N, -1, -1 )
00446                   WRKBL = MAX( WRKBL, M+M*ILAENV( 1, 'SORGLQ', ' ', M,
00447      $                    N, M, -1 ) )
00448                   WRKBL = MAX( WRKBL, 3*M+2*M*
00449      $                    ILAENV( 1, 'SGEBRD', ' ', M, M, -1, -1 ) )
00450                   WRKBL = MAX( WRKBL, 3*M+( M-1 )*
00451      $                    ILAENV( 1, 'SORGBR', 'P', M, M, M, -1 ) )
00452                   WRKBL = MAX( WRKBL, 3*M+M*
00453      $                    ILAENV( 1, 'SORGBR', 'Q', M, M, M, -1 ) )
00454                   WRKBL = MAX( WRKBL, BDSPAC )
00455                   MAXWRK = 2*M*M + WRKBL
00456                   MINWRK = MAX( 3*M+N, BDSPAC )
00457                   MAXWRK = MAX( MAXWRK, MINWRK )
00458                ELSE IF( WNTVS .AND. WNTUAS ) THEN
00459 *
00460 *                 Path 6t(N much larger than M, JOBU='S' or 'A',
00461 *                 JOBVT='S')
00462 *
00463                   WRKBL = M + M*ILAENV( 1, 'SGELQF', ' ', M, N, -1, -1 )
00464                   WRKBL = MAX( WRKBL, M+M*ILAENV( 1, 'SORGLQ', ' ', M,
00465      $                    N, M, -1 ) )
00466                   WRKBL = MAX( WRKBL, 3*M+2*M*
00467      $                    ILAENV( 1, 'SGEBRD', ' ', M, M, -1, -1 ) )
00468                   WRKBL = MAX( WRKBL, 3*M+( M-1 )*
00469      $                    ILAENV( 1, 'SORGBR', 'P', M, M, M, -1 ) )
00470                   WRKBL = MAX( WRKBL, 3*M+M*
00471      $                    ILAENV( 1, 'SORGBR', 'Q', M, M, M, -1 ) )
00472                   WRKBL = MAX( WRKBL, BDSPAC )
00473                   MAXWRK = M*M + WRKBL
00474                   MINWRK = MAX( 3*M+N, BDSPAC )
00475                ELSE IF( WNTVA .AND. WNTUN ) THEN
00476 *
00477 *                 Path 7t(N much larger than M, JOBU='N', JOBVT='A')
00478 *
00479                   WRKBL = M + M*ILAENV( 1, 'SGELQF', ' ', M, N, -1, -1 )
00480                   WRKBL = MAX( WRKBL, M+N*ILAENV( 1, 'SORGLQ', ' ', N,
00481      $                    N, M, -1 ) )
00482                   WRKBL = MAX( WRKBL, 3*M+2*M*
00483      $                    ILAENV( 1, 'SGEBRD', ' ', M, M, -1, -1 ) )
00484                   WRKBL = MAX( WRKBL, 3*M+( M-1 )*
00485      $                    ILAENV( 1, 'SORGBR', 'P', M, M, M, -1 ) )
00486                   WRKBL = MAX( WRKBL, BDSPAC )
00487                   MAXWRK = M*M + WRKBL
00488                   MINWRK = MAX( 3*M+N, BDSPAC )
00489                ELSE IF( WNTVA .AND. WNTUO ) THEN
00490 *
00491 *                 Path 8t(N much larger than M, JOBU='O', JOBVT='A')
00492 *
00493                   WRKBL = M + M*ILAENV( 1, 'SGELQF', ' ', M, N, -1, -1 )
00494                   WRKBL = MAX( WRKBL, M+N*ILAENV( 1, 'SORGLQ', ' ', N,
00495      $                    N, M, -1 ) )
00496                   WRKBL = MAX( WRKBL, 3*M+2*M*
00497      $                    ILAENV( 1, 'SGEBRD', ' ', M, M, -1, -1 ) )
00498                   WRKBL = MAX( WRKBL, 3*M+( M-1 )*
00499      $                    ILAENV( 1, 'SORGBR', 'P', M, M, M, -1 ) )
00500                   WRKBL = MAX( WRKBL, 3*M+M*
00501      $                    ILAENV( 1, 'SORGBR', 'Q', M, M, M, -1 ) )
00502                   WRKBL = MAX( WRKBL, BDSPAC )
00503                   MAXWRK = 2*M*M + WRKBL
00504                   MINWRK = MAX( 3*M+N, BDSPAC )
00505                ELSE IF( WNTVA .AND. WNTUAS ) THEN
00506 *
00507 *                 Path 9t(N much larger than M, JOBU='S' or 'A',
00508 *                 JOBVT='A')
00509 *
00510                   WRKBL = M + M*ILAENV( 1, 'SGELQF', ' ', M, N, -1, -1 )
00511                   WRKBL = MAX( WRKBL, M+N*ILAENV( 1, 'SORGLQ', ' ', N,
00512      $                    N, M, -1 ) )
00513                   WRKBL = MAX( WRKBL, 3*M+2*M*
00514      $                    ILAENV( 1, 'SGEBRD', ' ', M, M, -1, -1 ) )
00515                   WRKBL = MAX( WRKBL, 3*M+( M-1 )*
00516      $                    ILAENV( 1, 'SORGBR', 'P', M, M, M, -1 ) )
00517                   WRKBL = MAX( WRKBL, 3*M+M*
00518      $                    ILAENV( 1, 'SORGBR', 'Q', M, M, M, -1 ) )
00519                   WRKBL = MAX( WRKBL, BDSPAC )
00520                   MAXWRK = M*M + WRKBL
00521                   MINWRK = MAX( 3*M+N, BDSPAC )
00522                END IF
00523             ELSE
00524 *
00525 *              Path 10t(N greater than M, but not much larger)
00526 *
00527                MAXWRK = 3*M + ( M+N )*ILAENV( 1, 'SGEBRD', ' ', M, N,
00528      $                  -1, -1 )
00529                IF( WNTVS .OR. WNTVO )
00530      $            MAXWRK = MAX( MAXWRK, 3*M+M*
00531      $                     ILAENV( 1, 'SORGBR', 'P', M, N, M, -1 ) )
00532                IF( WNTVA )
00533      $            MAXWRK = MAX( MAXWRK, 3*M+N*
00534      $                     ILAENV( 1, 'SORGBR', 'P', N, N, M, -1 ) )
00535                IF( .NOT.WNTUN )
00536      $            MAXWRK = MAX( MAXWRK, 3*M+( M-1 )*
00537      $                     ILAENV( 1, 'SORGBR', 'Q', M, M, M, -1 ) )
00538                MAXWRK = MAX( MAXWRK, BDSPAC )
00539                MINWRK = MAX( 3*M+N, BDSPAC )
00540             END IF
00541          END IF
00542          MAXWRK = MAX( MAXWRK, MINWRK )
00543          WORK( 1 ) = MAXWRK
00544 *
00545          IF( LWORK.LT.MINWRK .AND. .NOT.LQUERY ) THEN
00546             INFO = -13
00547          END IF
00548       END IF
00549 *
00550       IF( INFO.NE.0 ) THEN
00551          CALL XERBLA( 'SGESVD', -INFO )
00552          RETURN
00553       ELSE IF( LQUERY ) THEN
00554          RETURN
00555       END IF
00556 *
00557 *     Quick return if possible
00558 *
00559       IF( M.EQ.0 .OR. N.EQ.0 ) THEN
00560          RETURN
00561       END IF
00562 *
00563 *     Get machine constants
00564 *
00565       EPS = SLAMCH( 'P' )
00566       SMLNUM = SQRT( SLAMCH( 'S' ) ) / EPS
00567       BIGNUM = ONE / SMLNUM
00568 *
00569 *     Scale A if max element outside range [SMLNUM,BIGNUM]
00570 *
00571       ANRM = SLANGE( 'M', M, N, A, LDA, DUM )
00572       ISCL = 0
00573       IF( ANRM.GT.ZERO .AND. ANRM.LT.SMLNUM ) THEN
00574          ISCL = 1
00575          CALL SLASCL( 'G', 0, 0, ANRM, SMLNUM, M, N, A, LDA, IERR )
00576       ELSE IF( ANRM.GT.BIGNUM ) THEN
00577          ISCL = 1
00578          CALL SLASCL( 'G', 0, 0, ANRM, BIGNUM, M, N, A, LDA, IERR )
00579       END IF
00580 *
00581       IF( M.GE.N ) THEN
00582 *
00583 *        A has at least as many rows as columns. If A has sufficiently
00584 *        more rows than columns, first reduce using the QR
00585 *        decomposition (if sufficient workspace available)
00586 *
00587          IF( M.GE.MNTHR ) THEN
00588 *
00589             IF( WNTUN ) THEN
00590 *
00591 *              Path 1 (M much larger than N, JOBU='N')
00592 *              No left singular vectors to be computed
00593 *
00594                ITAU = 1
00595                IWORK = ITAU + N
00596 *
00597 *              Compute A=Q*R
00598 *              (Workspace: need 2*N, prefer N+N*NB)
00599 *
00600                CALL SGEQRF( M, N, A, LDA, WORK( ITAU ), WORK( IWORK ),
00601      $                      LWORK-IWORK+1, IERR )
00602 *
00603 *              Zero out below R
00604 *
00605                CALL SLASET( 'L', N-1, N-1, ZERO, ZERO, A( 2, 1 ), LDA )
00606                IE = 1
00607                ITAUQ = IE + N
00608                ITAUP = ITAUQ + N
00609                IWORK = ITAUP + N
00610 *
00611 *              Bidiagonalize R in A
00612 *              (Workspace: need 4*N, prefer 3*N+2*N*NB)
00613 *
00614                CALL SGEBRD( N, N, A, LDA, S, WORK( IE ), WORK( ITAUQ ),
00615      $                      WORK( ITAUP ), WORK( IWORK ), LWORK-IWORK+1,
00616      $                      IERR )
00617                NCVT = 0
00618                IF( WNTVO .OR. WNTVAS ) THEN
00619 *
00620 *                 If right singular vectors desired, generate P'.
00621 *                 (Workspace: need 4*N-1, prefer 3*N+(N-1)*NB)
00622 *
00623                   CALL SORGBR( 'P', N, N, N, A, LDA, WORK( ITAUP ),
00624      $                         WORK( IWORK ), LWORK-IWORK+1, IERR )
00625                   NCVT = N
00626                END IF
00627                IWORK = IE + N
00628 *
00629 *              Perform bidiagonal QR iteration, computing right
00630 *              singular vectors of A in A if desired
00631 *              (Workspace: need BDSPAC)
00632 *
00633                CALL SBDSQR( 'U', N, NCVT, 0, 0, S, WORK( IE ), A, LDA,
00634      $                      DUM, 1, DUM, 1, WORK( IWORK ), INFO )
00635 *
00636 *              If right singular vectors desired in VT, copy them there
00637 *
00638                IF( WNTVAS )
00639      $            CALL SLACPY( 'F', N, N, A, LDA, VT, LDVT )
00640 *
00641             ELSE IF( WNTUO .AND. WNTVN ) THEN
00642 *
00643 *              Path 2 (M much larger than N, JOBU='O', JOBVT='N')
00644 *              N left singular vectors to be overwritten on A and
00645 *              no right singular vectors to be computed
00646 *
00647                IF( LWORK.GE.N*N+MAX( 4*N, BDSPAC ) ) THEN
00648 *
00649 *                 Sufficient workspace for a fast algorithm
00650 *
00651                   IR = 1
00652                   IF( LWORK.GE.MAX( WRKBL, LDA*N+N )+LDA*N ) THEN
00653 *
00654 *                    WORK(IU) is LDA by N, WORK(IR) is LDA by N
00655 *
00656                      LDWRKU = LDA
00657                      LDWRKR = LDA
00658                   ELSE IF( LWORK.GE.MAX( WRKBL, LDA*N+N )+N*N ) THEN
00659 *
00660 *                    WORK(IU) is LDA by N, WORK(IR) is N by N
00661 *
00662                      LDWRKU = LDA
00663                      LDWRKR = N
00664                   ELSE
00665 *
00666 *                    WORK(IU) is LDWRKU by N, WORK(IR) is N by N
00667 *
00668                      LDWRKU = ( LWORK-N*N-N ) / N
00669                      LDWRKR = N
00670                   END IF
00671                   ITAU = IR + LDWRKR*N
00672                   IWORK = ITAU + N
00673 *
00674 *                 Compute A=Q*R
00675 *                 (Workspace: need N*N+2*N, prefer N*N+N+N*NB)
00676 *
00677                   CALL SGEQRF( M, N, A, LDA, WORK( ITAU ),
00678      $                         WORK( IWORK ), LWORK-IWORK+1, IERR )
00679 *
00680 *                 Copy R to WORK(IR) and zero out below it
00681 *
00682                   CALL SLACPY( 'U', N, N, A, LDA, WORK( IR ), LDWRKR )
00683                   CALL SLASET( 'L', N-1, N-1, ZERO, ZERO, WORK( IR+1 ),
00684      $                         LDWRKR )
00685 *
00686 *                 Generate Q in A
00687 *                 (Workspace: need N*N+2*N, prefer N*N+N+N*NB)
00688 *
00689                   CALL SORGQR( M, N, N, A, LDA, WORK( ITAU ),
00690      $                         WORK( IWORK ), LWORK-IWORK+1, IERR )
00691                   IE = ITAU
00692                   ITAUQ = IE + N
00693                   ITAUP = ITAUQ + N
00694                   IWORK = ITAUP + N
00695 *
00696 *                 Bidiagonalize R in WORK(IR)
00697 *                 (Workspace: need N*N+4*N, prefer N*N+3*N+2*N*NB)
00698 *
00699                   CALL SGEBRD( N, N, WORK( IR ), LDWRKR, S, WORK( IE ),
00700      $                         WORK( ITAUQ ), WORK( ITAUP ),
00701      $                         WORK( IWORK ), LWORK-IWORK+1, IERR )
00702 *
00703 *                 Generate left vectors bidiagonalizing R
00704 *                 (Workspace: need N*N+4*N, prefer N*N+3*N+N*NB)
00705 *
00706                   CALL SORGBR( 'Q', N, N, N, WORK( IR ), LDWRKR,
00707      $                         WORK( ITAUQ ), WORK( IWORK ),
00708      $                         LWORK-IWORK+1, IERR )
00709                   IWORK = IE + N
00710 *
00711 *                 Perform bidiagonal QR iteration, computing left
00712 *                 singular vectors of R in WORK(IR)
00713 *                 (Workspace: need N*N+BDSPAC)
00714 *
00715                   CALL SBDSQR( 'U', N, 0, N, 0, S, WORK( IE ), DUM, 1,
00716      $                         WORK( IR ), LDWRKR, DUM, 1,
00717      $                         WORK( IWORK ), INFO )
00718                   IU = IE + N
00719 *
00720 *                 Multiply Q in A by left singular vectors of R in
00721 *                 WORK(IR), storing result in WORK(IU) and copying to A
00722 *                 (Workspace: need N*N+2*N, prefer N*N+M*N+N)
00723 *
00724                   DO 10 I = 1, M, LDWRKU
00725                      CHUNK = MIN( M-I+1, LDWRKU )
00726                      CALL SGEMM( 'N', 'N', CHUNK, N, N, ONE, A( I, 1 ),
00727      $                           LDA, WORK( IR ), LDWRKR, ZERO,
00728      $                           WORK( IU ), LDWRKU )
00729                      CALL SLACPY( 'F', CHUNK, N, WORK( IU ), LDWRKU,
00730      $                            A( I, 1 ), LDA )
00731    10             CONTINUE
00732 *
00733                ELSE
00734 *
00735 *                 Insufficient workspace for a fast algorithm
00736 *
00737                   IE = 1
00738                   ITAUQ = IE + N
00739                   ITAUP = ITAUQ + N
00740                   IWORK = ITAUP + N
00741 *
00742 *                 Bidiagonalize A
00743 *                 (Workspace: need 3*N+M, prefer 3*N+(M+N)*NB)
00744 *
00745                   CALL SGEBRD( M, N, A, LDA, S, WORK( IE ),
00746      $                         WORK( ITAUQ ), WORK( ITAUP ),
00747      $                         WORK( IWORK ), LWORK-IWORK+1, IERR )
00748 *
00749 *                 Generate left vectors bidiagonalizing A
00750 *                 (Workspace: need 4*N, prefer 3*N+N*NB)
00751 *
00752                   CALL SORGBR( 'Q', M, N, N, A, LDA, WORK( ITAUQ ),
00753      $                         WORK( IWORK ), LWORK-IWORK+1, IERR )
00754                   IWORK = IE + N
00755 *
00756 *                 Perform bidiagonal QR iteration, computing left
00757 *                 singular vectors of A in A
00758 *                 (Workspace: need BDSPAC)
00759 *
00760                   CALL SBDSQR( 'U', N, 0, M, 0, S, WORK( IE ), DUM, 1,
00761      $                         A, LDA, DUM, 1, WORK( IWORK ), INFO )
00762 *
00763                END IF
00764 *
00765             ELSE IF( WNTUO .AND. WNTVAS ) THEN
00766 *
00767 *              Path 3 (M much larger than N, JOBU='O', JOBVT='S' or 'A')
00768 *              N left singular vectors to be overwritten on A and
00769 *              N right singular vectors to be computed in VT
00770 *
00771                IF( LWORK.GE.N*N+MAX( 4*N, BDSPAC ) ) THEN
00772 *
00773 *                 Sufficient workspace for a fast algorithm
00774 *
00775                   IR = 1
00776                   IF( LWORK.GE.MAX( WRKBL, LDA*N+N )+LDA*N ) THEN
00777 *
00778 *                    WORK(IU) is LDA by N and WORK(IR) is LDA by N
00779 *
00780                      LDWRKU = LDA
00781                      LDWRKR = LDA
00782                   ELSE IF( LWORK.GE.MAX( WRKBL, LDA*N+N )+N*N ) THEN
00783 *
00784 *                    WORK(IU) is LDA by N and WORK(IR) is N by N
00785 *
00786                      LDWRKU = LDA
00787                      LDWRKR = N
00788                   ELSE
00789 *
00790 *                    WORK(IU) is LDWRKU by N and WORK(IR) is N by N
00791 *
00792                      LDWRKU = ( LWORK-N*N-N ) / N
00793                      LDWRKR = N
00794                   END IF
00795                   ITAU = IR + LDWRKR*N
00796                   IWORK = ITAU + N
00797 *
00798 *                 Compute A=Q*R
00799 *                 (Workspace: need N*N+2*N, prefer N*N+N+N*NB)
00800 *
00801                   CALL SGEQRF( M, N, A, LDA, WORK( ITAU ),
00802      $                         WORK( IWORK ), LWORK-IWORK+1, IERR )
00803 *
00804 *                 Copy R to VT, zeroing out below it
00805 *
00806                   CALL SLACPY( 'U', N, N, A, LDA, VT, LDVT )
00807                   IF( N.GT.1 )
00808      $               CALL SLASET( 'L', N-1, N-1, ZERO, ZERO,
00809      $                            VT( 2, 1 ), LDVT )
00810 *
00811 *                 Generate Q in A
00812 *                 (Workspace: need N*N+2*N, prefer N*N+N+N*NB)
00813 *
00814                   CALL SORGQR( M, N, N, A, LDA, WORK( ITAU ),
00815      $                         WORK( IWORK ), LWORK-IWORK+1, IERR )
00816                   IE = ITAU
00817                   ITAUQ = IE + N
00818                   ITAUP = ITAUQ + N
00819                   IWORK = ITAUP + N
00820 *
00821 *                 Bidiagonalize R in VT, copying result to WORK(IR)
00822 *                 (Workspace: need N*N+4*N, prefer N*N+3*N+2*N*NB)
00823 *
00824                   CALL SGEBRD( N, N, VT, LDVT, S, WORK( IE ),
00825      $                         WORK( ITAUQ ), WORK( ITAUP ),
00826      $                         WORK( IWORK ), LWORK-IWORK+1, IERR )
00827                   CALL SLACPY( 'L', N, N, VT, LDVT, WORK( IR ), LDWRKR )
00828 *
00829 *                 Generate left vectors bidiagonalizing R in WORK(IR)
00830 *                 (Workspace: need N*N+4*N, prefer N*N+3*N+N*NB)
00831 *
00832                   CALL SORGBR( 'Q', N, N, N, WORK( IR ), LDWRKR,
00833      $                         WORK( ITAUQ ), WORK( IWORK ),
00834      $                         LWORK-IWORK+1, IERR )
00835 *
00836 *                 Generate right vectors bidiagonalizing R in VT
00837 *                 (Workspace: need N*N+4*N-1, prefer N*N+3*N+(N-1)*NB)
00838 *
00839                   CALL SORGBR( 'P', N, N, N, VT, LDVT, WORK( ITAUP ),
00840      $                         WORK( IWORK ), LWORK-IWORK+1, IERR )
00841                   IWORK = IE + N
00842 *
00843 *                 Perform bidiagonal QR iteration, computing left
00844 *                 singular vectors of R in WORK(IR) and computing right
00845 *                 singular vectors of R in VT
00846 *                 (Workspace: need N*N+BDSPAC)
00847 *
00848                   CALL SBDSQR( 'U', N, N, N, 0, S, WORK( IE ), VT, LDVT,
00849      $                         WORK( IR ), LDWRKR, DUM, 1,
00850      $                         WORK( IWORK ), INFO )
00851                   IU = IE + N
00852 *
00853 *                 Multiply Q in A by left singular vectors of R in
00854 *                 WORK(IR), storing result in WORK(IU) and copying to A
00855 *                 (Workspace: need N*N+2*N, prefer N*N+M*N+N)
00856 *
00857                   DO 20 I = 1, M, LDWRKU
00858                      CHUNK = MIN( M-I+1, LDWRKU )
00859                      CALL SGEMM( 'N', 'N', CHUNK, N, N, ONE, A( I, 1 ),
00860      $                           LDA, WORK( IR ), LDWRKR, ZERO,
00861      $                           WORK( IU ), LDWRKU )
00862                      CALL SLACPY( 'F', CHUNK, N, WORK( IU ), LDWRKU,
00863      $                            A( I, 1 ), LDA )
00864    20             CONTINUE
00865 *
00866                ELSE
00867 *
00868 *                 Insufficient workspace for a fast algorithm
00869 *
00870                   ITAU = 1
00871                   IWORK = ITAU + N
00872 *
00873 *                 Compute A=Q*R
00874 *                 (Workspace: need 2*N, prefer N+N*NB)
00875 *
00876                   CALL SGEQRF( M, N, A, LDA, WORK( ITAU ),
00877      $                         WORK( IWORK ), LWORK-IWORK+1, IERR )
00878 *
00879 *                 Copy R to VT, zeroing out below it
00880 *
00881                   CALL SLACPY( 'U', N, N, A, LDA, VT, LDVT )
00882                   IF( N.GT.1 )
00883      $               CALL SLASET( 'L', N-1, N-1, ZERO, ZERO,
00884      $                            VT( 2, 1 ), LDVT )
00885 *
00886 *                 Generate Q in A
00887 *                 (Workspace: need 2*N, prefer N+N*NB)
00888 *
00889                   CALL SORGQR( M, N, N, A, LDA, WORK( ITAU ),
00890      $                         WORK( IWORK ), LWORK-IWORK+1, IERR )
00891                   IE = ITAU
00892                   ITAUQ = IE + N
00893                   ITAUP = ITAUQ + N
00894                   IWORK = ITAUP + N
00895 *
00896 *                 Bidiagonalize R in VT
00897 *                 (Workspace: need 4*N, prefer 3*N+2*N*NB)
00898 *
00899                   CALL SGEBRD( N, N, VT, LDVT, S, WORK( IE ),
00900      $                         WORK( ITAUQ ), WORK( ITAUP ),
00901      $                         WORK( IWORK ), LWORK-IWORK+1, IERR )
00902 *
00903 *                 Multiply Q in A by left vectors bidiagonalizing R
00904 *                 (Workspace: need 3*N+M, prefer 3*N+M*NB)
00905 *
00906                   CALL SORMBR( 'Q', 'R', 'N', M, N, N, VT, LDVT,
00907      $                         WORK( ITAUQ ), A, LDA, WORK( IWORK ),
00908      $                         LWORK-IWORK+1, IERR )
00909 *
00910 *                 Generate right vectors bidiagonalizing R in VT
00911 *                 (Workspace: need 4*N-1, prefer 3*N+(N-1)*NB)
00912 *
00913                   CALL SORGBR( 'P', N, N, N, VT, LDVT, WORK( ITAUP ),
00914      $                         WORK( IWORK ), LWORK-IWORK+1, IERR )
00915                   IWORK = IE + N
00916 *
00917 *                 Perform bidiagonal QR iteration, computing left
00918 *                 singular vectors of A in A and computing right
00919 *                 singular vectors of A in VT
00920 *                 (Workspace: need BDSPAC)
00921 *
00922                   CALL SBDSQR( 'U', N, N, M, 0, S, WORK( IE ), VT, LDVT,
00923      $                         A, LDA, DUM, 1, WORK( IWORK ), INFO )
00924 *
00925                END IF
00926 *
00927             ELSE IF( WNTUS ) THEN
00928 *
00929                IF( WNTVN ) THEN
00930 *
00931 *                 Path 4 (M much larger than N, JOBU='S', JOBVT='N')
00932 *                 N left singular vectors to be computed in U and
00933 *                 no right singular vectors to be computed
00934 *
00935                   IF( LWORK.GE.N*N+MAX( 4*N, BDSPAC ) ) THEN
00936 *
00937 *                    Sufficient workspace for a fast algorithm
00938 *
00939                      IR = 1
00940                      IF( LWORK.GE.WRKBL+LDA*N ) THEN
00941 *
00942 *                       WORK(IR) is LDA by N
00943 *
00944                         LDWRKR = LDA
00945                      ELSE
00946 *
00947 *                       WORK(IR) is N by N
00948 *
00949                         LDWRKR = N
00950                      END IF
00951                      ITAU = IR + LDWRKR*N
00952                      IWORK = ITAU + N
00953 *
00954 *                    Compute A=Q*R
00955 *                    (Workspace: need N*N+2*N, prefer N*N+N+N*NB)
00956 *
00957                      CALL SGEQRF( M, N, A, LDA, WORK( ITAU ),
00958      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
00959 *
00960 *                    Copy R to WORK(IR), zeroing out below it
00961 *
00962                      CALL SLACPY( 'U', N, N, A, LDA, WORK( IR ),
00963      $                            LDWRKR )
00964                      CALL SLASET( 'L', N-1, N-1, ZERO, ZERO,
00965      $                            WORK( IR+1 ), LDWRKR )
00966 *
00967 *                    Generate Q in A
00968 *                    (Workspace: need N*N+2*N, prefer N*N+N+N*NB)
00969 *
00970                      CALL SORGQR( M, N, N, A, LDA, WORK( ITAU ),
00971      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
00972                      IE = ITAU
00973                      ITAUQ = IE + N
00974                      ITAUP = ITAUQ + N
00975                      IWORK = ITAUP + N
00976 *
00977 *                    Bidiagonalize R in WORK(IR)
00978 *                    (Workspace: need N*N+4*N, prefer N*N+3*N+2*N*NB)
00979 *
00980                      CALL SGEBRD( N, N, WORK( IR ), LDWRKR, S,
00981      $                            WORK( IE ), WORK( ITAUQ ),
00982      $                            WORK( ITAUP ), WORK( IWORK ),
00983      $                            LWORK-IWORK+1, IERR )
00984 *
00985 *                    Generate left vectors bidiagonalizing R in WORK(IR)
00986 *                    (Workspace: need N*N+4*N, prefer N*N+3*N+N*NB)
00987 *
00988                      CALL SORGBR( 'Q', N, N, N, WORK( IR ), LDWRKR,
00989      $                            WORK( ITAUQ ), WORK( IWORK ),
00990      $                            LWORK-IWORK+1, IERR )
00991                      IWORK = IE + N
00992 *
00993 *                    Perform bidiagonal QR iteration, computing left
00994 *                    singular vectors of R in WORK(IR)
00995 *                    (Workspace: need N*N+BDSPAC)
00996 *
00997                      CALL SBDSQR( 'U', N, 0, N, 0, S, WORK( IE ), DUM,
00998      $                            1, WORK( IR ), LDWRKR, DUM, 1,
00999      $                            WORK( IWORK ), INFO )
01000 *
01001 *                    Multiply Q in A by left singular vectors of R in
01002 *                    WORK(IR), storing result in U
01003 *                    (Workspace: need N*N)
01004 *
01005                      CALL SGEMM( 'N', 'N', M, N, N, ONE, A, LDA,
01006      $                           WORK( IR ), LDWRKR, ZERO, U, LDU )
01007 *
01008                   ELSE
01009 *
01010 *                    Insufficient workspace for a fast algorithm
01011 *
01012                      ITAU = 1
01013                      IWORK = ITAU + N
01014 *
01015 *                    Compute A=Q*R, copying result to U
01016 *                    (Workspace: need 2*N, prefer N+N*NB)
01017 *
01018                      CALL SGEQRF( M, N, A, LDA, WORK( ITAU ),
01019      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
01020                      CALL SLACPY( 'L', M, N, A, LDA, U, LDU )
01021 *
01022 *                    Generate Q in U
01023 *                    (Workspace: need 2*N, prefer N+N*NB)
01024 *
01025                      CALL SORGQR( M, N, N, U, LDU, WORK( ITAU ),
01026      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
01027                      IE = ITAU
01028                      ITAUQ = IE + N
01029                      ITAUP = ITAUQ + N
01030                      IWORK = ITAUP + N
01031 *
01032 *                    Zero out below R in A
01033 *
01034                      CALL SLASET( 'L', N-1, N-1, ZERO, ZERO, A( 2, 1 ),
01035      $                            LDA )
01036 *
01037 *                    Bidiagonalize R in A
01038 *                    (Workspace: need 4*N, prefer 3*N+2*N*NB)
01039 *
01040                      CALL SGEBRD( N, N, A, LDA, S, WORK( IE ),
01041      $                            WORK( ITAUQ ), WORK( ITAUP ),
01042      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
01043 *
01044 *                    Multiply Q in U by left vectors bidiagonalizing R
01045 *                    (Workspace: need 3*N+M, prefer 3*N+M*NB)
01046 *
01047                      CALL SORMBR( 'Q', 'R', 'N', M, N, N, A, LDA,
01048      $                            WORK( ITAUQ ), U, LDU, WORK( IWORK ),
01049      $                            LWORK-IWORK+1, IERR )
01050                      IWORK = IE + N
01051 *
01052 *                    Perform bidiagonal QR iteration, computing left
01053 *                    singular vectors of A in U
01054 *                    (Workspace: need BDSPAC)
01055 *
01056                      CALL SBDSQR( 'U', N, 0, M, 0, S, WORK( IE ), DUM,
01057      $                            1, U, LDU, DUM, 1, WORK( IWORK ),
01058      $                            INFO )
01059 *
01060                   END IF
01061 *
01062                ELSE IF( WNTVO ) THEN
01063 *
01064 *                 Path 5 (M much larger than N, JOBU='S', JOBVT='O')
01065 *                 N left singular vectors to be computed in U and
01066 *                 N right singular vectors to be overwritten on A
01067 *
01068                   IF( LWORK.GE.2*N*N+MAX( 4*N, BDSPAC ) ) THEN
01069 *
01070 *                    Sufficient workspace for a fast algorithm
01071 *
01072                      IU = 1
01073                      IF( LWORK.GE.WRKBL+2*LDA*N ) THEN
01074 *
01075 *                       WORK(IU) is LDA by N and WORK(IR) is LDA by N
01076 *
01077                         LDWRKU = LDA
01078                         IR = IU + LDWRKU*N
01079                         LDWRKR = LDA
01080                      ELSE IF( LWORK.GE.WRKBL+( LDA+N )*N ) THEN
01081 *
01082 *                       WORK(IU) is LDA by N and WORK(IR) is N by N
01083 *
01084                         LDWRKU = LDA
01085                         IR = IU + LDWRKU*N
01086                         LDWRKR = N
01087                      ELSE
01088 *
01089 *                       WORK(IU) is N by N and WORK(IR) is N by N
01090 *
01091                         LDWRKU = N
01092                         IR = IU + LDWRKU*N
01093                         LDWRKR = N
01094                      END IF
01095                      ITAU = IR + LDWRKR*N
01096                      IWORK = ITAU + N
01097 *
01098 *                    Compute A=Q*R
01099 *                    (Workspace: need 2*N*N+2*N, prefer 2*N*N+N+N*NB)
01100 *
01101                      CALL SGEQRF( M, N, A, LDA, WORK( ITAU ),
01102      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
01103 *
01104 *                    Copy R to WORK(IU), zeroing out below it
01105 *
01106                      CALL SLACPY( 'U', N, N, A, LDA, WORK( IU ),
01107      $                            LDWRKU )
01108                      CALL SLASET( 'L', N-1, N-1, ZERO, ZERO,
01109      $                            WORK( IU+1 ), LDWRKU )
01110 *
01111 *                    Generate Q in A
01112 *                    (Workspace: need 2*N*N+2*N, prefer 2*N*N+N+N*NB)
01113 *
01114                      CALL SORGQR( M, N, N, A, LDA, WORK( ITAU ),
01115      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
01116                      IE = ITAU
01117                      ITAUQ = IE + N
01118                      ITAUP = ITAUQ + N
01119                      IWORK = ITAUP + N
01120 *
01121 *                    Bidiagonalize R in WORK(IU), copying result to
01122 *                    WORK(IR)
01123 *                    (Workspace: need 2*N*N+4*N,
01124 *                                prefer 2*N*N+3*N+2*N*NB)
01125 *
01126                      CALL SGEBRD( N, N, WORK( IU ), LDWRKU, S,
01127      $                            WORK( IE ), WORK( ITAUQ ),
01128      $                            WORK( ITAUP ), WORK( IWORK ),
01129      $                            LWORK-IWORK+1, IERR )
01130                      CALL SLACPY( 'U', N, N, WORK( IU ), LDWRKU,
01131      $                            WORK( IR ), LDWRKR )
01132 *
01133 *                    Generate left bidiagonalizing vectors in WORK(IU)
01134 *                    (Workspace: need 2*N*N+4*N, prefer 2*N*N+3*N+N*NB)
01135 *
01136                      CALL SORGBR( 'Q', N, N, N, WORK( IU ), LDWRKU,
01137      $                            WORK( ITAUQ ), WORK( IWORK ),
01138      $                            LWORK-IWORK+1, IERR )
01139 *
01140 *                    Generate right bidiagonalizing vectors in WORK(IR)
01141 *                    (Workspace: need 2*N*N+4*N-1,
01142 *                                prefer 2*N*N+3*N+(N-1)*NB)
01143 *
01144                      CALL SORGBR( 'P', N, N, N, WORK( IR ), LDWRKR,
01145      $                            WORK( ITAUP ), WORK( IWORK ),
01146      $                            LWORK-IWORK+1, IERR )
01147                      IWORK = IE + N
01148 *
01149 *                    Perform bidiagonal QR iteration, computing left
01150 *                    singular vectors of R in WORK(IU) and computing
01151 *                    right singular vectors of R in WORK(IR)
01152 *                    (Workspace: need 2*N*N+BDSPAC)
01153 *
01154                      CALL SBDSQR( 'U', N, N, N, 0, S, WORK( IE ),
01155      $                            WORK( IR ), LDWRKR, WORK( IU ),
01156      $                            LDWRKU, DUM, 1, WORK( IWORK ), INFO )
01157 *
01158 *                    Multiply Q in A by left singular vectors of R in
01159 *                    WORK(IU), storing result in U
01160 *                    (Workspace: need N*N)
01161 *
01162                      CALL SGEMM( 'N', 'N', M, N, N, ONE, A, LDA,
01163      $                           WORK( IU ), LDWRKU, ZERO, U, LDU )
01164 *
01165 *                    Copy right singular vectors of R to A
01166 *                    (Workspace: need N*N)
01167 *
01168                      CALL SLACPY( 'F', N, N, WORK( IR ), LDWRKR, A,
01169      $                            LDA )
01170 *
01171                   ELSE
01172 *
01173 *                    Insufficient workspace for a fast algorithm
01174 *
01175                      ITAU = 1
01176                      IWORK = ITAU + N
01177 *
01178 *                    Compute A=Q*R, copying result to U
01179 *                    (Workspace: need 2*N, prefer N+N*NB)
01180 *
01181                      CALL SGEQRF( M, N, A, LDA, WORK( ITAU ),
01182      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
01183                      CALL SLACPY( 'L', M, N, A, LDA, U, LDU )
01184 *
01185 *                    Generate Q in U
01186 *                    (Workspace: need 2*N, prefer N+N*NB)
01187 *
01188                      CALL SORGQR( M, N, N, U, LDU, WORK( ITAU ),
01189      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
01190                      IE = ITAU
01191                      ITAUQ = IE + N
01192                      ITAUP = ITAUQ + N
01193                      IWORK = ITAUP + N
01194 *
01195 *                    Zero out below R in A
01196 *
01197                      CALL SLASET( 'L', N-1, N-1, ZERO, ZERO, A( 2, 1 ),
01198      $                            LDA )
01199 *
01200 *                    Bidiagonalize R in A
01201 *                    (Workspace: need 4*N, prefer 3*N+2*N*NB)
01202 *
01203                      CALL SGEBRD( N, N, A, LDA, S, WORK( IE ),
01204      $                            WORK( ITAUQ ), WORK( ITAUP ),
01205      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
01206 *
01207 *                    Multiply Q in U by left vectors bidiagonalizing R
01208 *                    (Workspace: need 3*N+M, prefer 3*N+M*NB)
01209 *
01210                      CALL SORMBR( 'Q', 'R', 'N', M, N, N, A, LDA,
01211      $                            WORK( ITAUQ ), U, LDU, WORK( IWORK ),
01212      $                            LWORK-IWORK+1, IERR )
01213 *
01214 *                    Generate right vectors bidiagonalizing R in A
01215 *                    (Workspace: need 4*N-1, prefer 3*N+(N-1)*NB)
01216 *
01217                      CALL SORGBR( 'P', N, N, N, A, LDA, WORK( ITAUP ),
01218      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
01219                      IWORK = IE + N
01220 *
01221 *                    Perform bidiagonal QR iteration, computing left
01222 *                    singular vectors of A in U and computing right
01223 *                    singular vectors of A in A
01224 *                    (Workspace: need BDSPAC)
01225 *
01226                      CALL SBDSQR( 'U', N, N, M, 0, S, WORK( IE ), A,
01227      $                            LDA, U, LDU, DUM, 1, WORK( IWORK ),
01228      $                            INFO )
01229 *
01230                   END IF
01231 *
01232                ELSE IF( WNTVAS ) THEN
01233 *
01234 *                 Path 6 (M much larger than N, JOBU='S', JOBVT='S'
01235 *                         or 'A')
01236 *                 N left singular vectors to be computed in U and
01237 *                 N right singular vectors to be computed in VT
01238 *
01239                   IF( LWORK.GE.N*N+MAX( 4*N, BDSPAC ) ) THEN
01240 *
01241 *                    Sufficient workspace for a fast algorithm
01242 *
01243                      IU = 1
01244                      IF( LWORK.GE.WRKBL+LDA*N ) THEN
01245 *
01246 *                       WORK(IU) is LDA by N
01247 *
01248                         LDWRKU = LDA
01249                      ELSE
01250 *
01251 *                       WORK(IU) is N by N
01252 *
01253                         LDWRKU = N
01254                      END IF
01255                      ITAU = IU + LDWRKU*N
01256                      IWORK = ITAU + N
01257 *
01258 *                    Compute A=Q*R
01259 *                    (Workspace: need N*N+2*N, prefer N*N+N+N*NB)
01260 *
01261                      CALL SGEQRF( M, N, A, LDA, WORK( ITAU ),
01262      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
01263 *
01264 *                    Copy R to WORK(IU), zeroing out below it
01265 *
01266                      CALL SLACPY( 'U', N, N, A, LDA, WORK( IU ),
01267      $                            LDWRKU )
01268                      CALL SLASET( 'L', N-1, N-1, ZERO, ZERO,
01269      $                            WORK( IU+1 ), LDWRKU )
01270 *
01271 *                    Generate Q in A
01272 *                    (Workspace: need N*N+2*N, prefer N*N+N+N*NB)
01273 *
01274                      CALL SORGQR( M, N, N, A, LDA, WORK( ITAU ),
01275      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
01276                      IE = ITAU
01277                      ITAUQ = IE + N
01278                      ITAUP = ITAUQ + N
01279                      IWORK = ITAUP + N
01280 *
01281 *                    Bidiagonalize R in WORK(IU), copying result to VT
01282 *                    (Workspace: need N*N+4*N, prefer N*N+3*N+2*N*NB)
01283 *
01284                      CALL SGEBRD( N, N, WORK( IU ), LDWRKU, S,
01285      $                            WORK( IE ), WORK( ITAUQ ),
01286      $                            WORK( ITAUP ), WORK( IWORK ),
01287      $                            LWORK-IWORK+1, IERR )
01288                      CALL SLACPY( 'U', N, N, WORK( IU ), LDWRKU, VT,
01289      $                            LDVT )
01290 *
01291 *                    Generate left bidiagonalizing vectors in WORK(IU)
01292 *                    (Workspace: need N*N+4*N, prefer N*N+3*N+N*NB)
01293 *
01294                      CALL SORGBR( 'Q', N, N, N, WORK( IU ), LDWRKU,
01295      $                            WORK( ITAUQ ), WORK( IWORK ),
01296      $                            LWORK-IWORK+1, IERR )
01297 *
01298 *                    Generate right bidiagonalizing vectors in VT
01299 *                    (Workspace: need N*N+4*N-1,
01300 *                                prefer N*N+3*N+(N-1)*NB)
01301 *
01302                      CALL SORGBR( 'P', N, N, N, VT, LDVT, WORK( ITAUP ),
01303      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
01304                      IWORK = IE + N
01305 *
01306 *                    Perform bidiagonal QR iteration, computing left
01307 *                    singular vectors of R in WORK(IU) and computing
01308 *                    right singular vectors of R in VT
01309 *                    (Workspace: need N*N+BDSPAC)
01310 *
01311                      CALL SBDSQR( 'U', N, N, N, 0, S, WORK( IE ), VT,
01312      $                            LDVT, WORK( IU ), LDWRKU, DUM, 1,
01313      $                            WORK( IWORK ), INFO )
01314 *
01315 *                    Multiply Q in A by left singular vectors of R in
01316 *                    WORK(IU), storing result in U
01317 *                    (Workspace: need N*N)
01318 *
01319                      CALL SGEMM( 'N', 'N', M, N, N, ONE, A, LDA,
01320      $                           WORK( IU ), LDWRKU, ZERO, U, LDU )
01321 *
01322                   ELSE
01323 *
01324 *                    Insufficient workspace for a fast algorithm
01325 *
01326                      ITAU = 1
01327                      IWORK = ITAU + N
01328 *
01329 *                    Compute A=Q*R, copying result to U
01330 *                    (Workspace: need 2*N, prefer N+N*NB)
01331 *
01332                      CALL SGEQRF( M, N, A, LDA, WORK( ITAU ),
01333      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
01334                      CALL SLACPY( 'L', M, N, A, LDA, U, LDU )
01335 *
01336 *                    Generate Q in U
01337 *                    (Workspace: need 2*N, prefer N+N*NB)
01338 *
01339                      CALL SORGQR( M, N, N, U, LDU, WORK( ITAU ),
01340      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
01341 *
01342 *                    Copy R to VT, zeroing out below it
01343 *
01344                      CALL SLACPY( 'U', N, N, A, LDA, VT, LDVT )
01345                      IF( N.GT.1 )
01346      $                  CALL SLASET( 'L', N-1, N-1, ZERO, ZERO,
01347      $                               VT( 2, 1 ), LDVT )
01348                      IE = ITAU
01349                      ITAUQ = IE + N
01350                      ITAUP = ITAUQ + N
01351                      IWORK = ITAUP + N
01352 *
01353 *                    Bidiagonalize R in VT
01354 *                    (Workspace: need 4*N, prefer 3*N+2*N*NB)
01355 *
01356                      CALL SGEBRD( N, N, VT, LDVT, S, WORK( IE ),
01357      $                            WORK( ITAUQ ), WORK( ITAUP ),
01358      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
01359 *
01360 *                    Multiply Q in U by left bidiagonalizing vectors
01361 *                    in VT
01362 *                    (Workspace: need 3*N+M, prefer 3*N+M*NB)
01363 *
01364                      CALL SORMBR( 'Q', 'R', 'N', M, N, N, VT, LDVT,
01365      $                            WORK( ITAUQ ), U, LDU, WORK( IWORK ),
01366      $                            LWORK-IWORK+1, IERR )
01367 *
01368 *                    Generate right bidiagonalizing vectors in VT
01369 *                    (Workspace: need 4*N-1, prefer 3*N+(N-1)*NB)
01370 *
01371                      CALL SORGBR( 'P', N, N, N, VT, LDVT, WORK( ITAUP ),
01372      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
01373                      IWORK = IE + N
01374 *
01375 *                    Perform bidiagonal QR iteration, computing left
01376 *                    singular vectors of A in U and computing right
01377 *                    singular vectors of A in VT
01378 *                    (Workspace: need BDSPAC)
01379 *
01380                      CALL SBDSQR( 'U', N, N, M, 0, S, WORK( IE ), VT,
01381      $                            LDVT, U, LDU, DUM, 1, WORK( IWORK ),
01382      $                            INFO )
01383 *
01384                   END IF
01385 *
01386                END IF
01387 *
01388             ELSE IF( WNTUA ) THEN
01389 *
01390                IF( WNTVN ) THEN
01391 *
01392 *                 Path 7 (M much larger than N, JOBU='A', JOBVT='N')
01393 *                 M left singular vectors to be computed in U and
01394 *                 no right singular vectors to be computed
01395 *
01396                   IF( LWORK.GE.N*N+MAX( N+M, 4*N, BDSPAC ) ) THEN
01397 *
01398 *                    Sufficient workspace for a fast algorithm
01399 *
01400                      IR = 1
01401                      IF( LWORK.GE.WRKBL+LDA*N ) THEN
01402 *
01403 *                       WORK(IR) is LDA by N
01404 *
01405                         LDWRKR = LDA
01406                      ELSE
01407 *
01408 *                       WORK(IR) is N by N
01409 *
01410                         LDWRKR = N
01411                      END IF
01412                      ITAU = IR + LDWRKR*N
01413                      IWORK = ITAU + N
01414 *
01415 *                    Compute A=Q*R, copying result to U
01416 *                    (Workspace: need N*N+2*N, prefer N*N+N+N*NB)
01417 *
01418                      CALL SGEQRF( M, N, A, LDA, WORK( ITAU ),
01419      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
01420                      CALL SLACPY( 'L', M, N, A, LDA, U, LDU )
01421 *
01422 *                    Copy R to WORK(IR), zeroing out below it
01423 *
01424                      CALL SLACPY( 'U', N, N, A, LDA, WORK( IR ),
01425      $                            LDWRKR )
01426                      CALL SLASET( 'L', N-1, N-1, ZERO, ZERO,
01427      $                            WORK( IR+1 ), LDWRKR )
01428 *
01429 *                    Generate Q in U
01430 *                    (Workspace: need N*N+N+M, prefer N*N+N+M*NB)
01431 *
01432                      CALL SORGQR( M, M, N, U, LDU, WORK( ITAU ),
01433      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
01434                      IE = ITAU
01435                      ITAUQ = IE + N
01436                      ITAUP = ITAUQ + N
01437                      IWORK = ITAUP + N
01438 *
01439 *                    Bidiagonalize R in WORK(IR)
01440 *                    (Workspace: need N*N+4*N, prefer N*N+3*N+2*N*NB)
01441 *
01442                      CALL SGEBRD( N, N, WORK( IR ), LDWRKR, S,
01443      $                            WORK( IE ), WORK( ITAUQ ),
01444      $                            WORK( ITAUP ), WORK( IWORK ),
01445      $                            LWORK-IWORK+1, IERR )
01446 *
01447 *                    Generate left bidiagonalizing vectors in WORK(IR)
01448 *                    (Workspace: need N*N+4*N, prefer N*N+3*N+N*NB)
01449 *
01450                      CALL SORGBR( 'Q', N, N, N, WORK( IR ), LDWRKR,
01451      $                            WORK( ITAUQ ), WORK( IWORK ),
01452      $                            LWORK-IWORK+1, IERR )
01453                      IWORK = IE + N
01454 *
01455 *                    Perform bidiagonal QR iteration, computing left
01456 *                    singular vectors of R in WORK(IR)
01457 *                    (Workspace: need N*N+BDSPAC)
01458 *
01459                      CALL SBDSQR( 'U', N, 0, N, 0, S, WORK( IE ), DUM,
01460      $                            1, WORK( IR ), LDWRKR, DUM, 1,
01461      $                            WORK( IWORK ), INFO )
01462 *
01463 *                    Multiply Q in U by left singular vectors of R in
01464 *                    WORK(IR), storing result in A
01465 *                    (Workspace: need N*N)
01466 *
01467                      CALL SGEMM( 'N', 'N', M, N, N, ONE, U, LDU,
01468      $                           WORK( IR ), LDWRKR, ZERO, A, LDA )
01469 *
01470 *                    Copy left singular vectors of A from A to U
01471 *
01472                      CALL SLACPY( 'F', M, N, A, LDA, U, LDU )
01473 *
01474                   ELSE
01475 *
01476 *                    Insufficient workspace for a fast algorithm
01477 *
01478                      ITAU = 1
01479                      IWORK = ITAU + N
01480 *
01481 *                    Compute A=Q*R, copying result to U
01482 *                    (Workspace: need 2*N, prefer N+N*NB)
01483 *
01484                      CALL SGEQRF( M, N, A, LDA, WORK( ITAU ),
01485      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
01486                      CALL SLACPY( 'L', M, N, A, LDA, U, LDU )
01487 *
01488 *                    Generate Q in U
01489 *                    (Workspace: need N+M, prefer N+M*NB)
01490 *
01491                      CALL SORGQR( M, M, N, U, LDU, WORK( ITAU ),
01492      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
01493                      IE = ITAU
01494                      ITAUQ = IE + N
01495                      ITAUP = ITAUQ + N
01496                      IWORK = ITAUP + N
01497 *
01498 *                    Zero out below R in A
01499 *
01500                      CALL SLASET( 'L', N-1, N-1, ZERO, ZERO, A( 2, 1 ),
01501      $                            LDA )
01502 *
01503 *                    Bidiagonalize R in A
01504 *                    (Workspace: need 4*N, prefer 3*N+2*N*NB)
01505 *
01506                      CALL SGEBRD( N, N, A, LDA, S, WORK( IE ),
01507      $                            WORK( ITAUQ ), WORK( ITAUP ),
01508      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
01509 *
01510 *                    Multiply Q in U by left bidiagonalizing vectors
01511 *                    in A
01512 *                    (Workspace: need 3*N+M, prefer 3*N+M*NB)
01513 *
01514                      CALL SORMBR( 'Q', 'R', 'N', M, N, N, A, LDA,
01515      $                            WORK( ITAUQ ), U, LDU, WORK( IWORK ),
01516      $                            LWORK-IWORK+1, IERR )
01517                      IWORK = IE + N
01518 *
01519 *                    Perform bidiagonal QR iteration, computing left
01520 *                    singular vectors of A in U
01521 *                    (Workspace: need BDSPAC)
01522 *
01523                      CALL SBDSQR( 'U', N, 0, M, 0, S, WORK( IE ), DUM,
01524      $                            1, U, LDU, DUM, 1, WORK( IWORK ),
01525      $                            INFO )
01526 *
01527                   END IF
01528 *
01529                ELSE IF( WNTVO ) THEN
01530 *
01531 *                 Path 8 (M much larger than N, JOBU='A', JOBVT='O')
01532 *                 M left singular vectors to be computed in U and
01533 *                 N right singular vectors to be overwritten on A
01534 *
01535                   IF( LWORK.GE.2*N*N+MAX( N+M, 4*N, BDSPAC ) ) THEN
01536 *
01537 *                    Sufficient workspace for a fast algorithm
01538 *
01539                      IU = 1
01540                      IF( LWORK.GE.WRKBL+2*LDA*N ) THEN
01541 *
01542 *                       WORK(IU) is LDA by N and WORK(IR) is LDA by N
01543 *
01544                         LDWRKU = LDA
01545                         IR = IU + LDWRKU*N
01546                         LDWRKR = LDA
01547                      ELSE IF( LWORK.GE.WRKBL+( LDA+N )*N ) THEN
01548 *
01549 *                       WORK(IU) is LDA by N and WORK(IR) is N by N
01550 *
01551                         LDWRKU = LDA
01552                         IR = IU + LDWRKU*N
01553                         LDWRKR = N
01554                      ELSE
01555 *
01556 *                       WORK(IU) is N by N and WORK(IR) is N by N
01557 *
01558                         LDWRKU = N
01559                         IR = IU + LDWRKU*N
01560                         LDWRKR = N
01561                      END IF
01562                      ITAU = IR + LDWRKR*N
01563                      IWORK = ITAU + N
01564 *
01565 *                    Compute A=Q*R, copying result to U
01566 *                    (Workspace: need 2*N*N+2*N, prefer 2*N*N+N+N*NB)
01567 *
01568                      CALL SGEQRF( M, N, A, LDA, WORK( ITAU ),
01569      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
01570                      CALL SLACPY( 'L', M, N, A, LDA, U, LDU )
01571 *
01572 *                    Generate Q in U
01573 *                    (Workspace: need 2*N*N+N+M, prefer 2*N*N+N+M*NB)
01574 *
01575                      CALL SORGQR( M, M, N, U, LDU, WORK( ITAU ),
01576      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
01577 *
01578 *                    Copy R to WORK(IU), zeroing out below it
01579 *
01580                      CALL SLACPY( 'U', N, N, A, LDA, WORK( IU ),
01581      $                            LDWRKU )
01582                      CALL SLASET( 'L', N-1, N-1, ZERO, ZERO,
01583      $                            WORK( IU+1 ), LDWRKU )
01584                      IE = ITAU
01585                      ITAUQ = IE + N
01586                      ITAUP = ITAUQ + N
01587                      IWORK = ITAUP + N
01588 *
01589 *                    Bidiagonalize R in WORK(IU), copying result to
01590 *                    WORK(IR)
01591 *                    (Workspace: need 2*N*N+4*N,
01592 *                                prefer 2*N*N+3*N+2*N*NB)
01593 *
01594                      CALL SGEBRD( N, N, WORK( IU ), LDWRKU, S,
01595      $                            WORK( IE ), WORK( ITAUQ ),
01596      $                            WORK( ITAUP ), WORK( IWORK ),
01597      $                            LWORK-IWORK+1, IERR )
01598                      CALL SLACPY( 'U', N, N, WORK( IU ), LDWRKU,
01599      $                            WORK( IR ), LDWRKR )
01600 *
01601 *                    Generate left bidiagonalizing vectors in WORK(IU)
01602 *                    (Workspace: need 2*N*N+4*N, prefer 2*N*N+3*N+N*NB)
01603 *
01604                      CALL SORGBR( 'Q', N, N, N, WORK( IU ), LDWRKU,
01605      $                            WORK( ITAUQ ), WORK( IWORK ),
01606      $                            LWORK-IWORK+1, IERR )
01607 *
01608 *                    Generate right bidiagonalizing vectors in WORK(IR)
01609 *                    (Workspace: need 2*N*N+4*N-1,
01610 *                                prefer 2*N*N+3*N+(N-1)*NB)
01611 *
01612                      CALL SORGBR( 'P', N, N, N, WORK( IR ), LDWRKR,
01613      $                            WORK( ITAUP ), WORK( IWORK ),
01614      $                            LWORK-IWORK+1, IERR )
01615                      IWORK = IE + N
01616 *
01617 *                    Perform bidiagonal QR iteration, computing left
01618 *                    singular vectors of R in WORK(IU) and computing
01619 *                    right singular vectors of R in WORK(IR)
01620 *                    (Workspace: need 2*N*N+BDSPAC)
01621 *
01622                      CALL SBDSQR( 'U', N, N, N, 0, S, WORK( IE ),
01623      $                            WORK( IR ), LDWRKR, WORK( IU ),
01624      $                            LDWRKU, DUM, 1, WORK( IWORK ), INFO )
01625 *
01626 *                    Multiply Q in U by left singular vectors of R in
01627 *                    WORK(IU), storing result in A
01628 *                    (Workspace: need N*N)
01629 *
01630                      CALL SGEMM( 'N', 'N', M, N, N, ONE, U, LDU,
01631      $                           WORK( IU ), LDWRKU, ZERO, A, LDA )
01632 *
01633 *                    Copy left singular vectors of A from A to U
01634 *
01635                      CALL SLACPY( 'F', M, N, A, LDA, U, LDU )
01636 *
01637 *                    Copy right singular vectors of R from WORK(IR) to A
01638 *
01639                      CALL SLACPY( 'F', N, N, WORK( IR ), LDWRKR, A,
01640      $                            LDA )
01641 *
01642                   ELSE
01643 *
01644 *                    Insufficient workspace for a fast algorithm
01645 *
01646                      ITAU = 1
01647                      IWORK = ITAU + N
01648 *
01649 *                    Compute A=Q*R, copying result to U
01650 *                    (Workspace: need 2*N, prefer N+N*NB)
01651 *
01652                      CALL SGEQRF( M, N, A, LDA, WORK( ITAU ),
01653      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
01654                      CALL SLACPY( 'L', M, N, A, LDA, U, LDU )
01655 *
01656 *                    Generate Q in U
01657 *                    (Workspace: need N+M, prefer N+M*NB)
01658 *
01659                      CALL SORGQR( M, M, N, U, LDU, WORK( ITAU ),
01660      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
01661                      IE = ITAU
01662                      ITAUQ = IE + N
01663                      ITAUP = ITAUQ + N
01664                      IWORK = ITAUP + N
01665 *
01666 *                    Zero out below R in A
01667 *
01668                      CALL SLASET( 'L', N-1, N-1, ZERO, ZERO, A( 2, 1 ),
01669      $                            LDA )
01670 *
01671 *                    Bidiagonalize R in A
01672 *                    (Workspace: need 4*N, prefer 3*N+2*N*NB)
01673 *
01674                      CALL SGEBRD( N, N, A, LDA, S, WORK( IE ),
01675      $                            WORK( ITAUQ ), WORK( ITAUP ),
01676      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
01677 *
01678 *                    Multiply Q in U by left bidiagonalizing vectors
01679 *                    in A
01680 *                    (Workspace: need 3*N+M, prefer 3*N+M*NB)
01681 *
01682                      CALL SORMBR( 'Q', 'R', 'N', M, N, N, A, LDA,
01683      $                            WORK( ITAUQ ), U, LDU, WORK( IWORK ),
01684      $                            LWORK-IWORK+1, IERR )
01685 *
01686 *                    Generate right bidiagonalizing vectors in A
01687 *                    (Workspace: need 4*N-1, prefer 3*N+(N-1)*NB)
01688 *
01689                      CALL SORGBR( 'P', N, N, N, A, LDA, WORK( ITAUP ),
01690      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
01691                      IWORK = IE + N
01692 *
01693 *                    Perform bidiagonal QR iteration, computing left
01694 *                    singular vectors of A in U and computing right
01695 *                    singular vectors of A in A
01696 *                    (Workspace: need BDSPAC)
01697 *
01698                      CALL SBDSQR( 'U', N, N, M, 0, S, WORK( IE ), A,
01699      $                            LDA, U, LDU, DUM, 1, WORK( IWORK ),
01700      $                            INFO )
01701 *
01702                   END IF
01703 *
01704                ELSE IF( WNTVAS ) THEN
01705 *
01706 *                 Path 9 (M much larger than N, JOBU='A', JOBVT='S'
01707 *                         or 'A')
01708 *                 M left singular vectors to be computed in U and
01709 *                 N right singular vectors to be computed in VT
01710 *
01711                   IF( LWORK.GE.N*N+MAX( N+M, 4*N, BDSPAC ) ) THEN
01712 *
01713 *                    Sufficient workspace for a fast algorithm
01714 *
01715                      IU = 1
01716                      IF( LWORK.GE.WRKBL+LDA*N ) THEN
01717 *
01718 *                       WORK(IU) is LDA by N
01719 *
01720                         LDWRKU = LDA
01721                      ELSE
01722 *
01723 *                       WORK(IU) is N by N
01724 *
01725                         LDWRKU = N
01726                      END IF
01727                      ITAU = IU + LDWRKU*N
01728                      IWORK = ITAU + N
01729 *
01730 *                    Compute A=Q*R, copying result to U
01731 *                    (Workspace: need N*N+2*N, prefer N*N+N+N*NB)
01732 *
01733                      CALL SGEQRF( M, N, A, LDA, WORK( ITAU ),
01734      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
01735                      CALL SLACPY( 'L', M, N, A, LDA, U, LDU )
01736 *
01737 *                    Generate Q in U
01738 *                    (Workspace: need N*N+N+M, prefer N*N+N+M*NB)
01739 *
01740                      CALL SORGQR( M, M, N, U, LDU, WORK( ITAU ),
01741      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
01742 *
01743 *                    Copy R to WORK(IU), zeroing out below it
01744 *
01745                      CALL SLACPY( 'U', N, N, A, LDA, WORK( IU ),
01746      $                            LDWRKU )
01747                      CALL SLASET( 'L', N-1, N-1, ZERO, ZERO,
01748      $                            WORK( IU+1 ), LDWRKU )
01749                      IE = ITAU
01750                      ITAUQ = IE + N
01751                      ITAUP = ITAUQ + N
01752                      IWORK = ITAUP + N
01753 *
01754 *                    Bidiagonalize R in WORK(IU), copying result to VT
01755 *                    (Workspace: need N*N+4*N, prefer N*N+3*N+2*N*NB)
01756 *
01757                      CALL SGEBRD( N, N, WORK( IU ), LDWRKU, S,
01758      $                            WORK( IE ), WORK( ITAUQ ),
01759      $                            WORK( ITAUP ), WORK( IWORK ),
01760      $                            LWORK-IWORK+1, IERR )
01761                      CALL SLACPY( 'U', N, N, WORK( IU ), LDWRKU, VT,
01762      $                            LDVT )
01763 *
01764 *                    Generate left bidiagonalizing vectors in WORK(IU)
01765 *                    (Workspace: need N*N+4*N, prefer N*N+3*N+N*NB)
01766 *
01767                      CALL SORGBR( 'Q', N, N, N, WORK( IU ), LDWRKU,
01768      $                            WORK( ITAUQ ), WORK( IWORK ),
01769      $                            LWORK-IWORK+1, IERR )
01770 *
01771 *                    Generate right bidiagonalizing vectors in VT
01772 *                    (Workspace: need N*N+4*N-1,
01773 *                                prefer N*N+3*N+(N-1)*NB)
01774 *
01775                      CALL SORGBR( 'P', N, N, N, VT, LDVT, WORK( ITAUP ),
01776      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
01777                      IWORK = IE + N
01778 *
01779 *                    Perform bidiagonal QR iteration, computing left
01780 *                    singular vectors of R in WORK(IU) and computing
01781 *                    right singular vectors of R in VT
01782 *                    (Workspace: need N*N+BDSPAC)
01783 *
01784                      CALL SBDSQR( 'U', N, N, N, 0, S, WORK( IE ), VT,
01785      $                            LDVT, WORK( IU ), LDWRKU, DUM, 1,
01786      $                            WORK( IWORK ), INFO )
01787 *
01788 *                    Multiply Q in U by left singular vectors of R in
01789 *                    WORK(IU), storing result in A
01790 *                    (Workspace: need N*N)
01791 *
01792                      CALL SGEMM( 'N', 'N', M, N, N, ONE, U, LDU,
01793      $                           WORK( IU ), LDWRKU, ZERO, A, LDA )
01794 *
01795 *                    Copy left singular vectors of A from A to U
01796 *
01797                      CALL SLACPY( 'F', M, N, A, LDA, U, LDU )
01798 *
01799                   ELSE
01800 *
01801 *                    Insufficient workspace for a fast algorithm
01802 *
01803                      ITAU = 1
01804                      IWORK = ITAU + N
01805 *
01806 *                    Compute A=Q*R, copying result to U
01807 *                    (Workspace: need 2*N, prefer N+N*NB)
01808 *
01809                      CALL SGEQRF( M, N, A, LDA, WORK( ITAU ),
01810      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
01811                      CALL SLACPY( 'L', M, N, A, LDA, U, LDU )
01812 *
01813 *                    Generate Q in U
01814 *                    (Workspace: need N+M, prefer N+M*NB)
01815 *
01816                      CALL SORGQR( M, M, N, U, LDU, WORK( ITAU ),
01817      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
01818 *
01819 *                    Copy R from A to VT, zeroing out below it
01820 *
01821                      CALL SLACPY( 'U', N, N, A, LDA, VT, LDVT )
01822                      IF( N.GT.1 )
01823      $                  CALL SLASET( 'L', N-1, N-1, ZERO, ZERO,
01824      $                               VT( 2, 1 ), LDVT )
01825                      IE = ITAU
01826                      ITAUQ = IE + N
01827                      ITAUP = ITAUQ + N
01828                      IWORK = ITAUP + N
01829 *
01830 *                    Bidiagonalize R in VT
01831 *                    (Workspace: need 4*N, prefer 3*N+2*N*NB)
01832 *
01833                      CALL SGEBRD( N, N, VT, LDVT, S, WORK( IE ),
01834      $                            WORK( ITAUQ ), WORK( ITAUP ),
01835      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
01836 *
01837 *                    Multiply Q in U by left bidiagonalizing vectors
01838 *                    in VT
01839 *                    (Workspace: need 3*N+M, prefer 3*N+M*NB)
01840 *
01841                      CALL SORMBR( 'Q', 'R', 'N', M, N, N, VT, LDVT,
01842      $                            WORK( ITAUQ ), U, LDU, WORK( IWORK ),
01843      $                            LWORK-IWORK+1, IERR )
01844 *
01845 *                    Generate right bidiagonalizing vectors in VT
01846 *                    (Workspace: need 4*N-1, prefer 3*N+(N-1)*NB)
01847 *
01848                      CALL SORGBR( 'P', N, N, N, VT, LDVT, WORK( ITAUP ),
01849      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
01850                      IWORK = IE + N
01851 *
01852 *                    Perform bidiagonal QR iteration, computing left
01853 *                    singular vectors of A in U and computing right
01854 *                    singular vectors of A in VT
01855 *                    (Workspace: need BDSPAC)
01856 *
01857                      CALL SBDSQR( 'U', N, N, M, 0, S, WORK( IE ), VT,
01858      $                            LDVT, U, LDU, DUM, 1, WORK( IWORK ),
01859      $                            INFO )
01860 *
01861                   END IF
01862 *
01863                END IF
01864 *
01865             END IF
01866 *
01867          ELSE
01868 *
01869 *           M .LT. MNTHR
01870 *
01871 *           Path 10 (M at least N, but not much larger)
01872 *           Reduce to bidiagonal form without QR decomposition
01873 *
01874             IE = 1
01875             ITAUQ = IE + N
01876             ITAUP = ITAUQ + N
01877             IWORK = ITAUP + N
01878 *
01879 *           Bidiagonalize A
01880 *           (Workspace: need 3*N+M, prefer 3*N+(M+N)*NB)
01881 *
01882             CALL SGEBRD( M, N, A, LDA, S, WORK( IE ), WORK( ITAUQ ),
01883      $                   WORK( ITAUP ), WORK( IWORK ), LWORK-IWORK+1,
01884      $                   IERR )
01885             IF( WNTUAS ) THEN
01886 *
01887 *              If left singular vectors desired in U, copy result to U
01888 *              and generate left bidiagonalizing vectors in U
01889 *              (Workspace: need 3*N+NCU, prefer 3*N+NCU*NB)
01890 *
01891                CALL SLACPY( 'L', M, N, A, LDA, U, LDU )
01892                IF( WNTUS )
01893      $            NCU = N
01894                IF( WNTUA )
01895      $            NCU = M
01896                CALL SORGBR( 'Q', M, NCU, N, U, LDU, WORK( ITAUQ ),
01897      $                      WORK( IWORK ), LWORK-IWORK+1, IERR )
01898             END IF
01899             IF( WNTVAS ) THEN
01900 *
01901 *              If right singular vectors desired in VT, copy result to
01902 *              VT and generate right bidiagonalizing vectors in VT
01903 *              (Workspace: need 4*N-1, prefer 3*N+(N-1)*NB)
01904 *
01905                CALL SLACPY( 'U', N, N, A, LDA, VT, LDVT )
01906                CALL SORGBR( 'P', N, N, N, VT, LDVT, WORK( ITAUP ),
01907      $                      WORK( IWORK ), LWORK-IWORK+1, IERR )
01908             END IF
01909             IF( WNTUO ) THEN
01910 *
01911 *              If left singular vectors desired in A, generate left
01912 *              bidiagonalizing vectors in A
01913 *              (Workspace: need 4*N, prefer 3*N+N*NB)
01914 *
01915                CALL SORGBR( 'Q', M, N, N, A, LDA, WORK( ITAUQ ),
01916      $                      WORK( IWORK ), LWORK-IWORK+1, IERR )
01917             END IF
01918             IF( WNTVO ) THEN
01919 *
01920 *              If right singular vectors desired in A, generate right
01921 *              bidiagonalizing vectors in A
01922 *              (Workspace: need 4*N-1, prefer 3*N+(N-1)*NB)
01923 *
01924                CALL SORGBR( 'P', N, N, N, A, LDA, WORK( ITAUP ),
01925      $                      WORK( IWORK ), LWORK-IWORK+1, IERR )
01926             END IF
01927             IWORK = IE + N
01928             IF( WNTUAS .OR. WNTUO )
01929      $         NRU = M
01930             IF( WNTUN )
01931      $         NRU = 0
01932             IF( WNTVAS .OR. WNTVO )
01933      $         NCVT = N
01934             IF( WNTVN )
01935      $         NCVT = 0
01936             IF( ( .NOT.WNTUO ) .AND. ( .NOT.WNTVO ) ) THEN
01937 *
01938 *              Perform bidiagonal QR iteration, if desired, computing
01939 *              left singular vectors in U and computing right singular
01940 *              vectors in VT
01941 *              (Workspace: need BDSPAC)
01942 *
01943                CALL SBDSQR( 'U', N, NCVT, NRU, 0, S, WORK( IE ), VT,
01944      $                      LDVT, U, LDU, DUM, 1, WORK( IWORK ), INFO )
01945             ELSE IF( ( .NOT.WNTUO ) .AND. WNTVO ) THEN
01946 *
01947 *              Perform bidiagonal QR iteration, if desired, computing
01948 *              left singular vectors in U and computing right singular
01949 *              vectors in A
01950 *              (Workspace: need BDSPAC)
01951 *
01952                CALL SBDSQR( 'U', N, NCVT, NRU, 0, S, WORK( IE ), A, LDA,
01953      $                      U, LDU, DUM, 1, WORK( IWORK ), INFO )
01954             ELSE
01955 *
01956 *              Perform bidiagonal QR iteration, if desired, computing
01957 *              left singular vectors in A and computing right singular
01958 *              vectors in VT
01959 *              (Workspace: need BDSPAC)
01960 *
01961                CALL SBDSQR( 'U', N, NCVT, NRU, 0, S, WORK( IE ), VT,
01962      $                      LDVT, A, LDA, DUM, 1, WORK( IWORK ), INFO )
01963             END IF
01964 *
01965          END IF
01966 *
01967       ELSE
01968 *
01969 *        A has more columns than rows. If A has sufficiently more
01970 *        columns than rows, first reduce using the LQ decomposition (if
01971 *        sufficient workspace available)
01972 *
01973          IF( N.GE.MNTHR ) THEN
01974 *
01975             IF( WNTVN ) THEN
01976 *
01977 *              Path 1t(N much larger than M, JOBVT='N')
01978 *              No right singular vectors to be computed
01979 *
01980                ITAU = 1
01981                IWORK = ITAU + M
01982 *
01983 *              Compute A=L*Q
01984 *              (Workspace: need 2*M, prefer M+M*NB)
01985 *
01986                CALL SGELQF( M, N, A, LDA, WORK( ITAU ), WORK( IWORK ),
01987      $                      LWORK-IWORK+1, IERR )
01988 *
01989 *              Zero out above L
01990 *
01991                CALL SLASET( 'U', M-1, M-1, ZERO, ZERO, A( 1, 2 ), LDA )
01992                IE = 1
01993                ITAUQ = IE + M
01994                ITAUP = ITAUQ + M
01995                IWORK = ITAUP + M
01996 *
01997 *              Bidiagonalize L in A
01998 *              (Workspace: need 4*M, prefer 3*M+2*M*NB)
01999 *
02000                CALL SGEBRD( M, M, A, LDA, S, WORK( IE ), WORK( ITAUQ ),
02001      $                      WORK( ITAUP ), WORK( IWORK ), LWORK-IWORK+1,
02002      $                      IERR )
02003                IF( WNTUO .OR. WNTUAS ) THEN
02004 *
02005 *                 If left singular vectors desired, generate Q
02006 *                 (Workspace: need 4*M, prefer 3*M+M*NB)
02007 *
02008                   CALL SORGBR( 'Q', M, M, M, A, LDA, WORK( ITAUQ ),
02009      $                         WORK( IWORK ), LWORK-IWORK+1, IERR )
02010                END IF
02011                IWORK = IE + M
02012                NRU = 0
02013                IF( WNTUO .OR. WNTUAS )
02014      $            NRU = M
02015 *
02016 *              Perform bidiagonal QR iteration, computing left singular
02017 *              vectors of A in A if desired
02018 *              (Workspace: need BDSPAC)
02019 *
02020                CALL SBDSQR( 'U', M, 0, NRU, 0, S, WORK( IE ), DUM, 1, A,
02021      $                      LDA, DUM, 1, WORK( IWORK ), INFO )
02022 *
02023 *              If left singular vectors desired in U, copy them there
02024 *
02025                IF( WNTUAS )
02026      $            CALL SLACPY( 'F', M, M, A, LDA, U, LDU )
02027 *
02028             ELSE IF( WNTVO .AND. WNTUN ) THEN
02029 *
02030 *              Path 2t(N much larger than M, JOBU='N', JOBVT='O')
02031 *              M right singular vectors to be overwritten on A and
02032 *              no left singular vectors to be computed
02033 *
02034                IF( LWORK.GE.M*M+MAX( 4*M, BDSPAC ) ) THEN
02035 *
02036 *                 Sufficient workspace for a fast algorithm
02037 *
02038                   IR = 1
02039                   IF( LWORK.GE.MAX( WRKBL, LDA*N+M )+LDA*M ) THEN
02040 *
02041 *                    WORK(IU) is LDA by N and WORK(IR) is LDA by M
02042 *
02043                      LDWRKU = LDA
02044                      CHUNK = N
02045                      LDWRKR = LDA
02046                   ELSE IF( LWORK.GE.MAX( WRKBL, LDA*N+M )+M*M ) THEN
02047 *
02048 *                    WORK(IU) is LDA by N and WORK(IR) is M by M
02049 *
02050                      LDWRKU = LDA
02051                      CHUNK = N
02052                      LDWRKR = M
02053                   ELSE
02054 *
02055 *                    WORK(IU) is M by CHUNK and WORK(IR) is M by M
02056 *
02057                      LDWRKU = M
02058                      CHUNK = ( LWORK-M*M-M ) / M
02059                      LDWRKR = M
02060                   END IF
02061                   ITAU = IR + LDWRKR*M
02062                   IWORK = ITAU + M
02063 *
02064 *                 Compute A=L*Q
02065 *                 (Workspace: need M*M+2*M, prefer M*M+M+M*NB)
02066 *
02067                   CALL SGELQF( M, N, A, LDA, WORK( ITAU ),
02068      $                         WORK( IWORK ), LWORK-IWORK+1, IERR )
02069 *
02070 *                 Copy L to WORK(IR) and zero out above it
02071 *
02072                   CALL SLACPY( 'L', M, M, A, LDA, WORK( IR ), LDWRKR )
02073                   CALL SLASET( 'U', M-1, M-1, ZERO, ZERO,
02074      $                         WORK( IR+LDWRKR ), LDWRKR )
02075 *
02076 *                 Generate Q in A
02077 *                 (Workspace: need M*M+2*M, prefer M*M+M+M*NB)
02078 *
02079                   CALL SORGLQ( M, N, M, A, LDA, WORK( ITAU ),
02080      $                         WORK( IWORK ), LWORK-IWORK+1, IERR )
02081                   IE = ITAU
02082                   ITAUQ = IE + M
02083                   ITAUP = ITAUQ + M
02084                   IWORK = ITAUP + M
02085 *
02086 *                 Bidiagonalize L in WORK(IR)
02087 *                 (Workspace: need M*M+4*M, prefer M*M+3*M+2*M*NB)
02088 *
02089                   CALL SGEBRD( M, M, WORK( IR ), LDWRKR, S, WORK( IE ),
02090      $                         WORK( ITAUQ ), WORK( ITAUP ),
02091      $                         WORK( IWORK ), LWORK-IWORK+1, IERR )
02092 *
02093 *                 Generate right vectors bidiagonalizing L
02094 *                 (Workspace: need M*M+4*M-1, prefer M*M+3*M+(M-1)*NB)
02095 *
02096                   CALL SORGBR( 'P', M, M, M, WORK( IR ), LDWRKR,
02097      $                         WORK( ITAUP ), WORK( IWORK ),
02098      $                         LWORK-IWORK+1, IERR )
02099                   IWORK = IE + M
02100 *
02101 *                 Perform bidiagonal QR iteration, computing right
02102 *                 singular vectors of L in WORK(IR)
02103 *                 (Workspace: need M*M+BDSPAC)
02104 *
02105                   CALL SBDSQR( 'U', M, M, 0, 0, S, WORK( IE ),
02106      $                         WORK( IR ), LDWRKR, DUM, 1, DUM, 1,
02107      $                         WORK( IWORK ), INFO )
02108                   IU = IE + M
02109 *
02110 *                 Multiply right singular vectors of L in WORK(IR) by Q
02111 *                 in A, storing result in WORK(IU) and copying to A
02112 *                 (Workspace: need M*M+2*M, prefer M*M+M*N+M)
02113 *
02114                   DO 30 I = 1, N, CHUNK
02115                      BLK = MIN( N-I+1, CHUNK )
02116                      CALL SGEMM( 'N', 'N', M, BLK, M, ONE, WORK( IR ),
02117      $                           LDWRKR, A( 1, I ), LDA, ZERO,
02118      $                           WORK( IU ), LDWRKU )
02119                      CALL SLACPY( 'F', M, BLK, WORK( IU ), LDWRKU,
02120      $                            A( 1, I ), LDA )
02121    30             CONTINUE
02122 *
02123                ELSE
02124 *
02125 *                 Insufficient workspace for a fast algorithm
02126 *
02127                   IE = 1
02128                   ITAUQ = IE + M
02129                   ITAUP = ITAUQ + M
02130                   IWORK = ITAUP + M
02131 *
02132 *                 Bidiagonalize A
02133 *                 (Workspace: need 3*M+N, prefer 3*M+(M+N)*NB)
02134 *
02135                   CALL SGEBRD( M, N, A, LDA, S, WORK( IE ),
02136      $                         WORK( ITAUQ ), WORK( ITAUP ),
02137      $                         WORK( IWORK ), LWORK-IWORK+1, IERR )
02138 *
02139 *                 Generate right vectors bidiagonalizing A
02140 *                 (Workspace: need 4*M, prefer 3*M+M*NB)
02141 *
02142                   CALL SORGBR( 'P', M, N, M, A, LDA, WORK( ITAUP ),
02143      $                         WORK( IWORK ), LWORK-IWORK+1, IERR )
02144                   IWORK = IE + M
02145 *
02146 *                 Perform bidiagonal QR iteration, computing right
02147 *                 singular vectors of A in A
02148 *                 (Workspace: need BDSPAC)
02149 *
02150                   CALL SBDSQR( 'L', M, N, 0, 0, S, WORK( IE ), A, LDA,
02151      $                         DUM, 1, DUM, 1, WORK( IWORK ), INFO )
02152 *
02153                END IF
02154 *
02155             ELSE IF( WNTVO .AND. WNTUAS ) THEN
02156 *
02157 *              Path 3t(N much larger than M, JOBU='S' or 'A', JOBVT='O')
02158 *              M right singular vectors to be overwritten on A and
02159 *              M left singular vectors to be computed in U
02160 *
02161                IF( LWORK.GE.M*M+MAX( 4*M, BDSPAC ) ) THEN
02162 *
02163 *                 Sufficient workspace for a fast algorithm
02164 *
02165                   IR = 1
02166                   IF( LWORK.GE.MAX( WRKBL, LDA*N+M )+LDA*M ) THEN
02167 *
02168 *                    WORK(IU) is LDA by N and WORK(IR) is LDA by M
02169 *
02170                      LDWRKU = LDA
02171                      CHUNK = N
02172                      LDWRKR = LDA
02173                   ELSE IF( LWORK.GE.MAX( WRKBL, LDA*N+M )+M*M ) THEN
02174 *
02175 *                    WORK(IU) is LDA by N and WORK(IR) is M by M
02176 *
02177                      LDWRKU = LDA
02178                      CHUNK = N
02179                      LDWRKR = M
02180                   ELSE
02181 *
02182 *                    WORK(IU) is M by CHUNK and WORK(IR) is M by M
02183 *
02184                      LDWRKU = M
02185                      CHUNK = ( LWORK-M*M-M ) / M
02186                      LDWRKR = M
02187                   END IF
02188                   ITAU = IR + LDWRKR*M
02189                   IWORK = ITAU + M
02190 *
02191 *                 Compute A=L*Q
02192 *                 (Workspace: need M*M+2*M, prefer M*M+M+M*NB)
02193 *
02194                   CALL SGELQF( M, N, A, LDA, WORK( ITAU ),
02195      $                         WORK( IWORK ), LWORK-IWORK+1, IERR )
02196 *
02197 *                 Copy L to U, zeroing about above it
02198 *
02199                   CALL SLACPY( 'L', M, M, A, LDA, U, LDU )
02200                   CALL SLASET( 'U', M-1, M-1, ZERO, ZERO, U( 1, 2 ),
02201      $                         LDU )
02202 *
02203 *                 Generate Q in A
02204 *                 (Workspace: need M*M+2*M, prefer M*M+M+M*NB)
02205 *
02206                   CALL SORGLQ( M, N, M, A, LDA, WORK( ITAU ),
02207      $                         WORK( IWORK ), LWORK-IWORK+1, IERR )
02208                   IE = ITAU
02209                   ITAUQ = IE + M
02210                   ITAUP = ITAUQ + M
02211                   IWORK = ITAUP + M
02212 *
02213 *                 Bidiagonalize L in U, copying result to WORK(IR)
02214 *                 (Workspace: need M*M+4*M, prefer M*M+3*M+2*M*NB)
02215 *
02216                   CALL SGEBRD( M, M, U, LDU, S, WORK( IE ),
02217      $                         WORK( ITAUQ ), WORK( ITAUP ),
02218      $                         WORK( IWORK ), LWORK-IWORK+1, IERR )
02219                   CALL SLACPY( 'U', M, M, U, LDU, WORK( IR ), LDWRKR )
02220 *
02221 *                 Generate right vectors bidiagonalizing L in WORK(IR)
02222 *                 (Workspace: need M*M+4*M-1, prefer M*M+3*M+(M-1)*NB)
02223 *
02224                   CALL SORGBR( 'P', M, M, M, WORK( IR ), LDWRKR,
02225      $                         WORK( ITAUP ), WORK( IWORK ),
02226      $                         LWORK-IWORK+1, IERR )
02227 *
02228 *                 Generate left vectors bidiagonalizing L in U
02229 *                 (Workspace: need M*M+4*M, prefer M*M+3*M+M*NB)
02230 *
02231                   CALL SORGBR( 'Q', M, M, M, U, LDU, WORK( ITAUQ ),
02232      $                         WORK( IWORK ), LWORK-IWORK+1, IERR )
02233                   IWORK = IE + M
02234 *
02235 *                 Perform bidiagonal QR iteration, computing left
02236 *                 singular vectors of L in U, and computing right
02237 *                 singular vectors of L in WORK(IR)
02238 *                 (Workspace: need M*M+BDSPAC)
02239 *
02240                   CALL SBDSQR( 'U', M, M, M, 0, S, WORK( IE ),
02241      $                         WORK( IR ), LDWRKR, U, LDU, DUM, 1,
02242      $                         WORK( IWORK ), INFO )
02243                   IU = IE + M
02244 *
02245 *                 Multiply right singular vectors of L in WORK(IR) by Q
02246 *                 in A, storing result in WORK(IU) and copying to A
02247 *                 (Workspace: need M*M+2*M, prefer M*M+M*N+M))
02248 *
02249                   DO 40 I = 1, N, CHUNK
02250                      BLK = MIN( N-I+1, CHUNK )
02251                      CALL SGEMM( 'N', 'N', M, BLK, M, ONE, WORK( IR ),
02252      $                           LDWRKR, A( 1, I ), LDA, ZERO,
02253      $                           WORK( IU ), LDWRKU )
02254                      CALL SLACPY( 'F', M, BLK, WORK( IU ), LDWRKU,
02255      $                            A( 1, I ), LDA )
02256    40             CONTINUE
02257 *
02258                ELSE
02259 *
02260 *                 Insufficient workspace for a fast algorithm
02261 *
02262                   ITAU = 1
02263                   IWORK = ITAU + M
02264 *
02265 *                 Compute A=L*Q
02266 *                 (Workspace: need 2*M, prefer M+M*NB)
02267 *
02268                   CALL SGELQF( M, N, A, LDA, WORK( ITAU ),
02269      $                         WORK( IWORK ), LWORK-IWORK+1, IERR )
02270 *
02271 *                 Copy L to U, zeroing out above it
02272 *
02273                   CALL SLACPY( 'L', M, M, A, LDA, U, LDU )
02274                   CALL SLASET( 'U', M-1, M-1, ZERO, ZERO, U( 1, 2 ),
02275      $                         LDU )
02276 *
02277 *                 Generate Q in A
02278 *                 (Workspace: need 2*M, prefer M+M*NB)
02279 *
02280                   CALL SORGLQ( M, N, M, A, LDA, WORK( ITAU ),
02281      $                         WORK( IWORK ), LWORK-IWORK+1, IERR )
02282                   IE = ITAU
02283                   ITAUQ = IE + M
02284                   ITAUP = ITAUQ + M
02285                   IWORK = ITAUP + M
02286 *
02287 *                 Bidiagonalize L in U
02288 *                 (Workspace: need 4*M, prefer 3*M+2*M*NB)
02289 *
02290                   CALL SGEBRD( M, M, U, LDU, S, WORK( IE ),
02291      $                         WORK( ITAUQ ), WORK( ITAUP ),
02292      $                         WORK( IWORK ), LWORK-IWORK+1, IERR )
02293 *
02294 *                 Multiply right vectors bidiagonalizing L by Q in A
02295 *                 (Workspace: need 3*M+N, prefer 3*M+N*NB)
02296 *
02297                   CALL SORMBR( 'P', 'L', 'T', M, N, M, U, LDU,
02298      $                         WORK( ITAUP ), A, LDA, WORK( IWORK ),
02299      $                         LWORK-IWORK+1, IERR )
02300 *
02301 *                 Generate left vectors bidiagonalizing L in U
02302 *                 (Workspace: need 4*M, prefer 3*M+M*NB)
02303 *
02304                   CALL SORGBR( 'Q', M, M, M, U, LDU, WORK( ITAUQ ),
02305      $                         WORK( IWORK ), LWORK-IWORK+1, IERR )
02306                   IWORK = IE + M
02307 *
02308 *                 Perform bidiagonal QR iteration, computing left
02309 *                 singular vectors of A in U and computing right
02310 *                 singular vectors of A in A
02311 *                 (Workspace: need BDSPAC)
02312 *
02313                   CALL SBDSQR( 'U', M, N, M, 0, S, WORK( IE ), A, LDA,
02314      $                         U, LDU, DUM, 1, WORK( IWORK ), INFO )
02315 *
02316                END IF
02317 *
02318             ELSE IF( WNTVS ) THEN
02319 *
02320                IF( WNTUN ) THEN
02321 *
02322 *                 Path 4t(N much larger than M, JOBU='N', JOBVT='S')
02323 *                 M right singular vectors to be computed in VT and
02324 *                 no left singular vectors to be computed
02325 *
02326                   IF( LWORK.GE.M*M+MAX( 4*M, BDSPAC ) ) THEN
02327 *
02328 *                    Sufficient workspace for a fast algorithm
02329 *
02330                      IR = 1
02331                      IF( LWORK.GE.WRKBL+LDA*M ) THEN
02332 *
02333 *                       WORK(IR) is LDA by M
02334 *
02335                         LDWRKR = LDA
02336                      ELSE
02337 *
02338 *                       WORK(IR) is M by M
02339 *
02340                         LDWRKR = M
02341                      END IF
02342                      ITAU = IR + LDWRKR*M
02343                      IWORK = ITAU + M
02344 *
02345 *                    Compute A=L*Q
02346 *                    (Workspace: need M*M+2*M, prefer M*M+M+M*NB)
02347 *
02348                      CALL SGELQF( M, N, A, LDA, WORK( ITAU ),
02349      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
02350 *
02351 *                    Copy L to WORK(IR), zeroing out above it
02352 *
02353                      CALL SLACPY( 'L', M, M, A, LDA, WORK( IR ),
02354      $                            LDWRKR )
02355                      CALL SLASET( 'U', M-1, M-1, ZERO, ZERO,
02356      $                            WORK( IR+LDWRKR ), LDWRKR )
02357 *
02358 *                    Generate Q in A
02359 *                    (Workspace: need M*M+2*M, prefer M*M+M+M*NB)
02360 *
02361                      CALL SORGLQ( M, N, M, A, LDA, WORK( ITAU ),
02362      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
02363                      IE = ITAU
02364                      ITAUQ = IE + M
02365                      ITAUP = ITAUQ + M
02366                      IWORK = ITAUP + M
02367 *
02368 *                    Bidiagonalize L in WORK(IR)
02369 *                    (Workspace: need M*M+4*M, prefer M*M+3*M+2*M*NB)
02370 *
02371                      CALL SGEBRD( M, M, WORK( IR ), LDWRKR, S,
02372      $                            WORK( IE ), WORK( ITAUQ ),
02373      $                            WORK( ITAUP ), WORK( IWORK ),
02374      $                            LWORK-IWORK+1, IERR )
02375 *
02376 *                    Generate right vectors bidiagonalizing L in
02377 *                    WORK(IR)
02378 *                    (Workspace: need M*M+4*M, prefer M*M+3*M+(M-1)*NB)
02379 *
02380                      CALL SORGBR( 'P', M, M, M, WORK( IR ), LDWRKR,
02381      $                            WORK( ITAUP ), WORK( IWORK ),
02382      $                            LWORK-IWORK+1, IERR )
02383                      IWORK = IE + M
02384 *
02385 *                    Perform bidiagonal QR iteration, computing right
02386 *                    singular vectors of L in WORK(IR)
02387 *                    (Workspace: need M*M+BDSPAC)
02388 *
02389                      CALL SBDSQR( 'U', M, M, 0, 0, S, WORK( IE ),
02390      $                            WORK( IR ), LDWRKR, DUM, 1, DUM, 1,
02391      $                            WORK( IWORK ), INFO )
02392 *
02393 *                    Multiply right singular vectors of L in WORK(IR) by
02394 *                    Q in A, storing result in VT
02395 *                    (Workspace: need M*M)
02396 *
02397                      CALL SGEMM( 'N', 'N', M, N, M, ONE, WORK( IR ),
02398      $                           LDWRKR, A, LDA, ZERO, VT, LDVT )
02399 *
02400                   ELSE
02401 *
02402 *                    Insufficient workspace for a fast algorithm
02403 *
02404                      ITAU = 1
02405                      IWORK = ITAU + M
02406 *
02407 *                    Compute A=L*Q
02408 *                    (Workspace: need 2*M, prefer M+M*NB)
02409 *
02410                      CALL SGELQF( M, N, A, LDA, WORK( ITAU ),
02411      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
02412 *
02413 *                    Copy result to VT
02414 *
02415                      CALL SLACPY( 'U', M, N, A, LDA, VT, LDVT )
02416 *
02417 *                    Generate Q in VT
02418 *                    (Workspace: need 2*M, prefer M+M*NB)
02419 *
02420                      CALL SORGLQ( M, N, M, VT, LDVT, WORK( ITAU ),
02421      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
02422                      IE = ITAU
02423                      ITAUQ = IE + M
02424                      ITAUP = ITAUQ + M
02425                      IWORK = ITAUP + M
02426 *
02427 *                    Zero out above L in A
02428 *
02429                      CALL SLASET( 'U', M-1, M-1, ZERO, ZERO, A( 1, 2 ),
02430      $                            LDA )
02431 *
02432 *                    Bidiagonalize L in A
02433 *                    (Workspace: need 4*M, prefer 3*M+2*M*NB)
02434 *
02435                      CALL SGEBRD( M, M, A, LDA, S, WORK( IE ),
02436      $                            WORK( ITAUQ ), WORK( ITAUP ),
02437      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
02438 *
02439 *                    Multiply right vectors bidiagonalizing L by Q in VT
02440 *                    (Workspace: need 3*M+N, prefer 3*M+N*NB)
02441 *
02442                      CALL SORMBR( 'P', 'L', 'T', M, N, M, A, LDA,
02443      $                            WORK( ITAUP ), VT, LDVT,
02444      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
02445                      IWORK = IE + M
02446 *
02447 *                    Perform bidiagonal QR iteration, computing right
02448 *                    singular vectors of A in VT
02449 *                    (Workspace: need BDSPAC)
02450 *
02451                      CALL SBDSQR( 'U', M, N, 0, 0, S, WORK( IE ), VT,
02452      $                            LDVT, DUM, 1, DUM, 1, WORK( IWORK ),
02453      $                            INFO )
02454 *
02455                   END IF
02456 *
02457                ELSE IF( WNTUO ) THEN
02458 *
02459 *                 Path 5t(N much larger than M, JOBU='O', JOBVT='S')
02460 *                 M right singular vectors to be computed in VT and
02461 *                 M left singular vectors to be overwritten on A
02462 *
02463                   IF( LWORK.GE.2*M*M+MAX( 4*M, BDSPAC ) ) THEN
02464 *
02465 *                    Sufficient workspace for a fast algorithm
02466 *
02467                      IU = 1
02468                      IF( LWORK.GE.WRKBL+2*LDA*M ) THEN
02469 *
02470 *                       WORK(IU) is LDA by M and WORK(IR) is LDA by M
02471 *
02472                         LDWRKU = LDA
02473                         IR = IU + LDWRKU*M
02474                         LDWRKR = LDA
02475                      ELSE IF( LWORK.GE.WRKBL+( LDA+M )*M ) THEN
02476 *
02477 *                       WORK(IU) is LDA by M and WORK(IR) is M by M
02478 *
02479                         LDWRKU = LDA
02480                         IR = IU + LDWRKU*M
02481                         LDWRKR = M
02482                      ELSE
02483 *
02484 *                       WORK(IU) is M by M and WORK(IR) is M by M
02485 *
02486                         LDWRKU = M
02487                         IR = IU + LDWRKU*M
02488                         LDWRKR = M
02489                      END IF
02490                      ITAU = IR + LDWRKR*M
02491                      IWORK = ITAU + M
02492 *
02493 *                    Compute A=L*Q
02494 *                    (Workspace: need 2*M*M+2*M, prefer 2*M*M+M+M*NB)
02495 *
02496                      CALL SGELQF( M, N, A, LDA, WORK( ITAU ),
02497      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
02498 *
02499 *                    Copy L to WORK(IU), zeroing out below it
02500 *
02501                      CALL SLACPY( 'L', M, M, A, LDA, WORK( IU ),
02502      $                            LDWRKU )
02503                      CALL SLASET( 'U', M-1, M-1, ZERO, ZERO,
02504      $                            WORK( IU+LDWRKU ), LDWRKU )
02505 *
02506 *                    Generate Q in A
02507 *                    (Workspace: need 2*M*M+2*M, prefer 2*M*M+M+M*NB)
02508 *
02509                      CALL SORGLQ( M, N, M, A, LDA, WORK( ITAU ),
02510      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
02511                      IE = ITAU
02512                      ITAUQ = IE + M
02513                      ITAUP = ITAUQ + M
02514                      IWORK = ITAUP + M
02515 *
02516 *                    Bidiagonalize L in WORK(IU), copying result to
02517 *                    WORK(IR)
02518 *                    (Workspace: need 2*M*M+4*M,
02519 *                                prefer 2*M*M+3*M+2*M*NB)
02520 *
02521                      CALL SGEBRD( M, M, WORK( IU ), LDWRKU, S,
02522      $                            WORK( IE ), WORK( ITAUQ ),
02523      $                            WORK( ITAUP ), WORK( IWORK ),
02524      $                            LWORK-IWORK+1, IERR )
02525                      CALL SLACPY( 'L', M, M, WORK( IU ), LDWRKU,
02526      $                            WORK( IR ), LDWRKR )
02527 *
02528 *                    Generate right bidiagonalizing vectors in WORK(IU)
02529 *                    (Workspace: need 2*M*M+4*M-1,
02530 *                                prefer 2*M*M+3*M+(M-1)*NB)
02531 *
02532                      CALL SORGBR( 'P', M, M, M, WORK( IU ), LDWRKU,
02533      $                            WORK( ITAUP ), WORK( IWORK ),
02534      $                            LWORK-IWORK+1, IERR )
02535 *
02536 *                    Generate left bidiagonalizing vectors in WORK(IR)
02537 *                    (Workspace: need 2*M*M+4*M, prefer 2*M*M+3*M+M*NB)
02538 *
02539                      CALL SORGBR( 'Q', M, M, M, WORK( IR ), LDWRKR,
02540      $                            WORK( ITAUQ ), WORK( IWORK ),
02541      $                            LWORK-IWORK+1, IERR )
02542                      IWORK = IE + M
02543 *
02544 *                    Perform bidiagonal QR iteration, computing left
02545 *                    singular vectors of L in WORK(IR) and computing
02546 *                    right singular vectors of L in WORK(IU)
02547 *                    (Workspace: need 2*M*M+BDSPAC)
02548 *
02549                      CALL SBDSQR( 'U', M, M, M, 0, S, WORK( IE ),
02550      $                            WORK( IU ), LDWRKU, WORK( IR ),
02551      $                            LDWRKR, DUM, 1, WORK( IWORK ), INFO )
02552 *
02553 *                    Multiply right singular vectors of L in WORK(IU) by
02554 *                    Q in A, storing result in VT
02555 *                    (Workspace: need M*M)
02556 *
02557                      CALL SGEMM( 'N', 'N', M, N, M, ONE, WORK( IU ),
02558      $                           LDWRKU, A, LDA, ZERO, VT, LDVT )
02559 *
02560 *                    Copy left singular vectors of L to A
02561 *                    (Workspace: need M*M)
02562 *
02563                      CALL SLACPY( 'F', M, M, WORK( IR ), LDWRKR, A,
02564      $                            LDA )
02565 *
02566                   ELSE
02567 *
02568 *                    Insufficient workspace for a fast algorithm
02569 *
02570                      ITAU = 1
02571                      IWORK = ITAU + M
02572 *
02573 *                    Compute A=L*Q, copying result to VT
02574 *                    (Workspace: need 2*M, prefer M+M*NB)
02575 *
02576                      CALL SGELQF( M, N, A, LDA, WORK( ITAU ),
02577      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
02578                      CALL SLACPY( 'U', M, N, A, LDA, VT, LDVT )
02579 *
02580 *                    Generate Q in VT
02581 *                    (Workspace: need 2*M, prefer M+M*NB)
02582 *
02583                      CALL SORGLQ( M, N, M, VT, LDVT, WORK( ITAU ),
02584      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
02585                      IE = ITAU
02586                      ITAUQ = IE + M
02587                      ITAUP = ITAUQ + M
02588                      IWORK = ITAUP + M
02589 *
02590 *                    Zero out above L in A
02591 *
02592                      CALL SLASET( 'U', M-1, M-1, ZERO, ZERO, A( 1, 2 ),
02593      $                            LDA )
02594 *
02595 *                    Bidiagonalize L in A
02596 *                    (Workspace: need 4*M, prefer 3*M+2*M*NB)
02597 *
02598                      CALL SGEBRD( M, M, A, LDA, S, WORK( IE ),
02599      $                            WORK( ITAUQ ), WORK( ITAUP ),
02600      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
02601 *
02602 *                    Multiply right vectors bidiagonalizing L by Q in VT
02603 *                    (Workspace: need 3*M+N, prefer 3*M+N*NB)
02604 *
02605                      CALL SORMBR( 'P', 'L', 'T', M, N, M, A, LDA,
02606      $                            WORK( ITAUP ), VT, LDVT,
02607      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
02608 *
02609 *                    Generate left bidiagonalizing vectors of L in A
02610 *                    (Workspace: need 4*M, prefer 3*M+M*NB)
02611 *
02612                      CALL SORGBR( 'Q', M, M, M, A, LDA, WORK( ITAUQ ),
02613      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
02614                      IWORK = IE + M
02615 *
02616 *                    Perform bidiagonal QR iteration, compute left
02617 *                    singular vectors of A in A and compute right
02618 *                    singular vectors of A in VT
02619 *                    (Workspace: need BDSPAC)
02620 *
02621                      CALL SBDSQR( 'U', M, N, M, 0, S, WORK( IE ), VT,
02622      $                            LDVT, A, LDA, DUM, 1, WORK( IWORK ),
02623      $                            INFO )
02624 *
02625                   END IF
02626 *
02627                ELSE IF( WNTUAS ) THEN
02628 *
02629 *                 Path 6t(N much larger than M, JOBU='S' or 'A',
02630 *                         JOBVT='S')
02631 *                 M right singular vectors to be computed in VT and
02632 *                 M left singular vectors to be computed in U
02633 *
02634                   IF( LWORK.GE.M*M+MAX( 4*M, BDSPAC ) ) THEN
02635 *
02636 *                    Sufficient workspace for a fast algorithm
02637 *
02638                      IU = 1
02639                      IF( LWORK.GE.WRKBL+LDA*M ) THEN
02640 *
02641 *                       WORK(IU) is LDA by N
02642 *
02643                         LDWRKU = LDA
02644                      ELSE
02645 *
02646 *                       WORK(IU) is LDA by M
02647 *
02648                         LDWRKU = M
02649                      END IF
02650                      ITAU = IU + LDWRKU*M
02651                      IWORK = ITAU + M
02652 *
02653 *                    Compute A=L*Q
02654 *                    (Workspace: need M*M+2*M, prefer M*M+M+M*NB)
02655 *
02656                      CALL SGELQF( M, N, A, LDA, WORK( ITAU ),
02657      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
02658 *
02659 *                    Copy L to WORK(IU), zeroing out above it
02660 *
02661                      CALL SLACPY( 'L', M, M, A, LDA, WORK( IU ),
02662      $                            LDWRKU )
02663                      CALL SLASET( 'U', M-1, M-1, ZERO, ZERO,
02664      $                            WORK( IU+LDWRKU ), LDWRKU )
02665 *
02666 *                    Generate Q in A
02667 *                    (Workspace: need M*M+2*M, prefer M*M+M+M*NB)
02668 *
02669                      CALL SORGLQ( M, N, M, A, LDA, WORK( ITAU ),
02670      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
02671                      IE = ITAU
02672                      ITAUQ = IE + M
02673                      ITAUP = ITAUQ + M
02674                      IWORK = ITAUP + M
02675 *
02676 *                    Bidiagonalize L in WORK(IU), copying result to U
02677 *                    (Workspace: need M*M+4*M, prefer M*M+3*M+2*M*NB)
02678 *
02679                      CALL SGEBRD( M, M, WORK( IU ), LDWRKU, S,
02680      $                            WORK( IE ), WORK( ITAUQ ),
02681      $                            WORK( ITAUP ), WORK( IWORK ),
02682      $                            LWORK-IWORK+1, IERR )
02683                      CALL SLACPY( 'L', M, M, WORK( IU ), LDWRKU, U,
02684      $                            LDU )
02685 *
02686 *                    Generate right bidiagonalizing vectors in WORK(IU)
02687 *                    (Workspace: need M*M+4*M-1,
02688 *                                prefer M*M+3*M+(M-1)*NB)
02689 *
02690                      CALL SORGBR( 'P', M, M, M, WORK( IU ), LDWRKU,
02691      $                            WORK( ITAUP ), WORK( IWORK ),
02692      $                            LWORK-IWORK+1, IERR )
02693 *
02694 *                    Generate left bidiagonalizing vectors in U
02695 *                    (Workspace: need M*M+4*M, prefer M*M+3*M+M*NB)
02696 *
02697                      CALL SORGBR( 'Q', M, M, M, U, LDU, WORK( ITAUQ ),
02698      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
02699                      IWORK = IE + M
02700 *
02701 *                    Perform bidiagonal QR iteration, computing left
02702 *                    singular vectors of L in U and computing right
02703 *                    singular vectors of L in WORK(IU)
02704 *                    (Workspace: need M*M+BDSPAC)
02705 *
02706                      CALL SBDSQR( 'U', M, M, M, 0, S, WORK( IE ),
02707      $                            WORK( IU ), LDWRKU, U, LDU, DUM, 1,
02708      $                            WORK( IWORK ), INFO )
02709 *
02710 *                    Multiply right singular vectors of L in WORK(IU) by
02711 *                    Q in A, storing result in VT
02712 *                    (Workspace: need M*M)
02713 *
02714                      CALL SGEMM( 'N', 'N', M, N, M, ONE, WORK( IU ),
02715      $                           LDWRKU, A, LDA, ZERO, VT, LDVT )
02716 *
02717                   ELSE
02718 *
02719 *                    Insufficient workspace for a fast algorithm
02720 *
02721                      ITAU = 1
02722                      IWORK = ITAU + M
02723 *
02724 *                    Compute A=L*Q, copying result to VT
02725 *                    (Workspace: need 2*M, prefer M+M*NB)
02726 *
02727                      CALL SGELQF( M, N, A, LDA, WORK( ITAU ),
02728      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
02729                      CALL SLACPY( 'U', M, N, A, LDA, VT, LDVT )
02730 *
02731 *                    Generate Q in VT
02732 *                    (Workspace: need 2*M, prefer M+M*NB)
02733 *
02734                      CALL SORGLQ( M, N, M, VT, LDVT, WORK( ITAU ),
02735      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
02736 *
02737 *                    Copy L to U, zeroing out above it
02738 *
02739                      CALL SLACPY( 'L', M, M, A, LDA, U, LDU )
02740                      CALL SLASET( 'U', M-1, M-1, ZERO, ZERO, U( 1, 2 ),
02741      $                            LDU )
02742                      IE = ITAU
02743                      ITAUQ = IE + M
02744                      ITAUP = ITAUQ + M
02745                      IWORK = ITAUP + M
02746 *
02747 *                    Bidiagonalize L in U
02748 *                    (Workspace: need 4*M, prefer 3*M+2*M*NB)
02749 *
02750                      CALL SGEBRD( M, M, U, LDU, S, WORK( IE ),
02751      $                            WORK( ITAUQ ), WORK( ITAUP ),
02752      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
02753 *
02754 *                    Multiply right bidiagonalizing vectors in U by Q
02755 *                    in VT
02756 *                    (Workspace: need 3*M+N, prefer 3*M+N*NB)
02757 *
02758                      CALL SORMBR( 'P', 'L', 'T', M, N, M, U, LDU,
02759      $                            WORK( ITAUP ), VT, LDVT,
02760      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
02761 *
02762 *                    Generate left bidiagonalizing vectors in U
02763 *                    (Workspace: need 4*M, prefer 3*M+M*NB)
02764 *
02765                      CALL SORGBR( 'Q', M, M, M, U, LDU, WORK( ITAUQ ),
02766      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
02767                      IWORK = IE + M
02768 *
02769 *                    Perform bidiagonal QR iteration, computing left
02770 *                    singular vectors of A in U and computing right
02771 *                    singular vectors of A in VT
02772 *                    (Workspace: need BDSPAC)
02773 *
02774                      CALL SBDSQR( 'U', M, N, M, 0, S, WORK( IE ), VT,
02775      $                            LDVT, U, LDU, DUM, 1, WORK( IWORK ),
02776      $                            INFO )
02777 *
02778                   END IF
02779 *
02780                END IF
02781 *
02782             ELSE IF( WNTVA ) THEN
02783 *
02784                IF( WNTUN ) THEN
02785 *
02786 *                 Path 7t(N much larger than M, JOBU='N', JOBVT='A')
02787 *                 N right singular vectors to be computed in VT and
02788 *                 no left singular vectors to be computed
02789 *
02790                   IF( LWORK.GE.M*M+MAX( N+M, 4*M, BDSPAC ) ) THEN
02791 *
02792 *                    Sufficient workspace for a fast algorithm
02793 *
02794                      IR = 1
02795                      IF( LWORK.GE.WRKBL+LDA*M ) THEN
02796 *
02797 *                       WORK(IR) is LDA by M
02798 *
02799                         LDWRKR = LDA
02800                      ELSE
02801 *
02802 *                       WORK(IR) is M by M
02803 *
02804                         LDWRKR = M
02805                      END IF
02806                      ITAU = IR + LDWRKR*M
02807                      IWORK = ITAU + M
02808 *
02809 *                    Compute A=L*Q, copying result to VT
02810 *                    (Workspace: need M*M+2*M, prefer M*M+M+M*NB)
02811 *
02812                      CALL SGELQF( M, N, A, LDA, WORK( ITAU ),
02813      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
02814                      CALL SLACPY( 'U', M, N, A, LDA, VT, LDVT )
02815 *
02816 *                    Copy L to WORK(IR), zeroing out above it
02817 *
02818                      CALL SLACPY( 'L', M, M, A, LDA, WORK( IR ),
02819      $                            LDWRKR )
02820                      CALL SLASET( 'U', M-1, M-1, ZERO, ZERO,
02821      $                            WORK( IR+LDWRKR ), LDWRKR )
02822 *
02823 *                    Generate Q in VT
02824 *                    (Workspace: need M*M+M+N, prefer M*M+M+N*NB)
02825 *
02826                      CALL SORGLQ( N, N, M, VT, LDVT, WORK( ITAU ),
02827      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
02828                      IE = ITAU
02829                      ITAUQ = IE + M
02830                      ITAUP = ITAUQ + M
02831                      IWORK = ITAUP + M
02832 *
02833 *                    Bidiagonalize L in WORK(IR)
02834 *                    (Workspace: need M*M+4*M, prefer M*M+3*M+2*M*NB)
02835 *
02836                      CALL SGEBRD( M, M, WORK( IR ), LDWRKR, S,
02837      $                            WORK( IE ), WORK( ITAUQ ),
02838      $                            WORK( ITAUP ), WORK( IWORK ),
02839      $                            LWORK-IWORK+1, IERR )
02840 *
02841 *                    Generate right bidiagonalizing vectors in WORK(IR)
02842 *                    (Workspace: need M*M+4*M-1,
02843 *                                prefer M*M+3*M+(M-1)*NB)
02844 *
02845                      CALL SORGBR( 'P', M, M, M, WORK( IR ), LDWRKR,
02846      $                            WORK( ITAUP ), WORK( IWORK ),
02847      $                            LWORK-IWORK+1, IERR )
02848                      IWORK = IE + M
02849 *
02850 *                    Perform bidiagonal QR iteration, computing right
02851 *                    singular vectors of L in WORK(IR)
02852 *                    (Workspace: need M*M+BDSPAC)
02853 *
02854                      CALL SBDSQR( 'U', M, M, 0, 0, S, WORK( IE ),
02855      $                            WORK( IR ), LDWRKR, DUM, 1, DUM, 1,
02856      $                            WORK( IWORK ), INFO )
02857 *
02858 *                    Multiply right singular vectors of L in WORK(IR) by
02859 *                    Q in VT, storing result in A
02860 *                    (Workspace: need M*M)
02861 *
02862                      CALL SGEMM( 'N', 'N', M, N, M, ONE, WORK( IR ),
02863      $                           LDWRKR, VT, LDVT, ZERO, A, LDA )
02864 *
02865 *                    Copy right singular vectors of A from A to VT
02866 *
02867                      CALL SLACPY( 'F', M, N, A, LDA, VT, LDVT )
02868 *
02869                   ELSE
02870 *
02871 *                    Insufficient workspace for a fast algorithm
02872 *
02873                      ITAU = 1
02874                      IWORK = ITAU + M
02875 *
02876 *                    Compute A=L*Q, copying result to VT
02877 *                    (Workspace: need 2*M, prefer M+M*NB)
02878 *
02879                      CALL SGELQF( M, N, A, LDA, WORK( ITAU ),
02880      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
02881                      CALL SLACPY( 'U', M, N, A, LDA, VT, LDVT )
02882 *
02883 *                    Generate Q in VT
02884 *                    (Workspace: need M+N, prefer M+N*NB)
02885 *
02886                      CALL SORGLQ( N, N, M, VT, LDVT, WORK( ITAU ),
02887      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
02888                      IE = ITAU
02889                      ITAUQ = IE + M
02890                      ITAUP = ITAUQ + M
02891                      IWORK = ITAUP + M
02892 *
02893 *                    Zero out above L in A
02894 *
02895                      CALL SLASET( 'U', M-1, M-1, ZERO, ZERO, A( 1, 2 ),
02896      $                            LDA )
02897 *
02898 *                    Bidiagonalize L in A
02899 *                    (Workspace: need 4*M, prefer 3*M+2*M*NB)
02900 *
02901                      CALL SGEBRD( M, M, A, LDA, S, WORK( IE ),
02902      $                            WORK( ITAUQ ), WORK( ITAUP ),
02903      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
02904 *
02905 *                    Multiply right bidiagonalizing vectors in A by Q
02906 *                    in VT
02907 *                    (Workspace: need 3*M+N, prefer 3*M+N*NB)
02908 *
02909                      CALL SORMBR( 'P', 'L', 'T', M, N, M, A, LDA,
02910      $                            WORK( ITAUP ), VT, LDVT,
02911      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
02912                      IWORK = IE + M
02913 *
02914 *                    Perform bidiagonal QR iteration, computing right
02915 *                    singular vectors of A in VT
02916 *                    (Workspace: need BDSPAC)
02917 *
02918                      CALL SBDSQR( 'U', M, N, 0, 0, S, WORK( IE ), VT,
02919      $                            LDVT, DUM, 1, DUM, 1, WORK( IWORK ),
02920      $                            INFO )
02921 *
02922                   END IF
02923 *
02924                ELSE IF( WNTUO ) THEN
02925 *
02926 *                 Path 8t(N much larger than M, JOBU='O', JOBVT='A')
02927 *                 N right singular vectors to be computed in VT and
02928 *                 M left singular vectors to be overwritten on A
02929 *
02930                   IF( LWORK.GE.2*M*M+MAX( N+M, 4*M, BDSPAC ) ) THEN
02931 *
02932 *                    Sufficient workspace for a fast algorithm
02933 *
02934                      IU = 1
02935                      IF( LWORK.GE.WRKBL+2*LDA*M ) THEN
02936 *
02937 *                       WORK(IU) is LDA by M and WORK(IR) is LDA by M
02938 *
02939                         LDWRKU = LDA
02940                         IR = IU + LDWRKU*M
02941                         LDWRKR = LDA
02942                      ELSE IF( LWORK.GE.WRKBL+( LDA+M )*M ) THEN
02943 *
02944 *                       WORK(IU) is LDA by M and WORK(IR) is M by M
02945 *
02946                         LDWRKU = LDA
02947                         IR = IU + LDWRKU*M
02948                         LDWRKR = M
02949                      ELSE
02950 *
02951 *                       WORK(IU) is M by M and WORK(IR) is M by M
02952 *
02953                         LDWRKU = M
02954                         IR = IU + LDWRKU*M
02955                         LDWRKR = M
02956                      END IF
02957                      ITAU = IR + LDWRKR*M
02958                      IWORK = ITAU + M
02959 *
02960 *                    Compute A=L*Q, copying result to VT
02961 *                    (Workspace: need 2*M*M+2*M, prefer 2*M*M+M+M*NB)
02962 *
02963                      CALL SGELQF( M, N, A, LDA, WORK( ITAU ),
02964      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
02965                      CALL SLACPY( 'U', M, N, A, LDA, VT, LDVT )
02966 *
02967 *                    Generate Q in VT
02968 *                    (Workspace: need 2*M*M+M+N, prefer 2*M*M+M+N*NB)
02969 *
02970                      CALL SORGLQ( N, N, M, VT, LDVT, WORK( ITAU ),
02971      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
02972 *
02973 *                    Copy L to WORK(IU), zeroing out above it
02974 *
02975                      CALL SLACPY( 'L', M, M, A, LDA, WORK( IU ),
02976      $                            LDWRKU )
02977                      CALL SLASET( 'U', M-1, M-1, ZERO, ZERO,
02978      $                            WORK( IU+LDWRKU ), LDWRKU )
02979                      IE = ITAU
02980                      ITAUQ = IE + M
02981                      ITAUP = ITAUQ + M
02982                      IWORK = ITAUP + M
02983 *
02984 *                    Bidiagonalize L in WORK(IU), copying result to
02985 *                    WORK(IR)
02986 *                    (Workspace: need 2*M*M+4*M,
02987 *                                prefer 2*M*M+3*M+2*M*NB)
02988 *
02989                      CALL SGEBRD( M, M, WORK( IU ), LDWRKU, S,
02990      $                            WORK( IE ), WORK( ITAUQ ),
02991      $                            WORK( ITAUP ), WORK( IWORK ),
02992      $                            LWORK-IWORK+1, IERR )
02993                      CALL SLACPY( 'L', M, M, WORK( IU ), LDWRKU,
02994      $                            WORK( IR ), LDWRKR )
02995 *
02996 *                    Generate right bidiagonalizing vectors in WORK(IU)
02997 *                    (Workspace: need 2*M*M+4*M-1,
02998 *                                prefer 2*M*M+3*M+(M-1)*NB)
02999 *
03000                      CALL SORGBR( 'P', M, M, M, WORK( IU ), LDWRKU,
03001      $                            WORK( ITAUP ), WORK( IWORK ),
03002      $                            LWORK-IWORK+1, IERR )
03003 *
03004 *                    Generate left bidiagonalizing vectors in WORK(IR)
03005 *                    (Workspace: need 2*M*M+4*M, prefer 2*M*M+3*M+M*NB)
03006 *
03007                      CALL SORGBR( 'Q', M, M, M, WORK( IR ), LDWRKR,
03008      $                            WORK( ITAUQ ), WORK( IWORK ),
03009      $                            LWORK-IWORK+1, IERR )
03010                      IWORK = IE + M
03011 *
03012 *                    Perform bidiagonal QR iteration, computing left
03013 *                    singular vectors of L in WORK(IR) and computing
03014 *                    right singular vectors of L in WORK(IU)
03015 *                    (Workspace: need 2*M*M+BDSPAC)
03016 *
03017                      CALL SBDSQR( 'U', M, M, M, 0, S, WORK( IE ),
03018      $                            WORK( IU ), LDWRKU, WORK( IR ),
03019      $                            LDWRKR, DUM, 1, WORK( IWORK ), INFO )
03020 *
03021 *                    Multiply right singular vectors of L in WORK(IU) by
03022 *                    Q in VT, storing result in A
03023 *                    (Workspace: need M*M)
03024 *
03025                      CALL SGEMM( 'N', 'N', M, N, M, ONE, WORK( IU ),
03026      $                           LDWRKU, VT, LDVT, ZERO, A, LDA )
03027 *
03028 *                    Copy right singular vectors of A from A to VT
03029 *
03030                      CALL SLACPY( 'F', M, N, A, LDA, VT, LDVT )
03031 *
03032 *                    Copy left singular vectors of A from WORK(IR) to A
03033 *
03034                      CALL SLACPY( 'F', M, M, WORK( IR ), LDWRKR, A,
03035      $                            LDA )
03036 *
03037                   ELSE
03038 *
03039 *                    Insufficient workspace for a fast algorithm
03040 *
03041                      ITAU = 1
03042                      IWORK = ITAU + M
03043 *
03044 *                    Compute A=L*Q, copying result to VT
03045 *                    (Workspace: need 2*M, prefer M+M*NB)
03046 *
03047                      CALL SGELQF( M, N, A, LDA, WORK( ITAU ),
03048      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
03049                      CALL SLACPY( 'U', M, N, A, LDA, VT, LDVT )
03050 *
03051 *                    Generate Q in VT
03052 *                    (Workspace: need M+N, prefer M+N*NB)
03053 *
03054                      CALL SORGLQ( N, N, M, VT, LDVT, WORK( ITAU ),
03055      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
03056                      IE = ITAU
03057                      ITAUQ = IE + M
03058                      ITAUP = ITAUQ + M
03059                      IWORK = ITAUP + M
03060 *
03061 *                    Zero out above L in A
03062 *
03063                      CALL SLASET( 'U', M-1, M-1, ZERO, ZERO, A( 1, 2 ),
03064      $                            LDA )
03065 *
03066 *                    Bidiagonalize L in A
03067 *                    (Workspace: need 4*M, prefer 3*M+2*M*NB)
03068 *
03069                      CALL SGEBRD( M, M, A, LDA, S, WORK( IE ),
03070      $                            WORK( ITAUQ ), WORK( ITAUP ),
03071      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
03072 *
03073 *                    Multiply right bidiagonalizing vectors in A by Q
03074 *                    in VT
03075 *                    (Workspace: need 3*M+N, prefer 3*M+N*NB)
03076 *
03077                      CALL SORMBR( 'P', 'L', 'T', M, N, M, A, LDA,
03078      $                            WORK( ITAUP ), VT, LDVT,
03079      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
03080 *
03081 *                    Generate left bidiagonalizing vectors in A
03082 *                    (Workspace: need 4*M, prefer 3*M+M*NB)
03083 *
03084                      CALL SORGBR( 'Q', M, M, M, A, LDA, WORK( ITAUQ ),
03085      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
03086                      IWORK = IE + M
03087 *
03088 *                    Perform bidiagonal QR iteration, computing left
03089 *                    singular vectors of A in A and computing right
03090 *                    singular vectors of A in VT
03091 *                    (Workspace: need BDSPAC)
03092 *
03093                      CALL SBDSQR( 'U', M, N, M, 0, S, WORK( IE ), VT,
03094      $                            LDVT, A, LDA, DUM, 1, WORK( IWORK ),
03095      $                            INFO )
03096 *
03097                   END IF
03098 *
03099                ELSE IF( WNTUAS ) THEN
03100 *
03101 *                 Path 9t(N much larger than M, JOBU='S' or 'A',
03102 *                         JOBVT='A')
03103 *                 N right singular vectors to be computed in VT and
03104 *                 M left singular vectors to be computed in U
03105 *
03106                   IF( LWORK.GE.M*M+MAX( N+M, 4*M, BDSPAC ) ) THEN
03107 *
03108 *                    Sufficient workspace for a fast algorithm
03109 *
03110                      IU = 1
03111                      IF( LWORK.GE.WRKBL+LDA*M ) THEN
03112 *
03113 *                       WORK(IU) is LDA by M
03114 *
03115                         LDWRKU = LDA
03116                      ELSE
03117 *
03118 *                       WORK(IU) is M by M
03119 *
03120                         LDWRKU = M
03121                      END IF
03122                      ITAU = IU + LDWRKU*M
03123                      IWORK = ITAU + M
03124 *
03125 *                    Compute A=L*Q, copying result to VT
03126 *                    (Workspace: need M*M+2*M, prefer M*M+M+M*NB)
03127 *
03128                      CALL SGELQF( M, N, A, LDA, WORK( ITAU ),
03129      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
03130                      CALL SLACPY( 'U', M, N, A, LDA, VT, LDVT )
03131 *
03132 *                    Generate Q in VT
03133 *                    (Workspace: need M*M+M+N, prefer M*M+M+N*NB)
03134 *
03135                      CALL SORGLQ( N, N, M, VT, LDVT, WORK( ITAU ),
03136      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
03137 *
03138 *                    Copy L to WORK(IU), zeroing out above it
03139 *
03140                      CALL SLACPY( 'L', M, M, A, LDA, WORK( IU ),
03141      $                            LDWRKU )
03142                      CALL SLASET( 'U', M-1, M-1, ZERO, ZERO,
03143      $                            WORK( IU+LDWRKU ), LDWRKU )
03144                      IE = ITAU
03145                      ITAUQ = IE + M
03146                      ITAUP = ITAUQ + M
03147                      IWORK = ITAUP + M
03148 *
03149 *                    Bidiagonalize L in WORK(IU), copying result to U
03150 *                    (Workspace: need M*M+4*M, prefer M*M+3*M+2*M*NB)
03151 *
03152                      CALL SGEBRD( M, M, WORK( IU ), LDWRKU, S,
03153      $                            WORK( IE ), WORK( ITAUQ ),
03154      $                            WORK( ITAUP ), WORK( IWORK ),
03155      $                            LWORK-IWORK+1, IERR )
03156                      CALL SLACPY( 'L', M, M, WORK( IU ), LDWRKU, U,
03157      $                            LDU )
03158 *
03159 *                    Generate right bidiagonalizing vectors in WORK(IU)
03160 *                    (Workspace: need M*M+4*M, prefer M*M+3*M+(M-1)*NB)
03161 *
03162                      CALL SORGBR( 'P', M, M, M, WORK( IU ), LDWRKU,
03163      $                            WORK( ITAUP ), WORK( IWORK ),
03164      $                            LWORK-IWORK+1, IERR )
03165 *
03166 *                    Generate left bidiagonalizing vectors in U
03167 *                    (Workspace: need M*M+4*M, prefer M*M+3*M+M*NB)
03168 *
03169                      CALL SORGBR( 'Q', M, M, M, U, LDU, WORK( ITAUQ ),
03170      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
03171                      IWORK = IE + M
03172 *
03173 *                    Perform bidiagonal QR iteration, computing left
03174 *                    singular vectors of L in U and computing right
03175 *                    singular vectors of L in WORK(IU)
03176 *                    (Workspace: need M*M+BDSPAC)
03177 *
03178                      CALL SBDSQR( 'U', M, M, M, 0, S, WORK( IE ),
03179      $                            WORK( IU ), LDWRKU, U, LDU, DUM, 1,
03180      $                            WORK( IWORK ), INFO )
03181 *
03182 *                    Multiply right singular vectors of L in WORK(IU) by
03183 *                    Q in VT, storing result in A
03184 *                    (Workspace: need M*M)
03185 *
03186                      CALL SGEMM( 'N', 'N', M, N, M, ONE, WORK( IU ),
03187      $                           LDWRKU, VT, LDVT, ZERO, A, LDA )
03188 *
03189 *                    Copy right singular vectors of A from A to VT
03190 *
03191                      CALL SLACPY( 'F', M, N, A, LDA, VT, LDVT )
03192 *
03193                   ELSE
03194 *
03195 *                    Insufficient workspace for a fast algorithm
03196 *
03197                      ITAU = 1
03198                      IWORK = ITAU + M
03199 *
03200 *                    Compute A=L*Q, copying result to VT
03201 *                    (Workspace: need 2*M, prefer M+M*NB)
03202 *
03203                      CALL SGELQF( M, N, A, LDA, WORK( ITAU ),
03204      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
03205                      CALL SLACPY( 'U', M, N, A, LDA, VT, LDVT )
03206 *
03207 *                    Generate Q in VT
03208 *                    (Workspace: need M+N, prefer M+N*NB)
03209 *
03210                      CALL SORGLQ( N, N, M, VT, LDVT, WORK( ITAU ),
03211      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
03212 *
03213 *                    Copy L to U, zeroing out above it
03214 *
03215                      CALL SLACPY( 'L', M, M, A, LDA, U, LDU )
03216                      CALL SLASET( 'U', M-1, M-1, ZERO, ZERO, U( 1, 2 ),
03217      $                            LDU )
03218                      IE = ITAU
03219                      ITAUQ = IE + M
03220                      ITAUP = ITAUQ + M
03221                      IWORK = ITAUP + M
03222 *
03223 *                    Bidiagonalize L in U
03224 *                    (Workspace: need 4*M, prefer 3*M+2*M*NB)
03225 *
03226                      CALL SGEBRD( M, M, U, LDU, S, WORK( IE ),
03227      $                            WORK( ITAUQ ), WORK( ITAUP ),
03228      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
03229 *
03230 *                    Multiply right bidiagonalizing vectors in U by Q
03231 *                    in VT
03232 *                    (Workspace: need 3*M+N, prefer 3*M+N*NB)
03233 *
03234                      CALL SORMBR( 'P', 'L', 'T', M, N, M, U, LDU,
03235      $                            WORK( ITAUP ), VT, LDVT,
03236      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
03237 *
03238 *                    Generate left bidiagonalizing vectors in U
03239 *                    (Workspace: need 4*M, prefer 3*M+M*NB)
03240 *
03241                      CALL SORGBR( 'Q', M, M, M, U, LDU, WORK( ITAUQ ),
03242      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
03243                      IWORK = IE + M
03244 *
03245 *                    Perform bidiagonal QR iteration, computing left
03246 *                    singular vectors of A in U and computing right
03247 *                    singular vectors of A in VT
03248 *                    (Workspace: need BDSPAC)
03249 *
03250                      CALL SBDSQR( 'U', M, N, M, 0, S, WORK( IE ), VT,
03251      $                            LDVT, U, LDU, DUM, 1, WORK( IWORK ),
03252      $                            INFO )
03253 *
03254                   END IF
03255 *
03256                END IF
03257 *
03258             END IF
03259 *
03260          ELSE
03261 *
03262 *           N .LT. MNTHR
03263 *
03264 *           Path 10t(N greater than M, but not much larger)
03265 *           Reduce to bidiagonal form without LQ decomposition
03266 *
03267             IE = 1
03268             ITAUQ = IE + M
03269             ITAUP = ITAUQ + M
03270             IWORK = ITAUP + M
03271 *
03272 *           Bidiagonalize A
03273 *           (Workspace: need 3*M+N, prefer 3*M+(M+N)*NB)
03274 *
03275             CALL SGEBRD( M, N, A, LDA, S, WORK( IE ), WORK( ITAUQ ),
03276      $                   WORK( ITAUP ), WORK( IWORK ), LWORK-IWORK+1,
03277      $                   IERR )
03278             IF( WNTUAS ) THEN
03279 *
03280 *              If left singular vectors desired in U, copy result to U
03281 *              and generate left bidiagonalizing vectors in U
03282 *              (Workspace: need 4*M-1, prefer 3*M+(M-1)*NB)
03283 *
03284                CALL SLACPY( 'L', M, M, A, LDA, U, LDU )
03285                CALL SORGBR( 'Q', M, M, N, U, LDU, WORK( ITAUQ ),
03286      $                      WORK( IWORK ), LWORK-IWORK+1, IERR )
03287             END IF
03288             IF( WNTVAS ) THEN
03289 *
03290 *              If right singular vectors desired in VT, copy result to
03291 *              VT and generate right bidiagonalizing vectors in VT
03292 *              (Workspace: need 3*M+NRVT, prefer 3*M+NRVT*NB)
03293 *
03294                CALL SLACPY( 'U', M, N, A, LDA, VT, LDVT )
03295                IF( WNTVA )
03296      $            NRVT = N
03297                IF( WNTVS )
03298      $            NRVT = M
03299                CALL SORGBR( 'P', NRVT, N, M, VT, LDVT, WORK( ITAUP ),
03300      $                      WORK( IWORK ), LWORK-IWORK+1, IERR )
03301             END IF
03302             IF( WNTUO ) THEN
03303 *
03304 *              If left singular vectors desired in A, generate left
03305 *              bidiagonalizing vectors in A
03306 *              (Workspace: need 4*M-1, prefer 3*M+(M-1)*NB)
03307 *
03308                CALL SORGBR( 'Q', M, M, N, A, LDA, WORK( ITAUQ ),
03309      $                      WORK( IWORK ), LWORK-IWORK+1, IERR )
03310             END IF
03311             IF( WNTVO ) THEN
03312 *
03313 *              If right singular vectors desired in A, generate right
03314 *              bidiagonalizing vectors in A
03315 *              (Workspace: need 4*M, prefer 3*M+M*NB)
03316 *
03317                CALL SORGBR( 'P', M, N, M, A, LDA, WORK( ITAUP ),
03318      $                      WORK( IWORK ), LWORK-IWORK+1, IERR )
03319             END IF
03320             IWORK = IE + M
03321             IF( WNTUAS .OR. WNTUO )
03322      $         NRU = M
03323             IF( WNTUN )
03324      $         NRU = 0
03325             IF( WNTVAS .OR. WNTVO )
03326      $         NCVT = N
03327             IF( WNTVN )
03328      $         NCVT = 0
03329             IF( ( .NOT.WNTUO ) .AND. ( .NOT.WNTVO ) ) THEN
03330 *
03331 *              Perform bidiagonal QR iteration, if desired, computing
03332 *              left singular vectors in U and computing right singular
03333 *              vectors in VT
03334 *              (Workspace: need BDSPAC)
03335 *
03336                CALL SBDSQR( 'L', M, NCVT, NRU, 0, S, WORK( IE ), VT,
03337      $                      LDVT, U, LDU, DUM, 1, WORK( IWORK ), INFO )
03338             ELSE IF( ( .NOT.WNTUO ) .AND. WNTVO ) THEN
03339 *
03340 *              Perform bidiagonal QR iteration, if desired, computing
03341 *              left singular vectors in U and computing right singular
03342 *              vectors in A
03343 *              (Workspace: need BDSPAC)
03344 *
03345                CALL SBDSQR( 'L', M, NCVT, NRU, 0, S, WORK( IE ), A, LDA,
03346      $                      U, LDU, DUM, 1, WORK( IWORK ), INFO )
03347             ELSE
03348 *
03349 *              Perform bidiagonal QR iteration, if desired, computing
03350 *              left singular vectors in A and computing right singular
03351 *              vectors in VT
03352 *              (Workspace: need BDSPAC)
03353 *
03354                CALL SBDSQR( 'L', M, NCVT, NRU, 0, S, WORK( IE ), VT,
03355      $                      LDVT, A, LDA, DUM, 1, WORK( IWORK ), INFO )
03356             END IF
03357 *
03358          END IF
03359 *
03360       END IF
03361 *
03362 *     If SBDSQR failed to converge, copy unconverged superdiagonals
03363 *     to WORK( 2:MINMN )
03364 *
03365       IF( INFO.NE.0 ) THEN
03366          IF( IE.GT.2 ) THEN
03367             DO 50 I = 1, MINMN - 1
03368                WORK( I+1 ) = WORK( I+IE-1 )
03369    50       CONTINUE
03370          END IF
03371          IF( IE.LT.2 ) THEN
03372             DO 60 I = MINMN - 1, 1, -1
03373                WORK( I+1 ) = WORK( I+IE-1 )
03374    60       CONTINUE
03375          END IF
03376       END IF
03377 *
03378 *     Undo scaling if necessary
03379 *
03380       IF( ISCL.EQ.1 ) THEN
03381          IF( ANRM.GT.BIGNUM )
03382      $      CALL SLASCL( 'G', 0, 0, BIGNUM, ANRM, MINMN, 1, S, MINMN,
03383      $                   IERR )
03384          IF( INFO.NE.0 .AND. ANRM.GT.BIGNUM )
03385      $      CALL SLASCL( 'G', 0, 0, BIGNUM, ANRM, MINMN-1, 1, WORK( 2 ),
03386      $                   MINMN, IERR )
03387          IF( ANRM.LT.SMLNUM )
03388      $      CALL SLASCL( 'G', 0, 0, SMLNUM, ANRM, MINMN, 1, S, MINMN,
03389      $                   IERR )
03390          IF( INFO.NE.0 .AND. ANRM.LT.SMLNUM )
03391      $      CALL SLASCL( 'G', 0, 0, SMLNUM, ANRM, MINMN-1, 1, WORK( 2 ),
03392      $                   MINMN, IERR )
03393       END IF
03394 *
03395 *     Return optimal workspace in WORK(1)
03396 *
03397       WORK( 1 ) = MAXWRK
03398 *
03399       RETURN
03400 *
03401 *     End of SGESVD
03402 *
03403       END