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