rpl/lapack/lapack/zgesdd.f - diff

Return to zgesdd.f CVS log

Up to [local] / rpl / lapack / lapack

Diff for /rpl/lapack/lapack/zgesdd.f between versions 1.9 and 1.17

version 1.9, 2011/11/21 20:43:09	version 1.17, 2017/06/17 10:54:11
Line 2	Line 2
*	*
* =========== DOCUMENTATION ===========	* =========== DOCUMENTATION ===========
*	*
* Online html documentation available at	* Online html documentation available at
* http://www.netlib.org/lapack/explore-html/	* http://www.netlib.org/lapack/explore-html/
*	*
*> \htmlonly	*> \htmlonly
*> Download ZGESDD + dependencies	*> Download ZGESDD + dependencies
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zgesdd.f">	*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zgesdd.f">
*> [TGZ]</a>	*> [TGZ]</a>
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zgesdd.f">	*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zgesdd.f">
*> [ZIP]</a>	*> [ZIP]</a>
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zgesdd.f">	*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zgesdd.f">
*> [TXT]</a>	*> [TXT]</a>
*> \endhtmlonly	*> \endhtmlonly
*	*
* Definition:	* Definition:
* ===========	* ===========
*	*
* SUBROUTINE ZGESDD( JOBZ, M, N, A, LDA, S, U, LDU, VT, LDVT, WORK,	* SUBROUTINE ZGESDD( JOBZ, M, N, A, LDA, S, U, LDU, VT, LDVT,
* LWORK, RWORK, IWORK, INFO )	* WORK, LWORK, RWORK, IWORK, INFO )
*	*
* .. Scalar Arguments ..	* .. Scalar Arguments ..
* CHARACTER JOBZ	* CHARACTER JOBZ
* INTEGER INFO, LDA, LDU, LDVT, LWORK, M, N	* INTEGER INFO, LDA, LDU, LDVT, LWORK, M, N
Line 31	Line 31
* COMPLEX16 A( LDA, ), U( LDU, * ), VT( LDVT, * ),	* COMPLEX16 A( LDA, ), U( LDU, * ), VT( LDVT, * ),
* $ WORK( * )	* $ WORK( * )
* ..	* ..
*	*
*	*
*> \par Purpose:	*> \par Purpose:
* =============	* =============
Line 135	Line 135
*> \param[in] LDU	*> \param[in] LDU
*> \verbatim	*> \verbatim
*> LDU is INTEGER	*> LDU is INTEGER
*> The leading dimension of the array U. LDU >= 1; if	*> The leading dimension of the array U. LDU >= 1;
*> JOBZ = 'S' or 'A' or JOBZ = 'O' and M < N, LDU >= M.	*> if JOBZ = 'S' or 'A' or JOBZ = 'O' and M < N, LDU >= M.
*> \endverbatim	*> \endverbatim
*>	*>
*> \param[out] VT	*> \param[out] VT
Line 152	Line 152
*> \param[in] LDVT	*> \param[in] LDVT
*> \verbatim	*> \verbatim
*> LDVT is INTEGER	*> LDVT is INTEGER
*> The leading dimension of the array VT. LDVT >= 1; if	*> The leading dimension of the array VT. LDVT >= 1;
*> JOBZ = 'A' or JOBZ = 'O' and M >= N, LDVT >= N;	*> if JOBZ = 'A' or JOBZ = 'O' and M >= N, LDVT >= N;
*> if JOBZ = 'S', LDVT >= min(M,N).	*> if JOBZ = 'S', LDVT >= min(M,N).
*> \endverbatim	*> \endverbatim
*>	*>
Line 167	Line 167
*> \verbatim	*> \verbatim
*> LWORK is INTEGER	*> LWORK is INTEGER
*> The dimension of the array WORK. LWORK >= 1.	*> The dimension of the array WORK. LWORK >= 1.
> if JOBZ = 'N', LWORK >= 2min(M,N)+max(M,N).
*> if JOBZ = 'O',
> LWORK >= 2min(M,N)min(M,N)+2min(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	*> If LWORK = -1, a workspace query is assumed. The optimal
*> size for the WORK array is calculated and stored in WORK(1),	*> size for the WORK array is calculated and stored in WORK(1),
*> and no other work except argument checking is performed.	*> and no other work except argument checking is performed.
	*>
	*> Let mx = max(M,N) and mn = min(M,N).
	> If JOBZ = 'N', LWORK >= 2mn + mx.
	> If JOBZ = 'O', LWORK >= 2mnmn + 2mn + mx.
	> If JOBZ = 'S', LWORK >= mnmn + 3*mn.
	> If JOBZ = 'A', LWORK >= mnmn + 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	*> \endverbatim
*>	*>
*> \param[out] RWORK	*> \param[out] RWORK
*> \verbatim	*> \verbatim
*> RWORK is DOUBLE PRECISION array, dimension (MAX(1,LRWORK))	*> RWORK is DOUBLE PRECISION array, dimension (MAX(1,LRWORK))
> If JOBZ = 'N', LRWORK >= 5min(M,N).	*> Let mx = max(M,N) and mn = min(M,N).
*> Otherwise,	> If JOBZ = 'N', LRWORK >= 5mn (LAPACK <= 3.6 needs 7*mn);
> LRWORK >= min(M,N)max(5min(M,N)+7,2max(M,N)+2*min(M,N)+1)	> else if mx >> mn, LRWORK >= 5mnmn + 5mn;
	> else LRWORK >= max( 5mnmn + 5mn,
	> 2mxmn + 2mn*mn + mn ).
*> \endverbatim	*> \endverbatim
*>	*>
*> \param[out] IWORK	*> \param[out] IWORK
Line 203	Line 207
* Authors:	* Authors:
* ========	* ========
*	*
*> \author Univ. of Tennessee	*> \author Univ. of Tennessee
*> \author Univ. of California Berkeley	*> \author Univ. of California Berkeley
*> \author Univ. of Colorado Denver	*> \author Univ. of Colorado Denver
*> \author NAG Ltd.	*> \author NAG Ltd.
*	*
*> \date November 2011	*> \date June 2016
*	*
*> \ingroup complex16GEsing	*> \ingroup complex16GEsing
*	*
Line 219	Line 223
*> California at Berkeley, USA	*> California at Berkeley, USA
*>	*>
* =====================================================================	* =====================================================================
SUBROUTINE ZGESDD( JOBZ, M, N, A, LDA, S, U, LDU, VT, LDVT, WORK,	SUBROUTINE ZGESDD( JOBZ, M, N, A, LDA, S, U, LDU, VT, LDVT,
$ LWORK, RWORK, IWORK, INFO )	$ WORK, LWORK, RWORK, IWORK, INFO )
	implicit none
*	*
* -- LAPACK driver routine (version 3.4.0) --	* -- LAPACK driver routine (version 3.7.0) --
* -- LAPACK is a software package provided by Univ. of Tennessee, --	* -- LAPACK is a software package provided by Univ. of Tennessee, --
* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--	* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
* November 2011	* June 2016
*	*
* .. Scalar Arguments ..	* .. Scalar Arguments ..
CHARACTER JOBZ	CHARACTER JOBZ
Line 241	Line 246
* =====================================================================	* =====================================================================
*	*
* .. 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 250	Line 253
PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 )	PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 )
* ..	* ..
* .. Local Scalars ..	* .. Local Scalars ..
LOGICAL WNTQA, WNTQAS, WNTQN, WNTQO, WNTQS	LOGICAL LQUERY, WNTQA, WNTQAS, WNTQN, WNTQO, WNTQS
INTEGER BLK, CHUNK, I, IE, IERR, IL, IR, IRU, IRVT,	INTEGER BLK, CHUNK, I, IE, IERR, IL, IR, IRU, IRVT,
$ ISCL, ITAU, ITAUP, ITAUQ, IU, IVT, LDWKVT,	$ ISCL, ITAU, ITAUP, ITAUQ, IU, IVT, LDWKVT,
$ LDWRKL, LDWRKR, LDWRKU, MAXWRK, MINMN, MINWRK,	$ LDWRKL, LDWRKR, LDWRKU, MAXWRK, MINMN, MINWRK,
$ MNTHR1, MNTHR2, NRWORK, NWORK, WRKBL	$ MNTHR1, MNTHR2, NRWORK, NWORK, WRKBL
	INTEGER LWORK_ZGEBRD_MN, LWORK_ZGEBRD_MM,
	$ LWORK_ZGEBRD_NN, LWORK_ZGELQF_MN,
	$ LWORK_ZGEQRF_MN,
	$ LWORK_ZUNGBR_P_MN, LWORK_ZUNGBR_P_NN,
	$ LWORK_ZUNGBR_Q_MN, LWORK_ZUNGBR_Q_MM,
	$ LWORK_ZUNGLQ_MN, LWORK_ZUNGLQ_NN,
	$ LWORK_ZUNGQR_MM, LWORK_ZUNGQR_MN,
	$ LWORK_ZUNMBR_PRC_MM, LWORK_ZUNMBR_QLN_MM,
	$ LWORK_ZUNMBR_PRC_MN, LWORK_ZUNMBR_QLN_MN,
	$ LWORK_ZUNMBR_PRC_NN, LWORK_ZUNMBR_QLN_NN
DOUBLE PRECISION ANRM, BIGNUM, EPS, SMLNUM	DOUBLE PRECISION ANRM, BIGNUM, EPS, SMLNUM
* ..	* ..
* .. Local Arrays ..	* .. Local Arrays ..
INTEGER IDUM( 1 )	INTEGER IDUM( 1 )
DOUBLE PRECISION DUM( 1 )	DOUBLE PRECISION DUM( 1 )
	COMPLEX*16 CDUM( 1 )
* ..	* ..
* .. External Subroutines ..	* .. External Subroutines ..
EXTERNAL DBDSDC, DLASCL, XERBLA, ZGEBRD, ZGELQF, ZGEMM,	EXTERNAL DBDSDC, DLASCL, XERBLA, ZGEBRD, ZGELQF, ZGEMM,
Line 268	Line 282
* ..	* ..
* .. External Functions ..	* .. External Functions ..
LOGICAL LSAME	LOGICAL LSAME
INTEGER ILAENV
DOUBLE PRECISION DLAMCH, ZLANGE	DOUBLE PRECISION DLAMCH, ZLANGE
EXTERNAL LSAME, ILAENV, DLAMCH, ZLANGE	EXTERNAL LSAME, DLAMCH, ZLANGE
* ..	* ..
* .. Intrinsic Functions ..	* .. Intrinsic Functions ..
INTRINSIC INT, MAX, MIN, SQRT	INTRINSIC INT, MAX, MIN, SQRT
Line 279	Line 292
*	*
* Test the input arguments	* Test the input arguments
*	*
INFO = 0	INFO = 0
MINMN = MIN( M, N )	MINMN = MIN( M, N )
MNTHR1 = INT( MINMN*17.0D0 / 9.0D0 )	MNTHR1 = INT( MINMN*17.0D0 / 9.0D0 )
MNTHR2 = INT( MINMN*5.0D0 / 3.0D0 )	MNTHR2 = INT( MINMN*5.0D0 / 3.0D0 )
WNTQA = LSAME( JOBZ, 'A' )	WNTQA = LSAME( JOBZ, 'A' )
WNTQS = LSAME( JOBZ, 'S' )	WNTQS = LSAME( JOBZ, 'S' )
WNTQAS = WNTQA .OR. WNTQS	WNTQAS = WNTQA .OR. WNTQS
WNTQO = LSAME( JOBZ, 'O' )	WNTQO = LSAME( JOBZ, 'O' )
WNTQN = LSAME( JOBZ, 'N' )	WNTQN = LSAME( JOBZ, 'N' )
	LQUERY = ( LWORK.EQ.-1 )
MINWRK = 1	MINWRK = 1
MAXWRK = 1	MAXWRK = 1
*	*
Line 309	Line 323
END IF	END IF
*	*
* Compute workspace	* Compute workspace
* (Note: Comments in the code beginning "Workspace:" describe the	* Note: Comments in the code beginning "Workspace:" describe the
* minimal amount of workspace needed at that point in the code,	* minimal amount of workspace allocated at that point in the code,
* as well as the preferred amount for good performance.	* as well as the preferred amount for good performance.
* CWorkspace refers to complex workspace, and RWorkspace to	* CWorkspace refers to complex workspace, and RWorkspace to
* real workspace. NB refers to the optimal block size for the	* real workspace. NB refers to the optimal block size for the
* immediately following subroutine, as returned by ILAENV.)	* immediately following subroutine, as returned by ILAENV.)
*	*
IF( INFO.EQ.0 .AND. M.GT.0 .AND. N.GT.0 ) THEN	IF( INFO.EQ.0 ) THEN
IF( M.GE.N ) THEN	MINWRK = 1
	MAXWRK = 1
	IF( M.GE.N .AND. MINMN.GT.0 ) THEN
*	*
* There is no complex work space needed for bidiagonal SVD	* There is no complex work space needed for bidiagonal SVD
* The real work space needed for bidiagonal SVD is BDSPAC	* The real work space needed for bidiagonal SVD (dbdsdc) is
* for computing singular values and singular vectors; BDSPAN	* BDSPAC = 3NN + 4*N for singular values and vectors;
* for computing singular values only.	* BDSPAC = 4*N for singular values only;
* BDSPAC = 5NN + 7*N	* not including e, RU, and RVT matrices.
* BDSPAN = MAX(7N+4, 3N+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, 2N+2N*
$ 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, 2N+2N*	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, 2N+N
$ ILAENV( 1, 'ZUNMBR', 'QLN', N, N, N, -1 ) )
WRKBL = MAX( WRKBL, 2N+N
$ ILAENV( 1, 'ZUNMBR', 'PRC', N, N, N, -1 ) )
MAXWRK = MN + NN + WRKBL	MAXWRK = MN + NN + WRKBL
MINWRK = 2NN + 3*N	MINWRK = 2NN + 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, 2N+2N*	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, 2N+N
$ ILAENV( 1, 'ZUNMBR', 'QLN', N, N, N, -1 ) )
WRKBL = MAX( WRKBL, 2N+N
$ ILAENV( 1, 'ZUNMBR', 'PRC', N, N, N, -1 ) )
MAXWRK = N*N + WRKBL	MAXWRK = N*N + WRKBL
MINWRK = NN + 3N	MINWRK = NN + 3N
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, 2N+2N*	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, 2N+N
$ ILAENV( 1, 'ZUNMBR', 'QLN', N, N, N, -1 ) )
WRKBL = MAX( WRKBL, 2N+N
$ ILAENV( 1, 'ZUNMBR', 'PRC', N, N, N, -1 ) )
MAXWRK = N*N + WRKBL	MAXWRK = N*N + WRKBL
MINWRK = NN + 2N + M	MINWRK = NN + MAX( 3N, 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 = 2N + ( 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, 2N+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, 2N+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, 2N+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, 2N+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, 2N+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, 2N+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 = 2N + ( 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, 2N+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, 2N+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, 2N+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, 2N+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, 2N+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, 2N+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 = 3MM + 4*M for singular values and vectors;
* for computing singular values only.	* BDSPAC = 4*M for singular values only;
* BDSPAC = 5MM + 7*M	* not including e, RU, and RVT matrices.
* BDSPAN = MAX(7M+4, 3M+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, 2M+2M*
$ 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, 2M+2M*	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, 2M+M
$ ILAENV( 1, 'ZUNMBR', 'PRC', M, M, M, -1 ) )
WRKBL = MAX( WRKBL, 2M+M
$ ILAENV( 1, 'ZUNMBR', 'QLN', M, M, M, -1 ) )
MAXWRK = MN + MM + WRKBL	MAXWRK = MN + MM + WRKBL
MINWRK = 2MM + 3*M	MINWRK = 2MM + 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, 2M+2M*	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, 2M+M
$ ILAENV( 1, 'ZUNMBR', 'PRC', M, M, M, -1 ) )
WRKBL = MAX( WRKBL, 2M+M
$ ILAENV( 1, 'ZUNMBR', 'QLN', M, M, M, -1 ) )
MAXWRK = M*M + WRKBL	MAXWRK = M*M + WRKBL
MINWRK = MM + 3M	MINWRK = MM + 3M
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, 2M+2M*	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, 2M+M
$ ILAENV( 1, 'ZUNMBR', 'PRC', M, M, M, -1 ) )
WRKBL = MAX( WRKBL, 2M+M
$ ILAENV( 1, 'ZUNMBR', 'QLN', M, M, M, -1 ) )
MAXWRK = M*M + WRKBL	MAXWRK = M*M + WRKBL
MINWRK = MM + 2M + N	MINWRK = MM + MAX( 3M, 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 = 2M + ( 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, 2M+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, 2M+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, 2M+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, 2M+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, 2M+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, 2M+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 = 2M + ( 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, 2M+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, 2M+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, 2M+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, 2M+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, 2M+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, 2M+M	MAXWRK = MAX( MAXWRK, 2*M + LWORK_ZUNMBR_PRC_NN )
$ ILAENV( 1, 'ZUNGBR', 'QLN', M, M, N, -1 ) )
END IF	END IF
END IF	END IF
END IF	END IF
Line 554	Line 620
END IF	END IF
IF( INFO.EQ.0 ) THEN	IF( INFO.EQ.0 ) THEN
WORK( 1 ) = MAXWRK	WORK( 1 ) = MAXWRK
IF( LWORK.LT.MINWRK .AND. LWORK.NE.LQUERV )	IF( LWORK.LT.MINWRK .AND. .NOT. LQUERY ) THEN
$ INFO = -13	INFO = -12
	END IF
END IF	END IF
*	*
* Quick returns
*
IF( INFO.NE.0 ) THEN	IF( INFO.NE.0 ) THEN
CALL XERBLA( 'ZGESDD', -INFO )	CALL XERBLA( 'ZGESDD', -INFO )
RETURN	RETURN
	ELSE IF( LQUERY ) THEN
	RETURN
END IF	END IF
IF( LWORK.EQ.LQUERV )	*
$ RETURN	* Quick return if possible
	*
IF( M.EQ.0 .OR. N.EQ.0 ) THEN	IF( M.EQ.0 .OR. N.EQ.0 ) THEN
RETURN	RETURN
END IF	END IF
Line 598	Line 666
*	*
IF( WNTQN ) THEN	IF( WNTQN ) THEN
*	*
* Path 1 (M much larger than N, JOBZ='N')	* Path 1 (M >> N, JOBZ='N')
* No singular vectors to be computed	* No singular vectors to be computed
*	*
ITAU = 1	ITAU = 1
NWORK = ITAU + N	NWORK = ITAU + N
*	*
* Compute A=Q*R	* Compute A=Q*R
* (CWorkspace: need 2N, prefer N+NNB)	* 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 621	Line 690
NWORK = ITAUP + N	NWORK = ITAUP + N
*	*
* Bidiagonalize R in A	* Bidiagonalize R in A
* (CWorkspace: need 3N, prefer 2N+2NNB)	* CWorkspace: need 2*N [tauq, taup] + N [work]
* (RWorkspace: need N)	* CWorkspace: prefer 2N [tauq, taup] + 2N*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 630	Line 700
NRWORK = IE + N	NRWORK = IE + N
*	*
* Perform bidiagonal SVD, compute singular values only	* Perform bidiagonal SVD, compute singular values only
* (CWorkspace: 0)	* CWorkspace: need 0
* (RWorkspace: need BDSPAN)	* RWorkspace: need N [e] + BDSPAC
*	*
CALL DBDSDC( 'U', 'N', N, S, RWORK( IE ), DUM, 1, DUM, 1,	CALL DBDSDC( 'U', 'N', N, S, RWORK( IE ), DUM,1,DUM,1,
$ DUM, IDUM, RWORK( NRWORK ), IWORK, INFO )	$ DUM, IDUM, RWORK( NRWORK ), IWORK, INFO )
*	*
ELSE IF( WNTQO ) THEN	ELSE IF( WNTQO ) THEN
*	*
* Path 2 (M much larger than N, JOBZ='O')	* Path 2 (M >> N, JOBZ='O')
* N left singular vectors to be overwritten on A and	* N left singular vectors to be overwritten on A and
* N right singular vectors to be computed in VT	* N right singular vectors to be computed in VT
*	*
Line 648	Line 718
*	*
LDWRKU = N	LDWRKU = N
IR = IU + LDWRKU*N	IR = IU + LDWRKU*N
IF( LWORK.GE.MN+NN+3*N ) THEN	IF( LWORK .GE. MN + NN + 3*N ) THEN
*	*
* WORK(IR) is M by N	* WORK(IR) is M by N
*	*
LDWRKR = M	LDWRKR = M
ELSE	ELSE
LDWRKR = ( LWORK-NN-3N ) / N	LDWRKR = ( LWORK - NN - 3N ) / 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 NN+2N, prefer MN+N+NNB)	* CWorkspace: need NN [U] + NN [R] + N [tau] + N [work]
* (RWorkspace: 0)	* CWorkspace: prefer NN [U] + NN [R] + N [tau] + N*NB [work]
	* RWorkspace: need 0
*	*
CALL ZGEQRF( M, N, A, LDA, WORK( ITAU ), WORK( NWORK ),	CALL ZGEQRF( M, N, A, LDA, WORK( ITAU ), WORK( NWORK ),
$ LWORK-NWORK+1, IERR )	$ LWORK-NWORK+1, IERR )
Line 673	Line 744
$ LDWRKR )	$ LDWRKR )
*	*
* Generate Q in A	* Generate Q in A
* (CWorkspace: need 2N, prefer N+NNB)	* CWorkspace: need NN [U] + NN [R] + N [tau] + N [work]
* (RWorkspace: 0)	* CWorkspace: prefer NN [U] + NN [R] + N [tau] + N*NB [work]
	* RWorkspace: need 0
*	*
CALL ZUNGQR( M, N, N, A, LDA, WORK( ITAU ),	CALL ZUNGQR( M, N, N, A, LDA, WORK( ITAU ),
$ WORK( NWORK ), LWORK-NWORK+1, IERR )	$ WORK( NWORK ), LWORK-NWORK+1, IERR )
Line 684	Line 756
NWORK = ITAUP + N	NWORK = ITAUP + N
*	*
* Bidiagonalize R in WORK(IR)	* Bidiagonalize R in WORK(IR)
* (CWorkspace: need NN+3N, prefer MN+2N+2NNB)	* CWorkspace: need NN [U] + NN [R] + 2*N [tauq, taup] + N [work]
* (RWorkspace: need N)	* CWorkspace: prefer NN [U] + NN [R] + 2N [tauq, taup] + 2N*NB [work]
	* RWorkspace: need N [e]
*	*
CALL ZGEBRD( N, N, WORK( IR ), LDWRKR, S, RWORK( IE ),	CALL ZGEBRD( N, N, WORK( IR ), LDWRKR, S, RWORK( IE ),
$ WORK( ITAUQ ), WORK( ITAUP ), WORK( NWORK ),	$ WORK( ITAUQ ), WORK( ITAUP ), WORK( NWORK ),
Line 694	Line 767
* Perform bidiagonal SVD, computing left singular vectors	* Perform bidiagonal SVD, computing left singular vectors
* of R in WORK(IRU) and computing right singular vectors	* of R in WORK(IRU) and computing right singular vectors
* of R in WORK(IRVT)	* of R in WORK(IRVT)
* (CWorkspace: need 0)	* CWorkspace: need 0
* (RWorkspace: need BDSPAC)	* RWorkspace: need N [e] + NN [RU] + NN [RVT] + BDSPAC
*	*
IRU = IE + N	IRU = IE + N
IRVT = IRU + N*N	IRVT = IRU + N*N
Line 706	Line 779
*	*
* Copy real matrix RWORK(IRU) to complex matrix WORK(IU)	* Copy real matrix RWORK(IRU) to complex matrix WORK(IU)
* Overwrite WORK(IU) by the left singular vectors of R	* Overwrite WORK(IU) by the left singular vectors of R
* (CWorkspace: need 2NN+3N, prefer MN+NN+2N+N*NB)	* CWorkspace: need NN [U] + NN [R] + 2*N [tauq, taup] + N [work]
* (RWorkspace: 0)	* CWorkspace: prefer NN [U] + NN [R] + 2N [tauq, taup] + NNB [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 717	Line 791
*	*
* Copy real matrix RWORK(IRVT) to complex matrix VT	* Copy real matrix RWORK(IRVT) to complex matrix VT
* Overwrite VT by the right singular vectors of R	* Overwrite VT by the right singular vectors of R
* (CWorkspace: need NN+3N, prefer MN+2N+N*NB)	* CWorkspace: need NN [U] + NN [R] + 2*N [tauq, taup] + N [work]
* (RWorkspace: 0)	* CWorkspace: prefer NN [U] + NN [R] + 2N [tauq, taup] + NNB [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 727	Line 802
*	*
* Multiply Q in A by left singular vectors of R in	* Multiply Q in A by left singular vectors of R in
* WORK(IU), storing result in WORK(IR) and copying to A	* WORK(IU), storing result in WORK(IR) and copying to A
* (CWorkspace: need 2NN, prefer NN+MN)	* CWorkspace: need NN [U] + NN [R]
* (RWorkspace: 0)	* CWorkspace: prefer NN [U] + MN [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 741	Line 817
*	*
ELSE IF( WNTQS ) THEN	ELSE IF( WNTQS ) THEN
*	*
* Path 3 (M much larger than N, JOBZ='S')	* Path 3 (M >> N, JOBZ='S')
* N left singular vectors to be computed in U and	* N left singular vectors to be computed in U and
* N right singular vectors to be computed in VT	* N right singular vectors to be computed in VT
*	*
Line 754	Line 830
NWORK = ITAU + N	NWORK = ITAU + N
*	*
* Compute A=Q*R	* Compute A=Q*R
* (CWorkspace: need NN+2N, prefer NN+N+NNB)	* CWorkspace: need N*N [R] + N [tau] + N [work]
* (RWorkspace: 0)	* CWorkspace: prefer NN [R] + N [tau] + NNB [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 767	Line 844
$ LDWRKR )	$ LDWRKR )
*	*
* Generate Q in A	* Generate Q in A
* (CWorkspace: need 2N, prefer N+NNB)	* CWorkspace: need N*N [R] + N [tau] + N [work]
* (RWorkspace: 0)	* CWorkspace: prefer NN [R] + N [tau] + NNB [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 778	Line 856
NWORK = ITAUP + N	NWORK = ITAUP + N
*	*
* Bidiagonalize R in WORK(IR)	* Bidiagonalize R in WORK(IR)
* (CWorkspace: need NN+3N, prefer NN+2N+2NNB)	* CWorkspace: need NN [R] + 2N [tauq, taup] + N [work]
* (RWorkspace: need N)	* CWorkspace: prefer NN [R] + 2N [tauq, taup] + 2NNB [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 788	Line 867
* Perform bidiagonal SVD, computing left singular vectors	* Perform bidiagonal SVD, computing left singular vectors
* of bidiagonal matrix in RWORK(IRU) and computing right	* of bidiagonal matrix in RWORK(IRU) and computing right
* singular vectors of bidiagonal matrix in RWORK(IRVT)	* singular vectors of bidiagonal matrix in RWORK(IRVT)
* (CWorkspace: need 0)	* CWorkspace: need 0
* (RWorkspace: need BDSPAC)	* RWorkspace: need N [e] + NN [RU] + NN [RVT] + BDSPAC
*	*
IRU = IE + N	IRU = IE + N
IRVT = IRU + N*N	IRVT = IRU + N*N
Line 800	Line 879
*	*
* Copy real matrix RWORK(IRU) to complex matrix U	* Copy real matrix RWORK(IRU) to complex matrix U
* Overwrite U by left singular vectors of R	* Overwrite U by left singular vectors of R
* (CWorkspace: need NN+3N, prefer NN+2N+N*NB)	* CWorkspace: need NN [R] + 2N [tauq, taup] + N [work]
* (RWorkspace: 0)	* CWorkspace: prefer NN [R] + 2N [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 810	Line 890
*	*
* Copy real matrix RWORK(IRVT) to complex matrix VT	* Copy real matrix RWORK(IRVT) to complex matrix VT
* Overwrite VT by right singular vectors of R	* Overwrite VT by right singular vectors of R
* (CWorkspace: need NN+3N, prefer NN+2N+N*NB)	* CWorkspace: need NN [R] + 2N [tauq, taup] + N [work]
* (RWorkspace: 0)	* CWorkspace: prefer NN [R] + 2N [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 820	Line 901
*	*
* Multiply Q in A by left singular vectors of R in	* Multiply Q in A by left singular vectors of R in
* WORK(IR), storing result in U	* WORK(IR), storing result in U
* (CWorkspace: need N*N)	* CWorkspace: need N*N [R]
* (RWorkspace: 0)	* RWorkspace: need 0
*	*
CALL ZLACPY( 'F', N, N, U, LDU, WORK( IR ), LDWRKR )	CALL ZLACPY( 'F', N, N, U, LDU, WORK( IR ), LDWRKR )
CALL ZGEMM( 'N', 'N', M, N, N, CONE, A, LDA, WORK( IR ),	CALL ZGEMM( 'N', 'N', M, N, N, CONE, A, LDA, WORK( IR ),
Line 829	Line 910
*	*
ELSE IF( WNTQA ) THEN	ELSE IF( WNTQA ) THEN
*	*
* Path 4 (M much larger than N, JOBZ='A')	* Path 4 (M >> N, JOBZ='A')
* M left singular vectors to be computed in U and	* M left singular vectors to be computed in U and
* N right singular vectors to be computed in VT	* N right singular vectors to be computed in VT
*	*
Line 842	Line 923
NWORK = ITAU + N	NWORK = ITAU + N
*	*
* Compute A=Q*R, copying result to U	* Compute A=Q*R, copying result to U
* (CWorkspace: need 2N, prefer N+NNB)	* CWorkspace: need N*N [U] + N [tau] + N [work]
* (RWorkspace: 0)	* CWorkspace: prefer NN [U] + N [tau] + NNB [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 NN [U] + N [tau] + MNB [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 866	Line 949
NWORK = ITAUP + N	NWORK = ITAUP + N
*	*
* Bidiagonalize R in A	* Bidiagonalize R in A
* (CWorkspace: need 3N, prefer 2N+2NNB)	* CWorkspace: need NN [U] + 2N [tauq, taup] + N [work]
* (RWorkspace: need N)	* CWorkspace: prefer NN [U] + 2N [tauq, taup] + 2NNB [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 879	Line 963
* Perform bidiagonal SVD, computing left singular vectors	* Perform bidiagonal SVD, computing left singular vectors
* of bidiagonal matrix in RWORK(IRU) and computing right	* of bidiagonal matrix in RWORK(IRU) and computing right
* singular vectors of bidiagonal matrix in RWORK(IRVT)	* singular vectors of bidiagonal matrix in RWORK(IRVT)
* (CWorkspace: need 0)	* CWorkspace: need 0
* (RWorkspace: need BDSPAC)	* RWorkspace: need N [e] + NN [RU] + NN [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 888	Line 972
*	*
* Copy real matrix RWORK(IRU) to complex matrix WORK(IU)	* Copy real matrix RWORK(IRU) to complex matrix WORK(IU)
* Overwrite WORK(IU) by left singular vectors of R	* Overwrite WORK(IU) by left singular vectors of R
* (CWorkspace: need NN+3N, prefer NN+2N+N*NB)	* CWorkspace: need NN [U] + 2N [tauq, taup] + N [work]
* (RWorkspace: 0)	* CWorkspace: prefer NN [U] + 2N [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 899	Line 984
*	*
* Copy real matrix RWORK(IRVT) to complex matrix VT	* Copy real matrix RWORK(IRVT) to complex matrix VT
* Overwrite VT by right singular vectors of R	* Overwrite VT by right singular vectors of R
* (CWorkspace: need 3N, prefer 2N+N*NB)	* CWorkspace: need NN [U] + 2N [tauq, taup] + N [work]
* (RWorkspace: 0)	* CWorkspace: prefer NN [U] + 2N [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 909	Line 995
*	*
* Multiply Q in U by left singular vectors of R in	* Multiply Q in U by left singular vectors of R in
* WORK(IU), storing result in A	* WORK(IU), storing result in A
* (CWorkspace: need N*N)	* CWorkspace: need N*N [U]
* (RWorkspace: 0)	* RWorkspace: need 0
*	*
CALL ZGEMM( 'N', 'N', M, N, N, CONE, U, LDU, WORK( IU ),	CALL ZGEMM( 'N', 'N', M, N, N, CONE, U, LDU, WORK( IU ),
$ LDWRKU, CZERO, A, LDA )	$ LDWRKU, CZERO, A, LDA )
Line 925	Line 1011
*	*
* MNTHR2 <= M < MNTHR1	* MNTHR2 <= M < MNTHR1
*	*
* Path 5 (M much larger than N, but not as much as MNTHR1)	* Path 5 (M >> N, but not as much as MNTHR1)
* Reduce to bidiagonal form without QR decomposition, use	* Reduce to bidiagonal form without QR decomposition, use
* ZUNGBR and matrix multiplication to compute singular vectors	* ZUNGBR and matrix multiplication to compute singular vectors
*	*
Line 936	Line 1022
NWORK = ITAUP + N	NWORK = ITAUP + N
*	*
* Bidiagonalize A	* Bidiagonalize A
* (CWorkspace: need 2N+M, prefer 2N+(M+N)*NB)	* CWorkspace: need 2*N [tauq, taup] + M [work]
* (RWorkspace: need N)	* CWorkspace: prefer 2N [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 956	Line 1044
IRVT = IRU + N*N	IRVT = IRU + N*N
NRWORK = IRVT + N*N	NRWORK = IRVT + N*N
*	*
	* Path 5o (M >> N, JOBZ='O')
* Copy A to VT, generate P**H	* Copy A to VT, generate P**H
* (Cworkspace: need 2N, prefer N+NNB)	* CWorkspace: need 2*N [tauq, taup] + N [work]
* (Rworkspace: 0)	* CWorkspace: prefer 2N [tauq, taup] + NNB [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 2N, prefer N+NNB)	* CWorkspace: need 2*N [tauq, taup] + N [work]
* (RWorkspace: 0)	* CWorkspace: prefer 2N [tauq, taup] + NNB [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.MN+3N ) THEN	IF( LWORK .GE. MN + 3N ) THEN
*	*
* WORK( IU ) is M by N	* WORK( IU ) is M by N
*	*
Line 980	Line 1071
*	*
* WORK(IU) is LDWRKU by N	* WORK(IU) is LDWRKU by N
*	*
LDWRKU = ( LWORK-3*N ) / N	LDWRKU = ( LWORK - 3*N ) / N
END IF	END IF
NWORK = IU + LDWRKU*N	NWORK = IU + LDWRKU*N
*	*
* Perform bidiagonal SVD, computing left singular vectors	* Perform bidiagonal SVD, computing left singular vectors
* of bidiagonal matrix in RWORK(IRU) and computing right	* of bidiagonal matrix in RWORK(IRU) and computing right
* singular vectors of bidiagonal matrix in RWORK(IRVT)	* singular vectors of bidiagonal matrix in RWORK(IRVT)
* (CWorkspace: need 0)	* CWorkspace: need 0
* (RWorkspace: need BDSPAC)	* RWorkspace: need N [e] + NN [RU] + NN [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 996	Line 1087
*	*
* Multiply real matrix RWORK(IRVT) by P**H in VT,	* Multiply real matrix RWORK(IRVT) by P**H in VT,
* storing the result in WORK(IU), copying to VT	* storing the result in WORK(IU), copying to VT
* (Cworkspace: need 0)	* CWorkspace: need 2N [tauq, taup] + NN [U]
* (Rworkspace: need 3NN)	* RWorkspace: need N [e] + NN [RU] + NN [RVT] + 2NN [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 1005	Line 1096
*	*
* Multiply Q in A by real matrix RWORK(IRU), storing the	* Multiply Q in A by real matrix RWORK(IRU), storing the
* result in WORK(IU), copying to A	* result in WORK(IU), copying to A
* (CWorkspace: need NN, prefer MN)	* CWorkspace: need 2N [tauq, taup] + NN [U]
* (Rworkspace: need 3NN, prefer NN+2M*N)	* CWorkspace: prefer 2N [tauq, taup] + MN [U]
	* RWorkspace: need N [e] + NN [RU] + 2N*N [rwork]
	* RWorkspace: prefer N [e] + NN [RU] + 2MN [rwork] < N + 5NN since M < 2N here
*	*
NRWORK = IRVT	NRWORK = IRVT
DO 20 I = 1, M, LDWRKU	DO 20 I = 1, M, LDWRKU
Line 1019	Line 1112
*	*
ELSE IF( WNTQS ) THEN	ELSE IF( WNTQS ) THEN
*	*
	* Path 5s (M >> N, JOBZ='S')
* Copy A to VT, generate P**H	* Copy A to VT, generate P**H
* (Cworkspace: need 2N, prefer N+NNB)	* CWorkspace: need 2*N [tauq, taup] + N [work]
* (Rworkspace: 0)	* CWorkspace: prefer 2N [tauq, taup] + NNB [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 2N, prefer N+NNB)	* CWorkspace: need 2*N [tauq, taup] + N [work]
* (Rworkspace: 0)	* CWorkspace: prefer 2N [tauq, taup] + NNB [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 1038	Line 1134
* Perform bidiagonal SVD, computing left singular vectors	* Perform bidiagonal SVD, computing left singular vectors
* of bidiagonal matrix in RWORK(IRU) and computing right	* of bidiagonal matrix in RWORK(IRU) and computing right
* singular vectors of bidiagonal matrix in RWORK(IRVT)	* singular vectors of bidiagonal matrix in RWORK(IRVT)
* (CWorkspace: need 0)	* CWorkspace: need 0
* (RWorkspace: need BDSPAC)	* RWorkspace: need N [e] + NN [RU] + NN [RVT] + BDSPAC
*	*
IRU = NRWORK	IRU = NRWORK
IRVT = IRU + N*N	IRVT = IRU + N*N
Line 1050	Line 1146
*	*
* Multiply real matrix RWORK(IRVT) by P**H in VT,	* Multiply real matrix RWORK(IRVT) by P**H in VT,
* storing the result in A, copying to VT	* storing the result in A, copying to VT
* (Cworkspace: need 0)	* CWorkspace: need 0
* (Rworkspace: need 3NN)	* RWorkspace: need N [e] + NN [RU] + NN [RVT] + 2NN [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 1059	Line 1155
*	*
* Multiply Q in U by real matrix RWORK(IRU), storing the	* Multiply Q in U by real matrix RWORK(IRU), storing the
* result in A, copying to U	* result in A, copying to U
* (CWorkspace: need 0)	* CWorkspace: need 0
* (Rworkspace: need NN+2M*N)	* RWorkspace: need N [e] + NN [RU] + 2MN [rwork] < N + 5NN since M < 2N 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 1068	Line 1164
CALL ZLACPY( 'F', M, N, A, LDA, U, LDU )	CALL ZLACPY( 'F', M, N, A, LDA, U, LDU )
ELSE	ELSE
*	*
	* Path 5a (M >> N, JOBZ='A')
* Copy A to VT, generate P**H	* Copy A to VT, generate P**H
* (Cworkspace: need 2N, prefer N+NNB)	* CWorkspace: need 2*N [tauq, taup] + N [work]
* (Rworkspace: 0)	* CWorkspace: prefer 2N [tauq, taup] + NNB [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 2N, prefer N+NNB)	* CWorkspace: need 2*N [tauq, taup] + M [work]
* (Rworkspace: 0)	* CWorkspace: prefer 2N [tauq, taup] + MNB [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 1087	Line 1186
* Perform bidiagonal SVD, computing left singular vectors	* Perform bidiagonal SVD, computing left singular vectors
* of bidiagonal matrix in RWORK(IRU) and computing right	* of bidiagonal matrix in RWORK(IRU) and computing right
* singular vectors of bidiagonal matrix in RWORK(IRVT)	* singular vectors of bidiagonal matrix in RWORK(IRVT)
* (CWorkspace: need 0)	* CWorkspace: need 0
* (RWorkspace: need BDSPAC)	* RWorkspace: need N [e] + NN [RU] + NN [RVT] + BDSPAC
*	*
IRU = NRWORK	IRU = NRWORK
IRVT = IRU + N*N	IRVT = IRU + N*N
Line 1099	Line 1198
*	*
* Multiply real matrix RWORK(IRVT) by P**H in VT,	* Multiply real matrix RWORK(IRVT) by P**H in VT,
* storing the result in A, copying to VT	* storing the result in A, copying to VT
* (Cworkspace: need 0)	* CWorkspace: need 0
* (Rworkspace: need 3NN)	* RWorkspace: need N [e] + NN [RU] + NN [RVT] + 2NN [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 1108	Line 1207
*	*
* Multiply Q in U by real matrix RWORK(IRU), storing the	* Multiply Q in U by real matrix RWORK(IRU), storing the
* result in A, copying to U	* result in A, copying to U
* (CWorkspace: 0)	* CWorkspace: need 0
* (Rworkspace: need 3NN)	* RWorkspace: need N [e] + NN [RU] + 2MN [rwork] < N + 5NN since M < 2N 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 1121	Line 1220
*	*
* M .LT. MNTHR2	* M .LT. MNTHR2
*	*
* Path 6 (M at least N, but not much larger)	* Path 6 (M >= N, but not much larger)
* Reduce to bidiagonal form without QR decomposition	* Reduce to bidiagonal form without QR decomposition
* Use ZUNMBR to compute singular vectors	* Use ZUNMBR to compute singular vectors
*	*
Line 1132	Line 1231
NWORK = ITAUP + N	NWORK = ITAUP + N
*	*
* Bidiagonalize A	* Bidiagonalize A
* (CWorkspace: need 2N+M, prefer 2N+(M+N)*NB)	* CWorkspace: need 2*N [tauq, taup] + M [work]
* (RWorkspace: need N)	* CWorkspace: prefer 2N [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.MN+3N ) THEN	IF( LWORK .GE. MN + 3N ) THEN
*	*
* WORK( IU ) is M by N	* WORK( IU ) is M by N
*	*
Line 1160	Line 1261
*	*
* WORK( IU ) is LDWRKU by N	* WORK( IU ) is LDWRKU by N
*	*
LDWRKU = ( LWORK-3*N ) / N	LDWRKU = ( LWORK - 3*N ) / N
END IF	END IF
NWORK = IU + LDWRKU*N	NWORK = IU + LDWRKU*N
*	*
	* Path 6o (M >= N, JOBZ='O')
* Perform bidiagonal SVD, computing left singular vectors	* Perform bidiagonal SVD, computing left singular vectors
* of bidiagonal matrix in RWORK(IRU) and computing right	* of bidiagonal matrix in RWORK(IRU) and computing right
* singular vectors of bidiagonal matrix in RWORK(IRVT)	* singular vectors of bidiagonal matrix in RWORK(IRVT)
* (CWorkspace: need 0)	* CWorkspace: need 0
* (RWorkspace: need BDSPAC)	* RWorkspace: need N [e] + NN [RU] + NN [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 1176	Line 1278
*	*
* Copy real matrix RWORK(IRVT) to complex matrix VT	* Copy real matrix RWORK(IRVT) to complex matrix VT
* Overwrite VT by right singular vectors of A	* Overwrite VT by right singular vectors of A
* (Cworkspace: need 2N, prefer N+NNB)	* CWorkspace: need 2N [tauq, taup] + NN [U] + N [work]
* (Rworkspace: need 0)	* CWorkspace: prefer 2N [tauq, taup] + NN [U] + N*NB [work]
	* RWorkspace: need N [e] + NN [RU] + NN [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.MN+3N ) THEN	IF( LWORK .GE. MN + 3N ) 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 MN+2N, prefer MN+N+NNB)	* to A
* (Rworkspace: need 0)	* CWorkspace: need 2N [tauq, taup] + MN [U] + N [work]
	* CWorkspace: prefer 2N [tauq, taup] + MN [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 1202	Line 1307
CALL ZLACPY( 'F', M, N, WORK( IU ), LDWRKU, A, LDA )	CALL ZLACPY( 'F', M, N, WORK( IU ), LDWRKU, A, LDA )
ELSE	ELSE
*	*
	* Path 6o-slow
* Generate Q in A	* Generate Q in A
* (Cworkspace: need 2N, prefer N+NNB)	* CWorkspace: need 2N [tauq, taup] + NN [U] + N [work]
* (Rworkspace: need 0)	* CWorkspace: prefer 2N [tauq, taup] + NN [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 NN, prefer MN)	* CWorkspace: need 2N [tauq, taup] + NN [U]
* (Rworkspace: need 3NN, prefer NN+2M*N)	* CWorkspace: prefer 2N [tauq, taup] + MN [U]
	* RWorkspace: need N [e] + NN [RU] + 2N*N [rwork]
	* RWorkspace: prefer N [e] + NN [RU] + 2MN [rwork] < N + 5NN since M < 2N here
*	*
NRWORK = IRVT	NRWORK = IRVT
DO 30 I = 1, M, LDWRKU	DO 30 I = 1, M, LDWRKU
Line 1227	Line 1336
*	*
ELSE IF( WNTQS ) THEN	ELSE IF( WNTQS ) THEN
*	*
	* Path 6s (M >= N, JOBZ='S')
* Perform bidiagonal SVD, computing left singular vectors	* Perform bidiagonal SVD, computing left singular vectors
* of bidiagonal matrix in RWORK(IRU) and computing right	* of bidiagonal matrix in RWORK(IRU) and computing right
* singular vectors of bidiagonal matrix in RWORK(IRVT)	* singular vectors of bidiagonal matrix in RWORK(IRVT)
* (CWorkspace: need 0)	* CWorkspace: need 0
* (RWorkspace: need BDSPAC)	* RWorkspace: need N [e] + NN [RU] + NN [RVT] + BDSPAC
*	*
IRU = NRWORK	IRU = NRWORK
IRVT = IRU + N*N	IRVT = IRU + N*N
Line 1242	Line 1352
*	*
* Copy real matrix RWORK(IRU) to complex matrix U	* Copy real matrix RWORK(IRU) to complex matrix U
* Overwrite U by left singular vectors of A	* Overwrite U by left singular vectors of A
* (CWorkspace: need 3N, prefer 2N+N*NB)	* CWorkspace: need 2*N [tauq, taup] + N [work]
* (RWorkspace: 0)	* CWorkspace: prefer 2N [tauq, taup] + NNB [work]
	* RWorkspace: need N [e] + NN [RU] + NN [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 1253	Line 1364
*	*
* Copy real matrix RWORK(IRVT) to complex matrix VT	* Copy real matrix RWORK(IRVT) to complex matrix VT
* Overwrite VT by right singular vectors of A	* Overwrite VT by right singular vectors of A
* (CWorkspace: need 3N, prefer 2N+N*NB)	* CWorkspace: need 2*N [tauq, taup] + N [work]
* (RWorkspace: 0)	* CWorkspace: prefer 2N [tauq, taup] + NNB [work]
	* RWorkspace: need N [e] + NN [RU] + NN [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 1262	Line 1374
$ LWORK-NWORK+1, IERR )	$ LWORK-NWORK+1, IERR )
ELSE	ELSE
*	*
	* Path 6a (M >= N, JOBZ='A')
* Perform bidiagonal SVD, computing left singular vectors	* Perform bidiagonal SVD, computing left singular vectors
* of bidiagonal matrix in RWORK(IRU) and computing right	* of bidiagonal matrix in RWORK(IRU) and computing right
* singular vectors of bidiagonal matrix in RWORK(IRVT)	* singular vectors of bidiagonal matrix in RWORK(IRVT)
* (CWorkspace: need 0)	* CWorkspace: need 0
* (RWorkspace: need BDSPAC)	* RWorkspace: need N [e] + NN [RU] + NN [RVT] + BDSPAC
*	*
IRU = NRWORK	IRU = NRWORK
IRVT = IRU + N*N	IRVT = IRU + N*N
Line 1285	Line 1398
*	*
* Copy real matrix RWORK(IRU) to complex matrix U	* Copy real matrix RWORK(IRU) to complex matrix U
* Overwrite U by left singular vectors of A	* Overwrite U by left singular vectors of A
* (CWorkspace: need 2N+M, prefer 2N+M*NB)	* CWorkspace: need 2*N [tauq, taup] + M [work]
* (RWorkspace: 0)	* CWorkspace: prefer 2N [tauq, taup] + MNB [work]
	* RWorkspace: need N [e] + NN [RU] + NN [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 1295	Line 1409
*	*
* Copy real matrix RWORK(IRVT) to complex matrix VT	* Copy real matrix RWORK(IRVT) to complex matrix VT
* Overwrite VT by right singular vectors of A	* Overwrite VT by right singular vectors of A
* (CWorkspace: need 3N, prefer 2N+N*NB)	* CWorkspace: need 2*N [tauq, taup] + N [work]
* (RWorkspace: 0)	* CWorkspace: prefer 2N [tauq, taup] + NNB [work]
	* RWorkspace: need N [e] + NN [RU] + NN [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 1316	Line 1431
*	*
IF( WNTQN ) THEN	IF( WNTQN ) THEN
*	*
* Path 1t (N much larger than M, JOBZ='N')	* Path 1t (N >> M, JOBZ='N')
* No singular vectors to be computed	* No singular vectors to be computed
*	*
ITAU = 1	ITAU = 1
NWORK = ITAU + M	NWORK = ITAU + M
*	*
* Compute A=L*Q	* Compute A=L*Q
* (CWorkspace: need 2M, prefer M+MNB)	* 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 1339	Line 1455
NWORK = ITAUP + M	NWORK = ITAUP + M
*	*
* Bidiagonalize L in A	* Bidiagonalize L in A
* (CWorkspace: need 3M, prefer 2M+2MNB)	* CWorkspace: need 2*M [tauq, taup] + M [work]
* (RWorkspace: need M)	* CWorkspace: prefer 2M [tauq, taup] + 2M*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 1348	Line 1465
NRWORK = IE + M	NRWORK = IE + M
*	*
* Perform bidiagonal SVD, compute singular values only	* Perform bidiagonal SVD, compute singular values only
* (CWorkspace: 0)	* CWorkspace: need 0
* (RWorkspace: need BDSPAN)	* RWorkspace: need M [e] + BDSPAC
*	*
CALL DBDSDC( 'U', 'N', M, S, RWORK( IE ), DUM, 1, DUM, 1,	CALL DBDSDC( 'U', 'N', M, S, RWORK( IE ), DUM,1,DUM,1,
$ DUM, IDUM, RWORK( NRWORK ), IWORK, INFO )	$ DUM, IDUM, RWORK( NRWORK ), IWORK, INFO )
*	*
ELSE IF( WNTQO ) THEN	ELSE IF( WNTQO ) THEN
*	*
* Path 2t (N much larger than M, JOBZ='O')	* Path 2t (N >> M, JOBZ='O')
* M right singular vectors to be overwritten on A and	* M right singular vectors to be overwritten on A and
* M left singular vectors to be computed in U	* M left singular vectors to be computed in U
*	*
Line 1366	Line 1483
* WORK(IVT) is M by M	* WORK(IVT) is M by M
*	*
IL = IVT + LDWKVT*M	IL = IVT + LDWKVT*M
IF( LWORK.GE.MN+MM+3*M ) THEN	IF( LWORK .GE. MN + MM + 3*M ) THEN
*	*
* WORK(IL) M by N	* WORK(IL) M by N
*	*
Line 1377	Line 1494
* WORK(IL) is M by CHUNK	* WORK(IL) is M by CHUNK
*	*
LDWRKL = M	LDWRKL = M
CHUNK = ( LWORK-MM-3M ) / M	CHUNK = ( LWORK - MM - 3M ) / 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 2M, prefer M+MNB)	* CWorkspace: need MM [VT] + MM [L] + M [tau] + M [work]
* (RWorkspace: 0)	* CWorkspace: prefer MM [VT] + MM [L] + M [tau] + M*NB [work]
	* RWorkspace: need 0
*	*
CALL ZGELQF( M, N, A, LDA, WORK( ITAU ), WORK( NWORK ),	CALL ZGELQF( M, N, A, LDA, WORK( ITAU ), WORK( NWORK ),
$ LWORK-NWORK+1, IERR )	$ LWORK-NWORK+1, IERR )
Line 1396	Line 1514
$ WORK( IL+LDWRKL ), LDWRKL )	$ WORK( IL+LDWRKL ), LDWRKL )
*	*
* Generate Q in A	* Generate Q in A
* (CWorkspace: need MM+2M, prefer MM+M+MNB)	* CWorkspace: need MM [VT] + MM [L] + M [tau] + M [work]
* (RWorkspace: 0)	* CWorkspace: prefer MM [VT] + MM [L] + M [tau] + M*NB [work]
	* RWorkspace: need 0
*	*
CALL ZUNGLQ( M, N, M, A, LDA, WORK( ITAU ),	CALL ZUNGLQ( M, N, M, A, LDA, WORK( ITAU ),
$ WORK( NWORK ), LWORK-NWORK+1, IERR )	$ WORK( NWORK ), LWORK-NWORK+1, IERR )
Line 1407	Line 1526
NWORK = ITAUP + M	NWORK = ITAUP + M
*	*
* Bidiagonalize L in WORK(IL)	* Bidiagonalize L in WORK(IL)
* (CWorkspace: need MM+3M, prefer MM+2M+2MNB)	* CWorkspace: need MM [VT] + MM [L] + 2*M [tauq, taup] + M [work]
* (RWorkspace: need M)	* CWorkspace: prefer MM [VT] + MM [L] + 2M [tauq, taup] + 2M*NB [work]
	* RWorkspace: need M [e]
*	*
CALL ZGEBRD( M, M, WORK( IL ), LDWRKL, S, RWORK( IE ),	CALL ZGEBRD( M, M, WORK( IL ), LDWRKL, S, RWORK( IE ),
$ WORK( ITAUQ ), WORK( ITAUP ), WORK( NWORK ),	$ WORK( ITAUQ ), WORK( ITAUP ), WORK( NWORK ),
Line 1417	Line 1537
* Perform bidiagonal SVD, computing left singular vectors	* Perform bidiagonal SVD, computing left singular vectors
* of bidiagonal matrix in RWORK(IRU) and computing right	* of bidiagonal matrix in RWORK(IRU) and computing right
* singular vectors of bidiagonal matrix in RWORK(IRVT)	* singular vectors of bidiagonal matrix in RWORK(IRVT)
* (CWorkspace: need 0)	* CWorkspace: need 0
* (RWorkspace: need BDSPAC)	* RWorkspace: need M [e] + MM [RU] + MM [RVT] + BDSPAC
*	*
IRU = IE + M	IRU = IE + M
IRVT = IRU + M*M	IRVT = IRU + M*M
Line 1429	Line 1549
*	*
* Copy real matrix RWORK(IRU) to complex matrix WORK(IU)	* Copy real matrix RWORK(IRU) to complex matrix WORK(IU)
* Overwrite WORK(IU) by the left singular vectors of L	* Overwrite WORK(IU) by the left singular vectors of L
* (CWorkspace: need NN+3N, prefer MN+2N+N*NB)	* CWorkspace: need MM [VT] + MM [L] + 2*M [tauq, taup] + M [work]
* (RWorkspace: 0)	* CWorkspace: prefer MM [VT] + MM [L] + 2M [tauq, taup] + MNB [work]
	* RWorkspace: need 0
*	*
CALL ZLACP2( 'F', M, M, RWORK( IRU ), M, U, LDU )	CALL ZLACP2( 'F', M, M, RWORK( IRU ), M, U, LDU )
CALL ZUNMBR( 'Q', 'L', 'N', M, M, M, WORK( IL ), LDWRKL,	CALL ZUNMBR( 'Q', 'L', 'N', M, M, M, WORK( IL ), LDWRKL,
Line 1439	Line 1560
*	*
* Copy real matrix RWORK(IRVT) to complex matrix WORK(IVT)	* Copy real matrix RWORK(IRVT) to complex matrix WORK(IVT)
* Overwrite WORK(IVT) by the right singular vectors of L	* Overwrite WORK(IVT) by the right singular vectors of L
* (CWorkspace: need NN+3N, prefer MN+2N+N*NB)	* CWorkspace: need MM [VT] + MM [L] + 2*M [tauq, taup] + M [work]
* (RWorkspace: 0)	* CWorkspace: prefer MM [VT] + MM [L] + 2M [tauq, taup] + MNB [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 1450	Line 1572
*	*
* Multiply right singular vectors of L in WORK(IL) by Q	* Multiply right singular vectors of L in WORK(IL) by Q
* in A, storing result in WORK(IL) and copying to A	* in A, storing result in WORK(IL) and copying to A
* (CWorkspace: need 2MM, prefer MM+MN))	* CWorkspace: need MM [VT] + MM [L]
* (RWorkspace: 0)	* CWorkspace: prefer MM [VT] + MN [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 1464	Line 1587
*	*
ELSE IF( WNTQS ) THEN	ELSE IF( WNTQS ) THEN
*	*
* Path 3t (N much larger than M, JOBZ='S')	* Path 3t (N >> M, JOBZ='S')
* M right singular vectors to be computed in VT and	* M right singular vectors to be computed in VT and
* M left singular vectors to be computed in U	* M left singular vectors to be computed in U
*	*
IL = 1	IL = 1
*	*
Line 1477	Line 1600
NWORK = ITAU + M	NWORK = ITAU + M
*	*
* Compute A=L*Q	* Compute A=L*Q
* (CWorkspace: need 2M, prefer M+MNB)	* CWorkspace: need M*M [L] + M [tau] + M [work]
* (RWorkspace: 0)	* CWorkspace: prefer MM [L] + M [tau] + MNB [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 1490	Line 1614
$ WORK( IL+LDWRKL ), LDWRKL )	$ WORK( IL+LDWRKL ), LDWRKL )
*	*
* Generate Q in A	* Generate Q in A
* (CWorkspace: need MM+2M, prefer MM+M+MNB)	* CWorkspace: need M*M [L] + M [tau] + M [work]
* (RWorkspace: 0)	* CWorkspace: prefer MM [L] + M [tau] + MNB [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 1501	Line 1626
NWORK = ITAUP + M	NWORK = ITAUP + M
*	*
* Bidiagonalize L in WORK(IL)	* Bidiagonalize L in WORK(IL)
* (CWorkspace: need MM+3M, prefer MM+2M+2MNB)	* CWorkspace: need MM [L] + 2M [tauq, taup] + M [work]
* (RWorkspace: need M)	* CWorkspace: prefer MM [L] + 2M [tauq, taup] + 2MNB [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 1511	Line 1637
* Perform bidiagonal SVD, computing left singular vectors	* Perform bidiagonal SVD, computing left singular vectors
* of bidiagonal matrix in RWORK(IRU) and computing right	* of bidiagonal matrix in RWORK(IRU) and computing right
* singular vectors of bidiagonal matrix in RWORK(IRVT)	* singular vectors of bidiagonal matrix in RWORK(IRVT)
* (CWorkspace: need 0)	* CWorkspace: need 0
* (RWorkspace: need BDSPAC)	* RWorkspace: need M [e] + MM [RU] + MM [RVT] + BDSPAC
*	*
IRU = IE + M	IRU = IE + M
IRVT = IRU + M*M	IRVT = IRU + M*M
Line 1523	Line 1649
*	*
* Copy real matrix RWORK(IRU) to complex matrix U	* Copy real matrix RWORK(IRU) to complex matrix U
* Overwrite U by left singular vectors of L	* Overwrite U by left singular vectors of L
* (CWorkspace: need MM+3M, prefer MM+2M+M*NB)	* CWorkspace: need MM [L] + 2M [tauq, taup] + M [work]
* (RWorkspace: 0)	* CWorkspace: prefer MM [L] + 2M [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 1533	Line 1660
*	*
* Copy real matrix RWORK(IRVT) to complex matrix VT	* Copy real matrix RWORK(IRVT) to complex matrix VT
* Overwrite VT by left singular vectors of L	* Overwrite VT by left singular vectors of L
* (CWorkspace: need MM+3M, prefer MM+2M+M*NB)	* CWorkspace: need MM [L] + 2M [tauq, taup] + M [work]
* (RWorkspace: 0)	* CWorkspace: prefer MM [L] + 2M [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 1543	Line 1671
*	*
* Copy VT to WORK(IL), multiply right singular vectors of L	* Copy VT to WORK(IL), multiply right singular vectors of L
* in WORK(IL) by Q in A, storing result in VT	* in WORK(IL) by Q in A, storing result in VT
* (CWorkspace: need M*M)	* CWorkspace: need M*M [L]
* (RWorkspace: 0)	* RWorkspace: need 0
*	*
CALL ZLACPY( 'F', M, M, VT, LDVT, WORK( IL ), LDWRKL )	CALL ZLACPY( 'F', M, M, VT, LDVT, WORK( IL ), LDWRKL )
CALL ZGEMM( 'N', 'N', M, N, M, CONE, WORK( IL ), LDWRKL,	CALL ZGEMM( 'N', 'N', M, N, M, CONE, WORK( IL ), LDWRKL,
Line 1552	Line 1680
*	*
ELSE IF( WNTQA ) THEN	ELSE IF( WNTQA ) THEN
*	*
* Path 9t (N much larger than M, JOBZ='A')	* Path 4t (N >> M, JOBZ='A')
* N right singular vectors to be computed in VT and	* N right singular vectors to be computed in VT and
* M left singular vectors to be computed in U	* M left singular vectors to be computed in U
*	*
Line 1565	Line 1693
NWORK = ITAU + M	NWORK = ITAU + M
*	*
* Compute A=L*Q, copying result to VT	* Compute A=L*Q, copying result to VT
* (CWorkspace: need 2M, prefer M+MNB)	* CWorkspace: need M*M [VT] + M [tau] + M [work]
* (RWorkspace: 0)	* CWorkspace: prefer MM [VT] + M [tau] + MNB [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 MM [VT] + M [tau] + NNB [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 1589	Line 1719
NWORK = ITAUP + M	NWORK = ITAUP + M
*	*
* Bidiagonalize L in A	* Bidiagonalize L in A
* (CWorkspace: need MM+3M, prefer MM+2M+2MNB)	* CWorkspace: need MM [VT] + 2M [tauq, taup] + M [work]
* (RWorkspace: need M)	* CWorkspace: prefer MM [VT] + 2M [tauq, taup] + 2MNB [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 1599	Line 1730
* Perform bidiagonal SVD, computing left singular vectors	* Perform bidiagonal SVD, computing left singular vectors
* of bidiagonal matrix in RWORK(IRU) and computing right	* of bidiagonal matrix in RWORK(IRU) and computing right
* singular vectors of bidiagonal matrix in RWORK(IRVT)	* singular vectors of bidiagonal matrix in RWORK(IRVT)
* (CWorkspace: need 0)	* CWorkspace: need 0
* (RWorkspace: need BDSPAC)	* RWorkspace: need M [e] + MM [RU] + MM [RVT] + BDSPAC
*	*
IRU = IE + M	IRU = IE + M
IRVT = IRU + M*M	IRVT = IRU + M*M
Line 1611	Line 1742
*	*
* Copy real matrix RWORK(IRU) to complex matrix U	* Copy real matrix RWORK(IRU) to complex matrix U
* Overwrite U by left singular vectors of L	* Overwrite U by left singular vectors of L
* (CWorkspace: need 3M, prefer 2M+M*NB)	* CWorkspace: need MM [VT] + 2M [tauq, taup] + M [work]
* (RWorkspace: 0)	* CWorkspace: prefer MM [VT] + 2M [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 1621	Line 1753
*	*
* Copy real matrix RWORK(IRVT) to complex matrix WORK(IVT)	* Copy real matrix RWORK(IRVT) to complex matrix WORK(IVT)
* Overwrite WORK(IVT) by right singular vectors of L	* Overwrite WORK(IVT) by right singular vectors of L
* (CWorkspace: need MM+3M, prefer MM+2M+M*NB)	* CWorkspace: need MM [VT] + 2M [tauq, taup] + M [work]
* (RWorkspace: 0)	* CWorkspace: prefer MM [VT] + 2M [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 1632	Line 1765
*	*
* Multiply right singular vectors of L in WORK(IVT) by	* Multiply right singular vectors of L in WORK(IVT) by
* Q in VT, storing result in A	* Q in VT, storing result in A
* (CWorkspace: need M*M)	* CWorkspace: need M*M [VT]
* (RWorkspace: 0)	* RWorkspace: need 0
*	*
CALL ZGEMM( 'N', 'N', M, N, M, CONE, WORK( IVT ), LDWKVT,	CALL ZGEMM( 'N', 'N', M, N, M, CONE, WORK( IVT ), LDWKVT,
$ VT, LDVT, CZERO, A, LDA )	$ VT, LDVT, CZERO, A, LDA )
Line 1648	Line 1781
*	*
* MNTHR2 <= N < MNTHR1	* MNTHR2 <= N < MNTHR1
*	*
* Path 5t (N much larger than M, but not as much as MNTHR1)	* Path 5t (N >> M, but not as much as MNTHR1)
* Reduce to bidiagonal form without QR decomposition, use	* Reduce to bidiagonal form without QR decomposition, use
* ZUNGBR and matrix multiplication to compute singular vectors	* ZUNGBR and matrix multiplication to compute singular vectors
*	*
*
IE = 1	IE = 1
NRWORK = IE + M	NRWORK = IE + M
ITAUQ = 1	ITAUQ = 1
Line 1660	Line 1792
NWORK = ITAUP + M	NWORK = ITAUP + M
*	*
* Bidiagonalize A	* Bidiagonalize A
* (CWorkspace: need 2M+N, prefer 2M+(M+N)*NB)	* CWorkspace: need 2*M [tauq, taup] + N [work]
* (RWorkspace: M)	* CWorkspace: prefer 2M [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 1669	Line 1802
*	*
IF( WNTQN ) THEN	IF( WNTQN ) THEN
*	*
	* Path 5tn (N >> M, JOBZ='N')
* Compute singular values only	* Compute singular values only
* (Cworkspace: 0)	* CWorkspace: need 0
* (Rworkspace: need BDSPAN)	* RWorkspace: need M [e] + BDSPAC
*	*
CALL DBDSDC( 'L', 'N', M, S, RWORK( IE ), DUM, 1, DUM, 1,	CALL DBDSDC( 'L', 'N', M, S, RWORK( IE ), DUM,1,DUM,1,
$ DUM, IDUM, RWORK( NRWORK ), IWORK, INFO )	$ DUM, IDUM, RWORK( NRWORK ), IWORK, INFO )
ELSE IF( WNTQO ) THEN	ELSE IF( WNTQO ) THEN
IRVT = NRWORK	IRVT = NRWORK
Line 1681	Line 1815
NRWORK = IRU + M*M	NRWORK = IRU + M*M
IVT = NWORK	IVT = NWORK
*	*
	* Path 5to (N >> M, JOBZ='O')
* Copy A to U, generate Q	* Copy A to U, generate Q
* (Cworkspace: need 2M, prefer M+MNB)	* CWorkspace: need 2*M [tauq, taup] + M [work]
* (Rworkspace: 0)	* CWorkspace: prefer 2M [tauq, taup] + MNB [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 2M, prefer M+MNB)	* CWorkspace: need 2*M [tauq, taup] + M [work]
* (Rworkspace: 0)	* CWorkspace: prefer 2M [tauq, taup] + MNB [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.MN+3M ) THEN	IF( LWORK .GE. MN + 3M ) THEN
*	*
* WORK( IVT ) is M by N	* WORK( IVT ) is M by N
*	*
Line 1707	Line 1844
*	*
* WORK( IVT ) is M by CHUNK	* WORK( IVT ) is M by CHUNK
*	*
CHUNK = ( LWORK-3*M ) / M	CHUNK = ( LWORK - 3*M ) / M
NWORK = IVT + LDWKVT*CHUNK	NWORK = IVT + LDWKVT*CHUNK
END IF	END IF
*	*
* Perform bidiagonal SVD, computing left singular vectors	* Perform bidiagonal SVD, computing left singular vectors
* of bidiagonal matrix in RWORK(IRU) and computing right	* of bidiagonal matrix in RWORK(IRU) and computing right
* singular vectors of bidiagonal matrix in RWORK(IRVT)	* singular vectors of bidiagonal matrix in RWORK(IRVT)
* (CWorkspace: need 0)	* CWorkspace: need 0
* (RWorkspace: need BDSPAC)	* RWorkspace: need M [e] + MM [RVT] + MM [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 1723	Line 1860
*	*
* Multiply Q in U by real matrix RWORK(IRVT)	* Multiply Q in U by real matrix RWORK(IRVT)
* storing the result in WORK(IVT), copying to U	* storing the result in WORK(IVT), copying to U
* (Cworkspace: need 0)	* CWorkspace: need 2M [tauq, taup] + MM [VT]
* (Rworkspace: need 2MM)	* RWorkspace: need M [e] + MM [RVT] + MM [RU] + 2MM [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 1732	Line 1869
*	*
* Multiply RWORK(IRVT) by P**H in A, storing the	* Multiply RWORK(IRVT) by P**H in A, storing the
* result in WORK(IVT), copying to A	* result in WORK(IVT), copying to A
* (CWorkspace: need MM, prefer MN)	* CWorkspace: need 2M [tauq, taup] + MM [VT]
* (Rworkspace: need 2MM, prefer 2MN)	* CWorkspace: prefer 2M [tauq, taup] + MN [VT]
	* RWorkspace: need M [e] + MM [RVT] + 2M*M [rwork]
	* RWorkspace: prefer M [e] + MM [RVT] + 2MN [rwork] < M + 5MM since N < 2M here
*	*
NRWORK = IRU	NRWORK = IRU
DO 50 I = 1, N, CHUNK	DO 50 I = 1, N, CHUNK
Line 1745	Line 1884
50 CONTINUE	50 CONTINUE
ELSE IF( WNTQS ) THEN	ELSE IF( WNTQS ) THEN
*	*
	* Path 5ts (N >> M, JOBZ='S')
* Copy A to U, generate Q	* Copy A to U, generate Q
* (Cworkspace: need 2M, prefer M+MNB)	* CWorkspace: need 2*M [tauq, taup] + M [work]
* (Rworkspace: 0)	* CWorkspace: prefer 2M [tauq, taup] + MNB [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 2M, prefer M+MNB)	* CWorkspace: need 2*M [tauq, taup] + M [work]
* (Rworkspace: 0)	* CWorkspace: prefer 2M [tauq, taup] + MNB [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 1764	Line 1906
* Perform bidiagonal SVD, computing left singular vectors	* Perform bidiagonal SVD, computing left singular vectors
* of bidiagonal matrix in RWORK(IRU) and computing right	* of bidiagonal matrix in RWORK(IRU) and computing right
* singular vectors of bidiagonal matrix in RWORK(IRVT)	* singular vectors of bidiagonal matrix in RWORK(IRVT)
* (CWorkspace: need 0)	* CWorkspace: need 0
* (RWorkspace: need BDSPAC)	* RWorkspace: need M [e] + MM [RVT] + MM [RU] + BDSPAC
*	*
IRVT = NRWORK	IRVT = NRWORK
IRU = IRVT + M*M	IRU = IRVT + M*M
Line 1776	Line 1918
*	*
* Multiply Q in U by real matrix RWORK(IRU), storing the	* Multiply Q in U by real matrix RWORK(IRU), storing the
* result in A, copying to U	* result in A, copying to U
* (CWorkspace: need 0)	* CWorkspace: need 0
* (Rworkspace: need 3MM)	* RWorkspace: need M [e] + MM [RVT] + MM [RU] + 2MM [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 1785	Line 1927
*	*
* Multiply real matrix RWORK(IRVT) by P**H in VT,	* Multiply real matrix RWORK(IRVT) by P**H in VT,
* storing the result in A, copying to VT	* storing the result in A, copying to VT
* (Cworkspace: need 0)	* CWorkspace: need 0
* (Rworkspace: need MM+2M*N)	* RWorkspace: need M [e] + MM [RVT] + 2MN [rwork] < M + 5MM since N < 2M 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 1794	Line 1936
CALL ZLACPY( 'F', M, N, A, LDA, VT, LDVT )	CALL ZLACPY( 'F', M, N, A, LDA, VT, LDVT )
ELSE	ELSE
*	*
	* Path 5ta (N >> M, JOBZ='A')
* Copy A to U, generate Q	* Copy A to U, generate Q
* (Cworkspace: need 2M, prefer M+MNB)	* CWorkspace: need 2*M [tauq, taup] + M [work]
* (Rworkspace: 0)	* CWorkspace: prefer 2M [tauq, taup] + MNB [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 2M, prefer M+MNB)	* CWorkspace: need 2*M [tauq, taup] + N [work]
* (Rworkspace: 0)	* CWorkspace: prefer 2M [tauq, taup] + NNB [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 1813	Line 1958
* Perform bidiagonal SVD, computing left singular vectors	* Perform bidiagonal SVD, computing left singular vectors
* of bidiagonal matrix in RWORK(IRU) and computing right	* of bidiagonal matrix in RWORK(IRU) and computing right
* singular vectors of bidiagonal matrix in RWORK(IRVT)	* singular vectors of bidiagonal matrix in RWORK(IRVT)
* (CWorkspace: need 0)	* CWorkspace: need 0
* (RWorkspace: need BDSPAC)	* RWorkspace: need M [e] + MM [RVT] + MM [RU] + BDSPAC
*	*
IRVT = NRWORK	IRVT = NRWORK
IRU = IRVT + M*M	IRU = IRVT + M*M
Line 1825	Line 1970
*	*
* Multiply Q in U by real matrix RWORK(IRU), storing the	* Multiply Q in U by real matrix RWORK(IRU), storing the
* result in A, copying to U	* result in A, copying to U
* (CWorkspace: need 0)	* CWorkspace: need 0
* (Rworkspace: need 3MM)	* RWorkspace: need M [e] + MM [RVT] + MM [RU] + 2MM [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 1834	Line 1979
*	*
* Multiply real matrix RWORK(IRVT) by P**H in VT,	* Multiply real matrix RWORK(IRVT) by P**H in VT,
* storing the result in A, copying to VT	* storing the result in A, copying to VT
* (Cworkspace: need 0)	* CWorkspace: need 0
* (Rworkspace: need MM+2M*N)	* RWorkspace: need M [e] + MM [RVT] + 2MN [rwork] < M + 5MM since N < 2M 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 1846	Line 1992
*	*
* N .LT. MNTHR2	* N .LT. MNTHR2
*	*
* Path 6t (N greater than M, but not much larger)	* Path 6t (N > M, but not much larger)
* Reduce to bidiagonal form without LQ decomposition	* Reduce to bidiagonal form without LQ decomposition
* Use ZUNMBR to compute singular vectors	* Use ZUNMBR to compute singular vectors
*	*
Line 1857	Line 2003
NWORK = ITAUP + M	NWORK = ITAUP + M
*	*
* Bidiagonalize A	* Bidiagonalize A
* (CWorkspace: need 2M+N, prefer 2M+(M+N)*NB)	* CWorkspace: need 2*M [tauq, taup] + N [work]
* (RWorkspace: M)	* CWorkspace: prefer 2M [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.MN+3M ) THEN	IF( LWORK .GE. MN + 3M ) THEN
*	*
* WORK( IVT ) is M by N	* WORK( IVT ) is M by N
*	*
Line 1885	Line 2034
*	*
* WORK( IVT ) is M by CHUNK	* WORK( IVT ) is M by CHUNK
*	*
CHUNK = ( LWORK-3*M ) / M	CHUNK = ( LWORK - 3*M ) / M
NWORK = IVT + LDWKVT*CHUNK	NWORK = IVT + LDWKVT*CHUNK
END IF	END IF
*	*
* Perform bidiagonal SVD, computing left singular vectors	* Perform bidiagonal SVD, computing left singular vectors
* of bidiagonal matrix in RWORK(IRU) and computing right	* of bidiagonal matrix in RWORK(IRU) and computing right
* singular vectors of bidiagonal matrix in RWORK(IRVT)	* singular vectors of bidiagonal matrix in RWORK(IRVT)
* (CWorkspace: need 0)	* CWorkspace: need 0
* (RWorkspace: need BDSPAC)	* RWorkspace: need M [e] + MM [RVT] + MM [RU] + BDSPAC
*	*
IRVT = NRWORK	IRVT = NRWORK
IRU = IRVT + M*M	IRU = IRVT + M*M
Line 1904	Line 2053
*	*
* Copy real matrix RWORK(IRU) to complex matrix U	* Copy real matrix RWORK(IRU) to complex matrix U
* Overwrite U by left singular vectors of A	* Overwrite U by left singular vectors of A
* (Cworkspace: need 2M, prefer M+MNB)	* CWorkspace: need 2M [tauq, taup] + MM [VT] + M [work]
* (Rworkspace: need 0)	* CWorkspace: prefer 2M [tauq, taup] + MM [VT] + M*NB [work]
	* RWorkspace: need M [e] + MM [RVT] + MM [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.MN+3M ) THEN	IF( LWORK .GE. MN + 3M ) 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 MN+2M, prefer MN+M+MNB)	* copying to A
* (Rworkspace: need 0)	* CWorkspace: need 2M [tauq, taup] + MN [VT] + M [work]
	* CWorkspace: prefer 2M [tauq, taup] + MN [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 1928	Line 2080
CALL ZLACPY( 'F', M, N, WORK( IVT ), LDWKVT, A, LDA )	CALL ZLACPY( 'F', M, N, WORK( IVT ), LDWKVT, A, LDA )
ELSE	ELSE
*	*
	* Path 6to-slow
* Generate P**H in A	* Generate P**H in A
* (Cworkspace: need 2M, prefer M+MNB)	* CWorkspace: need 2M [tauq, taup] + MM [VT] + M [work]
* (Rworkspace: need 0)	* CWorkspace: prefer 2M [tauq, taup] + MM [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 MM, prefer MN)	* CWorkspace: need 2M [tauq, taup] + MM [VT]
* (Rworkspace: need 3MM, prefer MM+2M*N)	* CWorkspace: prefer 2M [tauq, taup] + MN [VT]
	* RWorkspace: need M [e] + MM [RVT] + 2M*M [rwork]
	* RWorkspace: prefer M [e] + MM [RVT] + 2MN [rwork] < M + 5MM since N < 2M here
*	*
NRWORK = IRU	NRWORK = IRU
DO 60 I = 1, N, CHUNK	DO 60 I = 1, N, CHUNK
Line 1952	Line 2108
END IF	END IF
ELSE IF( WNTQS ) THEN	ELSE IF( WNTQS ) THEN
*	*
	* Path 6ts (N > M, JOBZ='S')
* Perform bidiagonal SVD, computing left singular vectors	* Perform bidiagonal SVD, computing left singular vectors
* of bidiagonal matrix in RWORK(IRU) and computing right	* of bidiagonal matrix in RWORK(IRU) and computing right
* singular vectors of bidiagonal matrix in RWORK(IRVT)	* singular vectors of bidiagonal matrix in RWORK(IRVT)
* (CWorkspace: need 0)	* CWorkspace: need 0
* (RWorkspace: need BDSPAC)	* RWorkspace: need M [e] + MM [RVT] + MM [RU] + BDSPAC
*	*
IRVT = NRWORK	IRVT = NRWORK
IRU = IRVT + M*M	IRU = IRVT + M*M
Line 1967	Line 2124
*	*
* Copy real matrix RWORK(IRU) to complex matrix U	* Copy real matrix RWORK(IRU) to complex matrix U
* Overwrite U by left singular vectors of A	* Overwrite U by left singular vectors of A
* (CWorkspace: need 3M, prefer 2M+M*NB)	* CWorkspace: need 2*M [tauq, taup] + M [work]
* (RWorkspace: M*M)	* CWorkspace: prefer 2M [tauq, taup] + MNB [work]
	* RWorkspace: need M [e] + MM [RVT] + MM [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 1977	Line 2135
*	*
* Copy real matrix RWORK(IRVT) to complex matrix VT	* Copy real matrix RWORK(IRVT) to complex matrix VT
* Overwrite VT by right singular vectors of A	* Overwrite VT by right singular vectors of A
* (CWorkspace: need 3M, prefer 2M+M*NB)	* CWorkspace: need 2*M [tauq, taup] + M [work]
* (RWorkspace: M*M)	* CWorkspace: prefer 2M [tauq, taup] + MNB [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 1987	Line 2146
$ LWORK-NWORK+1, IERR )	$ LWORK-NWORK+1, IERR )
ELSE	ELSE
*	*
	* Path 6ta (N > M, JOBZ='A')
* Perform bidiagonal SVD, computing left singular vectors	* Perform bidiagonal SVD, computing left singular vectors
* of bidiagonal matrix in RWORK(IRU) and computing right	* of bidiagonal matrix in RWORK(IRU) and computing right
* singular vectors of bidiagonal matrix in RWORK(IRVT)	* singular vectors of bidiagonal matrix in RWORK(IRVT)
* (CWorkspace: need 0)	* CWorkspace: need 0
* (RWorkspace: need BDSPAC)	* RWorkspace: need M [e] + MM [RVT] + MM [RU] + BDSPAC
*	*
IRVT = NRWORK	IRVT = NRWORK
IRU = IRVT + M*M	IRU = IRVT + M*M
Line 2003	Line 2163
*	*
* Copy real matrix RWORK(IRU) to complex matrix U	* Copy real matrix RWORK(IRU) to complex matrix U
* Overwrite U by left singular vectors of A	* Overwrite U by left singular vectors of A
* (CWorkspace: need 3M, prefer 2M+M*NB)	* CWorkspace: need 2*M [tauq, taup] + M [work]
* (RWorkspace: M*M)	* CWorkspace: prefer 2M [tauq, taup] + MNB [work]
	* RWorkspace: need M [e] + MM [RVT] + MM [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 2017	Line 2178
*	*
* Copy real matrix RWORK(IRVT) to complex matrix VT	* Copy real matrix RWORK(IRVT) to complex matrix VT
* Overwrite VT by right singular vectors of A	* Overwrite VT by right singular vectors of A
* (CWorkspace: need 2M+N, prefer 2M+N*NB)	* CWorkspace: need 2*M [tauq, taup] + N [work]
* (RWorkspace: M*M)	* CWorkspace: prefer 2M [tauq, taup] + NNB [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,

CVSweb interface <joel.bertrand@systella.fr>

Removed from v.1.9
changed lines
	Added in v.1.17