Diff for /rpl/lapack/lapack/zgesdd.f between versions 1.3 and 1.20

version 1.3, 2010/08/06 15:28:52 version 1.20, 2023/08/07 08:39:19
Line 1 Line 1
       SUBROUTINE ZGESDD( JOBZ, M, N, A, LDA, S, U, LDU, VT, LDVT, WORK,  *> \brief \b ZGESDD
      $                   LWORK, RWORK, IWORK, INFO )  
 *  *
 *  -- LAPACK driver routine (version 3.2) --  *  =========== DOCUMENTATION ===========
   *
   * Online html documentation available at
   *            http://www.netlib.org/lapack/explore-html/
   *
   *> \htmlonly
   *> Download ZGESDD + dependencies
   *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zgesdd.f">
   *> [TGZ]</a>
   *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zgesdd.f">
   *> [ZIP]</a>
   *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zgesdd.f">
   *> [TXT]</a>
   *> \endhtmlonly
   *
   *  Definition:
   *  ===========
   *
   *       SUBROUTINE ZGESDD( JOBZ, M, N, A, LDA, S, U, LDU, VT, LDVT,
   *                          WORK, LWORK, RWORK, IWORK, INFO )
   *
   *       .. Scalar Arguments ..
   *       CHARACTER          JOBZ
   *       INTEGER            INFO, LDA, LDU, LDVT, LWORK, M, N
   *       ..
   *       .. Array Arguments ..
   *       INTEGER            IWORK( * )
   *       DOUBLE PRECISION   RWORK( * ), S( * )
   *       COMPLEX*16         A( LDA, * ), U( LDU, * ), VT( LDVT, * ),
   *      $                   WORK( * )
   *       ..
   *
   *
   *> \par Purpose:
   *  =============
   *>
   *> \verbatim
   *>
   *> ZGESDD computes the singular value decomposition (SVD) of a complex
   *> M-by-N matrix A, optionally computing the left and/or right singular
   *> vectors, by using divide-and-conquer method. The SVD is written
   *>
   *>      A = U * SIGMA * conjugate-transpose(V)
   *>
   *> where SIGMA is an M-by-N matrix which is zero except for its
   *> min(m,n) diagonal elements, U is an M-by-M unitary matrix, and
   *> V is an N-by-N unitary matrix.  The diagonal elements of SIGMA
   *> are the singular values of A; they are real and non-negative, and
   *> are returned in descending order.  The first min(m,n) columns of
   *> U and V are the left and right singular vectors of A.
   *>
   *> Note that the routine returns VT = V**H, not V.
   *>
   *> The divide and conquer algorithm makes very mild assumptions about
   *> floating point arithmetic. It will work on machines with a guard
   *> digit in add/subtract, or on those binary machines without guard
   *> digits which subtract like the Cray X-MP, Cray Y-MP, Cray C-90, or
   *> Cray-2. It could conceivably fail on hexadecimal or decimal machines
   *> without guard digits, but we know of none.
   *> \endverbatim
   *
   *  Arguments:
   *  ==========
   *
   *> \param[in] JOBZ
   *> \verbatim
   *>          JOBZ is CHARACTER*1
   *>          Specifies options for computing all or part of the matrix U:
   *>          = 'A':  all M columns of U and all N rows of V**H are
   *>                  returned in the arrays U and VT;
   *>          = 'S':  the first min(M,N) columns of U and the first
   *>                  min(M,N) rows of V**H are returned in the arrays U
   *>                  and VT;
   *>          = 'O':  If M >= N, the first N columns of U are overwritten
   *>                  in the array A and all rows of V**H are returned in
   *>                  the array VT;
   *>                  otherwise, all columns of U are returned in the
   *>                  array U and the first M rows of V**H are overwritten
   *>                  in the array A;
   *>          = 'N':  no columns of U or rows of V**H are computed.
   *> \endverbatim
   *>
   *> \param[in] M
   *> \verbatim
   *>          M is INTEGER
   *>          The number of rows of the input matrix A.  M >= 0.
   *> \endverbatim
   *>
   *> \param[in] N
   *> \verbatim
   *>          N is INTEGER
   *>          The number of columns of the input matrix A.  N >= 0.
   *> \endverbatim
   *>
   *> \param[in,out] A
   *> \verbatim
   *>          A is COMPLEX*16 array, dimension (LDA,N)
   *>          On entry, the M-by-N matrix A.
   *>          On exit,
   *>          if JOBZ = 'O',  A is overwritten with the first N columns
   *>                          of U (the left singular vectors, stored
   *>                          columnwise) if M >= N;
   *>                          A is overwritten with the first M rows
   *>                          of V**H (the right singular vectors, stored
   *>                          rowwise) otherwise.
   *>          if JOBZ .ne. 'O', the contents of A are destroyed.
   *> \endverbatim
   *>
   *> \param[in] LDA
   *> \verbatim
   *>          LDA is INTEGER
   *>          The leading dimension of the array A.  LDA >= max(1,M).
   *> \endverbatim
   *>
   *> \param[out] S
   *> \verbatim
   *>          S is DOUBLE PRECISION array, dimension (min(M,N))
   *>          The singular values of A, sorted so that S(i) >= S(i+1).
   *> \endverbatim
   *>
   *> \param[out] U
   *> \verbatim
   *>          U is COMPLEX*16 array, dimension (LDU,UCOL)
   *>          UCOL = M if JOBZ = 'A' or JOBZ = 'O' and M < N;
   *>          UCOL = min(M,N) if JOBZ = 'S'.
   *>          If JOBZ = 'A' or JOBZ = 'O' and M < N, U contains the M-by-M
   *>          unitary matrix U;
   *>          if JOBZ = 'S', U contains the first min(M,N) columns of U
   *>          (the left singular vectors, stored columnwise);
   *>          if JOBZ = 'O' and M >= N, or JOBZ = 'N', U is not referenced.
   *> \endverbatim
   *>
   *> \param[in] LDU
   *> \verbatim
   *>          LDU is INTEGER
   *>          The leading dimension of the array U.  LDU >= 1;
   *>          if JOBZ = 'S' or 'A' or JOBZ = 'O' and M < N, LDU >= M.
   *> \endverbatim
   *>
   *> \param[out] VT
   *> \verbatim
   *>          VT is COMPLEX*16 array, dimension (LDVT,N)
   *>          If JOBZ = 'A' or JOBZ = 'O' and M >= N, VT contains the
   *>          N-by-N unitary matrix V**H;
   *>          if JOBZ = 'S', VT contains the first min(M,N) rows of
   *>          V**H (the right singular vectors, stored rowwise);
   *>          if JOBZ = 'O' and M < N, or JOBZ = 'N', VT is not referenced.
   *> \endverbatim
   *>
   *> \param[in] LDVT
   *> \verbatim
   *>          LDVT is INTEGER
   *>          The leading dimension of the array VT.  LDVT >= 1;
   *>          if JOBZ = 'A' or JOBZ = 'O' and M >= N, LDVT >= N;
   *>          if JOBZ = 'S', LDVT >= min(M,N).
   *> \endverbatim
   *>
   *> \param[out] WORK
   *> \verbatim
   *>          WORK is COMPLEX*16 array, dimension (MAX(1,LWORK))
   *>          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
   *> \endverbatim
   *>
   *> \param[in] LWORK
   *> \verbatim
   *>          LWORK is INTEGER
   *>          The dimension of the array WORK. LWORK >= 1.
   *>          If LWORK = -1, a workspace query is assumed.  The optimal
   *>          size for the WORK array is calculated and stored in WORK(1),
   *>          and no other work except argument checking is performed.
   *>
   *>          Let mx = max(M,N) and mn = min(M,N).
   *>          If JOBZ = 'N', LWORK >= 2*mn + mx.
   *>          If JOBZ = 'O', LWORK >= 2*mn*mn + 2*mn + mx.
   *>          If JOBZ = 'S', LWORK >=   mn*mn + 3*mn.
   *>          If JOBZ = 'A', LWORK >=   mn*mn + 2*mn + mx.
   *>          These are not tight minimums in all cases; see comments inside code.
   *>          For good performance, LWORK should generally be larger;
   *>          a query is recommended.
   *> \endverbatim
   *>
   *> \param[out] RWORK
   *> \verbatim
   *>          RWORK is DOUBLE PRECISION array, dimension (MAX(1,LRWORK))
   *>          Let mx = max(M,N) and mn = min(M,N).
   *>          If JOBZ = 'N',    LRWORK >= 5*mn (LAPACK <= 3.6 needs 7*mn);
   *>          else if mx >> mn, LRWORK >= 5*mn*mn + 5*mn;
   *>          else              LRWORK >= max( 5*mn*mn + 5*mn,
   *>                                           2*mx*mn + 2*mn*mn + mn ).
   *> \endverbatim
   *>
   *> \param[out] IWORK
   *> \verbatim
   *>          IWORK is INTEGER array, dimension (8*min(M,N))
   *> \endverbatim
   *>
   *> \param[out] INFO
   *> \verbatim
   *>          INFO is INTEGER
   *>          <  0:  if INFO = -i, the i-th argument had an illegal value.
   *>          = -4:  if A had a NAN entry.
   *>          >  0:  The updating process of DBDSDC did not converge.
   *>          =  0:  successful exit.
   *> \endverbatim
   *
   *  Authors:
   *  ========
   *
   *> \author Univ. of Tennessee
   *> \author Univ. of California Berkeley
   *> \author Univ. of Colorado Denver
   *> \author NAG Ltd.
   *
   *> \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 --
 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --  *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--  *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
 *     November 2006  
 *     8-15-00:  Improve consistency of WS calculations (eca)  
 *  *
 *     .. Scalar Arguments ..  *     .. Scalar Arguments ..
       CHARACTER          JOBZ        CHARACTER          JOBZ
Line 18 Line 241
      $                   WORK( * )       $                   WORK( * )
 *     ..  *     ..
 *  *
 *  Purpose  
 *  =======  
 *  
 *  ZGESDD computes the singular value decomposition (SVD) of a complex  
 *  M-by-N matrix A, optionally computing the left and/or right singular  
 *  vectors, by using divide-and-conquer method. The SVD is written  
 *  
 *       A = U * SIGMA * conjugate-transpose(V)  
 *  
 *  where SIGMA is an M-by-N matrix which is zero except for its  
 *  min(m,n) diagonal elements, U is an M-by-M unitary matrix, and  
 *  V is an N-by-N unitary matrix.  The diagonal elements of SIGMA  
 *  are the singular values of A; they are real and non-negative, and  
 *  are returned in descending order.  The first min(m,n) columns of  
 *  U and V are the left and right singular vectors of A.  
 *  
 *  Note that the routine returns VT = V**H, not V.  
 *  
 *  The divide and conquer algorithm makes very mild assumptions about  
 *  floating point arithmetic. It will work on machines with a guard  
 *  digit in add/subtract, or on those binary machines without guard  
 *  digits which subtract like the Cray X-MP, Cray Y-MP, Cray C-90, or  
 *  Cray-2. It could conceivably fail on hexadecimal or decimal machines  
 *  without guard digits, but we know of none.  
 *  
 *  Arguments  
 *  =========  
 *  
 *  JOBZ    (input) CHARACTER*1  
 *          Specifies options for computing all or part of the matrix U:  
 *          = 'A':  all M columns of U and all N rows of V**H are  
 *                  returned in the arrays U and VT;  
 *          = 'S':  the first min(M,N) columns of U and the first  
 *                  min(M,N) rows of V**H are returned in the arrays U  
 *                  and VT;  
 *          = 'O':  If M >= N, the first N columns of U are overwritten  
 *                  in the array A and all rows of V**H are returned in  
 *                  the array VT;  
 *                  otherwise, all columns of U are returned in the  
 *                  array U and the first M rows of V**H are overwritten  
 *                  in the array A;  
 *          = 'N':  no columns of U or rows of V**H are computed.  
 *  
 *  M       (input) INTEGER  
 *          The number of rows of the input matrix A.  M >= 0.  
 *  
 *  N       (input) INTEGER  
 *          The number of columns of the input matrix A.  N >= 0.  
 *  
 *  A       (input/output) COMPLEX*16 array, dimension (LDA,N)  
 *          On entry, the M-by-N matrix A.  
 *          On exit,  
 *          if JOBZ = 'O',  A is overwritten with the first N columns  
 *                          of U (the left singular vectors, stored  
 *                          columnwise) if M >= N;  
 *                          A is overwritten with the first M rows  
 *                          of V**H (the right singular vectors, stored  
 *                          rowwise) otherwise.  
 *          if JOBZ .ne. 'O', the contents of A are destroyed.  
 *  
 *  LDA     (input) INTEGER  
 *          The leading dimension of the array A.  LDA >= max(1,M).  
 *  
 *  S       (output) DOUBLE PRECISION array, dimension (min(M,N))  
 *          The singular values of A, sorted so that S(i) >= S(i+1).  
 *  
 *  U       (output) COMPLEX*16 array, dimension (LDU,UCOL)  
 *          UCOL = M if JOBZ = 'A' or JOBZ = 'O' and M < N;  
 *          UCOL = min(M,N) if JOBZ = 'S'.  
 *          If JOBZ = 'A' or JOBZ = 'O' and M < N, U contains the M-by-M  
 *          unitary matrix U;  
 *          if JOBZ = 'S', U contains the first min(M,N) columns of U  
 *          (the left singular vectors, stored columnwise);  
 *          if JOBZ = 'O' and M >= N, or JOBZ = 'N', U is not referenced.  
 *  
 *  LDU     (input) INTEGER  
 *          The leading dimension of the array U.  LDU >= 1; if  
 *          JOBZ = 'S' or 'A' or JOBZ = 'O' and M < N, LDU >= M.  
 *  
 *  VT      (output) COMPLEX*16 array, dimension (LDVT,N)  
 *          If JOBZ = 'A' or JOBZ = 'O' and M >= N, VT contains the  
 *          N-by-N unitary matrix V**H;  
 *          if JOBZ = 'S', VT contains the first min(M,N) rows of  
 *          V**H (the right singular vectors, stored rowwise);  
 *          if JOBZ = 'O' and M < N, or JOBZ = 'N', VT is not referenced.  
 *  
 *  LDVT    (input) INTEGER  
 *          The leading dimension of the array VT.  LDVT >= 1; if  
 *          JOBZ = 'A' or JOBZ = 'O' and M >= N, LDVT >= N;  
 *          if JOBZ = 'S', LDVT >= min(M,N).  
 *  
 *  WORK    (workspace/output) COMPLEX*16 array, dimension (MAX(1,LWORK))  
 *          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.  
 *  
 *  LWORK   (input) INTEGER  
 *          The dimension of the array WORK. LWORK >= 1.  
 *          if JOBZ = 'N', LWORK >= 2*min(M,N)+max(M,N).  
 *          if JOBZ = 'O',  
 *                LWORK >= 2*min(M,N)*min(M,N)+2*min(M,N)+max(M,N).  
 *          if JOBZ = 'S' or 'A',  
 *                LWORK >= min(M,N)*min(M,N)+2*min(M,N)+max(M,N).  
 *          For good performance, LWORK should generally be larger.  
 *  
 *          If LWORK = -1, a workspace query is assumed.  The optimal  
 *          size for the WORK array is calculated and stored in WORK(1),  
 *          and no other work except argument checking is performed.  
 *  
 *  RWORK   (workspace) DOUBLE PRECISION array, dimension (MAX(1,LRWORK))  
 *          If JOBZ = 'N', LRWORK >= 5*min(M,N).  
 *          Otherwise, LRWORK >= 5*min(M,N)*min(M,N) + 7*min(M,N)  
 *  
 *  IWORK   (workspace) INTEGER array, dimension (8*min(M,N))  
 *  
 *  INFO    (output) INTEGER  
 *          = 0:  successful exit.  
 *          < 0:  if INFO = -i, the i-th argument had an illegal value.  
 *          > 0:  The updating process of DBDSDC did not converge.  
 *  
 *  Further Details  
 *  ===============  
 *  
 *  Based on contributions by  
 *     Ming Gu and Huan Ren, Computer Science Division, University of  
 *     California at Berkeley, USA  
 *  
 *  =====================================================================  *  =====================================================================
 *  *
 *     .. Parameters ..  *     .. Parameters ..
       INTEGER            LQUERV  
       PARAMETER          ( LQUERV = -1 )  
       COMPLEX*16         CZERO, CONE        COMPLEX*16         CZERO, CONE
       PARAMETER          ( CZERO = ( 0.0D+0, 0.0D+0 ),        PARAMETER          ( CZERO = ( 0.0D+0, 0.0D+0 ),
      $                   CONE = ( 1.0D+0, 0.0D+0 ) )       $                   CONE = ( 1.0D+0, 0.0D+0 ) )
Line 155 Line 251
       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 172 Line 279
      $                   ZLASET, ZUNGBR, ZUNGLQ, ZUNGQR, ZUNMBR       $                   ZLASET, ZUNGBR, ZUNGLQ, ZUNGQR, ZUNMBR
 *     ..  *     ..
 *     .. External Functions ..  *     .. External Functions ..
       LOGICAL            LSAME        LOGICAL            LSAME, DISNAN
       INTEGER            ILAENV        DOUBLE PRECISION   DLAMCH, ZLANGE, DROUNDUP_LWORK
       DOUBLE PRECISION   DLAMCH, ZLANGE        EXTERNAL           LSAME, DLAMCH, ZLANGE, DISNAN, 
       EXTERNAL           LSAME, ILAENV, DLAMCH, ZLANGE       $                   DROUNDUP_LWORK
 *     ..  *     ..
 *     .. Intrinsic Functions ..  *     .. Intrinsic Functions ..
       INTRINSIC          INT, MAX, MIN, SQRT        INTRINSIC          INT, MAX, MIN, SQRT
Line 184 Line 291
 *  *
 *     Test the input arguments  *     Test the input arguments
 *  *
       INFO = 0        INFO   = 0
       MINMN = MIN( M, N )        MINMN  = MIN( M, N )
       MNTHR1 = INT( MINMN*17.0D0 / 9.0D0 )        MNTHR1 = INT( MINMN*17.0D0 / 9.0D0 )
       MNTHR2 = INT( MINMN*5.0D0 / 3.0D0 )        MNTHR2 = INT( MINMN*5.0D0 / 3.0D0 )
       WNTQA = LSAME( JOBZ, 'A' )        WNTQA  = LSAME( JOBZ, 'A' )
       WNTQS = LSAME( JOBZ, 'S' )        WNTQS  = LSAME( JOBZ, 'S' )
       WNTQAS = WNTQA .OR. WNTQS        WNTQAS = WNTQA .OR. WNTQS
       WNTQO = LSAME( JOBZ, 'O' )        WNTQO  = LSAME( JOBZ, 'O' )
       WNTQN = LSAME( JOBZ, 'N' )        WNTQN  = LSAME( JOBZ, 'N' )
         LQUERY = ( LWORK.EQ.-1 )
       MINWRK = 1        MINWRK = 1
       MAXWRK = 1        MAXWRK = 1
 *  *
Line 214 Line 322
       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
          MAXWRK = MAX( MAXWRK, MINWRK )           MAXWRK = MAX( MAXWRK, MINWRK )
       END IF        END IF
       IF( INFO.EQ.0 ) THEN        IF( INFO.EQ.0 ) THEN
          WORK( 1 ) = MAXWRK           WORK( 1 ) = DROUNDUP_LWORK( 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 484 Line 646
 *     Scale A if max element outside range [SMLNUM,BIGNUM]  *     Scale A if max element outside range [SMLNUM,BIGNUM]
 *  *
       ANRM = ZLANGE( 'M', M, N, A, LDA, DUM )        ANRM = ZLANGE( 'M', M, N, A, LDA, DUM )
         IF( DISNAN( ANRM ) ) THEN
             INFO = -4
             RETURN
         END IF
       ISCL = 0        ISCL = 0
       IF( ANRM.GT.ZERO .AND. ANRM.LT.SMLNUM ) THEN        IF( ANRM.GT.ZERO .AND. ANRM.LT.SMLNUM ) THEN
          ISCL = 1           ISCL = 1
Line 503 Line 669
 *  *
             IF( WNTQN ) THEN              IF( WNTQN ) THEN
 *  *
 *              Path 1 (M much larger than N, JOBZ='N')  *              Path 1 (M >> N, JOBZ='N')
 *              No singular vectors to be computed  *              No singular vectors to be computed
 *  *
                ITAU = 1                 ITAU = 1
                NWORK = ITAU + N                 NWORK = ITAU + N
 *  *
 *              Compute A=Q*R  *              Compute A=Q*R
 *              (CWorkspace: need 2*N, prefer N+N*NB)  *              CWorkspace: need   N [tau] + N    [work]
 *              (RWorkspace: need 0)  *              CWorkspace: prefer N [tau] + N*NB [work]
   *              RWorkspace: need   0
 *  *
                CALL ZGEQRF( M, N, A, LDA, WORK( ITAU ), WORK( NWORK ),                 CALL ZGEQRF( M, N, A, LDA, WORK( ITAU ), WORK( NWORK ),
      $                      LWORK-NWORK+1, IERR )       $                      LWORK-NWORK+1, IERR )
Line 526 Line 693
                NWORK = ITAUP + N                 NWORK = ITAUP + N
 *  *
 *              Bidiagonalize R in A  *              Bidiagonalize R in A
 *              (CWorkspace: need 3*N, prefer 2*N+2*N*NB)  *              CWorkspace: need   2*N [tauq, taup] + N      [work]
 *              (RWorkspace: need N)  *              CWorkspace: prefer 2*N [tauq, taup] + 2*N*NB [work]
   *              RWorkspace: need   N [e]
 *  *
                CALL ZGEBRD( N, N, A, LDA, S, RWORK( IE ), WORK( ITAUQ ),                 CALL ZGEBRD( N, N, A, LDA, S, RWORK( IE ), WORK( ITAUQ ),
      $                      WORK( ITAUP ), WORK( NWORK ), LWORK-NWORK+1,       $                      WORK( ITAUP ), WORK( NWORK ), LWORK-NWORK+1,
Line 535 Line 703
                NRWORK = IE + N                 NRWORK = IE + N
 *  *
 *              Perform bidiagonal SVD, compute singular values only  *              Perform bidiagonal SVD, compute singular values only
 *              (CWorkspace: 0)  *              CWorkspace: need   0
 *              (RWorkspace: need BDSPAN)  *              RWorkspace: need   N [e] + BDSPAC
 *  *
                CALL DBDSDC( 'U', 'N', N, S, RWORK( IE ), DUM, 1, DUM, 1,                 CALL DBDSDC( 'U', 'N', N, S, RWORK( IE ), DUM,1,DUM,1,
      $                      DUM, IDUM, RWORK( NRWORK ), IWORK, INFO )       $                      DUM, IDUM, RWORK( NRWORK ), IWORK, INFO )
 *  *
             ELSE IF( WNTQO ) THEN              ELSE IF( WNTQO ) THEN
 *  *
 *              Path 2 (M much larger than N, JOBZ='O')  *              Path 2 (M >> N, JOBZ='O')
 *              N left singular vectors to be overwritten on A and  *              N left singular vectors to be overwritten on A and
 *              N right singular vectors to be computed in VT  *              N right singular vectors to be computed in VT
 *  *
Line 553 Line 721
 *  *
                LDWRKU = N                 LDWRKU = N
                IR = IU + LDWRKU*N                 IR = IU + LDWRKU*N
                IF( LWORK.GE.M*N+N*N+3*N ) THEN                 IF( LWORK .GE. M*N + N*N + 3*N ) THEN
 *  *
 *                 WORK(IR) is M by N  *                 WORK(IR) is M by N
 *  *
                   LDWRKR = M                    LDWRKR = M
                ELSE                 ELSE
                   LDWRKR = ( LWORK-N*N-3*N ) / N                    LDWRKR = ( LWORK - N*N - 3*N ) / N
                END IF                 END IF
                ITAU = IR + LDWRKR*N                 ITAU = IR + LDWRKR*N
                NWORK = ITAU + N                 NWORK = ITAU + N
 *  *
 *              Compute A=Q*R  *              Compute A=Q*R
 *              (CWorkspace: need N*N+2*N, prefer M*N+N+N*NB)  *              CWorkspace: need   N*N [U] + N*N [R] + N [tau] + N    [work]
 *              (RWorkspace: 0)  *              CWorkspace: prefer N*N [U] + N*N [R] + N [tau] + N*NB [work]
   *              RWorkspace: need   0
 *  *
                CALL ZGEQRF( M, N, A, LDA, WORK( ITAU ), WORK( NWORK ),                 CALL ZGEQRF( M, N, A, LDA, WORK( ITAU ), WORK( NWORK ),
      $                      LWORK-NWORK+1, IERR )       $                      LWORK-NWORK+1, IERR )
Line 578 Line 747
      $                      LDWRKR )       $                      LDWRKR )
 *  *
 *              Generate Q in A  *              Generate Q in A
 *              (CWorkspace: need 2*N, prefer N+N*NB)  *              CWorkspace: need   N*N [U] + N*N [R] + N [tau] + N    [work]
 *              (RWorkspace: 0)  *              CWorkspace: prefer N*N [U] + N*N [R] + N [tau] + N*NB [work]
   *              RWorkspace: need   0
 *  *
                CALL ZUNGQR( M, N, N, A, LDA, WORK( ITAU ),                 CALL ZUNGQR( M, N, N, A, LDA, WORK( ITAU ),
      $                      WORK( NWORK ), LWORK-NWORK+1, IERR )       $                      WORK( NWORK ), LWORK-NWORK+1, IERR )
Line 589 Line 759
                NWORK = ITAUP + N                 NWORK = ITAUP + N
 *  *
 *              Bidiagonalize R in WORK(IR)  *              Bidiagonalize R in WORK(IR)
 *              (CWorkspace: need N*N+3*N, prefer M*N+2*N+2*N*NB)  *              CWorkspace: need   N*N [U] + N*N [R] + 2*N [tauq, taup] + N      [work]
 *              (RWorkspace: need N)  *              CWorkspace: prefer N*N [U] + N*N [R] + 2*N [tauq, taup] + 2*N*NB [work]
   *              RWorkspace: need   N [e]
 *  *
                CALL ZGEBRD( N, N, WORK( IR ), LDWRKR, S, RWORK( IE ),                 CALL ZGEBRD( N, N, WORK( IR ), LDWRKR, S, RWORK( IE ),
      $                      WORK( ITAUQ ), WORK( ITAUP ), WORK( NWORK ),       $                      WORK( ITAUQ ), WORK( ITAUP ), WORK( NWORK ),
Line 599 Line 770
 *              Perform bidiagonal SVD, computing left singular vectors  *              Perform bidiagonal SVD, computing left singular vectors
 *              of R in WORK(IRU) and computing right singular vectors  *              of R in WORK(IRU) and computing right singular vectors
 *              of R in WORK(IRVT)  *              of R in WORK(IRVT)
 *              (CWorkspace: need 0)  *              CWorkspace: need   0
 *              (RWorkspace: need BDSPAC)  *              RWorkspace: need   N [e] + N*N [RU] + N*N [RVT] + BDSPAC
 *  *
                IRU = IE + N                 IRU = IE + N
                IRVT = IRU + N*N                 IRVT = IRU + N*N
Line 611 Line 782
 *  *
 *              Copy real matrix RWORK(IRU) to complex matrix WORK(IU)  *              Copy real matrix RWORK(IRU) to complex matrix WORK(IU)
 *              Overwrite WORK(IU) by the left singular vectors of R  *              Overwrite WORK(IU) by the left singular vectors of R
 *              (CWorkspace: need 2*N*N+3*N, prefer M*N+N*N+2*N+N*NB)  *              CWorkspace: need   N*N [U] + N*N [R] + 2*N [tauq, taup] + N    [work]
 *              (RWorkspace: 0)  *              CWorkspace: prefer N*N [U] + N*N [R] + 2*N [tauq, taup] + N*NB [work]
   *              RWorkspace: need   0
 *  *
                CALL ZLACP2( 'F', N, N, RWORK( IRU ), N, WORK( IU ),                 CALL ZLACP2( 'F', N, N, RWORK( IRU ), N, WORK( IU ),
      $                      LDWRKU )       $                      LDWRKU )
Line 622 Line 794
 *  *
 *              Copy real matrix RWORK(IRVT) to complex matrix VT  *              Copy real matrix RWORK(IRVT) to complex matrix VT
 *              Overwrite VT by the right singular vectors of R  *              Overwrite VT by the right singular vectors of R
 *              (CWorkspace: need N*N+3*N, prefer M*N+2*N+N*NB)  *              CWorkspace: need   N*N [U] + N*N [R] + 2*N [tauq, taup] + N    [work]
 *              (RWorkspace: 0)  *              CWorkspace: prefer N*N [U] + N*N [R] + 2*N [tauq, taup] + N*NB [work]
   *              RWorkspace: need   0
 *  *
                CALL ZLACP2( 'F', N, N, RWORK( IRVT ), N, VT, LDVT )                 CALL ZLACP2( 'F', N, N, RWORK( IRVT ), N, VT, LDVT )
                CALL ZUNMBR( 'P', 'R', 'C', N, N, N, WORK( IR ), LDWRKR,                 CALL ZUNMBR( 'P', 'R', 'C', N, N, N, WORK( IR ), LDWRKR,
Line 632 Line 805
 *  *
 *              Multiply Q in A by left singular vectors of R in  *              Multiply Q in A by left singular vectors of R in
 *              WORK(IU), storing result in WORK(IR) and copying to A  *              WORK(IU), storing result in WORK(IR) and copying to A
 *              (CWorkspace: need 2*N*N, prefer N*N+M*N)  *              CWorkspace: need   N*N [U] + N*N [R]
 *              (RWorkspace: 0)  *              CWorkspace: prefer N*N [U] + M*N [R]
   *              RWorkspace: need   0
 *  *
                DO 10 I = 1, M, LDWRKR                 DO 10 I = 1, M, LDWRKR
                   CHUNK = MIN( M-I+1, LDWRKR )                    CHUNK = MIN( M-I+1, LDWRKR )
Line 646 Line 820
 *  *
             ELSE IF( WNTQS ) THEN              ELSE IF( WNTQS ) THEN
 *  *
 *              Path 3 (M much larger than N, JOBZ='S')  *              Path 3 (M >> N, JOBZ='S')
 *              N left singular vectors to be computed in U and  *              N left singular vectors to be computed in U and
 *              N right singular vectors to be computed in VT  *              N right singular vectors to be computed in VT
 *  *
Line 659 Line 833
                NWORK = ITAU + N                 NWORK = ITAU + N
 *  *
 *              Compute A=Q*R  *              Compute A=Q*R
 *              (CWorkspace: need N*N+2*N, prefer N*N+N+N*NB)  *              CWorkspace: need   N*N [R] + N [tau] + N    [work]
 *              (RWorkspace: 0)  *              CWorkspace: prefer N*N [R] + N [tau] + N*NB [work]
   *              RWorkspace: need   0
 *  *
                CALL ZGEQRF( M, N, A, LDA, WORK( ITAU ), WORK( NWORK ),                 CALL ZGEQRF( M, N, A, LDA, WORK( ITAU ), WORK( NWORK ),
      $                      LWORK-NWORK+1, IERR )       $                      LWORK-NWORK+1, IERR )
Line 672 Line 847
      $                      LDWRKR )       $                      LDWRKR )
 *  *
 *              Generate Q in A  *              Generate Q in A
 *              (CWorkspace: need 2*N, prefer N+N*NB)  *              CWorkspace: need   N*N [R] + N [tau] + N    [work]
 *              (RWorkspace: 0)  *              CWorkspace: prefer N*N [R] + N [tau] + N*NB [work]
   *              RWorkspace: need   0
 *  *
                CALL ZUNGQR( M, N, N, A, LDA, WORK( ITAU ),                 CALL ZUNGQR( M, N, N, A, LDA, WORK( ITAU ),
      $                      WORK( NWORK ), LWORK-NWORK+1, IERR )       $                      WORK( NWORK ), LWORK-NWORK+1, IERR )
Line 683 Line 859
                NWORK = ITAUP + N                 NWORK = ITAUP + N
 *  *
 *              Bidiagonalize R in WORK(IR)  *              Bidiagonalize R in WORK(IR)
 *              (CWorkspace: need N*N+3*N, prefer N*N+2*N+2*N*NB)  *              CWorkspace: need   N*N [R] + 2*N [tauq, taup] + N      [work]
 *              (RWorkspace: need N)  *              CWorkspace: prefer N*N [R] + 2*N [tauq, taup] + 2*N*NB [work]
   *              RWorkspace: need   N [e]
 *  *
                CALL ZGEBRD( N, N, WORK( IR ), LDWRKR, S, RWORK( IE ),                 CALL ZGEBRD( N, N, WORK( IR ), LDWRKR, S, RWORK( IE ),
      $                      WORK( ITAUQ ), WORK( ITAUP ), WORK( NWORK ),       $                      WORK( ITAUQ ), WORK( ITAUP ), WORK( NWORK ),
Line 693 Line 870
 *              Perform bidiagonal SVD, computing left singular vectors  *              Perform bidiagonal SVD, computing left singular vectors
 *              of bidiagonal matrix in RWORK(IRU) and computing right  *              of bidiagonal matrix in RWORK(IRU) and computing right
 *              singular vectors of bidiagonal matrix in RWORK(IRVT)  *              singular vectors of bidiagonal matrix in RWORK(IRVT)
 *              (CWorkspace: need 0)  *              CWorkspace: need   0
 *              (RWorkspace: need BDSPAC)  *              RWorkspace: need   N [e] + N*N [RU] + N*N [RVT] + BDSPAC
 *  *
                IRU = IE + N                 IRU = IE + N
                IRVT = IRU + N*N                 IRVT = IRU + N*N
Line 705 Line 882
 *  *
 *              Copy real matrix RWORK(IRU) to complex matrix U  *              Copy real matrix RWORK(IRU) to complex matrix U
 *              Overwrite U by left singular vectors of R  *              Overwrite U by left singular vectors of R
 *              (CWorkspace: need N*N+3*N, prefer N*N+2*N+N*NB)  *              CWorkspace: need   N*N [R] + 2*N [tauq, taup] + N    [work]
 *              (RWorkspace: 0)  *              CWorkspace: prefer N*N [R] + 2*N [tauq, taup] + N*NB [work]
   *              RWorkspace: need   0
 *  *
                CALL ZLACP2( 'F', N, N, RWORK( IRU ), N, U, LDU )                 CALL ZLACP2( 'F', N, N, RWORK( IRU ), N, U, LDU )
                CALL ZUNMBR( 'Q', 'L', 'N', N, N, N, WORK( IR ), LDWRKR,                 CALL ZUNMBR( 'Q', 'L', 'N', N, N, N, WORK( IR ), LDWRKR,
Line 715 Line 893
 *  *
 *              Copy real matrix RWORK(IRVT) to complex matrix VT  *              Copy real matrix RWORK(IRVT) to complex matrix VT
 *              Overwrite VT by right singular vectors of R  *              Overwrite VT by right singular vectors of R
 *              (CWorkspace: need N*N+3*N, prefer N*N+2*N+N*NB)  *              CWorkspace: need   N*N [R] + 2*N [tauq, taup] + N    [work]
 *              (RWorkspace: 0)  *              CWorkspace: prefer N*N [R] + 2*N [tauq, taup] + N*NB [work]
   *              RWorkspace: need   0
 *  *
                CALL ZLACP2( 'F', N, N, RWORK( IRVT ), N, VT, LDVT )                 CALL ZLACP2( 'F', N, N, RWORK( IRVT ), N, VT, LDVT )
                CALL ZUNMBR( 'P', 'R', 'C', N, N, N, WORK( IR ), LDWRKR,                 CALL ZUNMBR( 'P', 'R', 'C', N, N, N, WORK( IR ), LDWRKR,
Line 725 Line 904
 *  *
 *              Multiply Q in A by left singular vectors of R in  *              Multiply Q in A by left singular vectors of R in
 *              WORK(IR), storing result in U  *              WORK(IR), storing result in U
 *              (CWorkspace: need N*N)  *              CWorkspace: need   N*N [R]
 *              (RWorkspace: 0)  *              RWorkspace: need   0
 *  *
                CALL ZLACPY( 'F', N, N, U, LDU, WORK( IR ), LDWRKR )                 CALL ZLACPY( 'F', N, N, U, LDU, WORK( IR ), LDWRKR )
                CALL ZGEMM( 'N', 'N', M, N, N, CONE, A, LDA, WORK( IR ),                 CALL ZGEMM( 'N', 'N', M, N, N, CONE, A, LDA, WORK( IR ),
Line 734 Line 913
 *  *
             ELSE IF( WNTQA ) THEN              ELSE IF( WNTQA ) THEN
 *  *
 *              Path 4 (M much larger than N, JOBZ='A')  *              Path 4 (M >> N, JOBZ='A')
 *              M left singular vectors to be computed in U and  *              M left singular vectors to be computed in U and
 *              N right singular vectors to be computed in VT  *              N right singular vectors to be computed in VT
 *  *
Line 747 Line 926
                NWORK = ITAU + N                 NWORK = ITAU + N
 *  *
 *              Compute A=Q*R, copying result to U  *              Compute A=Q*R, copying result to U
 *              (CWorkspace: need 2*N, prefer N+N*NB)  *              CWorkspace: need   N*N [U] + N [tau] + N    [work]
 *              (RWorkspace: 0)  *              CWorkspace: prefer N*N [U] + N [tau] + N*NB [work]
   *              RWorkspace: need   0
 *  *
                CALL ZGEQRF( M, N, A, LDA, WORK( ITAU ), WORK( NWORK ),                 CALL ZGEQRF( M, N, A, LDA, WORK( ITAU ), WORK( NWORK ),
      $                      LWORK-NWORK+1, IERR )       $                      LWORK-NWORK+1, IERR )
                CALL ZLACPY( 'L', M, N, A, LDA, U, LDU )                 CALL ZLACPY( 'L', M, N, A, LDA, U, LDU )
 *  *
 *              Generate Q in U  *              Generate Q in U
 *              (CWorkspace: need N+M, prefer N+M*NB)  *              CWorkspace: need   N*N [U] + N [tau] + M    [work]
 *              (RWorkspace: 0)  *              CWorkspace: prefer N*N [U] + N [tau] + M*NB [work]
   *              RWorkspace: need   0
 *  *
                CALL ZUNGQR( M, M, N, U, LDU, WORK( ITAU ),                 CALL ZUNGQR( M, M, N, U, LDU, WORK( ITAU ),
      $                      WORK( NWORK ), LWORK-NWORK+1, IERR )       $                      WORK( NWORK ), LWORK-NWORK+1, IERR )
Line 771 Line 952
                NWORK = ITAUP + N                 NWORK = ITAUP + N
 *  *
 *              Bidiagonalize R in A  *              Bidiagonalize R in A
 *              (CWorkspace: need 3*N, prefer 2*N+2*N*NB)  *              CWorkspace: need   N*N [U] + 2*N [tauq, taup] + N      [work]
 *              (RWorkspace: need N)  *              CWorkspace: prefer N*N [U] + 2*N [tauq, taup] + 2*N*NB [work]
   *              RWorkspace: need   N [e]
 *  *
                CALL ZGEBRD( N, N, A, LDA, S, RWORK( IE ), WORK( ITAUQ ),                 CALL ZGEBRD( N, N, A, LDA, S, RWORK( IE ), WORK( ITAUQ ),
      $                      WORK( ITAUP ), WORK( NWORK ), LWORK-NWORK+1,       $                      WORK( ITAUP ), WORK( NWORK ), LWORK-NWORK+1,
Line 784 Line 966
 *              Perform bidiagonal SVD, computing left singular vectors  *              Perform bidiagonal SVD, computing left singular vectors
 *              of bidiagonal matrix in RWORK(IRU) and computing right  *              of bidiagonal matrix in RWORK(IRU) and computing right
 *              singular vectors of bidiagonal matrix in RWORK(IRVT)  *              singular vectors of bidiagonal matrix in RWORK(IRVT)
 *              (CWorkspace: need 0)  *              CWorkspace: need   0
 *              (RWorkspace: need BDSPAC)  *              RWorkspace: need   N [e] + N*N [RU] + N*N [RVT] + BDSPAC
 *  *
                CALL DBDSDC( 'U', 'I', N, S, RWORK( IE ), RWORK( IRU ),                 CALL DBDSDC( 'U', 'I', N, S, RWORK( IE ), RWORK( IRU ),
      $                      N, RWORK( IRVT ), N, DUM, IDUM,       $                      N, RWORK( IRVT ), N, DUM, IDUM,
Line 793 Line 975
 *  *
 *              Copy real matrix RWORK(IRU) to complex matrix WORK(IU)  *              Copy real matrix RWORK(IRU) to complex matrix WORK(IU)
 *              Overwrite WORK(IU) by left singular vectors of R  *              Overwrite WORK(IU) by left singular vectors of R
 *              (CWorkspace: need N*N+3*N, prefer N*N+2*N+N*NB)  *              CWorkspace: need   N*N [U] + 2*N [tauq, taup] + N    [work]
 *              (RWorkspace: 0)  *              CWorkspace: prefer N*N [U] + 2*N [tauq, taup] + N*NB [work]
   *              RWorkspace: need   0
 *  *
                CALL ZLACP2( 'F', N, N, RWORK( IRU ), N, WORK( IU ),                 CALL ZLACP2( 'F', N, N, RWORK( IRU ), N, WORK( IU ),
      $                      LDWRKU )       $                      LDWRKU )
Line 804 Line 987
 *  *
 *              Copy real matrix RWORK(IRVT) to complex matrix VT  *              Copy real matrix RWORK(IRVT) to complex matrix VT
 *              Overwrite VT by right singular vectors of R  *              Overwrite VT by right singular vectors of R
 *              (CWorkspace: need 3*N, prefer 2*N+N*NB)  *              CWorkspace: need   N*N [U] + 2*N [tauq, taup] + N    [work]
 *              (RWorkspace: 0)  *              CWorkspace: prefer N*N [U] + 2*N [tauq, taup] + N*NB [work]
   *              RWorkspace: need   0
 *  *
                CALL ZLACP2( 'F', N, N, RWORK( IRVT ), N, VT, LDVT )                 CALL ZLACP2( 'F', N, N, RWORK( IRVT ), N, VT, LDVT )
                CALL ZUNMBR( 'P', 'R', 'C', N, N, N, A, LDA,                 CALL ZUNMBR( 'P', 'R', 'C', N, N, N, A, LDA,
Line 814 Line 998
 *  *
 *              Multiply Q in U by left singular vectors of R in  *              Multiply Q in U by left singular vectors of R in
 *              WORK(IU), storing result in A  *              WORK(IU), storing result in A
 *              (CWorkspace: need N*N)  *              CWorkspace: need   N*N [U]
 *              (RWorkspace: 0)  *              RWorkspace: need   0
 *  *
                CALL ZGEMM( 'N', 'N', M, N, N, CONE, U, LDU, WORK( IU ),                 CALL ZGEMM( 'N', 'N', M, N, N, CONE, U, LDU, WORK( IU ),
      $                     LDWRKU, CZERO, A, LDA )       $                     LDWRKU, CZERO, A, LDA )
Line 830 Line 1014
 *  *
 *           MNTHR2 <= M < MNTHR1  *           MNTHR2 <= M < MNTHR1
 *  *
 *           Path 5 (M much larger than N, but not as much as MNTHR1)  *           Path 5 (M >> N, but not as much as MNTHR1)
 *           Reduce to bidiagonal form without QR decomposition, use  *           Reduce to bidiagonal form without QR decomposition, use
 *           ZUNGBR and matrix multiplication to compute singular vectors  *           ZUNGBR and matrix multiplication to compute singular vectors
 *  *
Line 841 Line 1025
             NWORK = ITAUP + N              NWORK = ITAUP + N
 *  *
 *           Bidiagonalize A  *           Bidiagonalize A
 *           (CWorkspace: need 2*N+M, prefer 2*N+(M+N)*NB)  *           CWorkspace: need   2*N [tauq, taup] + M        [work]
 *           (RWorkspace: need N)  *           CWorkspace: prefer 2*N [tauq, taup] + (M+N)*NB [work]
   *           RWorkspace: need   N [e]
 *  *
             CALL ZGEBRD( M, N, A, LDA, S, RWORK( IE ), WORK( ITAUQ ),              CALL ZGEBRD( M, N, A, LDA, S, RWORK( IE ), WORK( ITAUQ ),
      $                   WORK( ITAUP ), WORK( NWORK ), LWORK-NWORK+1,       $                   WORK( ITAUP ), WORK( NWORK ), LWORK-NWORK+1,
      $                   IERR )       $                   IERR )
             IF( WNTQN ) THEN              IF( WNTQN ) THEN
 *  *
   *              Path 5n (M >> N, JOBZ='N')
 *              Compute singular values only  *              Compute singular values only
 *              (Cworkspace: 0)  *              CWorkspace: need   0
 *              (Rworkspace: need BDSPAN)  *              RWorkspace: need   N [e] + BDSPAC
 *  *
                CALL DBDSDC( 'U', 'N', N, S, RWORK( IE ), DUM, 1, DUM, 1,                 CALL DBDSDC( 'U', 'N', N, S, RWORK( IE ), DUM, 1,DUM,1,
      $                      DUM, IDUM, RWORK( NRWORK ), IWORK, INFO )       $                      DUM, IDUM, RWORK( NRWORK ), IWORK, INFO )
             ELSE IF( WNTQO ) THEN              ELSE IF( WNTQO ) THEN
                IU = NWORK                 IU = NWORK
Line 861 Line 1047
                IRVT = IRU + N*N                 IRVT = IRU + N*N
                NRWORK = IRVT + N*N                 NRWORK = IRVT + N*N
 *  *
   *              Path 5o (M >> N, JOBZ='O')
 *              Copy A to VT, generate P**H  *              Copy A to VT, generate P**H
 *              (Cworkspace: need 2*N, prefer N+N*NB)  *              CWorkspace: need   2*N [tauq, taup] + N    [work]
 *              (Rworkspace: 0)  *              CWorkspace: prefer 2*N [tauq, taup] + N*NB [work]
   *              RWorkspace: need   0
 *  *
                CALL ZLACPY( 'U', N, N, A, LDA, VT, LDVT )                 CALL ZLACPY( 'U', N, N, A, LDA, VT, LDVT )
                CALL ZUNGBR( 'P', N, N, N, VT, LDVT, WORK( ITAUP ),                 CALL ZUNGBR( 'P', N, N, N, VT, LDVT, WORK( ITAUP ),
      $                      WORK( NWORK ), LWORK-NWORK+1, IERR )       $                      WORK( NWORK ), LWORK-NWORK+1, IERR )
 *  *
 *              Generate Q in A  *              Generate Q in A
 *              (CWorkspace: need 2*N, prefer N+N*NB)  *              CWorkspace: need   2*N [tauq, taup] + N    [work]
 *              (RWorkspace: 0)  *              CWorkspace: prefer 2*N [tauq, taup] + N*NB [work]
   *              RWorkspace: need   0
 *  *
                CALL ZUNGBR( 'Q', M, N, N, A, LDA, WORK( ITAUQ ),                 CALL ZUNGBR( 'Q', M, N, N, A, LDA, WORK( ITAUQ ),
      $                      WORK( NWORK ), LWORK-NWORK+1, IERR )       $                      WORK( NWORK ), LWORK-NWORK+1, IERR )
 *  *
                IF( LWORK.GE.M*N+3*N ) THEN                 IF( LWORK .GE. M*N + 3*N ) THEN
 *  *
 *                 WORK( IU ) is M by N  *                 WORK( IU ) is M by N
 *  *
Line 885 Line 1074
 *  *
 *                 WORK(IU) is LDWRKU by N  *                 WORK(IU) is LDWRKU by N
 *  *
                   LDWRKU = ( LWORK-3*N ) / N                    LDWRKU = ( LWORK - 3*N ) / N
                END IF                 END IF
                NWORK = IU + LDWRKU*N                 NWORK = IU + LDWRKU*N
 *  *
 *              Perform bidiagonal SVD, computing left singular vectors  *              Perform bidiagonal SVD, computing left singular vectors
 *              of bidiagonal matrix in RWORK(IRU) and computing right  *              of bidiagonal matrix in RWORK(IRU) and computing right
 *              singular vectors of bidiagonal matrix in RWORK(IRVT)  *              singular vectors of bidiagonal matrix in RWORK(IRVT)
 *              (CWorkspace: need 0)  *              CWorkspace: need   0
 *              (RWorkspace: need BDSPAC)  *              RWorkspace: need   N [e] + N*N [RU] + N*N [RVT] + BDSPAC
 *  *
                CALL DBDSDC( 'U', 'I', N, S, RWORK( IE ), RWORK( IRU ),                 CALL DBDSDC( 'U', 'I', N, S, RWORK( IE ), RWORK( IRU ),
      $                      N, RWORK( IRVT ), N, DUM, IDUM,       $                      N, RWORK( IRVT ), N, DUM, IDUM,
Line 901 Line 1090
 *  *
 *              Multiply real matrix RWORK(IRVT) by P**H in VT,  *              Multiply real matrix RWORK(IRVT) by P**H in VT,
 *              storing the result in WORK(IU), copying to VT  *              storing the result in WORK(IU), copying to VT
 *              (Cworkspace: need 0)  *              CWorkspace: need   2*N [tauq, taup] + N*N [U]
 *              (Rworkspace: need 3*N*N)  *              RWorkspace: need   N [e] + N*N [RU] + N*N [RVT] + 2*N*N [rwork]
 *  *
                CALL ZLARCM( N, N, RWORK( IRVT ), N, VT, LDVT,                 CALL ZLARCM( N, N, RWORK( IRVT ), N, VT, LDVT,
      $                      WORK( IU ), LDWRKU, RWORK( NRWORK ) )       $                      WORK( IU ), LDWRKU, RWORK( NRWORK ) )
Line 910 Line 1099
 *  *
 *              Multiply Q in A by real matrix RWORK(IRU), storing the  *              Multiply Q in A by real matrix RWORK(IRU), storing the
 *              result in WORK(IU), copying to A  *              result in WORK(IU), copying to A
 *              (CWorkspace: need N*N, prefer M*N)  *              CWorkspace: need   2*N [tauq, taup] + N*N [U]
 *              (Rworkspace: need 3*N*N, prefer N*N+2*M*N)  *              CWorkspace: prefer 2*N [tauq, taup] + M*N [U]
   *              RWorkspace: need   N [e] + N*N [RU] + 2*N*N [rwork]
   *              RWorkspace: prefer N [e] + N*N [RU] + 2*M*N [rwork] < N + 5*N*N since M < 2*N here
 *  *
                NRWORK = IRVT                 NRWORK = IRVT
                DO 20 I = 1, M, LDWRKU                 DO 20 I = 1, M, LDWRKU
Line 924 Line 1115
 *  *
             ELSE IF( WNTQS ) THEN              ELSE IF( WNTQS ) THEN
 *  *
   *              Path 5s (M >> N, JOBZ='S')
 *              Copy A to VT, generate P**H  *              Copy A to VT, generate P**H
 *              (Cworkspace: need 2*N, prefer N+N*NB)  *              CWorkspace: need   2*N [tauq, taup] + N    [work]
 *              (Rworkspace: 0)  *              CWorkspace: prefer 2*N [tauq, taup] + N*NB [work]
   *              RWorkspace: need   0
 *  *
                CALL ZLACPY( 'U', N, N, A, LDA, VT, LDVT )                 CALL ZLACPY( 'U', N, N, A, LDA, VT, LDVT )
                CALL ZUNGBR( 'P', N, N, N, VT, LDVT, WORK( ITAUP ),                 CALL ZUNGBR( 'P', N, N, N, VT, LDVT, WORK( ITAUP ),
      $                      WORK( NWORK ), LWORK-NWORK+1, IERR )       $                      WORK( NWORK ), LWORK-NWORK+1, IERR )
 *  *
 *              Copy A to U, generate Q  *              Copy A to U, generate Q
 *              (Cworkspace: need 2*N, prefer N+N*NB)  *              CWorkspace: need   2*N [tauq, taup] + N    [work]
 *              (Rworkspace: 0)  *              CWorkspace: prefer 2*N [tauq, taup] + N*NB [work]
   *              RWorkspace: need   0
 *  *
                CALL ZLACPY( 'L', M, N, A, LDA, U, LDU )                 CALL ZLACPY( 'L', M, N, A, LDA, U, LDU )
                CALL ZUNGBR( 'Q', M, N, N, U, LDU, WORK( ITAUQ ),                 CALL ZUNGBR( 'Q', M, N, N, U, LDU, WORK( ITAUQ ),
Line 943 Line 1137
 *              Perform bidiagonal SVD, computing left singular vectors  *              Perform bidiagonal SVD, computing left singular vectors
 *              of bidiagonal matrix in RWORK(IRU) and computing right  *              of bidiagonal matrix in RWORK(IRU) and computing right
 *              singular vectors of bidiagonal matrix in RWORK(IRVT)  *              singular vectors of bidiagonal matrix in RWORK(IRVT)
 *              (CWorkspace: need 0)  *              CWorkspace: need   0
 *              (RWorkspace: need BDSPAC)  *              RWorkspace: need   N [e] + N*N [RU] + N*N [RVT] + BDSPAC
 *  *
                IRU = NRWORK                 IRU = NRWORK
                IRVT = IRU + N*N                 IRVT = IRU + N*N
Line 955 Line 1149
 *  *
 *              Multiply real matrix RWORK(IRVT) by P**H in VT,  *              Multiply real matrix RWORK(IRVT) by P**H in VT,
 *              storing the result in A, copying to VT  *              storing the result in A, copying to VT
 *              (Cworkspace: need 0)  *              CWorkspace: need   0
 *              (Rworkspace: need 3*N*N)  *              RWorkspace: need   N [e] + N*N [RU] + N*N [RVT] + 2*N*N [rwork]
 *  *
                CALL ZLARCM( N, N, RWORK( IRVT ), N, VT, LDVT, A, LDA,                 CALL ZLARCM( N, N, RWORK( IRVT ), N, VT, LDVT, A, LDA,
      $                      RWORK( NRWORK ) )       $                      RWORK( NRWORK ) )
Line 964 Line 1158
 *  *
 *              Multiply Q in U by real matrix RWORK(IRU), storing the  *              Multiply Q in U by real matrix RWORK(IRU), storing the
 *              result in A, copying to U  *              result in A, copying to U
 *              (CWorkspace: need 0)  *              CWorkspace: need   0
 *              (Rworkspace: need N*N+2*M*N)  *              RWorkspace: need   N [e] + N*N [RU] + 2*M*N [rwork] < N + 5*N*N since M < 2*N here
 *  *
                NRWORK = IRVT                 NRWORK = IRVT
                CALL ZLACRM( M, N, U, LDU, RWORK( IRU ), N, A, LDA,                 CALL ZLACRM( M, N, U, LDU, RWORK( IRU ), N, A, LDA,
Line 973 Line 1167
                CALL ZLACPY( 'F', M, N, A, LDA, U, LDU )                 CALL ZLACPY( 'F', M, N, A, LDA, U, LDU )
             ELSE              ELSE
 *  *
   *              Path 5a (M >> N, JOBZ='A')
 *              Copy A to VT, generate P**H  *              Copy A to VT, generate P**H
 *              (Cworkspace: need 2*N, prefer N+N*NB)  *              CWorkspace: need   2*N [tauq, taup] + N    [work]
 *              (Rworkspace: 0)  *              CWorkspace: prefer 2*N [tauq, taup] + N*NB [work]
   *              RWorkspace: need   0
 *  *
                CALL ZLACPY( 'U', N, N, A, LDA, VT, LDVT )                 CALL ZLACPY( 'U', N, N, A, LDA, VT, LDVT )
                CALL ZUNGBR( 'P', N, N, N, VT, LDVT, WORK( ITAUP ),                 CALL ZUNGBR( 'P', N, N, N, VT, LDVT, WORK( ITAUP ),
      $                      WORK( NWORK ), LWORK-NWORK+1, IERR )       $                      WORK( NWORK ), LWORK-NWORK+1, IERR )
 *  *
 *              Copy A to U, generate Q  *              Copy A to U, generate Q
 *              (Cworkspace: need 2*N, prefer N+N*NB)  *              CWorkspace: need   2*N [tauq, taup] + M    [work]
 *              (Rworkspace: 0)  *              CWorkspace: prefer 2*N [tauq, taup] + M*NB [work]
   *              RWorkspace: need   0
 *  *
                CALL ZLACPY( 'L', M, N, A, LDA, U, LDU )                 CALL ZLACPY( 'L', M, N, A, LDA, U, LDU )
                CALL ZUNGBR( 'Q', M, M, N, U, LDU, WORK( ITAUQ ),                 CALL ZUNGBR( 'Q', M, M, N, U, LDU, WORK( ITAUQ ),
Line 992 Line 1189
 *              Perform bidiagonal SVD, computing left singular vectors  *              Perform bidiagonal SVD, computing left singular vectors
 *              of bidiagonal matrix in RWORK(IRU) and computing right  *              of bidiagonal matrix in RWORK(IRU) and computing right
 *              singular vectors of bidiagonal matrix in RWORK(IRVT)  *              singular vectors of bidiagonal matrix in RWORK(IRVT)
 *              (CWorkspace: need 0)  *              CWorkspace: need   0
 *              (RWorkspace: need BDSPAC)  *              RWorkspace: need   N [e] + N*N [RU] + N*N [RVT] + BDSPAC
 *  *
                IRU = NRWORK                 IRU = NRWORK
                IRVT = IRU + N*N                 IRVT = IRU + N*N
Line 1004 Line 1201
 *  *
 *              Multiply real matrix RWORK(IRVT) by P**H in VT,  *              Multiply real matrix RWORK(IRVT) by P**H in VT,
 *              storing the result in A, copying to VT  *              storing the result in A, copying to VT
 *              (Cworkspace: need 0)  *              CWorkspace: need   0
 *              (Rworkspace: need 3*N*N)  *              RWorkspace: need   N [e] + N*N [RU] + N*N [RVT] + 2*N*N [rwork]
 *  *
                CALL ZLARCM( N, N, RWORK( IRVT ), N, VT, LDVT, A, LDA,                 CALL ZLARCM( N, N, RWORK( IRVT ), N, VT, LDVT, A, LDA,
      $                      RWORK( NRWORK ) )       $                      RWORK( NRWORK ) )
Line 1013 Line 1210
 *  *
 *              Multiply Q in U by real matrix RWORK(IRU), storing the  *              Multiply Q in U by real matrix RWORK(IRU), storing the
 *              result in A, copying to U  *              result in A, copying to U
 *              (CWorkspace: 0)  *              CWorkspace: need   0
 *              (Rworkspace: need 3*N*N)  *              RWorkspace: need   N [e] + N*N [RU] + 2*M*N [rwork] < N + 5*N*N since M < 2*N here
 *  *
                NRWORK = IRVT                 NRWORK = IRVT
                CALL ZLACRM( M, N, U, LDU, RWORK( IRU ), N, A, LDA,                 CALL ZLACRM( M, N, U, LDU, RWORK( IRU ), N, A, LDA,
Line 1026 Line 1223
 *  *
 *           M .LT. MNTHR2  *           M .LT. MNTHR2
 *  *
 *           Path 6 (M at least N, but not much larger)  *           Path 6 (M >= N, but not much larger)
 *           Reduce to bidiagonal form without QR decomposition  *           Reduce to bidiagonal form without QR decomposition
 *           Use ZUNMBR to compute singular vectors  *           Use ZUNMBR to compute singular vectors
 *  *
Line 1037 Line 1234
             NWORK = ITAUP + N              NWORK = ITAUP + N
 *  *
 *           Bidiagonalize A  *           Bidiagonalize A
 *           (CWorkspace: need 2*N+M, prefer 2*N+(M+N)*NB)  *           CWorkspace: need   2*N [tauq, taup] + M        [work]
 *           (RWorkspace: need N)  *           CWorkspace: prefer 2*N [tauq, taup] + (M+N)*NB [work]
   *           RWorkspace: need   N [e]
 *  *
             CALL ZGEBRD( M, N, A, LDA, S, RWORK( IE ), WORK( ITAUQ ),              CALL ZGEBRD( M, N, A, LDA, S, RWORK( IE ), WORK( ITAUQ ),
      $                   WORK( ITAUP ), WORK( NWORK ), LWORK-NWORK+1,       $                   WORK( ITAUP ), WORK( NWORK ), LWORK-NWORK+1,
      $                   IERR )       $                   IERR )
             IF( WNTQN ) THEN              IF( WNTQN ) THEN
 *  *
   *              Path 6n (M >= N, JOBZ='N')
 *              Compute singular values only  *              Compute singular values only
 *              (Cworkspace: 0)  *              CWorkspace: need   0
 *              (Rworkspace: need BDSPAN)  *              RWorkspace: need   N [e] + BDSPAC
 *  *
                CALL DBDSDC( 'U', 'N', N, S, RWORK( IE ), DUM, 1, DUM, 1,                 CALL DBDSDC( 'U', 'N', N, S, RWORK( IE ), DUM,1,DUM,1,
      $                      DUM, IDUM, RWORK( NRWORK ), IWORK, INFO )       $                      DUM, IDUM, RWORK( NRWORK ), IWORK, INFO )
             ELSE IF( WNTQO ) THEN              ELSE IF( WNTQO ) THEN
                IU = NWORK                 IU = NWORK
                IRU = NRWORK                 IRU = NRWORK
                IRVT = IRU + N*N                 IRVT = IRU + N*N
                NRWORK = IRVT + N*N                 NRWORK = IRVT + N*N
                IF( LWORK.GE.M*N+3*N ) THEN                 IF( LWORK .GE. M*N + 3*N ) THEN
 *  *
 *                 WORK( IU ) is M by N  *                 WORK( IU ) is M by N
 *  *
Line 1065 Line 1264
 *  *
 *                 WORK( IU ) is LDWRKU by N  *                 WORK( IU ) is LDWRKU by N
 *  *
                   LDWRKU = ( LWORK-3*N ) / N                    LDWRKU = ( LWORK - 3*N ) / N
                END IF                 END IF
                NWORK = IU + LDWRKU*N                 NWORK = IU + LDWRKU*N
 *  *
   *              Path 6o (M >= N, JOBZ='O')
 *              Perform bidiagonal SVD, computing left singular vectors  *              Perform bidiagonal SVD, computing left singular vectors
 *              of bidiagonal matrix in RWORK(IRU) and computing right  *              of bidiagonal matrix in RWORK(IRU) and computing right
 *              singular vectors of bidiagonal matrix in RWORK(IRVT)  *              singular vectors of bidiagonal matrix in RWORK(IRVT)
 *              (CWorkspace: need 0)  *              CWorkspace: need   0
 *              (RWorkspace: need BDSPAC)  *              RWorkspace: need   N [e] + N*N [RU] + N*N [RVT] + BDSPAC
 *  *
                CALL DBDSDC( 'U', 'I', N, S, RWORK( IE ), RWORK( IRU ),                 CALL DBDSDC( 'U', 'I', N, S, RWORK( IE ), RWORK( IRU ),
      $                      N, RWORK( IRVT ), N, DUM, IDUM,       $                      N, RWORK( IRVT ), N, DUM, IDUM,
Line 1081 Line 1281
 *  *
 *              Copy real matrix RWORK(IRVT) to complex matrix VT  *              Copy real matrix RWORK(IRVT) to complex matrix VT
 *              Overwrite VT by right singular vectors of A  *              Overwrite VT by right singular vectors of A
 *              (Cworkspace: need 2*N, prefer N+N*NB)  *              CWorkspace: need   2*N [tauq, taup] + N*N [U] + N    [work]
 *              (Rworkspace: need 0)  *              CWorkspace: prefer 2*N [tauq, taup] + N*N [U] + N*NB [work]
   *              RWorkspace: need   N [e] + N*N [RU] + N*N [RVT]
 *  *
                CALL ZLACP2( 'F', N, N, RWORK( IRVT ), N, VT, LDVT )                 CALL ZLACP2( 'F', N, N, RWORK( IRVT ), N, VT, LDVT )
                CALL ZUNMBR( 'P', 'R', 'C', N, N, N, A, LDA,                 CALL ZUNMBR( 'P', 'R', 'C', N, N, N, A, LDA,
      $                      WORK( ITAUP ), VT, LDVT, WORK( NWORK ),       $                      WORK( ITAUP ), VT, LDVT, WORK( NWORK ),
      $                      LWORK-NWORK+1, IERR )       $                      LWORK-NWORK+1, IERR )
 *  *
                IF( LWORK.GE.M*N+3*N ) THEN                 IF( LWORK .GE. M*N + 3*N ) THEN
 *  *
 *              Copy real matrix RWORK(IRU) to complex matrix WORK(IU)  *                 Path 6o-fast
 *              Overwrite WORK(IU) by left singular vectors of A, copying  *                 Copy real matrix RWORK(IRU) to complex matrix WORK(IU)
 *              to A  *                 Overwrite WORK(IU) by left singular vectors of A, copying
 *              (Cworkspace: need M*N+2*N, prefer M*N+N+N*NB)  *                 to A
 *              (Rworkspace: need 0)  *                 CWorkspace: need   2*N [tauq, taup] + M*N [U] + N    [work]
   *                 CWorkspace: prefer 2*N [tauq, taup] + M*N [U] + N*NB [work]
   *                 RWorkspace: need   N [e] + N*N [RU]
 *  *
                   CALL ZLASET( 'F', M, N, CZERO, CZERO, WORK( IU ),                    CALL ZLASET( 'F', M, N, CZERO, CZERO, WORK( IU ),
      $                         LDWRKU )       $                         LDWRKU )
Line 1107 Line 1310
                   CALL ZLACPY( 'F', M, N, WORK( IU ), LDWRKU, A, LDA )                    CALL ZLACPY( 'F', M, N, WORK( IU ), LDWRKU, A, LDA )
                ELSE                 ELSE
 *  *
   *                 Path 6o-slow
 *                 Generate Q in A  *                 Generate Q in A
 *                 (Cworkspace: need 2*N, prefer N+N*NB)  *                 CWorkspace: need   2*N [tauq, taup] + N*N [U] + N    [work]
 *                 (Rworkspace: need 0)  *                 CWorkspace: prefer 2*N [tauq, taup] + N*N [U] + N*NB [work]
   *                 RWorkspace: need   0
 *  *
                   CALL ZUNGBR( 'Q', M, N, N, A, LDA, WORK( ITAUQ ),                    CALL ZUNGBR( 'Q', M, N, N, A, LDA, WORK( ITAUQ ),
      $                         WORK( NWORK ), LWORK-NWORK+1, IERR )       $                         WORK( NWORK ), LWORK-NWORK+1, IERR )
 *  *
 *                 Multiply Q in A by real matrix RWORK(IRU), storing the  *                 Multiply Q in A by real matrix RWORK(IRU), storing the
 *                 result in WORK(IU), copying to A  *                 result in WORK(IU), copying to A
 *                 (CWorkspace: need N*N, prefer M*N)  *                 CWorkspace: need   2*N [tauq, taup] + N*N [U]
 *                 (Rworkspace: need 3*N*N, prefer N*N+2*M*N)  *                 CWorkspace: prefer 2*N [tauq, taup] + M*N [U]
   *                 RWorkspace: need   N [e] + N*N [RU] + 2*N*N [rwork]
   *                 RWorkspace: prefer N [e] + N*N [RU] + 2*M*N [rwork] < N + 5*N*N since M < 2*N here
 *  *
                   NRWORK = IRVT                    NRWORK = IRVT
                   DO 30 I = 1, M, LDWRKU                    DO 30 I = 1, M, LDWRKU
Line 1132 Line 1339
 *  *
             ELSE IF( WNTQS ) THEN              ELSE IF( WNTQS ) THEN
 *  *
   *              Path 6s (M >= N, JOBZ='S')
 *              Perform bidiagonal SVD, computing left singular vectors  *              Perform bidiagonal SVD, computing left singular vectors
 *              of bidiagonal matrix in RWORK(IRU) and computing right  *              of bidiagonal matrix in RWORK(IRU) and computing right
 *              singular vectors of bidiagonal matrix in RWORK(IRVT)  *              singular vectors of bidiagonal matrix in RWORK(IRVT)
 *              (CWorkspace: need 0)  *              CWorkspace: need   0
 *              (RWorkspace: need BDSPAC)  *              RWorkspace: need   N [e] + N*N [RU] + N*N [RVT] + BDSPAC
 *  *
                IRU = NRWORK                 IRU = NRWORK
                IRVT = IRU + N*N                 IRVT = IRU + N*N
Line 1147 Line 1355
 *  *
 *              Copy real matrix RWORK(IRU) to complex matrix U  *              Copy real matrix RWORK(IRU) to complex matrix U
 *              Overwrite U by left singular vectors of A  *              Overwrite U by left singular vectors of A
 *              (CWorkspace: need 3*N, prefer 2*N+N*NB)  *              CWorkspace: need   2*N [tauq, taup] + N    [work]
 *              (RWorkspace: 0)  *              CWorkspace: prefer 2*N [tauq, taup] + N*NB [work]
   *              RWorkspace: need   N [e] + N*N [RU] + N*N [RVT]
 *  *
                CALL ZLASET( 'F', M, N, CZERO, CZERO, U, LDU )                 CALL ZLASET( 'F', M, N, CZERO, CZERO, U, LDU )
                CALL ZLACP2( 'F', N, N, RWORK( IRU ), N, U, LDU )                 CALL ZLACP2( 'F', N, N, RWORK( IRU ), N, U, LDU )
Line 1158 Line 1367
 *  *
 *              Copy real matrix RWORK(IRVT) to complex matrix VT  *              Copy real matrix RWORK(IRVT) to complex matrix VT
 *              Overwrite VT by right singular vectors of A  *              Overwrite VT by right singular vectors of A
 *              (CWorkspace: need 3*N, prefer 2*N+N*NB)  *              CWorkspace: need   2*N [tauq, taup] + N    [work]
 *              (RWorkspace: 0)  *              CWorkspace: prefer 2*N [tauq, taup] + N*NB [work]
   *              RWorkspace: need   N [e] + N*N [RU] + N*N [RVT]
 *  *
                CALL ZLACP2( 'F', N, N, RWORK( IRVT ), N, VT, LDVT )                 CALL ZLACP2( 'F', N, N, RWORK( IRVT ), N, VT, LDVT )
                CALL ZUNMBR( 'P', 'R', 'C', N, N, N, A, LDA,                 CALL ZUNMBR( 'P', 'R', 'C', N, N, N, A, LDA,
Line 1167 Line 1377
      $                      LWORK-NWORK+1, IERR )       $                      LWORK-NWORK+1, IERR )
             ELSE              ELSE
 *  *
   *              Path 6a (M >= N, JOBZ='A')
 *              Perform bidiagonal SVD, computing left singular vectors  *              Perform bidiagonal SVD, computing left singular vectors
 *              of bidiagonal matrix in RWORK(IRU) and computing right  *              of bidiagonal matrix in RWORK(IRU) and computing right
 *              singular vectors of bidiagonal matrix in RWORK(IRVT)  *              singular vectors of bidiagonal matrix in RWORK(IRVT)
 *              (CWorkspace: need 0)  *              CWorkspace: need   0
 *              (RWorkspace: need BDSPAC)  *              RWorkspace: need   N [e] + N*N [RU] + N*N [RVT] + BDSPAC
 *  *
                IRU = NRWORK                 IRU = NRWORK
                IRVT = IRU + N*N                 IRVT = IRU + N*N
Line 1190 Line 1401
 *  *
 *              Copy real matrix RWORK(IRU) to complex matrix U  *              Copy real matrix RWORK(IRU) to complex matrix U
 *              Overwrite U by left singular vectors of A  *              Overwrite U by left singular vectors of A
 *              (CWorkspace: need 2*N+M, prefer 2*N+M*NB)  *              CWorkspace: need   2*N [tauq, taup] + M    [work]
 *              (RWorkspace: 0)  *              CWorkspace: prefer 2*N [tauq, taup] + M*NB [work]
   *              RWorkspace: need   N [e] + N*N [RU] + N*N [RVT]
 *  *
                CALL ZLACP2( 'F', N, N, RWORK( IRU ), N, U, LDU )                 CALL ZLACP2( 'F', N, N, RWORK( IRU ), N, U, LDU )
                CALL ZUNMBR( 'Q', 'L', 'N', M, M, N, A, LDA,                 CALL ZUNMBR( 'Q', 'L', 'N', M, M, N, A, LDA,
Line 1200 Line 1412
 *  *
 *              Copy real matrix RWORK(IRVT) to complex matrix VT  *              Copy real matrix RWORK(IRVT) to complex matrix VT
 *              Overwrite VT by right singular vectors of A  *              Overwrite VT by right singular vectors of A
 *              (CWorkspace: need 3*N, prefer 2*N+N*NB)  *              CWorkspace: need   2*N [tauq, taup] + N    [work]
 *              (RWorkspace: 0)  *              CWorkspace: prefer 2*N [tauq, taup] + N*NB [work]
   *              RWorkspace: need   N [e] + N*N [RU] + N*N [RVT]
 *  *
                CALL ZLACP2( 'F', N, N, RWORK( IRVT ), N, VT, LDVT )                 CALL ZLACP2( 'F', N, N, RWORK( IRVT ), N, VT, LDVT )
                CALL ZUNMBR( 'P', 'R', 'C', N, N, N, A, LDA,                 CALL ZUNMBR( 'P', 'R', 'C', N, N, N, A, LDA,
Line 1221 Line 1434
 *  *
             IF( WNTQN ) THEN              IF( WNTQN ) THEN
 *  *
 *              Path 1t (N much larger than M, JOBZ='N')  *              Path 1t (N >> M, JOBZ='N')
 *              No singular vectors to be computed  *              No singular vectors to be computed
 *  *
                ITAU = 1                 ITAU = 1
                NWORK = ITAU + M                 NWORK = ITAU + M
 *  *
 *              Compute A=L*Q  *              Compute A=L*Q
 *              (CWorkspace: need 2*M, prefer M+M*NB)  *              CWorkspace: need   M [tau] + M    [work]
 *              (RWorkspace: 0)  *              CWorkspace: prefer M [tau] + M*NB [work]
   *              RWorkspace: need   0
 *  *
                CALL ZGELQF( M, N, A, LDA, WORK( ITAU ), WORK( NWORK ),                 CALL ZGELQF( M, N, A, LDA, WORK( ITAU ), WORK( NWORK ),
      $                      LWORK-NWORK+1, IERR )       $                      LWORK-NWORK+1, IERR )
Line 1244 Line 1458
                NWORK = ITAUP + M                 NWORK = ITAUP + M
 *  *
 *              Bidiagonalize L in A  *              Bidiagonalize L in A
 *              (CWorkspace: need 3*M, prefer 2*M+2*M*NB)  *              CWorkspace: need   2*M [tauq, taup] + M      [work]
 *              (RWorkspace: need M)  *              CWorkspace: prefer 2*M [tauq, taup] + 2*M*NB [work]
   *              RWorkspace: need   M [e]
 *  *
                CALL ZGEBRD( M, M, A, LDA, S, RWORK( IE ), WORK( ITAUQ ),                 CALL ZGEBRD( M, M, A, LDA, S, RWORK( IE ), WORK( ITAUQ ),
      $                      WORK( ITAUP ), WORK( NWORK ), LWORK-NWORK+1,       $                      WORK( ITAUP ), WORK( NWORK ), LWORK-NWORK+1,
Line 1253 Line 1468
                NRWORK = IE + M                 NRWORK = IE + M
 *  *
 *              Perform bidiagonal SVD, compute singular values only  *              Perform bidiagonal SVD, compute singular values only
 *              (CWorkspace: 0)  *              CWorkspace: need   0
 *              (RWorkspace: need BDSPAN)  *              RWorkspace: need   M [e] + BDSPAC
 *  *
                CALL DBDSDC( 'U', 'N', M, S, RWORK( IE ), DUM, 1, DUM, 1,                 CALL DBDSDC( 'U', 'N', M, S, RWORK( IE ), DUM,1,DUM,1,
      $                      DUM, IDUM, RWORK( NRWORK ), IWORK, INFO )       $                      DUM, IDUM, RWORK( NRWORK ), IWORK, INFO )
 *  *
             ELSE IF( WNTQO ) THEN              ELSE IF( WNTQO ) THEN
 *  *
 *              Path 2t (N much larger than M, JOBZ='O')  *              Path 2t (N >> M, JOBZ='O')
 *              M right singular vectors to be overwritten on A and  *              M right singular vectors to be overwritten on A and
 *              M left singular vectors to be computed in U  *              M left singular vectors to be computed in U
 *  *
Line 1271 Line 1486
 *              WORK(IVT) is M by M  *              WORK(IVT) is M by M
 *  *
                IL = IVT + LDWKVT*M                 IL = IVT + LDWKVT*M
                IF( LWORK.GE.M*N+M*M+3*M ) THEN                 IF( LWORK .GE. M*N + M*M + 3*M ) THEN
 *  *
 *                 WORK(IL) M by N  *                 WORK(IL) M by N
 *  *
Line 1282 Line 1497
 *                 WORK(IL) is M by CHUNK  *                 WORK(IL) is M by CHUNK
 *  *
                   LDWRKL = M                    LDWRKL = M
                   CHUNK = ( LWORK-M*M-3*M ) / M                    CHUNK = ( LWORK - M*M - 3*M ) / M
                END IF                 END IF
                ITAU = IL + LDWRKL*CHUNK                 ITAU = IL + LDWRKL*CHUNK
                NWORK = ITAU + M                 NWORK = ITAU + M
 *  *
 *              Compute A=L*Q  *              Compute A=L*Q
 *              (CWorkspace: need 2*M, prefer M+M*NB)  *              CWorkspace: need   M*M [VT] + M*M [L] + M [tau] + M    [work]
 *              (RWorkspace: 0)  *              CWorkspace: prefer M*M [VT] + M*M [L] + M [tau] + M*NB [work]
   *              RWorkspace: need   0
 *  *
                CALL ZGELQF( M, N, A, LDA, WORK( ITAU ), WORK( NWORK ),                 CALL ZGELQF( M, N, A, LDA, WORK( ITAU ), WORK( NWORK ),
      $                      LWORK-NWORK+1, IERR )       $                      LWORK-NWORK+1, IERR )
Line 1301 Line 1517
      $                      WORK( IL+LDWRKL ), LDWRKL )       $                      WORK( IL+LDWRKL ), LDWRKL )
 *  *
 *              Generate Q in A  *              Generate Q in A
 *              (CWorkspace: need M*M+2*M, prefer M*M+M+M*NB)  *              CWorkspace: need   M*M [VT] + M*M [L] + M [tau] + M    [work]
 *              (RWorkspace: 0)  *              CWorkspace: prefer M*M [VT] + M*M [L] + M [tau] + M*NB [work]
   *              RWorkspace: need   0
 *  *
                CALL ZUNGLQ( M, N, M, A, LDA, WORK( ITAU ),                 CALL ZUNGLQ( M, N, M, A, LDA, WORK( ITAU ),
      $                      WORK( NWORK ), LWORK-NWORK+1, IERR )       $                      WORK( NWORK ), LWORK-NWORK+1, IERR )
Line 1312 Line 1529
                NWORK = ITAUP + M                 NWORK = ITAUP + M
 *  *
 *              Bidiagonalize L in WORK(IL)  *              Bidiagonalize L in WORK(IL)
 *              (CWorkspace: need M*M+3*M, prefer M*M+2*M+2*M*NB)  *              CWorkspace: need   M*M [VT] + M*M [L] + 2*M [tauq, taup] + M      [work]
 *              (RWorkspace: need M)  *              CWorkspace: prefer M*M [VT] + M*M [L] + 2*M [tauq, taup] + 2*M*NB [work]
   *              RWorkspace: need   M [e]
 *  *
                CALL ZGEBRD( M, M, WORK( IL ), LDWRKL, S, RWORK( IE ),                 CALL ZGEBRD( M, M, WORK( IL ), LDWRKL, S, RWORK( IE ),
      $                      WORK( ITAUQ ), WORK( ITAUP ), WORK( NWORK ),       $                      WORK( ITAUQ ), WORK( ITAUP ), WORK( NWORK ),
Line 1322 Line 1540
 *              Perform bidiagonal SVD, computing left singular vectors  *              Perform bidiagonal SVD, computing left singular vectors
 *              of bidiagonal matrix in RWORK(IRU) and computing right  *              of bidiagonal matrix in RWORK(IRU) and computing right
 *              singular vectors of bidiagonal matrix in RWORK(IRVT)  *              singular vectors of bidiagonal matrix in RWORK(IRVT)
 *              (CWorkspace: need 0)  *              CWorkspace: need   0
 *              (RWorkspace: need BDSPAC)  *              RWorkspace: need   M [e] + M*M [RU] + M*M [RVT] + BDSPAC
 *  *
                IRU = IE + M                 IRU = IE + M
                IRVT = IRU + M*M                 IRVT = IRU + M*M
Line 1334 Line 1552
 *  *
 *              Copy real matrix RWORK(IRU) to complex matrix WORK(IU)  *              Copy real matrix RWORK(IRU) to complex matrix WORK(IU)
 *              Overwrite WORK(IU) by the left singular vectors of L  *              Overwrite WORK(IU) by the left singular vectors of L
 *              (CWorkspace: need N*N+3*N, prefer M*N+2*N+N*NB)  *              CWorkspace: need   M*M [VT] + M*M [L] + 2*M [tauq, taup] + M    [work]
 *              (RWorkspace: 0)  *              CWorkspace: prefer M*M [VT] + M*M [L] + 2*M [tauq, taup] + M*NB [work]
   *              RWorkspace: need   0
 *  *
                CALL ZLACP2( 'F', M, M, RWORK( IRU ), M, U, LDU )                 CALL ZLACP2( 'F', M, M, RWORK( IRU ), M, U, LDU )
                CALL ZUNMBR( 'Q', 'L', 'N', M, M, M, WORK( IL ), LDWRKL,                 CALL ZUNMBR( 'Q', 'L', 'N', M, M, M, WORK( IL ), LDWRKL,
Line 1344 Line 1563
 *  *
 *              Copy real matrix RWORK(IRVT) to complex matrix WORK(IVT)  *              Copy real matrix RWORK(IRVT) to complex matrix WORK(IVT)
 *              Overwrite WORK(IVT) by the right singular vectors of L  *              Overwrite WORK(IVT) by the right singular vectors of L
 *              (CWorkspace: need N*N+3*N, prefer M*N+2*N+N*NB)  *              CWorkspace: need   M*M [VT] + M*M [L] + 2*M [tauq, taup] + M    [work]
 *              (RWorkspace: 0)  *              CWorkspace: prefer M*M [VT] + M*M [L] + 2*M [tauq, taup] + M*NB [work]
   *              RWorkspace: need   0
 *  *
                CALL ZLACP2( 'F', M, M, RWORK( IRVT ), M, WORK( IVT ),                 CALL ZLACP2( 'F', M, M, RWORK( IRVT ), M, WORK( IVT ),
      $                      LDWKVT )       $                      LDWKVT )
Line 1355 Line 1575
 *  *
 *              Multiply right singular vectors of L in WORK(IL) by Q  *              Multiply right singular vectors of L in WORK(IL) by Q
 *              in A, storing result in WORK(IL) and copying to A  *              in A, storing result in WORK(IL) and copying to A
 *              (CWorkspace: need 2*M*M, prefer M*M+M*N))  *              CWorkspace: need   M*M [VT] + M*M [L]
 *              (RWorkspace: 0)  *              CWorkspace: prefer M*M [VT] + M*N [L]
   *              RWorkspace: need   0
 *  *
                DO 40 I = 1, N, CHUNK                 DO 40 I = 1, N, CHUNK
                   BLK = MIN( N-I+1, CHUNK )                    BLK = MIN( N-I+1, CHUNK )
Line 1369 Line 1590
 *  *
             ELSE IF( WNTQS ) THEN              ELSE IF( WNTQS ) THEN
 *  *
 *             Path 3t (N much larger than M, JOBZ='S')  *              Path 3t (N >> M, JOBZ='S')
 *             M right singular vectors to be computed in VT and  *              M right singular vectors to be computed in VT and
 *             M left singular vectors to be computed in U  *              M left singular vectors to be computed in U
 *  *
                IL = 1                 IL = 1
 *  *
Line 1382 Line 1603
                NWORK = ITAU + M                 NWORK = ITAU + M
 *  *
 *              Compute A=L*Q  *              Compute A=L*Q
 *              (CWorkspace: need 2*M, prefer M+M*NB)  *              CWorkspace: need   M*M [L] + M [tau] + M    [work]
 *              (RWorkspace: 0)  *              CWorkspace: prefer M*M [L] + M [tau] + M*NB [work]
   *              RWorkspace: need   0
 *  *
                CALL ZGELQF( M, N, A, LDA, WORK( ITAU ), WORK( NWORK ),                 CALL ZGELQF( M, N, A, LDA, WORK( ITAU ), WORK( NWORK ),
      $                      LWORK-NWORK+1, IERR )       $                      LWORK-NWORK+1, IERR )
Line 1395 Line 1617
      $                      WORK( IL+LDWRKL ), LDWRKL )       $                      WORK( IL+LDWRKL ), LDWRKL )
 *  *
 *              Generate Q in A  *              Generate Q in A
 *              (CWorkspace: need M*M+2*M, prefer M*M+M+M*NB)  *              CWorkspace: need   M*M [L] + M [tau] + M    [work]
 *              (RWorkspace: 0)  *              CWorkspace: prefer M*M [L] + M [tau] + M*NB [work]
   *              RWorkspace: need   0
 *  *
                CALL ZUNGLQ( M, N, M, A, LDA, WORK( ITAU ),                 CALL ZUNGLQ( M, N, M, A, LDA, WORK( ITAU ),
      $                      WORK( NWORK ), LWORK-NWORK+1, IERR )       $                      WORK( NWORK ), LWORK-NWORK+1, IERR )
Line 1406 Line 1629
                NWORK = ITAUP + M                 NWORK = ITAUP + M
 *  *
 *              Bidiagonalize L in WORK(IL)  *              Bidiagonalize L in WORK(IL)
 *              (CWorkspace: need M*M+3*M, prefer M*M+2*M+2*M*NB)  *              CWorkspace: need   M*M [L] + 2*M [tauq, taup] + M      [work]
 *              (RWorkspace: need M)  *              CWorkspace: prefer M*M [L] + 2*M [tauq, taup] + 2*M*NB [work]
   *              RWorkspace: need   M [e]
 *  *
                CALL ZGEBRD( M, M, WORK( IL ), LDWRKL, S, RWORK( IE ),                 CALL ZGEBRD( M, M, WORK( IL ), LDWRKL, S, RWORK( IE ),
      $                      WORK( ITAUQ ), WORK( ITAUP ), WORK( NWORK ),       $                      WORK( ITAUQ ), WORK( ITAUP ), WORK( NWORK ),
Line 1416 Line 1640
 *              Perform bidiagonal SVD, computing left singular vectors  *              Perform bidiagonal SVD, computing left singular vectors
 *              of bidiagonal matrix in RWORK(IRU) and computing right  *              of bidiagonal matrix in RWORK(IRU) and computing right
 *              singular vectors of bidiagonal matrix in RWORK(IRVT)  *              singular vectors of bidiagonal matrix in RWORK(IRVT)
 *              (CWorkspace: need 0)  *              CWorkspace: need   0
 *              (RWorkspace: need BDSPAC)  *              RWorkspace: need   M [e] + M*M [RU] + M*M [RVT] + BDSPAC
 *  *
                IRU = IE + M                 IRU = IE + M
                IRVT = IRU + M*M                 IRVT = IRU + M*M
Line 1428 Line 1652
 *  *
 *              Copy real matrix RWORK(IRU) to complex matrix U  *              Copy real matrix RWORK(IRU) to complex matrix U
 *              Overwrite U by left singular vectors of L  *              Overwrite U by left singular vectors of L
 *              (CWorkspace: need M*M+3*M, prefer M*M+2*M+M*NB)  *              CWorkspace: need   M*M [L] + 2*M [tauq, taup] + M    [work]
 *              (RWorkspace: 0)  *              CWorkspace: prefer M*M [L] + 2*M [tauq, taup] + M*NB [work]
   *              RWorkspace: need   0
 *  *
                CALL ZLACP2( 'F', M, M, RWORK( IRU ), M, U, LDU )                 CALL ZLACP2( 'F', M, M, RWORK( IRU ), M, U, LDU )
                CALL ZUNMBR( 'Q', 'L', 'N', M, M, M, WORK( IL ), LDWRKL,                 CALL ZUNMBR( 'Q', 'L', 'N', M, M, M, WORK( IL ), LDWRKL,
Line 1438 Line 1663
 *  *
 *              Copy real matrix RWORK(IRVT) to complex matrix VT  *              Copy real matrix RWORK(IRVT) to complex matrix VT
 *              Overwrite VT by left singular vectors of L  *              Overwrite VT by left singular vectors of L
 *              (CWorkspace: need M*M+3*M, prefer M*M+2*M+M*NB)  *              CWorkspace: need   M*M [L] + 2*M [tauq, taup] + M    [work]
 *              (RWorkspace: 0)  *              CWorkspace: prefer M*M [L] + 2*M [tauq, taup] + M*NB [work]
   *              RWorkspace: need   0
 *  *
                CALL ZLACP2( 'F', M, M, RWORK( IRVT ), M, VT, LDVT )                 CALL ZLACP2( 'F', M, M, RWORK( IRVT ), M, VT, LDVT )
                CALL ZUNMBR( 'P', 'R', 'C', M, M, M, WORK( IL ), LDWRKL,                 CALL ZUNMBR( 'P', 'R', 'C', M, M, M, WORK( IL ), LDWRKL,
Line 1448 Line 1674
 *  *
 *              Copy VT to WORK(IL), multiply right singular vectors of L  *              Copy VT to WORK(IL), multiply right singular vectors of L
 *              in WORK(IL) by Q in A, storing result in VT  *              in WORK(IL) by Q in A, storing result in VT
 *              (CWorkspace: need M*M)  *              CWorkspace: need   M*M [L]
 *              (RWorkspace: 0)  *              RWorkspace: need   0
 *  *
                CALL ZLACPY( 'F', M, M, VT, LDVT, WORK( IL ), LDWRKL )                 CALL ZLACPY( 'F', M, M, VT, LDVT, WORK( IL ), LDWRKL )
                CALL ZGEMM( 'N', 'N', M, N, M, CONE, WORK( IL ), LDWRKL,                 CALL ZGEMM( 'N', 'N', M, N, M, CONE, WORK( IL ), LDWRKL,
Line 1457 Line 1683
 *  *
             ELSE IF( WNTQA ) THEN              ELSE IF( WNTQA ) THEN
 *  *
 *              Path 9t (N much larger than M, JOBZ='A')  *              Path 4t (N >> M, JOBZ='A')
 *              N right singular vectors to be computed in VT and  *              N right singular vectors to be computed in VT and
 *              M left singular vectors to be computed in U  *              M left singular vectors to be computed in U
 *  *
Line 1470 Line 1696
                NWORK = ITAU + M                 NWORK = ITAU + M
 *  *
 *              Compute A=L*Q, copying result to VT  *              Compute A=L*Q, copying result to VT
 *              (CWorkspace: need 2*M, prefer M+M*NB)  *              CWorkspace: need   M*M [VT] + M [tau] + M    [work]
 *              (RWorkspace: 0)  *              CWorkspace: prefer M*M [VT] + M [tau] + M*NB [work]
   *              RWorkspace: need   0
 *  *
                CALL ZGELQF( M, N, A, LDA, WORK( ITAU ), WORK( NWORK ),                 CALL ZGELQF( M, N, A, LDA, WORK( ITAU ), WORK( NWORK ),
      $                      LWORK-NWORK+1, IERR )       $                      LWORK-NWORK+1, IERR )
                CALL ZLACPY( 'U', M, N, A, LDA, VT, LDVT )                 CALL ZLACPY( 'U', M, N, A, LDA, VT, LDVT )
 *  *
 *              Generate Q in VT  *              Generate Q in VT
 *              (CWorkspace: need M+N, prefer M+N*NB)  *              CWorkspace: need   M*M [VT] + M [tau] + N    [work]
 *              (RWorkspace: 0)  *              CWorkspace: prefer M*M [VT] + M [tau] + N*NB [work]
   *              RWorkspace: need   0
 *  *
                CALL ZUNGLQ( N, N, M, VT, LDVT, WORK( ITAU ),                 CALL ZUNGLQ( N, N, M, VT, LDVT, WORK( ITAU ),
      $                      WORK( NWORK ), LWORK-NWORK+1, IERR )       $                      WORK( NWORK ), LWORK-NWORK+1, IERR )
Line 1494 Line 1722
                NWORK = ITAUP + M                 NWORK = ITAUP + M
 *  *
 *              Bidiagonalize L in A  *              Bidiagonalize L in A
 *              (CWorkspace: need M*M+3*M, prefer M*M+2*M+2*M*NB)  *              CWorkspace: need   M*M [VT] + 2*M [tauq, taup] + M      [work]
 *              (RWorkspace: need M)  *              CWorkspace: prefer M*M [VT] + 2*M [tauq, taup] + 2*M*NB [work]
   *              RWorkspace: need   M [e]
 *  *
                CALL ZGEBRD( M, M, A, LDA, S, RWORK( IE ), WORK( ITAUQ ),                 CALL ZGEBRD( M, M, A, LDA, S, RWORK( IE ), WORK( ITAUQ ),
      $                      WORK( ITAUP ), WORK( NWORK ), LWORK-NWORK+1,       $                      WORK( ITAUP ), WORK( NWORK ), LWORK-NWORK+1,
Line 1504 Line 1733
 *              Perform bidiagonal SVD, computing left singular vectors  *              Perform bidiagonal SVD, computing left singular vectors
 *              of bidiagonal matrix in RWORK(IRU) and computing right  *              of bidiagonal matrix in RWORK(IRU) and computing right
 *              singular vectors of bidiagonal matrix in RWORK(IRVT)  *              singular vectors of bidiagonal matrix in RWORK(IRVT)
 *              (CWorkspace: need 0)  *              CWorkspace: need   0
 *              (RWorkspace: need BDSPAC)  *              RWorkspace: need   M [e] + M*M [RU] + M*M [RVT] + BDSPAC
 *  *
                IRU = IE + M                 IRU = IE + M
                IRVT = IRU + M*M                 IRVT = IRU + M*M
Line 1516 Line 1745
 *  *
 *              Copy real matrix RWORK(IRU) to complex matrix U  *              Copy real matrix RWORK(IRU) to complex matrix U
 *              Overwrite U by left singular vectors of L  *              Overwrite U by left singular vectors of L
 *              (CWorkspace: need 3*M, prefer 2*M+M*NB)  *              CWorkspace: need   M*M [VT] + 2*M [tauq, taup] + M    [work]
 *              (RWorkspace: 0)  *              CWorkspace: prefer M*M [VT] + 2*M [tauq, taup] + M*NB [work]
   *              RWorkspace: need   0
 *  *
                CALL ZLACP2( 'F', M, M, RWORK( IRU ), M, U, LDU )                 CALL ZLACP2( 'F', M, M, RWORK( IRU ), M, U, LDU )
                CALL ZUNMBR( 'Q', 'L', 'N', M, M, M, A, LDA,                 CALL ZUNMBR( 'Q', 'L', 'N', M, M, M, A, LDA,
Line 1526 Line 1756
 *  *
 *              Copy real matrix RWORK(IRVT) to complex matrix WORK(IVT)  *              Copy real matrix RWORK(IRVT) to complex matrix WORK(IVT)
 *              Overwrite WORK(IVT) by right singular vectors of L  *              Overwrite WORK(IVT) by right singular vectors of L
 *              (CWorkspace: need M*M+3*M, prefer M*M+2*M+M*NB)  *              CWorkspace: need   M*M [VT] + 2*M [tauq, taup] + M    [work]
 *              (RWorkspace: 0)  *              CWorkspace: prefer M*M [VT] + 2*M [tauq, taup] + M*NB [work]
   *              RWorkspace: need   0
 *  *
                CALL ZLACP2( 'F', M, M, RWORK( IRVT ), M, WORK( IVT ),                 CALL ZLACP2( 'F', M, M, RWORK( IRVT ), M, WORK( IVT ),
      $                      LDWKVT )       $                      LDWKVT )
Line 1537 Line 1768
 *  *
 *              Multiply right singular vectors of L in WORK(IVT) by  *              Multiply right singular vectors of L in WORK(IVT) by
 *              Q in VT, storing result in A  *              Q in VT, storing result in A
 *              (CWorkspace: need M*M)  *              CWorkspace: need   M*M [VT]
 *              (RWorkspace: 0)  *              RWorkspace: need   0
 *  *
                CALL ZGEMM( 'N', 'N', M, N, M, CONE, WORK( IVT ), LDWKVT,                 CALL ZGEMM( 'N', 'N', M, N, M, CONE, WORK( IVT ), LDWKVT,
      $                     VT, LDVT, CZERO, A, LDA )       $                     VT, LDVT, CZERO, A, LDA )
Line 1553 Line 1784
 *  *
 *           MNTHR2 <= N < MNTHR1  *           MNTHR2 <= N < MNTHR1
 *  *
 *           Path 5t (N much larger than M, but not as much as MNTHR1)  *           Path 5t (N >> M, but not as much as MNTHR1)
 *           Reduce to bidiagonal form without QR decomposition, use  *           Reduce to bidiagonal form without QR decomposition, use
 *           ZUNGBR and matrix multiplication to compute singular vectors  *           ZUNGBR and matrix multiplication to compute singular vectors
 *  *
 *  
             IE = 1              IE = 1
             NRWORK = IE + M              NRWORK = IE + M
             ITAUQ = 1              ITAUQ = 1
Line 1565 Line 1795
             NWORK = ITAUP + M              NWORK = ITAUP + M
 *  *
 *           Bidiagonalize A  *           Bidiagonalize A
 *           (CWorkspace: need 2*M+N, prefer 2*M+(M+N)*NB)  *           CWorkspace: need   2*M [tauq, taup] + N        [work]
 *           (RWorkspace: M)  *           CWorkspace: prefer 2*M [tauq, taup] + (M+N)*NB [work]
   *           RWorkspace: need   M [e]
 *  *
             CALL ZGEBRD( M, N, A, LDA, S, RWORK( IE ), WORK( ITAUQ ),              CALL ZGEBRD( M, N, A, LDA, S, RWORK( IE ), WORK( ITAUQ ),
      $                   WORK( ITAUP ), WORK( NWORK ), LWORK-NWORK+1,       $                   WORK( ITAUP ), WORK( NWORK ), LWORK-NWORK+1,
Line 1574 Line 1805
 *  *
             IF( WNTQN ) THEN              IF( WNTQN ) THEN
 *  *
   *              Path 5tn (N >> M, JOBZ='N')
 *              Compute singular values only  *              Compute singular values only
 *              (Cworkspace: 0)  *              CWorkspace: need   0
 *              (Rworkspace: need BDSPAN)  *              RWorkspace: need   M [e] + BDSPAC
 *  *
                CALL DBDSDC( 'L', 'N', M, S, RWORK( IE ), DUM, 1, DUM, 1,                 CALL DBDSDC( 'L', 'N', M, S, RWORK( IE ), DUM,1,DUM,1,
      $                      DUM, IDUM, RWORK( NRWORK ), IWORK, INFO )       $                      DUM, IDUM, RWORK( NRWORK ), IWORK, INFO )
             ELSE IF( WNTQO ) THEN              ELSE IF( WNTQO ) THEN
                IRVT = NRWORK                 IRVT = NRWORK
Line 1586 Line 1818
                NRWORK = IRU + M*M                 NRWORK = IRU + M*M
                IVT = NWORK                 IVT = NWORK
 *  *
   *              Path 5to (N >> M, JOBZ='O')
 *              Copy A to U, generate Q  *              Copy A to U, generate Q
 *              (Cworkspace: need 2*M, prefer M+M*NB)  *              CWorkspace: need   2*M [tauq, taup] + M    [work]
 *              (Rworkspace: 0)  *              CWorkspace: prefer 2*M [tauq, taup] + M*NB [work]
   *              RWorkspace: need   0
 *  *
                CALL ZLACPY( 'L', M, M, A, LDA, U, LDU )                 CALL ZLACPY( 'L', M, M, A, LDA, U, LDU )
                CALL ZUNGBR( 'Q', M, M, N, U, LDU, WORK( ITAUQ ),                 CALL ZUNGBR( 'Q', M, M, N, U, LDU, WORK( ITAUQ ),
      $                      WORK( NWORK ), LWORK-NWORK+1, IERR )       $                      WORK( NWORK ), LWORK-NWORK+1, IERR )
 *  *
 *              Generate P**H in A  *              Generate P**H in A
 *              (Cworkspace: need 2*M, prefer M+M*NB)  *              CWorkspace: need   2*M [tauq, taup] + M    [work]
 *              (Rworkspace: 0)  *              CWorkspace: prefer 2*M [tauq, taup] + M*NB [work]
   *              RWorkspace: need   0
 *  *
                CALL ZUNGBR( 'P', M, N, M, A, LDA, WORK( ITAUP ),                 CALL ZUNGBR( 'P', M, N, M, A, LDA, WORK( ITAUP ),
      $                      WORK( NWORK ), LWORK-NWORK+1, IERR )       $                      WORK( NWORK ), LWORK-NWORK+1, IERR )
 *  *
                LDWKVT = M                 LDWKVT = M
                IF( LWORK.GE.M*N+3*M ) THEN                 IF( LWORK .GE. M*N + 3*M ) THEN
 *  *
 *                 WORK( IVT ) is M by N  *                 WORK( IVT ) is M by N
 *  *
Line 1612 Line 1847
 *  *
 *                 WORK( IVT ) is M by CHUNK  *                 WORK( IVT ) is M by CHUNK
 *  *
                   CHUNK = ( LWORK-3*M ) / M                    CHUNK = ( LWORK - 3*M ) / M
                   NWORK = IVT + LDWKVT*CHUNK                    NWORK = IVT + LDWKVT*CHUNK
                END IF                 END IF
 *  *
 *              Perform bidiagonal SVD, computing left singular vectors  *              Perform bidiagonal SVD, computing left singular vectors
 *              of bidiagonal matrix in RWORK(IRU) and computing right  *              of bidiagonal matrix in RWORK(IRU) and computing right
 *              singular vectors of bidiagonal matrix in RWORK(IRVT)  *              singular vectors of bidiagonal matrix in RWORK(IRVT)
 *              (CWorkspace: need 0)  *              CWorkspace: need   0
 *              (RWorkspace: need BDSPAC)  *              RWorkspace: need   M [e] + M*M [RVT] + M*M [RU] + BDSPAC
 *  *
                CALL DBDSDC( 'L', 'I', M, S, RWORK( IE ), RWORK( IRU ),                 CALL DBDSDC( 'L', 'I', M, S, RWORK( IE ), RWORK( IRU ),
      $                      M, RWORK( IRVT ), M, DUM, IDUM,       $                      M, RWORK( IRVT ), M, DUM, IDUM,
Line 1628 Line 1863
 *  *
 *              Multiply Q in U by real matrix RWORK(IRVT)  *              Multiply Q in U by real matrix RWORK(IRVT)
 *              storing the result in WORK(IVT), copying to U  *              storing the result in WORK(IVT), copying to U
 *              (Cworkspace: need 0)  *              CWorkspace: need   2*M [tauq, taup] + M*M [VT]
 *              (Rworkspace: need 2*M*M)  *              RWorkspace: need   M [e] + M*M [RVT] + M*M [RU] + 2*M*M [rwork]
 *  *
                CALL ZLACRM( M, M, U, LDU, RWORK( IRU ), M, WORK( IVT ),                 CALL ZLACRM( M, M, U, LDU, RWORK( IRU ), M, WORK( IVT ),
      $                      LDWKVT, RWORK( NRWORK ) )       $                      LDWKVT, RWORK( NRWORK ) )
Line 1637 Line 1872
 *  *
 *              Multiply RWORK(IRVT) by P**H in A, storing the  *              Multiply RWORK(IRVT) by P**H in A, storing the
 *              result in WORK(IVT), copying to A  *              result in WORK(IVT), copying to A
 *              (CWorkspace: need M*M, prefer M*N)  *              CWorkspace: need   2*M [tauq, taup] + M*M [VT]
 *              (Rworkspace: need 2*M*M, prefer 2*M*N)  *              CWorkspace: prefer 2*M [tauq, taup] + M*N [VT]
   *              RWorkspace: need   M [e] + M*M [RVT] + 2*M*M [rwork]
   *              RWorkspace: prefer M [e] + M*M [RVT] + 2*M*N [rwork] < M + 5*M*M since N < 2*M here
 *  *
                NRWORK = IRU                 NRWORK = IRU
                DO 50 I = 1, N, CHUNK                 DO 50 I = 1, N, CHUNK
Line 1650 Line 1887
    50          CONTINUE     50          CONTINUE
             ELSE IF( WNTQS ) THEN              ELSE IF( WNTQS ) THEN
 *  *
   *              Path 5ts (N >> M, JOBZ='S')
 *              Copy A to U, generate Q  *              Copy A to U, generate Q
 *              (Cworkspace: need 2*M, prefer M+M*NB)  *              CWorkspace: need   2*M [tauq, taup] + M    [work]
 *              (Rworkspace: 0)  *              CWorkspace: prefer 2*M [tauq, taup] + M*NB [work]
   *              RWorkspace: need   0
 *  *
                CALL ZLACPY( 'L', M, M, A, LDA, U, LDU )                 CALL ZLACPY( 'L', M, M, A, LDA, U, LDU )
                CALL ZUNGBR( 'Q', M, M, N, U, LDU, WORK( ITAUQ ),                 CALL ZUNGBR( 'Q', M, M, N, U, LDU, WORK( ITAUQ ),
      $                      WORK( NWORK ), LWORK-NWORK+1, IERR )       $                      WORK( NWORK ), LWORK-NWORK+1, IERR )
 *  *
 *              Copy A to VT, generate P**H  *              Copy A to VT, generate P**H
 *              (Cworkspace: need 2*M, prefer M+M*NB)  *              CWorkspace: need   2*M [tauq, taup] + M    [work]
 *              (Rworkspace: 0)  *              CWorkspace: prefer 2*M [tauq, taup] + M*NB [work]
   *              RWorkspace: need   0
 *  *
                CALL ZLACPY( 'U', M, N, A, LDA, VT, LDVT )                 CALL ZLACPY( 'U', M, N, A, LDA, VT, LDVT )
                CALL ZUNGBR( 'P', M, N, M, VT, LDVT, WORK( ITAUP ),                 CALL ZUNGBR( 'P', M, N, M, VT, LDVT, WORK( ITAUP ),
Line 1669 Line 1909
 *              Perform bidiagonal SVD, computing left singular vectors  *              Perform bidiagonal SVD, computing left singular vectors
 *              of bidiagonal matrix in RWORK(IRU) and computing right  *              of bidiagonal matrix in RWORK(IRU) and computing right
 *              singular vectors of bidiagonal matrix in RWORK(IRVT)  *              singular vectors of bidiagonal matrix in RWORK(IRVT)
 *              (CWorkspace: need 0)  *              CWorkspace: need   0
 *              (RWorkspace: need BDSPAC)  *              RWorkspace: need   M [e] + M*M [RVT] + M*M [RU] + BDSPAC
 *  *
                IRVT = NRWORK                 IRVT = NRWORK
                IRU = IRVT + M*M                 IRU = IRVT + M*M
Line 1681 Line 1921
 *  *
 *              Multiply Q in U by real matrix RWORK(IRU), storing the  *              Multiply Q in U by real matrix RWORK(IRU), storing the
 *              result in A, copying to U  *              result in A, copying to U
 *              (CWorkspace: need 0)  *              CWorkspace: need   0
 *              (Rworkspace: need 3*M*M)  *              RWorkspace: need   M [e] + M*M [RVT] + M*M [RU] + 2*M*M [rwork]
 *  *
                CALL ZLACRM( M, M, U, LDU, RWORK( IRU ), M, A, LDA,                 CALL ZLACRM( M, M, U, LDU, RWORK( IRU ), M, A, LDA,
      $                      RWORK( NRWORK ) )       $                      RWORK( NRWORK ) )
Line 1690 Line 1930
 *  *
 *              Multiply real matrix RWORK(IRVT) by P**H in VT,  *              Multiply real matrix RWORK(IRVT) by P**H in VT,
 *              storing the result in A, copying to VT  *              storing the result in A, copying to VT
 *              (Cworkspace: need 0)  *              CWorkspace: need   0
 *              (Rworkspace: need M*M+2*M*N)  *              RWorkspace: need   M [e] + M*M [RVT] + 2*M*N [rwork] < M + 5*M*M since N < 2*M here
 *  *
                NRWORK = IRU                 NRWORK = IRU
                CALL ZLARCM( M, N, RWORK( IRVT ), M, VT, LDVT, A, LDA,                 CALL ZLARCM( M, N, RWORK( IRVT ), M, VT, LDVT, A, LDA,
Line 1699 Line 1939
                CALL ZLACPY( 'F', M, N, A, LDA, VT, LDVT )                 CALL ZLACPY( 'F', M, N, A, LDA, VT, LDVT )
             ELSE              ELSE
 *  *
   *              Path 5ta (N >> M, JOBZ='A')
 *              Copy A to U, generate Q  *              Copy A to U, generate Q
 *              (Cworkspace: need 2*M, prefer M+M*NB)  *              CWorkspace: need   2*M [tauq, taup] + M    [work]
 *              (Rworkspace: 0)  *              CWorkspace: prefer 2*M [tauq, taup] + M*NB [work]
   *              RWorkspace: need   0
 *  *
                CALL ZLACPY( 'L', M, M, A, LDA, U, LDU )                 CALL ZLACPY( 'L', M, M, A, LDA, U, LDU )
                CALL ZUNGBR( 'Q', M, M, N, U, LDU, WORK( ITAUQ ),                 CALL ZUNGBR( 'Q', M, M, N, U, LDU, WORK( ITAUQ ),
      $                      WORK( NWORK ), LWORK-NWORK+1, IERR )       $                      WORK( NWORK ), LWORK-NWORK+1, IERR )
 *  *
 *              Copy A to VT, generate P**H  *              Copy A to VT, generate P**H
 *              (Cworkspace: need 2*M, prefer M+M*NB)  *              CWorkspace: need   2*M [tauq, taup] + N    [work]
 *              (Rworkspace: 0)  *              CWorkspace: prefer 2*M [tauq, taup] + N*NB [work]
   *              RWorkspace: need   0
 *  *
                CALL ZLACPY( 'U', M, N, A, LDA, VT, LDVT )                 CALL ZLACPY( 'U', M, N, A, LDA, VT, LDVT )
                CALL ZUNGBR( 'P', N, N, M, VT, LDVT, WORK( ITAUP ),                 CALL ZUNGBR( 'P', N, N, M, VT, LDVT, WORK( ITAUP ),
Line 1718 Line 1961
 *              Perform bidiagonal SVD, computing left singular vectors  *              Perform bidiagonal SVD, computing left singular vectors
 *              of bidiagonal matrix in RWORK(IRU) and computing right  *              of bidiagonal matrix in RWORK(IRU) and computing right
 *              singular vectors of bidiagonal matrix in RWORK(IRVT)  *              singular vectors of bidiagonal matrix in RWORK(IRVT)
 *              (CWorkspace: need 0)  *              CWorkspace: need   0
 *              (RWorkspace: need BDSPAC)  *              RWorkspace: need   M [e] + M*M [RVT] + M*M [RU] + BDSPAC
 *  *
                IRVT = NRWORK                 IRVT = NRWORK
                IRU = IRVT + M*M                 IRU = IRVT + M*M
Line 1730 Line 1973
 *  *
 *              Multiply Q in U by real matrix RWORK(IRU), storing the  *              Multiply Q in U by real matrix RWORK(IRU), storing the
 *              result in A, copying to U  *              result in A, copying to U
 *              (CWorkspace: need 0)  *              CWorkspace: need   0
 *              (Rworkspace: need 3*M*M)  *              RWorkspace: need   M [e] + M*M [RVT] + M*M [RU] + 2*M*M [rwork]
 *  *
                CALL ZLACRM( M, M, U, LDU, RWORK( IRU ), M, A, LDA,                 CALL ZLACRM( M, M, U, LDU, RWORK( IRU ), M, A, LDA,
      $                      RWORK( NRWORK ) )       $                      RWORK( NRWORK ) )
Line 1739 Line 1982
 *  *
 *              Multiply real matrix RWORK(IRVT) by P**H in VT,  *              Multiply real matrix RWORK(IRVT) by P**H in VT,
 *              storing the result in A, copying to VT  *              storing the result in A, copying to VT
 *              (Cworkspace: need 0)  *              CWorkspace: need   0
 *              (Rworkspace: need M*M+2*M*N)  *              RWorkspace: need   M [e] + M*M [RVT] + 2*M*N [rwork] < M + 5*M*M since N < 2*M here
 *  *
                  NRWORK = IRU
                CALL ZLARCM( M, N, RWORK( IRVT ), M, VT, LDVT, A, LDA,                 CALL ZLARCM( M, N, RWORK( IRVT ), M, VT, LDVT, A, LDA,
      $                      RWORK( NRWORK ) )       $                      RWORK( NRWORK ) )
                CALL ZLACPY( 'F', M, N, A, LDA, VT, LDVT )                 CALL ZLACPY( 'F', M, N, A, LDA, VT, LDVT )
Line 1751 Line 1995
 *  *
 *           N .LT. MNTHR2  *           N .LT. MNTHR2
 *  *
 *           Path 6t (N greater than M, but not much larger)  *           Path 6t (N > M, but not much larger)
 *           Reduce to bidiagonal form without LQ decomposition  *           Reduce to bidiagonal form without LQ decomposition
 *           Use ZUNMBR to compute singular vectors  *           Use ZUNMBR to compute singular vectors
 *  *
Line 1762 Line 2006
             NWORK = ITAUP + M              NWORK = ITAUP + M
 *  *
 *           Bidiagonalize A  *           Bidiagonalize A
 *           (CWorkspace: need 2*M+N, prefer 2*M+(M+N)*NB)  *           CWorkspace: need   2*M [tauq, taup] + N        [work]
 *           (RWorkspace: M)  *           CWorkspace: prefer 2*M [tauq, taup] + (M+N)*NB [work]
   *           RWorkspace: need   M [e]
 *  *
             CALL ZGEBRD( M, N, A, LDA, S, RWORK( IE ), WORK( ITAUQ ),              CALL ZGEBRD( M, N, A, LDA, S, RWORK( IE ), WORK( ITAUQ ),
      $                   WORK( ITAUP ), WORK( NWORK ), LWORK-NWORK+1,       $                   WORK( ITAUP ), WORK( NWORK ), LWORK-NWORK+1,
      $                   IERR )       $                   IERR )
             IF( WNTQN ) THEN              IF( WNTQN ) THEN
 *  *
   *              Path 6tn (N > M, JOBZ='N')
 *              Compute singular values only  *              Compute singular values only
 *              (Cworkspace: 0)  *              CWorkspace: need   0
 *              (Rworkspace: need BDSPAN)  *              RWorkspace: need   M [e] + BDSPAC
 *  *
                CALL DBDSDC( 'L', 'N', M, S, RWORK( IE ), DUM, 1, DUM, 1,                 CALL DBDSDC( 'L', 'N', M, S, RWORK( IE ), DUM,1,DUM,1,
      $                      DUM, IDUM, RWORK( NRWORK ), IWORK, INFO )       $                      DUM, IDUM, RWORK( NRWORK ), IWORK, INFO )
             ELSE IF( WNTQO ) THEN              ELSE IF( WNTQO ) THEN
   *              Path 6to (N > M, JOBZ='O')
                LDWKVT = M                 LDWKVT = M
                IVT = NWORK                 IVT = NWORK
                IF( LWORK.GE.M*N+3*M ) THEN                 IF( LWORK .GE. M*N + 3*M ) THEN
 *  *
 *                 WORK( IVT ) is M by N  *                 WORK( IVT ) is M by N
 *  *
Line 1790 Line 2037
 *  *
 *                 WORK( IVT ) is M by CHUNK  *                 WORK( IVT ) is M by CHUNK
 *  *
                   CHUNK = ( LWORK-3*M ) / M                    CHUNK = ( LWORK - 3*M ) / M
                   NWORK = IVT + LDWKVT*CHUNK                    NWORK = IVT + LDWKVT*CHUNK
                END IF                 END IF
 *  *
 *              Perform bidiagonal SVD, computing left singular vectors  *              Perform bidiagonal SVD, computing left singular vectors
 *              of bidiagonal matrix in RWORK(IRU) and computing right  *              of bidiagonal matrix in RWORK(IRU) and computing right
 *              singular vectors of bidiagonal matrix in RWORK(IRVT)  *              singular vectors of bidiagonal matrix in RWORK(IRVT)
 *              (CWorkspace: need 0)  *              CWorkspace: need   0
 *              (RWorkspace: need BDSPAC)  *              RWorkspace: need   M [e] + M*M [RVT] + M*M [RU] + BDSPAC
 *  *
                IRVT = NRWORK                 IRVT = NRWORK
                IRU = IRVT + M*M                 IRU = IRVT + M*M
Line 1809 Line 2056
 *  *
 *              Copy real matrix RWORK(IRU) to complex matrix U  *              Copy real matrix RWORK(IRU) to complex matrix U
 *              Overwrite U by left singular vectors of A  *              Overwrite U by left singular vectors of A
 *              (Cworkspace: need 2*M, prefer M+M*NB)  *              CWorkspace: need   2*M [tauq, taup] + M*M [VT] + M    [work]
 *              (Rworkspace: need 0)  *              CWorkspace: prefer 2*M [tauq, taup] + M*M [VT] + M*NB [work]
   *              RWorkspace: need   M [e] + M*M [RVT] + M*M [RU]
 *  *
                CALL ZLACP2( 'F', M, M, RWORK( IRU ), M, U, LDU )                 CALL ZLACP2( 'F', M, M, RWORK( IRU ), M, U, LDU )
                CALL ZUNMBR( 'Q', 'L', 'N', M, M, N, A, LDA,                 CALL ZUNMBR( 'Q', 'L', 'N', M, M, N, A, LDA,
      $                      WORK( ITAUQ ), U, LDU, WORK( NWORK ),       $                      WORK( ITAUQ ), U, LDU, WORK( NWORK ),
      $                      LWORK-NWORK+1, IERR )       $                      LWORK-NWORK+1, IERR )
 *  *
                IF( LWORK.GE.M*N+3*M ) THEN                 IF( LWORK .GE. M*N + 3*M ) THEN
 *  *
 *              Copy real matrix RWORK(IRVT) to complex matrix WORK(IVT)  *                 Path 6to-fast
 *              Overwrite WORK(IVT) by right singular vectors of A,  *                 Copy real matrix RWORK(IRVT) to complex matrix WORK(IVT)
 *              copying to A  *                 Overwrite WORK(IVT) by right singular vectors of A,
 *              (Cworkspace: need M*N+2*M, prefer M*N+M+M*NB)  *                 copying to A
 *              (Rworkspace: need 0)  *                 CWorkspace: need   2*M [tauq, taup] + M*N [VT] + M    [work]
   *                 CWorkspace: prefer 2*M [tauq, taup] + M*N [VT] + M*NB [work]
   *                 RWorkspace: need   M [e] + M*M [RVT]
 *  *
                   CALL ZLACP2( 'F', M, M, RWORK( IRVT ), M, WORK( IVT ),                    CALL ZLACP2( 'F', M, M, RWORK( IRVT ), M, WORK( IVT ),
      $                         LDWKVT )       $                         LDWKVT )
Line 1833 Line 2083
                   CALL ZLACPY( 'F', M, N, WORK( IVT ), LDWKVT, A, LDA )                    CALL ZLACPY( 'F', M, N, WORK( IVT ), LDWKVT, A, LDA )
                ELSE                 ELSE
 *  *
   *                 Path 6to-slow
 *                 Generate P**H in A  *                 Generate P**H in A
 *                 (Cworkspace: need 2*M, prefer M+M*NB)  *                 CWorkspace: need   2*M [tauq, taup] + M*M [VT] + M    [work]
 *                 (Rworkspace: need 0)  *                 CWorkspace: prefer 2*M [tauq, taup] + M*M [VT] + M*NB [work]
   *                 RWorkspace: need   0
 *  *
                   CALL ZUNGBR( 'P', M, N, M, A, LDA, WORK( ITAUP ),                    CALL ZUNGBR( 'P', M, N, M, A, LDA, WORK( ITAUP ),
      $                         WORK( NWORK ), LWORK-NWORK+1, IERR )       $                         WORK( NWORK ), LWORK-NWORK+1, IERR )
 *  *
 *                 Multiply Q in A by real matrix RWORK(IRU), storing the  *                 Multiply Q in A by real matrix RWORK(IRU), storing the
 *                 result in WORK(IU), copying to A  *                 result in WORK(IU), copying to A
 *                 (CWorkspace: need M*M, prefer M*N)  *                 CWorkspace: need   2*M [tauq, taup] + M*M [VT]
 *                 (Rworkspace: need 3*M*M, prefer M*M+2*M*N)  *                 CWorkspace: prefer 2*M [tauq, taup] + M*N [VT]
   *                 RWorkspace: need   M [e] + M*M [RVT] + 2*M*M [rwork]
   *                 RWorkspace: prefer M [e] + M*M [RVT] + 2*M*N [rwork] < M + 5*M*M since N < 2*M here
 *  *
                   NRWORK = IRU                    NRWORK = IRU
                   DO 60 I = 1, N, CHUNK                    DO 60 I = 1, N, CHUNK
Line 1857 Line 2111
                END IF                 END IF
             ELSE IF( WNTQS ) THEN              ELSE IF( WNTQS ) THEN
 *  *
   *              Path 6ts (N > M, JOBZ='S')
 *              Perform bidiagonal SVD, computing left singular vectors  *              Perform bidiagonal SVD, computing left singular vectors
 *              of bidiagonal matrix in RWORK(IRU) and computing right  *              of bidiagonal matrix in RWORK(IRU) and computing right
 *              singular vectors of bidiagonal matrix in RWORK(IRVT)  *              singular vectors of bidiagonal matrix in RWORK(IRVT)
 *              (CWorkspace: need 0)  *              CWorkspace: need   0
 *              (RWorkspace: need BDSPAC)  *              RWorkspace: need   M [e] + M*M [RVT] + M*M [RU] + BDSPAC
 *  *
                IRVT = NRWORK                 IRVT = NRWORK
                IRU = IRVT + M*M                 IRU = IRVT + M*M
Line 1872 Line 2127
 *  *
 *              Copy real matrix RWORK(IRU) to complex matrix U  *              Copy real matrix RWORK(IRU) to complex matrix U
 *              Overwrite U by left singular vectors of A  *              Overwrite U by left singular vectors of A
 *              (CWorkspace: need 3*M, prefer 2*M+M*NB)  *              CWorkspace: need   2*M [tauq, taup] + M    [work]
 *              (RWorkspace: M*M)  *              CWorkspace: prefer 2*M [tauq, taup] + M*NB [work]
   *              RWorkspace: need   M [e] + M*M [RVT] + M*M [RU]
 *  *
                CALL ZLACP2( 'F', M, M, RWORK( IRU ), M, U, LDU )                 CALL ZLACP2( 'F', M, M, RWORK( IRU ), M, U, LDU )
                CALL ZUNMBR( 'Q', 'L', 'N', M, M, N, A, LDA,                 CALL ZUNMBR( 'Q', 'L', 'N', M, M, N, A, LDA,
Line 1882 Line 2138
 *  *
 *              Copy real matrix RWORK(IRVT) to complex matrix VT  *              Copy real matrix RWORK(IRVT) to complex matrix VT
 *              Overwrite VT by right singular vectors of A  *              Overwrite VT by right singular vectors of A
 *              (CWorkspace: need 3*M, prefer 2*M+M*NB)  *              CWorkspace: need   2*M [tauq, taup] + M    [work]
 *              (RWorkspace: M*M)  *              CWorkspace: prefer 2*M [tauq, taup] + M*NB [work]
   *              RWorkspace: need   M [e] + M*M [RVT]
 *  *
                CALL ZLASET( 'F', M, N, CZERO, CZERO, VT, LDVT )                 CALL ZLASET( 'F', M, N, CZERO, CZERO, VT, LDVT )
                CALL ZLACP2( 'F', M, M, RWORK( IRVT ), M, VT, LDVT )                 CALL ZLACP2( 'F', M, M, RWORK( IRVT ), M, VT, LDVT )
Line 1892 Line 2149
      $                      LWORK-NWORK+1, IERR )       $                      LWORK-NWORK+1, IERR )
             ELSE              ELSE
 *  *
   *              Path 6ta (N > M, JOBZ='A')
 *              Perform bidiagonal SVD, computing left singular vectors  *              Perform bidiagonal SVD, computing left singular vectors
 *              of bidiagonal matrix in RWORK(IRU) and computing right  *              of bidiagonal matrix in RWORK(IRU) and computing right
 *              singular vectors of bidiagonal matrix in RWORK(IRVT)  *              singular vectors of bidiagonal matrix in RWORK(IRVT)
 *              (CWorkspace: need 0)  *              CWorkspace: need   0
 *              (RWorkspace: need BDSPAC)  *              RWorkspace: need   M [e] + M*M [RVT] + M*M [RU] + BDSPAC
 *  *
                IRVT = NRWORK                 IRVT = NRWORK
                IRU = IRVT + M*M                 IRU = IRVT + M*M
Line 1908 Line 2166
 *  *
 *              Copy real matrix RWORK(IRU) to complex matrix U  *              Copy real matrix RWORK(IRU) to complex matrix U
 *              Overwrite U by left singular vectors of A  *              Overwrite U by left singular vectors of A
 *              (CWorkspace: need 3*M, prefer 2*M+M*NB)  *              CWorkspace: need   2*M [tauq, taup] + M    [work]
 *              (RWorkspace: M*M)  *              CWorkspace: prefer 2*M [tauq, taup] + M*NB [work]
   *              RWorkspace: need   M [e] + M*M [RVT] + M*M [RU]
 *  *
                CALL ZLACP2( 'F', M, M, RWORK( IRU ), M, U, LDU )                 CALL ZLACP2( 'F', M, M, RWORK( IRU ), M, U, LDU )
                CALL ZUNMBR( 'Q', 'L', 'N', M, M, N, A, LDA,                 CALL ZUNMBR( 'Q', 'L', 'N', M, M, N, A, LDA,
Line 1922 Line 2181
 *  *
 *              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,
Line 1954 Line 2214
 *  *
 *     Return optimal workspace in WORK(1)  *     Return optimal workspace in WORK(1)
 *  *
       WORK( 1 ) = MAXWRK        WORK( 1 ) = DROUNDUP_LWORK( MAXWRK )
 *  *
       RETURN        RETURN
 *  *

Removed from v.1.3  
changed lines
  Added in v.1.20


CVSweb interface <joel.bertrand@systella.fr>