Diff for /rpl/lapack/lapack/zgesdd.f between versions 1.5 and 1.18

version 1.5, 2010/08/07 13:18:09 version 1.18, 2017/06/17 11:06:43
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.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..--
 *     June 2010  *     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 >= min(M,N)*max(5*min(M,N)+7,2*max(M,N)+2*min(M,N)+1)  
 *  
 *  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 156 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 174 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 185 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 215 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 460 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 504 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 527 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 536 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 554 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 579 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 590 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 600 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 612 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 623 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 633 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 647 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 660 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 673 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 684 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 694 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 706 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 716 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 726 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 735 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 748 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 772 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 785 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 794 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 805 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 815 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 831 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 842 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 862 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 886 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 902 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 911 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 925 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 944 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 956 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 965 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 974 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 993 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 1005 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 1014 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 1027 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 1038 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 1066 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 1082 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 1108 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 1133 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 1148 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 1159 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 1168 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 1191 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 1201 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 1222 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 1245 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 1254 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 1272 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 1283 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 1302 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 1313 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 1323 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 1335 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 1345 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 1356 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 1370 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 1383 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 1396 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 1407 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 1417 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 1429 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 1439 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 1449 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 1458 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 1471 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 1495 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 1505 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 1517 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 1527 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 1538 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 1554 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 1566 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 1575 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 1587 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 1613 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 1629 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 1638 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 1651 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 1670 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 1682 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 1691 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 1700 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 1719 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 1731 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 1740 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 1752 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 1763 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 1791 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 1810 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 1834 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 1858 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 1873 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 1883 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 1893 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 1909 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 1923 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.5  
changed lines
  Added in v.1.18


CVSweb interface <joel.bertrand@systella.fr>