Diff for /rpl/lapack/lapack/zgesdd.f between versions 1.1 and 1.19

version 1.1, 2010/01/26 15:22:46 version 1.19, 2018/05/29 07:18:15
Line 1 Line 1
       SUBROUTINE ZGESDD( JOBZ, M, N, A, LDA, S, U, LDU, VT, LDVT, WORK,  *> \brief \b ZGESDD
      $                   LWORK, RWORK, IWORK, INFO )  
 *  *
 *  -- LAPACK driver routine (version 3.2) --  *  =========== DOCUMENTATION ===========
   *
   * Online html documentation available at
   *            http://www.netlib.org/lapack/explore-html/
   *
   *> \htmlonly
   *> Download ZGESDD + dependencies
   *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zgesdd.f">
   *> [TGZ]</a>
   *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zgesdd.f">
   *> [ZIP]</a>
   *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zgesdd.f">
   *> [TXT]</a>
   *> \endhtmlonly
   *
   *  Definition:
   *  ===========
   *
   *       SUBROUTINE ZGESDD( JOBZ, M, N, A, LDA, S, U, LDU, VT, LDVT,
   *                          WORK, LWORK, RWORK, IWORK, INFO )
   *
   *       .. Scalar Arguments ..
   *       CHARACTER          JOBZ
   *       INTEGER            INFO, LDA, LDU, LDVT, LWORK, M, N
   *       ..
   *       .. Array Arguments ..
   *       INTEGER            IWORK( * )
   *       DOUBLE PRECISION   RWORK( * ), S( * )
   *       COMPLEX*16         A( LDA, * ), U( LDU, * ), VT( LDVT, * ),
   *      $                   WORK( * )
   *       ..
   *
   *
   *> \par Purpose:
   *  =============
   *>
   *> \verbatim
   *>
   *> ZGESDD computes the singular value decomposition (SVD) of a complex
   *> M-by-N matrix A, optionally computing the left and/or right singular
   *> vectors, by using divide-and-conquer method. The SVD is written
   *>
   *>      A = U * SIGMA * conjugate-transpose(V)
   *>
   *> where SIGMA is an M-by-N matrix which is zero except for its
   *> min(m,n) diagonal elements, U is an M-by-M unitary matrix, and
   *> V is an N-by-N unitary matrix.  The diagonal elements of SIGMA
   *> are the singular values of A; they are real and non-negative, and
   *> are returned in descending order.  The first min(m,n) columns of
   *> U and V are the left and right singular vectors of A.
   *>
   *> Note that the routine returns VT = V**H, not V.
   *>
   *> The divide and conquer algorithm makes very mild assumptions about
   *> floating point arithmetic. It will work on machines with a guard
   *> digit in add/subtract, or on those binary machines without guard
   *> digits which subtract like the Cray X-MP, Cray Y-MP, Cray C-90, or
   *> Cray-2. It could conceivably fail on hexadecimal or decimal machines
   *> without guard digits, but we know of none.
   *> \endverbatim
   *
   *  Arguments:
   *  ==========
   *
   *> \param[in] JOBZ
   *> \verbatim
   *>          JOBZ is CHARACTER*1
   *>          Specifies options for computing all or part of the matrix U:
   *>          = 'A':  all M columns of U and all N rows of V**H are
   *>                  returned in the arrays U and VT;
   *>          = 'S':  the first min(M,N) columns of U and the first
   *>                  min(M,N) rows of V**H are returned in the arrays U
   *>                  and VT;
   *>          = 'O':  If M >= N, the first N columns of U are overwritten
   *>                  in the array A and all rows of V**H are returned in
   *>                  the array VT;
   *>                  otherwise, all columns of U are returned in the
   *>                  array U and the first M rows of V**H are overwritten
   *>                  in the array A;
   *>          = 'N':  no columns of U or rows of V**H are computed.
   *> \endverbatim
   *>
   *> \param[in] M
   *> \verbatim
   *>          M is INTEGER
   *>          The number of rows of the input matrix A.  M >= 0.
   *> \endverbatim
   *>
   *> \param[in] N
   *> \verbatim
   *>          N is INTEGER
   *>          The number of columns of the input matrix A.  N >= 0.
   *> \endverbatim
   *>
   *> \param[in,out] A
   *> \verbatim
   *>          A is COMPLEX*16 array, dimension (LDA,N)
   *>          On entry, the M-by-N matrix A.
   *>          On exit,
   *>          if JOBZ = 'O',  A is overwritten with the first N columns
   *>                          of U (the left singular vectors, stored
   *>                          columnwise) if M >= N;
   *>                          A is overwritten with the first M rows
   *>                          of V**H (the right singular vectors, stored
   *>                          rowwise) otherwise.
   *>          if JOBZ .ne. 'O', the contents of A are destroyed.
   *> \endverbatim
   *>
   *> \param[in] LDA
   *> \verbatim
   *>          LDA is INTEGER
   *>          The leading dimension of the array A.  LDA >= max(1,M).
   *> \endverbatim
   *>
   *> \param[out] S
   *> \verbatim
   *>          S is DOUBLE PRECISION array, dimension (min(M,N))
   *>          The singular values of A, sorted so that S(i) >= S(i+1).
   *> \endverbatim
   *>
   *> \param[out] U
   *> \verbatim
   *>          U is COMPLEX*16 array, dimension (LDU,UCOL)
   *>          UCOL = M if JOBZ = 'A' or JOBZ = 'O' and M < N;
   *>          UCOL = min(M,N) if JOBZ = 'S'.
   *>          If JOBZ = 'A' or JOBZ = 'O' and M < N, U contains the M-by-M
   *>          unitary matrix U;
   *>          if JOBZ = 'S', U contains the first min(M,N) columns of U
   *>          (the left singular vectors, stored columnwise);
   *>          if JOBZ = 'O' and M >= N, or JOBZ = 'N', U is not referenced.
   *> \endverbatim
   *>
   *> \param[in] LDU
   *> \verbatim
   *>          LDU is INTEGER
   *>          The leading dimension of the array U.  LDU >= 1;
   *>          if JOBZ = 'S' or 'A' or JOBZ = 'O' and M < N, LDU >= M.
   *> \endverbatim
   *>
   *> \param[out] VT
   *> \verbatim
   *>          VT is COMPLEX*16 array, dimension (LDVT,N)
   *>          If JOBZ = 'A' or JOBZ = 'O' and M >= N, VT contains the
   *>          N-by-N unitary matrix V**H;
   *>          if JOBZ = 'S', VT contains the first min(M,N) rows of
   *>          V**H (the right singular vectors, stored rowwise);
   *>          if JOBZ = 'O' and M < N, or JOBZ = 'N', VT is not referenced.
   *> \endverbatim
   *>
   *> \param[in] LDVT
   *> \verbatim
   *>          LDVT is INTEGER
   *>          The leading dimension of the array VT.  LDVT >= 1;
   *>          if JOBZ = 'A' or JOBZ = 'O' and M >= N, LDVT >= N;
   *>          if JOBZ = 'S', LDVT >= min(M,N).
   *> \endverbatim
   *>
   *> \param[out] WORK
   *> \verbatim
   *>          WORK is COMPLEX*16 array, dimension (MAX(1,LWORK))
   *>          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
   *> \endverbatim
   *>
   *> \param[in] LWORK
   *> \verbatim
   *>          LWORK is INTEGER
   *>          The dimension of the array WORK. LWORK >= 1.
   *>          If LWORK = -1, a workspace query is assumed.  The optimal
   *>          size for the WORK array is calculated and stored in WORK(1),
   *>          and no other work except argument checking is performed.
   *>
   *>          Let mx = max(M,N) and mn = min(M,N).
   *>          If JOBZ = 'N', LWORK >= 2*mn + mx.
   *>          If JOBZ = 'O', LWORK >= 2*mn*mn + 2*mn + mx.
   *>          If JOBZ = 'S', LWORK >=   mn*mn + 3*mn.
   *>          If JOBZ = 'A', LWORK >=   mn*mn + 2*mn + mx.
   *>          These are not tight minimums in all cases; see comments inside code.
   *>          For good performance, LWORK should generally be larger;
   *>          a query is recommended.
   *> \endverbatim
   *>
   *> \param[out] RWORK
   *> \verbatim
   *>          RWORK is DOUBLE PRECISION array, dimension (MAX(1,LRWORK))
   *>          Let mx = max(M,N) and mn = min(M,N).
   *>          If JOBZ = 'N',    LRWORK >= 5*mn (LAPACK <= 3.6 needs 7*mn);
   *>          else if mx >> mn, LRWORK >= 5*mn*mn + 5*mn;
   *>          else              LRWORK >= max( 5*mn*mn + 5*mn,
   *>                                           2*mx*mn + 2*mn*mn + mn ).
   *> \endverbatim
   *>
   *> \param[out] IWORK
   *> \verbatim
   *>          IWORK is INTEGER array, dimension (8*min(M,N))
   *> \endverbatim
   *>
   *> \param[out] INFO
   *> \verbatim
   *>          INFO is INTEGER
   *>          = 0:  successful exit.
   *>          < 0:  if INFO = -i, the i-th argument had an illegal value.
   *>          > 0:  The updating process of DBDSDC did not converge.
   *> \endverbatim
   *
   *  Authors:
   *  ========
   *
   *> \author Univ. of Tennessee
   *> \author Univ. of California Berkeley
   *> \author Univ. of Colorado Denver
   *> \author NAG Ltd.
   *
   *> \date June 2016
   *
   *> \ingroup complex16GEsing
   *
   *> \par Contributors:
   *  ==================
   *>
   *>     Ming Gu and Huan Ren, Computer Science Division, University of
   *>     California at Berkeley, USA
   *>
   *  =====================================================================
         SUBROUTINE ZGESDD( JOBZ, M, N, A, LDA, S, U, LDU, VT, LDVT,
        $                   WORK, LWORK, RWORK, IWORK, INFO )
         implicit none
   *
   *  -- LAPACK driver routine (version 3.7.0) --
 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --  *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--  *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
 *     November 2006  *     June 2016
 *     8-15-00:  Improve consistency of WS calculations (eca)  
 *  *
 *     .. Scalar Arguments ..  *     .. Scalar Arguments ..
       CHARACTER          JOBZ        CHARACTER          JOBZ
Line 18 Line 243
      $                   WORK( * )       $                   WORK( * )
 *     ..  *     ..
 *  *
 *  Purpose  
 *  =======  
 *  
 *  ZGESDD computes the singular value decomposition (SVD) of a complex  
 *  M-by-N matrix A, optionally computing the left and/or right singular  
 *  vectors, by using divide-and-conquer method. The SVD is written  
 *  
 *       A = U * SIGMA * conjugate-transpose(V)  
 *  
 *  where SIGMA is an M-by-N matrix which is zero except for its  
 *  min(m,n) diagonal elements, U is an M-by-M unitary matrix, and  
 *  V is an N-by-N unitary matrix.  The diagonal elements of SIGMA  
 *  are the singular values of A; they are real and non-negative, and  
 *  are returned in descending order.  The first min(m,n) columns of  
 *  U and V are the left and right singular vectors of A.  
 *  
 *  Note that the routine returns VT = V**H, not V.  
 *  
 *  The divide and conquer algorithm makes very mild assumptions about  
 *  floating point arithmetic. It will work on machines with a guard  
 *  digit in add/subtract, or on those binary machines without guard  
 *  digits which subtract like the Cray X-MP, Cray Y-MP, Cray C-90, or  
 *  Cray-2. It could conceivably fail on hexadecimal or decimal machines  
 *  without guard digits, but we know of none.  
 *  
 *  Arguments  
 *  =========  
 *  
 *  JOBZ    (input) CHARACTER*1  
 *          Specifies options for computing all or part of the matrix U:  
 *          = 'A':  all M columns of U and all N rows of V**H are  
 *                  returned in the arrays U and VT;  
 *          = 'S':  the first min(M,N) columns of U and the first  
 *                  min(M,N) rows of V**H are returned in the arrays U  
 *                  and VT;  
 *          = 'O':  If M >= N, the first N columns of U are overwritten  
 *                  in the array A and all rows of V**H are returned in  
 *                  the array VT;  
 *                  otherwise, all columns of U are returned in the  
 *                  array U and the first M rows of V**H are overwritten  
 *                  in the array A;  
 *          = 'N':  no columns of U or rows of V**H are computed.  
 *  
 *  M       (input) INTEGER  
 *          The number of rows of the input matrix A.  M >= 0.  
 *  
 *  N       (input) INTEGER  
 *          The number of columns of the input matrix A.  N >= 0.  
 *  
 *  A       (input/output) COMPLEX*16 array, dimension (LDA,N)  
 *          On entry, the M-by-N matrix A.  
 *          On exit,  
 *          if JOBZ = 'O',  A is overwritten with the first N columns  
 *                          of U (the left singular vectors, stored  
 *                          columnwise) if M >= N;  
 *                          A is overwritten with the first M rows  
 *                          of V**H (the right singular vectors, stored  
 *                          rowwise) otherwise.  
 *          if JOBZ .ne. 'O', the contents of A are destroyed.  
 *  
 *  LDA     (input) INTEGER  
 *          The leading dimension of the array A.  LDA >= max(1,M).  
 *  
 *  S       (output) DOUBLE PRECISION array, dimension (min(M,N))  
 *          The singular values of A, sorted so that S(i) >= S(i+1).  
 *  
 *  U       (output) COMPLEX*16 array, dimension (LDU,UCOL)  
 *          UCOL = M if JOBZ = 'A' or JOBZ = 'O' and M < N;  
 *          UCOL = min(M,N) if JOBZ = 'S'.  
 *          If JOBZ = 'A' or JOBZ = 'O' and M < N, U contains the M-by-M  
 *          unitary matrix U;  
 *          if JOBZ = 'S', U contains the first min(M,N) columns of U  
 *          (the left singular vectors, stored columnwise);  
 *          if JOBZ = 'O' and M >= N, or JOBZ = 'N', U is not referenced.  
 *  
 *  LDU     (input) INTEGER  
 *          The leading dimension of the array U.  LDU >= 1; if  
 *          JOBZ = 'S' or 'A' or JOBZ = 'O' and M < N, LDU >= M.  
 *  
 *  VT      (output) COMPLEX*16 array, dimension (LDVT,N)  
 *          If JOBZ = 'A' or JOBZ = 'O' and M >= N, VT contains the  
 *          N-by-N unitary matrix V**H;  
 *          if JOBZ = 'S', VT contains the first min(M,N) rows of  
 *          V**H (the right singular vectors, stored rowwise);  
 *          if JOBZ = 'O' and M < N, or JOBZ = 'N', VT is not referenced.  
 *  
 *  LDVT    (input) INTEGER  
 *          The leading dimension of the array VT.  LDVT >= 1; if  
 *          JOBZ = 'A' or JOBZ = 'O' and M >= N, LDVT >= N;  
 *          if JOBZ = 'S', LDVT >= min(M,N).  
 *  
 *  WORK    (workspace/output) COMPLEX*16 array, dimension (MAX(1,LWORK))  
 *          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.  
 *  
 *  LWORK   (input) INTEGER  
 *          The dimension of the array WORK. LWORK >= 1.  
 *          if JOBZ = 'N', LWORK >= 2*min(M,N)+max(M,N).  
 *          if JOBZ = 'O',  
 *                LWORK >= 2*min(M,N)*min(M,N)+2*min(M,N)+max(M,N).  
 *          if JOBZ = 'S' or 'A',  
 *                LWORK >= min(M,N)*min(M,N)+2*min(M,N)+max(M,N).  
 *          For good performance, LWORK should generally be larger.  
 *  
 *          If LWORK = -1, a workspace query is assumed.  The optimal  
 *          size for the WORK array is calculated and stored in WORK(1),  
 *          and no other work except argument checking is performed.  
 *  
 *  RWORK   (workspace) DOUBLE PRECISION array, dimension (MAX(1,LRWORK))  
 *          If JOBZ = 'N', LRWORK >= 5*min(M,N).  
 *          Otherwise, LRWORK >= 5*min(M,N)*min(M,N) + 7*min(M,N)  
 *  
 *  IWORK   (workspace) INTEGER array, dimension (8*min(M,N))  
 *  
 *  INFO    (output) INTEGER  
 *          = 0:  successful exit.  
 *          < 0:  if INFO = -i, the i-th argument had an illegal value.  
 *          > 0:  The updating process of DBDSDC did not converge.  
 *  
 *  Further Details  
 *  ===============  
 *  
 *  Based on contributions by  
 *     Ming Gu and Huan Ren, Computer Science Division, University of  
 *     California at Berkeley, USA  
 *  
 *  =====================================================================  *  =====================================================================
 *  *
 *     .. Parameters ..  *     .. Parameters ..
       INTEGER            LQUERV  
       PARAMETER          ( LQUERV = -1 )  
       COMPLEX*16         CZERO, CONE        COMPLEX*16         CZERO, CONE
       PARAMETER          ( CZERO = ( 0.0D+0, 0.0D+0 ),        PARAMETER          ( CZERO = ( 0.0D+0, 0.0D+0 ),
      $                   CONE = ( 1.0D+0, 0.0D+0 ) )       $                   CONE = ( 1.0D+0, 0.0D+0 ) )
Line 155 Line 253
       PARAMETER          ( ZERO = 0.0D+0, ONE = 1.0D+0 )        PARAMETER          ( ZERO = 0.0D+0, ONE = 1.0D+0 )
 *     ..  *     ..
 *     .. Local Scalars ..  *     .. Local Scalars ..
       LOGICAL            WNTQA, WNTQAS, WNTQN, WNTQO, WNTQS        LOGICAL            LQUERY, WNTQA, WNTQAS, WNTQN, WNTQO, WNTQS
       INTEGER            BLK, CHUNK, I, IE, IERR, IL, IR, IRU, IRVT,        INTEGER            BLK, CHUNK, I, IE, IERR, IL, IR, IRU, IRVT,
      $                   ISCL, ITAU, ITAUP, ITAUQ, IU, IVT, LDWKVT,       $                   ISCL, ITAU, ITAUP, ITAUQ, IU, IVT, LDWKVT,
      $                   LDWRKL, LDWRKR, LDWRKU, MAXWRK, MINMN, MINWRK,       $                   LDWRKL, LDWRKR, LDWRKU, MAXWRK, MINMN, MINWRK,
      $                   MNTHR1, MNTHR2, NRWORK, NWORK, WRKBL       $                   MNTHR1, MNTHR2, NRWORK, NWORK, WRKBL
         INTEGER            LWORK_ZGEBRD_MN, LWORK_ZGEBRD_MM,
        $                   LWORK_ZGEBRD_NN, LWORK_ZGELQF_MN,
        $                   LWORK_ZGEQRF_MN,
        $                   LWORK_ZUNGBR_P_MN, LWORK_ZUNGBR_P_NN,
        $                   LWORK_ZUNGBR_Q_MN, LWORK_ZUNGBR_Q_MM,
        $                   LWORK_ZUNGLQ_MN, LWORK_ZUNGLQ_NN,
        $                   LWORK_ZUNGQR_MM, LWORK_ZUNGQR_MN,
        $                   LWORK_ZUNMBR_PRC_MM, LWORK_ZUNMBR_QLN_MM,
        $                   LWORK_ZUNMBR_PRC_MN, LWORK_ZUNMBR_QLN_MN,
        $                   LWORK_ZUNMBR_PRC_NN, LWORK_ZUNMBR_QLN_NN
       DOUBLE PRECISION   ANRM, BIGNUM, EPS, SMLNUM        DOUBLE PRECISION   ANRM, BIGNUM, EPS, SMLNUM
 *     ..  *     ..
 *     .. Local Arrays ..  *     .. Local Arrays ..
       INTEGER            IDUM( 1 )        INTEGER            IDUM( 1 )
       DOUBLE PRECISION   DUM( 1 )        DOUBLE PRECISION   DUM( 1 )
         COMPLEX*16         CDUM( 1 )
 *     ..  *     ..
 *     .. External Subroutines ..  *     .. External Subroutines ..
       EXTERNAL           DBDSDC, DLASCL, XERBLA, ZGEBRD, ZGELQF, ZGEMM,        EXTERNAL           DBDSDC, DLASCL, XERBLA, ZGEBRD, ZGELQF, ZGEMM,
Line 173 Line 282
 *     ..  *     ..
 *     .. External Functions ..  *     .. External Functions ..
       LOGICAL            LSAME        LOGICAL            LSAME
       INTEGER            ILAENV  
       DOUBLE PRECISION   DLAMCH, ZLANGE        DOUBLE PRECISION   DLAMCH, ZLANGE
       EXTERNAL           LSAME, ILAENV, DLAMCH, ZLANGE        EXTERNAL           LSAME, DLAMCH, ZLANGE
 *     ..  *     ..
 *     .. Intrinsic Functions ..  *     .. Intrinsic Functions ..
       INTRINSIC          INT, MAX, MIN, SQRT        INTRINSIC          INT, MAX, MIN, SQRT
Line 184 Line 292
 *  *
 *     Test the input arguments  *     Test the input arguments
 *  *
       INFO = 0        INFO   = 0
       MINMN = MIN( M, N )        MINMN  = MIN( M, N )
       MNTHR1 = INT( MINMN*17.0D0 / 9.0D0 )        MNTHR1 = INT( MINMN*17.0D0 / 9.0D0 )
       MNTHR2 = INT( MINMN*5.0D0 / 3.0D0 )        MNTHR2 = INT( MINMN*5.0D0 / 3.0D0 )
       WNTQA = LSAME( JOBZ, 'A' )        WNTQA  = LSAME( JOBZ, 'A' )
       WNTQS = LSAME( JOBZ, 'S' )        WNTQS  = LSAME( JOBZ, 'S' )
       WNTQAS = WNTQA .OR. WNTQS        WNTQAS = WNTQA .OR. WNTQS
       WNTQO = LSAME( JOBZ, 'O' )        WNTQO  = LSAME( JOBZ, 'O' )
       WNTQN = LSAME( JOBZ, 'N' )        WNTQN  = LSAME( JOBZ, 'N' )
         LQUERY = ( LWORK.EQ.-1 )
       MINWRK = 1        MINWRK = 1
       MAXWRK = 1        MAXWRK = 1
 *  *
Line 214 Line 323
       END IF        END IF
 *  *
 *     Compute workspace  *     Compute workspace
 *      (Note: Comments in the code beginning "Workspace:" describe the  *       Note: Comments in the code beginning "Workspace:" describe the
 *       minimal amount of workspace needed at that point in the code,  *       minimal amount of workspace allocated at that point in the code,
 *       as well as the preferred amount for good performance.  *       as well as the preferred amount for good performance.
 *       CWorkspace refers to complex workspace, and RWorkspace to  *       CWorkspace refers to complex workspace, and RWorkspace to
 *       real workspace. NB refers to the optimal block size for the  *       real workspace. NB refers to the optimal block size for the
 *       immediately following subroutine, as returned by ILAENV.)  *       immediately following subroutine, as returned by ILAENV.)
 *  *
       IF( INFO.EQ.0 .AND. M.GT.0 .AND. N.GT.0 ) THEN        IF( INFO.EQ.0 ) THEN
          IF( M.GE.N ) THEN           MINWRK = 1
            MAXWRK = 1
            IF( M.GE.N .AND. MINMN.GT.0 ) THEN
 *  *
 *           There is no complex work space needed for bidiagonal SVD  *           There is no complex work space needed for bidiagonal SVD
 *           The real work space needed for bidiagonal SVD is BDSPAC  *           The real work space needed for bidiagonal SVD (dbdsdc) is
 *           for computing singular values and singular vectors; BDSPAN  *           BDSPAC = 3*N*N + 4*N for singular values and vectors;
 *           for computing singular values only.  *           BDSPAC = 4*N         for singular values only;
 *           BDSPAC = 5*N*N + 7*N  *           not including e, RU, and RVT matrices.
 *           BDSPAN = MAX(7*N+4, 3*N+2+SMLSIZ*(SMLSIZ+8))  *
   *           Compute space preferred for each routine
               CALL ZGEBRD( M, N, CDUM(1), M, DUM(1), DUM(1), CDUM(1),
        $                   CDUM(1), CDUM(1), -1, IERR )
               LWORK_ZGEBRD_MN = INT( CDUM(1) )
   *
               CALL ZGEBRD( N, N, CDUM(1), N, DUM(1), DUM(1), CDUM(1),
        $                   CDUM(1), CDUM(1), -1, IERR )
               LWORK_ZGEBRD_NN = INT( CDUM(1) )
   *
               CALL ZGEQRF( M, N, CDUM(1), M, CDUM(1), CDUM(1), -1, IERR )
               LWORK_ZGEQRF_MN = INT( CDUM(1) )
   *
               CALL ZUNGBR( 'P', N, N, N, CDUM(1), N, CDUM(1), CDUM(1),
        $                   -1, IERR )
               LWORK_ZUNGBR_P_NN = INT( CDUM(1) )
   *
               CALL ZUNGBR( 'Q', M, M, N, CDUM(1), M, CDUM(1), CDUM(1),
        $                   -1, IERR )
               LWORK_ZUNGBR_Q_MM = INT( CDUM(1) )
   *
               CALL ZUNGBR( 'Q', M, N, N, CDUM(1), M, CDUM(1), CDUM(1),
        $                   -1, IERR )
               LWORK_ZUNGBR_Q_MN = INT( CDUM(1) )
   *
               CALL ZUNGQR( M, M, N, CDUM(1), M, CDUM(1), CDUM(1),
        $                   -1, IERR )
               LWORK_ZUNGQR_MM = INT( CDUM(1) )
   *
               CALL ZUNGQR( M, N, N, CDUM(1), M, CDUM(1), CDUM(1),
        $                   -1, IERR )
               LWORK_ZUNGQR_MN = INT( CDUM(1) )
   *
               CALL ZUNMBR( 'P', 'R', 'C', N, N, N, CDUM(1), N, CDUM(1),
        $                   CDUM(1), N, CDUM(1), -1, IERR )
               LWORK_ZUNMBR_PRC_NN = INT( CDUM(1) )
   *
               CALL ZUNMBR( 'Q', 'L', 'N', M, M, N, CDUM(1), M, CDUM(1),
        $                   CDUM(1), M, CDUM(1), -1, IERR )
               LWORK_ZUNMBR_QLN_MM = INT( CDUM(1) )
   *
               CALL ZUNMBR( 'Q', 'L', 'N', M, N, N, CDUM(1), M, CDUM(1),
        $                   CDUM(1), M, CDUM(1), -1, IERR )
               LWORK_ZUNMBR_QLN_MN = INT( CDUM(1) )
   *
               CALL ZUNMBR( 'Q', 'L', 'N', N, N, N, CDUM(1), N, CDUM(1),
        $                   CDUM(1), N, CDUM(1), -1, IERR )
               LWORK_ZUNMBR_QLN_NN = INT( CDUM(1) )
 *  *
             IF( M.GE.MNTHR1 ) THEN              IF( M.GE.MNTHR1 ) THEN
                IF( WNTQN ) THEN                 IF( WNTQN ) THEN
 *  *
 *                 Path 1 (M much larger than N, JOBZ='N')  *                 Path 1 (M >> N, JOBZ='N')
 *  *
                   MAXWRK = N + N*ILAENV( 1, 'ZGEQRF', ' ', M, N, -1,                    MAXWRK = N + LWORK_ZGEQRF_MN
      $                     -1 )                    MAXWRK = MAX( MAXWRK, 2*N + LWORK_ZGEBRD_NN )
                   MAXWRK = MAX( MAXWRK, 2*N+2*N*  
      $                     ILAENV( 1, 'ZGEBRD', ' ', N, N, -1, -1 ) )  
                   MINWRK = 3*N                    MINWRK = 3*N
                ELSE IF( WNTQO ) THEN                 ELSE IF( WNTQO ) THEN
 *  *
 *                 Path 2 (M much larger than N, JOBZ='O')  *                 Path 2 (M >> N, JOBZ='O')
 *  *
                   WRKBL = N + N*ILAENV( 1, 'ZGEQRF', ' ', M, N, -1, -1 )                    WRKBL = N + LWORK_ZGEQRF_MN
                   WRKBL = MAX( WRKBL, N+N*ILAENV( 1, 'ZUNGQR', ' ', M,                    WRKBL = MAX( WRKBL,   N + LWORK_ZUNGQR_MN )
      $                    N, N, -1 ) )                    WRKBL = MAX( WRKBL, 2*N + LWORK_ZGEBRD_NN )
                   WRKBL = MAX( WRKBL, 2*N+2*N*                    WRKBL = MAX( WRKBL, 2*N + LWORK_ZUNMBR_QLN_NN )
      $                    ILAENV( 1, 'ZGEBRD', ' ', N, N, -1, -1 ) )                    WRKBL = MAX( WRKBL, 2*N + LWORK_ZUNMBR_PRC_NN )
                   WRKBL = MAX( WRKBL, 2*N+N*  
      $                    ILAENV( 1, 'ZUNMBR', 'QLN', N, N, N, -1 ) )  
                   WRKBL = MAX( WRKBL, 2*N+N*  
      $                    ILAENV( 1, 'ZUNMBR', 'PRC', N, N, N, -1 ) )  
                   MAXWRK = M*N + N*N + WRKBL                    MAXWRK = M*N + N*N + WRKBL
                   MINWRK = 2*N*N + 3*N                    MINWRK = 2*N*N + 3*N
                ELSE IF( WNTQS ) THEN                 ELSE IF( WNTQS ) THEN
 *  *
 *                 Path 3 (M much larger than N, JOBZ='S')  *                 Path 3 (M >> N, JOBZ='S')
 *  *
                   WRKBL = N + N*ILAENV( 1, 'ZGEQRF', ' ', M, N, -1, -1 )                    WRKBL = N + LWORK_ZGEQRF_MN
                   WRKBL = MAX( WRKBL, N+N*ILAENV( 1, 'ZUNGQR', ' ', M,                    WRKBL = MAX( WRKBL,   N + LWORK_ZUNGQR_MN )
      $                    N, N, -1 ) )                    WRKBL = MAX( WRKBL, 2*N + LWORK_ZGEBRD_NN )
                   WRKBL = MAX( WRKBL, 2*N+2*N*                    WRKBL = MAX( WRKBL, 2*N + LWORK_ZUNMBR_QLN_NN )
      $                    ILAENV( 1, 'ZGEBRD', ' ', N, N, -1, -1 ) )                    WRKBL = MAX( WRKBL, 2*N + LWORK_ZUNMBR_PRC_NN )
                   WRKBL = MAX( WRKBL, 2*N+N*  
      $                    ILAENV( 1, 'ZUNMBR', 'QLN', N, N, N, -1 ) )  
                   WRKBL = MAX( WRKBL, 2*N+N*  
      $                    ILAENV( 1, 'ZUNMBR', 'PRC', N, N, N, -1 ) )  
                   MAXWRK = N*N + WRKBL                    MAXWRK = N*N + WRKBL
                   MINWRK = N*N + 3*N                    MINWRK = N*N + 3*N
                ELSE IF( WNTQA ) THEN                 ELSE IF( WNTQA ) THEN
 *  *
 *                 Path 4 (M much larger than N, JOBZ='A')  *                 Path 4 (M >> N, JOBZ='A')
 *  *
                   WRKBL = N + N*ILAENV( 1, 'ZGEQRF', ' ', M, N, -1, -1 )                    WRKBL = N + LWORK_ZGEQRF_MN
                   WRKBL = MAX( WRKBL, N+M*ILAENV( 1, 'ZUNGQR', ' ', M,                    WRKBL = MAX( WRKBL,   N + LWORK_ZUNGQR_MM )
      $                    M, N, -1 ) )                    WRKBL = MAX( WRKBL, 2*N + LWORK_ZGEBRD_NN )
                   WRKBL = MAX( WRKBL, 2*N+2*N*                    WRKBL = MAX( WRKBL, 2*N + LWORK_ZUNMBR_QLN_NN )
      $                    ILAENV( 1, 'ZGEBRD', ' ', N, N, -1, -1 ) )                    WRKBL = MAX( WRKBL, 2*N + LWORK_ZUNMBR_PRC_NN )
                   WRKBL = MAX( WRKBL, 2*N+N*  
      $                    ILAENV( 1, 'ZUNMBR', 'QLN', N, N, N, -1 ) )  
                   WRKBL = MAX( WRKBL, 2*N+N*  
      $                    ILAENV( 1, 'ZUNMBR', 'PRC', N, N, N, -1 ) )  
                   MAXWRK = N*N + WRKBL                    MAXWRK = N*N + WRKBL
                   MINWRK = N*N + 2*N + M                    MINWRK = N*N + MAX( 3*N, N + M )
                END IF                 END IF
             ELSE IF( M.GE.MNTHR2 ) THEN              ELSE IF( M.GE.MNTHR2 ) THEN
 *  *
 *              Path 5 (M much larger than N, but not as much as MNTHR1)  *              Path 5 (M >> N, but not as much as MNTHR1)
 *  *
                MAXWRK = 2*N + ( M+N )*ILAENV( 1, 'ZGEBRD', ' ', M, N,                 MAXWRK = 2*N + LWORK_ZGEBRD_MN
      $                  -1, -1 )  
                MINWRK = 2*N + M                 MINWRK = 2*N + M
                IF( WNTQO ) THEN                 IF( WNTQO ) THEN
                   MAXWRK = MAX( MAXWRK, 2*N+N*  *                 Path 5o (M >> N, JOBZ='O')
      $                     ILAENV( 1, 'ZUNGBR', 'P', N, N, N, -1 ) )                    MAXWRK = MAX( MAXWRK, 2*N + LWORK_ZUNGBR_P_NN )
                   MAXWRK = MAX( MAXWRK, 2*N+N*                    MAXWRK = MAX( MAXWRK, 2*N + LWORK_ZUNGBR_Q_MN )
      $                     ILAENV( 1, 'ZUNGBR', 'Q', M, N, N, -1 ) )  
                   MAXWRK = MAXWRK + M*N                    MAXWRK = MAXWRK + M*N
                   MINWRK = MINWRK + N*N                    MINWRK = MINWRK + N*N
                ELSE IF( WNTQS ) THEN                 ELSE IF( WNTQS ) THEN
                   MAXWRK = MAX( MAXWRK, 2*N+N*  *                 Path 5s (M >> N, JOBZ='S')
      $                     ILAENV( 1, 'ZUNGBR', 'P', N, N, N, -1 ) )                    MAXWRK = MAX( MAXWRK, 2*N + LWORK_ZUNGBR_P_NN )
                   MAXWRK = MAX( MAXWRK, 2*N+N*                    MAXWRK = MAX( MAXWRK, 2*N + LWORK_ZUNGBR_Q_MN )
      $                     ILAENV( 1, 'ZUNGBR', 'Q', M, N, N, -1 ) )  
                ELSE IF( WNTQA ) THEN                 ELSE IF( WNTQA ) THEN
                   MAXWRK = MAX( MAXWRK, 2*N+N*  *                 Path 5a (M >> N, JOBZ='A')
      $                     ILAENV( 1, 'ZUNGBR', 'P', N, N, N, -1 ) )                    MAXWRK = MAX( MAXWRK, 2*N + LWORK_ZUNGBR_P_NN )
                   MAXWRK = MAX( MAXWRK, 2*N+M*                    MAXWRK = MAX( MAXWRK, 2*N + LWORK_ZUNGBR_Q_MM )
      $                     ILAENV( 1, 'ZUNGBR', 'Q', M, M, N, -1 ) )  
                END IF                 END IF
             ELSE              ELSE
 *  *
 *              Path 6 (M at least N, but not much larger)  *              Path 6 (M >= N, but not much larger)
 *  *
                MAXWRK = 2*N + ( M+N )*ILAENV( 1, 'ZGEBRD', ' ', M, N,                 MAXWRK = 2*N + LWORK_ZGEBRD_MN
      $                  -1, -1 )  
                MINWRK = 2*N + M                 MINWRK = 2*N + M
                IF( WNTQO ) THEN                 IF( WNTQO ) THEN
                   MAXWRK = MAX( MAXWRK, 2*N+N*  *                 Path 6o (M >= N, JOBZ='O')
      $                     ILAENV( 1, 'ZUNMBR', 'PRC', N, N, N, -1 ) )                    MAXWRK = MAX( MAXWRK, 2*N + LWORK_ZUNMBR_PRC_NN )
                   MAXWRK = MAX( MAXWRK, 2*N+N*                    MAXWRK = MAX( MAXWRK, 2*N + LWORK_ZUNMBR_QLN_MN )
      $                     ILAENV( 1, 'ZUNMBR', 'QLN', M, N, N, -1 ) )  
                   MAXWRK = MAXWRK + M*N                    MAXWRK = MAXWRK + M*N
                   MINWRK = MINWRK + N*N                    MINWRK = MINWRK + N*N
                ELSE IF( WNTQS ) THEN                 ELSE IF( WNTQS ) THEN
                   MAXWRK = MAX( MAXWRK, 2*N+N*  *                 Path 6s (M >= N, JOBZ='S')
      $                     ILAENV( 1, 'ZUNMBR', 'PRC', N, N, N, -1 ) )                    MAXWRK = MAX( MAXWRK, 2*N + LWORK_ZUNMBR_QLN_MN )
                   MAXWRK = MAX( MAXWRK, 2*N+N*                    MAXWRK = MAX( MAXWRK, 2*N + LWORK_ZUNMBR_PRC_NN )
      $                     ILAENV( 1, 'ZUNMBR', 'QLN', M, N, N, -1 ) )  
                ELSE IF( WNTQA ) THEN                 ELSE IF( WNTQA ) THEN
                   MAXWRK = MAX( MAXWRK, 2*N+N*  *                 Path 6a (M >= N, JOBZ='A')
      $                     ILAENV( 1, 'ZUNGBR', 'PRC', N, N, N, -1 ) )                    MAXWRK = MAX( MAXWRK, 2*N + LWORK_ZUNMBR_QLN_MM )
                   MAXWRK = MAX( MAXWRK, 2*N+M*                    MAXWRK = MAX( MAXWRK, 2*N + LWORK_ZUNMBR_PRC_NN )
      $                     ILAENV( 1, 'ZUNGBR', 'QLN', M, M, N, -1 ) )  
                END IF                 END IF
             END IF              END IF
          ELSE           ELSE IF( MINMN.GT.0 ) THEN
 *  *
 *           There is no complex work space needed for bidiagonal SVD  *           There is no complex work space needed for bidiagonal SVD
 *           The real work space needed for bidiagonal SVD is BDSPAC  *           The real work space needed for bidiagonal SVD (dbdsdc) is
 *           for computing singular values and singular vectors; BDSPAN  *           BDSPAC = 3*M*M + 4*M for singular values and vectors;
 *           for computing singular values only.  *           BDSPAC = 4*M         for singular values only;
 *           BDSPAC = 5*M*M + 7*M  *           not including e, RU, and RVT matrices.
 *           BDSPAN = MAX(7*M+4, 3*M+2+SMLSIZ*(SMLSIZ+8))  *
   *           Compute space preferred for each routine
               CALL ZGEBRD( M, N, CDUM(1), M, DUM(1), DUM(1), CDUM(1),
        $                   CDUM(1), CDUM(1), -1, IERR )
               LWORK_ZGEBRD_MN = INT( CDUM(1) )
   *
               CALL ZGEBRD( M, M, CDUM(1), M, DUM(1), DUM(1), CDUM(1),
        $                   CDUM(1), CDUM(1), -1, IERR )
               LWORK_ZGEBRD_MM = INT( CDUM(1) )
   *
               CALL ZGELQF( M, N, CDUM(1), M, CDUM(1), CDUM(1), -1, IERR )
               LWORK_ZGELQF_MN = INT( CDUM(1) )
   *
               CALL ZUNGBR( 'P', M, N, M, CDUM(1), M, CDUM(1), CDUM(1),
        $                   -1, IERR )
               LWORK_ZUNGBR_P_MN = INT( CDUM(1) )
   *
               CALL ZUNGBR( 'P', N, N, M, CDUM(1), N, CDUM(1), CDUM(1),
        $                   -1, IERR )
               LWORK_ZUNGBR_P_NN = INT( CDUM(1) )
   *
               CALL ZUNGBR( 'Q', M, M, N, CDUM(1), M, CDUM(1), CDUM(1),
        $                   -1, IERR )
               LWORK_ZUNGBR_Q_MM = INT( CDUM(1) )
   *
               CALL ZUNGLQ( M, N, M, CDUM(1), M, CDUM(1), CDUM(1),
        $                   -1, IERR )
               LWORK_ZUNGLQ_MN = INT( CDUM(1) )
   *
               CALL ZUNGLQ( N, N, M, CDUM(1), N, CDUM(1), CDUM(1),
        $                   -1, IERR )
               LWORK_ZUNGLQ_NN = INT( CDUM(1) )
   *
               CALL ZUNMBR( 'P', 'R', 'C', M, M, M, CDUM(1), M, CDUM(1),
        $                   CDUM(1), M, CDUM(1), -1, IERR )
               LWORK_ZUNMBR_PRC_MM = INT( CDUM(1) )
   *
               CALL ZUNMBR( 'P', 'R', 'C', M, N, M, CDUM(1), M, CDUM(1),
        $                   CDUM(1), M, CDUM(1), -1, IERR )
               LWORK_ZUNMBR_PRC_MN = INT( CDUM(1) )
   *
               CALL ZUNMBR( 'P', 'R', 'C', N, N, M, CDUM(1), N, CDUM(1),
        $                   CDUM(1), N, CDUM(1), -1, IERR )
               LWORK_ZUNMBR_PRC_NN = INT( CDUM(1) )
   *
               CALL ZUNMBR( 'Q', 'L', 'N', M, M, M, CDUM(1), M, CDUM(1),
        $                   CDUM(1), M, CDUM(1), -1, IERR )
               LWORK_ZUNMBR_QLN_MM = INT( CDUM(1) )
 *  *
             IF( N.GE.MNTHR1 ) THEN              IF( N.GE.MNTHR1 ) THEN
                IF( WNTQN ) THEN                 IF( WNTQN ) THEN
 *  *
 *                 Path 1t (N much larger than M, JOBZ='N')  *                 Path 1t (N >> M, JOBZ='N')
 *  *
                   MAXWRK = M + M*ILAENV( 1, 'ZGELQF', ' ', M, N, -1,                    MAXWRK = M + LWORK_ZGELQF_MN
      $                     -1 )                    MAXWRK = MAX( MAXWRK, 2*M + LWORK_ZGEBRD_MM )
                   MAXWRK = MAX( MAXWRK, 2*M+2*M*  
      $                     ILAENV( 1, 'ZGEBRD', ' ', M, M, -1, -1 ) )  
                   MINWRK = 3*M                    MINWRK = 3*M
                ELSE IF( WNTQO ) THEN                 ELSE IF( WNTQO ) THEN
 *  *
 *                 Path 2t (N much larger than M, JOBZ='O')  *                 Path 2t (N >> M, JOBZ='O')
 *  *
                   WRKBL = M + M*ILAENV( 1, 'ZGELQF', ' ', M, N, -1, -1 )                    WRKBL = M + LWORK_ZGELQF_MN
                   WRKBL = MAX( WRKBL, M+M*ILAENV( 1, 'ZUNGLQ', ' ', M,                    WRKBL = MAX( WRKBL,   M + LWORK_ZUNGLQ_MN )
      $                    N, M, -1 ) )                    WRKBL = MAX( WRKBL, 2*M + LWORK_ZGEBRD_MM )
                   WRKBL = MAX( WRKBL, 2*M+2*M*                    WRKBL = MAX( WRKBL, 2*M + LWORK_ZUNMBR_QLN_MM )
      $                    ILAENV( 1, 'ZGEBRD', ' ', M, M, -1, -1 ) )                    WRKBL = MAX( WRKBL, 2*M + LWORK_ZUNMBR_PRC_MM )
                   WRKBL = MAX( WRKBL, 2*M+M*  
      $                    ILAENV( 1, 'ZUNMBR', 'PRC', M, M, M, -1 ) )  
                   WRKBL = MAX( WRKBL, 2*M+M*  
      $                    ILAENV( 1, 'ZUNMBR', 'QLN', M, M, M, -1 ) )  
                   MAXWRK = M*N + M*M + WRKBL                    MAXWRK = M*N + M*M + WRKBL
                   MINWRK = 2*M*M + 3*M                    MINWRK = 2*M*M + 3*M
                ELSE IF( WNTQS ) THEN                 ELSE IF( WNTQS ) THEN
 *  *
 *                 Path 3t (N much larger than M, JOBZ='S')  *                 Path 3t (N >> M, JOBZ='S')
 *  *
                   WRKBL = M + M*ILAENV( 1, 'ZGELQF', ' ', M, N, -1, -1 )                    WRKBL = M + LWORK_ZGELQF_MN
                   WRKBL = MAX( WRKBL, M+M*ILAENV( 1, 'ZUNGLQ', ' ', M,                    WRKBL = MAX( WRKBL,   M + LWORK_ZUNGLQ_MN )
      $                    N, M, -1 ) )                    WRKBL = MAX( WRKBL, 2*M + LWORK_ZGEBRD_MM )
                   WRKBL = MAX( WRKBL, 2*M+2*M*                    WRKBL = MAX( WRKBL, 2*M + LWORK_ZUNMBR_QLN_MM )
      $                    ILAENV( 1, 'ZGEBRD', ' ', M, M, -1, -1 ) )                    WRKBL = MAX( WRKBL, 2*M + LWORK_ZUNMBR_PRC_MM )
                   WRKBL = MAX( WRKBL, 2*M+M*  
      $                    ILAENV( 1, 'ZUNMBR', 'PRC', M, M, M, -1 ) )  
                   WRKBL = MAX( WRKBL, 2*M+M*  
      $                    ILAENV( 1, 'ZUNMBR', 'QLN', M, M, M, -1 ) )  
                   MAXWRK = M*M + WRKBL                    MAXWRK = M*M + WRKBL
                   MINWRK = M*M + 3*M                    MINWRK = M*M + 3*M
                ELSE IF( WNTQA ) THEN                 ELSE IF( WNTQA ) THEN
 *  *
 *                 Path 4t (N much larger than M, JOBZ='A')  *                 Path 4t (N >> M, JOBZ='A')
 *  *
                   WRKBL = M + M*ILAENV( 1, 'ZGELQF', ' ', M, N, -1, -1 )                    WRKBL = M + LWORK_ZGELQF_MN
                   WRKBL = MAX( WRKBL, M+N*ILAENV( 1, 'ZUNGLQ', ' ', N,                    WRKBL = MAX( WRKBL,   M + LWORK_ZUNGLQ_NN )
      $                    N, M, -1 ) )                    WRKBL = MAX( WRKBL, 2*M + LWORK_ZGEBRD_MM )
                   WRKBL = MAX( WRKBL, 2*M+2*M*                    WRKBL = MAX( WRKBL, 2*M + LWORK_ZUNMBR_QLN_MM )
      $                    ILAENV( 1, 'ZGEBRD', ' ', M, M, -1, -1 ) )                    WRKBL = MAX( WRKBL, 2*M + LWORK_ZUNMBR_PRC_MM )
                   WRKBL = MAX( WRKBL, 2*M+M*  
      $                    ILAENV( 1, 'ZUNMBR', 'PRC', M, M, M, -1 ) )  
                   WRKBL = MAX( WRKBL, 2*M+M*  
      $                    ILAENV( 1, 'ZUNMBR', 'QLN', M, M, M, -1 ) )  
                   MAXWRK = M*M + WRKBL                    MAXWRK = M*M + WRKBL
                   MINWRK = M*M + 2*M + N                    MINWRK = M*M + MAX( 3*M, M + N )
                END IF                 END IF
             ELSE IF( N.GE.MNTHR2 ) THEN              ELSE IF( N.GE.MNTHR2 ) THEN
 *  *
 *              Path 5t (N much larger than M, but not as much as MNTHR1)  *              Path 5t (N >> M, but not as much as MNTHR1)
 *  *
                MAXWRK = 2*M + ( M+N )*ILAENV( 1, 'ZGEBRD', ' ', M, N,                 MAXWRK = 2*M + LWORK_ZGEBRD_MN
      $                  -1, -1 )  
                MINWRK = 2*M + N                 MINWRK = 2*M + N
                IF( WNTQO ) THEN                 IF( WNTQO ) THEN
                   MAXWRK = MAX( MAXWRK, 2*M+M*  *                 Path 5to (N >> M, JOBZ='O')
      $                     ILAENV( 1, 'ZUNGBR', 'P', M, N, M, -1 ) )                    MAXWRK = MAX( MAXWRK, 2*M + LWORK_ZUNGBR_Q_MM )
                   MAXWRK = MAX( MAXWRK, 2*M+M*                    MAXWRK = MAX( MAXWRK, 2*M + LWORK_ZUNGBR_P_MN )
      $                     ILAENV( 1, 'ZUNGBR', 'Q', M, M, N, -1 ) )  
                   MAXWRK = MAXWRK + M*N                    MAXWRK = MAXWRK + M*N
                   MINWRK = MINWRK + M*M                    MINWRK = MINWRK + M*M
                ELSE IF( WNTQS ) THEN                 ELSE IF( WNTQS ) THEN
                   MAXWRK = MAX( MAXWRK, 2*M+M*  *                 Path 5ts (N >> M, JOBZ='S')
      $                     ILAENV( 1, 'ZUNGBR', 'P', M, N, M, -1 ) )                    MAXWRK = MAX( MAXWRK, 2*M + LWORK_ZUNGBR_Q_MM )
                   MAXWRK = MAX( MAXWRK, 2*M+M*                    MAXWRK = MAX( MAXWRK, 2*M + LWORK_ZUNGBR_P_MN )
      $                     ILAENV( 1, 'ZUNGBR', 'Q', M, M, N, -1 ) )  
                ELSE IF( WNTQA ) THEN                 ELSE IF( WNTQA ) THEN
                   MAXWRK = MAX( MAXWRK, 2*M+N*  *                 Path 5ta (N >> M, JOBZ='A')
      $                     ILAENV( 1, 'ZUNGBR', 'P', N, N, M, -1 ) )                    MAXWRK = MAX( MAXWRK, 2*M + LWORK_ZUNGBR_Q_MM )
                   MAXWRK = MAX( MAXWRK, 2*M+M*                    MAXWRK = MAX( MAXWRK, 2*M + LWORK_ZUNGBR_P_NN )
      $                     ILAENV( 1, 'ZUNGBR', 'Q', M, M, N, -1 ) )  
                END IF                 END IF
             ELSE              ELSE
 *  *
 *              Path 6t (N greater than M, but not much larger)  *              Path 6t (N > M, but not much larger)
 *  *
                MAXWRK = 2*M + ( M+N )*ILAENV( 1, 'ZGEBRD', ' ', M, N,                 MAXWRK = 2*M + LWORK_ZGEBRD_MN
      $                  -1, -1 )  
                MINWRK = 2*M + N                 MINWRK = 2*M + N
                IF( WNTQO ) THEN                 IF( WNTQO ) THEN
                   MAXWRK = MAX( MAXWRK, 2*M+M*  *                 Path 6to (N > M, JOBZ='O')
      $                     ILAENV( 1, 'ZUNMBR', 'PRC', M, N, M, -1 ) )                    MAXWRK = MAX( MAXWRK, 2*M + LWORK_ZUNMBR_QLN_MM )
                   MAXWRK = MAX( MAXWRK, 2*M+M*                    MAXWRK = MAX( MAXWRK, 2*M + LWORK_ZUNMBR_PRC_MN )
      $                     ILAENV( 1, 'ZUNMBR', 'QLN', M, M, N, -1 ) )  
                   MAXWRK = MAXWRK + M*N                    MAXWRK = MAXWRK + M*N
                   MINWRK = MINWRK + M*M                    MINWRK = MINWRK + M*M
                ELSE IF( WNTQS ) THEN                 ELSE IF( WNTQS ) THEN
                   MAXWRK = MAX( MAXWRK, 2*M+M*  *                 Path 6ts (N > M, JOBZ='S')
      $                     ILAENV( 1, 'ZUNGBR', 'PRC', M, N, M, -1 ) )                    MAXWRK = MAX( MAXWRK, 2*M + LWORK_ZUNMBR_QLN_MM )
                   MAXWRK = MAX( MAXWRK, 2*M+M*                    MAXWRK = MAX( MAXWRK, 2*M + LWORK_ZUNMBR_PRC_MN )
      $                     ILAENV( 1, 'ZUNGBR', 'QLN', M, M, N, -1 ) )  
                ELSE IF( WNTQA ) THEN                 ELSE IF( WNTQA ) THEN
                   MAXWRK = MAX( MAXWRK, 2*M+N*  *                 Path 6ta (N > M, JOBZ='A')
      $                     ILAENV( 1, 'ZUNGBR', 'PRC', N, N, M, -1 ) )                    MAXWRK = MAX( MAXWRK, 2*M + LWORK_ZUNMBR_QLN_MM )
                   MAXWRK = MAX( MAXWRK, 2*M+M*                    MAXWRK = MAX( MAXWRK, 2*M + LWORK_ZUNMBR_PRC_NN )
      $                     ILAENV( 1, 'ZUNGBR', 'QLN', M, M, N, -1 ) )  
                END IF                 END IF
             END IF              END IF
          END IF           END IF
Line 459 Line 620
       END IF        END IF
       IF( INFO.EQ.0 ) THEN        IF( INFO.EQ.0 ) THEN
          WORK( 1 ) = MAXWRK           WORK( 1 ) = MAXWRK
          IF( LWORK.LT.MINWRK .AND. LWORK.NE.LQUERV )           IF( LWORK.LT.MINWRK .AND. .NOT. LQUERY ) THEN
      $      INFO = -13              INFO = -12
            END IF
       END IF        END IF
 *  *
 *     Quick returns  
 *  
       IF( INFO.NE.0 ) THEN        IF( INFO.NE.0 ) THEN
          CALL XERBLA( 'ZGESDD', -INFO )           CALL XERBLA( 'ZGESDD', -INFO )
          RETURN           RETURN
         ELSE IF( LQUERY ) THEN
            RETURN
       END IF        END IF
       IF( LWORK.EQ.LQUERV )  *
      $   RETURN  *     Quick return if possible
   *
       IF( M.EQ.0 .OR. N.EQ.0 ) THEN        IF( M.EQ.0 .OR. N.EQ.0 ) THEN
          RETURN           RETURN
       END IF        END IF
Line 503 Line 666
 *  *
             IF( WNTQN ) THEN              IF( WNTQN ) THEN
 *  *
 *              Path 1 (M much larger than N, JOBZ='N')  *              Path 1 (M >> N, JOBZ='N')
 *              No singular vectors to be computed  *              No singular vectors to be computed
 *  *
                ITAU = 1                 ITAU = 1
                NWORK = ITAU + N                 NWORK = ITAU + N
 *  *
 *              Compute A=Q*R  *              Compute A=Q*R
 *              (CWorkspace: need 2*N, prefer N+N*NB)  *              CWorkspace: need   N [tau] + N    [work]
 *              (RWorkspace: need 0)  *              CWorkspace: prefer N [tau] + N*NB [work]
   *              RWorkspace: need   0
 *  *
                CALL ZGEQRF( M, N, A, LDA, WORK( ITAU ), WORK( NWORK ),                 CALL ZGEQRF( M, N, A, LDA, WORK( ITAU ), WORK( NWORK ),
      $                      LWORK-NWORK+1, IERR )       $                      LWORK-NWORK+1, IERR )
Line 526 Line 690
                NWORK = ITAUP + N                 NWORK = ITAUP + N
 *  *
 *              Bidiagonalize R in A  *              Bidiagonalize R in A
 *              (CWorkspace: need 3*N, prefer 2*N+2*N*NB)  *              CWorkspace: need   2*N [tauq, taup] + N      [work]
 *              (RWorkspace: need N)  *              CWorkspace: prefer 2*N [tauq, taup] + 2*N*NB [work]
   *              RWorkspace: need   N [e]
 *  *
                CALL ZGEBRD( N, N, A, LDA, S, RWORK( IE ), WORK( ITAUQ ),                 CALL ZGEBRD( N, N, A, LDA, S, RWORK( IE ), WORK( ITAUQ ),
      $                      WORK( ITAUP ), WORK( NWORK ), LWORK-NWORK+1,       $                      WORK( ITAUP ), WORK( NWORK ), LWORK-NWORK+1,
Line 535 Line 700
                NRWORK = IE + N                 NRWORK = IE + N
 *  *
 *              Perform bidiagonal SVD, compute singular values only  *              Perform bidiagonal SVD, compute singular values only
 *              (CWorkspace: 0)  *              CWorkspace: need   0
 *              (RWorkspace: need BDSPAN)  *              RWorkspace: need   N [e] + BDSPAC
 *  *
                CALL DBDSDC( 'U', 'N', N, S, RWORK( IE ), DUM, 1, DUM, 1,                 CALL DBDSDC( 'U', 'N', N, S, RWORK( IE ), DUM,1,DUM,1,
      $                      DUM, IDUM, RWORK( NRWORK ), IWORK, INFO )       $                      DUM, IDUM, RWORK( NRWORK ), IWORK, INFO )
 *  *
             ELSE IF( WNTQO ) THEN              ELSE IF( WNTQO ) THEN
 *  *
 *              Path 2 (M much larger than N, JOBZ='O')  *              Path 2 (M >> N, JOBZ='O')
 *              N left singular vectors to be overwritten on A and  *              N left singular vectors to be overwritten on A and
 *              N right singular vectors to be computed in VT  *              N right singular vectors to be computed in VT
 *  *
Line 553 Line 718
 *  *
                LDWRKU = N                 LDWRKU = N
                IR = IU + LDWRKU*N                 IR = IU + LDWRKU*N
                IF( LWORK.GE.M*N+N*N+3*N ) THEN                 IF( LWORK .GE. M*N + N*N + 3*N ) THEN
 *  *
 *                 WORK(IR) is M by N  *                 WORK(IR) is M by N
 *  *
                   LDWRKR = M                    LDWRKR = M
                ELSE                 ELSE
                   LDWRKR = ( LWORK-N*N-3*N ) / N                    LDWRKR = ( LWORK - N*N - 3*N ) / N
                END IF                 END IF
                ITAU = IR + LDWRKR*N                 ITAU = IR + LDWRKR*N
                NWORK = ITAU + N                 NWORK = ITAU + N
 *  *
 *              Compute A=Q*R  *              Compute A=Q*R
 *              (CWorkspace: need N*N+2*N, prefer M*N+N+N*NB)  *              CWorkspace: need   N*N [U] + N*N [R] + N [tau] + N    [work]
 *              (RWorkspace: 0)  *              CWorkspace: prefer N*N [U] + N*N [R] + N [tau] + N*NB [work]
   *              RWorkspace: need   0
 *  *
                CALL ZGEQRF( M, N, A, LDA, WORK( ITAU ), WORK( NWORK ),                 CALL ZGEQRF( M, N, A, LDA, WORK( ITAU ), WORK( NWORK ),
      $                      LWORK-NWORK+1, IERR )       $                      LWORK-NWORK+1, IERR )
Line 578 Line 744
      $                      LDWRKR )       $                      LDWRKR )
 *  *
 *              Generate Q in A  *              Generate Q in A
 *              (CWorkspace: need 2*N, prefer N+N*NB)  *              CWorkspace: need   N*N [U] + N*N [R] + N [tau] + N    [work]
 *              (RWorkspace: 0)  *              CWorkspace: prefer N*N [U] + N*N [R] + N [tau] + N*NB [work]
   *              RWorkspace: need   0
 *  *
                CALL ZUNGQR( M, N, N, A, LDA, WORK( ITAU ),                 CALL ZUNGQR( M, N, N, A, LDA, WORK( ITAU ),
      $                      WORK( NWORK ), LWORK-NWORK+1, IERR )       $                      WORK( NWORK ), LWORK-NWORK+1, IERR )
Line 589 Line 756
                NWORK = ITAUP + N                 NWORK = ITAUP + N
 *  *
 *              Bidiagonalize R in WORK(IR)  *              Bidiagonalize R in WORK(IR)
 *              (CWorkspace: need N*N+3*N, prefer M*N+2*N+2*N*NB)  *              CWorkspace: need   N*N [U] + N*N [R] + 2*N [tauq, taup] + N      [work]
 *              (RWorkspace: need N)  *              CWorkspace: prefer N*N [U] + N*N [R] + 2*N [tauq, taup] + 2*N*NB [work]
   *              RWorkspace: need   N [e]
 *  *
                CALL ZGEBRD( N, N, WORK( IR ), LDWRKR, S, RWORK( IE ),                 CALL ZGEBRD( N, N, WORK( IR ), LDWRKR, S, RWORK( IE ),
      $                      WORK( ITAUQ ), WORK( ITAUP ), WORK( NWORK ),       $                      WORK( ITAUQ ), WORK( ITAUP ), WORK( NWORK ),
Line 599 Line 767
 *              Perform bidiagonal SVD, computing left singular vectors  *              Perform bidiagonal SVD, computing left singular vectors
 *              of R in WORK(IRU) and computing right singular vectors  *              of R in WORK(IRU) and computing right singular vectors
 *              of R in WORK(IRVT)  *              of R in WORK(IRVT)
 *              (CWorkspace: need 0)  *              CWorkspace: need   0
 *              (RWorkspace: need BDSPAC)  *              RWorkspace: need   N [e] + N*N [RU] + N*N [RVT] + BDSPAC
 *  *
                IRU = IE + N                 IRU = IE + N
                IRVT = IRU + N*N                 IRVT = IRU + N*N
Line 611 Line 779
 *  *
 *              Copy real matrix RWORK(IRU) to complex matrix WORK(IU)  *              Copy real matrix RWORK(IRU) to complex matrix WORK(IU)
 *              Overwrite WORK(IU) by the left singular vectors of R  *              Overwrite WORK(IU) by the left singular vectors of R
 *              (CWorkspace: need 2*N*N+3*N, prefer M*N+N*N+2*N+N*NB)  *              CWorkspace: need   N*N [U] + N*N [R] + 2*N [tauq, taup] + N    [work]
 *              (RWorkspace: 0)  *              CWorkspace: prefer N*N [U] + N*N [R] + 2*N [tauq, taup] + N*NB [work]
   *              RWorkspace: need   0
 *  *
                CALL ZLACP2( 'F', N, N, RWORK( IRU ), N, WORK( IU ),                 CALL ZLACP2( 'F', N, N, RWORK( IRU ), N, WORK( IU ),
      $                      LDWRKU )       $                      LDWRKU )
Line 622 Line 791
 *  *
 *              Copy real matrix RWORK(IRVT) to complex matrix VT  *              Copy real matrix RWORK(IRVT) to complex matrix VT
 *              Overwrite VT by the right singular vectors of R  *              Overwrite VT by the right singular vectors of R
 *              (CWorkspace: need N*N+3*N, prefer M*N+2*N+N*NB)  *              CWorkspace: need   N*N [U] + N*N [R] + 2*N [tauq, taup] + N    [work]
 *              (RWorkspace: 0)  *              CWorkspace: prefer N*N [U] + N*N [R] + 2*N [tauq, taup] + N*NB [work]
   *              RWorkspace: need   0
 *  *
                CALL ZLACP2( 'F', N, N, RWORK( IRVT ), N, VT, LDVT )                 CALL ZLACP2( 'F', N, N, RWORK( IRVT ), N, VT, LDVT )
                CALL ZUNMBR( 'P', 'R', 'C', N, N, N, WORK( IR ), LDWRKR,                 CALL ZUNMBR( 'P', 'R', 'C', N, N, N, WORK( IR ), LDWRKR,
Line 632 Line 802
 *  *
 *              Multiply Q in A by left singular vectors of R in  *              Multiply Q in A by left singular vectors of R in
 *              WORK(IU), storing result in WORK(IR) and copying to A  *              WORK(IU), storing result in WORK(IR) and copying to A
 *              (CWorkspace: need 2*N*N, prefer N*N+M*N)  *              CWorkspace: need   N*N [U] + N*N [R]
 *              (RWorkspace: 0)  *              CWorkspace: prefer N*N [U] + M*N [R]
   *              RWorkspace: need   0
 *  *
                DO 10 I = 1, M, LDWRKR                 DO 10 I = 1, M, LDWRKR
                   CHUNK = MIN( M-I+1, LDWRKR )                    CHUNK = MIN( M-I+1, LDWRKR )
Line 646 Line 817
 *  *
             ELSE IF( WNTQS ) THEN              ELSE IF( WNTQS ) THEN
 *  *
 *              Path 3 (M much larger than N, JOBZ='S')  *              Path 3 (M >> N, JOBZ='S')
 *              N left singular vectors to be computed in U and  *              N left singular vectors to be computed in U and
 *              N right singular vectors to be computed in VT  *              N right singular vectors to be computed in VT
 *  *
Line 659 Line 830
                NWORK = ITAU + N                 NWORK = ITAU + N
 *  *
 *              Compute A=Q*R  *              Compute A=Q*R
 *              (CWorkspace: need N*N+2*N, prefer N*N+N+N*NB)  *              CWorkspace: need   N*N [R] + N [tau] + N    [work]
 *              (RWorkspace: 0)  *              CWorkspace: prefer N*N [R] + N [tau] + N*NB [work]
   *              RWorkspace: need   0
 *  *
                CALL ZGEQRF( M, N, A, LDA, WORK( ITAU ), WORK( NWORK ),                 CALL ZGEQRF( M, N, A, LDA, WORK( ITAU ), WORK( NWORK ),
      $                      LWORK-NWORK+1, IERR )       $                      LWORK-NWORK+1, IERR )
Line 672 Line 844
      $                      LDWRKR )       $                      LDWRKR )
 *  *
 *              Generate Q in A  *              Generate Q in A
 *              (CWorkspace: need 2*N, prefer N+N*NB)  *              CWorkspace: need   N*N [R] + N [tau] + N    [work]
 *              (RWorkspace: 0)  *              CWorkspace: prefer N*N [R] + N [tau] + N*NB [work]
   *              RWorkspace: need   0
 *  *
                CALL ZUNGQR( M, N, N, A, LDA, WORK( ITAU ),                 CALL ZUNGQR( M, N, N, A, LDA, WORK( ITAU ),
      $                      WORK( NWORK ), LWORK-NWORK+1, IERR )       $                      WORK( NWORK ), LWORK-NWORK+1, IERR )
Line 683 Line 856
                NWORK = ITAUP + N                 NWORK = ITAUP + N
 *  *
 *              Bidiagonalize R in WORK(IR)  *              Bidiagonalize R in WORK(IR)
 *              (CWorkspace: need N*N+3*N, prefer N*N+2*N+2*N*NB)  *              CWorkspace: need   N*N [R] + 2*N [tauq, taup] + N      [work]
 *              (RWorkspace: need N)  *              CWorkspace: prefer N*N [R] + 2*N [tauq, taup] + 2*N*NB [work]
   *              RWorkspace: need   N [e]
 *  *
                CALL ZGEBRD( N, N, WORK( IR ), LDWRKR, S, RWORK( IE ),                 CALL ZGEBRD( N, N, WORK( IR ), LDWRKR, S, RWORK( IE ),
      $                      WORK( ITAUQ ), WORK( ITAUP ), WORK( NWORK ),       $                      WORK( ITAUQ ), WORK( ITAUP ), WORK( NWORK ),
Line 693 Line 867
 *              Perform bidiagonal SVD, computing left singular vectors  *              Perform bidiagonal SVD, computing left singular vectors
 *              of bidiagonal matrix in RWORK(IRU) and computing right  *              of bidiagonal matrix in RWORK(IRU) and computing right
 *              singular vectors of bidiagonal matrix in RWORK(IRVT)  *              singular vectors of bidiagonal matrix in RWORK(IRVT)
 *              (CWorkspace: need 0)  *              CWorkspace: need   0
 *              (RWorkspace: need BDSPAC)  *              RWorkspace: need   N [e] + N*N [RU] + N*N [RVT] + BDSPAC
 *  *
                IRU = IE + N                 IRU = IE + N
                IRVT = IRU + N*N                 IRVT = IRU + N*N
Line 705 Line 879
 *  *
 *              Copy real matrix RWORK(IRU) to complex matrix U  *              Copy real matrix RWORK(IRU) to complex matrix U
 *              Overwrite U by left singular vectors of R  *              Overwrite U by left singular vectors of R
 *              (CWorkspace: need N*N+3*N, prefer N*N+2*N+N*NB)  *              CWorkspace: need   N*N [R] + 2*N [tauq, taup] + N    [work]
 *              (RWorkspace: 0)  *              CWorkspace: prefer N*N [R] + 2*N [tauq, taup] + N*NB [work]
   *              RWorkspace: need   0
 *  *
                CALL ZLACP2( 'F', N, N, RWORK( IRU ), N, U, LDU )                 CALL ZLACP2( 'F', N, N, RWORK( IRU ), N, U, LDU )
                CALL ZUNMBR( 'Q', 'L', 'N', N, N, N, WORK( IR ), LDWRKR,                 CALL ZUNMBR( 'Q', 'L', 'N', N, N, N, WORK( IR ), LDWRKR,
Line 715 Line 890
 *  *
 *              Copy real matrix RWORK(IRVT) to complex matrix VT  *              Copy real matrix RWORK(IRVT) to complex matrix VT
 *              Overwrite VT by right singular vectors of R  *              Overwrite VT by right singular vectors of R
 *              (CWorkspace: need N*N+3*N, prefer N*N+2*N+N*NB)  *              CWorkspace: need   N*N [R] + 2*N [tauq, taup] + N    [work]
 *              (RWorkspace: 0)  *              CWorkspace: prefer N*N [R] + 2*N [tauq, taup] + N*NB [work]
   *              RWorkspace: need   0
 *  *
                CALL ZLACP2( 'F', N, N, RWORK( IRVT ), N, VT, LDVT )                 CALL ZLACP2( 'F', N, N, RWORK( IRVT ), N, VT, LDVT )
                CALL ZUNMBR( 'P', 'R', 'C', N, N, N, WORK( IR ), LDWRKR,                 CALL ZUNMBR( 'P', 'R', 'C', N, N, N, WORK( IR ), LDWRKR,
Line 725 Line 901
 *  *
 *              Multiply Q in A by left singular vectors of R in  *              Multiply Q in A by left singular vectors of R in
 *              WORK(IR), storing result in U  *              WORK(IR), storing result in U
 *              (CWorkspace: need N*N)  *              CWorkspace: need   N*N [R]
 *              (RWorkspace: 0)  *              RWorkspace: need   0
 *  *
                CALL ZLACPY( 'F', N, N, U, LDU, WORK( IR ), LDWRKR )                 CALL ZLACPY( 'F', N, N, U, LDU, WORK( IR ), LDWRKR )
                CALL ZGEMM( 'N', 'N', M, N, N, CONE, A, LDA, WORK( IR ),                 CALL ZGEMM( 'N', 'N', M, N, N, CONE, A, LDA, WORK( IR ),
Line 734 Line 910
 *  *
             ELSE IF( WNTQA ) THEN              ELSE IF( WNTQA ) THEN
 *  *
 *              Path 4 (M much larger than N, JOBZ='A')  *              Path 4 (M >> N, JOBZ='A')
 *              M left singular vectors to be computed in U and  *              M left singular vectors to be computed in U and
 *              N right singular vectors to be computed in VT  *              N right singular vectors to be computed in VT
 *  *
Line 747 Line 923
                NWORK = ITAU + N                 NWORK = ITAU + N
 *  *
 *              Compute A=Q*R, copying result to U  *              Compute A=Q*R, copying result to U
 *              (CWorkspace: need 2*N, prefer N+N*NB)  *              CWorkspace: need   N*N [U] + N [tau] + N    [work]
 *              (RWorkspace: 0)  *              CWorkspace: prefer N*N [U] + N [tau] + N*NB [work]
   *              RWorkspace: need   0
 *  *
                CALL ZGEQRF( M, N, A, LDA, WORK( ITAU ), WORK( NWORK ),                 CALL ZGEQRF( M, N, A, LDA, WORK( ITAU ), WORK( NWORK ),
      $                      LWORK-NWORK+1, IERR )       $                      LWORK-NWORK+1, IERR )
                CALL ZLACPY( 'L', M, N, A, LDA, U, LDU )                 CALL ZLACPY( 'L', M, N, A, LDA, U, LDU )
 *  *
 *              Generate Q in U  *              Generate Q in U
 *              (CWorkspace: need N+M, prefer N+M*NB)  *              CWorkspace: need   N*N [U] + N [tau] + M    [work]
 *              (RWorkspace: 0)  *              CWorkspace: prefer N*N [U] + N [tau] + M*NB [work]
   *              RWorkspace: need   0
 *  *
                CALL ZUNGQR( M, M, N, U, LDU, WORK( ITAU ),                 CALL ZUNGQR( M, M, N, U, LDU, WORK( ITAU ),
      $                      WORK( NWORK ), LWORK-NWORK+1, IERR )       $                      WORK( NWORK ), LWORK-NWORK+1, IERR )
Line 771 Line 949
                NWORK = ITAUP + N                 NWORK = ITAUP + N
 *  *
 *              Bidiagonalize R in A  *              Bidiagonalize R in A
 *              (CWorkspace: need 3*N, prefer 2*N+2*N*NB)  *              CWorkspace: need   N*N [U] + 2*N [tauq, taup] + N      [work]
 *              (RWorkspace: need N)  *              CWorkspace: prefer N*N [U] + 2*N [tauq, taup] + 2*N*NB [work]
   *              RWorkspace: need   N [e]
 *  *
                CALL ZGEBRD( N, N, A, LDA, S, RWORK( IE ), WORK( ITAUQ ),                 CALL ZGEBRD( N, N, A, LDA, S, RWORK( IE ), WORK( ITAUQ ),
      $                      WORK( ITAUP ), WORK( NWORK ), LWORK-NWORK+1,       $                      WORK( ITAUP ), WORK( NWORK ), LWORK-NWORK+1,
Line 784 Line 963
 *              Perform bidiagonal SVD, computing left singular vectors  *              Perform bidiagonal SVD, computing left singular vectors
 *              of bidiagonal matrix in RWORK(IRU) and computing right  *              of bidiagonal matrix in RWORK(IRU) and computing right
 *              singular vectors of bidiagonal matrix in RWORK(IRVT)  *              singular vectors of bidiagonal matrix in RWORK(IRVT)
 *              (CWorkspace: need 0)  *              CWorkspace: need   0
 *              (RWorkspace: need BDSPAC)  *              RWorkspace: need   N [e] + N*N [RU] + N*N [RVT] + BDSPAC
 *  *
                CALL DBDSDC( 'U', 'I', N, S, RWORK( IE ), RWORK( IRU ),                 CALL DBDSDC( 'U', 'I', N, S, RWORK( IE ), RWORK( IRU ),
      $                      N, RWORK( IRVT ), N, DUM, IDUM,       $                      N, RWORK( IRVT ), N, DUM, IDUM,
Line 793 Line 972
 *  *
 *              Copy real matrix RWORK(IRU) to complex matrix WORK(IU)  *              Copy real matrix RWORK(IRU) to complex matrix WORK(IU)
 *              Overwrite WORK(IU) by left singular vectors of R  *              Overwrite WORK(IU) by left singular vectors of R
 *              (CWorkspace: need N*N+3*N, prefer N*N+2*N+N*NB)  *              CWorkspace: need   N*N [U] + 2*N [tauq, taup] + N    [work]
 *              (RWorkspace: 0)  *              CWorkspace: prefer N*N [U] + 2*N [tauq, taup] + N*NB [work]
   *              RWorkspace: need   0
 *  *
                CALL ZLACP2( 'F', N, N, RWORK( IRU ), N, WORK( IU ),                 CALL ZLACP2( 'F', N, N, RWORK( IRU ), N, WORK( IU ),
      $                      LDWRKU )       $                      LDWRKU )
Line 804 Line 984
 *  *
 *              Copy real matrix RWORK(IRVT) to complex matrix VT  *              Copy real matrix RWORK(IRVT) to complex matrix VT
 *              Overwrite VT by right singular vectors of R  *              Overwrite VT by right singular vectors of R
 *              (CWorkspace: need 3*N, prefer 2*N+N*NB)  *              CWorkspace: need   N*N [U] + 2*N [tauq, taup] + N    [work]
 *              (RWorkspace: 0)  *              CWorkspace: prefer N*N [U] + 2*N [tauq, taup] + N*NB [work]
   *              RWorkspace: need   0
 *  *
                CALL ZLACP2( 'F', N, N, RWORK( IRVT ), N, VT, LDVT )                 CALL ZLACP2( 'F', N, N, RWORK( IRVT ), N, VT, LDVT )
                CALL ZUNMBR( 'P', 'R', 'C', N, N, N, A, LDA,                 CALL ZUNMBR( 'P', 'R', 'C', N, N, N, A, LDA,
Line 814 Line 995
 *  *
 *              Multiply Q in U by left singular vectors of R in  *              Multiply Q in U by left singular vectors of R in
 *              WORK(IU), storing result in A  *              WORK(IU), storing result in A
 *              (CWorkspace: need N*N)  *              CWorkspace: need   N*N [U]
 *              (RWorkspace: 0)  *              RWorkspace: need   0
 *  *
                CALL ZGEMM( 'N', 'N', M, N, N, CONE, U, LDU, WORK( IU ),                 CALL ZGEMM( 'N', 'N', M, N, N, CONE, U, LDU, WORK( IU ),
      $                     LDWRKU, CZERO, A, LDA )       $                     LDWRKU, CZERO, A, LDA )
Line 830 Line 1011
 *  *
 *           MNTHR2 <= M < MNTHR1  *           MNTHR2 <= M < MNTHR1
 *  *
 *           Path 5 (M much larger than N, but not as much as MNTHR1)  *           Path 5 (M >> N, but not as much as MNTHR1)
 *           Reduce to bidiagonal form without QR decomposition, use  *           Reduce to bidiagonal form without QR decomposition, use
 *           ZUNGBR and matrix multiplication to compute singular vectors  *           ZUNGBR and matrix multiplication to compute singular vectors
 *  *
Line 841 Line 1022
             NWORK = ITAUP + N              NWORK = ITAUP + N
 *  *
 *           Bidiagonalize A  *           Bidiagonalize A
 *           (CWorkspace: need 2*N+M, prefer 2*N+(M+N)*NB)  *           CWorkspace: need   2*N [tauq, taup] + M        [work]
 *           (RWorkspace: need N)  *           CWorkspace: prefer 2*N [tauq, taup] + (M+N)*NB [work]
   *           RWorkspace: need   N [e]
 *  *
             CALL ZGEBRD( M, N, A, LDA, S, RWORK( IE ), WORK( ITAUQ ),              CALL ZGEBRD( M, N, A, LDA, S, RWORK( IE ), WORK( ITAUQ ),
      $                   WORK( ITAUP ), WORK( NWORK ), LWORK-NWORK+1,       $                   WORK( ITAUP ), WORK( NWORK ), LWORK-NWORK+1,
      $                   IERR )       $                   IERR )
             IF( WNTQN ) THEN              IF( WNTQN ) THEN
 *  *
   *              Path 5n (M >> N, JOBZ='N')
 *              Compute singular values only  *              Compute singular values only
 *              (Cworkspace: 0)  *              CWorkspace: need   0
 *              (Rworkspace: need BDSPAN)  *              RWorkspace: need   N [e] + BDSPAC
 *  *
                CALL DBDSDC( 'U', 'N', N, S, RWORK( IE ), DUM, 1, DUM, 1,                 CALL DBDSDC( 'U', 'N', N, S, RWORK( IE ), DUM, 1,DUM,1,
      $                      DUM, IDUM, RWORK( NRWORK ), IWORK, INFO )       $                      DUM, IDUM, RWORK( NRWORK ), IWORK, INFO )
             ELSE IF( WNTQO ) THEN              ELSE IF( WNTQO ) THEN
                IU = NWORK                 IU = NWORK
Line 861 Line 1044
                IRVT = IRU + N*N                 IRVT = IRU + N*N
                NRWORK = IRVT + N*N                 NRWORK = IRVT + N*N
 *  *
   *              Path 5o (M >> N, JOBZ='O')
 *              Copy A to VT, generate P**H  *              Copy A to VT, generate P**H
 *              (Cworkspace: need 2*N, prefer N+N*NB)  *              CWorkspace: need   2*N [tauq, taup] + N    [work]
 *              (Rworkspace: 0)  *              CWorkspace: prefer 2*N [tauq, taup] + N*NB [work]
   *              RWorkspace: need   0
 *  *
                CALL ZLACPY( 'U', N, N, A, LDA, VT, LDVT )                 CALL ZLACPY( 'U', N, N, A, LDA, VT, LDVT )
                CALL ZUNGBR( 'P', N, N, N, VT, LDVT, WORK( ITAUP ),                 CALL ZUNGBR( 'P', N, N, N, VT, LDVT, WORK( ITAUP ),
      $                      WORK( NWORK ), LWORK-NWORK+1, IERR )       $                      WORK( NWORK ), LWORK-NWORK+1, IERR )
 *  *
 *              Generate Q in A  *              Generate Q in A
 *              (CWorkspace: need 2*N, prefer N+N*NB)  *              CWorkspace: need   2*N [tauq, taup] + N    [work]
 *              (RWorkspace: 0)  *              CWorkspace: prefer 2*N [tauq, taup] + N*NB [work]
   *              RWorkspace: need   0
 *  *
                CALL ZUNGBR( 'Q', M, N, N, A, LDA, WORK( ITAUQ ),                 CALL ZUNGBR( 'Q', M, N, N, A, LDA, WORK( ITAUQ ),
      $                      WORK( NWORK ), LWORK-NWORK+1, IERR )       $                      WORK( NWORK ), LWORK-NWORK+1, IERR )
 *  *
                IF( LWORK.GE.M*N+3*N ) THEN                 IF( LWORK .GE. M*N + 3*N ) THEN
 *  *
 *                 WORK( IU ) is M by N  *                 WORK( IU ) is M by N
 *  *
Line 885 Line 1071
 *  *
 *                 WORK(IU) is LDWRKU by N  *                 WORK(IU) is LDWRKU by N
 *  *
                   LDWRKU = ( LWORK-3*N ) / N                    LDWRKU = ( LWORK - 3*N ) / N
                END IF                 END IF
                NWORK = IU + LDWRKU*N                 NWORK = IU + LDWRKU*N
 *  *
 *              Perform bidiagonal SVD, computing left singular vectors  *              Perform bidiagonal SVD, computing left singular vectors
 *              of bidiagonal matrix in RWORK(IRU) and computing right  *              of bidiagonal matrix in RWORK(IRU) and computing right
 *              singular vectors of bidiagonal matrix in RWORK(IRVT)  *              singular vectors of bidiagonal matrix in RWORK(IRVT)
 *              (CWorkspace: need 0)  *              CWorkspace: need   0
 *              (RWorkspace: need BDSPAC)  *              RWorkspace: need   N [e] + N*N [RU] + N*N [RVT] + BDSPAC
 *  *
                CALL DBDSDC( 'U', 'I', N, S, RWORK( IE ), RWORK( IRU ),                 CALL DBDSDC( 'U', 'I', N, S, RWORK( IE ), RWORK( IRU ),
      $                      N, RWORK( IRVT ), N, DUM, IDUM,       $                      N, RWORK( IRVT ), N, DUM, IDUM,
Line 901 Line 1087
 *  *
 *              Multiply real matrix RWORK(IRVT) by P**H in VT,  *              Multiply real matrix RWORK(IRVT) by P**H in VT,
 *              storing the result in WORK(IU), copying to VT  *              storing the result in WORK(IU), copying to VT
 *              (Cworkspace: need 0)  *              CWorkspace: need   2*N [tauq, taup] + N*N [U]
 *              (Rworkspace: need 3*N*N)  *              RWorkspace: need   N [e] + N*N [RU] + N*N [RVT] + 2*N*N [rwork]
 *  *
                CALL ZLARCM( N, N, RWORK( IRVT ), N, VT, LDVT,                 CALL ZLARCM( N, N, RWORK( IRVT ), N, VT, LDVT,
      $                      WORK( IU ), LDWRKU, RWORK( NRWORK ) )       $                      WORK( IU ), LDWRKU, RWORK( NRWORK ) )
Line 910 Line 1096
 *  *
 *              Multiply Q in A by real matrix RWORK(IRU), storing the  *              Multiply Q in A by real matrix RWORK(IRU), storing the
 *              result in WORK(IU), copying to A  *              result in WORK(IU), copying to A
 *              (CWorkspace: need N*N, prefer M*N)  *              CWorkspace: need   2*N [tauq, taup] + N*N [U]
 *              (Rworkspace: need 3*N*N, prefer N*N+2*M*N)  *              CWorkspace: prefer 2*N [tauq, taup] + M*N [U]
   *              RWorkspace: need   N [e] + N*N [RU] + 2*N*N [rwork]
   *              RWorkspace: prefer N [e] + N*N [RU] + 2*M*N [rwork] < N + 5*N*N since M < 2*N here
 *  *
                NRWORK = IRVT                 NRWORK = IRVT
                DO 20 I = 1, M, LDWRKU                 DO 20 I = 1, M, LDWRKU
Line 924 Line 1112
 *  *
             ELSE IF( WNTQS ) THEN              ELSE IF( WNTQS ) THEN
 *  *
   *              Path 5s (M >> N, JOBZ='S')
 *              Copy A to VT, generate P**H  *              Copy A to VT, generate P**H
 *              (Cworkspace: need 2*N, prefer N+N*NB)  *              CWorkspace: need   2*N [tauq, taup] + N    [work]
 *              (Rworkspace: 0)  *              CWorkspace: prefer 2*N [tauq, taup] + N*NB [work]
   *              RWorkspace: need   0
 *  *
                CALL ZLACPY( 'U', N, N, A, LDA, VT, LDVT )                 CALL ZLACPY( 'U', N, N, A, LDA, VT, LDVT )
                CALL ZUNGBR( 'P', N, N, N, VT, LDVT, WORK( ITAUP ),                 CALL ZUNGBR( 'P', N, N, N, VT, LDVT, WORK( ITAUP ),
      $                      WORK( NWORK ), LWORK-NWORK+1, IERR )       $                      WORK( NWORK ), LWORK-NWORK+1, IERR )
 *  *
 *              Copy A to U, generate Q  *              Copy A to U, generate Q
 *              (Cworkspace: need 2*N, prefer N+N*NB)  *              CWorkspace: need   2*N [tauq, taup] + N    [work]
 *              (Rworkspace: 0)  *              CWorkspace: prefer 2*N [tauq, taup] + N*NB [work]
   *              RWorkspace: need   0
 *  *
                CALL ZLACPY( 'L', M, N, A, LDA, U, LDU )                 CALL ZLACPY( 'L', M, N, A, LDA, U, LDU )
                CALL ZUNGBR( 'Q', M, N, N, U, LDU, WORK( ITAUQ ),                 CALL ZUNGBR( 'Q', M, N, N, U, LDU, WORK( ITAUQ ),
Line 943 Line 1134
 *              Perform bidiagonal SVD, computing left singular vectors  *              Perform bidiagonal SVD, computing left singular vectors
 *              of bidiagonal matrix in RWORK(IRU) and computing right  *              of bidiagonal matrix in RWORK(IRU) and computing right
 *              singular vectors of bidiagonal matrix in RWORK(IRVT)  *              singular vectors of bidiagonal matrix in RWORK(IRVT)
 *              (CWorkspace: need 0)  *              CWorkspace: need   0
 *              (RWorkspace: need BDSPAC)  *              RWorkspace: need   N [e] + N*N [RU] + N*N [RVT] + BDSPAC
 *  *
                IRU = NRWORK                 IRU = NRWORK
                IRVT = IRU + N*N                 IRVT = IRU + N*N
Line 955 Line 1146
 *  *
 *              Multiply real matrix RWORK(IRVT) by P**H in VT,  *              Multiply real matrix RWORK(IRVT) by P**H in VT,
 *              storing the result in A, copying to VT  *              storing the result in A, copying to VT
 *              (Cworkspace: need 0)  *              CWorkspace: need   0
 *              (Rworkspace: need 3*N*N)  *              RWorkspace: need   N [e] + N*N [RU] + N*N [RVT] + 2*N*N [rwork]
 *  *
                CALL ZLARCM( N, N, RWORK( IRVT ), N, VT, LDVT, A, LDA,                 CALL ZLARCM( N, N, RWORK( IRVT ), N, VT, LDVT, A, LDA,
      $                      RWORK( NRWORK ) )       $                      RWORK( NRWORK ) )
Line 964 Line 1155
 *  *
 *              Multiply Q in U by real matrix RWORK(IRU), storing the  *              Multiply Q in U by real matrix RWORK(IRU), storing the
 *              result in A, copying to U  *              result in A, copying to U
 *              (CWorkspace: need 0)  *              CWorkspace: need   0
 *              (Rworkspace: need N*N+2*M*N)  *              RWorkspace: need   N [e] + N*N [RU] + 2*M*N [rwork] < N + 5*N*N since M < 2*N here
 *  *
                NRWORK = IRVT                 NRWORK = IRVT
                CALL ZLACRM( M, N, U, LDU, RWORK( IRU ), N, A, LDA,                 CALL ZLACRM( M, N, U, LDU, RWORK( IRU ), N, A, LDA,
Line 973 Line 1164
                CALL ZLACPY( 'F', M, N, A, LDA, U, LDU )                 CALL ZLACPY( 'F', M, N, A, LDA, U, LDU )
             ELSE              ELSE
 *  *
   *              Path 5a (M >> N, JOBZ='A')
 *              Copy A to VT, generate P**H  *              Copy A to VT, generate P**H
 *              (Cworkspace: need 2*N, prefer N+N*NB)  *              CWorkspace: need   2*N [tauq, taup] + N    [work]
 *              (Rworkspace: 0)  *              CWorkspace: prefer 2*N [tauq, taup] + N*NB [work]
   *              RWorkspace: need   0
 *  *
                CALL ZLACPY( 'U', N, N, A, LDA, VT, LDVT )                 CALL ZLACPY( 'U', N, N, A, LDA, VT, LDVT )
                CALL ZUNGBR( 'P', N, N, N, VT, LDVT, WORK( ITAUP ),                 CALL ZUNGBR( 'P', N, N, N, VT, LDVT, WORK( ITAUP ),
      $                      WORK( NWORK ), LWORK-NWORK+1, IERR )       $                      WORK( NWORK ), LWORK-NWORK+1, IERR )
 *  *
 *              Copy A to U, generate Q  *              Copy A to U, generate Q
 *              (Cworkspace: need 2*N, prefer N+N*NB)  *              CWorkspace: need   2*N [tauq, taup] + M    [work]
 *              (Rworkspace: 0)  *              CWorkspace: prefer 2*N [tauq, taup] + M*NB [work]
   *              RWorkspace: need   0
 *  *
                CALL ZLACPY( 'L', M, N, A, LDA, U, LDU )                 CALL ZLACPY( 'L', M, N, A, LDA, U, LDU )
                CALL ZUNGBR( 'Q', M, M, N, U, LDU, WORK( ITAUQ ),                 CALL ZUNGBR( 'Q', M, M, N, U, LDU, WORK( ITAUQ ),
Line 992 Line 1186
 *              Perform bidiagonal SVD, computing left singular vectors  *              Perform bidiagonal SVD, computing left singular vectors
 *              of bidiagonal matrix in RWORK(IRU) and computing right  *              of bidiagonal matrix in RWORK(IRU) and computing right
 *              singular vectors of bidiagonal matrix in RWORK(IRVT)  *              singular vectors of bidiagonal matrix in RWORK(IRVT)
 *              (CWorkspace: need 0)  *              CWorkspace: need   0
 *              (RWorkspace: need BDSPAC)  *              RWorkspace: need   N [e] + N*N [RU] + N*N [RVT] + BDSPAC
 *  *
                IRU = NRWORK                 IRU = NRWORK
                IRVT = IRU + N*N                 IRVT = IRU + N*N
Line 1004 Line 1198
 *  *
 *              Multiply real matrix RWORK(IRVT) by P**H in VT,  *              Multiply real matrix RWORK(IRVT) by P**H in VT,
 *              storing the result in A, copying to VT  *              storing the result in A, copying to VT
 *              (Cworkspace: need 0)  *              CWorkspace: need   0
 *              (Rworkspace: need 3*N*N)  *              RWorkspace: need   N [e] + N*N [RU] + N*N [RVT] + 2*N*N [rwork]
 *  *
                CALL ZLARCM( N, N, RWORK( IRVT ), N, VT, LDVT, A, LDA,                 CALL ZLARCM( N, N, RWORK( IRVT ), N, VT, LDVT, A, LDA,
      $                      RWORK( NRWORK ) )       $                      RWORK( NRWORK ) )
Line 1013 Line 1207
 *  *
 *              Multiply Q in U by real matrix RWORK(IRU), storing the  *              Multiply Q in U by real matrix RWORK(IRU), storing the
 *              result in A, copying to U  *              result in A, copying to U
 *              (CWorkspace: 0)  *              CWorkspace: need   0
 *              (Rworkspace: need 3*N*N)  *              RWorkspace: need   N [e] + N*N [RU] + 2*M*N [rwork] < N + 5*N*N since M < 2*N here
 *  *
                NRWORK = IRVT                 NRWORK = IRVT
                CALL ZLACRM( M, N, U, LDU, RWORK( IRU ), N, A, LDA,                 CALL ZLACRM( M, N, U, LDU, RWORK( IRU ), N, A, LDA,
Line 1026 Line 1220
 *  *
 *           M .LT. MNTHR2  *           M .LT. MNTHR2
 *  *
 *           Path 6 (M at least N, but not much larger)  *           Path 6 (M >= N, but not much larger)
 *           Reduce to bidiagonal form without QR decomposition  *           Reduce to bidiagonal form without QR decomposition
 *           Use ZUNMBR to compute singular vectors  *           Use ZUNMBR to compute singular vectors
 *  *
Line 1037 Line 1231
             NWORK = ITAUP + N              NWORK = ITAUP + N
 *  *
 *           Bidiagonalize A  *           Bidiagonalize A
 *           (CWorkspace: need 2*N+M, prefer 2*N+(M+N)*NB)  *           CWorkspace: need   2*N [tauq, taup] + M        [work]
 *           (RWorkspace: need N)  *           CWorkspace: prefer 2*N [tauq, taup] + (M+N)*NB [work]
   *           RWorkspace: need   N [e]
 *  *
             CALL ZGEBRD( M, N, A, LDA, S, RWORK( IE ), WORK( ITAUQ ),              CALL ZGEBRD( M, N, A, LDA, S, RWORK( IE ), WORK( ITAUQ ),
      $                   WORK( ITAUP ), WORK( NWORK ), LWORK-NWORK+1,       $                   WORK( ITAUP ), WORK( NWORK ), LWORK-NWORK+1,
      $                   IERR )       $                   IERR )
             IF( WNTQN ) THEN              IF( WNTQN ) THEN
 *  *
   *              Path 6n (M >= N, JOBZ='N')
 *              Compute singular values only  *              Compute singular values only
 *              (Cworkspace: 0)  *              CWorkspace: need   0
 *              (Rworkspace: need BDSPAN)  *              RWorkspace: need   N [e] + BDSPAC
 *  *
                CALL DBDSDC( 'U', 'N', N, S, RWORK( IE ), DUM, 1, DUM, 1,                 CALL DBDSDC( 'U', 'N', N, S, RWORK( IE ), DUM,1,DUM,1,
      $                      DUM, IDUM, RWORK( NRWORK ), IWORK, INFO )       $                      DUM, IDUM, RWORK( NRWORK ), IWORK, INFO )
             ELSE IF( WNTQO ) THEN              ELSE IF( WNTQO ) THEN
                IU = NWORK                 IU = NWORK
                IRU = NRWORK                 IRU = NRWORK
                IRVT = IRU + N*N                 IRVT = IRU + N*N
                NRWORK = IRVT + N*N                 NRWORK = IRVT + N*N
                IF( LWORK.GE.M*N+3*N ) THEN                 IF( LWORK .GE. M*N + 3*N ) THEN
 *  *
 *                 WORK( IU ) is M by N  *                 WORK( IU ) is M by N
 *  *
Line 1065 Line 1261
 *  *
 *                 WORK( IU ) is LDWRKU by N  *                 WORK( IU ) is LDWRKU by N
 *  *
                   LDWRKU = ( LWORK-3*N ) / N                    LDWRKU = ( LWORK - 3*N ) / N
                END IF                 END IF
                NWORK = IU + LDWRKU*N                 NWORK = IU + LDWRKU*N
 *  *
   *              Path 6o (M >= N, JOBZ='O')
 *              Perform bidiagonal SVD, computing left singular vectors  *              Perform bidiagonal SVD, computing left singular vectors
 *              of bidiagonal matrix in RWORK(IRU) and computing right  *              of bidiagonal matrix in RWORK(IRU) and computing right
 *              singular vectors of bidiagonal matrix in RWORK(IRVT)  *              singular vectors of bidiagonal matrix in RWORK(IRVT)
 *              (CWorkspace: need 0)  *              CWorkspace: need   0
 *              (RWorkspace: need BDSPAC)  *              RWorkspace: need   N [e] + N*N [RU] + N*N [RVT] + BDSPAC
 *  *
                CALL DBDSDC( 'U', 'I', N, S, RWORK( IE ), RWORK( IRU ),                 CALL DBDSDC( 'U', 'I', N, S, RWORK( IE ), RWORK( IRU ),
      $                      N, RWORK( IRVT ), N, DUM, IDUM,       $                      N, RWORK( IRVT ), N, DUM, IDUM,
Line 1081 Line 1278
 *  *
 *              Copy real matrix RWORK(IRVT) to complex matrix VT  *              Copy real matrix RWORK(IRVT) to complex matrix VT
 *              Overwrite VT by right singular vectors of A  *              Overwrite VT by right singular vectors of A
 *              (Cworkspace: need 2*N, prefer N+N*NB)  *              CWorkspace: need   2*N [tauq, taup] + N*N [U] + N    [work]
 *              (Rworkspace: need 0)  *              CWorkspace: prefer 2*N [tauq, taup] + N*N [U] + N*NB [work]
   *              RWorkspace: need   N [e] + N*N [RU] + N*N [RVT]
 *  *
                CALL ZLACP2( 'F', N, N, RWORK( IRVT ), N, VT, LDVT )                 CALL ZLACP2( 'F', N, N, RWORK( IRVT ), N, VT, LDVT )
                CALL ZUNMBR( 'P', 'R', 'C', N, N, N, A, LDA,                 CALL ZUNMBR( 'P', 'R', 'C', N, N, N, A, LDA,
      $                      WORK( ITAUP ), VT, LDVT, WORK( NWORK ),       $                      WORK( ITAUP ), VT, LDVT, WORK( NWORK ),
      $                      LWORK-NWORK+1, IERR )       $                      LWORK-NWORK+1, IERR )
 *  *
                IF( LWORK.GE.M*N+3*N ) THEN                 IF( LWORK .GE. M*N + 3*N ) THEN
 *  *
 *              Copy real matrix RWORK(IRU) to complex matrix WORK(IU)  *                 Path 6o-fast
 *              Overwrite WORK(IU) by left singular vectors of A, copying  *                 Copy real matrix RWORK(IRU) to complex matrix WORK(IU)
 *              to A  *                 Overwrite WORK(IU) by left singular vectors of A, copying
 *              (Cworkspace: need M*N+2*N, prefer M*N+N+N*NB)  *                 to A
 *              (Rworkspace: need 0)  *                 CWorkspace: need   2*N [tauq, taup] + M*N [U] + N    [work]
   *                 CWorkspace: prefer 2*N [tauq, taup] + M*N [U] + N*NB [work]
   *                 RWorkspace: need   N [e] + N*N [RU]
 *  *
                   CALL ZLASET( 'F', M, N, CZERO, CZERO, WORK( IU ),                    CALL ZLASET( 'F', M, N, CZERO, CZERO, WORK( IU ),
      $                         LDWRKU )       $                         LDWRKU )
Line 1107 Line 1307
                   CALL ZLACPY( 'F', M, N, WORK( IU ), LDWRKU, A, LDA )                    CALL ZLACPY( 'F', M, N, WORK( IU ), LDWRKU, A, LDA )
                ELSE                 ELSE
 *  *
   *                 Path 6o-slow
 *                 Generate Q in A  *                 Generate Q in A
 *                 (Cworkspace: need 2*N, prefer N+N*NB)  *                 CWorkspace: need   2*N [tauq, taup] + N*N [U] + N    [work]
 *                 (Rworkspace: need 0)  *                 CWorkspace: prefer 2*N [tauq, taup] + N*N [U] + N*NB [work]
   *                 RWorkspace: need   0
 *  *
                   CALL ZUNGBR( 'Q', M, N, N, A, LDA, WORK( ITAUQ ),                    CALL ZUNGBR( 'Q', M, N, N, A, LDA, WORK( ITAUQ ),
      $                         WORK( NWORK ), LWORK-NWORK+1, IERR )       $                         WORK( NWORK ), LWORK-NWORK+1, IERR )
 *  *
 *                 Multiply Q in A by real matrix RWORK(IRU), storing the  *                 Multiply Q in A by real matrix RWORK(IRU), storing the
 *                 result in WORK(IU), copying to A  *                 result in WORK(IU), copying to A
 *                 (CWorkspace: need N*N, prefer M*N)  *                 CWorkspace: need   2*N [tauq, taup] + N*N [U]
 *                 (Rworkspace: need 3*N*N, prefer N*N+2*M*N)  *                 CWorkspace: prefer 2*N [tauq, taup] + M*N [U]
   *                 RWorkspace: need   N [e] + N*N [RU] + 2*N*N [rwork]
   *                 RWorkspace: prefer N [e] + N*N [RU] + 2*M*N [rwork] < N + 5*N*N since M < 2*N here
 *  *
                   NRWORK = IRVT                    NRWORK = IRVT
                   DO 30 I = 1, M, LDWRKU                    DO 30 I = 1, M, LDWRKU
Line 1132 Line 1336
 *  *
             ELSE IF( WNTQS ) THEN              ELSE IF( WNTQS ) THEN
 *  *
   *              Path 6s (M >= N, JOBZ='S')
 *              Perform bidiagonal SVD, computing left singular vectors  *              Perform bidiagonal SVD, computing left singular vectors
 *              of bidiagonal matrix in RWORK(IRU) and computing right  *              of bidiagonal matrix in RWORK(IRU) and computing right
 *              singular vectors of bidiagonal matrix in RWORK(IRVT)  *              singular vectors of bidiagonal matrix in RWORK(IRVT)
 *              (CWorkspace: need 0)  *              CWorkspace: need   0
 *              (RWorkspace: need BDSPAC)  *              RWorkspace: need   N [e] + N*N [RU] + N*N [RVT] + BDSPAC
 *  *
                IRU = NRWORK                 IRU = NRWORK
                IRVT = IRU + N*N                 IRVT = IRU + N*N
Line 1147 Line 1352
 *  *
 *              Copy real matrix RWORK(IRU) to complex matrix U  *              Copy real matrix RWORK(IRU) to complex matrix U
 *              Overwrite U by left singular vectors of A  *              Overwrite U by left singular vectors of A
 *              (CWorkspace: need 3*N, prefer 2*N+N*NB)  *              CWorkspace: need   2*N [tauq, taup] + N    [work]
 *              (RWorkspace: 0)  *              CWorkspace: prefer 2*N [tauq, taup] + N*NB [work]
   *              RWorkspace: need   N [e] + N*N [RU] + N*N [RVT]
 *  *
                CALL ZLASET( 'F', M, N, CZERO, CZERO, U, LDU )                 CALL ZLASET( 'F', M, N, CZERO, CZERO, U, LDU )
                CALL ZLACP2( 'F', N, N, RWORK( IRU ), N, U, LDU )                 CALL ZLACP2( 'F', N, N, RWORK( IRU ), N, U, LDU )
Line 1158 Line 1364
 *  *
 *              Copy real matrix RWORK(IRVT) to complex matrix VT  *              Copy real matrix RWORK(IRVT) to complex matrix VT
 *              Overwrite VT by right singular vectors of A  *              Overwrite VT by right singular vectors of A
 *              (CWorkspace: need 3*N, prefer 2*N+N*NB)  *              CWorkspace: need   2*N [tauq, taup] + N    [work]
 *              (RWorkspace: 0)  *              CWorkspace: prefer 2*N [tauq, taup] + N*NB [work]
   *              RWorkspace: need   N [e] + N*N [RU] + N*N [RVT]
 *  *
                CALL ZLACP2( 'F', N, N, RWORK( IRVT ), N, VT, LDVT )                 CALL ZLACP2( 'F', N, N, RWORK( IRVT ), N, VT, LDVT )
                CALL ZUNMBR( 'P', 'R', 'C', N, N, N, A, LDA,                 CALL ZUNMBR( 'P', 'R', 'C', N, N, N, A, LDA,
Line 1167 Line 1374
      $                      LWORK-NWORK+1, IERR )       $                      LWORK-NWORK+1, IERR )
             ELSE              ELSE
 *  *
   *              Path 6a (M >= N, JOBZ='A')
 *              Perform bidiagonal SVD, computing left singular vectors  *              Perform bidiagonal SVD, computing left singular vectors
 *              of bidiagonal matrix in RWORK(IRU) and computing right  *              of bidiagonal matrix in RWORK(IRU) and computing right
 *              singular vectors of bidiagonal matrix in RWORK(IRVT)  *              singular vectors of bidiagonal matrix in RWORK(IRVT)
 *              (CWorkspace: need 0)  *              CWorkspace: need   0
 *              (RWorkspace: need BDSPAC)  *              RWorkspace: need   N [e] + N*N [RU] + N*N [RVT] + BDSPAC
 *  *
                IRU = NRWORK                 IRU = NRWORK
                IRVT = IRU + N*N                 IRVT = IRU + N*N
Line 1190 Line 1398
 *  *
 *              Copy real matrix RWORK(IRU) to complex matrix U  *              Copy real matrix RWORK(IRU) to complex matrix U
 *              Overwrite U by left singular vectors of A  *              Overwrite U by left singular vectors of A
 *              (CWorkspace: need 2*N+M, prefer 2*N+M*NB)  *              CWorkspace: need   2*N [tauq, taup] + M    [work]
 *              (RWorkspace: 0)  *              CWorkspace: prefer 2*N [tauq, taup] + M*NB [work]
   *              RWorkspace: need   N [e] + N*N [RU] + N*N [RVT]
 *  *
                CALL ZLACP2( 'F', N, N, RWORK( IRU ), N, U, LDU )                 CALL ZLACP2( 'F', N, N, RWORK( IRU ), N, U, LDU )
                CALL ZUNMBR( 'Q', 'L', 'N', M, M, N, A, LDA,                 CALL ZUNMBR( 'Q', 'L', 'N', M, M, N, A, LDA,
Line 1200 Line 1409
 *  *
 *              Copy real matrix RWORK(IRVT) to complex matrix VT  *              Copy real matrix RWORK(IRVT) to complex matrix VT
 *              Overwrite VT by right singular vectors of A  *              Overwrite VT by right singular vectors of A
 *              (CWorkspace: need 3*N, prefer 2*N+N*NB)  *              CWorkspace: need   2*N [tauq, taup] + N    [work]
 *              (RWorkspace: 0)  *              CWorkspace: prefer 2*N [tauq, taup] + N*NB [work]
   *              RWorkspace: need   N [e] + N*N [RU] + N*N [RVT]
 *  *
                CALL ZLACP2( 'F', N, N, RWORK( IRVT ), N, VT, LDVT )                 CALL ZLACP2( 'F', N, N, RWORK( IRVT ), N, VT, LDVT )
                CALL ZUNMBR( 'P', 'R', 'C', N, N, N, A, LDA,                 CALL ZUNMBR( 'P', 'R', 'C', N, N, N, A, LDA,
Line 1221 Line 1431
 *  *
             IF( WNTQN ) THEN              IF( WNTQN ) THEN
 *  *
 *              Path 1t (N much larger than M, JOBZ='N')  *              Path 1t (N >> M, JOBZ='N')
 *              No singular vectors to be computed  *              No singular vectors to be computed
 *  *
                ITAU = 1                 ITAU = 1
                NWORK = ITAU + M                 NWORK = ITAU + M
 *  *
 *              Compute A=L*Q  *              Compute A=L*Q
 *              (CWorkspace: need 2*M, prefer M+M*NB)  *              CWorkspace: need   M [tau] + M    [work]
 *              (RWorkspace: 0)  *              CWorkspace: prefer M [tau] + M*NB [work]
   *              RWorkspace: need   0
 *  *
                CALL ZGELQF( M, N, A, LDA, WORK( ITAU ), WORK( NWORK ),                 CALL ZGELQF( M, N, A, LDA, WORK( ITAU ), WORK( NWORK ),
      $                      LWORK-NWORK+1, IERR )       $                      LWORK-NWORK+1, IERR )
Line 1244 Line 1455
                NWORK = ITAUP + M                 NWORK = ITAUP + M
 *  *
 *              Bidiagonalize L in A  *              Bidiagonalize L in A
 *              (CWorkspace: need 3*M, prefer 2*M+2*M*NB)  *              CWorkspace: need   2*M [tauq, taup] + M      [work]
 *              (RWorkspace: need M)  *              CWorkspace: prefer 2*M [tauq, taup] + 2*M*NB [work]
   *              RWorkspace: need   M [e]
 *  *
                CALL ZGEBRD( M, M, A, LDA, S, RWORK( IE ), WORK( ITAUQ ),                 CALL ZGEBRD( M, M, A, LDA, S, RWORK( IE ), WORK( ITAUQ ),
      $                      WORK( ITAUP ), WORK( NWORK ), LWORK-NWORK+1,       $                      WORK( ITAUP ), WORK( NWORK ), LWORK-NWORK+1,
Line 1253 Line 1465
                NRWORK = IE + M                 NRWORK = IE + M
 *  *
 *              Perform bidiagonal SVD, compute singular values only  *              Perform bidiagonal SVD, compute singular values only
 *              (CWorkspace: 0)  *              CWorkspace: need   0
 *              (RWorkspace: need BDSPAN)  *              RWorkspace: need   M [e] + BDSPAC
 *  *
                CALL DBDSDC( 'U', 'N', M, S, RWORK( IE ), DUM, 1, DUM, 1,                 CALL DBDSDC( 'U', 'N', M, S, RWORK( IE ), DUM,1,DUM,1,
      $                      DUM, IDUM, RWORK( NRWORK ), IWORK, INFO )       $                      DUM, IDUM, RWORK( NRWORK ), IWORK, INFO )
 *  *
             ELSE IF( WNTQO ) THEN              ELSE IF( WNTQO ) THEN
 *  *
 *              Path 2t (N much larger than M, JOBZ='O')  *              Path 2t (N >> M, JOBZ='O')
 *              M right singular vectors to be overwritten on A and  *              M right singular vectors to be overwritten on A and
 *              M left singular vectors to be computed in U  *              M left singular vectors to be computed in U
 *  *
Line 1271 Line 1483
 *              WORK(IVT) is M by M  *              WORK(IVT) is M by M
 *  *
                IL = IVT + LDWKVT*M                 IL = IVT + LDWKVT*M
                IF( LWORK.GE.M*N+M*M+3*M ) THEN                 IF( LWORK .GE. M*N + M*M + 3*M ) THEN
 *  *
 *                 WORK(IL) M by N  *                 WORK(IL) M by N
 *  *
Line 1282 Line 1494
 *                 WORK(IL) is M by CHUNK  *                 WORK(IL) is M by CHUNK
 *  *
                   LDWRKL = M                    LDWRKL = M
                   CHUNK = ( LWORK-M*M-3*M ) / M                    CHUNK = ( LWORK - M*M - 3*M ) / M
                END IF                 END IF
                ITAU = IL + LDWRKL*CHUNK                 ITAU = IL + LDWRKL*CHUNK
                NWORK = ITAU + M                 NWORK = ITAU + M
 *  *
 *              Compute A=L*Q  *              Compute A=L*Q
 *              (CWorkspace: need 2*M, prefer M+M*NB)  *              CWorkspace: need   M*M [VT] + M*M [L] + M [tau] + M    [work]
 *              (RWorkspace: 0)  *              CWorkspace: prefer M*M [VT] + M*M [L] + M [tau] + M*NB [work]
   *              RWorkspace: need   0
 *  *
                CALL ZGELQF( M, N, A, LDA, WORK( ITAU ), WORK( NWORK ),                 CALL ZGELQF( M, N, A, LDA, WORK( ITAU ), WORK( NWORK ),
      $                      LWORK-NWORK+1, IERR )       $                      LWORK-NWORK+1, IERR )
Line 1301 Line 1514
      $                      WORK( IL+LDWRKL ), LDWRKL )       $                      WORK( IL+LDWRKL ), LDWRKL )
 *  *
 *              Generate Q in A  *              Generate Q in A
 *              (CWorkspace: need M*M+2*M, prefer M*M+M+M*NB)  *              CWorkspace: need   M*M [VT] + M*M [L] + M [tau] + M    [work]
 *              (RWorkspace: 0)  *              CWorkspace: prefer M*M [VT] + M*M [L] + M [tau] + M*NB [work]
   *              RWorkspace: need   0
 *  *
                CALL ZUNGLQ( M, N, M, A, LDA, WORK( ITAU ),                 CALL ZUNGLQ( M, N, M, A, LDA, WORK( ITAU ),
      $                      WORK( NWORK ), LWORK-NWORK+1, IERR )       $                      WORK( NWORK ), LWORK-NWORK+1, IERR )
Line 1312 Line 1526
                NWORK = ITAUP + M                 NWORK = ITAUP + M
 *  *
 *              Bidiagonalize L in WORK(IL)  *              Bidiagonalize L in WORK(IL)
 *              (CWorkspace: need M*M+3*M, prefer M*M+2*M+2*M*NB)  *              CWorkspace: need   M*M [VT] + M*M [L] + 2*M [tauq, taup] + M      [work]
 *              (RWorkspace: need M)  *              CWorkspace: prefer M*M [VT] + M*M [L] + 2*M [tauq, taup] + 2*M*NB [work]
   *              RWorkspace: need   M [e]
 *  *
                CALL ZGEBRD( M, M, WORK( IL ), LDWRKL, S, RWORK( IE ),                 CALL ZGEBRD( M, M, WORK( IL ), LDWRKL, S, RWORK( IE ),
      $                      WORK( ITAUQ ), WORK( ITAUP ), WORK( NWORK ),       $                      WORK( ITAUQ ), WORK( ITAUP ), WORK( NWORK ),
Line 1322 Line 1537
 *              Perform bidiagonal SVD, computing left singular vectors  *              Perform bidiagonal SVD, computing left singular vectors
 *              of bidiagonal matrix in RWORK(IRU) and computing right  *              of bidiagonal matrix in RWORK(IRU) and computing right
 *              singular vectors of bidiagonal matrix in RWORK(IRVT)  *              singular vectors of bidiagonal matrix in RWORK(IRVT)
 *              (CWorkspace: need 0)  *              CWorkspace: need   0
 *              (RWorkspace: need BDSPAC)  *              RWorkspace: need   M [e] + M*M [RU] + M*M [RVT] + BDSPAC
 *  *
                IRU = IE + M                 IRU = IE + M
                IRVT = IRU + M*M                 IRVT = IRU + M*M
Line 1334 Line 1549
 *  *
 *              Copy real matrix RWORK(IRU) to complex matrix WORK(IU)  *              Copy real matrix RWORK(IRU) to complex matrix WORK(IU)
 *              Overwrite WORK(IU) by the left singular vectors of L  *              Overwrite WORK(IU) by the left singular vectors of L
 *              (CWorkspace: need N*N+3*N, prefer M*N+2*N+N*NB)  *              CWorkspace: need   M*M [VT] + M*M [L] + 2*M [tauq, taup] + M    [work]
 *              (RWorkspace: 0)  *              CWorkspace: prefer M*M [VT] + M*M [L] + 2*M [tauq, taup] + M*NB [work]
   *              RWorkspace: need   0
 *  *
                CALL ZLACP2( 'F', M, M, RWORK( IRU ), M, U, LDU )                 CALL ZLACP2( 'F', M, M, RWORK( IRU ), M, U, LDU )
                CALL ZUNMBR( 'Q', 'L', 'N', M, M, M, WORK( IL ), LDWRKL,                 CALL ZUNMBR( 'Q', 'L', 'N', M, M, M, WORK( IL ), LDWRKL,
Line 1344 Line 1560
 *  *
 *              Copy real matrix RWORK(IRVT) to complex matrix WORK(IVT)  *              Copy real matrix RWORK(IRVT) to complex matrix WORK(IVT)
 *              Overwrite WORK(IVT) by the right singular vectors of L  *              Overwrite WORK(IVT) by the right singular vectors of L
 *              (CWorkspace: need N*N+3*N, prefer M*N+2*N+N*NB)  *              CWorkspace: need   M*M [VT] + M*M [L] + 2*M [tauq, taup] + M    [work]
 *              (RWorkspace: 0)  *              CWorkspace: prefer M*M [VT] + M*M [L] + 2*M [tauq, taup] + M*NB [work]
   *              RWorkspace: need   0
 *  *
                CALL ZLACP2( 'F', M, M, RWORK( IRVT ), M, WORK( IVT ),                 CALL ZLACP2( 'F', M, M, RWORK( IRVT ), M, WORK( IVT ),
      $                      LDWKVT )       $                      LDWKVT )
Line 1355 Line 1572
 *  *
 *              Multiply right singular vectors of L in WORK(IL) by Q  *              Multiply right singular vectors of L in WORK(IL) by Q
 *              in A, storing result in WORK(IL) and copying to A  *              in A, storing result in WORK(IL) and copying to A
 *              (CWorkspace: need 2*M*M, prefer M*M+M*N))  *              CWorkspace: need   M*M [VT] + M*M [L]
 *              (RWorkspace: 0)  *              CWorkspace: prefer M*M [VT] + M*N [L]
   *              RWorkspace: need   0
 *  *
                DO 40 I = 1, N, CHUNK                 DO 40 I = 1, N, CHUNK
                   BLK = MIN( N-I+1, CHUNK )                    BLK = MIN( N-I+1, CHUNK )
Line 1369 Line 1587
 *  *
             ELSE IF( WNTQS ) THEN              ELSE IF( WNTQS ) THEN
 *  *
 *             Path 3t (N much larger than M, JOBZ='S')  *              Path 3t (N >> M, JOBZ='S')
 *             M right singular vectors to be computed in VT and  *              M right singular vectors to be computed in VT and
 *             M left singular vectors to be computed in U  *              M left singular vectors to be computed in U
 *  *
                IL = 1                 IL = 1
 *  *
Line 1382 Line 1600
                NWORK = ITAU + M                 NWORK = ITAU + M
 *  *
 *              Compute A=L*Q  *              Compute A=L*Q
 *              (CWorkspace: need 2*M, prefer M+M*NB)  *              CWorkspace: need   M*M [L] + M [tau] + M    [work]
 *              (RWorkspace: 0)  *              CWorkspace: prefer M*M [L] + M [tau] + M*NB [work]
   *              RWorkspace: need   0
 *  *
                CALL ZGELQF( M, N, A, LDA, WORK( ITAU ), WORK( NWORK ),                 CALL ZGELQF( M, N, A, LDA, WORK( ITAU ), WORK( NWORK ),
      $                      LWORK-NWORK+1, IERR )       $                      LWORK-NWORK+1, IERR )
Line 1395 Line 1614
      $                      WORK( IL+LDWRKL ), LDWRKL )       $                      WORK( IL+LDWRKL ), LDWRKL )
 *  *
 *              Generate Q in A  *              Generate Q in A
 *              (CWorkspace: need M*M+2*M, prefer M*M+M+M*NB)  *              CWorkspace: need   M*M [L] + M [tau] + M    [work]
 *              (RWorkspace: 0)  *              CWorkspace: prefer M*M [L] + M [tau] + M*NB [work]
   *              RWorkspace: need   0
 *  *
                CALL ZUNGLQ( M, N, M, A, LDA, WORK( ITAU ),                 CALL ZUNGLQ( M, N, M, A, LDA, WORK( ITAU ),
      $                      WORK( NWORK ), LWORK-NWORK+1, IERR )       $                      WORK( NWORK ), LWORK-NWORK+1, IERR )
Line 1406 Line 1626
                NWORK = ITAUP + M                 NWORK = ITAUP + M
 *  *
 *              Bidiagonalize L in WORK(IL)  *              Bidiagonalize L in WORK(IL)
 *              (CWorkspace: need M*M+3*M, prefer M*M+2*M+2*M*NB)  *              CWorkspace: need   M*M [L] + 2*M [tauq, taup] + M      [work]
 *              (RWorkspace: need M)  *              CWorkspace: prefer M*M [L] + 2*M [tauq, taup] + 2*M*NB [work]
   *              RWorkspace: need   M [e]
 *  *
                CALL ZGEBRD( M, M, WORK( IL ), LDWRKL, S, RWORK( IE ),                 CALL ZGEBRD( M, M, WORK( IL ), LDWRKL, S, RWORK( IE ),
      $                      WORK( ITAUQ ), WORK( ITAUP ), WORK( NWORK ),       $                      WORK( ITAUQ ), WORK( ITAUP ), WORK( NWORK ),
Line 1416 Line 1637
 *              Perform bidiagonal SVD, computing left singular vectors  *              Perform bidiagonal SVD, computing left singular vectors
 *              of bidiagonal matrix in RWORK(IRU) and computing right  *              of bidiagonal matrix in RWORK(IRU) and computing right
 *              singular vectors of bidiagonal matrix in RWORK(IRVT)  *              singular vectors of bidiagonal matrix in RWORK(IRVT)
 *              (CWorkspace: need 0)  *              CWorkspace: need   0
 *              (RWorkspace: need BDSPAC)  *              RWorkspace: need   M [e] + M*M [RU] + M*M [RVT] + BDSPAC
 *  *
                IRU = IE + M                 IRU = IE + M
                IRVT = IRU + M*M                 IRVT = IRU + M*M
Line 1428 Line 1649
 *  *
 *              Copy real matrix RWORK(IRU) to complex matrix U  *              Copy real matrix RWORK(IRU) to complex matrix U
 *              Overwrite U by left singular vectors of L  *              Overwrite U by left singular vectors of L
 *              (CWorkspace: need M*M+3*M, prefer M*M+2*M+M*NB)  *              CWorkspace: need   M*M [L] + 2*M [tauq, taup] + M    [work]
 *              (RWorkspace: 0)  *              CWorkspace: prefer M*M [L] + 2*M [tauq, taup] + M*NB [work]
   *              RWorkspace: need   0
 *  *
                CALL ZLACP2( 'F', M, M, RWORK( IRU ), M, U, LDU )                 CALL ZLACP2( 'F', M, M, RWORK( IRU ), M, U, LDU )
                CALL ZUNMBR( 'Q', 'L', 'N', M, M, M, WORK( IL ), LDWRKL,                 CALL ZUNMBR( 'Q', 'L', 'N', M, M, M, WORK( IL ), LDWRKL,
Line 1438 Line 1660
 *  *
 *              Copy real matrix RWORK(IRVT) to complex matrix VT  *              Copy real matrix RWORK(IRVT) to complex matrix VT
 *              Overwrite VT by left singular vectors of L  *              Overwrite VT by left singular vectors of L
 *              (CWorkspace: need M*M+3*M, prefer M*M+2*M+M*NB)  *              CWorkspace: need   M*M [L] + 2*M [tauq, taup] + M    [work]
 *              (RWorkspace: 0)  *              CWorkspace: prefer M*M [L] + 2*M [tauq, taup] + M*NB [work]
   *              RWorkspace: need   0
 *  *
                CALL ZLACP2( 'F', M, M, RWORK( IRVT ), M, VT, LDVT )                 CALL ZLACP2( 'F', M, M, RWORK( IRVT ), M, VT, LDVT )
                CALL ZUNMBR( 'P', 'R', 'C', M, M, M, WORK( IL ), LDWRKL,                 CALL ZUNMBR( 'P', 'R', 'C', M, M, M, WORK( IL ), LDWRKL,
Line 1448 Line 1671
 *  *
 *              Copy VT to WORK(IL), multiply right singular vectors of L  *              Copy VT to WORK(IL), multiply right singular vectors of L
 *              in WORK(IL) by Q in A, storing result in VT  *              in WORK(IL) by Q in A, storing result in VT
 *              (CWorkspace: need M*M)  *              CWorkspace: need   M*M [L]
 *              (RWorkspace: 0)  *              RWorkspace: need   0
 *  *
                CALL ZLACPY( 'F', M, M, VT, LDVT, WORK( IL ), LDWRKL )                 CALL ZLACPY( 'F', M, M, VT, LDVT, WORK( IL ), LDWRKL )
                CALL ZGEMM( 'N', 'N', M, N, M, CONE, WORK( IL ), LDWRKL,                 CALL ZGEMM( 'N', 'N', M, N, M, CONE, WORK( IL ), LDWRKL,
Line 1457 Line 1680
 *  *
             ELSE IF( WNTQA ) THEN              ELSE IF( WNTQA ) THEN
 *  *
 *              Path 9t (N much larger than M, JOBZ='A')  *              Path 4t (N >> M, JOBZ='A')
 *              N right singular vectors to be computed in VT and  *              N right singular vectors to be computed in VT and
 *              M left singular vectors to be computed in U  *              M left singular vectors to be computed in U
 *  *
Line 1470 Line 1693
                NWORK = ITAU + M                 NWORK = ITAU + M
 *  *
 *              Compute A=L*Q, copying result to VT  *              Compute A=L*Q, copying result to VT
 *              (CWorkspace: need 2*M, prefer M+M*NB)  *              CWorkspace: need   M*M [VT] + M [tau] + M    [work]
 *              (RWorkspace: 0)  *              CWorkspace: prefer M*M [VT] + M [tau] + M*NB [work]
   *              RWorkspace: need   0
 *  *
                CALL ZGELQF( M, N, A, LDA, WORK( ITAU ), WORK( NWORK ),                 CALL ZGELQF( M, N, A, LDA, WORK( ITAU ), WORK( NWORK ),
      $                      LWORK-NWORK+1, IERR )       $                      LWORK-NWORK+1, IERR )
                CALL ZLACPY( 'U', M, N, A, LDA, VT, LDVT )                 CALL ZLACPY( 'U', M, N, A, LDA, VT, LDVT )
 *  *
 *              Generate Q in VT  *              Generate Q in VT
 *              (CWorkspace: need M+N, prefer M+N*NB)  *              CWorkspace: need   M*M [VT] + M [tau] + N    [work]
 *              (RWorkspace: 0)  *              CWorkspace: prefer M*M [VT] + M [tau] + N*NB [work]
   *              RWorkspace: need   0
 *  *
                CALL ZUNGLQ( N, N, M, VT, LDVT, WORK( ITAU ),                 CALL ZUNGLQ( N, N, M, VT, LDVT, WORK( ITAU ),
      $                      WORK( NWORK ), LWORK-NWORK+1, IERR )       $                      WORK( NWORK ), LWORK-NWORK+1, IERR )
Line 1494 Line 1719
                NWORK = ITAUP + M                 NWORK = ITAUP + M
 *  *
 *              Bidiagonalize L in A  *              Bidiagonalize L in A
 *              (CWorkspace: need M*M+3*M, prefer M*M+2*M+2*M*NB)  *              CWorkspace: need   M*M [VT] + 2*M [tauq, taup] + M      [work]
 *              (RWorkspace: need M)  *              CWorkspace: prefer M*M [VT] + 2*M [tauq, taup] + 2*M*NB [work]
   *              RWorkspace: need   M [e]
 *  *
                CALL ZGEBRD( M, M, A, LDA, S, RWORK( IE ), WORK( ITAUQ ),                 CALL ZGEBRD( M, M, A, LDA, S, RWORK( IE ), WORK( ITAUQ ),
      $                      WORK( ITAUP ), WORK( NWORK ), LWORK-NWORK+1,       $                      WORK( ITAUP ), WORK( NWORK ), LWORK-NWORK+1,
Line 1504 Line 1730
 *              Perform bidiagonal SVD, computing left singular vectors  *              Perform bidiagonal SVD, computing left singular vectors
 *              of bidiagonal matrix in RWORK(IRU) and computing right  *              of bidiagonal matrix in RWORK(IRU) and computing right
 *              singular vectors of bidiagonal matrix in RWORK(IRVT)  *              singular vectors of bidiagonal matrix in RWORK(IRVT)
 *              (CWorkspace: need 0)  *              CWorkspace: need   0
 *              (RWorkspace: need BDSPAC)  *              RWorkspace: need   M [e] + M*M [RU] + M*M [RVT] + BDSPAC
 *  *
                IRU = IE + M                 IRU = IE + M
                IRVT = IRU + M*M                 IRVT = IRU + M*M
Line 1516 Line 1742
 *  *
 *              Copy real matrix RWORK(IRU) to complex matrix U  *              Copy real matrix RWORK(IRU) to complex matrix U
 *              Overwrite U by left singular vectors of L  *              Overwrite U by left singular vectors of L
 *              (CWorkspace: need 3*M, prefer 2*M+M*NB)  *              CWorkspace: need   M*M [VT] + 2*M [tauq, taup] + M    [work]
 *              (RWorkspace: 0)  *              CWorkspace: prefer M*M [VT] + 2*M [tauq, taup] + M*NB [work]
   *              RWorkspace: need   0
 *  *
                CALL ZLACP2( 'F', M, M, RWORK( IRU ), M, U, LDU )                 CALL ZLACP2( 'F', M, M, RWORK( IRU ), M, U, LDU )
                CALL ZUNMBR( 'Q', 'L', 'N', M, M, M, A, LDA,                 CALL ZUNMBR( 'Q', 'L', 'N', M, M, M, A, LDA,
Line 1526 Line 1753
 *  *
 *              Copy real matrix RWORK(IRVT) to complex matrix WORK(IVT)  *              Copy real matrix RWORK(IRVT) to complex matrix WORK(IVT)
 *              Overwrite WORK(IVT) by right singular vectors of L  *              Overwrite WORK(IVT) by right singular vectors of L
 *              (CWorkspace: need M*M+3*M, prefer M*M+2*M+M*NB)  *              CWorkspace: need   M*M [VT] + 2*M [tauq, taup] + M    [work]
 *              (RWorkspace: 0)  *              CWorkspace: prefer M*M [VT] + 2*M [tauq, taup] + M*NB [work]
   *              RWorkspace: need   0
 *  *
                CALL ZLACP2( 'F', M, M, RWORK( IRVT ), M, WORK( IVT ),                 CALL ZLACP2( 'F', M, M, RWORK( IRVT ), M, WORK( IVT ),
      $                      LDWKVT )       $                      LDWKVT )
Line 1537 Line 1765
 *  *
 *              Multiply right singular vectors of L in WORK(IVT) by  *              Multiply right singular vectors of L in WORK(IVT) by
 *              Q in VT, storing result in A  *              Q in VT, storing result in A
 *              (CWorkspace: need M*M)  *              CWorkspace: need   M*M [VT]
 *              (RWorkspace: 0)  *              RWorkspace: need   0
 *  *
                CALL ZGEMM( 'N', 'N', M, N, M, CONE, WORK( IVT ), LDWKVT,                 CALL ZGEMM( 'N', 'N', M, N, M, CONE, WORK( IVT ), LDWKVT,
      $                     VT, LDVT, CZERO, A, LDA )       $                     VT, LDVT, CZERO, A, LDA )
Line 1553 Line 1781
 *  *
 *           MNTHR2 <= N < MNTHR1  *           MNTHR2 <= N < MNTHR1
 *  *
 *           Path 5t (N much larger than M, but not as much as MNTHR1)  *           Path 5t (N >> M, but not as much as MNTHR1)
 *           Reduce to bidiagonal form without QR decomposition, use  *           Reduce to bidiagonal form without QR decomposition, use
 *           ZUNGBR and matrix multiplication to compute singular vectors  *           ZUNGBR and matrix multiplication to compute singular vectors
 *  *
 *  
             IE = 1              IE = 1
             NRWORK = IE + M              NRWORK = IE + M
             ITAUQ = 1              ITAUQ = 1
Line 1565 Line 1792
             NWORK = ITAUP + M              NWORK = ITAUP + M
 *  *
 *           Bidiagonalize A  *           Bidiagonalize A
 *           (CWorkspace: need 2*M+N, prefer 2*M+(M+N)*NB)  *           CWorkspace: need   2*M [tauq, taup] + N        [work]
 *           (RWorkspace: M)  *           CWorkspace: prefer 2*M [tauq, taup] + (M+N)*NB [work]
   *           RWorkspace: need   M [e]
 *  *
             CALL ZGEBRD( M, N, A, LDA, S, RWORK( IE ), WORK( ITAUQ ),              CALL ZGEBRD( M, N, A, LDA, S, RWORK( IE ), WORK( ITAUQ ),
      $                   WORK( ITAUP ), WORK( NWORK ), LWORK-NWORK+1,       $                   WORK( ITAUP ), WORK( NWORK ), LWORK-NWORK+1,
Line 1574 Line 1802
 *  *
             IF( WNTQN ) THEN              IF( WNTQN ) THEN
 *  *
   *              Path 5tn (N >> M, JOBZ='N')
 *              Compute singular values only  *              Compute singular values only
 *              (Cworkspace: 0)  *              CWorkspace: need   0
 *              (Rworkspace: need BDSPAN)  *              RWorkspace: need   M [e] + BDSPAC
 *  *
                CALL DBDSDC( 'L', 'N', M, S, RWORK( IE ), DUM, 1, DUM, 1,                 CALL DBDSDC( 'L', 'N', M, S, RWORK( IE ), DUM,1,DUM,1,
      $                      DUM, IDUM, RWORK( NRWORK ), IWORK, INFO )       $                      DUM, IDUM, RWORK( NRWORK ), IWORK, INFO )
             ELSE IF( WNTQO ) THEN              ELSE IF( WNTQO ) THEN
                IRVT = NRWORK                 IRVT = NRWORK
Line 1586 Line 1815
                NRWORK = IRU + M*M                 NRWORK = IRU + M*M
                IVT = NWORK                 IVT = NWORK
 *  *
   *              Path 5to (N >> M, JOBZ='O')
 *              Copy A to U, generate Q  *              Copy A to U, generate Q
 *              (Cworkspace: need 2*M, prefer M+M*NB)  *              CWorkspace: need   2*M [tauq, taup] + M    [work]
 *              (Rworkspace: 0)  *              CWorkspace: prefer 2*M [tauq, taup] + M*NB [work]
   *              RWorkspace: need   0
 *  *
                CALL ZLACPY( 'L', M, M, A, LDA, U, LDU )                 CALL ZLACPY( 'L', M, M, A, LDA, U, LDU )
                CALL ZUNGBR( 'Q', M, M, N, U, LDU, WORK( ITAUQ ),                 CALL ZUNGBR( 'Q', M, M, N, U, LDU, WORK( ITAUQ ),
      $                      WORK( NWORK ), LWORK-NWORK+1, IERR )       $                      WORK( NWORK ), LWORK-NWORK+1, IERR )
 *  *
 *              Generate P**H in A  *              Generate P**H in A
 *              (Cworkspace: need 2*M, prefer M+M*NB)  *              CWorkspace: need   2*M [tauq, taup] + M    [work]
 *              (Rworkspace: 0)  *              CWorkspace: prefer 2*M [tauq, taup] + M*NB [work]
   *              RWorkspace: need   0
 *  *
                CALL ZUNGBR( 'P', M, N, M, A, LDA, WORK( ITAUP ),                 CALL ZUNGBR( 'P', M, N, M, A, LDA, WORK( ITAUP ),
      $                      WORK( NWORK ), LWORK-NWORK+1, IERR )       $                      WORK( NWORK ), LWORK-NWORK+1, IERR )
 *  *
                LDWKVT = M                 LDWKVT = M
                IF( LWORK.GE.M*N+3*M ) THEN                 IF( LWORK .GE. M*N + 3*M ) THEN
 *  *
 *                 WORK( IVT ) is M by N  *                 WORK( IVT ) is M by N
 *  *
Line 1612 Line 1844
 *  *
 *                 WORK( IVT ) is M by CHUNK  *                 WORK( IVT ) is M by CHUNK
 *  *
                   CHUNK = ( LWORK-3*M ) / M                    CHUNK = ( LWORK - 3*M ) / M
                   NWORK = IVT + LDWKVT*CHUNK                    NWORK = IVT + LDWKVT*CHUNK
                END IF                 END IF
 *  *
 *              Perform bidiagonal SVD, computing left singular vectors  *              Perform bidiagonal SVD, computing left singular vectors
 *              of bidiagonal matrix in RWORK(IRU) and computing right  *              of bidiagonal matrix in RWORK(IRU) and computing right
 *              singular vectors of bidiagonal matrix in RWORK(IRVT)  *              singular vectors of bidiagonal matrix in RWORK(IRVT)
 *              (CWorkspace: need 0)  *              CWorkspace: need   0
 *              (RWorkspace: need BDSPAC)  *              RWorkspace: need   M [e] + M*M [RVT] + M*M [RU] + BDSPAC
 *  *
                CALL DBDSDC( 'L', 'I', M, S, RWORK( IE ), RWORK( IRU ),                 CALL DBDSDC( 'L', 'I', M, S, RWORK( IE ), RWORK( IRU ),
      $                      M, RWORK( IRVT ), M, DUM, IDUM,       $                      M, RWORK( IRVT ), M, DUM, IDUM,
Line 1628 Line 1860
 *  *
 *              Multiply Q in U by real matrix RWORK(IRVT)  *              Multiply Q in U by real matrix RWORK(IRVT)
 *              storing the result in WORK(IVT), copying to U  *              storing the result in WORK(IVT), copying to U
 *              (Cworkspace: need 0)  *              CWorkspace: need   2*M [tauq, taup] + M*M [VT]
 *              (Rworkspace: need 2*M*M)  *              RWorkspace: need   M [e] + M*M [RVT] + M*M [RU] + 2*M*M [rwork]
 *  *
                CALL ZLACRM( M, M, U, LDU, RWORK( IRU ), M, WORK( IVT ),                 CALL ZLACRM( M, M, U, LDU, RWORK( IRU ), M, WORK( IVT ),
      $                      LDWKVT, RWORK( NRWORK ) )       $                      LDWKVT, RWORK( NRWORK ) )
Line 1637 Line 1869
 *  *
 *              Multiply RWORK(IRVT) by P**H in A, storing the  *              Multiply RWORK(IRVT) by P**H in A, storing the
 *              result in WORK(IVT), copying to A  *              result in WORK(IVT), copying to A
 *              (CWorkspace: need M*M, prefer M*N)  *              CWorkspace: need   2*M [tauq, taup] + M*M [VT]
 *              (Rworkspace: need 2*M*M, prefer 2*M*N)  *              CWorkspace: prefer 2*M [tauq, taup] + M*N [VT]
   *              RWorkspace: need   M [e] + M*M [RVT] + 2*M*M [rwork]
   *              RWorkspace: prefer M [e] + M*M [RVT] + 2*M*N [rwork] < M + 5*M*M since N < 2*M here
 *  *
                NRWORK = IRU                 NRWORK = IRU
                DO 50 I = 1, N, CHUNK                 DO 50 I = 1, N, CHUNK
Line 1650 Line 1884
    50          CONTINUE     50          CONTINUE
             ELSE IF( WNTQS ) THEN              ELSE IF( WNTQS ) THEN
 *  *
   *              Path 5ts (N >> M, JOBZ='S')
 *              Copy A to U, generate Q  *              Copy A to U, generate Q
 *              (Cworkspace: need 2*M, prefer M+M*NB)  *              CWorkspace: need   2*M [tauq, taup] + M    [work]
 *              (Rworkspace: 0)  *              CWorkspace: prefer 2*M [tauq, taup] + M*NB [work]
   *              RWorkspace: need   0
 *  *
                CALL ZLACPY( 'L', M, M, A, LDA, U, LDU )                 CALL ZLACPY( 'L', M, M, A, LDA, U, LDU )
                CALL ZUNGBR( 'Q', M, M, N, U, LDU, WORK( ITAUQ ),                 CALL ZUNGBR( 'Q', M, M, N, U, LDU, WORK( ITAUQ ),
      $                      WORK( NWORK ), LWORK-NWORK+1, IERR )       $                      WORK( NWORK ), LWORK-NWORK+1, IERR )
 *  *
 *              Copy A to VT, generate P**H  *              Copy A to VT, generate P**H
 *              (Cworkspace: need 2*M, prefer M+M*NB)  *              CWorkspace: need   2*M [tauq, taup] + M    [work]
 *              (Rworkspace: 0)  *              CWorkspace: prefer 2*M [tauq, taup] + M*NB [work]
   *              RWorkspace: need   0
 *  *
                CALL ZLACPY( 'U', M, N, A, LDA, VT, LDVT )                 CALL ZLACPY( 'U', M, N, A, LDA, VT, LDVT )
                CALL ZUNGBR( 'P', M, N, M, VT, LDVT, WORK( ITAUP ),                 CALL ZUNGBR( 'P', M, N, M, VT, LDVT, WORK( ITAUP ),
Line 1669 Line 1906
 *              Perform bidiagonal SVD, computing left singular vectors  *              Perform bidiagonal SVD, computing left singular vectors
 *              of bidiagonal matrix in RWORK(IRU) and computing right  *              of bidiagonal matrix in RWORK(IRU) and computing right
 *              singular vectors of bidiagonal matrix in RWORK(IRVT)  *              singular vectors of bidiagonal matrix in RWORK(IRVT)
 *              (CWorkspace: need 0)  *              CWorkspace: need   0
 *              (RWorkspace: need BDSPAC)  *              RWorkspace: need   M [e] + M*M [RVT] + M*M [RU] + BDSPAC
 *  *
                IRVT = NRWORK                 IRVT = NRWORK
                IRU = IRVT + M*M                 IRU = IRVT + M*M
Line 1681 Line 1918
 *  *
 *              Multiply Q in U by real matrix RWORK(IRU), storing the  *              Multiply Q in U by real matrix RWORK(IRU), storing the
 *              result in A, copying to U  *              result in A, copying to U
 *              (CWorkspace: need 0)  *              CWorkspace: need   0
 *              (Rworkspace: need 3*M*M)  *              RWorkspace: need   M [e] + M*M [RVT] + M*M [RU] + 2*M*M [rwork]
 *  *
                CALL ZLACRM( M, M, U, LDU, RWORK( IRU ), M, A, LDA,                 CALL ZLACRM( M, M, U, LDU, RWORK( IRU ), M, A, LDA,
      $                      RWORK( NRWORK ) )       $                      RWORK( NRWORK ) )
Line 1690 Line 1927
 *  *
 *              Multiply real matrix RWORK(IRVT) by P**H in VT,  *              Multiply real matrix RWORK(IRVT) by P**H in VT,
 *              storing the result in A, copying to VT  *              storing the result in A, copying to VT
 *              (Cworkspace: need 0)  *              CWorkspace: need   0
 *              (Rworkspace: need M*M+2*M*N)  *              RWorkspace: need   M [e] + M*M [RVT] + 2*M*N [rwork] < M + 5*M*M since N < 2*M here
 *  *
                NRWORK = IRU                 NRWORK = IRU
                CALL ZLARCM( M, N, RWORK( IRVT ), M, VT, LDVT, A, LDA,                 CALL ZLARCM( M, N, RWORK( IRVT ), M, VT, LDVT, A, LDA,
Line 1699 Line 1936
                CALL ZLACPY( 'F', M, N, A, LDA, VT, LDVT )                 CALL ZLACPY( 'F', M, N, A, LDA, VT, LDVT )
             ELSE              ELSE
 *  *
   *              Path 5ta (N >> M, JOBZ='A')
 *              Copy A to U, generate Q  *              Copy A to U, generate Q
 *              (Cworkspace: need 2*M, prefer M+M*NB)  *              CWorkspace: need   2*M [tauq, taup] + M    [work]
 *              (Rworkspace: 0)  *              CWorkspace: prefer 2*M [tauq, taup] + M*NB [work]
   *              RWorkspace: need   0
 *  *
                CALL ZLACPY( 'L', M, M, A, LDA, U, LDU )                 CALL ZLACPY( 'L', M, M, A, LDA, U, LDU )
                CALL ZUNGBR( 'Q', M, M, N, U, LDU, WORK( ITAUQ ),                 CALL ZUNGBR( 'Q', M, M, N, U, LDU, WORK( ITAUQ ),
      $                      WORK( NWORK ), LWORK-NWORK+1, IERR )       $                      WORK( NWORK ), LWORK-NWORK+1, IERR )
 *  *
 *              Copy A to VT, generate P**H  *              Copy A to VT, generate P**H
 *              (Cworkspace: need 2*M, prefer M+M*NB)  *              CWorkspace: need   2*M [tauq, taup] + N    [work]
 *              (Rworkspace: 0)  *              CWorkspace: prefer 2*M [tauq, taup] + N*NB [work]
   *              RWorkspace: need   0
 *  *
                CALL ZLACPY( 'U', M, N, A, LDA, VT, LDVT )                 CALL ZLACPY( 'U', M, N, A, LDA, VT, LDVT )
                CALL ZUNGBR( 'P', N, N, M, VT, LDVT, WORK( ITAUP ),                 CALL ZUNGBR( 'P', N, N, M, VT, LDVT, WORK( ITAUP ),
Line 1718 Line 1958
 *              Perform bidiagonal SVD, computing left singular vectors  *              Perform bidiagonal SVD, computing left singular vectors
 *              of bidiagonal matrix in RWORK(IRU) and computing right  *              of bidiagonal matrix in RWORK(IRU) and computing right
 *              singular vectors of bidiagonal matrix in RWORK(IRVT)  *              singular vectors of bidiagonal matrix in RWORK(IRVT)
 *              (CWorkspace: need 0)  *              CWorkspace: need   0
 *              (RWorkspace: need BDSPAC)  *              RWorkspace: need   M [e] + M*M [RVT] + M*M [RU] + BDSPAC
 *  *
                IRVT = NRWORK                 IRVT = NRWORK
                IRU = IRVT + M*M                 IRU = IRVT + M*M
Line 1730 Line 1970
 *  *
 *              Multiply Q in U by real matrix RWORK(IRU), storing the  *              Multiply Q in U by real matrix RWORK(IRU), storing the
 *              result in A, copying to U  *              result in A, copying to U
 *              (CWorkspace: need 0)  *              CWorkspace: need   0
 *              (Rworkspace: need 3*M*M)  *              RWorkspace: need   M [e] + M*M [RVT] + M*M [RU] + 2*M*M [rwork]
 *  *
                CALL ZLACRM( M, M, U, LDU, RWORK( IRU ), M, A, LDA,                 CALL ZLACRM( M, M, U, LDU, RWORK( IRU ), M, A, LDA,
      $                      RWORK( NRWORK ) )       $                      RWORK( NRWORK ) )
Line 1739 Line 1979
 *  *
 *              Multiply real matrix RWORK(IRVT) by P**H in VT,  *              Multiply real matrix RWORK(IRVT) by P**H in VT,
 *              storing the result in A, copying to VT  *              storing the result in A, copying to VT
 *              (Cworkspace: need 0)  *              CWorkspace: need   0
 *              (Rworkspace: need M*M+2*M*N)  *              RWorkspace: need   M [e] + M*M [RVT] + 2*M*N [rwork] < M + 5*M*M since N < 2*M here
 *  *
                  NRWORK = IRU
                CALL ZLARCM( M, N, RWORK( IRVT ), M, VT, LDVT, A, LDA,                 CALL ZLARCM( M, N, RWORK( IRVT ), M, VT, LDVT, A, LDA,
      $                      RWORK( NRWORK ) )       $                      RWORK( NRWORK ) )
                CALL ZLACPY( 'F', M, N, A, LDA, VT, LDVT )                 CALL ZLACPY( 'F', M, N, A, LDA, VT, LDVT )
Line 1751 Line 1992
 *  *
 *           N .LT. MNTHR2  *           N .LT. MNTHR2
 *  *
 *           Path 6t (N greater than M, but not much larger)  *           Path 6t (N > M, but not much larger)
 *           Reduce to bidiagonal form without LQ decomposition  *           Reduce to bidiagonal form without LQ decomposition
 *           Use ZUNMBR to compute singular vectors  *           Use ZUNMBR to compute singular vectors
 *  *
Line 1762 Line 2003
             NWORK = ITAUP + M              NWORK = ITAUP + M
 *  *
 *           Bidiagonalize A  *           Bidiagonalize A
 *           (CWorkspace: need 2*M+N, prefer 2*M+(M+N)*NB)  *           CWorkspace: need   2*M [tauq, taup] + N        [work]
 *           (RWorkspace: M)  *           CWorkspace: prefer 2*M [tauq, taup] + (M+N)*NB [work]
   *           RWorkspace: need   M [e]
 *  *
             CALL ZGEBRD( M, N, A, LDA, S, RWORK( IE ), WORK( ITAUQ ),              CALL ZGEBRD( M, N, A, LDA, S, RWORK( IE ), WORK( ITAUQ ),
      $                   WORK( ITAUP ), WORK( NWORK ), LWORK-NWORK+1,       $                   WORK( ITAUP ), WORK( NWORK ), LWORK-NWORK+1,
      $                   IERR )       $                   IERR )
             IF( WNTQN ) THEN              IF( WNTQN ) THEN
 *  *
   *              Path 6tn (N > M, JOBZ='N')
 *              Compute singular values only  *              Compute singular values only
 *              (Cworkspace: 0)  *              CWorkspace: need   0
 *              (Rworkspace: need BDSPAN)  *              RWorkspace: need   M [e] + BDSPAC
 *  *
                CALL DBDSDC( 'L', 'N', M, S, RWORK( IE ), DUM, 1, DUM, 1,                 CALL DBDSDC( 'L', 'N', M, S, RWORK( IE ), DUM,1,DUM,1,
      $                      DUM, IDUM, RWORK( NRWORK ), IWORK, INFO )       $                      DUM, IDUM, RWORK( NRWORK ), IWORK, INFO )
             ELSE IF( WNTQO ) THEN              ELSE IF( WNTQO ) THEN
   *              Path 6to (N > M, JOBZ='O')
                LDWKVT = M                 LDWKVT = M
                IVT = NWORK                 IVT = NWORK
                IF( LWORK.GE.M*N+3*M ) THEN                 IF( LWORK .GE. M*N + 3*M ) THEN
 *  *
 *                 WORK( IVT ) is M by N  *                 WORK( IVT ) is M by N
 *  *
Line 1790 Line 2034
 *  *
 *                 WORK( IVT ) is M by CHUNK  *                 WORK( IVT ) is M by CHUNK
 *  *
                   CHUNK = ( LWORK-3*M ) / M                    CHUNK = ( LWORK - 3*M ) / M
                   NWORK = IVT + LDWKVT*CHUNK                    NWORK = IVT + LDWKVT*CHUNK
                END IF                 END IF
 *  *
 *              Perform bidiagonal SVD, computing left singular vectors  *              Perform bidiagonal SVD, computing left singular vectors
 *              of bidiagonal matrix in RWORK(IRU) and computing right  *              of bidiagonal matrix in RWORK(IRU) and computing right
 *              singular vectors of bidiagonal matrix in RWORK(IRVT)  *              singular vectors of bidiagonal matrix in RWORK(IRVT)
 *              (CWorkspace: need 0)  *              CWorkspace: need   0
 *              (RWorkspace: need BDSPAC)  *              RWorkspace: need   M [e] + M*M [RVT] + M*M [RU] + BDSPAC
 *  *
                IRVT = NRWORK                 IRVT = NRWORK
                IRU = IRVT + M*M                 IRU = IRVT + M*M
Line 1809 Line 2053
 *  *
 *              Copy real matrix RWORK(IRU) to complex matrix U  *              Copy real matrix RWORK(IRU) to complex matrix U
 *              Overwrite U by left singular vectors of A  *              Overwrite U by left singular vectors of A
 *              (Cworkspace: need 2*M, prefer M+M*NB)  *              CWorkspace: need   2*M [tauq, taup] + M*M [VT] + M    [work]
 *              (Rworkspace: need 0)  *              CWorkspace: prefer 2*M [tauq, taup] + M*M [VT] + M*NB [work]
   *              RWorkspace: need   M [e] + M*M [RVT] + M*M [RU]
 *  *
                CALL ZLACP2( 'F', M, M, RWORK( IRU ), M, U, LDU )                 CALL ZLACP2( 'F', M, M, RWORK( IRU ), M, U, LDU )
                CALL ZUNMBR( 'Q', 'L', 'N', M, M, N, A, LDA,                 CALL ZUNMBR( 'Q', 'L', 'N', M, M, N, A, LDA,
      $                      WORK( ITAUQ ), U, LDU, WORK( NWORK ),       $                      WORK( ITAUQ ), U, LDU, WORK( NWORK ),
      $                      LWORK-NWORK+1, IERR )       $                      LWORK-NWORK+1, IERR )
 *  *
                IF( LWORK.GE.M*N+3*M ) THEN                 IF( LWORK .GE. M*N + 3*M ) THEN
 *  *
 *              Copy real matrix RWORK(IRVT) to complex matrix WORK(IVT)  *                 Path 6to-fast
 *              Overwrite WORK(IVT) by right singular vectors of A,  *                 Copy real matrix RWORK(IRVT) to complex matrix WORK(IVT)
 *              copying to A  *                 Overwrite WORK(IVT) by right singular vectors of A,
 *              (Cworkspace: need M*N+2*M, prefer M*N+M+M*NB)  *                 copying to A
 *              (Rworkspace: need 0)  *                 CWorkspace: need   2*M [tauq, taup] + M*N [VT] + M    [work]
   *                 CWorkspace: prefer 2*M [tauq, taup] + M*N [VT] + M*NB [work]
   *                 RWorkspace: need   M [e] + M*M [RVT]
 *  *
                   CALL ZLACP2( 'F', M, M, RWORK( IRVT ), M, WORK( IVT ),                    CALL ZLACP2( 'F', M, M, RWORK( IRVT ), M, WORK( IVT ),
      $                         LDWKVT )       $                         LDWKVT )
Line 1833 Line 2080
                   CALL ZLACPY( 'F', M, N, WORK( IVT ), LDWKVT, A, LDA )                    CALL ZLACPY( 'F', M, N, WORK( IVT ), LDWKVT, A, LDA )
                ELSE                 ELSE
 *  *
   *                 Path 6to-slow
 *                 Generate P**H in A  *                 Generate P**H in A
 *                 (Cworkspace: need 2*M, prefer M+M*NB)  *                 CWorkspace: need   2*M [tauq, taup] + M*M [VT] + M    [work]
 *                 (Rworkspace: need 0)  *                 CWorkspace: prefer 2*M [tauq, taup] + M*M [VT] + M*NB [work]
   *                 RWorkspace: need   0
 *  *
                   CALL ZUNGBR( 'P', M, N, M, A, LDA, WORK( ITAUP ),                    CALL ZUNGBR( 'P', M, N, M, A, LDA, WORK( ITAUP ),
      $                         WORK( NWORK ), LWORK-NWORK+1, IERR )       $                         WORK( NWORK ), LWORK-NWORK+1, IERR )
 *  *
 *                 Multiply Q in A by real matrix RWORK(IRU), storing the  *                 Multiply Q in A by real matrix RWORK(IRU), storing the
 *                 result in WORK(IU), copying to A  *                 result in WORK(IU), copying to A
 *                 (CWorkspace: need M*M, prefer M*N)  *                 CWorkspace: need   2*M [tauq, taup] + M*M [VT]
 *                 (Rworkspace: need 3*M*M, prefer M*M+2*M*N)  *                 CWorkspace: prefer 2*M [tauq, taup] + M*N [VT]
   *                 RWorkspace: need   M [e] + M*M [RVT] + 2*M*M [rwork]
   *                 RWorkspace: prefer M [e] + M*M [RVT] + 2*M*N [rwork] < M + 5*M*M since N < 2*M here
 *  *
                   NRWORK = IRU                    NRWORK = IRU
                   DO 60 I = 1, N, CHUNK                    DO 60 I = 1, N, CHUNK
Line 1857 Line 2108
                END IF                 END IF
             ELSE IF( WNTQS ) THEN              ELSE IF( WNTQS ) THEN
 *  *
   *              Path 6ts (N > M, JOBZ='S')
 *              Perform bidiagonal SVD, computing left singular vectors  *              Perform bidiagonal SVD, computing left singular vectors
 *              of bidiagonal matrix in RWORK(IRU) and computing right  *              of bidiagonal matrix in RWORK(IRU) and computing right
 *              singular vectors of bidiagonal matrix in RWORK(IRVT)  *              singular vectors of bidiagonal matrix in RWORK(IRVT)
 *              (CWorkspace: need 0)  *              CWorkspace: need   0
 *              (RWorkspace: need BDSPAC)  *              RWorkspace: need   M [e] + M*M [RVT] + M*M [RU] + BDSPAC
 *  *
                IRVT = NRWORK                 IRVT = NRWORK
                IRU = IRVT + M*M                 IRU = IRVT + M*M
Line 1872 Line 2124
 *  *
 *              Copy real matrix RWORK(IRU) to complex matrix U  *              Copy real matrix RWORK(IRU) to complex matrix U
 *              Overwrite U by left singular vectors of A  *              Overwrite U by left singular vectors of A
 *              (CWorkspace: need 3*M, prefer 2*M+M*NB)  *              CWorkspace: need   2*M [tauq, taup] + M    [work]
 *              (RWorkspace: M*M)  *              CWorkspace: prefer 2*M [tauq, taup] + M*NB [work]
   *              RWorkspace: need   M [e] + M*M [RVT] + M*M [RU]
 *  *
                CALL ZLACP2( 'F', M, M, RWORK( IRU ), M, U, LDU )                 CALL ZLACP2( 'F', M, M, RWORK( IRU ), M, U, LDU )
                CALL ZUNMBR( 'Q', 'L', 'N', M, M, N, A, LDA,                 CALL ZUNMBR( 'Q', 'L', 'N', M, M, N, A, LDA,
Line 1882 Line 2135
 *  *
 *              Copy real matrix RWORK(IRVT) to complex matrix VT  *              Copy real matrix RWORK(IRVT) to complex matrix VT
 *              Overwrite VT by right singular vectors of A  *              Overwrite VT by right singular vectors of A
 *              (CWorkspace: need 3*M, prefer 2*M+M*NB)  *              CWorkspace: need   2*M [tauq, taup] + M    [work]
 *              (RWorkspace: M*M)  *              CWorkspace: prefer 2*M [tauq, taup] + M*NB [work]
   *              RWorkspace: need   M [e] + M*M [RVT]
 *  *
                CALL ZLASET( 'F', M, N, CZERO, CZERO, VT, LDVT )                 CALL ZLASET( 'F', M, N, CZERO, CZERO, VT, LDVT )
                CALL ZLACP2( 'F', M, M, RWORK( IRVT ), M, VT, LDVT )                 CALL ZLACP2( 'F', M, M, RWORK( IRVT ), M, VT, LDVT )
Line 1892 Line 2146
      $                      LWORK-NWORK+1, IERR )       $                      LWORK-NWORK+1, IERR )
             ELSE              ELSE
 *  *
   *              Path 6ta (N > M, JOBZ='A')
 *              Perform bidiagonal SVD, computing left singular vectors  *              Perform bidiagonal SVD, computing left singular vectors
 *              of bidiagonal matrix in RWORK(IRU) and computing right  *              of bidiagonal matrix in RWORK(IRU) and computing right
 *              singular vectors of bidiagonal matrix in RWORK(IRVT)  *              singular vectors of bidiagonal matrix in RWORK(IRVT)
 *              (CWorkspace: need 0)  *              CWorkspace: need   0
 *              (RWorkspace: need BDSPAC)  *              RWorkspace: need   M [e] + M*M [RVT] + M*M [RU] + BDSPAC
 *  *
                IRVT = NRWORK                 IRVT = NRWORK
                IRU = IRVT + M*M                 IRU = IRVT + M*M
Line 1908 Line 2163
 *  *
 *              Copy real matrix RWORK(IRU) to complex matrix U  *              Copy real matrix RWORK(IRU) to complex matrix U
 *              Overwrite U by left singular vectors of A  *              Overwrite U by left singular vectors of A
 *              (CWorkspace: need 3*M, prefer 2*M+M*NB)  *              CWorkspace: need   2*M [tauq, taup] + M    [work]
 *              (RWorkspace: M*M)  *              CWorkspace: prefer 2*M [tauq, taup] + M*NB [work]
   *              RWorkspace: need   M [e] + M*M [RVT] + M*M [RU]
 *  *
                CALL ZLACP2( 'F', M, M, RWORK( IRU ), M, U, LDU )                 CALL ZLACP2( 'F', M, M, RWORK( IRU ), M, U, LDU )
                CALL ZUNMBR( 'Q', 'L', 'N', M, M, N, A, LDA,                 CALL ZUNMBR( 'Q', 'L', 'N', M, M, N, A, LDA,
Line 1922 Line 2178
 *  *
 *              Copy real matrix RWORK(IRVT) to complex matrix VT  *              Copy real matrix RWORK(IRVT) to complex matrix VT
 *              Overwrite VT by right singular vectors of A  *              Overwrite VT by right singular vectors of A
 *              (CWorkspace: need 2*M+N, prefer 2*M+N*NB)  *              CWorkspace: need   2*M [tauq, taup] + N    [work]
 *              (RWorkspace: M*M)  *              CWorkspace: prefer 2*M [tauq, taup] + N*NB [work]
   *              RWorkspace: need   M [e] + M*M [RVT]
 *  *
                CALL ZLACP2( 'F', M, M, RWORK( IRVT ), M, VT, LDVT )                 CALL ZLACP2( 'F', M, M, RWORK( IRVT ), M, VT, LDVT )
                CALL ZUNMBR( 'P', 'R', 'C', N, N, M, A, LDA,                 CALL ZUNMBR( 'P', 'R', 'C', N, N, M, A, LDA,

Removed from v.1.1  
changed lines
  Added in v.1.19


CVSweb interface <joel.bertrand@systella.fr>