LAPACK  3.4.1
LAPACK: Linear Algebra PACKage
dgesdd.f
Go to the documentation of this file.
00001 *> \brief \b DGESDD
00002 *
00003 *  =========== DOCUMENTATION ===========
00004 *
00005 * Online html documentation available at 
00006 *            http://www.netlib.org/lapack/explore-html/ 
00007 *
00008 *> \htmlonly
00009 *> Download DGESDD + dependencies 
00010 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dgesdd.f"> 
00011 *> [TGZ]</a> 
00012 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dgesdd.f"> 
00013 *> [ZIP]</a> 
00014 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dgesdd.f"> 
00015 *> [TXT]</a>
00016 *> \endhtmlonly 
00017 *
00018 *  Definition:
00019 *  ===========
00020 *
00021 *       SUBROUTINE DGESDD( JOBZ, M, N, A, LDA, S, U, LDU, VT, LDVT, WORK,
00022 *                          LWORK, IWORK, INFO )
00023 * 
00024 *       .. Scalar Arguments ..
00025 *       CHARACTER          JOBZ
00026 *       INTEGER            INFO, LDA, LDU, LDVT, LWORK, M, N
00027 *       ..
00028 *       .. Array Arguments ..
00029 *       INTEGER            IWORK( * )
00030 *       DOUBLE PRECISION   A( LDA, * ), S( * ), U( LDU, * ),
00031 *      $                   VT( LDVT, * ), WORK( * )
00032 *       ..
00033 *  
00034 *
00035 *> \par Purpose:
00036 *  =============
00037 *>
00038 *> \verbatim
00039 *>
00040 *> DGESDD computes the singular value decomposition (SVD) of a real
00041 *> M-by-N matrix A, optionally computing the left and right singular
00042 *> vectors.  If singular vectors are desired, it uses a
00043 *> divide-and-conquer algorithm.
00044 *>
00045 *> The SVD is written
00046 *>
00047 *>      A = U * SIGMA * transpose(V)
00048 *>
00049 *> where SIGMA is an M-by-N matrix which is zero except for its
00050 *> min(m,n) diagonal elements, U is an M-by-M orthogonal matrix, and
00051 *> V is an N-by-N orthogonal matrix.  The diagonal elements of SIGMA
00052 *> are the singular values of A; they are real and non-negative, and
00053 *> are returned in descending order.  The first min(m,n) columns of
00054 *> U and V are the left and right singular vectors of A.
00055 *>
00056 *> Note that the routine returns VT = V**T, not V.
00057 *>
00058 *> The divide and conquer algorithm makes very mild assumptions about
00059 *> floating point arithmetic. It will work on machines with a guard
00060 *> digit in add/subtract, or on those binary machines without guard
00061 *> digits which subtract like the Cray X-MP, Cray Y-MP, Cray C-90, or
00062 *> Cray-2. It could conceivably fail on hexadecimal or decimal machines
00063 *> without guard digits, but we know of none.
00064 *> \endverbatim
00065 *
00066 *  Arguments:
00067 *  ==========
00068 *
00069 *> \param[in] JOBZ
00070 *> \verbatim
00071 *>          JOBZ is CHARACTER*1
00072 *>          Specifies options for computing all or part of the matrix U:
00073 *>          = 'A':  all M columns of U and all N rows of V**T are
00074 *>                  returned in the arrays U and VT;
00075 *>          = 'S':  the first min(M,N) columns of U and the first
00076 *>                  min(M,N) rows of V**T are returned in the arrays U
00077 *>                  and VT;
00078 *>          = 'O':  If M >= N, the first N columns of U are overwritten
00079 *>                  on the array A and all rows of V**T are returned in
00080 *>                  the array VT;
00081 *>                  otherwise, all columns of U are returned in the
00082 *>                  array U and the first M rows of V**T are overwritten
00083 *>                  in the array A;
00084 *>          = 'N':  no columns of U or rows of V**T are computed.
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 JOBZ = 'O',  A is overwritten with the first N columns
00105 *>                          of U (the left singular vectors, stored
00106 *>                          columnwise) if M >= N;
00107 *>                          A is overwritten with the first M rows
00108 *>                          of V**T (the right singular vectors, stored
00109 *>                          rowwise) otherwise.
00110 *>          if JOBZ .ne. 'O', the contents of A are destroyed.
00111 *> \endverbatim
00112 *>
00113 *> \param[in] LDA
00114 *> \verbatim
00115 *>          LDA is INTEGER
00116 *>          The leading dimension of the array A.  LDA >= max(1,M).
00117 *> \endverbatim
00118 *>
00119 *> \param[out] S
00120 *> \verbatim
00121 *>          S is DOUBLE PRECISION array, dimension (min(M,N))
00122 *>          The singular values of A, sorted so that S(i) >= S(i+1).
00123 *> \endverbatim
00124 *>
00125 *> \param[out] U
00126 *> \verbatim
00127 *>          U is DOUBLE PRECISION array, dimension (LDU,UCOL)
00128 *>          UCOL = M if JOBZ = 'A' or JOBZ = 'O' and M < N;
00129 *>          UCOL = min(M,N) if JOBZ = 'S'.
00130 *>          If JOBZ = 'A' or JOBZ = 'O' and M < N, U contains the M-by-M
00131 *>          orthogonal matrix U;
00132 *>          if JOBZ = 'S', U contains the first min(M,N) columns of U
00133 *>          (the left singular vectors, stored columnwise);
00134 *>          if JOBZ = 'O' and M >= N, or JOBZ = 'N', U is not referenced.
00135 *> \endverbatim
00136 *>
00137 *> \param[in] LDU
00138 *> \verbatim
00139 *>          LDU is INTEGER
00140 *>          The leading dimension of the array U.  LDU >= 1; if
00141 *>          JOBZ = 'S' or 'A' or JOBZ = 'O' and M < N, LDU >= M.
00142 *> \endverbatim
00143 *>
00144 *> \param[out] VT
00145 *> \verbatim
00146 *>          VT is DOUBLE PRECISION array, dimension (LDVT,N)
00147 *>          If JOBZ = 'A' or JOBZ = 'O' and M >= N, VT contains the
00148 *>          N-by-N orthogonal matrix V**T;
00149 *>          if JOBZ = 'S', VT contains the first min(M,N) rows of
00150 *>          V**T (the right singular vectors, stored rowwise);
00151 *>          if JOBZ = 'O' and M < N, or JOBZ = 'N', VT is not referenced.
00152 *> \endverbatim
00153 *>
00154 *> \param[in] LDVT
00155 *> \verbatim
00156 *>          LDVT is INTEGER
00157 *>          The leading dimension of the array VT.  LDVT >= 1; if
00158 *>          JOBZ = 'A' or JOBZ = 'O' and M >= N, LDVT >= N;
00159 *>          if JOBZ = 'S', LDVT >= min(M,N).
00160 *> \endverbatim
00161 *>
00162 *> \param[out] WORK
00163 *> \verbatim
00164 *>          WORK is DOUBLE PRECISION array, dimension (MAX(1,LWORK))
00165 *>          On exit, if INFO = 0, WORK(1) returns the optimal LWORK;
00166 *> \endverbatim
00167 *>
00168 *> \param[in] LWORK
00169 *> \verbatim
00170 *>          LWORK is INTEGER
00171 *>          The dimension of the array WORK. LWORK >= 1.
00172 *>          If JOBZ = 'N',
00173 *>            LWORK >= 3*min(M,N) + max(max(M,N),7*min(M,N)).
00174 *>          If JOBZ = 'O',
00175 *>            LWORK >= 3*min(M,N) + 
00176 *>                     max(max(M,N),5*min(M,N)*min(M,N)+4*min(M,N)).
00177 *>          If JOBZ = 'S' or 'A'
00178 *>            LWORK >= 3*min(M,N) +
00179 *>                     max(max(M,N),4*min(M,N)*min(M,N)+4*min(M,N)).
00180 *>          For good performance, LWORK should generally be larger.
00181 *>          If LWORK = -1 but other input arguments are legal, WORK(1)
00182 *>          returns the optimal LWORK.
00183 *> \endverbatim
00184 *>
00185 *> \param[out] IWORK
00186 *> \verbatim
00187 *>          IWORK is INTEGER array, dimension (8*min(M,N))
00188 *> \endverbatim
00189 *>
00190 *> \param[out] INFO
00191 *> \verbatim
00192 *>          INFO is INTEGER
00193 *>          = 0:  successful exit.
00194 *>          < 0:  if INFO = -i, the i-th argument had an illegal value.
00195 *>          > 0:  DBDSDC did not converge, updating process failed.
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 November 2011
00207 *
00208 *> \ingroup doubleGEsing
00209 *
00210 *> \par Contributors:
00211 *  ==================
00212 *>
00213 *>     Ming Gu and Huan Ren, Computer Science Division, University of
00214 *>     California at Berkeley, USA
00215 *>
00216 *  =====================================================================
00217       SUBROUTINE DGESDD( JOBZ, M, N, A, LDA, S, U, LDU, VT, LDVT, WORK,
00218      $                   LWORK, IWORK, INFO )
00219 *
00220 *  -- LAPACK driver routine (version 3.4.0) --
00221 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
00222 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
00223 *     November 2011
00224 *
00225 *     .. Scalar Arguments ..
00226       CHARACTER          JOBZ
00227       INTEGER            INFO, LDA, LDU, LDVT, LWORK, M, N
00228 *     ..
00229 *     .. Array Arguments ..
00230       INTEGER            IWORK( * )
00231       DOUBLE PRECISION   A( LDA, * ), S( * ), U( LDU, * ),
00232      $                   VT( LDVT, * ), WORK( * )
00233 *     ..
00234 *
00235 *  =====================================================================
00236 *
00237 *     .. Parameters ..
00238       DOUBLE PRECISION   ZERO, ONE
00239       PARAMETER          ( ZERO = 0.0D0, ONE = 1.0D0 )
00240 *     ..
00241 *     .. Local Scalars ..
00242       LOGICAL            LQUERY, WNTQA, WNTQAS, WNTQN, WNTQO, WNTQS
00243       INTEGER            BDSPAC, BLK, CHUNK, I, IE, IERR, IL,
00244      $                   IR, ISCL, ITAU, ITAUP, ITAUQ, IU, IVT, LDWKVT,
00245      $                   LDWRKL, LDWRKR, LDWRKU, MAXWRK, MINMN, MINWRK,
00246      $                   MNTHR, NWORK, WRKBL
00247       DOUBLE PRECISION   ANRM, BIGNUM, EPS, SMLNUM
00248 *     ..
00249 *     .. Local Arrays ..
00250       INTEGER            IDUM( 1 )
00251       DOUBLE PRECISION   DUM( 1 )
00252 *     ..
00253 *     .. External Subroutines ..
00254       EXTERNAL           DBDSDC, DGEBRD, DGELQF, DGEMM, DGEQRF, DLACPY,
00255      $                   DLASCL, DLASET, DORGBR, DORGLQ, DORGQR, DORMBR,
00256      $                   XERBLA
00257 *     ..
00258 *     .. External Functions ..
00259       LOGICAL            LSAME
00260       INTEGER            ILAENV
00261       DOUBLE PRECISION   DLAMCH, DLANGE
00262       EXTERNAL           DLAMCH, DLANGE, ILAENV, LSAME
00263 *     ..
00264 *     .. Intrinsic Functions ..
00265       INTRINSIC          INT, MAX, MIN, SQRT
00266 *     ..
00267 *     .. Executable Statements ..
00268 *
00269 *     Test the input arguments
00270 *
00271       INFO = 0
00272       MINMN = MIN( M, N )
00273       WNTQA = LSAME( JOBZ, 'A' )
00274       WNTQS = LSAME( JOBZ, 'S' )
00275       WNTQAS = WNTQA .OR. WNTQS
00276       WNTQO = LSAME( JOBZ, 'O' )
00277       WNTQN = LSAME( JOBZ, 'N' )
00278       LQUERY = ( LWORK.EQ.-1 )
00279 *
00280       IF( .NOT.( WNTQA .OR. WNTQS .OR. WNTQO .OR. WNTQN ) ) THEN
00281          INFO = -1
00282       ELSE IF( M.LT.0 ) THEN
00283          INFO = -2
00284       ELSE IF( N.LT.0 ) THEN
00285          INFO = -3
00286       ELSE IF( LDA.LT.MAX( 1, M ) ) THEN
00287          INFO = -5
00288       ELSE IF( LDU.LT.1 .OR. ( WNTQAS .AND. LDU.LT.M ) .OR.
00289      $         ( WNTQO .AND. M.LT.N .AND. LDU.LT.M ) ) THEN
00290          INFO = -8
00291       ELSE IF( LDVT.LT.1 .OR. ( WNTQA .AND. LDVT.LT.N ) .OR.
00292      $         ( WNTQS .AND. LDVT.LT.MINMN ) .OR.
00293      $         ( WNTQO .AND. M.GE.N .AND. LDVT.LT.N ) ) THEN
00294          INFO = -10
00295       END IF
00296 *
00297 *     Compute workspace
00298 *      (Note: Comments in the code beginning "Workspace:" describe the
00299 *       minimal amount of workspace needed at that point in the code,
00300 *       as well as the preferred amount for good performance.
00301 *       NB refers to the optimal block size for the immediately
00302 *       following subroutine, as returned by ILAENV.)
00303 *
00304       IF( INFO.EQ.0 ) THEN
00305          MINWRK = 1
00306          MAXWRK = 1
00307          IF( M.GE.N .AND. MINMN.GT.0 ) THEN
00308 *
00309 *           Compute space needed for DBDSDC
00310 *
00311             MNTHR = INT( MINMN*11.0D0 / 6.0D0 )
00312             IF( WNTQN ) THEN
00313                BDSPAC = 7*N
00314             ELSE
00315                BDSPAC = 3*N*N + 4*N
00316             END IF
00317             IF( M.GE.MNTHR ) THEN
00318                IF( WNTQN ) THEN
00319 *
00320 *                 Path 1 (M much larger than N, JOBZ='N')
00321 *
00322                   WRKBL = N + N*ILAENV( 1, 'DGEQRF', ' ', M, N, -1,
00323      $                    -1 )
00324                   WRKBL = MAX( WRKBL, 3*N+2*N*
00325      $                    ILAENV( 1, 'DGEBRD', ' ', N, N, -1, -1 ) )
00326                   MAXWRK = MAX( WRKBL, BDSPAC+N )
00327                   MINWRK = BDSPAC + N
00328                ELSE IF( WNTQO ) THEN
00329 *
00330 *                 Path 2 (M much larger than N, JOBZ='O')
00331 *
00332                   WRKBL = N + N*ILAENV( 1, 'DGEQRF', ' ', M, N, -1, -1 )
00333                   WRKBL = MAX( WRKBL, N+N*ILAENV( 1, 'DORGQR', ' ', M,
00334      $                    N, N, -1 ) )
00335                   WRKBL = MAX( WRKBL, 3*N+2*N*
00336      $                    ILAENV( 1, 'DGEBRD', ' ', N, N, -1, -1 ) )
00337                   WRKBL = MAX( WRKBL, 3*N+N*
00338      $                    ILAENV( 1, 'DORMBR', 'QLN', N, N, N, -1 ) )
00339                   WRKBL = MAX( WRKBL, 3*N+N*
00340      $                    ILAENV( 1, 'DORMBR', 'PRT', N, N, N, -1 ) )
00341                   WRKBL = MAX( WRKBL, BDSPAC+3*N )
00342                   MAXWRK = WRKBL + 2*N*N
00343                   MINWRK = BDSPAC + 2*N*N + 3*N
00344                ELSE IF( WNTQS ) THEN
00345 *
00346 *                 Path 3 (M much larger than N, JOBZ='S')
00347 *
00348                   WRKBL = N + N*ILAENV( 1, 'DGEQRF', ' ', M, N, -1, -1 )
00349                   WRKBL = MAX( WRKBL, N+N*ILAENV( 1, 'DORGQR', ' ', M,
00350      $                    N, N, -1 ) )
00351                   WRKBL = MAX( WRKBL, 3*N+2*N*
00352      $                    ILAENV( 1, 'DGEBRD', ' ', N, N, -1, -1 ) )
00353                   WRKBL = MAX( WRKBL, 3*N+N*
00354      $                    ILAENV( 1, 'DORMBR', 'QLN', N, N, N, -1 ) )
00355                   WRKBL = MAX( WRKBL, 3*N+N*
00356      $                    ILAENV( 1, 'DORMBR', 'PRT', N, N, N, -1 ) )
00357                   WRKBL = MAX( WRKBL, BDSPAC+3*N )
00358                   MAXWRK = WRKBL + N*N
00359                   MINWRK = BDSPAC + N*N + 3*N
00360                ELSE IF( WNTQA ) THEN
00361 *
00362 *                 Path 4 (M much larger than N, JOBZ='A')
00363 *
00364                   WRKBL = N + N*ILAENV( 1, 'DGEQRF', ' ', M, N, -1, -1 )
00365                   WRKBL = MAX( WRKBL, N+M*ILAENV( 1, 'DORGQR', ' ', M,
00366      $                    M, N, -1 ) )
00367                   WRKBL = MAX( WRKBL, 3*N+2*N*
00368      $                    ILAENV( 1, 'DGEBRD', ' ', N, N, -1, -1 ) )
00369                   WRKBL = MAX( WRKBL, 3*N+N*
00370      $                    ILAENV( 1, 'DORMBR', 'QLN', N, N, N, -1 ) )
00371                   WRKBL = MAX( WRKBL, 3*N+N*
00372      $                    ILAENV( 1, 'DORMBR', 'PRT', N, N, N, -1 ) )
00373                   WRKBL = MAX( WRKBL, BDSPAC+3*N )
00374                   MAXWRK = WRKBL + N*N
00375                   MINWRK = BDSPAC + N*N + 3*N
00376                END IF
00377             ELSE
00378 *
00379 *              Path 5 (M at least N, but not much larger)
00380 *
00381                WRKBL = 3*N + ( M+N )*ILAENV( 1, 'DGEBRD', ' ', M, N, -1,
00382      $                 -1 )
00383                IF( WNTQN ) THEN
00384                   MAXWRK = MAX( WRKBL, BDSPAC+3*N )
00385                   MINWRK = 3*N + MAX( M, BDSPAC )
00386                ELSE IF( WNTQO ) THEN
00387                   WRKBL = MAX( WRKBL, 3*N+N*
00388      $                    ILAENV( 1, 'DORMBR', 'QLN', M, N, N, -1 ) )
00389                   WRKBL = MAX( WRKBL, 3*N+N*
00390      $                    ILAENV( 1, 'DORMBR', 'PRT', N, N, N, -1 ) )
00391                   WRKBL = MAX( WRKBL, BDSPAC+3*N )
00392                   MAXWRK = WRKBL + M*N
00393                   MINWRK = 3*N + MAX( M, N*N+BDSPAC )
00394                ELSE IF( WNTQS ) THEN
00395                   WRKBL = MAX( WRKBL, 3*N+N*
00396      $                    ILAENV( 1, 'DORMBR', 'QLN', M, N, N, -1 ) )
00397                   WRKBL = MAX( WRKBL, 3*N+N*
00398      $                    ILAENV( 1, 'DORMBR', 'PRT', N, N, N, -1 ) )
00399                   MAXWRK = MAX( WRKBL, BDSPAC+3*N )
00400                   MINWRK = 3*N + MAX( M, BDSPAC )
00401                ELSE IF( WNTQA ) THEN
00402                   WRKBL = MAX( WRKBL, 3*N+M*
00403      $                    ILAENV( 1, 'DORMBR', 'QLN', M, M, N, -1 ) )
00404                   WRKBL = MAX( WRKBL, 3*N+N*
00405      $                    ILAENV( 1, 'DORMBR', 'PRT', N, N, N, -1 ) )
00406                   MAXWRK = MAX( MAXWRK, BDSPAC+3*N )
00407                   MINWRK = 3*N + MAX( M, BDSPAC )
00408                END IF
00409             END IF
00410          ELSE IF( MINMN.GT.0 ) THEN
00411 *
00412 *           Compute space needed for DBDSDC
00413 *
00414             MNTHR = INT( MINMN*11.0D0 / 6.0D0 )
00415             IF( WNTQN ) THEN
00416                BDSPAC = 7*M
00417             ELSE
00418                BDSPAC = 3*M*M + 4*M
00419             END IF
00420             IF( N.GE.MNTHR ) THEN
00421                IF( WNTQN ) THEN
00422 *
00423 *                 Path 1t (N much larger than M, JOBZ='N')
00424 *
00425                   WRKBL = M + M*ILAENV( 1, 'DGELQF', ' ', M, N, -1,
00426      $                    -1 )
00427                   WRKBL = MAX( WRKBL, 3*M+2*M*
00428      $                    ILAENV( 1, 'DGEBRD', ' ', M, M, -1, -1 ) )
00429                   MAXWRK = MAX( WRKBL, BDSPAC+M )
00430                   MINWRK = BDSPAC + M
00431                ELSE IF( WNTQO ) THEN
00432 *
00433 *                 Path 2t (N much larger than M, JOBZ='O')
00434 *
00435                   WRKBL = M + M*ILAENV( 1, 'DGELQF', ' ', M, N, -1, -1 )
00436                   WRKBL = MAX( WRKBL, M+M*ILAENV( 1, 'DORGLQ', ' ', M,
00437      $                    N, M, -1 ) )
00438                   WRKBL = MAX( WRKBL, 3*M+2*M*
00439      $                    ILAENV( 1, 'DGEBRD', ' ', M, M, -1, -1 ) )
00440                   WRKBL = MAX( WRKBL, 3*M+M*
00441      $                    ILAENV( 1, 'DORMBR', 'QLN', M, M, M, -1 ) )
00442                   WRKBL = MAX( WRKBL, 3*M+M*
00443      $                    ILAENV( 1, 'DORMBR', 'PRT', M, M, M, -1 ) )
00444                   WRKBL = MAX( WRKBL, BDSPAC+3*M )
00445                   MAXWRK = WRKBL + 2*M*M
00446                   MINWRK = BDSPAC + 2*M*M + 3*M
00447                ELSE IF( WNTQS ) THEN
00448 *
00449 *                 Path 3t (N much larger than M, JOBZ='S')
00450 *
00451                   WRKBL = M + M*ILAENV( 1, 'DGELQF', ' ', M, N, -1, -1 )
00452                   WRKBL = MAX( WRKBL, M+M*ILAENV( 1, 'DORGLQ', ' ', M,
00453      $                    N, M, -1 ) )
00454                   WRKBL = MAX( WRKBL, 3*M+2*M*
00455      $                    ILAENV( 1, 'DGEBRD', ' ', M, M, -1, -1 ) )
00456                   WRKBL = MAX( WRKBL, 3*M+M*
00457      $                    ILAENV( 1, 'DORMBR', 'QLN', M, M, M, -1 ) )
00458                   WRKBL = MAX( WRKBL, 3*M+M*
00459      $                    ILAENV( 1, 'DORMBR', 'PRT', M, M, M, -1 ) )
00460                   WRKBL = MAX( WRKBL, BDSPAC+3*M )
00461                   MAXWRK = WRKBL + M*M
00462                   MINWRK = BDSPAC + M*M + 3*M
00463                ELSE IF( WNTQA ) THEN
00464 *
00465 *                 Path 4t (N much larger than M, JOBZ='A')
00466 *
00467                   WRKBL = M + M*ILAENV( 1, 'DGELQF', ' ', M, N, -1, -1 )
00468                   WRKBL = MAX( WRKBL, M+N*ILAENV( 1, 'DORGLQ', ' ', N,
00469      $                    N, M, -1 ) )
00470                   WRKBL = MAX( WRKBL, 3*M+2*M*
00471      $                    ILAENV( 1, 'DGEBRD', ' ', M, M, -1, -1 ) )
00472                   WRKBL = MAX( WRKBL, 3*M+M*
00473      $                    ILAENV( 1, 'DORMBR', 'QLN', M, M, M, -1 ) )
00474                   WRKBL = MAX( WRKBL, 3*M+M*
00475      $                    ILAENV( 1, 'DORMBR', 'PRT', M, M, M, -1 ) )
00476                   WRKBL = MAX( WRKBL, BDSPAC+3*M )
00477                   MAXWRK = WRKBL + M*M
00478                   MINWRK = BDSPAC + M*M + 3*M
00479                END IF
00480             ELSE
00481 *
00482 *              Path 5t (N greater than M, but not much larger)
00483 *
00484                WRKBL = 3*M + ( M+N )*ILAENV( 1, 'DGEBRD', ' ', M, N, -1,
00485      $                 -1 )
00486                IF( WNTQN ) THEN
00487                   MAXWRK = MAX( WRKBL, BDSPAC+3*M )
00488                   MINWRK = 3*M + MAX( N, BDSPAC )
00489                ELSE IF( WNTQO ) THEN
00490                   WRKBL = MAX( WRKBL, 3*M+M*
00491      $                    ILAENV( 1, 'DORMBR', 'QLN', M, M, N, -1 ) )
00492                   WRKBL = MAX( WRKBL, 3*M+M*
00493      $                    ILAENV( 1, 'DORMBR', 'PRT', M, N, M, -1 ) )
00494                   WRKBL = MAX( WRKBL, BDSPAC+3*M )
00495                   MAXWRK = WRKBL + M*N
00496                   MINWRK = 3*M + MAX( N, M*M+BDSPAC )
00497                ELSE IF( WNTQS ) THEN
00498                   WRKBL = MAX( WRKBL, 3*M+M*
00499      $                    ILAENV( 1, 'DORMBR', 'QLN', M, M, N, -1 ) )
00500                   WRKBL = MAX( WRKBL, 3*M+M*
00501      $                    ILAENV( 1, 'DORMBR', 'PRT', M, N, M, -1 ) )
00502                   MAXWRK = MAX( WRKBL, BDSPAC+3*M )
00503                   MINWRK = 3*M + MAX( N, BDSPAC )
00504                ELSE IF( WNTQA ) THEN
00505                   WRKBL = MAX( WRKBL, 3*M+M*
00506      $                    ILAENV( 1, 'DORMBR', 'QLN', M, M, N, -1 ) )
00507                   WRKBL = MAX( WRKBL, 3*M+M*
00508      $                    ILAENV( 1, 'DORMBR', 'PRT', N, N, M, -1 ) )
00509                   MAXWRK = MAX( WRKBL, BDSPAC+3*M )
00510                   MINWRK = 3*M + MAX( N, BDSPAC )
00511                END IF
00512             END IF
00513          END IF
00514          MAXWRK = MAX( MAXWRK, MINWRK )
00515          WORK( 1 ) = MAXWRK
00516 *
00517          IF( LWORK.LT.MINWRK .AND. .NOT.LQUERY ) THEN
00518             INFO = -12
00519          END IF
00520       END IF
00521 *
00522       IF( INFO.NE.0 ) THEN
00523          CALL XERBLA( 'DGESDD', -INFO )
00524          RETURN
00525       ELSE IF( LQUERY ) THEN
00526          RETURN
00527       END IF
00528 *
00529 *     Quick return if possible
00530 *
00531       IF( M.EQ.0 .OR. N.EQ.0 ) THEN
00532          RETURN
00533       END IF
00534 *
00535 *     Get machine constants
00536 *
00537       EPS = DLAMCH( 'P' )
00538       SMLNUM = SQRT( DLAMCH( 'S' ) ) / EPS
00539       BIGNUM = ONE / SMLNUM
00540 *
00541 *     Scale A if max element outside range [SMLNUM,BIGNUM]
00542 *
00543       ANRM = DLANGE( 'M', M, N, A, LDA, DUM )
00544       ISCL = 0
00545       IF( ANRM.GT.ZERO .AND. ANRM.LT.SMLNUM ) THEN
00546          ISCL = 1
00547          CALL DLASCL( 'G', 0, 0, ANRM, SMLNUM, M, N, A, LDA, IERR )
00548       ELSE IF( ANRM.GT.BIGNUM ) THEN
00549          ISCL = 1
00550          CALL DLASCL( 'G', 0, 0, ANRM, BIGNUM, M, N, A, LDA, IERR )
00551       END IF
00552 *
00553       IF( M.GE.N ) THEN
00554 *
00555 *        A has at least as many rows as columns. If A has sufficiently
00556 *        more rows than columns, first reduce using the QR
00557 *        decomposition (if sufficient workspace available)
00558 *
00559          IF( M.GE.MNTHR ) THEN
00560 *
00561             IF( WNTQN ) THEN
00562 *
00563 *              Path 1 (M much larger than N, JOBZ='N')
00564 *              No singular vectors to be computed
00565 *
00566                ITAU = 1
00567                NWORK = ITAU + N
00568 *
00569 *              Compute A=Q*R
00570 *              (Workspace: need 2*N, prefer N+N*NB)
00571 *
00572                CALL DGEQRF( M, N, A, LDA, WORK( ITAU ), WORK( NWORK ),
00573      $                      LWORK-NWORK+1, IERR )
00574 *
00575 *              Zero out below R
00576 *
00577                CALL DLASET( 'L', N-1, N-1, ZERO, ZERO, A( 2, 1 ), LDA )
00578                IE = 1
00579                ITAUQ = IE + N
00580                ITAUP = ITAUQ + N
00581                NWORK = ITAUP + N
00582 *
00583 *              Bidiagonalize R in A
00584 *              (Workspace: need 4*N, prefer 3*N+2*N*NB)
00585 *
00586                CALL DGEBRD( N, N, A, LDA, S, WORK( IE ), WORK( ITAUQ ),
00587      $                      WORK( ITAUP ), WORK( NWORK ), LWORK-NWORK+1,
00588      $                      IERR )
00589                NWORK = IE + N
00590 *
00591 *              Perform bidiagonal SVD, computing singular values only
00592 *              (Workspace: need N+BDSPAC)
00593 *
00594                CALL DBDSDC( 'U', 'N', N, S, WORK( IE ), DUM, 1, DUM, 1,
00595      $                      DUM, IDUM, WORK( NWORK ), IWORK, INFO )
00596 *
00597             ELSE IF( WNTQO ) THEN
00598 *
00599 *              Path 2 (M much larger than N, JOBZ = 'O')
00600 *              N left singular vectors to be overwritten on A and
00601 *              N right singular vectors to be computed in VT
00602 *
00603                IR = 1
00604 *
00605 *              WORK(IR) is LDWRKR by N
00606 *
00607                IF( LWORK.GE.LDA*N+N*N+3*N+BDSPAC ) THEN
00608                   LDWRKR = LDA
00609                ELSE
00610                   LDWRKR = ( LWORK-N*N-3*N-BDSPAC ) / N
00611                END IF
00612                ITAU = IR + LDWRKR*N
00613                NWORK = ITAU + N
00614 *
00615 *              Compute A=Q*R
00616 *              (Workspace: need N*N+2*N, prefer N*N+N+N*NB)
00617 *
00618                CALL DGEQRF( M, N, A, LDA, WORK( ITAU ), WORK( NWORK ),
00619      $                      LWORK-NWORK+1, IERR )
00620 *
00621 *              Copy R to WORK(IR), zeroing out below it
00622 *
00623                CALL DLACPY( 'U', N, N, A, LDA, WORK( IR ), LDWRKR )
00624                CALL DLASET( 'L', N-1, N-1, ZERO, ZERO, WORK( IR+1 ),
00625      $                      LDWRKR )
00626 *
00627 *              Generate Q in A
00628 *              (Workspace: need N*N+2*N, prefer N*N+N+N*NB)
00629 *
00630                CALL DORGQR( M, N, N, A, LDA, WORK( ITAU ),
00631      $                      WORK( NWORK ), LWORK-NWORK+1, IERR )
00632                IE = ITAU
00633                ITAUQ = IE + N
00634                ITAUP = ITAUQ + N
00635                NWORK = ITAUP + N
00636 *
00637 *              Bidiagonalize R in VT, copying result to WORK(IR)
00638 *              (Workspace: need N*N+4*N, prefer N*N+3*N+2*N*NB)
00639 *
00640                CALL DGEBRD( N, N, WORK( IR ), LDWRKR, S, WORK( IE ),
00641      $                      WORK( ITAUQ ), WORK( ITAUP ), WORK( NWORK ),
00642      $                      LWORK-NWORK+1, IERR )
00643 *
00644 *              WORK(IU) is N by N
00645 *
00646                IU = NWORK
00647                NWORK = IU + N*N
00648 *
00649 *              Perform bidiagonal SVD, computing left singular vectors
00650 *              of bidiagonal matrix in WORK(IU) and computing right
00651 *              singular vectors of bidiagonal matrix in VT
00652 *              (Workspace: need N+N*N+BDSPAC)
00653 *
00654                CALL DBDSDC( 'U', 'I', N, S, WORK( IE ), WORK( IU ), N,
00655      $                      VT, LDVT, DUM, IDUM, WORK( NWORK ), IWORK,
00656      $                      INFO )
00657 *
00658 *              Overwrite WORK(IU) by left singular vectors of R
00659 *              and VT by right singular vectors of R
00660 *              (Workspace: need 2*N*N+3*N, prefer 2*N*N+2*N+N*NB)
00661 *
00662                CALL DORMBR( 'Q', 'L', 'N', N, N, N, WORK( IR ), LDWRKR,
00663      $                      WORK( ITAUQ ), WORK( IU ), N, WORK( NWORK ),
00664      $                      LWORK-NWORK+1, IERR )
00665                CALL DORMBR( 'P', 'R', 'T', N, N, N, WORK( IR ), LDWRKR,
00666      $                      WORK( ITAUP ), VT, LDVT, WORK( NWORK ),
00667      $                      LWORK-NWORK+1, IERR )
00668 *
00669 *              Multiply Q in A by left singular vectors of R in
00670 *              WORK(IU), storing result in WORK(IR) and copying to A
00671 *              (Workspace: need 2*N*N, prefer N*N+M*N)
00672 *
00673                DO 10 I = 1, M, LDWRKR
00674                   CHUNK = MIN( M-I+1, LDWRKR )
00675                   CALL DGEMM( 'N', 'N', CHUNK, N, N, ONE, A( I, 1 ),
00676      $                        LDA, WORK( IU ), N, ZERO, WORK( IR ),
00677      $                        LDWRKR )
00678                   CALL DLACPY( 'F', CHUNK, N, WORK( IR ), LDWRKR,
00679      $                         A( I, 1 ), LDA )
00680    10          CONTINUE
00681 *
00682             ELSE IF( WNTQS ) THEN
00683 *
00684 *              Path 3 (M much larger than N, JOBZ='S')
00685 *              N left singular vectors to be computed in U and
00686 *              N right singular vectors to be computed in VT
00687 *
00688                IR = 1
00689 *
00690 *              WORK(IR) is N by N
00691 *
00692                LDWRKR = N
00693                ITAU = IR + LDWRKR*N
00694                NWORK = ITAU + N
00695 *
00696 *              Compute A=Q*R
00697 *              (Workspace: need N*N+2*N, prefer N*N+N+N*NB)
00698 *
00699                CALL DGEQRF( M, N, A, LDA, WORK( ITAU ), WORK( NWORK ),
00700      $                      LWORK-NWORK+1, IERR )
00701 *
00702 *              Copy R to WORK(IR), zeroing out below it
00703 *
00704                CALL DLACPY( 'U', N, N, A, LDA, WORK( IR ), LDWRKR )
00705                CALL DLASET( 'L', N-1, N-1, ZERO, ZERO, WORK( IR+1 ),
00706      $                      LDWRKR )
00707 *
00708 *              Generate Q in A
00709 *              (Workspace: need N*N+2*N, prefer N*N+N+N*NB)
00710 *
00711                CALL DORGQR( M, N, N, A, LDA, WORK( ITAU ),
00712      $                      WORK( NWORK ), LWORK-NWORK+1, IERR )
00713                IE = ITAU
00714                ITAUQ = IE + N
00715                ITAUP = ITAUQ + N
00716                NWORK = ITAUP + N
00717 *
00718 *              Bidiagonalize R in WORK(IR)
00719 *              (Workspace: need N*N+4*N, prefer N*N+3*N+2*N*NB)
00720 *
00721                CALL DGEBRD( N, N, WORK( IR ), LDWRKR, S, WORK( IE ),
00722      $                      WORK( ITAUQ ), WORK( ITAUP ), WORK( NWORK ),
00723      $                      LWORK-NWORK+1, IERR )
00724 *
00725 *              Perform bidiagonal SVD, computing left singular vectors
00726 *              of bidiagoal matrix in U and computing right singular
00727 *              vectors of bidiagonal matrix in VT
00728 *              (Workspace: need N+BDSPAC)
00729 *
00730                CALL DBDSDC( 'U', 'I', N, S, WORK( IE ), U, LDU, VT,
00731      $                      LDVT, DUM, IDUM, WORK( NWORK ), IWORK,
00732      $                      INFO )
00733 *
00734 *              Overwrite U by left singular vectors of R and VT
00735 *              by right singular vectors of R
00736 *              (Workspace: need N*N+3*N, prefer N*N+2*N+N*NB)
00737 *
00738                CALL DORMBR( 'Q', 'L', 'N', N, N, N, WORK( IR ), LDWRKR,
00739      $                      WORK( ITAUQ ), U, LDU, WORK( NWORK ),
00740      $                      LWORK-NWORK+1, IERR )
00741 *
00742                CALL DORMBR( 'P', 'R', 'T', N, N, N, WORK( IR ), LDWRKR,
00743      $                      WORK( ITAUP ), VT, LDVT, WORK( NWORK ),
00744      $                      LWORK-NWORK+1, IERR )
00745 *
00746 *              Multiply Q in A by left singular vectors of R in
00747 *              WORK(IR), storing result in U
00748 *              (Workspace: need N*N)
00749 *
00750                CALL DLACPY( 'F', N, N, U, LDU, WORK( IR ), LDWRKR )
00751                CALL DGEMM( 'N', 'N', M, N, N, ONE, A, LDA, WORK( IR ),
00752      $                     LDWRKR, ZERO, U, LDU )
00753 *
00754             ELSE IF( WNTQA ) THEN
00755 *
00756 *              Path 4 (M much larger than N, JOBZ='A')
00757 *              M left singular vectors to be computed in U and
00758 *              N right singular vectors to be computed in VT
00759 *
00760                IU = 1
00761 *
00762 *              WORK(IU) is N by N
00763 *
00764                LDWRKU = N
00765                ITAU = IU + LDWRKU*N
00766                NWORK = ITAU + N
00767 *
00768 *              Compute A=Q*R, copying result to U
00769 *              (Workspace: need N*N+2*N, prefer N*N+N+N*NB)
00770 *
00771                CALL DGEQRF( M, N, A, LDA, WORK( ITAU ), WORK( NWORK ),
00772      $                      LWORK-NWORK+1, IERR )
00773                CALL DLACPY( 'L', M, N, A, LDA, U, LDU )
00774 *
00775 *              Generate Q in U
00776 *              (Workspace: need N*N+2*N, prefer N*N+N+N*NB)
00777                CALL DORGQR( M, M, N, U, LDU, WORK( ITAU ),
00778      $                      WORK( NWORK ), LWORK-NWORK+1, IERR )
00779 *
00780 *              Produce R in A, zeroing out other entries
00781 *
00782                CALL DLASET( 'L', N-1, N-1, ZERO, ZERO, A( 2, 1 ), LDA )
00783                IE = ITAU
00784                ITAUQ = IE + N
00785                ITAUP = ITAUQ + N
00786                NWORK = ITAUP + N
00787 *
00788 *              Bidiagonalize R in A
00789 *              (Workspace: need N*N+4*N, prefer N*N+3*N+2*N*NB)
00790 *
00791                CALL DGEBRD( N, N, A, LDA, S, WORK( IE ), WORK( ITAUQ ),
00792      $                      WORK( ITAUP ), WORK( NWORK ), LWORK-NWORK+1,
00793      $                      IERR )
00794 *
00795 *              Perform bidiagonal SVD, computing left singular vectors
00796 *              of bidiagonal matrix in WORK(IU) and computing right
00797 *              singular vectors of bidiagonal matrix in VT
00798 *              (Workspace: need N+N*N+BDSPAC)
00799 *
00800                CALL DBDSDC( 'U', 'I', N, S, WORK( IE ), WORK( IU ), N,
00801      $                      VT, LDVT, DUM, IDUM, WORK( NWORK ), IWORK,
00802      $                      INFO )
00803 *
00804 *              Overwrite WORK(IU) by left singular vectors of R and VT
00805 *              by right singular vectors of R
00806 *              (Workspace: need N*N+3*N, prefer N*N+2*N+N*NB)
00807 *
00808                CALL DORMBR( 'Q', 'L', 'N', N, N, N, A, LDA,
00809      $                      WORK( ITAUQ ), WORK( IU ), LDWRKU,
00810      $                      WORK( NWORK ), LWORK-NWORK+1, IERR )
00811                CALL DORMBR( 'P', 'R', 'T', N, N, N, A, LDA,
00812      $                      WORK( ITAUP ), VT, LDVT, WORK( NWORK ),
00813      $                      LWORK-NWORK+1, IERR )
00814 *
00815 *              Multiply Q in U by left singular vectors of R in
00816 *              WORK(IU), storing result in A
00817 *              (Workspace: need N*N)
00818 *
00819                CALL DGEMM( 'N', 'N', M, N, N, ONE, U, LDU, WORK( IU ),
00820      $                     LDWRKU, ZERO, A, LDA )
00821 *
00822 *              Copy left singular vectors of A from A to U
00823 *
00824                CALL DLACPY( 'F', M, N, A, LDA, U, LDU )
00825 *
00826             END IF
00827 *
00828          ELSE
00829 *
00830 *           M .LT. MNTHR
00831 *
00832 *           Path 5 (M at least N, but not much larger)
00833 *           Reduce to bidiagonal form without QR decomposition
00834 *
00835             IE = 1
00836             ITAUQ = IE + N
00837             ITAUP = ITAUQ + N
00838             NWORK = ITAUP + N
00839 *
00840 *           Bidiagonalize A
00841 *           (Workspace: need 3*N+M, prefer 3*N+(M+N)*NB)
00842 *
00843             CALL DGEBRD( M, N, A, LDA, S, WORK( IE ), WORK( ITAUQ ),
00844      $                   WORK( ITAUP ), WORK( NWORK ), LWORK-NWORK+1,
00845      $                   IERR )
00846             IF( WNTQN ) THEN
00847 *
00848 *              Perform bidiagonal SVD, only computing singular values
00849 *              (Workspace: need N+BDSPAC)
00850 *
00851                CALL DBDSDC( 'U', 'N', N, S, WORK( IE ), DUM, 1, DUM, 1,
00852      $                      DUM, IDUM, WORK( NWORK ), IWORK, INFO )
00853             ELSE IF( WNTQO ) THEN
00854                IU = NWORK
00855                IF( LWORK.GE.M*N+3*N+BDSPAC ) THEN
00856 *
00857 *                 WORK( IU ) is M by N
00858 *
00859                   LDWRKU = M
00860                   NWORK = IU + LDWRKU*N
00861                   CALL DLASET( 'F', M, N, ZERO, ZERO, WORK( IU ),
00862      $                         LDWRKU )
00863                ELSE
00864 *
00865 *                 WORK( IU ) is N by N
00866 *
00867                   LDWRKU = N
00868                   NWORK = IU + LDWRKU*N
00869 *
00870 *                 WORK(IR) is LDWRKR by N
00871 *
00872                   IR = NWORK
00873                   LDWRKR = ( LWORK-N*N-3*N ) / N
00874                END IF
00875                NWORK = IU + LDWRKU*N
00876 *
00877 *              Perform bidiagonal SVD, computing left singular vectors
00878 *              of bidiagonal matrix in WORK(IU) and computing right
00879 *              singular vectors of bidiagonal matrix in VT
00880 *              (Workspace: need N+N*N+BDSPAC)
00881 *
00882                CALL DBDSDC( 'U', 'I', N, S, WORK( IE ), WORK( IU ),
00883      $                      LDWRKU, VT, LDVT, DUM, IDUM, WORK( NWORK ),
00884      $                      IWORK, INFO )
00885 *
00886 *              Overwrite VT by right singular vectors of A
00887 *              (Workspace: need N*N+2*N, prefer N*N+N+N*NB)
00888 *
00889                CALL DORMBR( 'P', 'R', 'T', N, N, N, A, LDA,
00890      $                      WORK( ITAUP ), VT, LDVT, WORK( NWORK ),
00891      $                      LWORK-NWORK+1, IERR )
00892 *
00893                IF( LWORK.GE.M*N+3*N+BDSPAC ) THEN
00894 *
00895 *                 Overwrite WORK(IU) by left singular vectors of A
00896 *                 (Workspace: need N*N+2*N, prefer N*N+N+N*NB)
00897 *
00898                   CALL DORMBR( 'Q', 'L', 'N', M, N, N, A, LDA,
00899      $                         WORK( ITAUQ ), WORK( IU ), LDWRKU,
00900      $                         WORK( NWORK ), LWORK-NWORK+1, IERR )
00901 *
00902 *                 Copy left singular vectors of A from WORK(IU) to A
00903 *
00904                   CALL DLACPY( 'F', M, N, WORK( IU ), LDWRKU, A, LDA )
00905                ELSE
00906 *
00907 *                 Generate Q in A
00908 *                 (Workspace: need N*N+2*N, prefer N*N+N+N*NB)
00909 *
00910                   CALL DORGBR( 'Q', M, N, N, A, LDA, WORK( ITAUQ ),
00911      $                         WORK( NWORK ), LWORK-NWORK+1, IERR )
00912 *
00913 *                 Multiply Q in A by left singular vectors of
00914 *                 bidiagonal matrix in WORK(IU), storing result in
00915 *                 WORK(IR) and copying to A
00916 *                 (Workspace: need 2*N*N, prefer N*N+M*N)
00917 *
00918                   DO 20 I = 1, M, LDWRKR
00919                      CHUNK = MIN( M-I+1, LDWRKR )
00920                      CALL DGEMM( 'N', 'N', CHUNK, N, N, ONE, A( I, 1 ),
00921      $                           LDA, WORK( IU ), LDWRKU, ZERO,
00922      $                           WORK( IR ), LDWRKR )
00923                      CALL DLACPY( 'F', CHUNK, N, WORK( IR ), LDWRKR,
00924      $                            A( I, 1 ), LDA )
00925    20             CONTINUE
00926                END IF
00927 *
00928             ELSE IF( WNTQS ) THEN
00929 *
00930 *              Perform bidiagonal SVD, computing left singular vectors
00931 *              of bidiagonal matrix in U and computing right singular
00932 *              vectors of bidiagonal matrix in VT
00933 *              (Workspace: need N+BDSPAC)
00934 *
00935                CALL DLASET( 'F', M, N, ZERO, ZERO, U, LDU )
00936                CALL DBDSDC( 'U', 'I', N, S, WORK( IE ), U, LDU, VT,
00937      $                      LDVT, DUM, IDUM, WORK( NWORK ), IWORK,
00938      $                      INFO )
00939 *
00940 *              Overwrite U by left singular vectors of A and VT
00941 *              by right singular vectors of A
00942 *              (Workspace: need 3*N, prefer 2*N+N*NB)
00943 *
00944                CALL DORMBR( 'Q', 'L', 'N', M, N, N, A, LDA,
00945      $                      WORK( ITAUQ ), U, LDU, WORK( NWORK ),
00946      $                      LWORK-NWORK+1, IERR )
00947                CALL DORMBR( 'P', 'R', 'T', N, N, N, A, LDA,
00948      $                      WORK( ITAUP ), VT, LDVT, WORK( NWORK ),
00949      $                      LWORK-NWORK+1, IERR )
00950             ELSE IF( WNTQA ) THEN
00951 *
00952 *              Perform bidiagonal SVD, computing left singular vectors
00953 *              of bidiagonal matrix in U and computing right singular
00954 *              vectors of bidiagonal matrix in VT
00955 *              (Workspace: need N+BDSPAC)
00956 *
00957                CALL DLASET( 'F', M, M, ZERO, ZERO, U, LDU )
00958                CALL DBDSDC( 'U', 'I', N, S, WORK( IE ), U, LDU, VT,
00959      $                      LDVT, DUM, IDUM, WORK( NWORK ), IWORK,
00960      $                      INFO )
00961 *
00962 *              Set the right corner of U to identity matrix
00963 *
00964                IF( M.GT.N ) THEN
00965                   CALL DLASET( 'F', M-N, M-N, ZERO, ONE, U( N+1, N+1 ),
00966      $                         LDU )
00967                END IF
00968 *
00969 *              Overwrite U by left singular vectors of A and VT
00970 *              by right singular vectors of A
00971 *              (Workspace: need N*N+2*N+M, prefer N*N+2*N+M*NB)
00972 *
00973                CALL DORMBR( 'Q', 'L', 'N', M, M, N, A, LDA,
00974      $                      WORK( ITAUQ ), U, LDU, WORK( NWORK ),
00975      $                      LWORK-NWORK+1, IERR )
00976                CALL DORMBR( 'P', 'R', 'T', N, N, M, A, LDA,
00977      $                      WORK( ITAUP ), VT, LDVT, WORK( NWORK ),
00978      $                      LWORK-NWORK+1, IERR )
00979             END IF
00980 *
00981          END IF
00982 *
00983       ELSE
00984 *
00985 *        A has more columns than rows. If A has sufficiently more
00986 *        columns than rows, first reduce using the LQ decomposition (if
00987 *        sufficient workspace available)
00988 *
00989          IF( N.GE.MNTHR ) THEN
00990 *
00991             IF( WNTQN ) THEN
00992 *
00993 *              Path 1t (N much larger than M, JOBZ='N')
00994 *              No singular vectors to be computed
00995 *
00996                ITAU = 1
00997                NWORK = ITAU + M
00998 *
00999 *              Compute A=L*Q
01000 *              (Workspace: need 2*M, prefer M+M*NB)
01001 *
01002                CALL DGELQF( M, N, A, LDA, WORK( ITAU ), WORK( NWORK ),
01003      $                      LWORK-NWORK+1, IERR )
01004 *
01005 *              Zero out above L
01006 *
01007                CALL DLASET( 'U', M-1, M-1, ZERO, ZERO, A( 1, 2 ), LDA )
01008                IE = 1
01009                ITAUQ = IE + M
01010                ITAUP = ITAUQ + M
01011                NWORK = ITAUP + M
01012 *
01013 *              Bidiagonalize L in A
01014 *              (Workspace: need 4*M, prefer 3*M+2*M*NB)
01015 *
01016                CALL DGEBRD( M, M, A, LDA, S, WORK( IE ), WORK( ITAUQ ),
01017      $                      WORK( ITAUP ), WORK( NWORK ), LWORK-NWORK+1,
01018      $                      IERR )
01019                NWORK = IE + M
01020 *
01021 *              Perform bidiagonal SVD, computing singular values only
01022 *              (Workspace: need M+BDSPAC)
01023 *
01024                CALL DBDSDC( 'U', 'N', M, S, WORK( IE ), DUM, 1, DUM, 1,
01025      $                      DUM, IDUM, WORK( NWORK ), IWORK, INFO )
01026 *
01027             ELSE IF( WNTQO ) THEN
01028 *
01029 *              Path 2t (N much larger than M, JOBZ='O')
01030 *              M right singular vectors to be overwritten on A and
01031 *              M left singular vectors to be computed in U
01032 *
01033                IVT = 1
01034 *
01035 *              IVT is M by M
01036 *
01037                IL = IVT + M*M
01038                IF( LWORK.GE.M*N+M*M+3*M+BDSPAC ) THEN
01039 *
01040 *                 WORK(IL) is M by N
01041 *
01042                   LDWRKL = M
01043                   CHUNK = N
01044                ELSE
01045                   LDWRKL = M
01046                   CHUNK = ( LWORK-M*M ) / M
01047                END IF
01048                ITAU = IL + LDWRKL*M
01049                NWORK = ITAU + M
01050 *
01051 *              Compute A=L*Q
01052 *              (Workspace: need M*M+2*M, prefer M*M+M+M*NB)
01053 *
01054                CALL DGELQF( M, N, A, LDA, WORK( ITAU ), WORK( NWORK ),
01055      $                      LWORK-NWORK+1, IERR )
01056 *
01057 *              Copy L to WORK(IL), zeroing about above it
01058 *
01059                CALL DLACPY( 'L', M, M, A, LDA, WORK( IL ), LDWRKL )
01060                CALL DLASET( 'U', M-1, M-1, ZERO, ZERO,
01061      $                      WORK( IL+LDWRKL ), LDWRKL )
01062 *
01063 *              Generate Q in A
01064 *              (Workspace: need M*M+2*M, prefer M*M+M+M*NB)
01065 *
01066                CALL DORGLQ( M, N, M, A, LDA, WORK( ITAU ),
01067      $                      WORK( NWORK ), LWORK-NWORK+1, IERR )
01068                IE = ITAU
01069                ITAUQ = IE + M
01070                ITAUP = ITAUQ + M
01071                NWORK = ITAUP + M
01072 *
01073 *              Bidiagonalize L in WORK(IL)
01074 *              (Workspace: need M*M+4*M, prefer M*M+3*M+2*M*NB)
01075 *
01076                CALL DGEBRD( M, M, WORK( IL ), LDWRKL, S, WORK( IE ),
01077      $                      WORK( ITAUQ ), WORK( ITAUP ), WORK( NWORK ),
01078      $                      LWORK-NWORK+1, IERR )
01079 *
01080 *              Perform bidiagonal SVD, computing left singular vectors
01081 *              of bidiagonal matrix in U, and computing right singular
01082 *              vectors of bidiagonal matrix in WORK(IVT)
01083 *              (Workspace: need M+M*M+BDSPAC)
01084 *
01085                CALL DBDSDC( 'U', 'I', M, S, WORK( IE ), U, LDU,
01086      $                      WORK( IVT ), M, DUM, IDUM, WORK( NWORK ),
01087      $                      IWORK, INFO )
01088 *
01089 *              Overwrite U by left singular vectors of L and WORK(IVT)
01090 *              by right singular vectors of L
01091 *              (Workspace: need 2*M*M+3*M, prefer 2*M*M+2*M+M*NB)
01092 *
01093                CALL DORMBR( 'Q', 'L', 'N', M, M, M, WORK( IL ), LDWRKL,
01094      $                      WORK( ITAUQ ), U, LDU, WORK( NWORK ),
01095      $                      LWORK-NWORK+1, IERR )
01096                CALL DORMBR( 'P', 'R', 'T', M, M, M, WORK( IL ), LDWRKL,
01097      $                      WORK( ITAUP ), WORK( IVT ), M,
01098      $                      WORK( NWORK ), LWORK-NWORK+1, IERR )
01099 *
01100 *              Multiply right singular vectors of L in WORK(IVT) by Q
01101 *              in A, storing result in WORK(IL) and copying to A
01102 *              (Workspace: need 2*M*M, prefer M*M+M*N)
01103 *
01104                DO 30 I = 1, N, CHUNK
01105                   BLK = MIN( N-I+1, CHUNK )
01106                   CALL DGEMM( 'N', 'N', M, BLK, M, ONE, WORK( IVT ), M,
01107      $                        A( 1, I ), LDA, ZERO, WORK( IL ), LDWRKL )
01108                   CALL DLACPY( 'F', M, BLK, WORK( IL ), LDWRKL,
01109      $                         A( 1, I ), LDA )
01110    30          CONTINUE
01111 *
01112             ELSE IF( WNTQS ) THEN
01113 *
01114 *              Path 3t (N much larger than M, JOBZ='S')
01115 *              M right singular vectors to be computed in VT and
01116 *              M left singular vectors to be computed in U
01117 *
01118                IL = 1
01119 *
01120 *              WORK(IL) is M by M
01121 *
01122                LDWRKL = M
01123                ITAU = IL + LDWRKL*M
01124                NWORK = ITAU + M
01125 *
01126 *              Compute A=L*Q
01127 *              (Workspace: need M*M+2*M, prefer M*M+M+M*NB)
01128 *
01129                CALL DGELQF( M, N, A, LDA, WORK( ITAU ), WORK( NWORK ),
01130      $                      LWORK-NWORK+1, IERR )
01131 *
01132 *              Copy L to WORK(IL), zeroing out above it
01133 *
01134                CALL DLACPY( 'L', M, M, A, LDA, WORK( IL ), LDWRKL )
01135                CALL DLASET( 'U', M-1, M-1, ZERO, ZERO,
01136      $                      WORK( IL+LDWRKL ), LDWRKL )
01137 *
01138 *              Generate Q in A
01139 *              (Workspace: need M*M+2*M, prefer M*M+M+M*NB)
01140 *
01141                CALL DORGLQ( M, N, M, A, LDA, WORK( ITAU ),
01142      $                      WORK( NWORK ), LWORK-NWORK+1, IERR )
01143                IE = ITAU
01144                ITAUQ = IE + M
01145                ITAUP = ITAUQ + M
01146                NWORK = ITAUP + M
01147 *
01148 *              Bidiagonalize L in WORK(IU), copying result to U
01149 *              (Workspace: need M*M+4*M, prefer M*M+3*M+2*M*NB)
01150 *
01151                CALL DGEBRD( M, M, WORK( IL ), LDWRKL, S, WORK( IE ),
01152      $                      WORK( ITAUQ ), WORK( ITAUP ), WORK( NWORK ),
01153      $                      LWORK-NWORK+1, IERR )
01154 *
01155 *              Perform bidiagonal SVD, computing left singular vectors
01156 *              of bidiagonal matrix in U and computing right singular
01157 *              vectors of bidiagonal matrix in VT
01158 *              (Workspace: need M+BDSPAC)
01159 *
01160                CALL DBDSDC( 'U', 'I', M, S, WORK( IE ), U, LDU, VT,
01161      $                      LDVT, DUM, IDUM, WORK( NWORK ), IWORK,
01162      $                      INFO )
01163 *
01164 *              Overwrite U by left singular vectors of L and VT
01165 *              by right singular vectors of L
01166 *              (Workspace: need M*M+3*M, prefer M*M+2*M+M*NB)
01167 *
01168                CALL DORMBR( 'Q', 'L', 'N', M, M, M, WORK( IL ), LDWRKL,
01169      $                      WORK( ITAUQ ), U, LDU, WORK( NWORK ),
01170      $                      LWORK-NWORK+1, IERR )
01171                CALL DORMBR( 'P', 'R', 'T', M, M, M, WORK( IL ), LDWRKL,
01172      $                      WORK( ITAUP ), VT, LDVT, WORK( NWORK ),
01173      $                      LWORK-NWORK+1, IERR )
01174 *
01175 *              Multiply right singular vectors of L in WORK(IL) by
01176 *              Q in A, storing result in VT
01177 *              (Workspace: need M*M)
01178 *
01179                CALL DLACPY( 'F', M, M, VT, LDVT, WORK( IL ), LDWRKL )
01180                CALL DGEMM( 'N', 'N', M, N, M, ONE, WORK( IL ), LDWRKL,
01181      $                     A, LDA, ZERO, VT, LDVT )
01182 *
01183             ELSE IF( WNTQA ) THEN
01184 *
01185 *              Path 4t (N much larger than M, JOBZ='A')
01186 *              N right singular vectors to be computed in VT and
01187 *              M left singular vectors to be computed in U
01188 *
01189                IVT = 1
01190 *
01191 *              WORK(IVT) is M by M
01192 *
01193                LDWKVT = M
01194                ITAU = IVT + LDWKVT*M
01195                NWORK = ITAU + M
01196 *
01197 *              Compute A=L*Q, copying result to VT
01198 *              (Workspace: need M*M+2*M, prefer M*M+M+M*NB)
01199 *
01200                CALL DGELQF( M, N, A, LDA, WORK( ITAU ), WORK( NWORK ),
01201      $                      LWORK-NWORK+1, IERR )
01202                CALL DLACPY( 'U', M, N, A, LDA, VT, LDVT )
01203 *
01204 *              Generate Q in VT
01205 *              (Workspace: need M*M+2*M, prefer M*M+M+M*NB)
01206 *
01207                CALL DORGLQ( N, N, M, VT, LDVT, WORK( ITAU ),
01208      $                      WORK( NWORK ), LWORK-NWORK+1, IERR )
01209 *
01210 *              Produce L in A, zeroing out other entries
01211 *
01212                CALL DLASET( 'U', M-1, M-1, ZERO, ZERO, A( 1, 2 ), LDA )
01213                IE = ITAU
01214                ITAUQ = IE + M
01215                ITAUP = ITAUQ + M
01216                NWORK = ITAUP + M
01217 *
01218 *              Bidiagonalize L in A
01219 *              (Workspace: need M*M+4*M, prefer M*M+3*M+2*M*NB)
01220 *
01221                CALL DGEBRD( M, M, A, LDA, S, WORK( IE ), WORK( ITAUQ ),
01222      $                      WORK( ITAUP ), WORK( NWORK ), LWORK-NWORK+1,
01223      $                      IERR )
01224 *
01225 *              Perform bidiagonal SVD, computing left singular vectors
01226 *              of bidiagonal matrix in U and computing right singular
01227 *              vectors of bidiagonal matrix in WORK(IVT)
01228 *              (Workspace: need M+M*M+BDSPAC)
01229 *
01230                CALL DBDSDC( 'U', 'I', M, S, WORK( IE ), U, LDU,
01231      $                      WORK( IVT ), LDWKVT, DUM, IDUM,
01232      $                      WORK( NWORK ), IWORK, INFO )
01233 *
01234 *              Overwrite U by left singular vectors of L and WORK(IVT)
01235 *              by right singular vectors of L
01236 *              (Workspace: need M*M+3*M, prefer M*M+2*M+M*NB)
01237 *
01238                CALL DORMBR( 'Q', 'L', 'N', M, M, M, A, LDA,
01239      $                      WORK( ITAUQ ), U, LDU, WORK( NWORK ),
01240      $                      LWORK-NWORK+1, IERR )
01241                CALL DORMBR( 'P', 'R', 'T', M, M, M, A, LDA,
01242      $                      WORK( ITAUP ), WORK( IVT ), LDWKVT,
01243      $                      WORK( NWORK ), LWORK-NWORK+1, IERR )
01244 *
01245 *              Multiply right singular vectors of L in WORK(IVT) by
01246 *              Q in VT, storing result in A
01247 *              (Workspace: need M*M)
01248 *
01249                CALL DGEMM( 'N', 'N', M, N, M, ONE, WORK( IVT ), LDWKVT,
01250      $                     VT, LDVT, ZERO, A, LDA )
01251 *
01252 *              Copy right singular vectors of A from A to VT
01253 *
01254                CALL DLACPY( 'F', M, N, A, LDA, VT, LDVT )
01255 *
01256             END IF
01257 *
01258          ELSE
01259 *
01260 *           N .LT. MNTHR
01261 *
01262 *           Path 5t (N greater than M, but not much larger)
01263 *           Reduce to bidiagonal form without LQ decomposition
01264 *
01265             IE = 1
01266             ITAUQ = IE + M
01267             ITAUP = ITAUQ + M
01268             NWORK = ITAUP + M
01269 *
01270 *           Bidiagonalize A
01271 *           (Workspace: need 3*M+N, prefer 3*M+(M+N)*NB)
01272 *
01273             CALL DGEBRD( M, N, A, LDA, S, WORK( IE ), WORK( ITAUQ ),
01274      $                   WORK( ITAUP ), WORK( NWORK ), LWORK-NWORK+1,
01275      $                   IERR )
01276             IF( WNTQN ) THEN
01277 *
01278 *              Perform bidiagonal SVD, only computing singular values
01279 *              (Workspace: need M+BDSPAC)
01280 *
01281                CALL DBDSDC( 'L', 'N', M, S, WORK( IE ), DUM, 1, DUM, 1,
01282      $                      DUM, IDUM, WORK( NWORK ), IWORK, INFO )
01283             ELSE IF( WNTQO ) THEN
01284                LDWKVT = M
01285                IVT = NWORK
01286                IF( LWORK.GE.M*N+3*M+BDSPAC ) THEN
01287 *
01288 *                 WORK( IVT ) is M by N
01289 *
01290                   CALL DLASET( 'F', M, N, ZERO, ZERO, WORK( IVT ),
01291      $                         LDWKVT )
01292                   NWORK = IVT + LDWKVT*N
01293                ELSE
01294 *
01295 *                 WORK( IVT ) is M by M
01296 *
01297                   NWORK = IVT + LDWKVT*M
01298                   IL = NWORK
01299 *
01300 *                 WORK(IL) is M by CHUNK
01301 *
01302                   CHUNK = ( LWORK-M*M-3*M ) / M
01303                END IF
01304 *
01305 *              Perform bidiagonal SVD, computing left singular vectors
01306 *              of bidiagonal matrix in U and computing right singular
01307 *              vectors of bidiagonal matrix in WORK(IVT)
01308 *              (Workspace: need M*M+BDSPAC)
01309 *
01310                CALL DBDSDC( 'L', 'I', M, S, WORK( IE ), U, LDU,
01311      $                      WORK( IVT ), LDWKVT, DUM, IDUM,
01312      $                      WORK( NWORK ), IWORK, INFO )
01313 *
01314 *              Overwrite U by left singular vectors of A
01315 *              (Workspace: need M*M+2*M, prefer M*M+M+M*NB)
01316 *
01317                CALL DORMBR( 'Q', 'L', 'N', M, M, N, A, LDA,
01318      $                      WORK( ITAUQ ), U, LDU, WORK( NWORK ),
01319      $                      LWORK-NWORK+1, IERR )
01320 *
01321                IF( LWORK.GE.M*N+3*M+BDSPAC ) THEN
01322 *
01323 *                 Overwrite WORK(IVT) by left singular vectors of A
01324 *                 (Workspace: need M*M+2*M, prefer M*M+M+M*NB)
01325 *
01326                   CALL DORMBR( 'P', 'R', 'T', M, N, M, A, LDA,
01327      $                         WORK( ITAUP ), WORK( IVT ), LDWKVT,
01328      $                         WORK( NWORK ), LWORK-NWORK+1, IERR )
01329 *
01330 *                 Copy right singular vectors of A from WORK(IVT) to A
01331 *
01332                   CALL DLACPY( 'F', M, N, WORK( IVT ), LDWKVT, A, LDA )
01333                ELSE
01334 *
01335 *                 Generate P**T in A
01336 *                 (Workspace: need M*M+2*M, prefer M*M+M+M*NB)
01337 *
01338                   CALL DORGBR( 'P', M, N, M, A, LDA, WORK( ITAUP ),
01339      $                         WORK( NWORK ), LWORK-NWORK+1, IERR )
01340 *
01341 *                 Multiply Q in A by right singular vectors of
01342 *                 bidiagonal matrix in WORK(IVT), storing result in
01343 *                 WORK(IL) and copying to A
01344 *                 (Workspace: need 2*M*M, prefer M*M+M*N)
01345 *
01346                   DO 40 I = 1, N, CHUNK
01347                      BLK = MIN( N-I+1, CHUNK )
01348                      CALL DGEMM( 'N', 'N', M, BLK, M, ONE, WORK( IVT ),
01349      $                           LDWKVT, A( 1, I ), LDA, ZERO,
01350      $                           WORK( IL ), M )
01351                      CALL DLACPY( 'F', M, BLK, WORK( IL ), M, A( 1, I ),
01352      $                            LDA )
01353    40             CONTINUE
01354                END IF
01355             ELSE IF( WNTQS ) THEN
01356 *
01357 *              Perform bidiagonal SVD, computing left singular vectors
01358 *              of bidiagonal matrix in U and computing right singular
01359 *              vectors of bidiagonal matrix in VT
01360 *              (Workspace: need M+BDSPAC)
01361 *
01362                CALL DLASET( 'F', M, N, ZERO, ZERO, VT, LDVT )
01363                CALL DBDSDC( 'L', 'I', M, S, WORK( IE ), U, LDU, VT,
01364      $                      LDVT, DUM, IDUM, WORK( NWORK ), IWORK,
01365      $                      INFO )
01366 *
01367 *              Overwrite U by left singular vectors of A and VT
01368 *              by right singular vectors of A
01369 *              (Workspace: need 3*M, prefer 2*M+M*NB)
01370 *
01371                CALL DORMBR( 'Q', 'L', 'N', M, M, N, A, LDA,
01372      $                      WORK( ITAUQ ), U, LDU, WORK( NWORK ),
01373      $                      LWORK-NWORK+1, IERR )
01374                CALL DORMBR( 'P', 'R', 'T', M, N, M, A, LDA,
01375      $                      WORK( ITAUP ), VT, LDVT, WORK( NWORK ),
01376      $                      LWORK-NWORK+1, IERR )
01377             ELSE IF( WNTQA ) THEN
01378 *
01379 *              Perform bidiagonal SVD, computing left singular vectors
01380 *              of bidiagonal matrix in U and computing right singular
01381 *              vectors of bidiagonal matrix in VT
01382 *              (Workspace: need M+BDSPAC)
01383 *
01384                CALL DLASET( 'F', N, N, ZERO, ZERO, VT, LDVT )
01385                CALL DBDSDC( 'L', 'I', M, S, WORK( IE ), U, LDU, VT,
01386      $                      LDVT, DUM, IDUM, WORK( NWORK ), IWORK,
01387      $                      INFO )
01388 *
01389 *              Set the right corner of VT to identity matrix
01390 *
01391                IF( N.GT.M ) THEN
01392                   CALL DLASET( 'F', N-M, N-M, ZERO, ONE, VT( M+1, M+1 ),
01393      $                         LDVT )
01394                END IF
01395 *
01396 *              Overwrite U by left singular vectors of A and VT
01397 *              by right singular vectors of A
01398 *              (Workspace: need 2*M+N, prefer 2*M+N*NB)
01399 *
01400                CALL DORMBR( 'Q', 'L', 'N', M, M, N, A, LDA,
01401      $                      WORK( ITAUQ ), U, LDU, WORK( NWORK ),
01402      $                      LWORK-NWORK+1, IERR )
01403                CALL DORMBR( 'P', 'R', 'T', N, N, M, A, LDA,
01404      $                      WORK( ITAUP ), VT, LDVT, WORK( NWORK ),
01405      $                      LWORK-NWORK+1, IERR )
01406             END IF
01407 *
01408          END IF
01409 *
01410       END IF
01411 *
01412 *     Undo scaling if necessary
01413 *
01414       IF( ISCL.EQ.1 ) THEN
01415          IF( ANRM.GT.BIGNUM )
01416      $      CALL DLASCL( 'G', 0, 0, BIGNUM, ANRM, MINMN, 1, S, MINMN,
01417      $                   IERR )
01418          IF( ANRM.LT.SMLNUM )
01419      $      CALL DLASCL( 'G', 0, 0, SMLNUM, ANRM, MINMN, 1, S, MINMN,
01420      $                   IERR )
01421       END IF
01422 *
01423 *     Return optimal workspace in WORK(1)
01424 *
01425       WORK( 1 ) = MAXWRK
01426 *
01427       RETURN
01428 *
01429 *     End of DGESDD
01430 *
01431       END
 All Files Functions