LAPACK: SRC/zgesdd.f Source File

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