LAPACK  3.4.1
LAPACK: Linear Algebra PACKage
cgesdd.f
Go to the documentation of this file.
00001 *> \brief \b CGESDD
00002 *
00003 *  =========== DOCUMENTATION ===========
00004 *
00005 * Online html documentation available at 
00006 *            http://www.netlib.org/lapack/explore-html/ 
00007 *
00008 *> \htmlonly
00009 *> Download CGESDD + dependencies 
00010 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/cgesdd.f"> 
00011 *> [TGZ]</a> 
00012 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/cgesdd.f"> 
00013 *> [ZIP]</a> 
00014 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/cgesdd.f"> 
00015 *> [TXT]</a>
00016 *> \endhtmlonly 
00017 *
00018 *  Definition:
00019 *  ===========
00020 *
00021 *       SUBROUTINE CGESDD( JOBZ, M, N, A, LDA, S, U, LDU, VT, LDVT,
00022 *                          WORK, LWORK, RWORK, 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 *       REAL               RWORK( * ), S( * )
00031 *       COMPLEX            A( LDA, * ), U( LDU, * ), VT( LDVT, * ),
00032 *      $                   WORK( * )
00033 *       ..
00034 *  
00035 *
00036 *> \par Purpose:
00037 *  =============
00038 *>
00039 *> \verbatim
00040 *>
00041 *> CGESDD computes the singular value decomposition (SVD) of a complex
00042 *> M-by-N matrix A, optionally computing the left and/or right singular
00043 *> vectors, by using divide-and-conquer method. The SVD is written
00044 *>
00045 *>      A = U * SIGMA * conjugate-transpose(V)
00046 *>
00047 *> where SIGMA is an M-by-N matrix which is zero except for its
00048 *> min(m,n) diagonal elements, U is an M-by-M unitary matrix, and
00049 *> V is an N-by-N unitary matrix.  The diagonal elements of SIGMA
00050 *> are the singular values of A; they are real and non-negative, and
00051 *> are returned in descending order.  The first min(m,n) columns of
00052 *> U and V are the left and right singular vectors of A.
00053 *>
00054 *> Note that the routine returns VT = V**H, not V.
00055 *>
00056 *> The divide and conquer algorithm makes very mild assumptions about
00057 *> floating point arithmetic. It will work on machines with a guard
00058 *> digit in add/subtract, or on those binary machines without guard
00059 *> digits which subtract like the Cray X-MP, Cray Y-MP, Cray C-90, or
00060 *> Cray-2. It could conceivably fail on hexadecimal or decimal machines
00061 *> without guard digits, but we know of none.
00062 *> \endverbatim
00063 *
00064 *  Arguments:
00065 *  ==========
00066 *
00067 *> \param[in] JOBZ
00068 *> \verbatim
00069 *>          JOBZ is CHARACTER*1
00070 *>          Specifies options for computing all or part of the matrix U:
00071 *>          = 'A':  all M columns of U and all N rows of V**H are
00072 *>                  returned in the arrays U and VT;
00073 *>          = 'S':  the first min(M,N) columns of U and the first
00074 *>                  min(M,N) rows of V**H are returned in the arrays U
00075 *>                  and VT;
00076 *>          = 'O':  If M >= N, the first N columns of U are overwritten
00077 *>                  in the array A and all rows of V**H are returned in
00078 *>                  the array VT;
00079 *>                  otherwise, all columns of U are returned in the
00080 *>                  array U and the first M rows of V**H are overwritten
00081 *>                  in the array A;
00082 *>          = 'N':  no columns of U or rows of V**H are computed.
00083 *> \endverbatim
00084 *>
00085 *> \param[in] M
00086 *> \verbatim
00087 *>          M is INTEGER
00088 *>          The number of rows of the input matrix A.  M >= 0.
00089 *> \endverbatim
00090 *>
00091 *> \param[in] N
00092 *> \verbatim
00093 *>          N is INTEGER
00094 *>          The number of columns of the input matrix A.  N >= 0.
00095 *> \endverbatim
00096 *>
00097 *> \param[in,out] A
00098 *> \verbatim
00099 *>          A is COMPLEX array, dimension (LDA,N)
00100 *>          On entry, the M-by-N matrix A.
00101 *>          On exit,
00102 *>          if JOBZ = 'O',  A is overwritten with the first N columns
00103 *>                          of U (the left singular vectors, stored
00104 *>                          columnwise) if M >= N;
00105 *>                          A is overwritten with the first M rows
00106 *>                          of V**H (the right singular vectors, stored
00107 *>                          rowwise) otherwise.
00108 *>          if JOBZ .ne. 'O', the contents of A are destroyed.
00109 *> \endverbatim
00110 *>
00111 *> \param[in] LDA
00112 *> \verbatim
00113 *>          LDA is INTEGER
00114 *>          The leading dimension of the array A.  LDA >= max(1,M).
00115 *> \endverbatim
00116 *>
00117 *> \param[out] S
00118 *> \verbatim
00119 *>          S is REAL array, dimension (min(M,N))
00120 *>          The singular values of A, sorted so that S(i) >= S(i+1).
00121 *> \endverbatim
00122 *>
00123 *> \param[out] U
00124 *> \verbatim
00125 *>          U is COMPLEX array, dimension (LDU,UCOL)
00126 *>          UCOL = M if JOBZ = 'A' or JOBZ = 'O' and M < N;
00127 *>          UCOL = min(M,N) if JOBZ = 'S'.
00128 *>          If JOBZ = 'A' or JOBZ = 'O' and M < N, U contains the M-by-M
00129 *>          unitary matrix U;
00130 *>          if JOBZ = 'S', U contains the first min(M,N) columns of U
00131 *>          (the left singular vectors, stored columnwise);
00132 *>          if JOBZ = 'O' and M >= N, or JOBZ = 'N', U is not referenced.
00133 *> \endverbatim
00134 *>
00135 *> \param[in] LDU
00136 *> \verbatim
00137 *>          LDU is INTEGER
00138 *>          The leading dimension of the array U.  LDU >= 1; if
00139 *>          JOBZ = 'S' or 'A' or JOBZ = 'O' and M < N, LDU >= M.
00140 *> \endverbatim
00141 *>
00142 *> \param[out] VT
00143 *> \verbatim
00144 *>          VT is COMPLEX array, dimension (LDVT,N)
00145 *>          If JOBZ = 'A' or JOBZ = 'O' and M >= N, VT contains the
00146 *>          N-by-N unitary matrix V**H;
00147 *>          if JOBZ = 'S', VT contains the first min(M,N) rows of
00148 *>          V**H (the right singular vectors, stored rowwise);
00149 *>          if JOBZ = 'O' and M < N, or JOBZ = 'N', VT is not referenced.
00150 *> \endverbatim
00151 *>
00152 *> \param[in] LDVT
00153 *> \verbatim
00154 *>          LDVT is INTEGER
00155 *>          The leading dimension of the array VT.  LDVT >= 1; if
00156 *>          JOBZ = 'A' or JOBZ = 'O' and M >= N, LDVT >= N;
00157 *>          if JOBZ = 'S', LDVT >= min(M,N).
00158 *> \endverbatim
00159 *>
00160 *> \param[out] WORK
00161 *> \verbatim
00162 *>          WORK is COMPLEX array, dimension (MAX(1,LWORK))
00163 *>          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
00164 *> \endverbatim
00165 *>
00166 *> \param[in] LWORK
00167 *> \verbatim
00168 *>          LWORK is INTEGER
00169 *>          The dimension of the array WORK. LWORK >= 1.
00170 *>          if JOBZ = 'N', LWORK >= 2*min(M,N)+max(M,N).
00171 *>          if JOBZ = 'O',
00172 *>                LWORK >= 2*min(M,N)*min(M,N)+2*min(M,N)+max(M,N).
00173 *>          if JOBZ = 'S' or 'A',
00174 *>                LWORK >= min(M,N)*min(M,N)+2*min(M,N)+max(M,N).
00175 *>          For good performance, LWORK should generally be larger.
00176 *>
00177 *>          If LWORK = -1, a workspace query is assumed.  The optimal
00178 *>          size for the WORK array is calculated and stored in WORK(1),
00179 *>          and no other work except argument checking is performed.
00180 *> \endverbatim
00181 *>
00182 *> \param[out] RWORK
00183 *> \verbatim
00184 *>          RWORK is REAL array, dimension (MAX(1,LRWORK))
00185 *>          If JOBZ = 'N', LRWORK >= 5*min(M,N).
00186 *>          Otherwise, 
00187 *>          LRWORK >= min(M,N)*max(5*min(M,N)+7,2*max(M,N)+2*min(M,N)+1)
00188 *> \endverbatim
00189 *>
00190 *> \param[out] IWORK
00191 *> \verbatim
00192 *>          IWORK is INTEGER array, dimension (8*min(M,N))
00193 *> \endverbatim
00194 *>
00195 *> \param[out] INFO
00196 *> \verbatim
00197 *>          INFO is INTEGER
00198 *>          = 0:  successful exit.
00199 *>          < 0:  if INFO = -i, the i-th argument had an illegal value.
00200 *>          > 0:  The updating process of SBDSDC did not converge.
00201 *> \endverbatim
00202 *
00203 *  Authors:
00204 *  ========
00205 *
00206 *> \author Univ. of Tennessee 
00207 *> \author Univ. of California Berkeley 
00208 *> \author Univ. of Colorado Denver 
00209 *> \author NAG Ltd. 
00210 *
00211 *> \date November 2011
00212 *
00213 *> \ingroup complexGEsing
00214 *
00215 *> \par Contributors:
00216 *  ==================
00217 *>
00218 *>     Ming Gu and Huan Ren, Computer Science Division, University of
00219 *>     California at Berkeley, USA
00220 *>
00221 *  =====================================================================
00222       SUBROUTINE CGESDD( JOBZ, M, N, A, LDA, S, U, LDU, VT, LDVT,
00223      $                   WORK, LWORK, RWORK, IWORK, INFO )
00224 *
00225 *  -- LAPACK driver routine (version 3.4.0) --
00226 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
00227 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
00228 *     November 2011
00229 *
00230 *     .. Scalar Arguments ..
00231       CHARACTER          JOBZ
00232       INTEGER            INFO, LDA, LDU, LDVT, LWORK, M, N
00233 *     ..
00234 *     .. Array Arguments ..
00235       INTEGER            IWORK( * )
00236       REAL               RWORK( * ), S( * )
00237       COMPLEX            A( LDA, * ), U( LDU, * ), VT( LDVT, * ),
00238      $                   WORK( * )
00239 *     ..
00240 *
00241 *  =====================================================================
00242 *
00243 *     .. Parameters ..
00244       INTEGER            LQUERV
00245       PARAMETER          ( LQUERV = -1 )
00246       COMPLEX            CZERO, CONE
00247       PARAMETER          ( CZERO = ( 0.0E+0, 0.0E+0 ),
00248      $                   CONE = ( 1.0E+0, 0.0E+0 ) )
00249       REAL               ZERO, ONE
00250       PARAMETER          ( ZERO = 0.0E+0, ONE = 1.0E+0 )
00251 *     ..
00252 *     .. Local Scalars ..
00253       LOGICAL            WNTQA, WNTQAS, WNTQN, WNTQO, WNTQS
00254       INTEGER            BLK, CHUNK, I, IE, IERR, IL, IR, IRU, IRVT,
00255      $                   ISCL, ITAU, ITAUP, ITAUQ, IU, IVT, LDWKVT,
00256      $                   LDWRKL, LDWRKR, LDWRKU, MAXWRK, MINMN, MINWRK,
00257      $                   MNTHR1, MNTHR2, NRWORK, NWORK, WRKBL
00258       REAL               ANRM, BIGNUM, EPS, SMLNUM
00259 *     ..
00260 *     .. Local Arrays ..
00261       INTEGER            IDUM( 1 )
00262       REAL               DUM( 1 )
00263 *     ..
00264 *     .. External Subroutines ..
00265       EXTERNAL           CGEBRD, CGELQF, CGEMM, CGEQRF, CLACP2, CLACPY,
00266      $                   CLACRM, CLARCM, CLASCL, CLASET, CUNGBR, CUNGLQ,
00267      $                   CUNGQR, CUNMBR, SBDSDC, SLASCL, XERBLA
00268 *     ..
00269 *     .. External Functions ..
00270       LOGICAL            LSAME
00271       INTEGER            ILAENV
00272       REAL               CLANGE, SLAMCH
00273       EXTERNAL           CLANGE, SLAMCH, ILAENV, LSAME
00274 *     ..
00275 *     .. Intrinsic Functions ..
00276       INTRINSIC          INT, MAX, MIN, SQRT
00277 *     ..
00278 *     .. Executable Statements ..
00279 *
00280 *     Test the input arguments
00281 *
00282       INFO = 0
00283       MINMN = MIN( M, N )
00284       MNTHR1 = INT( MINMN*17.0 / 9.0 )
00285       MNTHR2 = INT( MINMN*5.0 / 3.0 )
00286       WNTQA = LSAME( JOBZ, 'A' )
00287       WNTQS = LSAME( JOBZ, 'S' )
00288       WNTQAS = WNTQA .OR. WNTQS
00289       WNTQO = LSAME( JOBZ, 'O' )
00290       WNTQN = LSAME( JOBZ, 'N' )
00291       MINWRK = 1
00292       MAXWRK = 1
00293 *
00294       IF( .NOT.( WNTQA .OR. WNTQS .OR. WNTQO .OR. WNTQN ) ) THEN
00295          INFO = -1
00296       ELSE IF( M.LT.0 ) THEN
00297          INFO = -2
00298       ELSE IF( N.LT.0 ) THEN
00299          INFO = -3
00300       ELSE IF( LDA.LT.MAX( 1, M ) ) THEN
00301          INFO = -5
00302       ELSE IF( LDU.LT.1 .OR. ( WNTQAS .AND. LDU.LT.M ) .OR.
00303      $         ( WNTQO .AND. M.LT.N .AND. LDU.LT.M ) ) THEN
00304          INFO = -8
00305       ELSE IF( LDVT.LT.1 .OR. ( WNTQA .AND. LDVT.LT.N ) .OR.
00306      $         ( WNTQS .AND. LDVT.LT.MINMN ) .OR.
00307      $         ( WNTQO .AND. M.GE.N .AND. LDVT.LT.N ) ) THEN
00308          INFO = -10
00309       END IF
00310 *
00311 *     Compute workspace
00312 *      (Note: Comments in the code beginning "Workspace:" describe the
00313 *       minimal amount of workspace needed at that point in the code,
00314 *       as well as the preferred amount for good performance.
00315 *       CWorkspace refers to complex workspace, and RWorkspace to
00316 *       real workspace. NB refers to the optimal block size for the
00317 *       immediately following subroutine, as returned by ILAENV.)
00318 *
00319       IF( INFO.EQ.0 .AND. M.GT.0 .AND. N.GT.0 ) THEN
00320          IF( M.GE.N ) THEN
00321 *
00322 *           There is no complex work space needed for bidiagonal SVD
00323 *           The real work space needed for bidiagonal SVD is BDSPAC
00324 *           for computing singular values and singular vectors; BDSPAN
00325 *           for computing singular values only.
00326 *           BDSPAC = 5*N*N + 7*N
00327 *           BDSPAN = MAX(7*N+4, 3*N+2+SMLSIZ*(SMLSIZ+8))
00328 *
00329             IF( M.GE.MNTHR1 ) THEN
00330                IF( WNTQN ) THEN
00331 *
00332 *                 Path 1 (M much larger than N, JOBZ='N')
00333 *
00334                   MAXWRK = N + N*ILAENV( 1, 'CGEQRF', ' ', M, N, -1,
00335      $                     -1 )
00336                   MAXWRK = MAX( MAXWRK, 2*N+2*N*
00337      $                     ILAENV( 1, 'CGEBRD', ' ', N, N, -1, -1 ) )
00338                   MINWRK = 3*N
00339                ELSE IF( WNTQO ) THEN
00340 *
00341 *                 Path 2 (M much larger than N, JOBZ='O')
00342 *
00343                   WRKBL = N + N*ILAENV( 1, 'CGEQRF', ' ', M, N, -1, -1 )
00344                   WRKBL = MAX( WRKBL, N+N*ILAENV( 1, 'CUNGQR', ' ', M,
00345      $                    N, N, -1 ) )
00346                   WRKBL = MAX( WRKBL, 2*N+2*N*
00347      $                    ILAENV( 1, 'CGEBRD', ' ', N, N, -1, -1 ) )
00348                   WRKBL = MAX( WRKBL, 2*N+N*
00349      $                    ILAENV( 1, 'CUNMBR', 'QLN', N, N, N, -1 ) )
00350                   WRKBL = MAX( WRKBL, 2*N+N*
00351      $                    ILAENV( 1, 'CUNMBR', 'PRC', N, N, N, -1 ) )
00352                   MAXWRK = M*N + N*N + WRKBL
00353                   MINWRK = 2*N*N + 3*N
00354                ELSE IF( WNTQS ) THEN
00355 *
00356 *                 Path 3 (M much larger than N, JOBZ='S')
00357 *
00358                   WRKBL = N + N*ILAENV( 1, 'CGEQRF', ' ', M, N, -1, -1 )
00359                   WRKBL = MAX( WRKBL, N+N*ILAENV( 1, 'CUNGQR', ' ', M,
00360      $                    N, N, -1 ) )
00361                   WRKBL = MAX( WRKBL, 2*N+2*N*
00362      $                    ILAENV( 1, 'CGEBRD', ' ', N, N, -1, -1 ) )
00363                   WRKBL = MAX( WRKBL, 2*N+N*
00364      $                    ILAENV( 1, 'CUNMBR', 'QLN', N, N, N, -1 ) )
00365                   WRKBL = MAX( WRKBL, 2*N+N*
00366      $                    ILAENV( 1, 'CUNMBR', 'PRC', N, N, N, -1 ) )
00367                   MAXWRK = N*N + WRKBL
00368                   MINWRK = N*N + 3*N
00369                ELSE IF( WNTQA ) THEN
00370 *
00371 *                 Path 4 (M much larger than N, JOBZ='A')
00372 *
00373                   WRKBL = N + N*ILAENV( 1, 'CGEQRF', ' ', M, N, -1, -1 )
00374                   WRKBL = MAX( WRKBL, N+M*ILAENV( 1, 'CUNGQR', ' ', M,
00375      $                    M, N, -1 ) )
00376                   WRKBL = MAX( WRKBL, 2*N+2*N*
00377      $                    ILAENV( 1, 'CGEBRD', ' ', N, N, -1, -1 ) )
00378                   WRKBL = MAX( WRKBL, 2*N+N*
00379      $                    ILAENV( 1, 'CUNMBR', 'QLN', N, N, N, -1 ) )
00380                   WRKBL = MAX( WRKBL, 2*N+N*
00381      $                    ILAENV( 1, 'CUNMBR', 'PRC', N, N, N, -1 ) )
00382                   MAXWRK = N*N + WRKBL
00383                   MINWRK = N*N + 2*N + M
00384                END IF
00385             ELSE IF( M.GE.MNTHR2 ) THEN
00386 *
00387 *              Path 5 (M much larger than N, but not as much as MNTHR1)
00388 *
00389                MAXWRK = 2*N + ( M+N )*ILAENV( 1, 'CGEBRD', ' ', M, N,
00390      $                  -1, -1 )
00391                MINWRK = 2*N + M
00392                IF( WNTQO ) THEN
00393                   MAXWRK = MAX( MAXWRK, 2*N+N*
00394      $                     ILAENV( 1, 'CUNGBR', 'P', N, N, N, -1 ) )
00395                   MAXWRK = MAX( MAXWRK, 2*N+N*
00396      $                     ILAENV( 1, 'CUNGBR', 'Q', M, N, N, -1 ) )
00397                   MAXWRK = MAXWRK + M*N
00398                   MINWRK = MINWRK + N*N
00399                ELSE IF( WNTQS ) THEN
00400                   MAXWRK = MAX( MAXWRK, 2*N+N*
00401      $                     ILAENV( 1, 'CUNGBR', 'P', N, N, N, -1 ) )
00402                   MAXWRK = MAX( MAXWRK, 2*N+N*
00403      $                     ILAENV( 1, 'CUNGBR', 'Q', M, N, N, -1 ) )
00404                ELSE IF( WNTQA ) THEN
00405                   MAXWRK = MAX( MAXWRK, 2*N+N*
00406      $                     ILAENV( 1, 'CUNGBR', 'P', N, N, N, -1 ) )
00407                   MAXWRK = MAX( MAXWRK, 2*N+M*
00408      $                     ILAENV( 1, 'CUNGBR', 'Q', M, M, N, -1 ) )
00409                END IF
00410             ELSE
00411 *
00412 *              Path 6 (M at least N, but not much larger)
00413 *
00414                MAXWRK = 2*N + ( M+N )*ILAENV( 1, 'CGEBRD', ' ', M, N,
00415      $                  -1, -1 )
00416                MINWRK = 2*N + M
00417                IF( WNTQO ) THEN
00418                   MAXWRK = MAX( MAXWRK, 2*N+N*
00419      $                     ILAENV( 1, 'CUNMBR', 'PRC', N, N, N, -1 ) )
00420                   MAXWRK = MAX( MAXWRK, 2*N+N*
00421      $                     ILAENV( 1, 'CUNMBR', 'QLN', M, N, N, -1 ) )
00422                   MAXWRK = MAXWRK + M*N
00423                   MINWRK = MINWRK + N*N
00424                ELSE IF( WNTQS ) THEN
00425                   MAXWRK = MAX( MAXWRK, 2*N+N*
00426      $                     ILAENV( 1, 'CUNMBR', 'PRC', N, N, N, -1 ) )
00427                   MAXWRK = MAX( MAXWRK, 2*N+N*
00428      $                     ILAENV( 1, 'CUNMBR', 'QLN', M, N, N, -1 ) )
00429                ELSE IF( WNTQA ) THEN
00430                   MAXWRK = MAX( MAXWRK, 2*N+N*
00431      $                     ILAENV( 1, 'CUNGBR', 'PRC', N, N, N, -1 ) )
00432                   MAXWRK = MAX( MAXWRK, 2*N+M*
00433      $                     ILAENV( 1, 'CUNGBR', 'QLN', M, M, N, -1 ) )
00434                END IF
00435             END IF
00436          ELSE
00437 *
00438 *           There is no complex work space needed for bidiagonal SVD
00439 *           The real work space needed for bidiagonal SVD is BDSPAC
00440 *           for computing singular values and singular vectors; BDSPAN
00441 *           for computing singular values only.
00442 *           BDSPAC = 5*M*M + 7*M
00443 *           BDSPAN = MAX(7*M+4, 3*M+2+SMLSIZ*(SMLSIZ+8))
00444 *
00445             IF( N.GE.MNTHR1 ) THEN
00446                IF( WNTQN ) THEN
00447 *
00448 *                 Path 1t (N much larger than M, JOBZ='N')
00449 *
00450                   MAXWRK = M + M*ILAENV( 1, 'CGELQF', ' ', M, N, -1,
00451      $                     -1 )
00452                   MAXWRK = MAX( MAXWRK, 2*M+2*M*
00453      $                     ILAENV( 1, 'CGEBRD', ' ', M, M, -1, -1 ) )
00454                   MINWRK = 3*M
00455                ELSE IF( WNTQO ) THEN
00456 *
00457 *                 Path 2t (N much larger than M, JOBZ='O')
00458 *
00459                   WRKBL = M + M*ILAENV( 1, 'CGELQF', ' ', M, N, -1, -1 )
00460                   WRKBL = MAX( WRKBL, M+M*ILAENV( 1, 'CUNGLQ', ' ', M,
00461      $                    N, M, -1 ) )
00462                   WRKBL = MAX( WRKBL, 2*M+2*M*
00463      $                    ILAENV( 1, 'CGEBRD', ' ', M, M, -1, -1 ) )
00464                   WRKBL = MAX( WRKBL, 2*M+M*
00465      $                    ILAENV( 1, 'CUNMBR', 'PRC', M, M, M, -1 ) )
00466                   WRKBL = MAX( WRKBL, 2*M+M*
00467      $                    ILAENV( 1, 'CUNMBR', 'QLN', M, M, M, -1 ) )
00468                   MAXWRK = M*N + M*M + WRKBL
00469                   MINWRK = 2*M*M + 3*M
00470                ELSE IF( WNTQS ) THEN
00471 *
00472 *                 Path 3t (N much larger than M, JOBZ='S')
00473 *
00474                   WRKBL = M + M*ILAENV( 1, 'CGELQF', ' ', M, N, -1, -1 )
00475                   WRKBL = MAX( WRKBL, M+M*ILAENV( 1, 'CUNGLQ', ' ', M,
00476      $                    N, M, -1 ) )
00477                   WRKBL = MAX( WRKBL, 2*M+2*M*
00478      $                    ILAENV( 1, 'CGEBRD', ' ', M, M, -1, -1 ) )
00479                   WRKBL = MAX( WRKBL, 2*M+M*
00480      $                    ILAENV( 1, 'CUNMBR', 'PRC', M, M, M, -1 ) )
00481                   WRKBL = MAX( WRKBL, 2*M+M*
00482      $                    ILAENV( 1, 'CUNMBR', 'QLN', M, M, M, -1 ) )
00483                   MAXWRK = M*M + WRKBL
00484                   MINWRK = M*M + 3*M
00485                ELSE IF( WNTQA ) THEN
00486 *
00487 *                 Path 4t (N much larger than M, JOBZ='A')
00488 *
00489                   WRKBL = M + M*ILAENV( 1, 'CGELQF', ' ', M, N, -1, -1 )
00490                   WRKBL = MAX( WRKBL, M+N*ILAENV( 1, 'CUNGLQ', ' ', N,
00491      $                    N, M, -1 ) )
00492                   WRKBL = MAX( WRKBL, 2*M+2*M*
00493      $                    ILAENV( 1, 'CGEBRD', ' ', M, M, -1, -1 ) )
00494                   WRKBL = MAX( WRKBL, 2*M+M*
00495      $                    ILAENV( 1, 'CUNMBR', 'PRC', M, M, M, -1 ) )
00496                   WRKBL = MAX( WRKBL, 2*M+M*
00497      $                    ILAENV( 1, 'CUNMBR', 'QLN', M, M, M, -1 ) )
00498                   MAXWRK = M*M + WRKBL
00499                   MINWRK = M*M + 2*M + N
00500                END IF
00501             ELSE IF( N.GE.MNTHR2 ) THEN
00502 *
00503 *              Path 5t (N much larger than M, but not as much as MNTHR1)
00504 *
00505                MAXWRK = 2*M + ( M+N )*ILAENV( 1, 'CGEBRD', ' ', M, N,
00506      $                  -1, -1 )
00507                MINWRK = 2*M + N
00508                IF( WNTQO ) THEN
00509                   MAXWRK = MAX( MAXWRK, 2*M+M*
00510      $                     ILAENV( 1, 'CUNGBR', 'P', M, N, M, -1 ) )
00511                   MAXWRK = MAX( MAXWRK, 2*M+M*
00512      $                     ILAENV( 1, 'CUNGBR', 'Q', M, M, N, -1 ) )
00513                   MAXWRK = MAXWRK + M*N
00514                   MINWRK = MINWRK + M*M
00515                ELSE IF( WNTQS ) THEN
00516                   MAXWRK = MAX( MAXWRK, 2*M+M*
00517      $                     ILAENV( 1, 'CUNGBR', 'P', M, N, M, -1 ) )
00518                   MAXWRK = MAX( MAXWRK, 2*M+M*
00519      $                     ILAENV( 1, 'CUNGBR', 'Q', M, M, N, -1 ) )
00520                ELSE IF( WNTQA ) THEN
00521                   MAXWRK = MAX( MAXWRK, 2*M+N*
00522      $                     ILAENV( 1, 'CUNGBR', 'P', N, N, M, -1 ) )
00523                   MAXWRK = MAX( MAXWRK, 2*M+M*
00524      $                     ILAENV( 1, 'CUNGBR', 'Q', M, M, N, -1 ) )
00525                END IF
00526             ELSE
00527 *
00528 *              Path 6t (N greater than M, but not much larger)
00529 *
00530                MAXWRK = 2*M + ( M+N )*ILAENV( 1, 'CGEBRD', ' ', M, N,
00531      $                  -1, -1 )
00532                MINWRK = 2*M + N
00533                IF( WNTQO ) THEN
00534                   MAXWRK = MAX( MAXWRK, 2*M+M*
00535      $                     ILAENV( 1, 'CUNMBR', 'PRC', M, N, M, -1 ) )
00536                   MAXWRK = MAX( MAXWRK, 2*M+M*
00537      $                     ILAENV( 1, 'CUNMBR', 'QLN', M, M, N, -1 ) )
00538                   MAXWRK = MAXWRK + M*N
00539                   MINWRK = MINWRK + M*M
00540                ELSE IF( WNTQS ) THEN
00541                   MAXWRK = MAX( MAXWRK, 2*M+M*
00542      $                     ILAENV( 1, 'CUNGBR', 'PRC', M, N, M, -1 ) )
00543                   MAXWRK = MAX( MAXWRK, 2*M+M*
00544      $                     ILAENV( 1, 'CUNGBR', 'QLN', M, M, N, -1 ) )
00545                ELSE IF( WNTQA ) THEN
00546                   MAXWRK = MAX( MAXWRK, 2*M+N*
00547      $                     ILAENV( 1, 'CUNGBR', 'PRC', N, N, M, -1 ) )
00548                   MAXWRK = MAX( MAXWRK, 2*M+M*
00549      $                     ILAENV( 1, 'CUNGBR', 'QLN', M, M, N, -1 ) )
00550                END IF
00551             END IF
00552          END IF
00553          MAXWRK = MAX( MAXWRK, MINWRK )
00554       END IF
00555       IF( INFO.EQ.0 ) THEN
00556          WORK( 1 ) = MAXWRK
00557          IF( LWORK.LT.MINWRK .AND. LWORK.NE.LQUERV )
00558      $      INFO = -13
00559       END IF
00560 *
00561 *     Quick returns
00562 *
00563       IF( INFO.NE.0 ) THEN
00564          CALL XERBLA( 'CGESDD', -INFO )
00565          RETURN
00566       END IF
00567       IF( LWORK.EQ.LQUERV )
00568      $   RETURN
00569       IF( M.EQ.0 .OR. N.EQ.0 ) THEN
00570          RETURN
00571       END IF
00572 *
00573 *     Get machine constants
00574 *
00575       EPS = SLAMCH( 'P' )
00576       SMLNUM = SQRT( SLAMCH( 'S' ) ) / EPS
00577       BIGNUM = ONE / SMLNUM
00578 *
00579 *     Scale A if max element outside range [SMLNUM,BIGNUM]
00580 *
00581       ANRM = CLANGE( 'M', M, N, A, LDA, DUM )
00582       ISCL = 0
00583       IF( ANRM.GT.ZERO .AND. ANRM.LT.SMLNUM ) THEN
00584          ISCL = 1
00585          CALL CLASCL( 'G', 0, 0, ANRM, SMLNUM, M, N, A, LDA, IERR )
00586       ELSE IF( ANRM.GT.BIGNUM ) THEN
00587          ISCL = 1
00588          CALL CLASCL( 'G', 0, 0, ANRM, BIGNUM, M, N, A, LDA, IERR )
00589       END IF
00590 *
00591       IF( M.GE.N ) THEN
00592 *
00593 *        A has at least as many rows as columns. If A has sufficiently
00594 *        more rows than columns, first reduce using the QR
00595 *        decomposition (if sufficient workspace available)
00596 *
00597          IF( M.GE.MNTHR1 ) THEN
00598 *
00599             IF( WNTQN ) THEN
00600 *
00601 *              Path 1 (M much larger than N, JOBZ='N')
00602 *              No singular vectors to be computed
00603 *
00604                ITAU = 1
00605                NWORK = ITAU + N
00606 *
00607 *              Compute A=Q*R
00608 *              (CWorkspace: need 2*N, prefer N+N*NB)
00609 *              (RWorkspace: need 0)
00610 *
00611                CALL CGEQRF( M, N, A, LDA, WORK( ITAU ), WORK( NWORK ),
00612      $                      LWORK-NWORK+1, IERR )
00613 *
00614 *              Zero out below R
00615 *
00616                CALL CLASET( 'L', N-1, N-1, CZERO, CZERO, A( 2, 1 ),
00617      $                      LDA )
00618                IE = 1
00619                ITAUQ = 1
00620                ITAUP = ITAUQ + N
00621                NWORK = ITAUP + N
00622 *
00623 *              Bidiagonalize R in A
00624 *              (CWorkspace: need 3*N, prefer 2*N+2*N*NB)
00625 *              (RWorkspace: need N)
00626 *
00627                CALL CGEBRD( N, N, A, LDA, S, RWORK( IE ), WORK( ITAUQ ),
00628      $                      WORK( ITAUP ), WORK( NWORK ), LWORK-NWORK+1,
00629      $                      IERR )
00630                NRWORK = IE + N
00631 *
00632 *              Perform bidiagonal SVD, compute singular values only
00633 *              (CWorkspace: 0)
00634 *              (RWorkspace: need BDSPAN)
00635 *
00636                CALL SBDSDC( 'U', 'N', N, S, RWORK( IE ), DUM, 1, DUM, 1,
00637      $                      DUM, IDUM, RWORK( NRWORK ), IWORK, INFO )
00638 *
00639             ELSE IF( WNTQO ) THEN
00640 *
00641 *              Path 2 (M much larger than N, JOBZ='O')
00642 *              N left singular vectors to be overwritten on A and
00643 *              N right singular vectors to be computed in VT
00644 *
00645                IU = 1
00646 *
00647 *              WORK(IU) is N by N
00648 *
00649                LDWRKU = N
00650                IR = IU + LDWRKU*N
00651                IF( LWORK.GE.M*N+N*N+3*N ) THEN
00652 *
00653 *                 WORK(IR) is M by N
00654 *
00655                   LDWRKR = M
00656                ELSE
00657                   LDWRKR = ( LWORK-N*N-3*N ) / N
00658                END IF
00659                ITAU = IR + LDWRKR*N
00660                NWORK = ITAU + N
00661 *
00662 *              Compute A=Q*R
00663 *              (CWorkspace: need N*N+2*N, prefer M*N+N+N*NB)
00664 *              (RWorkspace: 0)
00665 *
00666                CALL CGEQRF( M, N, A, LDA, WORK( ITAU ), WORK( NWORK ),
00667      $                      LWORK-NWORK+1, IERR )
00668 *
00669 *              Copy R to WORK( IR ), zeroing out below it
00670 *
00671                CALL CLACPY( 'U', N, N, A, LDA, WORK( IR ), LDWRKR )
00672                CALL CLASET( 'L', N-1, N-1, CZERO, CZERO, WORK( IR+1 ),
00673      $                      LDWRKR )
00674 *
00675 *              Generate Q in A
00676 *              (CWorkspace: need 2*N, prefer N+N*NB)
00677 *              (RWorkspace: 0)
00678 *
00679                CALL CUNGQR( M, N, N, A, LDA, WORK( ITAU ),
00680      $                      WORK( NWORK ), LWORK-NWORK+1, IERR )
00681                IE = 1
00682                ITAUQ = ITAU
00683                ITAUP = ITAUQ + N
00684                NWORK = ITAUP + N
00685 *
00686 *              Bidiagonalize R in WORK(IR)
00687 *              (CWorkspace: need N*N+3*N, prefer M*N+2*N+2*N*NB)
00688 *              (RWorkspace: need N)
00689 *
00690                CALL CGEBRD( N, N, WORK( IR ), LDWRKR, S, RWORK( IE ),
00691      $                      WORK( ITAUQ ), WORK( ITAUP ), WORK( NWORK ),
00692      $                      LWORK-NWORK+1, IERR )
00693 *
00694 *              Perform bidiagonal SVD, computing left singular vectors
00695 *              of R in WORK(IRU) and computing right singular vectors
00696 *              of R in WORK(IRVT)
00697 *              (CWorkspace: need 0)
00698 *              (RWorkspace: need BDSPAC)
00699 *
00700                IRU = IE + N
00701                IRVT = IRU + N*N
00702                NRWORK = IRVT + N*N
00703                CALL SBDSDC( 'U', 'I', N, S, RWORK( IE ), RWORK( IRU ),
00704      $                      N, RWORK( IRVT ), N, DUM, IDUM,
00705      $                      RWORK( NRWORK ), IWORK, INFO )
00706 *
00707 *              Copy real matrix RWORK(IRU) to complex matrix WORK(IU)
00708 *              Overwrite WORK(IU) by the left singular vectors of R
00709 *              (CWorkspace: need 2*N*N+3*N, prefer M*N+N*N+2*N+N*NB)
00710 *              (RWorkspace: 0)
00711 *
00712                CALL CLACP2( 'F', N, N, RWORK( IRU ), N, WORK( IU ),
00713      $                      LDWRKU )
00714                CALL CUNMBR( 'Q', 'L', 'N', N, N, N, WORK( IR ), LDWRKR,
00715      $                      WORK( ITAUQ ), WORK( IU ), LDWRKU,
00716      $                      WORK( NWORK ), LWORK-NWORK+1, IERR )
00717 *
00718 *              Copy real matrix RWORK(IRVT) to complex matrix VT
00719 *              Overwrite VT by the right singular vectors of R
00720 *              (CWorkspace: need N*N+3*N, prefer M*N+2*N+N*NB)
00721 *              (RWorkspace: 0)
00722 *
00723                CALL CLACP2( 'F', N, N, RWORK( IRVT ), N, VT, LDVT )
00724                CALL CUNMBR( 'P', 'R', 'C', N, N, N, WORK( IR ), LDWRKR,
00725      $                      WORK( ITAUP ), VT, LDVT, WORK( NWORK ),
00726      $                      LWORK-NWORK+1, IERR )
00727 *
00728 *              Multiply Q in A by left singular vectors of R in
00729 *              WORK(IU), storing result in WORK(IR) and copying to A
00730 *              (CWorkspace: need 2*N*N, prefer N*N+M*N)
00731 *              (RWorkspace: 0)
00732 *
00733                DO 10 I = 1, M, LDWRKR
00734                   CHUNK = MIN( M-I+1, LDWRKR )
00735                   CALL CGEMM( 'N', 'N', CHUNK, N, N, CONE, A( I, 1 ),
00736      $                        LDA, WORK( IU ), LDWRKU, CZERO,
00737      $                        WORK( IR ), LDWRKR )
00738                   CALL CLACPY( 'F', CHUNK, N, WORK( IR ), LDWRKR,
00739      $                         A( I, 1 ), LDA )
00740    10          CONTINUE
00741 *
00742             ELSE IF( WNTQS ) THEN
00743 *
00744 *              Path 3 (M much larger than N, JOBZ='S')
00745 *              N left singular vectors to be computed in U and
00746 *              N right singular vectors to be computed in VT
00747 *
00748                IR = 1
00749 *
00750 *              WORK(IR) is N by N
00751 *
00752                LDWRKR = N
00753                ITAU = IR + LDWRKR*N
00754                NWORK = ITAU + N
00755 *
00756 *              Compute A=Q*R
00757 *              (CWorkspace: need N*N+2*N, prefer N*N+N+N*NB)
00758 *              (RWorkspace: 0)
00759 *
00760                CALL CGEQRF( M, N, A, LDA, WORK( ITAU ), WORK( NWORK ),
00761      $                      LWORK-NWORK+1, IERR )
00762 *
00763 *              Copy R to WORK(IR), zeroing out below it
00764 *
00765                CALL CLACPY( 'U', N, N, A, LDA, WORK( IR ), LDWRKR )
00766                CALL CLASET( 'L', N-1, N-1, CZERO, CZERO, WORK( IR+1 ),
00767      $                      LDWRKR )
00768 *
00769 *              Generate Q in A
00770 *              (CWorkspace: need 2*N, prefer N+N*NB)
00771 *              (RWorkspace: 0)
00772 *
00773                CALL CUNGQR( M, N, N, A, LDA, WORK( ITAU ),
00774      $                      WORK( NWORK ), LWORK-NWORK+1, IERR )
00775                IE = 1
00776                ITAUQ = ITAU
00777                ITAUP = ITAUQ + N
00778                NWORK = ITAUP + N
00779 *
00780 *              Bidiagonalize R in WORK(IR)
00781 *              (CWorkspace: need N*N+3*N, prefer N*N+2*N+2*N*NB)
00782 *              (RWorkspace: need N)
00783 *
00784                CALL CGEBRD( N, N, WORK( IR ), LDWRKR, S, RWORK( IE ),
00785      $                      WORK( ITAUQ ), WORK( ITAUP ), WORK( NWORK ),
00786      $                      LWORK-NWORK+1, IERR )
00787 *
00788 *              Perform bidiagonal SVD, computing left singular vectors
00789 *              of bidiagonal matrix in RWORK(IRU) and computing right
00790 *              singular vectors of bidiagonal matrix in RWORK(IRVT)
00791 *              (CWorkspace: need 0)
00792 *              (RWorkspace: need BDSPAC)
00793 *
00794                IRU = IE + N
00795                IRVT = IRU + N*N
00796                NRWORK = IRVT + N*N
00797                CALL SBDSDC( 'U', 'I', N, S, RWORK( IE ), RWORK( IRU ),
00798      $                      N, RWORK( IRVT ), N, DUM, IDUM,
00799      $                      RWORK( NRWORK ), IWORK, INFO )
00800 *
00801 *              Copy real matrix RWORK(IRU) to complex matrix U
00802 *              Overwrite U by left singular vectors of R
00803 *              (CWorkspace: need N*N+3*N, prefer N*N+2*N+N*NB)
00804 *              (RWorkspace: 0)
00805 *
00806                CALL CLACP2( 'F', N, N, RWORK( IRU ), N, U, LDU )
00807                CALL CUNMBR( 'Q', 'L', 'N', N, N, N, WORK( IR ), LDWRKR,
00808      $                      WORK( ITAUQ ), U, LDU, WORK( NWORK ),
00809      $                      LWORK-NWORK+1, IERR )
00810 *
00811 *              Copy real matrix RWORK(IRVT) to complex matrix VT
00812 *              Overwrite VT by right singular vectors of R
00813 *              (CWorkspace: need N*N+3*N, prefer N*N+2*N+N*NB)
00814 *              (RWorkspace: 0)
00815 *
00816                CALL CLACP2( 'F', N, N, RWORK( IRVT ), N, VT, LDVT )
00817                CALL CUNMBR( 'P', 'R', 'C', N, N, N, WORK( IR ), LDWRKR,
00818      $                      WORK( ITAUP ), VT, LDVT, WORK( NWORK ),
00819      $                      LWORK-NWORK+1, IERR )
00820 *
00821 *              Multiply Q in A by left singular vectors of R in
00822 *              WORK(IR), storing result in U
00823 *              (CWorkspace: need N*N)
00824 *              (RWorkspace: 0)
00825 *
00826                CALL CLACPY( 'F', N, N, U, LDU, WORK( IR ), LDWRKR )
00827                CALL CGEMM( 'N', 'N', M, N, N, CONE, A, LDA, WORK( IR ),
00828      $                     LDWRKR, CZERO, U, LDU )
00829 *
00830             ELSE IF( WNTQA ) THEN
00831 *
00832 *              Path 4 (M much larger than N, JOBZ='A')
00833 *              M left singular vectors to be computed in U and
00834 *              N right singular vectors to be computed in VT
00835 *
00836                IU = 1
00837 *
00838 *              WORK(IU) is N by N
00839 *
00840                LDWRKU = N
00841                ITAU = IU + LDWRKU*N
00842                NWORK = ITAU + N
00843 *
00844 *              Compute A=Q*R, copying result to U
00845 *              (CWorkspace: need 2*N, prefer N+N*NB)
00846 *              (RWorkspace: 0)
00847 *
00848                CALL CGEQRF( M, N, A, LDA, WORK( ITAU ), WORK( NWORK ),
00849      $                      LWORK-NWORK+1, IERR )
00850                CALL CLACPY( 'L', M, N, A, LDA, U, LDU )
00851 *
00852 *              Generate Q in U
00853 *              (CWorkspace: need N+M, prefer N+M*NB)
00854 *              (RWorkspace: 0)
00855 *
00856                CALL CUNGQR( M, M, N, U, LDU, WORK( ITAU ),
00857      $                      WORK( NWORK ), LWORK-NWORK+1, IERR )
00858 *
00859 *              Produce R in A, zeroing out below it
00860 *
00861                CALL CLASET( 'L', N-1, N-1, CZERO, CZERO, A( 2, 1 ),
00862      $                      LDA )
00863                IE = 1
00864                ITAUQ = ITAU
00865                ITAUP = ITAUQ + N
00866                NWORK = ITAUP + N
00867 *
00868 *              Bidiagonalize R in A
00869 *              (CWorkspace: need 3*N, prefer 2*N+2*N*NB)
00870 *              (RWorkspace: need N)
00871 *
00872                CALL CGEBRD( N, N, A, LDA, S, RWORK( IE ), WORK( ITAUQ ),
00873      $                      WORK( ITAUP ), WORK( NWORK ), LWORK-NWORK+1,
00874      $                      IERR )
00875                IRU = IE + N
00876                IRVT = IRU + N*N
00877                NRWORK = IRVT + N*N
00878 *
00879 *              Perform bidiagonal SVD, computing left singular vectors
00880 *              of bidiagonal matrix in RWORK(IRU) and computing right
00881 *              singular vectors of bidiagonal matrix in RWORK(IRVT)
00882 *              (CWorkspace: need 0)
00883 *              (RWorkspace: need BDSPAC)
00884 *
00885                CALL SBDSDC( 'U', 'I', N, S, RWORK( IE ), RWORK( IRU ),
00886      $                      N, RWORK( IRVT ), N, DUM, IDUM,
00887      $                      RWORK( NRWORK ), IWORK, INFO )
00888 *
00889 *              Copy real matrix RWORK(IRU) to complex matrix WORK(IU)
00890 *              Overwrite WORK(IU) by left singular vectors of R
00891 *              (CWorkspace: need N*N+3*N, prefer N*N+2*N+N*NB)
00892 *              (RWorkspace: 0)
00893 *
00894                CALL CLACP2( 'F', N, N, RWORK( IRU ), N, WORK( IU ),
00895      $                      LDWRKU )
00896                CALL CUNMBR( 'Q', 'L', 'N', N, N, N, A, LDA,
00897      $                      WORK( ITAUQ ), WORK( IU ), LDWRKU,
00898      $                      WORK( NWORK ), LWORK-NWORK+1, IERR )
00899 *
00900 *              Copy real matrix RWORK(IRVT) to complex matrix VT
00901 *              Overwrite VT by right singular vectors of R
00902 *              (CWorkspace: need 3*N, prefer 2*N+N*NB)
00903 *              (RWorkspace: 0)
00904 *
00905                CALL CLACP2( 'F', N, N, RWORK( IRVT ), N, VT, LDVT )
00906                CALL CUNMBR( 'P', 'R', 'C', N, N, N, A, LDA,
00907      $                      WORK( ITAUP ), VT, LDVT, WORK( NWORK ),
00908      $                      LWORK-NWORK+1, IERR )
00909 *
00910 *              Multiply Q in U by left singular vectors of R in
00911 *              WORK(IU), storing result in A
00912 *              (CWorkspace: need N*N)
00913 *              (RWorkspace: 0)
00914 *
00915                CALL CGEMM( 'N', 'N', M, N, N, CONE, U, LDU, WORK( IU ),
00916      $                     LDWRKU, CZERO, A, LDA )
00917 *
00918 *              Copy left singular vectors of A from A to U
00919 *
00920                CALL CLACPY( 'F', M, N, A, LDA, U, LDU )
00921 *
00922             END IF
00923 *
00924          ELSE IF( M.GE.MNTHR2 ) THEN
00925 *
00926 *           MNTHR2 <= M < MNTHR1
00927 *
00928 *           Path 5 (M much larger than N, but not as much as MNTHR1)
00929 *           Reduce to bidiagonal form without QR decomposition, use
00930 *           CUNGBR and matrix multiplication to compute singular vectors
00931 *
00932             IE = 1
00933             NRWORK = IE + N
00934             ITAUQ = 1
00935             ITAUP = ITAUQ + N
00936             NWORK = ITAUP + N
00937 *
00938 *           Bidiagonalize A
00939 *           (CWorkspace: need 2*N+M, prefer 2*N+(M+N)*NB)
00940 *           (RWorkspace: need N)
00941 *
00942             CALL CGEBRD( M, N, A, LDA, S, RWORK( IE ), WORK( ITAUQ ),
00943      $                   WORK( ITAUP ), WORK( NWORK ), LWORK-NWORK+1,
00944      $                   IERR )
00945             IF( WNTQN ) THEN
00946 *
00947 *              Compute singular values only
00948 *              (Cworkspace: 0)
00949 *              (Rworkspace: need BDSPAN)
00950 *
00951                CALL SBDSDC( 'U', 'N', N, S, RWORK( IE ), DUM, 1, DUM, 1,
00952      $                      DUM, IDUM, RWORK( NRWORK ), IWORK, INFO )
00953             ELSE IF( WNTQO ) THEN
00954                IU = NWORK
00955                IRU = NRWORK
00956                IRVT = IRU + N*N
00957                NRWORK = IRVT + N*N
00958 *
00959 *              Copy A to VT, generate P**H
00960 *              (Cworkspace: need 2*N, prefer N+N*NB)
00961 *              (Rworkspace: 0)
00962 *
00963                CALL CLACPY( 'U', N, N, A, LDA, VT, LDVT )
00964                CALL CUNGBR( 'P', N, N, N, VT, LDVT, WORK( ITAUP ),
00965      $                      WORK( NWORK ), LWORK-NWORK+1, IERR )
00966 *
00967 *              Generate Q in A
00968 *              (CWorkspace: need 2*N, prefer N+N*NB)
00969 *              (RWorkspace: 0)
00970 *
00971                CALL CUNGBR( 'Q', M, N, N, A, LDA, WORK( ITAUQ ),
00972      $                      WORK( NWORK ), LWORK-NWORK+1, IERR )
00973 *
00974                IF( LWORK.GE.M*N+3*N ) THEN
00975 *
00976 *                 WORK( IU ) is M by N
00977 *
00978                   LDWRKU = M
00979                ELSE
00980 *
00981 *                 WORK(IU) is LDWRKU by N
00982 *
00983                   LDWRKU = ( LWORK-3*N ) / N
00984                END IF
00985                NWORK = IU + LDWRKU*N
00986 *
00987 *              Perform bidiagonal SVD, computing left singular vectors
00988 *              of bidiagonal matrix in RWORK(IRU) and computing right
00989 *              singular vectors of bidiagonal matrix in RWORK(IRVT)
00990 *              (CWorkspace: need 0)
00991 *              (RWorkspace: need BDSPAC)
00992 *
00993                CALL SBDSDC( 'U', 'I', N, S, RWORK( IE ), RWORK( IRU ),
00994      $                      N, RWORK( IRVT ), N, DUM, IDUM,
00995      $                      RWORK( NRWORK ), IWORK, INFO )
00996 *
00997 *              Multiply real matrix RWORK(IRVT) by P**H in VT,
00998 *              storing the result in WORK(IU), copying to VT
00999 *              (Cworkspace: need 0)
01000 *              (Rworkspace: need 3*N*N)
01001 *
01002                CALL CLARCM( N, N, RWORK( IRVT ), N, VT, LDVT,
01003      $                      WORK( IU ), LDWRKU, RWORK( NRWORK ) )
01004                CALL CLACPY( 'F', N, N, WORK( IU ), LDWRKU, VT, LDVT )
01005 *
01006 *              Multiply Q in A by real matrix RWORK(IRU), storing the
01007 *              result in WORK(IU), copying to A
01008 *              (CWorkspace: need N*N, prefer M*N)
01009 *              (Rworkspace: need 3*N*N, prefer N*N+2*M*N)
01010 *
01011                NRWORK = IRVT
01012                DO 20 I = 1, M, LDWRKU
01013                   CHUNK = MIN( M-I+1, LDWRKU )
01014                   CALL CLACRM( CHUNK, N, A( I, 1 ), LDA, RWORK( IRU ),
01015      $                         N, WORK( IU ), LDWRKU, RWORK( NRWORK ) )
01016                   CALL CLACPY( 'F', CHUNK, N, WORK( IU ), LDWRKU,
01017      $                         A( I, 1 ), LDA )
01018    20          CONTINUE
01019 *
01020             ELSE IF( WNTQS ) THEN
01021 *
01022 *              Copy A to VT, generate P**H
01023 *              (Cworkspace: need 2*N, prefer N+N*NB)
01024 *              (Rworkspace: 0)
01025 *
01026                CALL CLACPY( 'U', N, N, A, LDA, VT, LDVT )
01027                CALL CUNGBR( 'P', N, N, N, VT, LDVT, WORK( ITAUP ),
01028      $                      WORK( NWORK ), LWORK-NWORK+1, IERR )
01029 *
01030 *              Copy A to U, generate Q
01031 *              (Cworkspace: need 2*N, prefer N+N*NB)
01032 *              (Rworkspace: 0)
01033 *
01034                CALL CLACPY( 'L', M, N, A, LDA, U, LDU )
01035                CALL CUNGBR( 'Q', M, N, N, U, LDU, WORK( ITAUQ ),
01036      $                      WORK( NWORK ), LWORK-NWORK+1, IERR )
01037 *
01038 *              Perform bidiagonal SVD, computing left singular vectors
01039 *              of bidiagonal matrix in RWORK(IRU) and computing right
01040 *              singular vectors of bidiagonal matrix in RWORK(IRVT)
01041 *              (CWorkspace: need 0)
01042 *              (RWorkspace: need BDSPAC)
01043 *
01044                IRU = NRWORK
01045                IRVT = IRU + N*N
01046                NRWORK = IRVT + N*N
01047                CALL SBDSDC( 'U', 'I', N, S, RWORK( IE ), RWORK( IRU ),
01048      $                      N, RWORK( IRVT ), N, DUM, IDUM,
01049      $                      RWORK( NRWORK ), IWORK, INFO )
01050 *
01051 *              Multiply real matrix RWORK(IRVT) by P**H in VT,
01052 *              storing the result in A, copying to VT
01053 *              (Cworkspace: need 0)
01054 *              (Rworkspace: need 3*N*N)
01055 *
01056                CALL CLARCM( N, N, RWORK( IRVT ), N, VT, LDVT, A, LDA,
01057      $                      RWORK( NRWORK ) )
01058                CALL CLACPY( 'F', N, N, A, LDA, VT, LDVT )
01059 *
01060 *              Multiply Q in U by real matrix RWORK(IRU), storing the
01061 *              result in A, copying to U
01062 *              (CWorkspace: need 0)
01063 *              (Rworkspace: need N*N+2*M*N)
01064 *
01065                NRWORK = IRVT
01066                CALL CLACRM( M, N, U, LDU, RWORK( IRU ), N, A, LDA,
01067      $                      RWORK( NRWORK ) )
01068                CALL CLACPY( 'F', M, N, A, LDA, U, LDU )
01069             ELSE
01070 *
01071 *              Copy A to VT, generate P**H
01072 *              (Cworkspace: need 2*N, prefer N+N*NB)
01073 *              (Rworkspace: 0)
01074 *
01075                CALL CLACPY( 'U', N, N, A, LDA, VT, LDVT )
01076                CALL CUNGBR( 'P', N, N, N, VT, LDVT, WORK( ITAUP ),
01077      $                      WORK( NWORK ), LWORK-NWORK+1, IERR )
01078 *
01079 *              Copy A to U, generate Q
01080 *              (Cworkspace: need 2*N, prefer N+N*NB)
01081 *              (Rworkspace: 0)
01082 *
01083                CALL CLACPY( 'L', M, N, A, LDA, U, LDU )
01084                CALL CUNGBR( 'Q', M, M, N, U, LDU, WORK( ITAUQ ),
01085      $                      WORK( NWORK ), LWORK-NWORK+1, IERR )
01086 *
01087 *              Perform bidiagonal SVD, computing left singular vectors
01088 *              of bidiagonal matrix in RWORK(IRU) and computing right
01089 *              singular vectors of bidiagonal matrix in RWORK(IRVT)
01090 *              (CWorkspace: need 0)
01091 *              (RWorkspace: need BDSPAC)
01092 *
01093                IRU = NRWORK
01094                IRVT = IRU + N*N
01095                NRWORK = IRVT + N*N
01096                CALL SBDSDC( 'U', 'I', N, S, RWORK( IE ), RWORK( IRU ),
01097      $                      N, RWORK( IRVT ), N, DUM, IDUM,
01098      $                      RWORK( NRWORK ), IWORK, INFO )
01099 *
01100 *              Multiply real matrix RWORK(IRVT) by P**H in VT,
01101 *              storing the result in A, copying to VT
01102 *              (Cworkspace: need 0)
01103 *              (Rworkspace: need 3*N*N)
01104 *
01105                CALL CLARCM( N, N, RWORK( IRVT ), N, VT, LDVT, A, LDA,
01106      $                      RWORK( NRWORK ) )
01107                CALL CLACPY( 'F', N, N, A, LDA, VT, LDVT )
01108 *
01109 *              Multiply Q in U by real matrix RWORK(IRU), storing the
01110 *              result in A, copying to U
01111 *              (CWorkspace: 0)
01112 *              (Rworkspace: need 3*N*N)
01113 *
01114                NRWORK = IRVT
01115                CALL CLACRM( M, N, U, LDU, RWORK( IRU ), N, A, LDA,
01116      $                      RWORK( NRWORK ) )
01117                CALL CLACPY( 'F', M, N, A, LDA, U, LDU )
01118             END IF
01119 *
01120          ELSE
01121 *
01122 *           M .LT. MNTHR2
01123 *
01124 *           Path 6 (M at least N, but not much larger)
01125 *           Reduce to bidiagonal form without QR decomposition
01126 *           Use CUNMBR to compute singular vectors
01127 *
01128             IE = 1
01129             NRWORK = IE + N
01130             ITAUQ = 1
01131             ITAUP = ITAUQ + N
01132             NWORK = ITAUP + N
01133 *
01134 *           Bidiagonalize A
01135 *           (CWorkspace: need 2*N+M, prefer 2*N+(M+N)*NB)
01136 *           (RWorkspace: need N)
01137 *
01138             CALL CGEBRD( M, N, A, LDA, S, RWORK( IE ), WORK( ITAUQ ),
01139      $                   WORK( ITAUP ), WORK( NWORK ), LWORK-NWORK+1,
01140      $                   IERR )
01141             IF( WNTQN ) THEN
01142 *
01143 *              Compute singular values only
01144 *              (Cworkspace: 0)
01145 *              (Rworkspace: need BDSPAN)
01146 *
01147                CALL SBDSDC( 'U', 'N', N, S, RWORK( IE ), DUM, 1, DUM, 1,
01148      $                      DUM, IDUM, RWORK( NRWORK ), IWORK, INFO )
01149             ELSE IF( WNTQO ) THEN
01150                IU = NWORK
01151                IRU = NRWORK
01152                IRVT = IRU + N*N
01153                NRWORK = IRVT + N*N
01154                IF( LWORK.GE.M*N+3*N ) THEN
01155 *
01156 *                 WORK( IU ) is M by N
01157 *
01158                   LDWRKU = M
01159                ELSE
01160 *
01161 *                 WORK( IU ) is LDWRKU by N
01162 *
01163                   LDWRKU = ( LWORK-3*N ) / N
01164                END IF
01165                NWORK = IU + LDWRKU*N
01166 *
01167 *              Perform bidiagonal SVD, computing left singular vectors
01168 *              of bidiagonal matrix in RWORK(IRU) and computing right
01169 *              singular vectors of bidiagonal matrix in RWORK(IRVT)
01170 *              (CWorkspace: need 0)
01171 *              (RWorkspace: need BDSPAC)
01172 *
01173                CALL SBDSDC( 'U', 'I', N, S, RWORK( IE ), RWORK( IRU ),
01174      $                      N, RWORK( IRVT ), N, DUM, IDUM,
01175      $                      RWORK( NRWORK ), IWORK, INFO )
01176 *
01177 *              Copy real matrix RWORK(IRVT) to complex matrix VT
01178 *              Overwrite VT by right singular vectors of A
01179 *              (Cworkspace: need 2*N, prefer N+N*NB)
01180 *              (Rworkspace: need 0)
01181 *
01182                CALL CLACP2( 'F', N, N, RWORK( IRVT ), N, VT, LDVT )
01183                CALL CUNMBR( 'P', 'R', 'C', N, N, N, A, LDA,
01184      $                      WORK( ITAUP ), VT, LDVT, WORK( NWORK ),
01185      $                      LWORK-NWORK+1, IERR )
01186 *
01187                IF( LWORK.GE.M*N+3*N ) THEN
01188 *
01189 *              Copy real matrix RWORK(IRU) to complex matrix WORK(IU)
01190 *              Overwrite WORK(IU) by left singular vectors of A, copying
01191 *              to A
01192 *              (Cworkspace: need M*N+2*N, prefer M*N+N+N*NB)
01193 *              (Rworkspace: need 0)
01194 *
01195                   CALL CLASET( 'F', M, N, CZERO, CZERO, WORK( IU ),
01196      $                         LDWRKU )
01197                   CALL CLACP2( 'F', N, N, RWORK( IRU ), N, WORK( IU ),
01198      $                         LDWRKU )
01199                   CALL CUNMBR( 'Q', 'L', 'N', M, N, N, A, LDA,
01200      $                         WORK( ITAUQ ), WORK( IU ), LDWRKU,
01201      $                         WORK( NWORK ), LWORK-NWORK+1, IERR )
01202                   CALL CLACPY( 'F', M, N, WORK( IU ), LDWRKU, A, LDA )
01203                ELSE
01204 *
01205 *                 Generate Q in A
01206 *                 (Cworkspace: need 2*N, prefer N+N*NB)
01207 *                 (Rworkspace: need 0)
01208 *
01209                   CALL CUNGBR( 'Q', M, N, N, A, LDA, WORK( ITAUQ ),
01210      $                         WORK( NWORK ), LWORK-NWORK+1, IERR )
01211 *
01212 *                 Multiply Q in A by real matrix RWORK(IRU), storing the
01213 *                 result in WORK(IU), copying to A
01214 *                 (CWorkspace: need N*N, prefer M*N)
01215 *                 (Rworkspace: need 3*N*N, prefer N*N+2*M*N)
01216 *
01217                   NRWORK = IRVT
01218                   DO 30 I = 1, M, LDWRKU
01219                      CHUNK = MIN( M-I+1, LDWRKU )
01220                      CALL CLACRM( CHUNK, N, A( I, 1 ), LDA,
01221      $                            RWORK( IRU ), N, WORK( IU ), LDWRKU,
01222      $                            RWORK( NRWORK ) )
01223                      CALL CLACPY( 'F', CHUNK, N, WORK( IU ), LDWRKU,
01224      $                            A( I, 1 ), LDA )
01225    30             CONTINUE
01226                END IF
01227 *
01228             ELSE IF( WNTQS ) THEN
01229 *
01230 *              Perform bidiagonal SVD, computing left singular vectors
01231 *              of bidiagonal matrix in RWORK(IRU) and computing right
01232 *              singular vectors of bidiagonal matrix in RWORK(IRVT)
01233 *              (CWorkspace: need 0)
01234 *              (RWorkspace: need BDSPAC)
01235 *
01236                IRU = NRWORK
01237                IRVT = IRU + N*N
01238                NRWORK = IRVT + N*N
01239                CALL SBDSDC( 'U', 'I', N, S, RWORK( IE ), RWORK( IRU ),
01240      $                      N, RWORK( IRVT ), N, DUM, IDUM,
01241      $                      RWORK( NRWORK ), IWORK, INFO )
01242 *
01243 *              Copy real matrix RWORK(IRU) to complex matrix U
01244 *              Overwrite U by left singular vectors of A
01245 *              (CWorkspace: need 3*N, prefer 2*N+N*NB)
01246 *              (RWorkspace: 0)
01247 *
01248                CALL CLASET( 'F', M, N, CZERO, CZERO, U, LDU )
01249                CALL CLACP2( 'F', N, N, RWORK( IRU ), N, U, LDU )
01250                CALL CUNMBR( 'Q', 'L', 'N', M, N, N, A, LDA,
01251      $                      WORK( ITAUQ ), U, LDU, WORK( NWORK ),
01252      $                      LWORK-NWORK+1, IERR )
01253 *
01254 *              Copy real matrix RWORK(IRVT) to complex matrix VT
01255 *              Overwrite VT by right singular vectors of A
01256 *              (CWorkspace: need 3*N, prefer 2*N+N*NB)
01257 *              (RWorkspace: 0)
01258 *
01259                CALL CLACP2( 'F', N, N, RWORK( IRVT ), N, VT, LDVT )
01260                CALL CUNMBR( 'P', 'R', 'C', N, N, N, A, LDA,
01261      $                      WORK( ITAUP ), VT, LDVT, WORK( NWORK ),
01262      $                      LWORK-NWORK+1, IERR )
01263             ELSE
01264 *
01265 *              Perform bidiagonal SVD, computing left singular vectors
01266 *              of bidiagonal matrix in RWORK(IRU) and computing right
01267 *              singular vectors of bidiagonal matrix in RWORK(IRVT)
01268 *              (CWorkspace: need 0)
01269 *              (RWorkspace: need BDSPAC)
01270 *
01271                IRU = NRWORK
01272                IRVT = IRU + N*N
01273                NRWORK = IRVT + N*N
01274                CALL SBDSDC( 'U', 'I', N, S, RWORK( IE ), RWORK( IRU ),
01275      $                      N, RWORK( IRVT ), N, DUM, IDUM,
01276      $                      RWORK( NRWORK ), IWORK, INFO )
01277 *
01278 *              Set the right corner of U to identity matrix
01279 *
01280                CALL CLASET( 'F', M, M, CZERO, CZERO, U, LDU )
01281                IF( M.GT.N ) THEN
01282                   CALL CLASET( 'F', M-N, M-N, CZERO, CONE,
01283      $                         U( N+1, N+1 ), LDU )
01284                END IF
01285 *
01286 *              Copy real matrix RWORK(IRU) to complex matrix U
01287 *              Overwrite U by left singular vectors of A
01288 *              (CWorkspace: need 2*N+M, prefer 2*N+M*NB)
01289 *              (RWorkspace: 0)
01290 *
01291                CALL CLACP2( 'F', N, N, RWORK( IRU ), N, U, LDU )
01292                CALL CUNMBR( 'Q', 'L', 'N', M, M, N, A, LDA,
01293      $                      WORK( ITAUQ ), U, LDU, WORK( NWORK ),
01294      $                      LWORK-NWORK+1, IERR )
01295 *
01296 *              Copy real matrix RWORK(IRVT) to complex matrix VT
01297 *              Overwrite VT by right singular vectors of A
01298 *              (CWorkspace: need 3*N, prefer 2*N+N*NB)
01299 *              (RWorkspace: 0)
01300 *
01301                CALL CLACP2( 'F', N, N, RWORK( IRVT ), N, VT, LDVT )
01302                CALL CUNMBR( 'P', 'R', 'C', N, N, N, A, LDA,
01303      $                      WORK( ITAUP ), VT, LDVT, WORK( NWORK ),
01304      $                      LWORK-NWORK+1, IERR )
01305             END IF
01306 *
01307          END IF
01308 *
01309       ELSE
01310 *
01311 *        A has more columns than rows. If A has sufficiently more
01312 *        columns than rows, first reduce using the LQ decomposition (if
01313 *        sufficient workspace available)
01314 *
01315          IF( N.GE.MNTHR1 ) THEN
01316 *
01317             IF( WNTQN ) THEN
01318 *
01319 *              Path 1t (N much larger than M, JOBZ='N')
01320 *              No singular vectors to be computed
01321 *
01322                ITAU = 1
01323                NWORK = ITAU + M
01324 *
01325 *              Compute A=L*Q
01326 *              (CWorkspace: need 2*M, prefer M+M*NB)
01327 *              (RWorkspace: 0)
01328 *
01329                CALL CGELQF( M, N, A, LDA, WORK( ITAU ), WORK( NWORK ),
01330      $                      LWORK-NWORK+1, IERR )
01331 *
01332 *              Zero out above L
01333 *
01334                CALL CLASET( 'U', M-1, M-1, CZERO, CZERO, A( 1, 2 ),
01335      $                      LDA )
01336                IE = 1
01337                ITAUQ = 1
01338                ITAUP = ITAUQ + M
01339                NWORK = ITAUP + M
01340 *
01341 *              Bidiagonalize L in A
01342 *              (CWorkspace: need 3*M, prefer 2*M+2*M*NB)
01343 *              (RWorkspace: need M)
01344 *
01345                CALL CGEBRD( M, M, A, LDA, S, RWORK( IE ), WORK( ITAUQ ),
01346      $                      WORK( ITAUP ), WORK( NWORK ), LWORK-NWORK+1,
01347      $                      IERR )
01348                NRWORK = IE + M
01349 *
01350 *              Perform bidiagonal SVD, compute singular values only
01351 *              (CWorkspace: 0)
01352 *              (RWorkspace: need BDSPAN)
01353 *
01354                CALL SBDSDC( 'U', 'N', M, S, RWORK( IE ), DUM, 1, DUM, 1,
01355      $                      DUM, IDUM, RWORK( NRWORK ), IWORK, INFO )
01356 *
01357             ELSE IF( WNTQO ) THEN
01358 *
01359 *              Path 2t (N much larger than M, JOBZ='O')
01360 *              M right singular vectors to be overwritten on A and
01361 *              M left singular vectors to be computed in U
01362 *
01363                IVT = 1
01364                LDWKVT = M
01365 *
01366 *              WORK(IVT) is M by M
01367 *
01368                IL = IVT + LDWKVT*M
01369                IF( LWORK.GE.M*N+M*M+3*M ) THEN
01370 *
01371 *                 WORK(IL) M by N
01372 *
01373                   LDWRKL = M
01374                   CHUNK = N
01375                ELSE
01376 *
01377 *                 WORK(IL) is M by CHUNK
01378 *
01379                   LDWRKL = M
01380                   CHUNK = ( LWORK-M*M-3*M ) / M
01381                END IF
01382                ITAU = IL + LDWRKL*CHUNK
01383                NWORK = ITAU + M
01384 *
01385 *              Compute A=L*Q
01386 *              (CWorkspace: need 2*M, prefer M+M*NB)
01387 *              (RWorkspace: 0)
01388 *
01389                CALL CGELQF( M, N, A, LDA, WORK( ITAU ), WORK( NWORK ),
01390      $                      LWORK-NWORK+1, IERR )
01391 *
01392 *              Copy L to WORK(IL), zeroing about above it
01393 *
01394                CALL CLACPY( 'L', M, M, A, LDA, WORK( IL ), LDWRKL )
01395                CALL CLASET( 'U', M-1, M-1, CZERO, CZERO,
01396      $                      WORK( IL+LDWRKL ), LDWRKL )
01397 *
01398 *              Generate Q in A
01399 *              (CWorkspace: need M*M+2*M, prefer M*M+M+M*NB)
01400 *              (RWorkspace: 0)
01401 *
01402                CALL CUNGLQ( M, N, M, A, LDA, WORK( ITAU ),
01403      $                      WORK( NWORK ), LWORK-NWORK+1, IERR )
01404                IE = 1
01405                ITAUQ = ITAU
01406                ITAUP = ITAUQ + M
01407                NWORK = ITAUP + M
01408 *
01409 *              Bidiagonalize L in WORK(IL)
01410 *              (CWorkspace: need M*M+3*M, prefer M*M+2*M+2*M*NB)
01411 *              (RWorkspace: need M)
01412 *
01413                CALL CGEBRD( M, M, WORK( IL ), LDWRKL, S, RWORK( IE ),
01414      $                      WORK( ITAUQ ), WORK( ITAUP ), WORK( NWORK ),
01415      $                      LWORK-NWORK+1, IERR )
01416 *
01417 *              Perform bidiagonal SVD, computing left singular vectors
01418 *              of bidiagonal matrix in RWORK(IRU) and computing right
01419 *              singular vectors of bidiagonal matrix in RWORK(IRVT)
01420 *              (CWorkspace: need 0)
01421 *              (RWorkspace: need BDSPAC)
01422 *
01423                IRU = IE + M
01424                IRVT = IRU + M*M
01425                NRWORK = IRVT + M*M
01426                CALL SBDSDC( 'U', 'I', M, S, RWORK( IE ), RWORK( IRU ),
01427      $                      M, RWORK( IRVT ), M, DUM, IDUM,
01428      $                      RWORK( NRWORK ), IWORK, INFO )
01429 *
01430 *              Copy real matrix RWORK(IRU) to complex matrix WORK(IU)
01431 *              Overwrite WORK(IU) by the left singular vectors of L
01432 *              (CWorkspace: need N*N+3*N, prefer M*N+2*N+N*NB)
01433 *              (RWorkspace: 0)
01434 *
01435                CALL CLACP2( 'F', M, M, RWORK( IRU ), M, U, LDU )
01436                CALL CUNMBR( 'Q', 'L', 'N', M, M, M, WORK( IL ), LDWRKL,
01437      $                      WORK( ITAUQ ), U, LDU, WORK( NWORK ),
01438      $                      LWORK-NWORK+1, IERR )
01439 *
01440 *              Copy real matrix RWORK(IRVT) to complex matrix WORK(IVT)
01441 *              Overwrite WORK(IVT) by the right singular vectors of L
01442 *              (CWorkspace: need N*N+3*N, prefer M*N+2*N+N*NB)
01443 *              (RWorkspace: 0)
01444 *
01445                CALL CLACP2( 'F', M, M, RWORK( IRVT ), M, WORK( IVT ),
01446      $                      LDWKVT )
01447                CALL CUNMBR( 'P', 'R', 'C', M, M, M, WORK( IL ), LDWRKL,
01448      $                      WORK( ITAUP ), WORK( IVT ), LDWKVT,
01449      $                      WORK( NWORK ), LWORK-NWORK+1, IERR )
01450 *
01451 *              Multiply right singular vectors of L in WORK(IL) by Q
01452 *              in A, storing result in WORK(IL) and copying to A
01453 *              (CWorkspace: need 2*M*M, prefer M*M+M*N))
01454 *              (RWorkspace: 0)
01455 *
01456                DO 40 I = 1, N, CHUNK
01457                   BLK = MIN( N-I+1, CHUNK )
01458                   CALL CGEMM( 'N', 'N', M, BLK, M, CONE, WORK( IVT ), M,
01459      $                        A( 1, I ), LDA, CZERO, WORK( IL ),
01460      $                        LDWRKL )
01461                   CALL CLACPY( 'F', M, BLK, WORK( IL ), LDWRKL,
01462      $                         A( 1, I ), LDA )
01463    40          CONTINUE
01464 *
01465             ELSE IF( WNTQS ) THEN
01466 *
01467 *             Path 3t (N much larger than M, JOBZ='S')
01468 *             M right singular vectors to be computed in VT and
01469 *             M left singular vectors to be computed in U
01470 *
01471                IL = 1
01472 *
01473 *              WORK(IL) is M by M
01474 *
01475                LDWRKL = M
01476                ITAU = IL + LDWRKL*M
01477                NWORK = ITAU + M
01478 *
01479 *              Compute A=L*Q
01480 *              (CWorkspace: need 2*M, prefer M+M*NB)
01481 *              (RWorkspace: 0)
01482 *
01483                CALL CGELQF( M, N, A, LDA, WORK( ITAU ), WORK( NWORK ),
01484      $                      LWORK-NWORK+1, IERR )
01485 *
01486 *              Copy L to WORK(IL), zeroing out above it
01487 *
01488                CALL CLACPY( 'L', M, M, A, LDA, WORK( IL ), LDWRKL )
01489                CALL CLASET( 'U', M-1, M-1, CZERO, CZERO,
01490      $                      WORK( IL+LDWRKL ), LDWRKL )
01491 *
01492 *              Generate Q in A
01493 *              (CWorkspace: need M*M+2*M, prefer M*M+M+M*NB)
01494 *              (RWorkspace: 0)
01495 *
01496                CALL CUNGLQ( M, N, M, A, LDA, WORK( ITAU ),
01497      $                      WORK( NWORK ), LWORK-NWORK+1, IERR )
01498                IE = 1
01499                ITAUQ = ITAU
01500                ITAUP = ITAUQ + M
01501                NWORK = ITAUP + M
01502 *
01503 *              Bidiagonalize L in WORK(IL)
01504 *              (CWorkspace: need M*M+3*M, prefer M*M+2*M+2*M*NB)
01505 *              (RWorkspace: need M)
01506 *
01507                CALL CGEBRD( M, M, WORK( IL ), LDWRKL, S, RWORK( IE ),
01508      $                      WORK( ITAUQ ), WORK( ITAUP ), WORK( NWORK ),
01509      $                      LWORK-NWORK+1, IERR )
01510 *
01511 *              Perform bidiagonal SVD, computing left singular vectors
01512 *              of bidiagonal matrix in RWORK(IRU) and computing right
01513 *              singular vectors of bidiagonal matrix in RWORK(IRVT)
01514 *              (CWorkspace: need 0)
01515 *              (RWorkspace: need BDSPAC)
01516 *
01517                IRU = IE + M
01518                IRVT = IRU + M*M
01519                NRWORK = IRVT + M*M
01520                CALL SBDSDC( 'U', 'I', M, S, RWORK( IE ), RWORK( IRU ),
01521      $                      M, RWORK( IRVT ), M, DUM, IDUM,
01522      $                      RWORK( NRWORK ), IWORK, INFO )
01523 *
01524 *              Copy real matrix RWORK(IRU) to complex matrix U
01525 *              Overwrite U by left singular vectors of L
01526 *              (CWorkspace: need M*M+3*M, prefer M*M+2*M+M*NB)
01527 *              (RWorkspace: 0)
01528 *
01529                CALL CLACP2( 'F', M, M, RWORK( IRU ), M, U, LDU )
01530                CALL CUNMBR( 'Q', 'L', 'N', M, M, M, WORK( IL ), LDWRKL,
01531      $                      WORK( ITAUQ ), U, LDU, WORK( NWORK ),
01532      $                      LWORK-NWORK+1, IERR )
01533 *
01534 *              Copy real matrix RWORK(IRVT) to complex matrix VT
01535 *              Overwrite VT by left singular vectors of L
01536 *              (CWorkspace: need M*M+3*M, prefer M*M+2*M+M*NB)
01537 *              (RWorkspace: 0)
01538 *
01539                CALL CLACP2( 'F', M, M, RWORK( IRVT ), M, VT, LDVT )
01540                CALL CUNMBR( 'P', 'R', 'C', M, M, M, WORK( IL ), LDWRKL,
01541      $                      WORK( ITAUP ), VT, LDVT, WORK( NWORK ),
01542      $                      LWORK-NWORK+1, IERR )
01543 *
01544 *              Copy VT to WORK(IL), multiply right singular vectors of L
01545 *              in WORK(IL) by Q in A, storing result in VT
01546 *              (CWorkspace: need M*M)
01547 *              (RWorkspace: 0)
01548 *
01549                CALL CLACPY( 'F', M, M, VT, LDVT, WORK( IL ), LDWRKL )
01550                CALL CGEMM( 'N', 'N', M, N, M, CONE, WORK( IL ), LDWRKL,
01551      $                     A, LDA, CZERO, VT, LDVT )
01552 *
01553             ELSE IF( WNTQA ) THEN
01554 *
01555 *              Path 9t (N much larger than M, JOBZ='A')
01556 *              N right singular vectors to be computed in VT and
01557 *              M left singular vectors to be computed in U
01558 *
01559                IVT = 1
01560 *
01561 *              WORK(IVT) is M by M
01562 *
01563                LDWKVT = M
01564                ITAU = IVT + LDWKVT*M
01565                NWORK = ITAU + M
01566 *
01567 *              Compute A=L*Q, copying result to VT
01568 *              (CWorkspace: need 2*M, prefer M+M*NB)
01569 *              (RWorkspace: 0)
01570 *
01571                CALL CGELQF( M, N, A, LDA, WORK( ITAU ), WORK( NWORK ),
01572      $                      LWORK-NWORK+1, IERR )
01573                CALL CLACPY( 'U', M, N, A, LDA, VT, LDVT )
01574 *
01575 *              Generate Q in VT
01576 *              (CWorkspace: need M+N, prefer M+N*NB)
01577 *              (RWorkspace: 0)
01578 *
01579                CALL CUNGLQ( N, N, M, VT, LDVT, WORK( ITAU ),
01580      $                      WORK( NWORK ), LWORK-NWORK+1, IERR )
01581 *
01582 *              Produce L in A, zeroing out above it
01583 *
01584                CALL CLASET( 'U', M-1, M-1, CZERO, CZERO, A( 1, 2 ),
01585      $                      LDA )
01586                IE = 1
01587                ITAUQ = ITAU
01588                ITAUP = ITAUQ + M
01589                NWORK = ITAUP + M
01590 *
01591 *              Bidiagonalize L in A
01592 *              (CWorkspace: need M*M+3*M, prefer M*M+2*M+2*M*NB)
01593 *              (RWorkspace: need M)
01594 *
01595                CALL CGEBRD( M, M, A, LDA, S, RWORK( IE ), WORK( ITAUQ ),
01596      $                      WORK( ITAUP ), WORK( NWORK ), LWORK-NWORK+1,
01597      $                      IERR )
01598 *
01599 *              Perform bidiagonal SVD, computing left singular vectors
01600 *              of bidiagonal matrix in RWORK(IRU) and computing right
01601 *              singular vectors of bidiagonal matrix in RWORK(IRVT)
01602 *              (CWorkspace: need 0)
01603 *              (RWorkspace: need BDSPAC)
01604 *
01605                IRU = IE + M
01606                IRVT = IRU + M*M
01607                NRWORK = IRVT + M*M
01608                CALL SBDSDC( 'U', 'I', M, S, RWORK( IE ), RWORK( IRU ),
01609      $                      M, RWORK( IRVT ), M, DUM, IDUM,
01610      $                      RWORK( NRWORK ), IWORK, INFO )
01611 *
01612 *              Copy real matrix RWORK(IRU) to complex matrix U
01613 *              Overwrite U by left singular vectors of L
01614 *              (CWorkspace: need 3*M, prefer 2*M+M*NB)
01615 *              (RWorkspace: 0)
01616 *
01617                CALL CLACP2( 'F', M, M, RWORK( IRU ), M, U, LDU )
01618                CALL CUNMBR( 'Q', 'L', 'N', M, M, M, A, LDA,
01619      $                      WORK( ITAUQ ), U, LDU, WORK( NWORK ),
01620      $                      LWORK-NWORK+1, IERR )
01621 *
01622 *              Copy real matrix RWORK(IRVT) to complex matrix WORK(IVT)
01623 *              Overwrite WORK(IVT) by right singular vectors of L
01624 *              (CWorkspace: need M*M+3*M, prefer M*M+2*M+M*NB)
01625 *              (RWorkspace: 0)
01626 *
01627                CALL CLACP2( 'F', M, M, RWORK( IRVT ), M, WORK( IVT ),
01628      $                      LDWKVT )
01629                CALL CUNMBR( 'P', 'R', 'C', M, M, M, A, LDA,
01630      $                      WORK( ITAUP ), WORK( IVT ), LDWKVT,
01631      $                      WORK( NWORK ), LWORK-NWORK+1, IERR )
01632 *
01633 *              Multiply right singular vectors of L in WORK(IVT) by
01634 *              Q in VT, storing result in A
01635 *              (CWorkspace: need M*M)
01636 *              (RWorkspace: 0)
01637 *
01638                CALL CGEMM( 'N', 'N', M, N, M, CONE, WORK( IVT ),
01639      $                     LDWKVT, VT, LDVT, CZERO, A, LDA )
01640 *
01641 *              Copy right singular vectors of A from A to VT
01642 *
01643                CALL CLACPY( 'F', M, N, A, LDA, VT, LDVT )
01644 *
01645             END IF
01646 *
01647          ELSE IF( N.GE.MNTHR2 ) THEN
01648 *
01649 *           MNTHR2 <= N < MNTHR1
01650 *
01651 *           Path 5t (N much larger than M, but not as much as MNTHR1)
01652 *           Reduce to bidiagonal form without QR decomposition, use
01653 *           CUNGBR and matrix multiplication to compute singular vectors
01654 *
01655 *
01656             IE = 1
01657             NRWORK = IE + M
01658             ITAUQ = 1
01659             ITAUP = ITAUQ + M
01660             NWORK = ITAUP + M
01661 *
01662 *           Bidiagonalize A
01663 *           (CWorkspace: need 2*M+N, prefer 2*M+(M+N)*NB)
01664 *           (RWorkspace: M)
01665 *
01666             CALL CGEBRD( M, N, A, LDA, S, RWORK( IE ), WORK( ITAUQ ),
01667      $                   WORK( ITAUP ), WORK( NWORK ), LWORK-NWORK+1,
01668      $                   IERR )
01669 *
01670             IF( WNTQN ) THEN
01671 *
01672 *              Compute singular values only
01673 *              (Cworkspace: 0)
01674 *              (Rworkspace: need BDSPAN)
01675 *
01676                CALL SBDSDC( 'L', 'N', M, S, RWORK( IE ), DUM, 1, DUM, 1,
01677      $                      DUM, IDUM, RWORK( NRWORK ), IWORK, INFO )
01678             ELSE IF( WNTQO ) THEN
01679                IRVT = NRWORK
01680                IRU = IRVT + M*M
01681                NRWORK = IRU + M*M
01682                IVT = NWORK
01683 *
01684 *              Copy A to U, generate Q
01685 *              (Cworkspace: need 2*M, prefer M+M*NB)
01686 *              (Rworkspace: 0)
01687 *
01688                CALL CLACPY( 'L', M, M, A, LDA, U, LDU )
01689                CALL CUNGBR( 'Q', M, M, N, U, LDU, WORK( ITAUQ ),
01690      $                      WORK( NWORK ), LWORK-NWORK+1, IERR )
01691 *
01692 *              Generate P**H in A
01693 *              (Cworkspace: need 2*M, prefer M+M*NB)
01694 *              (Rworkspace: 0)
01695 *
01696                CALL CUNGBR( 'P', M, N, M, A, LDA, WORK( ITAUP ),
01697      $                      WORK( NWORK ), LWORK-NWORK+1, IERR )
01698 *
01699                LDWKVT = M
01700                IF( LWORK.GE.M*N+3*M ) THEN
01701 *
01702 *                 WORK( IVT ) is M by N
01703 *
01704                   NWORK = IVT + LDWKVT*N
01705                   CHUNK = N
01706                ELSE
01707 *
01708 *                 WORK( IVT ) is M by CHUNK
01709 *
01710                   CHUNK = ( LWORK-3*M ) / M
01711                   NWORK = IVT + LDWKVT*CHUNK
01712                END IF
01713 *
01714 *              Perform bidiagonal SVD, computing left singular vectors
01715 *              of bidiagonal matrix in RWORK(IRU) and computing right
01716 *              singular vectors of bidiagonal matrix in RWORK(IRVT)
01717 *              (CWorkspace: need 0)
01718 *              (RWorkspace: need BDSPAC)
01719 *
01720                CALL SBDSDC( 'L', 'I', M, S, RWORK( IE ), RWORK( IRU ),
01721      $                      M, RWORK( IRVT ), M, DUM, IDUM,
01722      $                      RWORK( NRWORK ), IWORK, INFO )
01723 *
01724 *              Multiply Q in U by real matrix RWORK(IRVT)
01725 *              storing the result in WORK(IVT), copying to U
01726 *              (Cworkspace: need 0)
01727 *              (Rworkspace: need 2*M*M)
01728 *
01729                CALL CLACRM( M, M, U, LDU, RWORK( IRU ), M, WORK( IVT ),
01730      $                      LDWKVT, RWORK( NRWORK ) )
01731                CALL CLACPY( 'F', M, M, WORK( IVT ), LDWKVT, U, LDU )
01732 *
01733 *              Multiply RWORK(IRVT) by P**H in A, storing the
01734 *              result in WORK(IVT), copying to A
01735 *              (CWorkspace: need M*M, prefer M*N)
01736 *              (Rworkspace: need 2*M*M, prefer 2*M*N)
01737 *
01738                NRWORK = IRU
01739                DO 50 I = 1, N, CHUNK
01740                   BLK = MIN( N-I+1, CHUNK )
01741                   CALL CLARCM( M, BLK, RWORK( IRVT ), M, A( 1, I ), LDA,
01742      $                         WORK( IVT ), LDWKVT, RWORK( NRWORK ) )
01743                   CALL CLACPY( 'F', M, BLK, WORK( IVT ), LDWKVT,
01744      $                         A( 1, I ), LDA )
01745    50          CONTINUE
01746             ELSE IF( WNTQS ) THEN
01747 *
01748 *              Copy A to U, generate Q
01749 *              (Cworkspace: need 2*M, prefer M+M*NB)
01750 *              (Rworkspace: 0)
01751 *
01752                CALL CLACPY( 'L', M, M, A, LDA, U, LDU )
01753                CALL CUNGBR( 'Q', M, M, N, U, LDU, WORK( ITAUQ ),
01754      $                      WORK( NWORK ), LWORK-NWORK+1, IERR )
01755 *
01756 *              Copy A to VT, generate P**H
01757 *              (Cworkspace: need 2*M, prefer M+M*NB)
01758 *              (Rworkspace: 0)
01759 *
01760                CALL CLACPY( 'U', M, N, A, LDA, VT, LDVT )
01761                CALL CUNGBR( 'P', M, N, M, VT, LDVT, WORK( ITAUP ),
01762      $                      WORK( NWORK ), LWORK-NWORK+1, IERR )
01763 *
01764 *              Perform bidiagonal SVD, computing left singular vectors
01765 *              of bidiagonal matrix in RWORK(IRU) and computing right
01766 *              singular vectors of bidiagonal matrix in RWORK(IRVT)
01767 *              (CWorkspace: need 0)
01768 *              (RWorkspace: need BDSPAC)
01769 *
01770                IRVT = NRWORK
01771                IRU = IRVT + M*M
01772                NRWORK = IRU + M*M
01773                CALL SBDSDC( 'L', 'I', M, S, RWORK( IE ), RWORK( IRU ),
01774      $                      M, RWORK( IRVT ), M, DUM, IDUM,
01775      $                      RWORK( NRWORK ), IWORK, INFO )
01776 *
01777 *              Multiply Q in U by real matrix RWORK(IRU), storing the
01778 *              result in A, copying to U
01779 *              (CWorkspace: need 0)
01780 *              (Rworkspace: need 3*M*M)
01781 *
01782                CALL CLACRM( M, M, U, LDU, RWORK( IRU ), M, A, LDA,
01783      $                      RWORK( NRWORK ) )
01784                CALL CLACPY( 'F', M, M, A, LDA, U, LDU )
01785 *
01786 *              Multiply real matrix RWORK(IRVT) by P**H in VT,
01787 *              storing the result in A, copying to VT
01788 *              (Cworkspace: need 0)
01789 *              (Rworkspace: need M*M+2*M*N)
01790 *
01791                NRWORK = IRU
01792                CALL CLARCM( M, N, RWORK( IRVT ), M, VT, LDVT, A, LDA,
01793      $                      RWORK( NRWORK ) )
01794                CALL CLACPY( 'F', M, N, A, LDA, VT, LDVT )
01795             ELSE
01796 *
01797 *              Copy A to U, generate Q
01798 *              (Cworkspace: need 2*M, prefer M+M*NB)
01799 *              (Rworkspace: 0)
01800 *
01801                CALL CLACPY( 'L', M, M, A, LDA, U, LDU )
01802                CALL CUNGBR( 'Q', M, M, N, U, LDU, WORK( ITAUQ ),
01803      $                      WORK( NWORK ), LWORK-NWORK+1, IERR )
01804 *
01805 *              Copy A to VT, generate P**H
01806 *              (Cworkspace: need 2*M, prefer M+M*NB)
01807 *              (Rworkspace: 0)
01808 *
01809                CALL CLACPY( 'U', M, N, A, LDA, VT, LDVT )
01810                CALL CUNGBR( 'P', N, N, M, VT, LDVT, WORK( ITAUP ),
01811      $                      WORK( NWORK ), LWORK-NWORK+1, IERR )
01812 *
01813 *              Perform bidiagonal SVD, computing left singular vectors
01814 *              of bidiagonal matrix in RWORK(IRU) and computing right
01815 *              singular vectors of bidiagonal matrix in RWORK(IRVT)
01816 *              (CWorkspace: need 0)
01817 *              (RWorkspace: need BDSPAC)
01818 *
01819                IRVT = NRWORK
01820                IRU = IRVT + M*M
01821                NRWORK = IRU + M*M
01822                CALL SBDSDC( 'L', 'I', M, S, RWORK( IE ), RWORK( IRU ),
01823      $                      M, RWORK( IRVT ), M, DUM, IDUM,
01824      $                      RWORK( NRWORK ), IWORK, INFO )
01825 *
01826 *              Multiply Q in U by real matrix RWORK(IRU), storing the
01827 *              result in A, copying to U
01828 *              (CWorkspace: need 0)
01829 *              (Rworkspace: need 3*M*M)
01830 *
01831                CALL CLACRM( M, M, U, LDU, RWORK( IRU ), M, A, LDA,
01832      $                      RWORK( NRWORK ) )
01833                CALL CLACPY( 'F', M, M, A, LDA, U, LDU )
01834 *
01835 *              Multiply real matrix RWORK(IRVT) by P**H in VT,
01836 *              storing the result in A, copying to VT
01837 *              (Cworkspace: need 0)
01838 *              (Rworkspace: need M*M+2*M*N)
01839 *
01840                CALL CLARCM( M, N, RWORK( IRVT ), M, VT, LDVT, A, LDA,
01841      $                      RWORK( NRWORK ) )
01842                CALL CLACPY( 'F', M, N, A, LDA, VT, LDVT )
01843             END IF
01844 *
01845          ELSE
01846 *
01847 *           N .LT. MNTHR2
01848 *
01849 *           Path 6t (N greater than M, but not much larger)
01850 *           Reduce to bidiagonal form without LQ decomposition
01851 *           Use CUNMBR to compute singular vectors
01852 *
01853             IE = 1
01854             NRWORK = IE + M
01855             ITAUQ = 1
01856             ITAUP = ITAUQ + M
01857             NWORK = ITAUP + M
01858 *
01859 *           Bidiagonalize A
01860 *           (CWorkspace: need 2*M+N, prefer 2*M+(M+N)*NB)
01861 *           (RWorkspace: M)
01862 *
01863             CALL CGEBRD( M, N, A, LDA, S, RWORK( IE ), WORK( ITAUQ ),
01864      $                   WORK( ITAUP ), WORK( NWORK ), LWORK-NWORK+1,
01865      $                   IERR )
01866             IF( WNTQN ) THEN
01867 *
01868 *              Compute singular values only
01869 *              (Cworkspace: 0)
01870 *              (Rworkspace: need BDSPAN)
01871 *
01872                CALL SBDSDC( 'L', 'N', M, S, RWORK( IE ), DUM, 1, DUM, 1,
01873      $                      DUM, IDUM, RWORK( NRWORK ), IWORK, INFO )
01874             ELSE IF( WNTQO ) THEN
01875                LDWKVT = M
01876                IVT = NWORK
01877                IF( LWORK.GE.M*N+3*M ) THEN
01878 *
01879 *                 WORK( IVT ) is M by N
01880 *
01881                   CALL CLASET( 'F', M, N, CZERO, CZERO, WORK( IVT ),
01882      $                         LDWKVT )
01883                   NWORK = IVT + LDWKVT*N
01884                ELSE
01885 *
01886 *                 WORK( IVT ) is M by CHUNK
01887 *
01888                   CHUNK = ( LWORK-3*M ) / M
01889                   NWORK = IVT + LDWKVT*CHUNK
01890                END IF
01891 *
01892 *              Perform bidiagonal SVD, computing left singular vectors
01893 *              of bidiagonal matrix in RWORK(IRU) and computing right
01894 *              singular vectors of bidiagonal matrix in RWORK(IRVT)
01895 *              (CWorkspace: need 0)
01896 *              (RWorkspace: need BDSPAC)
01897 *
01898                IRVT = NRWORK
01899                IRU = IRVT + M*M
01900                NRWORK = IRU + M*M
01901                CALL SBDSDC( 'L', 'I', M, S, RWORK( IE ), RWORK( IRU ),
01902      $                      M, RWORK( IRVT ), M, DUM, IDUM,
01903      $                      RWORK( NRWORK ), IWORK, INFO )
01904 *
01905 *              Copy real matrix RWORK(IRU) to complex matrix U
01906 *              Overwrite U by left singular vectors of A
01907 *              (Cworkspace: need 2*M, prefer M+M*NB)
01908 *              (Rworkspace: need 0)
01909 *
01910                CALL CLACP2( 'F', M, M, RWORK( IRU ), M, U, LDU )
01911                CALL CUNMBR( 'Q', 'L', 'N', M, M, N, A, LDA,
01912      $                      WORK( ITAUQ ), U, LDU, WORK( NWORK ),
01913      $                      LWORK-NWORK+1, IERR )
01914 *
01915                IF( LWORK.GE.M*N+3*M ) THEN
01916 *
01917 *              Copy real matrix RWORK(IRVT) to complex matrix WORK(IVT)
01918 *              Overwrite WORK(IVT) by right singular vectors of A,
01919 *              copying to A
01920 *              (Cworkspace: need M*N+2*M, prefer M*N+M+M*NB)
01921 *              (Rworkspace: need 0)
01922 *
01923                   CALL CLACP2( 'F', M, M, RWORK( IRVT ), M, WORK( IVT ),
01924      $                         LDWKVT )
01925                   CALL CUNMBR( 'P', 'R', 'C', M, N, M, A, LDA,
01926      $                         WORK( ITAUP ), WORK( IVT ), LDWKVT,
01927      $                         WORK( NWORK ), LWORK-NWORK+1, IERR )
01928                   CALL CLACPY( 'F', M, N, WORK( IVT ), LDWKVT, A, LDA )
01929                ELSE
01930 *
01931 *                 Generate P**H in A
01932 *                 (Cworkspace: need 2*M, prefer M+M*NB)
01933 *                 (Rworkspace: need 0)
01934 *
01935                   CALL CUNGBR( 'P', M, N, M, A, LDA, WORK( ITAUP ),
01936      $                         WORK( NWORK ), LWORK-NWORK+1, IERR )
01937 *
01938 *                 Multiply Q in A by real matrix RWORK(IRU), storing the
01939 *                 result in WORK(IU), copying to A
01940 *                 (CWorkspace: need M*M, prefer M*N)
01941 *                 (Rworkspace: need 3*M*M, prefer M*M+2*M*N)
01942 *
01943                   NRWORK = IRU
01944                   DO 60 I = 1, N, CHUNK
01945                      BLK = MIN( N-I+1, CHUNK )
01946                      CALL CLARCM( M, BLK, RWORK( IRVT ), M, A( 1, I ),
01947      $                            LDA, WORK( IVT ), LDWKVT,
01948      $                            RWORK( NRWORK ) )
01949                      CALL CLACPY( 'F', M, BLK, WORK( IVT ), LDWKVT,
01950      $                            A( 1, I ), LDA )
01951    60             CONTINUE
01952                END IF
01953             ELSE IF( WNTQS ) THEN
01954 *
01955 *              Perform bidiagonal SVD, computing left singular vectors
01956 *              of bidiagonal matrix in RWORK(IRU) and computing right
01957 *              singular vectors of bidiagonal matrix in RWORK(IRVT)
01958 *              (CWorkspace: need 0)
01959 *              (RWorkspace: need BDSPAC)
01960 *
01961                IRVT = NRWORK
01962                IRU = IRVT + M*M
01963                NRWORK = IRU + M*M
01964                CALL SBDSDC( 'L', 'I', M, S, RWORK( IE ), RWORK( IRU ),
01965      $                      M, RWORK( IRVT ), M, DUM, IDUM,
01966      $                      RWORK( NRWORK ), IWORK, INFO )
01967 *
01968 *              Copy real matrix RWORK(IRU) to complex matrix U
01969 *              Overwrite U by left singular vectors of A
01970 *              (CWorkspace: need 3*M, prefer 2*M+M*NB)
01971 *              (RWorkspace: M*M)
01972 *
01973                CALL CLACP2( 'F', M, M, RWORK( IRU ), M, U, LDU )
01974                CALL CUNMBR( 'Q', 'L', 'N', M, M, N, A, LDA,
01975      $                      WORK( ITAUQ ), U, LDU, WORK( NWORK ),
01976      $                      LWORK-NWORK+1, IERR )
01977 *
01978 *              Copy real matrix RWORK(IRVT) to complex matrix VT
01979 *              Overwrite VT by right singular vectors of A
01980 *              (CWorkspace: need 3*M, prefer 2*M+M*NB)
01981 *              (RWorkspace: M*M)
01982 *
01983                CALL CLASET( 'F', M, N, CZERO, CZERO, VT, LDVT )
01984                CALL CLACP2( 'F', M, M, RWORK( IRVT ), M, VT, LDVT )
01985                CALL CUNMBR( 'P', 'R', 'C', M, N, M, A, LDA,
01986      $                      WORK( ITAUP ), VT, LDVT, WORK( NWORK ),
01987      $                      LWORK-NWORK+1, IERR )
01988             ELSE
01989 *
01990 *              Perform bidiagonal SVD, computing left singular vectors
01991 *              of bidiagonal matrix in RWORK(IRU) and computing right
01992 *              singular vectors of bidiagonal matrix in RWORK(IRVT)
01993 *              (CWorkspace: need 0)
01994 *              (RWorkspace: need BDSPAC)
01995 *
01996                IRVT = NRWORK
01997                IRU = IRVT + M*M
01998                NRWORK = IRU + M*M
01999 *
02000                CALL SBDSDC( 'L', 'I', M, S, RWORK( IE ), RWORK( IRU ),
02001      $                      M, RWORK( IRVT ), M, DUM, IDUM,
02002      $                      RWORK( NRWORK ), IWORK, INFO )
02003 *
02004 *              Copy real matrix RWORK(IRU) to complex matrix U
02005 *              Overwrite U by left singular vectors of A
02006 *              (CWorkspace: need 3*M, prefer 2*M+M*NB)
02007 *              (RWorkspace: M*M)
02008 *
02009                CALL CLACP2( 'F', M, M, RWORK( IRU ), M, U, LDU )
02010                CALL CUNMBR( 'Q', 'L', 'N', M, M, N, A, LDA,
02011      $                      WORK( ITAUQ ), U, LDU, WORK( NWORK ),
02012      $                      LWORK-NWORK+1, IERR )
02013 *
02014 *              Set all of VT to identity matrix
02015 *
02016                CALL CLASET( 'F', N, N, CZERO, CONE, VT, LDVT )
02017 *
02018 *              Copy real matrix RWORK(IRVT) to complex matrix VT
02019 *              Overwrite VT by right singular vectors of A
02020 *              (CWorkspace: need 2*M+N, prefer 2*M+N*NB)
02021 *              (RWorkspace: M*M)
02022 *
02023                CALL CLACP2( 'F', M, M, RWORK( IRVT ), M, VT, LDVT )
02024                CALL CUNMBR( 'P', 'R', 'C', N, N, M, A, LDA,
02025      $                      WORK( ITAUP ), VT, LDVT, WORK( NWORK ),
02026      $                      LWORK-NWORK+1, IERR )
02027             END IF
02028 *
02029          END IF
02030 *
02031       END IF
02032 *
02033 *     Undo scaling if necessary
02034 *
02035       IF( ISCL.EQ.1 ) THEN
02036          IF( ANRM.GT.BIGNUM )
02037      $      CALL SLASCL( 'G', 0, 0, BIGNUM, ANRM, MINMN, 1, S, MINMN,
02038      $                   IERR )
02039          IF( INFO.NE.0 .AND. ANRM.GT.BIGNUM )
02040      $      CALL SLASCL( 'G', 0, 0, BIGNUM, ANRM, MINMN-1, 1,
02041      $                   RWORK( IE ), MINMN, IERR )
02042          IF( ANRM.LT.SMLNUM )
02043      $      CALL SLASCL( 'G', 0, 0, SMLNUM, ANRM, MINMN, 1, S, MINMN,
02044      $                   IERR )
02045          IF( INFO.NE.0 .AND. ANRM.LT.SMLNUM )
02046      $      CALL SLASCL( 'G', 0, 0, SMLNUM, ANRM, MINMN-1, 1,
02047      $                   RWORK( IE ), MINMN, IERR )
02048       END IF
02049 *
02050 *     Return optimal workspace in WORK(1)
02051 *
02052       WORK( 1 ) = MAXWRK
02053 *
02054       RETURN
02055 *
02056 *     End of CGESDD
02057 *
02058       END
 All Files Functions