LAPACK: SRC/dgesvd.f Source File

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