rpl/lapack/lapack/dgesvd.f - diff

Return to dgesvd.f CVS log

Up to [local] / rpl / lapack / lapack

Diff for /rpl/lapack/lapack/dgesvd.f between versions 1.3 and 1.17

version 1.3, 2010/08/06 15:28:37	version 1.17, 2017/06/17 10:53:49
Line 1	Line 1
SUBROUTINE DGESVD( JOBU, JOBVT, M, N, A, LDA, S, U, LDU, VT, LDVT,	*> \brief <b> DGESVD computes the singular value decomposition (SVD) for GE matrices</b>
$ WORK, LWORK, INFO )
*	*
* -- LAPACK driver routine (version 3.2) --	* =========== DOCUMENTATION ===========
	*
	* Online html documentation available at
	* http://www.netlib.org/lapack/explore-html/
	*
	*> \htmlonly
	*> Download DGESVD + dependencies
	*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dgesvd.f">
	*> [TGZ]</a>
	*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dgesvd.f">
	*> [ZIP]</a>
	*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dgesvd.f">
	*> [TXT]</a>
	*> \endhtmlonly
	*
	* Definition:
	* ===========
	*
	* SUBROUTINE DGESVD( JOBU, JOBVT, M, N, A, LDA, S, U, LDU, VT, LDVT,
	* WORK, LWORK, INFO )
	*
	* .. Scalar Arguments ..
	* CHARACTER JOBU, JOBVT
	* INTEGER INFO, LDA, LDU, LDVT, LWORK, M, N
	* ..
	* .. Array Arguments ..
	* DOUBLE PRECISION A( LDA, * ), S( * ), U( LDU, * ),
	* $ VT( LDVT, * ), WORK( * )
	* ..
	*
	*
	*> \par Purpose:
	* =============
	*>
	*> \verbatim
	*>
	*> DGESVD computes the singular value decomposition (SVD) of a real
	*> M-by-N matrix A, optionally computing the left and/or right singular
	*> vectors. The SVD is written
	*>
	> A = U SIGMA * 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 orthogonal matrix, and
	*> V is an N-by-N orthogonal 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 V*T, not V.
	*> \endverbatim
	*
	* Arguments:
	* ==========
	*
	*> \param[in] JOBU
	*> \verbatim
	> JOBU is CHARACTER1
	*> Specifies options for computing all or part of the matrix U:
	*> = 'A': all M columns of U are returned in array U:
	*> = 'S': the first min(m,n) columns of U (the left singular
	*> vectors) are returned in the array U;
	*> = 'O': the first min(m,n) columns of U (the left singular
	*> vectors) are overwritten on the array A;
	*> = 'N': no columns of U (no left singular vectors) are
	*> computed.
	*> \endverbatim
	*>
	*> \param[in] JOBVT
	*> \verbatim
	> JOBVT is CHARACTER1
	*> Specifies options for computing all or part of the matrix
	> V*T:
	> = 'A': all N rows of V*T are returned in the array VT;
	> = 'S': the first min(m,n) rows of V*T (the right singular
	*> vectors) are returned in the array VT;
	> = 'O': the first min(m,n) rows of V*T (the right singular
	*> vectors) are overwritten on the array A;
	> = 'N': no rows of V*T (no right singular vectors) are
	*> computed.
	*>
	*> JOBVT and JOBU cannot both be 'O'.
	*> \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 DOUBLE PRECISION array, dimension (LDA,N)
	*> On entry, the M-by-N matrix A.
	*> On exit,
	*> if JOBU = 'O', A is overwritten with the first min(m,n)
	*> columns of U (the left singular vectors,
	*> stored columnwise);
	*> if JOBVT = 'O', A is overwritten with the first min(m,n)
	> rows of V*T (the right singular vectors,
	*> stored rowwise);
	*> if JOBU .ne. 'O' and JOBVT .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 DOUBLE PRECISION array, dimension (LDU,UCOL)
	*> (LDU,M) if JOBU = 'A' or (LDU,min(M,N)) if JOBU = 'S'.
	*> If JOBU = 'A', U contains the M-by-M orthogonal matrix U;
	*> if JOBU = 'S', U contains the first min(m,n) columns of U
	*> (the left singular vectors, stored columnwise);
	*> if JOBU = 'N' or 'O', U is not referenced.
	*> \endverbatim
	*>
	*> \param[in] LDU
	*> \verbatim
	*> LDU is INTEGER
	*> The leading dimension of the array U. LDU >= 1; if
	*> JOBU = 'S' or 'A', LDU >= M.
	*> \endverbatim
	*>
	*> \param[out] VT
	*> \verbatim
	*> VT is DOUBLE PRECISION array, dimension (LDVT,N)
	*> If JOBVT = 'A', VT contains the N-by-N orthogonal matrix
	> V*T;
	*> if JOBVT = 'S', VT contains the first min(m,n) rows of
	> V*T (the right singular vectors, stored rowwise);
	*> if JOBVT = 'N' or 'O', VT is not referenced.
	*> \endverbatim
	*>
	*> \param[in] LDVT
	*> \verbatim
	*> LDVT is INTEGER
	*> The leading dimension of the array VT. LDVT >= 1; if
	*> JOBVT = 'A', LDVT >= N; if JOBVT = 'S', LDVT >= min(M,N).
	*> \endverbatim
	*>
	*> \param[out] WORK
	*> \verbatim
	*> WORK is DOUBLE PRECISION array, dimension (MAX(1,LWORK))
	*> On exit, if INFO = 0, WORK(1) returns the optimal LWORK;
	*> if INFO > 0, WORK(2:MIN(M,N)) contains the unconverged
	*> superdiagonal elements of an upper bidiagonal matrix B
	*> whose diagonal is in S (not necessarily sorted). B
	> satisfies A = U B * VT, so it has the same singular values
	*> as A, and singular vectors related by U and VT.
	*> \endverbatim
	*>
	*> \param[in] LWORK
	*> \verbatim
	*> LWORK is INTEGER
	*> The dimension of the array WORK.
	> LWORK >= MAX(1,5MIN(M,N)) for the paths (see comments inside code):
	*> - PATH 1 (M much larger than N, JOBU='N')
	*> - PATH 1t (N much larger than M, JOBVT='N')
	> LWORK >= MAX(1,3MIN(M,N) + MAX(M,N),5*MIN(M,N)) for the other paths
	*> For good performance, LWORK should generally be larger.
	*>
	*> If LWORK = -1, then a workspace query is assumed; the routine
	*> only calculates the optimal size of the WORK array, returns
	*> this value as the first entry of the WORK array, and no error
	*> message related to LWORK is issued by XERBLA.
	*> \endverbatim
	*>
	*> \param[out] INFO
	*> \verbatim
	*> INFO is INTEGER
	*> = 0: successful exit.
	*> < 0: if INFO = -i, the i-th argument had an illegal value.
	*> > 0: if DBDSQR did not converge, INFO specifies how many
	*> superdiagonals of an intermediate bidiagonal form B
	*> did not converge to zero. See the description of WORK
	*> above for details.
	*> \endverbatim
	*
	* Authors:
	* ========
	*
	*> \author Univ. of Tennessee
	*> \author Univ. of California Berkeley
	*> \author Univ. of Colorado Denver
	*> \author NAG Ltd.
	*
	*> \date April 2012
	*
	*> \ingroup doubleGEsing
	*
	* =====================================================================
	SUBROUTINE DGESVD( JOBU, JOBVT, M, N, A, LDA, S, U, LDU,
	$ VT, LDVT, WORK, LWORK, INFO )
	*
	* -- LAPACK driver routine (version 3.7.0) --
* -- LAPACK is a software package provided by Univ. of Tennessee, --	* -- LAPACK is a software package provided by Univ. of Tennessee, --
* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--	* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
* November 2006	* April 2012
*	*
* .. Scalar Arguments ..	* .. Scalar Arguments ..
CHARACTER JOBU, JOBVT	CHARACTER JOBU, JOBVT
Line 15	Line 225
$ VT( LDVT, * ), WORK( * )	$ VT( LDVT, * ), WORK( * )
* ..	* ..
*	*
* Purpose
* =======
*
* DGESVD computes the singular value decomposition (SVD) of a real
* M-by-N matrix A, optionally computing the left and/or right singular
* vectors. The SVD is written
*
* A = U * SIGMA * 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 orthogonal matrix, and
* V is an N-by-N orthogonal 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 V**T, not V.
*
* Arguments
* =========
*
* JOBU (input) CHARACTER*1
* Specifies options for computing all or part of the matrix U:
* = 'A': all M columns of U are returned in array U:
* = 'S': the first min(m,n) columns of U (the left singular
* vectors) are returned in the array U;
* = 'O': the first min(m,n) columns of U (the left singular
* vectors) are overwritten on the array A;
* = 'N': no columns of U (no left singular vectors) are
* computed.
*
* JOBVT (input) CHARACTER*1
* Specifies options for computing all or part of the matrix
* V**T:
* = 'A': all N rows of V**T are returned in the array VT;
* = 'S': the first min(m,n) rows of V**T (the right singular
* vectors) are returned in the array VT;
* = 'O': the first min(m,n) rows of V**T (the right singular
* vectors) are overwritten on the array A;
* = 'N': no rows of V**T (no right singular vectors) are
* computed.
*
* JOBVT and JOBU cannot both be 'O'.
*
* 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) DOUBLE PRECISION array, dimension (LDA,N)
* On entry, the M-by-N matrix A.
* On exit,
* if JOBU = 'O', A is overwritten with the first min(m,n)
* columns of U (the left singular vectors,
* stored columnwise);
* if JOBVT = 'O', A is overwritten with the first min(m,n)
* rows of V**T (the right singular vectors,
* stored rowwise);
* if JOBU .ne. 'O' and JOBVT .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) DOUBLE PRECISION array, dimension (LDU,UCOL)
* (LDU,M) if JOBU = 'A' or (LDU,min(M,N)) if JOBU = 'S'.
* If JOBU = 'A', U contains the M-by-M orthogonal matrix U;
* if JOBU = 'S', U contains the first min(m,n) columns of U
* (the left singular vectors, stored columnwise);
* if JOBU = 'N' or 'O', U is not referenced.
*
* LDU (input) INTEGER
* The leading dimension of the array U. LDU >= 1; if
* JOBU = 'S' or 'A', LDU >= M.
*
* VT (output) DOUBLE PRECISION array, dimension (LDVT,N)
* If JOBVT = 'A', VT contains the N-by-N orthogonal matrix
* V**T;
* if JOBVT = 'S', VT contains the first min(m,n) rows of
* V**T (the right singular vectors, stored rowwise);
* if JOBVT = 'N' or 'O', VT is not referenced.
*
* LDVT (input) INTEGER
* The leading dimension of the array VT. LDVT >= 1; if
* JOBVT = 'A', LDVT >= N; if JOBVT = 'S', LDVT >= min(M,N).
*
* WORK (workspace/output) DOUBLE PRECISION array, dimension (MAX(1,LWORK))
* On exit, if INFO = 0, WORK(1) returns the optimal LWORK;
* if INFO > 0, WORK(2:MIN(M,N)) contains the unconverged
* superdiagonal elements of an upper bidiagonal matrix B
* whose diagonal is in S (not necessarily sorted). B
* satisfies A = U * B * VT, so it has the same singular values
* as A, and singular vectors related by U and VT.
*
* LWORK (input) INTEGER
* The dimension of the array WORK.
* LWORK >= MAX(1,3MIN(M,N)+MAX(M,N),5MIN(M,N)).
* For good performance, LWORK should generally be larger.
*
* If LWORK = -1, then a workspace query is assumed; the routine
* only calculates the optimal size of the WORK array, returns
* this value as the first entry of the WORK array, and no error
* message related to LWORK is issued by XERBLA.
*
* INFO (output) INTEGER
* = 0: successful exit.
* < 0: if INFO = -i, the i-th argument had an illegal value.
* > 0: if DBDSQR did not converge, INFO specifies how many
* superdiagonals of an intermediate bidiagonal form B
* did not converge to zero. See the description of WORK
* above for details.
*
* =====================================================================	* =====================================================================
*	*
* .. Parameters ..	* .. Parameters ..
Line 144	Line 238
$ ITAU, ITAUP, ITAUQ, IU, IWORK, LDWRKR, LDWRKU,	$ ITAU, ITAUP, ITAUQ, IU, IWORK, LDWRKR, LDWRKU,
$ MAXWRK, MINMN, MINWRK, MNTHR, NCU, NCVT, NRU,	$ MAXWRK, MINMN, MINWRK, MNTHR, NCU, NCVT, NRU,
$ NRVT, WRKBL	$ NRVT, WRKBL
	INTEGER LWORK_DGEQRF, LWORK_DORGQR_N, LWORK_DORGQR_M,
	$ LWORK_DGEBRD, LWORK_DORGBR_P, LWORK_DORGBR_Q,
	$ LWORK_DGELQF, LWORK_DORGLQ_N, LWORK_DORGLQ_M
DOUBLE PRECISION ANRM, BIGNUM, EPS, SMLNUM	DOUBLE PRECISION ANRM, BIGNUM, EPS, SMLNUM
* ..	* ..
* .. Local Arrays ..	* .. Local Arrays ..
Line 215	Line 312
*	*
MNTHR = ILAENV( 6, 'DGESVD', JOBU // JOBVT, M, N, 0, 0 )	MNTHR = ILAENV( 6, 'DGESVD', JOBU // JOBVT, M, N, 0, 0 )
BDSPAC = 5*N	BDSPAC = 5*N
	* Compute space needed for DGEQRF
	CALL DGEQRF( M, N, A, LDA, DUM(1), DUM(1), -1, IERR )
	LWORK_DGEQRF = INT( DUM(1) )
	* Compute space needed for DORGQR
	CALL DORGQR( M, N, N, A, LDA, DUM(1), DUM(1), -1, IERR )
	LWORK_DORGQR_N = INT( DUM(1) )
	CALL DORGQR( M, M, N, A, LDA, DUM(1), DUM(1), -1, IERR )
	LWORK_DORGQR_M = INT( DUM(1) )
	* Compute space needed for DGEBRD
	CALL DGEBRD( N, N, A, LDA, S, DUM(1), DUM(1),
	$ DUM(1), DUM(1), -1, IERR )
	LWORK_DGEBRD = INT( DUM(1) )
	* Compute space needed for DORGBR P
	CALL DORGBR( 'P', N, N, N, A, LDA, DUM(1),
	$ DUM(1), -1, IERR )
	LWORK_DORGBR_P = INT( DUM(1) )
	* Compute space needed for DORGBR Q
	CALL DORGBR( 'Q', N, N, N, A, LDA, DUM(1),
	$ DUM(1), -1, IERR )
	LWORK_DORGBR_Q = INT( DUM(1) )
	*
IF( M.GE.MNTHR ) THEN	IF( M.GE.MNTHR ) THEN
IF( WNTUN ) THEN	IF( WNTUN ) THEN
*	*
* Path 1 (M much larger than N, JOBU='N')	* Path 1 (M much larger than N, JOBU='N')
*	*
MAXWRK = N + N*ILAENV( 1, 'DGEQRF', ' ', M, N, -1,	MAXWRK = N + LWORK_DGEQRF
$ -1 )	MAXWRK = MAX( MAXWRK, 3*N + LWORK_DGEBRD )
MAXWRK = MAX( MAXWRK, 3N+2N*
$ ILAENV( 1, 'DGEBRD', ' ', N, N, -1, -1 ) )
IF( WNTVO .OR. WNTVAS )	IF( WNTVO .OR. WNTVAS )
$ MAXWRK = MAX( MAXWRK, 3N+( N-1 )	$ MAXWRK = MAX( MAXWRK, 3*N + LWORK_DORGBR_P )
$ ILAENV( 1, 'DORGBR', 'P', N, N, N, -1 ) )
MAXWRK = MAX( MAXWRK, BDSPAC )	MAXWRK = MAX( MAXWRK, BDSPAC )
MINWRK = MAX( 4*N, BDSPAC )	MINWRK = MAX( 4*N, BDSPAC )
ELSE IF( WNTUO .AND. WNTVN ) THEN	ELSE IF( WNTUO .AND. WNTVN ) THEN
*	*
* Path 2 (M much larger than N, JOBU='O', JOBVT='N')	* Path 2 (M much larger than N, JOBU='O', JOBVT='N')
*	*
WRKBL = N + N*ILAENV( 1, 'DGEQRF', ' ', M, N, -1, -1 )	WRKBL = N + LWORK_DGEQRF
WRKBL = MAX( WRKBL, N+N*ILAENV( 1, 'DORGQR', ' ', M,	WRKBL = MAX( WRKBL, N + LWORK_DORGQR_N )
$ N, N, -1 ) )	WRKBL = MAX( WRKBL, 3*N + LWORK_DGEBRD )
WRKBL = MAX( WRKBL, 3N+2N*	WRKBL = MAX( WRKBL, 3*N + LWORK_DORGBR_Q )
$ ILAENV( 1, 'DGEBRD', ' ', N, N, -1, -1 ) )
WRKBL = MAX( WRKBL, 3N+N
$ ILAENV( 1, 'DORGBR', 'Q', N, N, N, -1 ) )
WRKBL = MAX( WRKBL, BDSPAC )	WRKBL = MAX( WRKBL, BDSPAC )
MAXWRK = MAX( NN+WRKBL, NN+M*N+N )	MAXWRK = MAX( NN + WRKBL, NN + M*N + N )
MINWRK = MAX( 3*N+M, BDSPAC )	MINWRK = MAX( 3*N + M, BDSPAC )
ELSE IF( WNTUO .AND. WNTVAS ) THEN	ELSE IF( WNTUO .AND. WNTVAS ) THEN
*	*
* Path 3 (M much larger than N, JOBU='O', JOBVT='S' or	* Path 3 (M much larger than N, JOBU='O', JOBVT='S' or
* 'A')	* 'A')
*	*
WRKBL = N + N*ILAENV( 1, 'DGEQRF', ' ', M, N, -1, -1 )	WRKBL = N + LWORK_DGEQRF
WRKBL = MAX( WRKBL, N+N*ILAENV( 1, 'DORGQR', ' ', M,	WRKBL = MAX( WRKBL, N + LWORK_DORGQR_N )
$ N, N, -1 ) )	WRKBL = MAX( WRKBL, 3*N + LWORK_DGEBRD )
WRKBL = MAX( WRKBL, 3N+2N*	WRKBL = MAX( WRKBL, 3*N + LWORK_DORGBR_Q )
$ ILAENV( 1, 'DGEBRD', ' ', N, N, -1, -1 ) )	WRKBL = MAX( WRKBL, 3*N + LWORK_DORGBR_P )
WRKBL = MAX( WRKBL, 3N+N
$ ILAENV( 1, 'DORGBR', 'Q', N, N, N, -1 ) )
WRKBL = MAX( WRKBL, 3N+( N-1 )
$ ILAENV( 1, 'DORGBR', 'P', N, N, N, -1 ) )
WRKBL = MAX( WRKBL, BDSPAC )	WRKBL = MAX( WRKBL, BDSPAC )
MAXWRK = MAX( NN+WRKBL, NN+M*N+N )	MAXWRK = MAX( NN + WRKBL, NN + M*N + N )
MINWRK = MAX( 3*N+M, BDSPAC )	MINWRK = MAX( 3*N + M, BDSPAC )
ELSE IF( WNTUS .AND. WNTVN ) THEN	ELSE IF( WNTUS .AND. WNTVN ) THEN
*	*
* Path 4 (M much larger than N, JOBU='S', JOBVT='N')	* Path 4 (M much larger than N, JOBU='S', JOBVT='N')
*	*
WRKBL = N + N*ILAENV( 1, 'DGEQRF', ' ', M, N, -1, -1 )	WRKBL = N + LWORK_DGEQRF
WRKBL = MAX( WRKBL, N+N*ILAENV( 1, 'DORGQR', ' ', M,	WRKBL = MAX( WRKBL, N + LWORK_DORGQR_N )
$ N, N, -1 ) )	WRKBL = MAX( WRKBL, 3*N + LWORK_DGEBRD )
WRKBL = MAX( WRKBL, 3N+2N*	WRKBL = MAX( WRKBL, 3*N + LWORK_DORGBR_Q )
$ ILAENV( 1, 'DGEBRD', ' ', N, N, -1, -1 ) )
WRKBL = MAX( WRKBL, 3N+N
$ ILAENV( 1, 'DORGBR', 'Q', N, N, N, -1 ) )
WRKBL = MAX( WRKBL, BDSPAC )	WRKBL = MAX( WRKBL, BDSPAC )
MAXWRK = N*N + WRKBL	MAXWRK = N*N + WRKBL
MINWRK = MAX( 3*N+M, BDSPAC )	MINWRK = MAX( 3*N + M, BDSPAC )
ELSE IF( WNTUS .AND. WNTVO ) THEN	ELSE IF( WNTUS .AND. WNTVO ) THEN
*	*
* Path 5 (M much larger than N, JOBU='S', JOBVT='O')	* Path 5 (M much larger than N, JOBU='S', JOBVT='O')
*	*
WRKBL = N + N*ILAENV( 1, 'DGEQRF', ' ', M, N, -1, -1 )	WRKBL = N + LWORK_DGEQRF
WRKBL = MAX( WRKBL, N+N*ILAENV( 1, 'DORGQR', ' ', M,	WRKBL = MAX( WRKBL, N + LWORK_DORGQR_N )
$ N, N, -1 ) )	WRKBL = MAX( WRKBL, 3*N + LWORK_DGEBRD )
WRKBL = MAX( WRKBL, 3N+2N*	WRKBL = MAX( WRKBL, 3*N + LWORK_DORGBR_Q )
$ ILAENV( 1, 'DGEBRD', ' ', N, N, -1, -1 ) )	WRKBL = MAX( WRKBL, 3*N + LWORK_DORGBR_P )
WRKBL = MAX( WRKBL, 3N+N
$ ILAENV( 1, 'DORGBR', 'Q', N, N, N, -1 ) )
WRKBL = MAX( WRKBL, 3N+( N-1 )
$ ILAENV( 1, 'DORGBR', 'P', N, N, N, -1 ) )
WRKBL = MAX( WRKBL, BDSPAC )	WRKBL = MAX( WRKBL, BDSPAC )
MAXWRK = 2NN + WRKBL	MAXWRK = 2NN + WRKBL
MINWRK = MAX( 3*N+M, BDSPAC )	MINWRK = MAX( 3*N + M, BDSPAC )
ELSE IF( WNTUS .AND. WNTVAS ) THEN	ELSE IF( WNTUS .AND. WNTVAS ) THEN
*	*
* Path 6 (M much larger than N, JOBU='S', JOBVT='S' or	* Path 6 (M much larger than N, JOBU='S', JOBVT='S' or
* 'A')	* 'A')
*	*
WRKBL = N + N*ILAENV( 1, 'DGEQRF', ' ', M, N, -1, -1 )	WRKBL = N + LWORK_DGEQRF
WRKBL = MAX( WRKBL, N+N*ILAENV( 1, 'DORGQR', ' ', M,	WRKBL = MAX( WRKBL, N + LWORK_DORGQR_N )
$ N, N, -1 ) )	WRKBL = MAX( WRKBL, 3*N + LWORK_DGEBRD )
WRKBL = MAX( WRKBL, 3N+2N*	WRKBL = MAX( WRKBL, 3*N + LWORK_DORGBR_Q )
$ ILAENV( 1, 'DGEBRD', ' ', N, N, -1, -1 ) )	WRKBL = MAX( WRKBL, 3*N + LWORK_DORGBR_P )
WRKBL = MAX( WRKBL, 3N+N
$ ILAENV( 1, 'DORGBR', 'Q', N, N, N, -1 ) )
WRKBL = MAX( WRKBL, 3N+( N-1 )
$ ILAENV( 1, 'DORGBR', 'P', N, N, N, -1 ) )
WRKBL = MAX( WRKBL, BDSPAC )	WRKBL = MAX( WRKBL, BDSPAC )
MAXWRK = N*N + WRKBL	MAXWRK = N*N + WRKBL
MINWRK = MAX( 3*N+M, BDSPAC )	MINWRK = MAX( 3*N + M, BDSPAC )
ELSE IF( WNTUA .AND. WNTVN ) THEN	ELSE IF( WNTUA .AND. WNTVN ) THEN
*	*
* Path 7 (M much larger than N, JOBU='A', JOBVT='N')	* Path 7 (M much larger than N, JOBU='A', JOBVT='N')
*	*
WRKBL = N + N*ILAENV( 1, 'DGEQRF', ' ', M, N, -1, -1 )	WRKBL = N + LWORK_DGEQRF
WRKBL = MAX( WRKBL, N+M*ILAENV( 1, 'DORGQR', ' ', M,	WRKBL = MAX( WRKBL, N + LWORK_DORGQR_M )
$ M, N, -1 ) )	WRKBL = MAX( WRKBL, 3*N + LWORK_DGEBRD )
WRKBL = MAX( WRKBL, 3N+2N*	WRKBL = MAX( WRKBL, 3*N + LWORK_DORGBR_Q )
$ ILAENV( 1, 'DGEBRD', ' ', N, N, -1, -1 ) )
WRKBL = MAX( WRKBL, 3N+N
$ ILAENV( 1, 'DORGBR', 'Q', N, N, N, -1 ) )
WRKBL = MAX( WRKBL, BDSPAC )	WRKBL = MAX( WRKBL, BDSPAC )
MAXWRK = N*N + WRKBL	MAXWRK = N*N + WRKBL
MINWRK = MAX( 3*N+M, BDSPAC )	MINWRK = MAX( 3*N + M, BDSPAC )
ELSE IF( WNTUA .AND. WNTVO ) THEN	ELSE IF( WNTUA .AND. WNTVO ) THEN
*	*
* Path 8 (M much larger than N, JOBU='A', JOBVT='O')	* Path 8 (M much larger than N, JOBU='A', JOBVT='O')
*	*
WRKBL = N + N*ILAENV( 1, 'DGEQRF', ' ', M, N, -1, -1 )	WRKBL = N + LWORK_DGEQRF
WRKBL = MAX( WRKBL, N+M*ILAENV( 1, 'DORGQR', ' ', M,	WRKBL = MAX( WRKBL, N + LWORK_DORGQR_M )
$ M, N, -1 ) )	WRKBL = MAX( WRKBL, 3*N + LWORK_DGEBRD )
WRKBL = MAX( WRKBL, 3N+2N*	WRKBL = MAX( WRKBL, 3*N + LWORK_DORGBR_Q )
$ ILAENV( 1, 'DGEBRD', ' ', N, N, -1, -1 ) )	WRKBL = MAX( WRKBL, 3*N + LWORK_DORGBR_P )
WRKBL = MAX( WRKBL, 3N+N
$ ILAENV( 1, 'DORGBR', 'Q', N, N, N, -1 ) )
WRKBL = MAX( WRKBL, 3N+( N-1 )
$ ILAENV( 1, 'DORGBR', 'P', N, N, N, -1 ) )
WRKBL = MAX( WRKBL, BDSPAC )	WRKBL = MAX( WRKBL, BDSPAC )
MAXWRK = 2NN + WRKBL	MAXWRK = 2NN + WRKBL
MINWRK = MAX( 3*N+M, BDSPAC )	MINWRK = MAX( 3*N + M, BDSPAC )
ELSE IF( WNTUA .AND. WNTVAS ) THEN	ELSE IF( WNTUA .AND. WNTVAS ) THEN
*	*
* Path 9 (M much larger than N, JOBU='A', JOBVT='S' or	* Path 9 (M much larger than N, JOBU='A', JOBVT='S' or
* 'A')	* 'A')
*	*
WRKBL = N + N*ILAENV( 1, 'DGEQRF', ' ', M, N, -1, -1 )	WRKBL = N + LWORK_DGEQRF
WRKBL = MAX( WRKBL, N+M*ILAENV( 1, 'DORGQR', ' ', M,	WRKBL = MAX( WRKBL, N + LWORK_DORGQR_M )
$ M, N, -1 ) )	WRKBL = MAX( WRKBL, 3*N + LWORK_DGEBRD )
WRKBL = MAX( WRKBL, 3N+2N*	WRKBL = MAX( WRKBL, 3*N + LWORK_DORGBR_Q )
$ ILAENV( 1, 'DGEBRD', ' ', N, N, -1, -1 ) )	WRKBL = MAX( WRKBL, 3*N + LWORK_DORGBR_P )
WRKBL = MAX( WRKBL, 3N+N
$ ILAENV( 1, 'DORGBR', 'Q', N, N, N, -1 ) )
WRKBL = MAX( WRKBL, 3N+( N-1 )
$ ILAENV( 1, 'DORGBR', 'P', N, N, N, -1 ) )
WRKBL = MAX( WRKBL, BDSPAC )	WRKBL = MAX( WRKBL, BDSPAC )
MAXWRK = N*N + WRKBL	MAXWRK = N*N + WRKBL
MINWRK = MAX( 3*N+M, BDSPAC )	MINWRK = MAX( 3*N + M, BDSPAC )
END IF	END IF
ELSE	ELSE
*	*
* Path 10 (M at least N, but not much larger)	* Path 10 (M at least N, but not much larger)
*	*
MAXWRK = 3N + ( M+N )ILAENV( 1, 'DGEBRD', ' ', M, N,	CALL DGEBRD( M, N, A, LDA, S, DUM(1), DUM(1),
$ -1, -1 )	$ DUM(1), DUM(1), -1, IERR )
IF( WNTUS .OR. WNTUO )	LWORK_DGEBRD = INT( DUM(1) )
$ MAXWRK = MAX( MAXWRK, 3N+N	MAXWRK = 3*N + LWORK_DGEBRD
$ ILAENV( 1, 'DORGBR', 'Q', M, N, N, -1 ) )	IF( WNTUS .OR. WNTUO ) THEN
IF( WNTUA )	CALL DORGBR( 'Q', M, N, N, A, LDA, DUM(1),
$ MAXWRK = MAX( MAXWRK, 3N+M	$ DUM(1), -1, IERR )
$ ILAENV( 1, 'DORGBR', 'Q', M, M, N, -1 ) )	LWORK_DORGBR_Q = INT( DUM(1) )
IF( .NOT.WNTVN )	MAXWRK = MAX( MAXWRK, 3*N + LWORK_DORGBR_Q )
$ MAXWRK = MAX( MAXWRK, 3N+( N-1 )	END IF
$ ILAENV( 1, 'DORGBR', 'P', N, N, N, -1 ) )	IF( WNTUA ) THEN
	CALL DORGBR( 'Q', M, M, N, A, LDA, DUM(1),
	$ DUM(1), -1, IERR )
	LWORK_DORGBR_Q = INT( DUM(1) )
	MAXWRK = MAX( MAXWRK, 3*N + LWORK_DORGBR_Q )
	END IF
	IF( .NOT.WNTVN ) THEN
	MAXWRK = MAX( MAXWRK, 3*N + LWORK_DORGBR_P )
	END IF
MAXWRK = MAX( MAXWRK, BDSPAC )	MAXWRK = MAX( MAXWRK, BDSPAC )
MINWRK = MAX( 3*N+M, BDSPAC )	MINWRK = MAX( 3*N + M, BDSPAC )
END IF	END IF
ELSE IF( MINMN.GT.0 ) THEN	ELSE IF( MINMN.GT.0 ) THEN
*	*
Line 379	Line 473
*	*
MNTHR = ILAENV( 6, 'DGESVD', JOBU // JOBVT, M, N, 0, 0 )	MNTHR = ILAENV( 6, 'DGESVD', JOBU // JOBVT, M, N, 0, 0 )
BDSPAC = 5*M	BDSPAC = 5*M
	* Compute space needed for DGELQF
	CALL DGELQF( M, N, A, LDA, DUM(1), DUM(1), -1, IERR )
	LWORK_DGELQF = INT( DUM(1) )
	* Compute space needed for DORGLQ
	CALL DORGLQ( N, N, M, DUM(1), N, DUM(1), DUM(1), -1, IERR )
	LWORK_DORGLQ_N = INT( DUM(1) )
	CALL DORGLQ( M, N, M, A, LDA, DUM(1), DUM(1), -1, IERR )
	LWORK_DORGLQ_M = INT( DUM(1) )
	* Compute space needed for DGEBRD
	CALL DGEBRD( M, M, A, LDA, S, DUM(1), DUM(1),
	$ DUM(1), DUM(1), -1, IERR )
	LWORK_DGEBRD = INT( DUM(1) )
	* Compute space needed for DORGBR P
	CALL DORGBR( 'P', M, M, M, A, N, DUM(1),
	$ DUM(1), -1, IERR )
	LWORK_DORGBR_P = INT( DUM(1) )
	* Compute space needed for DORGBR Q
	CALL DORGBR( 'Q', M, M, M, A, N, DUM(1),
	$ DUM(1), -1, IERR )
	LWORK_DORGBR_Q = INT( DUM(1) )
IF( N.GE.MNTHR ) THEN	IF( N.GE.MNTHR ) THEN
IF( WNTVN ) THEN	IF( WNTVN ) THEN
*	*
* Path 1t(N much larger than M, JOBVT='N')	* Path 1t(N much larger than M, JOBVT='N')
*	*
MAXWRK = M + M*ILAENV( 1, 'DGELQF', ' ', M, N, -1,	MAXWRK = M + LWORK_DGELQF
$ -1 )	MAXWRK = MAX( MAXWRK, 3*M + LWORK_DGEBRD )
MAXWRK = MAX( MAXWRK, 3M+2M*
$ ILAENV( 1, 'DGEBRD', ' ', M, M, -1, -1 ) )
IF( WNTUO .OR. WNTUAS )	IF( WNTUO .OR. WNTUAS )
$ MAXWRK = MAX( MAXWRK, 3M+M	$ MAXWRK = MAX( MAXWRK, 3*M + LWORK_DORGBR_Q )
$ ILAENV( 1, 'DORGBR', 'Q', M, M, M, -1 ) )
MAXWRK = MAX( MAXWRK, BDSPAC )	MAXWRK = MAX( MAXWRK, BDSPAC )
MINWRK = MAX( 4*M, BDSPAC )	MINWRK = MAX( 4*M, BDSPAC )
ELSE IF( WNTVO .AND. WNTUN ) THEN	ELSE IF( WNTVO .AND. WNTUN ) THEN
*	*
* Path 2t(N much larger than M, JOBU='N', JOBVT='O')	* Path 2t(N much larger than M, JOBU='N', JOBVT='O')
*	*
WRKBL = M + M*ILAENV( 1, 'DGELQF', ' ', M, N, -1, -1 )	WRKBL = M + LWORK_DGELQF
WRKBL = MAX( WRKBL, M+M*ILAENV( 1, 'DORGLQ', ' ', M,	WRKBL = MAX( WRKBL, M + LWORK_DORGLQ_M )
$ N, M, -1 ) )	WRKBL = MAX( WRKBL, 3*M + LWORK_DGEBRD )
WRKBL = MAX( WRKBL, 3M+2M*	WRKBL = MAX( WRKBL, 3*M + LWORK_DORGBR_P )
$ ILAENV( 1, 'DGEBRD', ' ', M, M, -1, -1 ) )
WRKBL = MAX( WRKBL, 3M+( M-1 )
$ ILAENV( 1, 'DORGBR', 'P', M, M, M, -1 ) )
WRKBL = MAX( WRKBL, BDSPAC )	WRKBL = MAX( WRKBL, BDSPAC )
MAXWRK = MAX( MM+WRKBL, MM+M*N+M )	MAXWRK = MAX( MM + WRKBL, MM + M*N + M )
MINWRK = MAX( 3*M+N, BDSPAC )	MINWRK = MAX( 3*M + N, BDSPAC )
ELSE IF( WNTVO .AND. WNTUAS ) THEN	ELSE IF( WNTVO .AND. WNTUAS ) THEN
*	*
* Path 3t(N much larger than M, JOBU='S' or 'A',	* Path 3t(N much larger than M, JOBU='S' or 'A',
* JOBVT='O')	* JOBVT='O')
*	*
WRKBL = M + M*ILAENV( 1, 'DGELQF', ' ', M, N, -1, -1 )	WRKBL = M + LWORK_DGELQF
WRKBL = MAX( WRKBL, M+M*ILAENV( 1, 'DORGLQ', ' ', M,	WRKBL = MAX( WRKBL, M + LWORK_DORGLQ_M )
$ N, M, -1 ) )	WRKBL = MAX( WRKBL, 3*M + LWORK_DGEBRD )
WRKBL = MAX( WRKBL, 3M+2M*	WRKBL = MAX( WRKBL, 3*M + LWORK_DORGBR_P )
$ ILAENV( 1, 'DGEBRD', ' ', M, M, -1, -1 ) )	WRKBL = MAX( WRKBL, 3*M + LWORK_DORGBR_Q )
WRKBL = MAX( WRKBL, 3M+( M-1 )
$ ILAENV( 1, 'DORGBR', 'P', M, M, M, -1 ) )
WRKBL = MAX( WRKBL, 3M+M
$ ILAENV( 1, 'DORGBR', 'Q', M, M, M, -1 ) )
WRKBL = MAX( WRKBL, BDSPAC )	WRKBL = MAX( WRKBL, BDSPAC )
MAXWRK = MAX( MM+WRKBL, MM+M*N+M )	MAXWRK = MAX( MM + WRKBL, MM + M*N + M )
MINWRK = MAX( 3*M+N, BDSPAC )	MINWRK = MAX( 3*M + N, BDSPAC )
ELSE IF( WNTVS .AND. WNTUN ) THEN	ELSE IF( WNTVS .AND. WNTUN ) THEN
*	*
* Path 4t(N much larger than M, JOBU='N', JOBVT='S')	* Path 4t(N much larger than M, JOBU='N', JOBVT='S')
*	*
WRKBL = M + M*ILAENV( 1, 'DGELQF', ' ', M, N, -1, -1 )	WRKBL = M + LWORK_DGELQF
WRKBL = MAX( WRKBL, M+M*ILAENV( 1, 'DORGLQ', ' ', M,	WRKBL = MAX( WRKBL, M + LWORK_DORGLQ_M )
$ N, M, -1 ) )	WRKBL = MAX( WRKBL, 3*M + LWORK_DGEBRD )
WRKBL = MAX( WRKBL, 3M+2M*	WRKBL = MAX( WRKBL, 3*M + LWORK_DORGBR_P )
$ ILAENV( 1, 'DGEBRD', ' ', M, M, -1, -1 ) )
WRKBL = MAX( WRKBL, 3M+( M-1 )
$ ILAENV( 1, 'DORGBR', 'P', M, M, M, -1 ) )
WRKBL = MAX( WRKBL, BDSPAC )	WRKBL = MAX( WRKBL, BDSPAC )
MAXWRK = M*M + WRKBL	MAXWRK = M*M + WRKBL
MINWRK = MAX( 3*M+N, BDSPAC )	MINWRK = MAX( 3*M + N, BDSPAC )
ELSE IF( WNTVS .AND. WNTUO ) THEN	ELSE IF( WNTVS .AND. WNTUO ) THEN
*	*
* Path 5t(N much larger than M, JOBU='O', JOBVT='S')	* Path 5t(N much larger than M, JOBU='O', JOBVT='S')
*	*
WRKBL = M + M*ILAENV( 1, 'DGELQF', ' ', M, N, -1, -1 )	WRKBL = M + LWORK_DGELQF
WRKBL = MAX( WRKBL, M+M*ILAENV( 1, 'DORGLQ', ' ', M,	WRKBL = MAX( WRKBL, M + LWORK_DORGLQ_M )
$ N, M, -1 ) )	WRKBL = MAX( WRKBL, 3*M + LWORK_DGEBRD )
WRKBL = MAX( WRKBL, 3M+2M*	WRKBL = MAX( WRKBL, 3*M + LWORK_DORGBR_P )
$ ILAENV( 1, 'DGEBRD', ' ', M, M, -1, -1 ) )	WRKBL = MAX( WRKBL, 3*M + LWORK_DORGBR_Q )
WRKBL = MAX( WRKBL, 3M+( M-1 )
$ ILAENV( 1, 'DORGBR', 'P', M, M, M, -1 ) )
WRKBL = MAX( WRKBL, 3M+M
$ ILAENV( 1, 'DORGBR', 'Q', M, M, M, -1 ) )
WRKBL = MAX( WRKBL, BDSPAC )	WRKBL = MAX( WRKBL, BDSPAC )
MAXWRK = 2MM + WRKBL	MAXWRK = 2MM + WRKBL
MINWRK = MAX( 3*M+N, BDSPAC )	MINWRK = MAX( 3*M + N, BDSPAC )
ELSE IF( WNTVS .AND. WNTUAS ) THEN	ELSE IF( WNTVS .AND. WNTUAS ) THEN
*	*
* Path 6t(N much larger than M, JOBU='S' or 'A',	* Path 6t(N much larger than M, JOBU='S' or 'A',
* JOBVT='S')	* JOBVT='S')
*	*
WRKBL = M + M*ILAENV( 1, 'DGELQF', ' ', M, N, -1, -1 )	WRKBL = M + LWORK_DGELQF
WRKBL = MAX( WRKBL, M+M*ILAENV( 1, 'DORGLQ', ' ', M,	WRKBL = MAX( WRKBL, M + LWORK_DORGLQ_M )
$ N, M, -1 ) )	WRKBL = MAX( WRKBL, 3*M + LWORK_DGEBRD )
WRKBL = MAX( WRKBL, 3M+2M*	WRKBL = MAX( WRKBL, 3*M + LWORK_DORGBR_P )
$ ILAENV( 1, 'DGEBRD', ' ', M, M, -1, -1 ) )	WRKBL = MAX( WRKBL, 3*M + LWORK_DORGBR_Q )
WRKBL = MAX( WRKBL, 3M+( M-1 )
$ ILAENV( 1, 'DORGBR', 'P', M, M, M, -1 ) )
WRKBL = MAX( WRKBL, 3M+M
$ ILAENV( 1, 'DORGBR', 'Q', M, M, M, -1 ) )
WRKBL = MAX( WRKBL, BDSPAC )	WRKBL = MAX( WRKBL, BDSPAC )
MAXWRK = M*M + WRKBL	MAXWRK = M*M + WRKBL
MINWRK = MAX( 3*M+N, BDSPAC )	MINWRK = MAX( 3*M + N, BDSPAC )
ELSE IF( WNTVA .AND. WNTUN ) THEN	ELSE IF( WNTVA .AND. WNTUN ) THEN
*	*
* Path 7t(N much larger than M, JOBU='N', JOBVT='A')	* Path 7t(N much larger than M, JOBU='N', JOBVT='A')
*	*
WRKBL = M + M*ILAENV( 1, 'DGELQF', ' ', M, N, -1, -1 )	WRKBL = M + LWORK_DGELQF
WRKBL = MAX( WRKBL, M+N*ILAENV( 1, 'DORGLQ', ' ', N,	WRKBL = MAX( WRKBL, M + LWORK_DORGLQ_N )
$ N, M, -1 ) )	WRKBL = MAX( WRKBL, 3*M + LWORK_DGEBRD )
WRKBL = MAX( WRKBL, 3M+2M*	WRKBL = MAX( WRKBL, 3*M + LWORK_DORGBR_P )
$ ILAENV( 1, 'DGEBRD', ' ', M, M, -1, -1 ) )
WRKBL = MAX( WRKBL, 3M+( M-1 )
$ ILAENV( 1, 'DORGBR', 'P', M, M, M, -1 ) )
WRKBL = MAX( WRKBL, BDSPAC )	WRKBL = MAX( WRKBL, BDSPAC )
MAXWRK = M*M + WRKBL	MAXWRK = M*M + WRKBL
MINWRK = MAX( 3*M+N, BDSPAC )	MINWRK = MAX( 3*M + N, BDSPAC )
ELSE IF( WNTVA .AND. WNTUO ) THEN	ELSE IF( WNTVA .AND. WNTUO ) THEN
*	*
* Path 8t(N much larger than M, JOBU='O', JOBVT='A')	* Path 8t(N much larger than M, JOBU='O', JOBVT='A')
*	*
WRKBL = M + M*ILAENV( 1, 'DGELQF', ' ', M, N, -1, -1 )	WRKBL = M + LWORK_DGELQF
WRKBL = MAX( WRKBL, M+N*ILAENV( 1, 'DORGLQ', ' ', N,	WRKBL = MAX( WRKBL, M + LWORK_DORGLQ_N )
$ N, M, -1 ) )	WRKBL = MAX( WRKBL, 3*M + LWORK_DGEBRD )
WRKBL = MAX( WRKBL, 3M+2M*	WRKBL = MAX( WRKBL, 3*M + LWORK_DORGBR_P )
$ ILAENV( 1, 'DGEBRD', ' ', M, M, -1, -1 ) )	WRKBL = MAX( WRKBL, 3*M + LWORK_DORGBR_Q )
WRKBL = MAX( WRKBL, 3M+( M-1 )
$ ILAENV( 1, 'DORGBR', 'P', M, M, M, -1 ) )
WRKBL = MAX( WRKBL, 3M+M
$ ILAENV( 1, 'DORGBR', 'Q', M, M, M, -1 ) )
WRKBL = MAX( WRKBL, BDSPAC )	WRKBL = MAX( WRKBL, BDSPAC )
MAXWRK = 2MM + WRKBL	MAXWRK = 2MM + WRKBL
MINWRK = MAX( 3*M+N, BDSPAC )	MINWRK = MAX( 3*M + N, BDSPAC )
ELSE IF( WNTVA .AND. WNTUAS ) THEN	ELSE IF( WNTVA .AND. WNTUAS ) THEN
*	*
* Path 9t(N much larger than M, JOBU='S' or 'A',	* Path 9t(N much larger than M, JOBU='S' or 'A',
* JOBVT='A')	* JOBVT='A')
*	*
WRKBL = M + M*ILAENV( 1, 'DGELQF', ' ', M, N, -1, -1 )	WRKBL = M + LWORK_DGELQF
WRKBL = MAX( WRKBL, M+N*ILAENV( 1, 'DORGLQ', ' ', N,	WRKBL = MAX( WRKBL, M + LWORK_DORGLQ_N )
$ N, M, -1 ) )	WRKBL = MAX( WRKBL, 3*M + LWORK_DGEBRD )
WRKBL = MAX( WRKBL, 3M+2M*	WRKBL = MAX( WRKBL, 3*M + LWORK_DORGBR_P )
$ ILAENV( 1, 'DGEBRD', ' ', M, M, -1, -1 ) )	WRKBL = MAX( WRKBL, 3*M + LWORK_DORGBR_Q )
WRKBL = MAX( WRKBL, 3M+( M-1 )
$ ILAENV( 1, 'DORGBR', 'P', M, M, M, -1 ) )
WRKBL = MAX( WRKBL, 3M+M
$ ILAENV( 1, 'DORGBR', 'Q', M, M, M, -1 ) )
WRKBL = MAX( WRKBL, BDSPAC )	WRKBL = MAX( WRKBL, BDSPAC )
MAXWRK = M*M + WRKBL	MAXWRK = M*M + WRKBL
MINWRK = MAX( 3*M+N, BDSPAC )	MINWRK = MAX( 3*M + N, BDSPAC )
END IF	END IF
ELSE	ELSE
*	*
* Path 10t(N greater than M, but not much larger)	* Path 10t(N greater than M, but not much larger)
*	*
MAXWRK = 3M + ( M+N )ILAENV( 1, 'DGEBRD', ' ', M, N,	CALL DGEBRD( M, N, A, LDA, S, DUM(1), DUM(1),
$ -1, -1 )	$ DUM(1), DUM(1), -1, IERR )
IF( WNTVS .OR. WNTVO )	LWORK_DGEBRD = INT( DUM(1) )
$ MAXWRK = MAX( MAXWRK, 3M+M	MAXWRK = 3*M + LWORK_DGEBRD
$ ILAENV( 1, 'DORGBR', 'P', M, N, M, -1 ) )	IF( WNTVS .OR. WNTVO ) THEN
IF( WNTVA )	* Compute space needed for DORGBR P
$ MAXWRK = MAX( MAXWRK, 3M+N	CALL DORGBR( 'P', M, N, M, A, N, DUM(1),
$ ILAENV( 1, 'DORGBR', 'P', N, N, M, -1 ) )	$ DUM(1), -1, IERR )
IF( .NOT.WNTUN )	LWORK_DORGBR_P = INT( DUM(1) )
$ MAXWRK = MAX( MAXWRK, 3M+( M-1 )	MAXWRK = MAX( MAXWRK, 3*M + LWORK_DORGBR_P )
$ ILAENV( 1, 'DORGBR', 'Q', M, M, M, -1 ) )	END IF
	IF( WNTVA ) THEN
	CALL DORGBR( 'P', N, N, M, A, N, DUM(1),
	$ DUM(1), -1, IERR )
	LWORK_DORGBR_P = INT( DUM(1) )
	MAXWRK = MAX( MAXWRK, 3*M + LWORK_DORGBR_P )
	END IF
	IF( .NOT.WNTUN ) THEN
	MAXWRK = MAX( MAXWRK, 3*M + LWORK_DORGBR_Q )
	END IF
MAXWRK = MAX( MAXWRK, BDSPAC )	MAXWRK = MAX( MAXWRK, BDSPAC )
MINWRK = MAX( 3*M+N, BDSPAC )	MINWRK = MAX( 3*M + N, BDSPAC )
END IF	END IF
END IF	END IF
MAXWRK = MAX( MAXWRK, MINWRK )	MAXWRK = MAX( MAXWRK, MINWRK )
Line 594	Line 685
IWORK = ITAU + N	IWORK = ITAU + N
*	*
* Compute A=Q*R	* Compute A=Q*R
* (Workspace: need 2N, prefer N+NNB)	* (Workspace: need 2N, prefer N + NNB)
*	*
CALL DGEQRF( M, N, A, LDA, WORK( ITAU ), WORK( IWORK ),	CALL DGEQRF( M, N, A, LDA, WORK( ITAU ), WORK( IWORK ),
$ LWORK-IWORK+1, IERR )	$ LWORK-IWORK+1, IERR )
*	*
* Zero out below R	* Zero out below R
*	*
CALL DLASET( 'L', N-1, N-1, ZERO, ZERO, A( 2, 1 ), LDA )	IF( N .GT. 1 ) THEN
	CALL DLASET( 'L', N-1, N-1, ZERO, ZERO, A( 2, 1 ),
	$ LDA )
	END IF
IE = 1	IE = 1
ITAUQ = IE + N	ITAUQ = IE + N
ITAUP = ITAUQ + N	ITAUP = ITAUQ + N
IWORK = ITAUP + N	IWORK = ITAUP + N
*	*
* Bidiagonalize R in A	* Bidiagonalize R in A
* (Workspace: need 4N, prefer 3N+2NNB)	* (Workspace: need 4N, prefer 3N + 2NNB)
*	*
CALL DGEBRD( N, N, A, LDA, S, WORK( IE ), WORK( ITAUQ ),	CALL DGEBRD( N, N, A, LDA, S, WORK( IE ), WORK( ITAUQ ),
$ WORK( ITAUP ), WORK( IWORK ), LWORK-IWORK+1,	$ WORK( ITAUP ), WORK( IWORK ), LWORK-IWORK+1,
Line 617	Line 711
IF( WNTVO .OR. WNTVAS ) THEN	IF( WNTVO .OR. WNTVAS ) THEN
*	*
* If right singular vectors desired, generate P'.	* If right singular vectors desired, generate P'.
* (Workspace: need 4N-1, prefer 3N+(N-1)*NB)	* (Workspace: need 4N-1, prefer 3N + (N-1)*NB)
*	*
CALL DORGBR( 'P', N, N, N, A, LDA, WORK( ITAUP ),	CALL DORGBR( 'P', N, N, N, A, LDA, WORK( ITAUP ),
$ WORK( IWORK ), LWORK-IWORK+1, IERR )	$ WORK( IWORK ), LWORK-IWORK+1, IERR )
Line 648	Line 742
* Sufficient workspace for a fast algorithm	* Sufficient workspace for a fast algorithm
*	*
IR = 1	IR = 1
IF( LWORK.GE.MAX( WRKBL, LDAN+N )+LDAN ) THEN	IF( LWORK.GE.MAX( WRKBL, LDAN + N ) + LDAN ) THEN
*	*
* WORK(IU) is LDA by N, WORK(IR) is LDA by N	* WORK(IU) is LDA by N, WORK(IR) is LDA by N
*	*
LDWRKU = LDA	LDWRKU = LDA
LDWRKR = LDA	LDWRKR = LDA
ELSE IF( LWORK.GE.MAX( WRKBL, LDAN+N )+NN ) THEN	ELSE IF( LWORK.GE.MAX( WRKBL, LDAN + N ) + NN ) THEN
*	*
* WORK(IU) is LDA by N, WORK(IR) is N by N	* WORK(IU) is LDA by N, WORK(IR) is N by N
*	*
Line 671	Line 765
IWORK = ITAU + N	IWORK = ITAU + N
*	*
* Compute A=Q*R	* Compute A=Q*R
* (Workspace: need NN+2N, prefer NN+N+NNB)	* (Workspace: need NN + 2N, prefer NN + N + NNB)
*	*
CALL DGEQRF( M, N, A, LDA, WORK( ITAU ),	CALL DGEQRF( M, N, A, LDA, WORK( ITAU ),
$ WORK( IWORK ), LWORK-IWORK+1, IERR )	$ WORK( IWORK ), LWORK-IWORK+1, IERR )
Line 683	Line 777
$ LDWRKR )	$ LDWRKR )
*	*
* Generate Q in A	* Generate Q in A
* (Workspace: need NN+2N, prefer NN+N+NNB)	* (Workspace: need NN + 2N, prefer NN + N + NNB)
*	*
CALL DORGQR( M, N, N, A, LDA, WORK( ITAU ),	CALL DORGQR( M, N, N, A, LDA, WORK( ITAU ),
$ WORK( IWORK ), LWORK-IWORK+1, IERR )	$ WORK( IWORK ), LWORK-IWORK+1, IERR )
Line 693	Line 787
IWORK = ITAUP + N	IWORK = ITAUP + N
*	*
* Bidiagonalize R in WORK(IR)	* Bidiagonalize R in WORK(IR)
* (Workspace: need NN+4N, prefer NN+3N+2NNB)	* (Workspace: need NN + 4N, prefer NN + 3N + 2NNB)
*	*
CALL DGEBRD( N, N, WORK( IR ), LDWRKR, S, WORK( IE ),	CALL DGEBRD( N, N, WORK( IR ), LDWRKR, S, WORK( IE ),
$ WORK( ITAUQ ), WORK( ITAUP ),	$ WORK( ITAUQ ), WORK( ITAUP ),
$ WORK( IWORK ), LWORK-IWORK+1, IERR )	$ WORK( IWORK ), LWORK-IWORK+1, IERR )
*	*
* Generate left vectors bidiagonalizing R	* Generate left vectors bidiagonalizing R
* (Workspace: need NN+4N, prefer NN+3N+N*NB)	* (Workspace: need NN + 4N, prefer NN + 3N + N*NB)
*	*
CALL DORGBR( 'Q', N, N, N, WORK( IR ), LDWRKR,	CALL DORGBR( 'Q', N, N, N, WORK( IR ), LDWRKR,
$ WORK( ITAUQ ), WORK( IWORK ),	$ WORK( ITAUQ ), WORK( IWORK ),
Line 709	Line 803
*	*
* Perform bidiagonal QR iteration, computing left	* Perform bidiagonal QR iteration, computing left
* singular vectors of R in WORK(IR)	* singular vectors of R in WORK(IR)
* (Workspace: need N*N+BDSPAC)	* (Workspace: need N*N + BDSPAC)
*	*
CALL DBDSQR( 'U', N, 0, N, 0, S, WORK( IE ), DUM, 1,	CALL DBDSQR( 'U', N, 0, N, 0, S, WORK( IE ), DUM, 1,
$ WORK( IR ), LDWRKR, DUM, 1,	$ WORK( IR ), LDWRKR, DUM, 1,
Line 718	Line 812
*	*
* 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 WORK(IU) and copying to A	* WORK(IR), storing result in WORK(IU) and copying to A
* (Workspace: need NN+2N, prefer NN+MN+N)	* (Workspace: need NN + 2N, prefer NN + MN + N)
*	*
DO 10 I = 1, M, LDWRKU	DO 10 I = 1, M, LDWRKU
CHUNK = MIN( M-I+1, LDWRKU )	CHUNK = MIN( M-I+1, LDWRKU )
Line 739	Line 833
IWORK = ITAUP + N	IWORK = ITAUP + N
*	*
* Bidiagonalize A	* Bidiagonalize A
* (Workspace: need 3N+M, prefer 3N+(M+N)*NB)	* (Workspace: need 3N + M, prefer 3N + (M + N)*NB)
*	*
CALL DGEBRD( M, N, A, LDA, S, WORK( IE ),	CALL DGEBRD( M, N, A, LDA, S, WORK( IE ),
$ WORK( ITAUQ ), WORK( ITAUP ),	$ WORK( ITAUQ ), WORK( ITAUP ),
$ WORK( IWORK ), LWORK-IWORK+1, IERR )	$ WORK( IWORK ), LWORK-IWORK+1, IERR )
*	*
* Generate left vectors bidiagonalizing A	* Generate left vectors bidiagonalizing A
* (Workspace: need 4N, prefer 3N+N*NB)	* (Workspace: need 4N, prefer 3N + N*NB)
*	*
CALL DORGBR( 'Q', M, N, N, A, LDA, WORK( ITAUQ ),	CALL DORGBR( 'Q', M, N, N, A, LDA, WORK( ITAUQ ),
$ WORK( IWORK ), LWORK-IWORK+1, IERR )	$ WORK( IWORK ), LWORK-IWORK+1, IERR )
Line 772	Line 866
* Sufficient workspace for a fast algorithm	* Sufficient workspace for a fast algorithm
*	*
IR = 1	IR = 1
IF( LWORK.GE.MAX( WRKBL, LDAN+N )+LDAN ) THEN	IF( LWORK.GE.MAX( WRKBL, LDAN + N ) + LDAN ) THEN
*	*
* WORK(IU) is LDA by N and WORK(IR) is LDA by N	* WORK(IU) is LDA by N and WORK(IR) is LDA by N
*	*
LDWRKU = LDA	LDWRKU = LDA
LDWRKR = LDA	LDWRKR = LDA
ELSE IF( LWORK.GE.MAX( WRKBL, LDAN+N )+NN ) THEN	ELSE IF( LWORK.GE.MAX( WRKBL, LDAN + N ) + NN ) THEN
*	*
* WORK(IU) is LDA by N and WORK(IR) is N by N	* WORK(IU) is LDA by N and WORK(IR) is N by N
*	*
Line 795	Line 889
IWORK = ITAU + N	IWORK = ITAU + N
*	*
* Compute A=Q*R	* Compute A=Q*R
* (Workspace: need NN+2N, prefer NN+N+NNB)	* (Workspace: need NN + 2N, prefer NN + N + NNB)
*	*
CALL DGEQRF( M, N, A, LDA, WORK( ITAU ),	CALL DGEQRF( M, N, A, LDA, WORK( ITAU ),
$ WORK( IWORK ), LWORK-IWORK+1, IERR )	$ WORK( IWORK ), LWORK-IWORK+1, IERR )
Line 808	Line 902
$ VT( 2, 1 ), LDVT )	$ VT( 2, 1 ), LDVT )
*	*
* Generate Q in A	* Generate Q in A
* (Workspace: need NN+2N, prefer NN+N+NNB)	* (Workspace: need NN + 2N, prefer NN + N + NNB)
*	*
CALL DORGQR( M, N, N, A, LDA, WORK( ITAU ),	CALL DORGQR( M, N, N, A, LDA, WORK( ITAU ),
$ WORK( IWORK ), LWORK-IWORK+1, IERR )	$ WORK( IWORK ), LWORK-IWORK+1, IERR )
Line 818	Line 912
IWORK = ITAUP + N	IWORK = ITAUP + N
*	*
* Bidiagonalize R in VT, copying result to WORK(IR)	* Bidiagonalize R in VT, copying result to WORK(IR)
* (Workspace: need NN+4N, prefer NN+3N+2NNB)	* (Workspace: need NN + 4N, prefer NN + 3N + 2NNB)
*	*
CALL DGEBRD( N, N, VT, LDVT, S, WORK( IE ),	CALL DGEBRD( N, N, VT, LDVT, S, WORK( IE ),
$ WORK( ITAUQ ), WORK( ITAUP ),	$ WORK( ITAUQ ), WORK( ITAUP ),
Line 826	Line 920
CALL DLACPY( 'L', N, N, VT, LDVT, WORK( IR ), LDWRKR )	CALL DLACPY( 'L', N, N, VT, LDVT, WORK( IR ), LDWRKR )
*	*
* Generate left vectors bidiagonalizing R in WORK(IR)	* Generate left vectors bidiagonalizing R in WORK(IR)
* (Workspace: need NN+4N, prefer NN+3N+N*NB)	* (Workspace: need NN + 4N, prefer NN + 3N + N*NB)
*	*
CALL DORGBR( 'Q', N, N, N, WORK( IR ), LDWRKR,	CALL DORGBR( 'Q', N, N, N, WORK( IR ), LDWRKR,
$ WORK( ITAUQ ), WORK( IWORK ),	$ WORK( ITAUQ ), WORK( IWORK ),
$ LWORK-IWORK+1, IERR )	$ LWORK-IWORK+1, IERR )
*	*
* Generate right vectors bidiagonalizing R in VT	* Generate right vectors bidiagonalizing R in VT
* (Workspace: need NN+4N-1, prefer NN+3N+(N-1)*NB)	* (Workspace: need NN + 4N-1, prefer NN + 3N + (N-1)*NB)
*	*
CALL DORGBR( 'P', N, N, N, VT, LDVT, WORK( ITAUP ),	CALL DORGBR( 'P', N, N, N, VT, LDVT, WORK( ITAUP ),
$ WORK( IWORK ), LWORK-IWORK+1, IERR )	$ WORK( IWORK ), LWORK-IWORK+1, IERR )
Line 842	Line 936
* Perform bidiagonal QR iteration, computing left	* Perform bidiagonal QR iteration, computing left
* singular vectors of R in WORK(IR) and computing right	* singular vectors of R in WORK(IR) and computing right
* singular vectors of R in VT	* singular vectors of R in VT
* (Workspace: need N*N+BDSPAC)	* (Workspace: need N*N + BDSPAC)
*	*
CALL DBDSQR( 'U', N, N, N, 0, S, WORK( IE ), VT, LDVT,	CALL DBDSQR( 'U', N, N, N, 0, S, WORK( IE ), VT, LDVT,
$ WORK( IR ), LDWRKR, DUM, 1,	$ WORK( IR ), LDWRKR, DUM, 1,
Line 851	Line 945
*	*
* 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 WORK(IU) and copying to A	* WORK(IR), storing result in WORK(IU) and copying to A
* (Workspace: need NN+2N, prefer NN+MN+N)	* (Workspace: need NN + 2N, prefer NN + MN + N)
*	*
DO 20 I = 1, M, LDWRKU	DO 20 I = 1, M, LDWRKU
CHUNK = MIN( M-I+1, LDWRKU )	CHUNK = MIN( M-I+1, LDWRKU )
Line 870	Line 964
IWORK = ITAU + N	IWORK = ITAU + N
*	*
* Compute A=Q*R	* Compute A=Q*R
* (Workspace: need 2N, prefer N+NNB)	* (Workspace: need 2N, prefer N + NNB)
*	*
CALL DGEQRF( M, N, A, LDA, WORK( ITAU ),	CALL DGEQRF( M, N, A, LDA, WORK( ITAU ),
$ WORK( IWORK ), LWORK-IWORK+1, IERR )	$ WORK( IWORK ), LWORK-IWORK+1, IERR )
Line 883	Line 977
$ VT( 2, 1 ), LDVT )	$ VT( 2, 1 ), LDVT )
*	*
* Generate Q in A	* Generate Q in A
* (Workspace: need 2N, prefer N+NNB)	* (Workspace: need 2N, prefer N + NNB)
*	*
CALL DORGQR( M, N, N, A, LDA, WORK( ITAU ),	CALL DORGQR( M, N, N, A, LDA, WORK( ITAU ),
$ WORK( IWORK ), LWORK-IWORK+1, IERR )	$ WORK( IWORK ), LWORK-IWORK+1, IERR )
Line 893	Line 987
IWORK = ITAUP + N	IWORK = ITAUP + N
*	*
* Bidiagonalize R in VT	* Bidiagonalize R in VT
* (Workspace: need 4N, prefer 3N+2NNB)	* (Workspace: need 4N, prefer 3N + 2NNB)
*	*
CALL DGEBRD( N, N, VT, LDVT, S, WORK( IE ),	CALL DGEBRD( N, N, VT, LDVT, S, WORK( IE ),
$ WORK( ITAUQ ), WORK( ITAUP ),	$ WORK( ITAUQ ), WORK( ITAUP ),
$ WORK( IWORK ), LWORK-IWORK+1, IERR )	$ WORK( IWORK ), LWORK-IWORK+1, IERR )
*	*
* Multiply Q in A by left vectors bidiagonalizing R	* Multiply Q in A by left vectors bidiagonalizing R
* (Workspace: need 3N+M, prefer 3N+M*NB)	* (Workspace: need 3N + M, prefer 3N + M*NB)
*	*
CALL DORMBR( 'Q', 'R', 'N', M, N, N, VT, LDVT,	CALL DORMBR( 'Q', 'R', 'N', M, N, N, VT, LDVT,
$ WORK( ITAUQ ), A, LDA, WORK( IWORK ),	$ WORK( ITAUQ ), A, LDA, WORK( IWORK ),
$ LWORK-IWORK+1, IERR )	$ LWORK-IWORK+1, IERR )
*	*
* Generate right vectors bidiagonalizing R in VT	* Generate right vectors bidiagonalizing R in VT
* (Workspace: need 4N-1, prefer 3N+(N-1)*NB)	* (Workspace: need 4N-1, prefer 3N + (N-1)*NB)
*	*
CALL DORGBR( 'P', N, N, N, VT, LDVT, WORK( ITAUP ),	CALL DORGBR( 'P', N, N, N, VT, LDVT, WORK( ITAUP ),
$ WORK( IWORK ), LWORK-IWORK+1, IERR )	$ WORK( IWORK ), LWORK-IWORK+1, IERR )
Line 951	Line 1045
IWORK = ITAU + N	IWORK = ITAU + N
*	*
* Compute A=Q*R	* Compute A=Q*R
* (Workspace: need NN+2N, prefer NN+N+NNB)	* (Workspace: need NN + 2N, prefer NN + N + NNB)
*	*
CALL DGEQRF( M, N, A, LDA, WORK( ITAU ),	CALL DGEQRF( M, N, A, LDA, WORK( ITAU ),
$ WORK( IWORK ), LWORK-IWORK+1, IERR )	$ WORK( IWORK ), LWORK-IWORK+1, IERR )
Line 964	Line 1058
$ WORK( IR+1 ), LDWRKR )	$ WORK( IR+1 ), LDWRKR )
*	*
* Generate Q in A	* Generate Q in A
* (Workspace: need NN+2N, prefer NN+N+NNB)	* (Workspace: need NN + 2N, prefer NN + N + NNB)
*	*
CALL DORGQR( M, N, N, A, LDA, WORK( ITAU ),	CALL DORGQR( M, N, N, A, LDA, WORK( ITAU ),
$ WORK( IWORK ), LWORK-IWORK+1, IERR )	$ WORK( IWORK ), LWORK-IWORK+1, IERR )
Line 974	Line 1068
IWORK = ITAUP + N	IWORK = ITAUP + N
*	*
* Bidiagonalize R in WORK(IR)	* Bidiagonalize R in WORK(IR)
* (Workspace: need NN+4N, prefer NN+3N+2NNB)	* (Workspace: need NN + 4N, prefer NN + 3N + 2NNB)
*	*
CALL DGEBRD( N, N, WORK( IR ), LDWRKR, S,	CALL DGEBRD( N, N, WORK( IR ), LDWRKR, S,
$ WORK( IE ), WORK( ITAUQ ),	$ WORK( IE ), WORK( ITAUQ ),
Line 982	Line 1076
$ LWORK-IWORK+1, IERR )	$ LWORK-IWORK+1, IERR )
*	*
* Generate left vectors bidiagonalizing R in WORK(IR)	* Generate left vectors bidiagonalizing R in WORK(IR)
* (Workspace: need NN+4N, prefer NN+3N+N*NB)	* (Workspace: need NN + 4N, prefer NN + 3N + N*NB)
*	*
CALL DORGBR( 'Q', N, N, N, WORK( IR ), LDWRKR,	CALL DORGBR( 'Q', N, N, N, WORK( IR ), LDWRKR,
$ WORK( ITAUQ ), WORK( IWORK ),	$ WORK( ITAUQ ), WORK( IWORK ),
Line 991	Line 1085
*	*
* Perform bidiagonal QR iteration, computing left	* Perform bidiagonal QR iteration, computing left
* singular vectors of R in WORK(IR)	* singular vectors of R in WORK(IR)
* (Workspace: need N*N+BDSPAC)	* (Workspace: need N*N + BDSPAC)
*	*
CALL DBDSQR( 'U', N, 0, N, 0, S, WORK( IE ), DUM,	CALL DBDSQR( 'U', N, 0, N, 0, S, WORK( IE ), DUM,
$ 1, WORK( IR ), LDWRKR, DUM, 1,	$ 1, WORK( IR ), LDWRKR, DUM, 1,
Line 1012	Line 1106
IWORK = ITAU + N	IWORK = ITAU + N
*	*
* Compute A=Q*R, copying result to U	* Compute A=Q*R, copying result to U
* (Workspace: need 2N, prefer N+NNB)	* (Workspace: need 2N, prefer N + NNB)
*	*
CALL DGEQRF( M, N, A, LDA, WORK( ITAU ),	CALL DGEQRF( M, N, A, LDA, WORK( ITAU ),
$ WORK( IWORK ), LWORK-IWORK+1, IERR )	$ WORK( IWORK ), LWORK-IWORK+1, IERR )
CALL DLACPY( 'L', M, N, A, LDA, U, LDU )	CALL DLACPY( 'L', M, N, A, LDA, U, LDU )
*	*
* Generate Q in U	* Generate Q in U
* (Workspace: need 2N, prefer N+NNB)	* (Workspace: need 2N, prefer N + NNB)
*	*
CALL DORGQR( M, N, N, U, LDU, WORK( ITAU ),	CALL DORGQR( M, N, N, U, LDU, WORK( ITAU ),
$ WORK( IWORK ), LWORK-IWORK+1, IERR )	$ WORK( IWORK ), LWORK-IWORK+1, IERR )
Line 1030	Line 1124
*	*
* Zero out below R in A	* Zero out below R in A
*	*
CALL DLASET( 'L', N-1, N-1, ZERO, ZERO, A( 2, 1 ),	IF( N .GT. 1 ) THEN
$ LDA )	CALL DLASET( 'L', N-1, N-1, ZERO, ZERO,
	$ A( 2, 1 ), LDA )
	END IF
*	*
* Bidiagonalize R in A	* Bidiagonalize R in A
* (Workspace: need 4N, prefer 3N+2NNB)	* (Workspace: need 4N, prefer 3N + 2NNB)
*	*
CALL DGEBRD( N, N, A, LDA, S, WORK( IE ),	CALL DGEBRD( N, N, A, LDA, S, WORK( IE ),
$ WORK( ITAUQ ), WORK( ITAUP ),	$ WORK( ITAUQ ), WORK( ITAUP ),
$ WORK( IWORK ), LWORK-IWORK+1, IERR )	$ WORK( IWORK ), LWORK-IWORK+1, IERR )
*	*
* Multiply Q in U by left vectors bidiagonalizing R	* Multiply Q in U by left vectors bidiagonalizing R
* (Workspace: need 3N+M, prefer 3N+M*NB)	* (Workspace: need 3N + M, prefer 3N + M*NB)
*	*
CALL DORMBR( 'Q', 'R', 'N', M, N, N, A, LDA,	CALL DORMBR( 'Q', 'R', 'N', M, N, N, A, LDA,
$ WORK( ITAUQ ), U, LDU, WORK( IWORK ),	$ WORK( ITAUQ ), U, LDU, WORK( IWORK ),
Line 1076	Line 1172
LDWRKU = LDA	LDWRKU = LDA
IR = IU + LDWRKU*N	IR = IU + LDWRKU*N
LDWRKR = LDA	LDWRKR = LDA
ELSE IF( LWORK.GE.WRKBL+( LDA+N )*N ) THEN	ELSE IF( LWORK.GE.WRKBL+( LDA + N )*N ) THEN
*	*
* WORK(IU) is LDA by N and WORK(IR) is N by N	* WORK(IU) is LDA by N and WORK(IR) is N by N
*	*
Line 1095	Line 1191
IWORK = ITAU + N	IWORK = ITAU + N
*	*
* Compute A=Q*R	* Compute A=Q*R
* (Workspace: need 2NN+2N, prefer 2NN+N+NNB)	* (Workspace: need 2NN + 2N, prefer 2NN + N + NNB)
*	*
CALL DGEQRF( M, N, A, LDA, WORK( ITAU ),	CALL DGEQRF( M, N, A, LDA, WORK( ITAU ),
$ WORK( IWORK ), LWORK-IWORK+1, IERR )	$ WORK( IWORK ), LWORK-IWORK+1, IERR )
Line 1108	Line 1204
$ WORK( IU+1 ), LDWRKU )	$ WORK( IU+1 ), LDWRKU )
*	*
* Generate Q in A	* Generate Q in A
* (Workspace: need 2NN+2N, prefer 2NN+N+NNB)	* (Workspace: need 2NN + 2N, prefer 2NN + N + NNB)
*	*
CALL DORGQR( M, N, N, A, LDA, WORK( ITAU ),	CALL DORGQR( M, N, N, A, LDA, WORK( ITAU ),
$ WORK( IWORK ), LWORK-IWORK+1, IERR )	$ WORK( IWORK ), LWORK-IWORK+1, IERR )
Line 1119	Line 1215
*	*
* Bidiagonalize R in WORK(IU), copying result to	* Bidiagonalize R in WORK(IU), copying result to
* WORK(IR)	* WORK(IR)
* (Workspace: need 2NN+4*N,	* (Workspace: need 2NN + 4*N,
* prefer 2NN+3N+2N*NB)	* prefer 2NN+3N+2N*NB)
*	*
CALL DGEBRD( N, N, WORK( IU ), LDWRKU, S,	CALL DGEBRD( N, N, WORK( IU ), LDWRKU, S,
Line 1130	Line 1226
$ WORK( IR ), LDWRKR )	$ WORK( IR ), LDWRKR )
*	*
* Generate left bidiagonalizing vectors in WORK(IU)	* Generate left bidiagonalizing vectors in WORK(IU)
* (Workspace: need 2NN+4N, prefer 2NN+3N+N*NB)	* (Workspace: need 2NN + 4N, prefer 2NN + 3N + N*NB)
*	*
CALL DORGBR( 'Q', N, N, N, WORK( IU ), LDWRKU,	CALL DORGBR( 'Q', N, N, N, WORK( IU ), LDWRKU,
$ WORK( ITAUQ ), WORK( IWORK ),	$ WORK( ITAUQ ), WORK( IWORK ),
$ LWORK-IWORK+1, IERR )	$ LWORK-IWORK+1, IERR )
*	*
* Generate right bidiagonalizing vectors in WORK(IR)	* Generate right bidiagonalizing vectors in WORK(IR)
* (Workspace: need 2NN+4*N-1,	* (Workspace: need 2NN + 4*N-1,
* prefer 2NN+3N+(N-1)NB)	* prefer 2NN+3N+(N-1)NB)
*	*
CALL DORGBR( 'P', N, N, N, WORK( IR ), LDWRKR,	CALL DORGBR( 'P', N, N, N, WORK( IR ), LDWRKR,
Line 1148	Line 1244
* Perform bidiagonal QR iteration, computing left	* Perform bidiagonal QR iteration, computing left
* singular vectors of R in WORK(IU) and computing	* singular vectors of R in WORK(IU) and computing
* right singular vectors of R in WORK(IR)	* right singular vectors of R in WORK(IR)
* (Workspace: need 2NN+BDSPAC)	* (Workspace: need 2NN + BDSPAC)
*	*
CALL DBDSQR( 'U', N, N, N, 0, S, WORK( IE ),	CALL DBDSQR( 'U', N, N, N, 0, S, WORK( IE ),
$ WORK( IR ), LDWRKR, WORK( IU ),	$ WORK( IR ), LDWRKR, WORK( IU ),
Line 1175	Line 1271
IWORK = ITAU + N	IWORK = ITAU + N
*	*
* Compute A=Q*R, copying result to U	* Compute A=Q*R, copying result to U
* (Workspace: need 2N, prefer N+NNB)	* (Workspace: need 2N, prefer N + NNB)
*	*
CALL DGEQRF( M, N, A, LDA, WORK( ITAU ),	CALL DGEQRF( M, N, A, LDA, WORK( ITAU ),
$ WORK( IWORK ), LWORK-IWORK+1, IERR )	$ WORK( IWORK ), LWORK-IWORK+1, IERR )
CALL DLACPY( 'L', M, N, A, LDA, U, LDU )	CALL DLACPY( 'L', M, N, A, LDA, U, LDU )
*	*
* Generate Q in U	* Generate Q in U
* (Workspace: need 2N, prefer N+NNB)	* (Workspace: need 2N, prefer N + NNB)
*	*
CALL DORGQR( M, N, N, U, LDU, WORK( ITAU ),	CALL DORGQR( M, N, N, U, LDU, WORK( ITAU ),
$ WORK( IWORK ), LWORK-IWORK+1, IERR )	$ WORK( IWORK ), LWORK-IWORK+1, IERR )
Line 1193	Line 1289
*	*
* Zero out below R in A	* Zero out below R in A
*	*
CALL DLASET( 'L', N-1, N-1, ZERO, ZERO, A( 2, 1 ),	IF( N .GT. 1 ) THEN
$ LDA )	CALL DLASET( 'L', N-1, N-1, ZERO, ZERO,
	$ A( 2, 1 ), LDA )
	END IF
*	*
* Bidiagonalize R in A	* Bidiagonalize R in A
* (Workspace: need 4N, prefer 3N+2NNB)	* (Workspace: need 4N, prefer 3N + 2NNB)
*	*
CALL DGEBRD( N, N, A, LDA, S, WORK( IE ),	CALL DGEBRD( N, N, A, LDA, S, WORK( IE ),
$ WORK( ITAUQ ), WORK( ITAUP ),	$ WORK( ITAUQ ), WORK( ITAUP ),
$ WORK( IWORK ), LWORK-IWORK+1, IERR )	$ WORK( IWORK ), LWORK-IWORK+1, IERR )
*	*
* Multiply Q in U by left vectors bidiagonalizing R	* Multiply Q in U by left vectors bidiagonalizing R
* (Workspace: need 3N+M, prefer 3N+M*NB)	* (Workspace: need 3N + M, prefer 3N + M*NB)
*	*
CALL DORMBR( 'Q', 'R', 'N', M, N, N, A, LDA,	CALL DORMBR( 'Q', 'R', 'N', M, N, N, A, LDA,
$ WORK( ITAUQ ), U, LDU, WORK( IWORK ),	$ WORK( ITAUQ ), U, LDU, WORK( IWORK ),
$ LWORK-IWORK+1, IERR )	$ LWORK-IWORK+1, IERR )
*	*
* Generate right vectors bidiagonalizing R in A	* Generate right vectors bidiagonalizing R in A
* (Workspace: need 4N-1, prefer 3N+(N-1)*NB)	* (Workspace: need 4N-1, prefer 3N + (N-1)*NB)
*	*
CALL DORGBR( 'P', N, N, N, A, LDA, WORK( ITAUP ),	CALL DORGBR( 'P', N, N, N, A, LDA, WORK( ITAUP ),
$ WORK( IWORK ), LWORK-IWORK+1, IERR )	$ WORK( IWORK ), LWORK-IWORK+1, IERR )
Line 1255	Line 1353
IWORK = ITAU + N	IWORK = ITAU + N
*	*
* Compute A=Q*R	* Compute A=Q*R
* (Workspace: need NN+2N, prefer NN+N+NNB)	* (Workspace: need NN + 2N, prefer NN + N + NNB)
*	*
CALL DGEQRF( M, N, A, LDA, WORK( ITAU ),	CALL DGEQRF( M, N, A, LDA, WORK( ITAU ),
$ WORK( IWORK ), LWORK-IWORK+1, IERR )	$ WORK( IWORK ), LWORK-IWORK+1, IERR )
Line 1268	Line 1366
$ WORK( IU+1 ), LDWRKU )	$ WORK( IU+1 ), LDWRKU )
*	*
* Generate Q in A	* Generate Q in A
* (Workspace: need NN+2N, prefer NN+N+NNB)	* (Workspace: need NN + 2N, prefer NN + N + NNB)
*	*
CALL DORGQR( M, N, N, A, LDA, WORK( ITAU ),	CALL DORGQR( M, N, N, A, LDA, WORK( ITAU ),
$ WORK( IWORK ), LWORK-IWORK+1, IERR )	$ WORK( IWORK ), LWORK-IWORK+1, IERR )
Line 1278	Line 1376
IWORK = ITAUP + N	IWORK = ITAUP + N
*	*
* Bidiagonalize R in WORK(IU), copying result to VT	* Bidiagonalize R in WORK(IU), copying result to VT
* (Workspace: need NN+4N, prefer NN+3N+2NNB)	* (Workspace: need NN + 4N, prefer NN + 3N + 2NNB)
*	*
CALL DGEBRD( N, N, WORK( IU ), LDWRKU, S,	CALL DGEBRD( N, N, WORK( IU ), LDWRKU, S,
$ WORK( IE ), WORK( ITAUQ ),	$ WORK( IE ), WORK( ITAUQ ),
Line 1288	Line 1386
$ LDVT )	$ LDVT )
*	*
* Generate left bidiagonalizing vectors in WORK(IU)	* Generate left bidiagonalizing vectors in WORK(IU)
* (Workspace: need NN+4N, prefer NN+3N+N*NB)	* (Workspace: need NN + 4N, prefer NN + 3N + N*NB)
*	*
CALL DORGBR( 'Q', N, N, N, WORK( IU ), LDWRKU,	CALL DORGBR( 'Q', N, N, N, WORK( IU ), LDWRKU,
$ WORK( ITAUQ ), WORK( IWORK ),	$ WORK( ITAUQ ), WORK( IWORK ),
$ LWORK-IWORK+1, IERR )	$ LWORK-IWORK+1, IERR )
*	*
* Generate right bidiagonalizing vectors in VT	* Generate right bidiagonalizing vectors in VT
* (Workspace: need NN+4N-1,	* (Workspace: need NN + 4N-1,
* prefer NN+3N+(N-1)*NB)	* prefer NN+3N+(N-1)*NB)
*	*
CALL DORGBR( 'P', N, N, N, VT, LDVT, WORK( ITAUP ),	CALL DORGBR( 'P', N, N, N, VT, LDVT, WORK( ITAUP ),
Line 1305	Line 1403
* Perform bidiagonal QR iteration, computing left	* Perform bidiagonal QR iteration, computing left
* singular vectors of R in WORK(IU) and computing	* singular vectors of R in WORK(IU) and computing
* right singular vectors of R in VT	* right singular vectors of R in VT
* (Workspace: need N*N+BDSPAC)	* (Workspace: need N*N + BDSPAC)
*	*
CALL DBDSQR( 'U', N, N, N, 0, S, WORK( IE ), VT,	CALL DBDSQR( 'U', N, N, N, 0, S, WORK( IE ), VT,
$ LDVT, WORK( IU ), LDWRKU, DUM, 1,	$ LDVT, WORK( IU ), LDWRKU, DUM, 1,
Line 1326	Line 1424
IWORK = ITAU + N	IWORK = ITAU + N
*	*
* Compute A=Q*R, copying result to U	* Compute A=Q*R, copying result to U
* (Workspace: need 2N, prefer N+NNB)	* (Workspace: need 2N, prefer N + NNB)
*	*
CALL DGEQRF( M, N, A, LDA, WORK( ITAU ),	CALL DGEQRF( M, N, A, LDA, WORK( ITAU ),
$ WORK( IWORK ), LWORK-IWORK+1, IERR )	$ WORK( IWORK ), LWORK-IWORK+1, IERR )
CALL DLACPY( 'L', M, N, A, LDA, U, LDU )	CALL DLACPY( 'L', M, N, A, LDA, U, LDU )
*	*
* Generate Q in U	* Generate Q in U
* (Workspace: need 2N, prefer N+NNB)	* (Workspace: need 2N, prefer N + NNB)
*	*
CALL DORGQR( M, N, N, U, LDU, WORK( ITAU ),	CALL DORGQR( M, N, N, U, LDU, WORK( ITAU ),
$ WORK( IWORK ), LWORK-IWORK+1, IERR )	$ WORK( IWORK ), LWORK-IWORK+1, IERR )
Line 1350	Line 1448
IWORK = ITAUP + N	IWORK = ITAUP + N
*	*
* Bidiagonalize R in VT	* Bidiagonalize R in VT
* (Workspace: need 4N, prefer 3N+2NNB)	* (Workspace: need 4N, prefer 3N + 2NNB)
*	*
CALL DGEBRD( N, N, VT, LDVT, S, WORK( IE ),	CALL DGEBRD( N, N, VT, LDVT, S, WORK( IE ),
$ WORK( ITAUQ ), WORK( ITAUP ),	$ WORK( ITAUQ ), WORK( ITAUP ),
Line 1358	Line 1456
*	*
* Multiply Q in U by left bidiagonalizing vectors	* Multiply Q in U by left bidiagonalizing vectors
* in VT	* in VT
* (Workspace: need 3N+M, prefer 3N+M*NB)	* (Workspace: need 3N + M, prefer 3N + M*NB)
*	*
CALL DORMBR( 'Q', 'R', 'N', M, N, N, VT, LDVT,	CALL DORMBR( 'Q', 'R', 'N', M, N, N, VT, LDVT,
$ WORK( ITAUQ ), U, LDU, WORK( IWORK ),	$ WORK( ITAUQ ), U, LDU, WORK( IWORK ),
$ LWORK-IWORK+1, IERR )	$ LWORK-IWORK+1, IERR )
*	*
* Generate right bidiagonalizing vectors in VT	* Generate right bidiagonalizing vectors in VT
* (Workspace: need 4N-1, prefer 3N+(N-1)*NB)	* (Workspace: need 4N-1, prefer 3N + (N-1)*NB)
*	*
CALL DORGBR( 'P', N, N, N, VT, LDVT, WORK( ITAUP ),	CALL DORGBR( 'P', N, N, N, VT, LDVT, WORK( ITAUP ),
$ WORK( IWORK ), LWORK-IWORK+1, IERR )	$ WORK( IWORK ), LWORK-IWORK+1, IERR )
Line 1412	Line 1510
IWORK = ITAU + N	IWORK = ITAU + N
*	*
* Compute A=Q*R, copying result to U	* Compute A=Q*R, copying result to U
* (Workspace: need NN+2N, prefer NN+N+NNB)	* (Workspace: need NN + 2N, prefer NN + N + NNB)
*	*
CALL DGEQRF( M, N, A, LDA, WORK( ITAU ),	CALL DGEQRF( M, N, A, LDA, WORK( ITAU ),
$ WORK( IWORK ), LWORK-IWORK+1, IERR )	$ WORK( IWORK ), LWORK-IWORK+1, IERR )
Line 1426	Line 1524
$ WORK( IR+1 ), LDWRKR )	$ WORK( IR+1 ), LDWRKR )
*	*
* Generate Q in U	* Generate Q in U
* (Workspace: need NN+N+M, prefer NN+N+M*NB)	* (Workspace: need NN + N + M, prefer NN + N + M*NB)
*	*
CALL DORGQR( M, M, N, U, LDU, WORK( ITAU ),	CALL DORGQR( M, M, N, U, LDU, WORK( ITAU ),
$ WORK( IWORK ), LWORK-IWORK+1, IERR )	$ WORK( IWORK ), LWORK-IWORK+1, IERR )
Line 1436	Line 1534
IWORK = ITAUP + N	IWORK = ITAUP + N
*	*
* Bidiagonalize R in WORK(IR)	* Bidiagonalize R in WORK(IR)
* (Workspace: need NN+4N, prefer NN+3N+2NNB)	* (Workspace: need NN + 4N, prefer NN + 3N + 2NNB)
*	*
CALL DGEBRD( N, N, WORK( IR ), LDWRKR, S,	CALL DGEBRD( N, N, WORK( IR ), LDWRKR, S,
$ WORK( IE ), WORK( ITAUQ ),	$ WORK( IE ), WORK( ITAUQ ),
Line 1444	Line 1542
$ LWORK-IWORK+1, IERR )	$ LWORK-IWORK+1, IERR )
*	*
* Generate left bidiagonalizing vectors in WORK(IR)	* Generate left bidiagonalizing vectors in WORK(IR)
* (Workspace: need NN+4N, prefer NN+3N+N*NB)	* (Workspace: need NN + 4N, prefer NN + 3N + N*NB)
*	*
CALL DORGBR( 'Q', N, N, N, WORK( IR ), LDWRKR,	CALL DORGBR( 'Q', N, N, N, WORK( IR ), LDWRKR,
$ WORK( ITAUQ ), WORK( IWORK ),	$ WORK( ITAUQ ), WORK( IWORK ),
Line 1453	Line 1551
*	*
* Perform bidiagonal QR iteration, computing left	* Perform bidiagonal QR iteration, computing left
* singular vectors of R in WORK(IR)	* singular vectors of R in WORK(IR)
* (Workspace: need N*N+BDSPAC)	* (Workspace: need N*N + BDSPAC)
*	*
CALL DBDSQR( 'U', N, 0, N, 0, S, WORK( IE ), DUM,	CALL DBDSQR( 'U', N, 0, N, 0, S, WORK( IE ), DUM,
$ 1, WORK( IR ), LDWRKR, DUM, 1,	$ 1, WORK( IR ), LDWRKR, DUM, 1,
Line 1478	Line 1576
IWORK = ITAU + N	IWORK = ITAU + N
*	*
* Compute A=Q*R, copying result to U	* Compute A=Q*R, copying result to U
* (Workspace: need 2N, prefer N+NNB)	* (Workspace: need 2N, prefer N + NNB)
*	*
CALL DGEQRF( M, N, A, LDA, WORK( ITAU ),	CALL DGEQRF( M, N, A, LDA, WORK( ITAU ),
$ WORK( IWORK ), LWORK-IWORK+1, IERR )	$ WORK( IWORK ), LWORK-IWORK+1, IERR )
CALL DLACPY( 'L', M, N, A, LDA, U, LDU )	CALL DLACPY( 'L', M, N, A, LDA, U, LDU )
*	*
* Generate Q in U	* Generate Q in U
* (Workspace: need N+M, prefer N+M*NB)	* (Workspace: need N + M, prefer N + M*NB)
*	*
CALL DORGQR( M, M, N, U, LDU, WORK( ITAU ),	CALL DORGQR( M, M, N, U, LDU, WORK( ITAU ),
$ WORK( IWORK ), LWORK-IWORK+1, IERR )	$ WORK( IWORK ), LWORK-IWORK+1, IERR )
Line 1496	Line 1594
*	*
* Zero out below R in A	* Zero out below R in A
*	*
CALL DLASET( 'L', N-1, N-1, ZERO, ZERO, A( 2, 1 ),	IF( N .GT. 1 ) THEN
$ LDA )	CALL DLASET( 'L', N-1, N-1, ZERO, ZERO,
	$ A( 2, 1 ), LDA )
	END IF
*	*
* Bidiagonalize R in A	* Bidiagonalize R in A
* (Workspace: need 4N, prefer 3N+2NNB)	* (Workspace: need 4N, prefer 3N + 2NNB)
*	*
CALL DGEBRD( N, N, A, LDA, S, WORK( IE ),	CALL DGEBRD( N, N, A, LDA, S, WORK( IE ),
$ WORK( ITAUQ ), WORK( ITAUP ),	$ WORK( ITAUQ ), WORK( ITAUP ),
Line 1508	Line 1608
*	*
* Multiply Q in U by left bidiagonalizing vectors	* Multiply Q in U by left bidiagonalizing vectors
* in A	* in A
* (Workspace: need 3N+M, prefer 3N+M*NB)	* (Workspace: need 3N + M, prefer 3N + M*NB)
*	*
CALL DORMBR( 'Q', 'R', 'N', M, N, N, A, LDA,	CALL DORMBR( 'Q', 'R', 'N', M, N, N, A, LDA,
$ WORK( ITAUQ ), U, LDU, WORK( IWORK ),	$ WORK( ITAUQ ), U, LDU, WORK( IWORK ),
Line 1543	Line 1643
LDWRKU = LDA	LDWRKU = LDA
IR = IU + LDWRKU*N	IR = IU + LDWRKU*N
LDWRKR = LDA	LDWRKR = LDA
ELSE IF( LWORK.GE.WRKBL+( LDA+N )*N ) THEN	ELSE IF( LWORK.GE.WRKBL+( LDA + N )*N ) THEN
*	*
* WORK(IU) is LDA by N and WORK(IR) is N by N	* WORK(IU) is LDA by N and WORK(IR) is N by N
*	*
Line 1562	Line 1662
IWORK = ITAU + N	IWORK = ITAU + N
*	*
* Compute A=Q*R, copying result to U	* Compute A=Q*R, copying result to U
* (Workspace: need 2NN+2N, prefer 2NN+N+NNB)	* (Workspace: need 2NN + 2N, prefer 2NN + N + NNB)
*	*
CALL DGEQRF( M, N, A, LDA, WORK( ITAU ),	CALL DGEQRF( M, N, A, LDA, WORK( ITAU ),
$ WORK( IWORK ), LWORK-IWORK+1, IERR )	$ WORK( IWORK ), LWORK-IWORK+1, IERR )
CALL DLACPY( 'L', M, N, A, LDA, U, LDU )	CALL DLACPY( 'L', M, N, A, LDA, U, LDU )
*	*
* Generate Q in U	* Generate Q in U
* (Workspace: need 2NN+N+M, prefer 2NN+N+M*NB)	* (Workspace: need 2NN + N + M, prefer 2NN + N + M*NB)
*	*
CALL DORGQR( M, M, N, U, LDU, WORK( ITAU ),	CALL DORGQR( M, M, N, U, LDU, WORK( ITAU ),
$ WORK( IWORK ), LWORK-IWORK+1, IERR )	$ WORK( IWORK ), LWORK-IWORK+1, IERR )
Line 1587	Line 1687
*	*
* Bidiagonalize R in WORK(IU), copying result to	* Bidiagonalize R in WORK(IU), copying result to
* WORK(IR)	* WORK(IR)
* (Workspace: need 2NN+4*N,	* (Workspace: need 2NN + 4*N,
* prefer 2NN+3N+2N*NB)	* prefer 2NN+3N+2N*NB)
*	*
CALL DGEBRD( N, N, WORK( IU ), LDWRKU, S,	CALL DGEBRD( N, N, WORK( IU ), LDWRKU, S,
Line 1598	Line 1698
$ WORK( IR ), LDWRKR )	$ WORK( IR ), LDWRKR )
*	*
* Generate left bidiagonalizing vectors in WORK(IU)	* Generate left bidiagonalizing vectors in WORK(IU)
* (Workspace: need 2NN+4N, prefer 2NN+3N+N*NB)	* (Workspace: need 2NN + 4N, prefer 2NN + 3N + N*NB)
*	*
CALL DORGBR( 'Q', N, N, N, WORK( IU ), LDWRKU,	CALL DORGBR( 'Q', N, N, N, WORK( IU ), LDWRKU,
$ WORK( ITAUQ ), WORK( IWORK ),	$ WORK( ITAUQ ), WORK( IWORK ),
$ LWORK-IWORK+1, IERR )	$ LWORK-IWORK+1, IERR )
*	*
* Generate right bidiagonalizing vectors in WORK(IR)	* Generate right bidiagonalizing vectors in WORK(IR)
* (Workspace: need 2NN+4*N-1,	* (Workspace: need 2NN + 4*N-1,
* prefer 2NN+3N+(N-1)NB)	* prefer 2NN+3N+(N-1)NB)
*	*
CALL DORGBR( 'P', N, N, N, WORK( IR ), LDWRKR,	CALL DORGBR( 'P', N, N, N, WORK( IR ), LDWRKR,
Line 1616	Line 1716
* Perform bidiagonal QR iteration, computing left	* Perform bidiagonal QR iteration, computing left
* singular vectors of R in WORK(IU) and computing	* singular vectors of R in WORK(IU) and computing
* right singular vectors of R in WORK(IR)	* right singular vectors of R in WORK(IR)
* (Workspace: need 2NN+BDSPAC)	* (Workspace: need 2NN + BDSPAC)
*	*
CALL DBDSQR( 'U', N, N, N, 0, S, WORK( IE ),	CALL DBDSQR( 'U', N, N, N, 0, S, WORK( IE ),
$ WORK( IR ), LDWRKR, WORK( IU ),	$ WORK( IR ), LDWRKR, WORK( IU ),
Line 1646	Line 1746
IWORK = ITAU + N	IWORK = ITAU + N
*	*
* Compute A=Q*R, copying result to U	* Compute A=Q*R, copying result to U
* (Workspace: need 2N, prefer N+NNB)	* (Workspace: need 2N, prefer N + NNB)
*	*
CALL DGEQRF( M, N, A, LDA, WORK( ITAU ),	CALL DGEQRF( M, N, A, LDA, WORK( ITAU ),
$ WORK( IWORK ), LWORK-IWORK+1, IERR )	$ WORK( IWORK ), LWORK-IWORK+1, IERR )
CALL DLACPY( 'L', M, N, A, LDA, U, LDU )	CALL DLACPY( 'L', M, N, A, LDA, U, LDU )
*	*
* Generate Q in U	* Generate Q in U
* (Workspace: need N+M, prefer N+M*NB)	* (Workspace: need N + M, prefer N + M*NB)
*	*
CALL DORGQR( M, M, N, U, LDU, WORK( ITAU ),	CALL DORGQR( M, M, N, U, LDU, WORK( ITAU ),
$ WORK( IWORK ), LWORK-IWORK+1, IERR )	$ WORK( IWORK ), LWORK-IWORK+1, IERR )
Line 1664	Line 1764
*	*
* Zero out below R in A	* Zero out below R in A
*	*
CALL DLASET( 'L', N-1, N-1, ZERO, ZERO, A( 2, 1 ),	IF( N .GT. 1 ) THEN
$ LDA )	CALL DLASET( 'L', N-1, N-1, ZERO, ZERO,
	$ A( 2, 1 ), LDA )
	END IF
*	*
* Bidiagonalize R in A	* Bidiagonalize R in A
* (Workspace: need 4N, prefer 3N+2NNB)	* (Workspace: need 4N, prefer 3N + 2NNB)
*	*
CALL DGEBRD( N, N, A, LDA, S, WORK( IE ),	CALL DGEBRD( N, N, A, LDA, S, WORK( IE ),
$ WORK( ITAUQ ), WORK( ITAUP ),	$ WORK( ITAUQ ), WORK( ITAUP ),
Line 1676	Line 1778
*	*
* Multiply Q in U by left bidiagonalizing vectors	* Multiply Q in U by left bidiagonalizing vectors
* in A	* in A
* (Workspace: need 3N+M, prefer 3N+M*NB)	* (Workspace: need 3N + M, prefer 3N + M*NB)
*	*
CALL DORMBR( 'Q', 'R', 'N', M, N, N, A, LDA,	CALL DORMBR( 'Q', 'R', 'N', M, N, N, A, LDA,
$ WORK( ITAUQ ), U, LDU, WORK( IWORK ),	$ WORK( ITAUQ ), U, LDU, WORK( IWORK ),
$ LWORK-IWORK+1, IERR )	$ LWORK-IWORK+1, IERR )
*	*
* Generate right bidiagonalizing vectors in A	* Generate right bidiagonalizing vectors in A
* (Workspace: need 4N-1, prefer 3N+(N-1)*NB)	* (Workspace: need 4N-1, prefer 3N + (N-1)*NB)
*	*
CALL DORGBR( 'P', N, N, N, A, LDA, WORK( ITAUP ),	CALL DORGBR( 'P', N, N, N, A, LDA, WORK( ITAUP ),
$ WORK( IWORK ), LWORK-IWORK+1, IERR )	$ WORK( IWORK ), LWORK-IWORK+1, IERR )
Line 1727	Line 1829
IWORK = ITAU + N	IWORK = ITAU + N
*	*
* Compute A=Q*R, copying result to U	* Compute A=Q*R, copying result to U
* (Workspace: need NN+2N, prefer NN+N+NNB)	* (Workspace: need NN + 2N, prefer NN + N + NNB)
*	*
CALL DGEQRF( M, N, A, LDA, WORK( ITAU ),	CALL DGEQRF( M, N, A, LDA, WORK( ITAU ),
$ WORK( IWORK ), LWORK-IWORK+1, IERR )	$ WORK( IWORK ), LWORK-IWORK+1, IERR )
CALL DLACPY( 'L', M, N, A, LDA, U, LDU )	CALL DLACPY( 'L', M, N, A, LDA, U, LDU )
*	*
* Generate Q in U	* Generate Q in U
* (Workspace: need NN+N+M, prefer NN+N+M*NB)	* (Workspace: need NN + N + M, prefer NN + N + M*NB)
*	*
CALL DORGQR( M, M, N, U, LDU, WORK( ITAU ),	CALL DORGQR( M, M, N, U, LDU, WORK( ITAU ),
$ WORK( IWORK ), LWORK-IWORK+1, IERR )	$ WORK( IWORK ), LWORK-IWORK+1, IERR )
Line 1751	Line 1853
IWORK = ITAUP + N	IWORK = ITAUP + N
*	*
* Bidiagonalize R in WORK(IU), copying result to VT	* Bidiagonalize R in WORK(IU), copying result to VT
* (Workspace: need NN+4N, prefer NN+3N+2NNB)	* (Workspace: need NN + 4N, prefer NN + 3N + 2NNB)
*	*
CALL DGEBRD( N, N, WORK( IU ), LDWRKU, S,	CALL DGEBRD( N, N, WORK( IU ), LDWRKU, S,
$ WORK( IE ), WORK( ITAUQ ),	$ WORK( IE ), WORK( ITAUQ ),
Line 1761	Line 1863
$ LDVT )	$ LDVT )
*	*
* Generate left bidiagonalizing vectors in WORK(IU)	* Generate left bidiagonalizing vectors in WORK(IU)
* (Workspace: need NN+4N, prefer NN+3N+N*NB)	* (Workspace: need NN + 4N, prefer NN + 3N + N*NB)
*	*
CALL DORGBR( 'Q', N, N, N, WORK( IU ), LDWRKU,	CALL DORGBR( 'Q', N, N, N, WORK( IU ), LDWRKU,
$ WORK( ITAUQ ), WORK( IWORK ),	$ WORK( ITAUQ ), WORK( IWORK ),
$ LWORK-IWORK+1, IERR )	$ LWORK-IWORK+1, IERR )
*	*
* Generate right bidiagonalizing vectors in VT	* Generate right bidiagonalizing vectors in VT
* (Workspace: need NN+4N-1,	* (Workspace: need NN + 4N-1,
* prefer NN+3N+(N-1)*NB)	* prefer NN+3N+(N-1)*NB)
*	*
CALL DORGBR( 'P', N, N, N, VT, LDVT, WORK( ITAUP ),	CALL DORGBR( 'P', N, N, N, VT, LDVT, WORK( ITAUP ),
Line 1778	Line 1880
* Perform bidiagonal QR iteration, computing left	* Perform bidiagonal QR iteration, computing left
* singular vectors of R in WORK(IU) and computing	* singular vectors of R in WORK(IU) and computing
* right singular vectors of R in VT	* right singular vectors of R in VT
* (Workspace: need N*N+BDSPAC)	* (Workspace: need N*N + BDSPAC)
*	*
CALL DBDSQR( 'U', N, N, N, 0, S, WORK( IE ), VT,	CALL DBDSQR( 'U', N, N, N, 0, S, WORK( IE ), VT,
$ LDVT, WORK( IU ), LDWRKU, DUM, 1,	$ LDVT, WORK( IU ), LDWRKU, DUM, 1,
Line 1803	Line 1905
IWORK = ITAU + N	IWORK = ITAU + N
*	*
* Compute A=Q*R, copying result to U	* Compute A=Q*R, copying result to U
* (Workspace: need 2N, prefer N+NNB)	* (Workspace: need 2N, prefer N + NNB)
*	*
CALL DGEQRF( M, N, A, LDA, WORK( ITAU ),	CALL DGEQRF( M, N, A, LDA, WORK( ITAU ),
$ WORK( IWORK ), LWORK-IWORK+1, IERR )	$ WORK( IWORK ), LWORK-IWORK+1, IERR )
CALL DLACPY( 'L', M, N, A, LDA, U, LDU )	CALL DLACPY( 'L', M, N, A, LDA, U, LDU )
*	*
* Generate Q in U	* Generate Q in U
* (Workspace: need N+M, prefer N+M*NB)	* (Workspace: need N + M, prefer N + M*NB)
*	*
CALL DORGQR( M, M, N, U, LDU, WORK( ITAU ),	CALL DORGQR( M, M, N, U, LDU, WORK( ITAU ),
$ WORK( IWORK ), LWORK-IWORK+1, IERR )	$ WORK( IWORK ), LWORK-IWORK+1, IERR )
Line 1827	Line 1929
IWORK = ITAUP + N	IWORK = ITAUP + N
*	*
* Bidiagonalize R in VT	* Bidiagonalize R in VT
* (Workspace: need 4N, prefer 3N+2NNB)	* (Workspace: need 4N, prefer 3N + 2NNB)
*	*
CALL DGEBRD( N, N, VT, LDVT, S, WORK( IE ),	CALL DGEBRD( N, N, VT, LDVT, S, WORK( IE ),
$ WORK( ITAUQ ), WORK( ITAUP ),	$ WORK( ITAUQ ), WORK( ITAUP ),
Line 1835	Line 1937
*	*
* Multiply Q in U by left bidiagonalizing vectors	* Multiply Q in U by left bidiagonalizing vectors
* in VT	* in VT
* (Workspace: need 3N+M, prefer 3N+M*NB)	* (Workspace: need 3N + M, prefer 3N + M*NB)
*	*
CALL DORMBR( 'Q', 'R', 'N', M, N, N, VT, LDVT,	CALL DORMBR( 'Q', 'R', 'N', M, N, N, VT, LDVT,
$ WORK( ITAUQ ), U, LDU, WORK( IWORK ),	$ WORK( ITAUQ ), U, LDU, WORK( IWORK ),
$ LWORK-IWORK+1, IERR )	$ LWORK-IWORK+1, IERR )
*	*
* Generate right bidiagonalizing vectors in VT	* Generate right bidiagonalizing vectors in VT
* (Workspace: need 4N-1, prefer 3N+(N-1)*NB)	* (Workspace: need 4N-1, prefer 3N + (N-1)*NB)
*	*
CALL DORGBR( 'P', N, N, N, VT, LDVT, WORK( ITAUP ),	CALL DORGBR( 'P', N, N, N, VT, LDVT, WORK( ITAUP ),
$ WORK( IWORK ), LWORK-IWORK+1, IERR )	$ WORK( IWORK ), LWORK-IWORK+1, IERR )
Line 1876	Line 1978
IWORK = ITAUP + N	IWORK = ITAUP + N
*	*
* Bidiagonalize A	* Bidiagonalize A
* (Workspace: need 3N+M, prefer 3N+(M+N)*NB)	* (Workspace: need 3N + M, prefer 3N + (M + N)*NB)
*	*
CALL DGEBRD( M, N, A, LDA, S, WORK( IE ), WORK( ITAUQ ),	CALL DGEBRD( M, N, A, LDA, S, WORK( IE ), WORK( ITAUQ ),
$ WORK( ITAUP ), WORK( IWORK ), LWORK-IWORK+1,	$ WORK( ITAUP ), WORK( IWORK ), LWORK-IWORK+1,
Line 1885	Line 1987
*	*
* If left singular vectors desired in U, copy result to U	* If left singular vectors desired in U, copy result to U
* and generate left bidiagonalizing vectors in U	* and generate left bidiagonalizing vectors in U
* (Workspace: need 3N+NCU, prefer 3N+NCU*NB)	* (Workspace: need 3N + NCU, prefer 3N + NCU*NB)
*	*
CALL DLACPY( 'L', M, N, A, LDA, U, LDU )	CALL DLACPY( 'L', M, N, A, LDA, U, LDU )
IF( WNTUS )	IF( WNTUS )
Line 1899	Line 2001
*	*
* If right singular vectors desired in VT, copy result to	* If right singular vectors desired in VT, copy result to
* VT and generate right bidiagonalizing vectors in VT	* VT and generate right bidiagonalizing vectors in VT
* (Workspace: need 4N-1, prefer 3N+(N-1)*NB)	* (Workspace: need 4N-1, prefer 3N + (N-1)*NB)
*	*
CALL DLACPY( 'U', N, N, A, LDA, VT, LDVT )	CALL DLACPY( 'U', N, N, A, LDA, VT, LDVT )
CALL DORGBR( 'P', N, N, N, VT, LDVT, WORK( ITAUP ),	CALL DORGBR( 'P', N, N, N, VT, LDVT, WORK( ITAUP ),
Line 1909	Line 2011
*	*
* If left singular vectors desired in A, generate left	* If left singular vectors desired in A, generate left
* bidiagonalizing vectors in A	* bidiagonalizing vectors in A
* (Workspace: need 4N, prefer 3N+N*NB)	* (Workspace: need 4N, prefer 3N + N*NB)
*	*
CALL DORGBR( 'Q', M, N, N, A, LDA, WORK( ITAUQ ),	CALL DORGBR( 'Q', M, N, N, A, LDA, WORK( ITAUQ ),
$ WORK( IWORK ), LWORK-IWORK+1, IERR )	$ WORK( IWORK ), LWORK-IWORK+1, IERR )
Line 1918	Line 2020
*	*
* If right singular vectors desired in A, generate right	* If right singular vectors desired in A, generate right
* bidiagonalizing vectors in A	* bidiagonalizing vectors in A
* (Workspace: need 4N-1, prefer 3N+(N-1)*NB)	* (Workspace: need 4N-1, prefer 3N + (N-1)*NB)
*	*
CALL DORGBR( 'P', N, N, N, A, LDA, WORK( ITAUP ),	CALL DORGBR( 'P', N, N, N, A, LDA, WORK( ITAUP ),
$ WORK( IWORK ), LWORK-IWORK+1, IERR )	$ WORK( IWORK ), LWORK-IWORK+1, IERR )
Line 1980	Line 2082
IWORK = ITAU + M	IWORK = ITAU + M
*	*
* Compute A=L*Q	* Compute A=L*Q
* (Workspace: need 2M, prefer M+MNB)	* (Workspace: need 2M, prefer M + MNB)
*	*
CALL DGELQF( M, N, A, LDA, WORK( ITAU ), WORK( IWORK ),	CALL DGELQF( M, N, A, LDA, WORK( ITAU ), WORK( IWORK ),
$ LWORK-IWORK+1, IERR )	$ LWORK-IWORK+1, IERR )
Line 1994	Line 2096
IWORK = ITAUP + M	IWORK = ITAUP + M
*	*
* Bidiagonalize L in A	* Bidiagonalize L in A
* (Workspace: need 4M, prefer 3M+2MNB)	* (Workspace: need 4M, prefer 3M + 2MNB)
*	*
CALL DGEBRD( M, M, A, LDA, S, WORK( IE ), WORK( ITAUQ ),	CALL DGEBRD( M, M, A, LDA, S, WORK( IE ), WORK( ITAUQ ),
$ WORK( ITAUP ), WORK( IWORK ), LWORK-IWORK+1,	$ WORK( ITAUP ), WORK( IWORK ), LWORK-IWORK+1,
Line 2002	Line 2104
IF( WNTUO .OR. WNTUAS ) THEN	IF( WNTUO .OR. WNTUAS ) THEN
*	*
* If left singular vectors desired, generate Q	* If left singular vectors desired, generate Q
* (Workspace: need 4M, prefer 3M+M*NB)	* (Workspace: need 4M, prefer 3M + M*NB)
*	*
CALL DORGBR( 'Q', M, M, M, A, LDA, WORK( ITAUQ ),	CALL DORGBR( 'Q', M, M, M, A, LDA, WORK( ITAUQ ),
$ WORK( IWORK ), LWORK-IWORK+1, IERR )	$ WORK( IWORK ), LWORK-IWORK+1, IERR )
Line 2035	Line 2137
* Sufficient workspace for a fast algorithm	* Sufficient workspace for a fast algorithm
*	*
IR = 1	IR = 1
IF( LWORK.GE.MAX( WRKBL, LDAN+M )+LDAM ) THEN	IF( LWORK.GE.MAX( WRKBL, LDAN + M ) + LDAM ) THEN
*	*
* WORK(IU) is LDA by N and WORK(IR) is LDA by M	* WORK(IU) is LDA by N and WORK(IR) is LDA by M
*	*
LDWRKU = LDA	LDWRKU = LDA
CHUNK = N	CHUNK = N
LDWRKR = LDA	LDWRKR = LDA
ELSE IF( LWORK.GE.MAX( WRKBL, LDAN+M )+MM ) THEN	ELSE IF( LWORK.GE.MAX( WRKBL, LDAN + M ) + MM ) THEN
*	*
* WORK(IU) is LDA by N and WORK(IR) is M by M	* WORK(IU) is LDA by N and WORK(IR) is M by M
*	*
Line 2061	Line 2163
IWORK = ITAU + M	IWORK = ITAU + M
*	*
* Compute A=L*Q	* Compute A=L*Q
* (Workspace: need MM+2M, prefer MM+M+MNB)	* (Workspace: need MM + 2M, prefer MM + M + MNB)
*	*
CALL DGELQF( M, N, A, LDA, WORK( ITAU ),	CALL DGELQF( M, N, A, LDA, WORK( ITAU ),
$ WORK( IWORK ), LWORK-IWORK+1, IERR )	$ WORK( IWORK ), LWORK-IWORK+1, IERR )
Line 2073	Line 2175
$ WORK( IR+LDWRKR ), LDWRKR )	$ WORK( IR+LDWRKR ), LDWRKR )
*	*
* Generate Q in A	* Generate Q in A
* (Workspace: need MM+2M, prefer MM+M+MNB)	* (Workspace: need MM + 2M, prefer MM + M + MNB)
*	*
CALL DORGLQ( M, N, M, A, LDA, WORK( ITAU ),	CALL DORGLQ( M, N, M, A, LDA, WORK( ITAU ),
$ WORK( IWORK ), LWORK-IWORK+1, IERR )	$ WORK( IWORK ), LWORK-IWORK+1, IERR )
Line 2083	Line 2185
IWORK = ITAUP + M	IWORK = ITAUP + M
*	*
* Bidiagonalize L in WORK(IR)	* Bidiagonalize L in WORK(IR)
* (Workspace: need MM+4M, prefer MM+3M+2MNB)	* (Workspace: need MM + 4M, prefer MM + 3M + 2MNB)
*	*
CALL DGEBRD( M, M, WORK( IR ), LDWRKR, S, WORK( IE ),	CALL DGEBRD( M, M, WORK( IR ), LDWRKR, S, WORK( IE ),
$ WORK( ITAUQ ), WORK( ITAUP ),	$ WORK( ITAUQ ), WORK( ITAUP ),
$ WORK( IWORK ), LWORK-IWORK+1, IERR )	$ WORK( IWORK ), LWORK-IWORK+1, IERR )
*	*
* Generate right vectors bidiagonalizing L	* Generate right vectors bidiagonalizing L
* (Workspace: need MM+4M-1, prefer MM+3M+(M-1)*NB)	* (Workspace: need MM + 4M-1, prefer MM + 3M + (M-1)*NB)
*	*
CALL DORGBR( 'P', M, M, M, WORK( IR ), LDWRKR,	CALL DORGBR( 'P', M, M, M, WORK( IR ), LDWRKR,
$ WORK( ITAUP ), WORK( IWORK ),	$ WORK( ITAUP ), WORK( IWORK ),
Line 2099	Line 2201
*	*
* Perform bidiagonal QR iteration, computing right	* Perform bidiagonal QR iteration, computing right
* singular vectors of L in WORK(IR)	* singular vectors of L in WORK(IR)
* (Workspace: need M*M+BDSPAC)	* (Workspace: need M*M + BDSPAC)
*	*
CALL DBDSQR( 'U', M, M, 0, 0, S, WORK( IE ),	CALL DBDSQR( 'U', M, M, 0, 0, S, WORK( IE ),
$ WORK( IR ), LDWRKR, DUM, 1, DUM, 1,	$ WORK( IR ), LDWRKR, DUM, 1, DUM, 1,
Line 2108	Line 2210
*	*
* Multiply right singular vectors of L in WORK(IR) by Q	* Multiply right singular vectors of L in WORK(IR) by Q
* in A, storing result in WORK(IU) and copying to A	* in A, storing result in WORK(IU) and copying to A
* (Workspace: need MM+2M, prefer MM+MN+M)	* (Workspace: need MM + 2M, prefer MM + MN + M)
*	*
DO 30 I = 1, N, CHUNK	DO 30 I = 1, N, CHUNK
BLK = MIN( N-I+1, CHUNK )	BLK = MIN( N-I+1, CHUNK )
Line 2129	Line 2231
IWORK = ITAUP + M	IWORK = ITAUP + M
*	*
* Bidiagonalize A	* Bidiagonalize A
* (Workspace: need 3M+N, prefer 3M+(M+N)*NB)	* (Workspace: need 3M + N, prefer 3M + (M + N)*NB)
*	*
CALL DGEBRD( M, N, A, LDA, S, WORK( IE ),	CALL DGEBRD( M, N, A, LDA, S, WORK( IE ),
$ WORK( ITAUQ ), WORK( ITAUP ),	$ WORK( ITAUQ ), WORK( ITAUP ),
$ WORK( IWORK ), LWORK-IWORK+1, IERR )	$ WORK( IWORK ), LWORK-IWORK+1, IERR )
*	*
* Generate right vectors bidiagonalizing A	* Generate right vectors bidiagonalizing A
* (Workspace: need 4M, prefer 3M+M*NB)	* (Workspace: need 4M, prefer 3M + M*NB)
*	*
CALL DORGBR( 'P', M, N, M, A, LDA, WORK( ITAUP ),	CALL DORGBR( 'P', M, N, M, A, LDA, WORK( ITAUP ),
$ WORK( IWORK ), LWORK-IWORK+1, IERR )	$ WORK( IWORK ), LWORK-IWORK+1, IERR )
Line 2162	Line 2264
* Sufficient workspace for a fast algorithm	* Sufficient workspace for a fast algorithm
*	*
IR = 1	IR = 1
IF( LWORK.GE.MAX( WRKBL, LDAN+M )+LDAM ) THEN	IF( LWORK.GE.MAX( WRKBL, LDAN + M ) + LDAM ) THEN
*	*
* WORK(IU) is LDA by N and WORK(IR) is LDA by M	* WORK(IU) is LDA by N and WORK(IR) is LDA by M
*	*
LDWRKU = LDA	LDWRKU = LDA
CHUNK = N	CHUNK = N
LDWRKR = LDA	LDWRKR = LDA
ELSE IF( LWORK.GE.MAX( WRKBL, LDAN+M )+MM ) THEN	ELSE IF( LWORK.GE.MAX( WRKBL, LDAN + M ) + MM ) THEN
*	*
* WORK(IU) is LDA by N and WORK(IR) is M by M	* WORK(IU) is LDA by N and WORK(IR) is M by M
*	*
Line 2188	Line 2290
IWORK = ITAU + M	IWORK = ITAU + M
*	*
* Compute A=L*Q	* Compute A=L*Q
* (Workspace: need MM+2M, prefer MM+M+MNB)	* (Workspace: need MM + 2M, prefer MM + M + MNB)
*	*
CALL DGELQF( M, N, A, LDA, WORK( ITAU ),	CALL DGELQF( M, N, A, LDA, WORK( ITAU ),
$ WORK( IWORK ), LWORK-IWORK+1, IERR )	$ WORK( IWORK ), LWORK-IWORK+1, IERR )
Line 2200	Line 2302
$ LDU )	$ LDU )
*	*
* Generate Q in A	* Generate Q in A
* (Workspace: need MM+2M, prefer MM+M+MNB)	* (Workspace: need MM + 2M, prefer MM + M + MNB)
*	*
CALL DORGLQ( M, N, M, A, LDA, WORK( ITAU ),	CALL DORGLQ( M, N, M, A, LDA, WORK( ITAU ),
$ WORK( IWORK ), LWORK-IWORK+1, IERR )	$ WORK( IWORK ), LWORK-IWORK+1, IERR )
Line 2210	Line 2312
IWORK = ITAUP + M	IWORK = ITAUP + M
*	*
* Bidiagonalize L in U, copying result to WORK(IR)	* Bidiagonalize L in U, copying result to WORK(IR)
* (Workspace: need MM+4M, prefer MM+3M+2MNB)	* (Workspace: need MM + 4M, prefer MM + 3M + 2MNB)
*	*
CALL DGEBRD( M, M, U, LDU, S, WORK( IE ),	CALL DGEBRD( M, M, U, LDU, S, WORK( IE ),
$ WORK( ITAUQ ), WORK( ITAUP ),	$ WORK( ITAUQ ), WORK( ITAUP ),
Line 2218	Line 2320
CALL DLACPY( 'U', M, M, U, LDU, WORK( IR ), LDWRKR )	CALL DLACPY( 'U', M, M, U, LDU, WORK( IR ), LDWRKR )
*	*
* Generate right vectors bidiagonalizing L in WORK(IR)	* Generate right vectors bidiagonalizing L in WORK(IR)
* (Workspace: need MM+4M-1, prefer MM+3M+(M-1)*NB)	* (Workspace: need MM + 4M-1, prefer MM + 3M + (M-1)*NB)
*	*
CALL DORGBR( 'P', M, M, M, WORK( IR ), LDWRKR,	CALL DORGBR( 'P', M, M, M, WORK( IR ), LDWRKR,
$ WORK( ITAUP ), WORK( IWORK ),	$ WORK( ITAUP ), WORK( IWORK ),
$ LWORK-IWORK+1, IERR )	$ LWORK-IWORK+1, IERR )
*	*
* Generate left vectors bidiagonalizing L in U	* Generate left vectors bidiagonalizing L in U
* (Workspace: need MM+4M, prefer MM+3M+M*NB)	* (Workspace: need MM + 4M, prefer MM + 3M + M*NB)
*	*
CALL DORGBR( 'Q', M, M, M, U, LDU, WORK( ITAUQ ),	CALL DORGBR( 'Q', M, M, M, U, LDU, WORK( ITAUQ ),
$ WORK( IWORK ), LWORK-IWORK+1, IERR )	$ WORK( IWORK ), LWORK-IWORK+1, IERR )
Line 2234	Line 2336
* Perform bidiagonal QR iteration, computing left	* Perform bidiagonal QR iteration, computing left
* singular vectors of L in U, and computing right	* singular vectors of L in U, and computing right
* singular vectors of L in WORK(IR)	* singular vectors of L in WORK(IR)
* (Workspace: need M*M+BDSPAC)	* (Workspace: need M*M + BDSPAC)
*	*
CALL DBDSQR( 'U', M, M, M, 0, S, WORK( IE ),	CALL DBDSQR( 'U', M, M, M, 0, S, WORK( IE ),
$ WORK( IR ), LDWRKR, U, LDU, DUM, 1,	$ WORK( IR ), LDWRKR, U, LDU, DUM, 1,
Line 2243	Line 2345
*	*
* Multiply right singular vectors of L in WORK(IR) by Q	* Multiply right singular vectors of L in WORK(IR) by Q
* in A, storing result in WORK(IU) and copying to A	* in A, storing result in WORK(IU) and copying to A
* (Workspace: need MM+2M, prefer MM+MN+M))	* (Workspace: need MM + 2M, prefer MM + MN + M))
*	*
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 2262	Line 2364
IWORK = ITAU + M	IWORK = ITAU + M
*	*
* Compute A=L*Q	* Compute A=L*Q
* (Workspace: need 2M, prefer M+MNB)	* (Workspace: need 2M, prefer M + MNB)
*	*
CALL DGELQF( M, N, A, LDA, WORK( ITAU ),	CALL DGELQF( M, N, A, LDA, WORK( ITAU ),
$ WORK( IWORK ), LWORK-IWORK+1, IERR )	$ WORK( IWORK ), LWORK-IWORK+1, IERR )
Line 2274	Line 2376
$ LDU )	$ LDU )
*	*
* Generate Q in A	* Generate Q in A
* (Workspace: need 2M, prefer M+MNB)	* (Workspace: need 2M, prefer M + MNB)
*	*
CALL DORGLQ( M, N, M, A, LDA, WORK( ITAU ),	CALL DORGLQ( M, N, M, A, LDA, WORK( ITAU ),
$ WORK( IWORK ), LWORK-IWORK+1, IERR )	$ WORK( IWORK ), LWORK-IWORK+1, IERR )
Line 2284	Line 2386
IWORK = ITAUP + M	IWORK = ITAUP + M
*	*
* Bidiagonalize L in U	* Bidiagonalize L in U
* (Workspace: need 4M, prefer 3M+2MNB)	* (Workspace: need 4M, prefer 3M + 2MNB)
*	*
CALL DGEBRD( M, M, U, LDU, S, WORK( IE ),	CALL DGEBRD( M, M, U, LDU, S, WORK( IE ),
$ WORK( ITAUQ ), WORK( ITAUP ),	$ WORK( ITAUQ ), WORK( ITAUP ),
$ WORK( IWORK ), LWORK-IWORK+1, IERR )	$ WORK( IWORK ), LWORK-IWORK+1, IERR )
*	*
* Multiply right vectors bidiagonalizing L by Q in A	* Multiply right vectors bidiagonalizing L by Q in A
* (Workspace: need 3M+N, prefer 3M+N*NB)	* (Workspace: need 3M + N, prefer 3M + N*NB)
*	*
CALL DORMBR( 'P', 'L', 'T', M, N, M, U, LDU,	CALL DORMBR( 'P', 'L', 'T', M, N, M, U, LDU,
$ WORK( ITAUP ), A, LDA, WORK( IWORK ),	$ WORK( ITAUP ), A, LDA, WORK( IWORK ),
$ LWORK-IWORK+1, IERR )	$ LWORK-IWORK+1, IERR )
*	*
* Generate left vectors bidiagonalizing L in U	* Generate left vectors bidiagonalizing L in U
* (Workspace: need 4M, prefer 3M+M*NB)	* (Workspace: need 4M, prefer 3M + M*NB)
*	*
CALL DORGBR( 'Q', M, M, M, U, LDU, WORK( ITAUQ ),	CALL DORGBR( 'Q', M, M, M, U, LDU, WORK( ITAUQ ),
$ WORK( IWORK ), LWORK-IWORK+1, IERR )	$ WORK( IWORK ), LWORK-IWORK+1, IERR )
Line 2342	Line 2444
IWORK = ITAU + M	IWORK = ITAU + M
*	*
* Compute A=L*Q	* Compute A=L*Q
* (Workspace: need MM+2M, prefer MM+M+MNB)	* (Workspace: need MM + 2M, prefer MM + M + MNB)
*	*
CALL DGELQF( M, N, A, LDA, WORK( ITAU ),	CALL DGELQF( M, N, A, LDA, WORK( ITAU ),
$ WORK( IWORK ), LWORK-IWORK+1, IERR )	$ WORK( IWORK ), LWORK-IWORK+1, IERR )
Line 2355	Line 2457
$ WORK( IR+LDWRKR ), LDWRKR )	$ WORK( IR+LDWRKR ), LDWRKR )
*	*
* Generate Q in A	* Generate Q in A
* (Workspace: need MM+2M, prefer MM+M+MNB)	* (Workspace: need MM + 2M, prefer MM + M + MNB)
*	*
CALL DORGLQ( M, N, M, A, LDA, WORK( ITAU ),	CALL DORGLQ( M, N, M, A, LDA, WORK( ITAU ),
$ WORK( IWORK ), LWORK-IWORK+1, IERR )	$ WORK( IWORK ), LWORK-IWORK+1, IERR )
Line 2365	Line 2467
IWORK = ITAUP + M	IWORK = ITAUP + M
*	*
* Bidiagonalize L in WORK(IR)	* Bidiagonalize L in WORK(IR)
* (Workspace: need MM+4M, prefer MM+3M+2MNB)	* (Workspace: need MM + 4M, prefer MM + 3M + 2MNB)
*	*
CALL DGEBRD( M, M, WORK( IR ), LDWRKR, S,	CALL DGEBRD( M, M, WORK( IR ), LDWRKR, S,
$ WORK( IE ), WORK( ITAUQ ),	$ WORK( IE ), WORK( ITAUQ ),
Line 2374	Line 2476
*	*
* Generate right vectors bidiagonalizing L in	* Generate right vectors bidiagonalizing L in
* WORK(IR)	* WORK(IR)
* (Workspace: need MM+4M, prefer MM+3M+(M-1)*NB)	* (Workspace: need MM + 4M, prefer MM + 3M + (M-1)*NB)
*	*
CALL DORGBR( 'P', M, M, M, WORK( IR ), LDWRKR,	CALL DORGBR( 'P', M, M, M, WORK( IR ), LDWRKR,
$ WORK( ITAUP ), WORK( IWORK ),	$ WORK( ITAUP ), WORK( IWORK ),
Line 2383	Line 2485
*	*
* Perform bidiagonal QR iteration, computing right	* Perform bidiagonal QR iteration, computing right
* singular vectors of L in WORK(IR)	* singular vectors of L in WORK(IR)
* (Workspace: need M*M+BDSPAC)	* (Workspace: need M*M + BDSPAC)
*	*
CALL DBDSQR( 'U', M, M, 0, 0, S, WORK( IE ),	CALL DBDSQR( 'U', M, M, 0, 0, S, WORK( IE ),
$ WORK( IR ), LDWRKR, DUM, 1, DUM, 1,	$ WORK( IR ), LDWRKR, DUM, 1, DUM, 1,
Line 2404	Line 2506
IWORK = ITAU + M	IWORK = ITAU + M
*	*
* Compute A=L*Q	* Compute A=L*Q
* (Workspace: need 2M, prefer M+MNB)	* (Workspace: need 2M, prefer M + MNB)
*	*
CALL DGELQF( M, N, A, LDA, WORK( ITAU ),	CALL DGELQF( M, N, A, LDA, WORK( ITAU ),
$ WORK( IWORK ), LWORK-IWORK+1, IERR )	$ WORK( IWORK ), LWORK-IWORK+1, IERR )
Line 2414	Line 2516
CALL DLACPY( 'U', M, N, A, LDA, VT, LDVT )	CALL DLACPY( 'U', M, N, A, LDA, VT, LDVT )
*	*
* Generate Q in VT	* Generate Q in VT
* (Workspace: need 2M, prefer M+MNB)	* (Workspace: need 2M, prefer M + MNB)
*	*
CALL DORGLQ( M, N, M, VT, LDVT, WORK( ITAU ),	CALL DORGLQ( M, N, M, VT, LDVT, WORK( ITAU ),
$ WORK( IWORK ), LWORK-IWORK+1, IERR )	$ WORK( IWORK ), LWORK-IWORK+1, IERR )
Line 2429	Line 2531
$ LDA )	$ LDA )
*	*
* Bidiagonalize L in A	* Bidiagonalize L in A
* (Workspace: need 4M, prefer 3M+2MNB)	* (Workspace: need 4M, prefer 3M + 2MNB)
*	*
CALL DGEBRD( M, M, A, LDA, S, WORK( IE ),	CALL DGEBRD( M, M, A, LDA, S, WORK( IE ),
$ WORK( ITAUQ ), WORK( ITAUP ),	$ WORK( ITAUQ ), WORK( ITAUP ),
$ WORK( IWORK ), LWORK-IWORK+1, IERR )	$ WORK( IWORK ), LWORK-IWORK+1, IERR )
*	*
* Multiply right vectors bidiagonalizing L by Q in VT	* Multiply right vectors bidiagonalizing L by Q in VT
* (Workspace: need 3M+N, prefer 3M+N*NB)	* (Workspace: need 3M + N, prefer 3M + N*NB)
*	*
CALL DORMBR( 'P', 'L', 'T', M, N, M, A, LDA,	CALL DORMBR( 'P', 'L', 'T', M, N, M, A, LDA,
$ WORK( ITAUP ), VT, LDVT,	$ WORK( ITAUP ), VT, LDVT,
Line 2471	Line 2573
LDWRKU = LDA	LDWRKU = LDA
IR = IU + LDWRKU*M	IR = IU + LDWRKU*M
LDWRKR = LDA	LDWRKR = LDA
ELSE IF( LWORK.GE.WRKBL+( LDA+M )*M ) THEN	ELSE IF( LWORK.GE.WRKBL+( LDA + M )*M ) THEN
*	*
* WORK(IU) is LDA by M and WORK(IR) is M by M	* WORK(IU) is LDA by M and WORK(IR) is M by M
*	*
Line 2490	Line 2592
IWORK = ITAU + M	IWORK = ITAU + M
*	*
* Compute A=L*Q	* Compute A=L*Q
* (Workspace: need 2MM+2M, prefer 2MM+M+MNB)	* (Workspace: need 2MM + 2M, prefer 2MM + M + MNB)
*	*
CALL DGELQF( M, N, A, LDA, WORK( ITAU ),	CALL DGELQF( M, N, A, LDA, WORK( ITAU ),
$ WORK( IWORK ), LWORK-IWORK+1, IERR )	$ WORK( IWORK ), LWORK-IWORK+1, IERR )
Line 2503	Line 2605
$ WORK( IU+LDWRKU ), LDWRKU )	$ WORK( IU+LDWRKU ), LDWRKU )
*	*
* Generate Q in A	* Generate Q in A
* (Workspace: need 2MM+2M, prefer 2MM+M+MNB)	* (Workspace: need 2MM + 2M, prefer 2MM + M + MNB)
*	*
CALL DORGLQ( M, N, M, A, LDA, WORK( ITAU ),	CALL DORGLQ( M, N, M, A, LDA, WORK( ITAU ),
$ WORK( IWORK ), LWORK-IWORK+1, IERR )	$ WORK( IWORK ), LWORK-IWORK+1, IERR )
Line 2514	Line 2616
*	*
* Bidiagonalize L in WORK(IU), copying result to	* Bidiagonalize L in WORK(IU), copying result to
* WORK(IR)	* WORK(IR)
* (Workspace: need 2MM+4*M,	* (Workspace: need 2MM + 4*M,
* prefer 2MM+3M+2M*NB)	* prefer 2MM+3M+2M*NB)
*	*
CALL DGEBRD( M, M, WORK( IU ), LDWRKU, S,	CALL DGEBRD( M, M, WORK( IU ), LDWRKU, S,
Line 2525	Line 2627
$ WORK( IR ), LDWRKR )	$ WORK( IR ), LDWRKR )
*	*
* Generate right bidiagonalizing vectors in WORK(IU)	* Generate right bidiagonalizing vectors in WORK(IU)
* (Workspace: need 2MM+4*M-1,	* (Workspace: need 2MM + 4*M-1,
* prefer 2MM+3M+(M-1)NB)	* prefer 2MM+3M+(M-1)NB)
*	*
CALL DORGBR( 'P', M, M, M, WORK( IU ), LDWRKU,	CALL DORGBR( 'P', M, M, M, WORK( IU ), LDWRKU,
Line 2533	Line 2635
$ LWORK-IWORK+1, IERR )	$ LWORK-IWORK+1, IERR )
*	*
* Generate left bidiagonalizing vectors in WORK(IR)	* Generate left bidiagonalizing vectors in WORK(IR)
* (Workspace: need 2MM+4M, prefer 2MM+3M+M*NB)	* (Workspace: need 2MM + 4M, prefer 2MM + 3M + M*NB)
*	*
CALL DORGBR( 'Q', M, M, M, WORK( IR ), LDWRKR,	CALL DORGBR( 'Q', M, M, M, WORK( IR ), LDWRKR,
$ WORK( ITAUQ ), WORK( IWORK ),	$ WORK( ITAUQ ), WORK( IWORK ),
Line 2543	Line 2645
* Perform bidiagonal QR iteration, computing left	* Perform bidiagonal QR iteration, computing left
* singular vectors of L in WORK(IR) and computing	* singular vectors of L in WORK(IR) and computing
* right singular vectors of L in WORK(IU)	* right singular vectors of L in WORK(IU)
* (Workspace: need 2MM+BDSPAC)	* (Workspace: need 2MM + BDSPAC)
*	*
CALL DBDSQR( 'U', M, M, M, 0, S, WORK( IE ),	CALL DBDSQR( 'U', M, M, M, 0, S, WORK( IE ),
$ WORK( IU ), LDWRKU, WORK( IR ),	$ WORK( IU ), LDWRKU, WORK( IR ),
Line 2570	Line 2672
IWORK = ITAU + M	IWORK = ITAU + M
*	*
* Compute A=L*Q, copying result to VT	* Compute A=L*Q, copying result to VT
* (Workspace: need 2M, prefer M+MNB)	* (Workspace: need 2M, prefer M + MNB)
*	*
CALL DGELQF( M, N, A, LDA, WORK( ITAU ),	CALL DGELQF( M, N, A, LDA, WORK( ITAU ),
$ WORK( IWORK ), LWORK-IWORK+1, IERR )	$ WORK( IWORK ), LWORK-IWORK+1, IERR )
CALL DLACPY( 'U', M, N, A, LDA, VT, LDVT )	CALL DLACPY( 'U', M, N, A, LDA, VT, LDVT )
*	*
* Generate Q in VT	* Generate Q in VT
* (Workspace: need 2M, prefer M+MNB)	* (Workspace: need 2M, prefer M + MNB)
*	*
CALL DORGLQ( M, N, M, VT, LDVT, WORK( ITAU ),	CALL DORGLQ( M, N, M, VT, LDVT, WORK( ITAU ),
$ WORK( IWORK ), LWORK-IWORK+1, IERR )	$ WORK( IWORK ), LWORK-IWORK+1, IERR )
Line 2592	Line 2694
$ LDA )	$ LDA )
*	*
* Bidiagonalize L in A	* Bidiagonalize L in A
* (Workspace: need 4M, prefer 3M+2MNB)	* (Workspace: need 4M, prefer 3M + 2MNB)
*	*
CALL DGEBRD( M, M, A, LDA, S, WORK( IE ),	CALL DGEBRD( M, M, A, LDA, S, WORK( IE ),
$ WORK( ITAUQ ), WORK( ITAUP ),	$ WORK( ITAUQ ), WORK( ITAUP ),
$ WORK( IWORK ), LWORK-IWORK+1, IERR )	$ WORK( IWORK ), LWORK-IWORK+1, IERR )
*	*
* Multiply right vectors bidiagonalizing L by Q in VT	* Multiply right vectors bidiagonalizing L by Q in VT
* (Workspace: need 3M+N, prefer 3M+N*NB)	* (Workspace: need 3M + N, prefer 3M + N*NB)
*	*
CALL DORMBR( 'P', 'L', 'T', M, N, M, A, LDA,	CALL DORMBR( 'P', 'L', 'T', M, N, M, A, LDA,
$ WORK( ITAUP ), VT, LDVT,	$ WORK( ITAUP ), VT, LDVT,
$ WORK( IWORK ), LWORK-IWORK+1, IERR )	$ WORK( IWORK ), LWORK-IWORK+1, IERR )
*	*
* Generate left bidiagonalizing vectors of L in A	* Generate left bidiagonalizing vectors of L in A
* (Workspace: need 4M, prefer 3M+M*NB)	* (Workspace: need 4M, prefer 3M + M*NB)
*	*
CALL DORGBR( 'Q', M, M, M, A, LDA, WORK( ITAUQ ),	CALL DORGBR( 'Q', M, M, M, A, LDA, WORK( ITAUQ ),
$ WORK( IWORK ), LWORK-IWORK+1, IERR )	$ WORK( IWORK ), LWORK-IWORK+1, IERR )
Line 2650	Line 2752
IWORK = ITAU + M	IWORK = ITAU + M
*	*
* Compute A=L*Q	* Compute A=L*Q
* (Workspace: need MM+2M, prefer MM+M+MNB)	* (Workspace: need MM + 2M, prefer MM + M + MNB)
*	*
CALL DGELQF( M, N, A, LDA, WORK( ITAU ),	CALL DGELQF( M, N, A, LDA, WORK( ITAU ),
$ WORK( IWORK ), LWORK-IWORK+1, IERR )	$ WORK( IWORK ), LWORK-IWORK+1, IERR )
Line 2663	Line 2765
$ WORK( IU+LDWRKU ), LDWRKU )	$ WORK( IU+LDWRKU ), LDWRKU )
*	*
* Generate Q in A	* Generate Q in A
* (Workspace: need MM+2M, prefer MM+M+MNB)	* (Workspace: need MM + 2M, prefer MM + M + MNB)
*	*
CALL DORGLQ( M, N, M, A, LDA, WORK( ITAU ),	CALL DORGLQ( M, N, M, A, LDA, WORK( ITAU ),
$ WORK( IWORK ), LWORK-IWORK+1, IERR )	$ WORK( IWORK ), LWORK-IWORK+1, IERR )
Line 2673	Line 2775
IWORK = ITAUP + M	IWORK = ITAUP + M
*	*
* Bidiagonalize L in WORK(IU), copying result to U	* Bidiagonalize L in WORK(IU), copying result to U
* (Workspace: need MM+4M, prefer MM+3M+2MNB)	* (Workspace: need MM + 4M, prefer MM + 3M + 2MNB)
*	*
CALL DGEBRD( M, M, WORK( IU ), LDWRKU, S,	CALL DGEBRD( M, M, WORK( IU ), LDWRKU, S,
$ WORK( IE ), WORK( ITAUQ ),	$ WORK( IE ), WORK( ITAUQ ),
Line 2683	Line 2785
$ LDU )	$ LDU )
*	*
* Generate right bidiagonalizing vectors in WORK(IU)	* Generate right bidiagonalizing vectors in WORK(IU)
* (Workspace: need MM+4M-1,	* (Workspace: need MM + 4M-1,
* prefer MM+3M+(M-1)*NB)	* prefer MM+3M+(M-1)*NB)
*	*
CALL DORGBR( 'P', M, M, M, WORK( IU ), LDWRKU,	CALL DORGBR( 'P', M, M, M, WORK( IU ), LDWRKU,
Line 2691	Line 2793
$ LWORK-IWORK+1, IERR )	$ LWORK-IWORK+1, IERR )
*	*
* Generate left bidiagonalizing vectors in U	* Generate left bidiagonalizing vectors in U
* (Workspace: need MM+4M, prefer MM+3M+M*NB)	* (Workspace: need MM + 4M, prefer MM + 3M + M*NB)
*	*
CALL DORGBR( 'Q', M, M, M, U, LDU, WORK( ITAUQ ),	CALL DORGBR( 'Q', M, M, M, U, LDU, WORK( ITAUQ ),
$ WORK( IWORK ), LWORK-IWORK+1, IERR )	$ WORK( IWORK ), LWORK-IWORK+1, IERR )
Line 2700	Line 2802
* Perform bidiagonal QR iteration, computing left	* Perform bidiagonal QR iteration, computing left
* singular vectors of L in U and computing right	* singular vectors of L in U and computing right
* singular vectors of L in WORK(IU)	* singular vectors of L in WORK(IU)
* (Workspace: need M*M+BDSPAC)	* (Workspace: need M*M + BDSPAC)
*	*
CALL DBDSQR( 'U', M, M, M, 0, S, WORK( IE ),	CALL DBDSQR( 'U', M, M, M, 0, S, WORK( IE ),
$ WORK( IU ), LDWRKU, U, LDU, DUM, 1,	$ WORK( IU ), LDWRKU, U, LDU, DUM, 1,
Line 2721	Line 2823
IWORK = ITAU + M	IWORK = ITAU + M
*	*
* Compute A=L*Q, copying result to VT	* Compute A=L*Q, copying result to VT
* (Workspace: need 2M, prefer M+MNB)	* (Workspace: need 2M, prefer M + MNB)
*	*
CALL DGELQF( M, N, A, LDA, WORK( ITAU ),	CALL DGELQF( M, N, A, LDA, WORK( ITAU ),
$ WORK( IWORK ), LWORK-IWORK+1, IERR )	$ WORK( IWORK ), LWORK-IWORK+1, IERR )
CALL DLACPY( 'U', M, N, A, LDA, VT, LDVT )	CALL DLACPY( 'U', M, N, A, LDA, VT, LDVT )
*	*
* Generate Q in VT	* Generate Q in VT
* (Workspace: need 2M, prefer M+MNB)	* (Workspace: need 2M, prefer M + MNB)
*	*
CALL DORGLQ( M, N, M, VT, LDVT, WORK( ITAU ),	CALL DORGLQ( M, N, M, VT, LDVT, WORK( ITAU ),
$ WORK( IWORK ), LWORK-IWORK+1, IERR )	$ WORK( IWORK ), LWORK-IWORK+1, IERR )
Line 2744	Line 2846
IWORK = ITAUP + M	IWORK = ITAUP + M
*	*
* Bidiagonalize L in U	* Bidiagonalize L in U
* (Workspace: need 4M, prefer 3M+2MNB)	* (Workspace: need 4M, prefer 3M + 2MNB)
*	*
CALL DGEBRD( M, M, U, LDU, S, WORK( IE ),	CALL DGEBRD( M, M, U, LDU, S, WORK( IE ),
$ WORK( ITAUQ ), WORK( ITAUP ),	$ WORK( ITAUQ ), WORK( ITAUP ),
Line 2752	Line 2854
*	*
* Multiply right bidiagonalizing vectors in U by Q	* Multiply right bidiagonalizing vectors in U by Q
* in VT	* in VT
* (Workspace: need 3M+N, prefer 3M+N*NB)	* (Workspace: need 3M + N, prefer 3M + N*NB)
*	*
CALL DORMBR( 'P', 'L', 'T', M, N, M, U, LDU,	CALL DORMBR( 'P', 'L', 'T', M, N, M, U, LDU,
$ WORK( ITAUP ), VT, LDVT,	$ WORK( ITAUP ), VT, LDVT,
$ WORK( IWORK ), LWORK-IWORK+1, IERR )	$ WORK( IWORK ), LWORK-IWORK+1, IERR )
*	*
* Generate left bidiagonalizing vectors in U	* Generate left bidiagonalizing vectors in U
* (Workspace: need 4M, prefer 3M+M*NB)	* (Workspace: need 4M, prefer 3M + M*NB)
*	*
CALL DORGBR( 'Q', M, M, M, U, LDU, WORK( ITAUQ ),	CALL DORGBR( 'Q', M, M, M, U, LDU, WORK( ITAUQ ),
$ WORK( IWORK ), LWORK-IWORK+1, IERR )	$ WORK( IWORK ), LWORK-IWORK+1, IERR )
Line 2786	Line 2888
* N right singular vectors to be computed in VT and	* N right singular vectors to be computed in VT and
* no left singular vectors to be computed	* no left singular vectors to be computed
*	*
IF( LWORK.GE.MM+MAX( N+M, 4M, BDSPAC ) ) THEN	IF( LWORK.GE.MM+MAX( N + M, 4M, BDSPAC ) ) THEN
*	*
* Sufficient workspace for a fast algorithm	* Sufficient workspace for a fast algorithm
*	*
Line 2806	Line 2908
IWORK = ITAU + M	IWORK = ITAU + M
*	*
* Compute A=L*Q, copying result to VT	* Compute A=L*Q, copying result to VT
* (Workspace: need MM+2M, prefer MM+M+MNB)	* (Workspace: need MM + 2M, prefer MM + M + MNB)
*	*
CALL DGELQF( M, N, A, LDA, WORK( ITAU ),	CALL DGELQF( M, N, A, LDA, WORK( ITAU ),
$ WORK( IWORK ), LWORK-IWORK+1, IERR )	$ WORK( IWORK ), LWORK-IWORK+1, IERR )
Line 2820	Line 2922
$ WORK( IR+LDWRKR ), LDWRKR )	$ WORK( IR+LDWRKR ), LDWRKR )
*	*
* Generate Q in VT	* Generate Q in VT
* (Workspace: need MM+M+N, prefer MM+M+N*NB)	* (Workspace: need MM + M + N, prefer MM + M + N*NB)
*	*
CALL DORGLQ( N, N, M, VT, LDVT, WORK( ITAU ),	CALL DORGLQ( N, N, M, VT, LDVT, WORK( ITAU ),
$ WORK( IWORK ), LWORK-IWORK+1, IERR )	$ WORK( IWORK ), LWORK-IWORK+1, IERR )
Line 2830	Line 2932
IWORK = ITAUP + M	IWORK = ITAUP + M
*	*
* Bidiagonalize L in WORK(IR)	* Bidiagonalize L in WORK(IR)
* (Workspace: need MM+4M, prefer MM+3M+2MNB)	* (Workspace: need MM + 4M, prefer MM + 3M + 2MNB)
*	*
CALL DGEBRD( M, M, WORK( IR ), LDWRKR, S,	CALL DGEBRD( M, M, WORK( IR ), LDWRKR, S,
$ WORK( IE ), WORK( ITAUQ ),	$ WORK( IE ), WORK( ITAUQ ),
Line 2838	Line 2940
$ LWORK-IWORK+1, IERR )	$ LWORK-IWORK+1, IERR )
*	*
* Generate right bidiagonalizing vectors in WORK(IR)	* Generate right bidiagonalizing vectors in WORK(IR)
* (Workspace: need MM+4M-1,	* (Workspace: need MM + 4M-1,
* prefer MM+3M+(M-1)*NB)	* prefer MM+3M+(M-1)*NB)
*	*
CALL DORGBR( 'P', M, M, M, WORK( IR ), LDWRKR,	CALL DORGBR( 'P', M, M, M, WORK( IR ), LDWRKR,
Line 2848	Line 2950
*	*
* Perform bidiagonal QR iteration, computing right	* Perform bidiagonal QR iteration, computing right
* singular vectors of L in WORK(IR)	* singular vectors of L in WORK(IR)
* (Workspace: need M*M+BDSPAC)	* (Workspace: need M*M + BDSPAC)
*	*
CALL DBDSQR( 'U', M, M, 0, 0, S, WORK( IE ),	CALL DBDSQR( 'U', M, M, 0, 0, S, WORK( IE ),
$ WORK( IR ), LDWRKR, DUM, 1, DUM, 1,	$ WORK( IR ), LDWRKR, DUM, 1, DUM, 1,
Line 2873	Line 2975
IWORK = ITAU + M	IWORK = ITAU + M
*	*
* Compute A=L*Q, copying result to VT	* Compute A=L*Q, copying result to VT
* (Workspace: need 2M, prefer M+MNB)	* (Workspace: need 2M, prefer M + MNB)
*	*
CALL DGELQF( M, N, A, LDA, WORK( ITAU ),	CALL DGELQF( M, N, A, LDA, WORK( ITAU ),
$ WORK( IWORK ), LWORK-IWORK+1, IERR )	$ WORK( IWORK ), LWORK-IWORK+1, IERR )
CALL DLACPY( 'U', M, N, A, LDA, VT, LDVT )	CALL DLACPY( 'U', M, N, A, LDA, VT, LDVT )
*	*
* Generate Q in VT	* Generate Q in VT
* (Workspace: need M+N, prefer M+N*NB)	* (Workspace: need M + N, prefer M + N*NB)
*	*
CALL DORGLQ( N, N, M, VT, LDVT, WORK( ITAU ),	CALL DORGLQ( N, N, M, VT, LDVT, WORK( ITAU ),
$ WORK( IWORK ), LWORK-IWORK+1, IERR )	$ WORK( IWORK ), LWORK-IWORK+1, IERR )
Line 2895	Line 2997
$ LDA )	$ LDA )
*	*
* Bidiagonalize L in A	* Bidiagonalize L in A
* (Workspace: need 4M, prefer 3M+2MNB)	* (Workspace: need 4M, prefer 3M + 2MNB)
*	*
CALL DGEBRD( M, M, A, LDA, S, WORK( IE ),	CALL DGEBRD( M, M, A, LDA, S, WORK( IE ),
$ WORK( ITAUQ ), WORK( ITAUP ),	$ WORK( ITAUQ ), WORK( ITAUP ),
Line 2903	Line 3005
*	*
* Multiply right bidiagonalizing vectors in A by Q	* Multiply right bidiagonalizing vectors in A by Q
* in VT	* in VT
* (Workspace: need 3M+N, prefer 3M+N*NB)	* (Workspace: need 3M + N, prefer 3M + N*NB)
*	*
CALL DORMBR( 'P', 'L', 'T', M, N, M, A, LDA,	CALL DORMBR( 'P', 'L', 'T', M, N, M, A, LDA,
$ WORK( ITAUP ), VT, LDVT,	$ WORK( ITAUP ), VT, LDVT,
Line 2926	Line 3028
* 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 overwritten on A	* M left singular vectors to be overwritten on A
*	*
IF( LWORK.GE.2MM+MAX( N+M, 4*M, BDSPAC ) ) THEN	IF( LWORK.GE.2MM+MAX( N + M, 4*M, BDSPAC ) ) THEN
*	*
* Sufficient workspace for a fast algorithm	* Sufficient workspace for a fast algorithm
*	*
Line 2938	Line 3040
LDWRKU = LDA	LDWRKU = LDA
IR = IU + LDWRKU*M	IR = IU + LDWRKU*M
LDWRKR = LDA	LDWRKR = LDA
ELSE IF( LWORK.GE.WRKBL+( LDA+M )*M ) THEN	ELSE IF( LWORK.GE.WRKBL+( LDA + M )*M ) THEN
*	*
* WORK(IU) is LDA by M and WORK(IR) is M by M	* WORK(IU) is LDA by M and WORK(IR) is M by M
*	*
Line 2957	Line 3059
IWORK = ITAU + M	IWORK = ITAU + M
*	*
* Compute A=L*Q, copying result to VT	* Compute A=L*Q, copying result to VT
* (Workspace: need 2MM+2M, prefer 2MM+M+MNB)	* (Workspace: need 2MM + 2M, prefer 2MM + M + MNB)
*	*
CALL DGELQF( M, N, A, LDA, WORK( ITAU ),	CALL DGELQF( M, N, A, LDA, WORK( ITAU ),
$ WORK( IWORK ), LWORK-IWORK+1, IERR )	$ WORK( IWORK ), LWORK-IWORK+1, IERR )
CALL DLACPY( 'U', M, N, A, LDA, VT, LDVT )	CALL DLACPY( 'U', M, N, A, LDA, VT, LDVT )
*	*
* Generate Q in VT	* Generate Q in VT
* (Workspace: need 2MM+M+N, prefer 2MM+M+N*NB)	* (Workspace: need 2MM + M + N, prefer 2MM + M + N*NB)
*	*
CALL DORGLQ( N, N, M, VT, LDVT, WORK( ITAU ),	CALL DORGLQ( N, N, M, VT, LDVT, WORK( ITAU ),
$ WORK( IWORK ), LWORK-IWORK+1, IERR )	$ WORK( IWORK ), LWORK-IWORK+1, IERR )
Line 2982	Line 3084
*	*
* Bidiagonalize L in WORK(IU), copying result to	* Bidiagonalize L in WORK(IU), copying result to
* WORK(IR)	* WORK(IR)
* (Workspace: need 2MM+4*M,	* (Workspace: need 2MM + 4*M,
* prefer 2MM+3M+2M*NB)	* prefer 2MM+3M+2M*NB)
*	*
CALL DGEBRD( M, M, WORK( IU ), LDWRKU, S,	CALL DGEBRD( M, M, WORK( IU ), LDWRKU, S,
Line 2993	Line 3095
$ WORK( IR ), LDWRKR )	$ WORK( IR ), LDWRKR )
*	*
* Generate right bidiagonalizing vectors in WORK(IU)	* Generate right bidiagonalizing vectors in WORK(IU)
* (Workspace: need 2MM+4*M-1,	* (Workspace: need 2MM + 4*M-1,
* prefer 2MM+3M+(M-1)NB)	* prefer 2MM+3M+(M-1)NB)
*	*
CALL DORGBR( 'P', M, M, M, WORK( IU ), LDWRKU,	CALL DORGBR( 'P', M, M, M, WORK( IU ), LDWRKU,
Line 3001	Line 3103
$ LWORK-IWORK+1, IERR )	$ LWORK-IWORK+1, IERR )
*	*
* Generate left bidiagonalizing vectors in WORK(IR)	* Generate left bidiagonalizing vectors in WORK(IR)
* (Workspace: need 2MM+4M, prefer 2MM+3M+M*NB)	* (Workspace: need 2MM + 4M, prefer 2MM + 3M + M*NB)
*	*
CALL DORGBR( 'Q', M, M, M, WORK( IR ), LDWRKR,	CALL DORGBR( 'Q', M, M, M, WORK( IR ), LDWRKR,
$ WORK( ITAUQ ), WORK( IWORK ),	$ WORK( ITAUQ ), WORK( IWORK ),
Line 3011	Line 3113
* Perform bidiagonal QR iteration, computing left	* Perform bidiagonal QR iteration, computing left
* singular vectors of L in WORK(IR) and computing	* singular vectors of L in WORK(IR) and computing
* right singular vectors of L in WORK(IU)	* right singular vectors of L in WORK(IU)
* (Workspace: need 2MM+BDSPAC)	* (Workspace: need 2MM + BDSPAC)
*	*
CALL DBDSQR( 'U', M, M, M, 0, S, WORK( IE ),	CALL DBDSQR( 'U', M, M, M, 0, S, WORK( IE ),
$ WORK( IU ), LDWRKU, WORK( IR ),	$ WORK( IU ), LDWRKU, WORK( IR ),
Line 3041	Line 3143
IWORK = ITAU + M	IWORK = ITAU + M
*	*
* Compute A=L*Q, copying result to VT	* Compute A=L*Q, copying result to VT
* (Workspace: need 2M, prefer M+MNB)	* (Workspace: need 2M, prefer M + MNB)
*	*
CALL DGELQF( M, N, A, LDA, WORK( ITAU ),	CALL DGELQF( M, N, A, LDA, WORK( ITAU ),
$ WORK( IWORK ), LWORK-IWORK+1, IERR )	$ WORK( IWORK ), LWORK-IWORK+1, IERR )
CALL DLACPY( 'U', M, N, A, LDA, VT, LDVT )	CALL DLACPY( 'U', M, N, A, LDA, VT, LDVT )
*	*
* Generate Q in VT	* Generate Q in VT
* (Workspace: need M+N, prefer M+N*NB)	* (Workspace: need M + N, prefer M + N*NB)
*	*
CALL DORGLQ( N, N, M, VT, LDVT, WORK( ITAU ),	CALL DORGLQ( N, N, M, VT, LDVT, WORK( ITAU ),
$ WORK( IWORK ), LWORK-IWORK+1, IERR )	$ WORK( IWORK ), LWORK-IWORK+1, IERR )
Line 3063	Line 3165
$ LDA )	$ LDA )
*	*
* Bidiagonalize L in A	* Bidiagonalize L in A
* (Workspace: need 4M, prefer 3M+2MNB)	* (Workspace: need 4M, prefer 3M + 2MNB)
*	*
CALL DGEBRD( M, M, A, LDA, S, WORK( IE ),	CALL DGEBRD( M, M, A, LDA, S, WORK( IE ),
$ WORK( ITAUQ ), WORK( ITAUP ),	$ WORK( ITAUQ ), WORK( ITAUP ),
Line 3071	Line 3173
*	*
* Multiply right bidiagonalizing vectors in A by Q	* Multiply right bidiagonalizing vectors in A by Q
* in VT	* in VT
* (Workspace: need 3M+N, prefer 3M+N*NB)	* (Workspace: need 3M + N, prefer 3M + N*NB)
*	*
CALL DORMBR( 'P', 'L', 'T', M, N, M, A, LDA,	CALL DORMBR( 'P', 'L', 'T', M, N, M, A, LDA,
$ WORK( ITAUP ), VT, LDVT,	$ WORK( ITAUP ), VT, LDVT,
$ WORK( IWORK ), LWORK-IWORK+1, IERR )	$ WORK( IWORK ), LWORK-IWORK+1, IERR )
*	*
* Generate left bidiagonalizing vectors in A	* Generate left bidiagonalizing vectors in A
* (Workspace: need 4M, prefer 3M+M*NB)	* (Workspace: need 4M, prefer 3M + M*NB)
*	*
CALL DORGBR( 'Q', M, M, M, A, LDA, WORK( ITAUQ ),	CALL DORGBR( 'Q', M, M, M, A, LDA, WORK( ITAUQ ),
$ WORK( IWORK ), LWORK-IWORK+1, IERR )	$ WORK( IWORK ), LWORK-IWORK+1, IERR )
Line 3102	Line 3204
* 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
*	*
IF( LWORK.GE.MM+MAX( N+M, 4M, BDSPAC ) ) THEN	IF( LWORK.GE.MM+MAX( N + M, 4M, BDSPAC ) ) THEN
*	*
* Sufficient workspace for a fast algorithm	* Sufficient workspace for a fast algorithm
*	*
Line 3122	Line 3224
IWORK = ITAU + M	IWORK = ITAU + M
*	*
* Compute A=L*Q, copying result to VT	* Compute A=L*Q, copying result to VT
* (Workspace: need MM+2M, prefer MM+M+MNB)	* (Workspace: need MM + 2M, prefer MM + M + MNB)
*	*
CALL DGELQF( M, N, A, LDA, WORK( ITAU ),	CALL DGELQF( M, N, A, LDA, WORK( ITAU ),
$ WORK( IWORK ), LWORK-IWORK+1, IERR )	$ WORK( IWORK ), LWORK-IWORK+1, IERR )
CALL DLACPY( 'U', M, N, A, LDA, VT, LDVT )	CALL DLACPY( 'U', M, N, A, LDA, VT, LDVT )
*	*
* Generate Q in VT	* Generate Q in VT
* (Workspace: need MM+M+N, prefer MM+M+N*NB)	* (Workspace: need MM + M + N, prefer MM + M + N*NB)
*	*
CALL DORGLQ( N, N, M, VT, LDVT, WORK( ITAU ),	CALL DORGLQ( N, N, M, VT, LDVT, WORK( ITAU ),
$ WORK( IWORK ), LWORK-IWORK+1, IERR )	$ WORK( IWORK ), LWORK-IWORK+1, IERR )
Line 3146	Line 3248
IWORK = ITAUP + M	IWORK = ITAUP + M
*	*
* Bidiagonalize L in WORK(IU), copying result to U	* Bidiagonalize L in WORK(IU), copying result to U
* (Workspace: need MM+4M, prefer MM+3M+2MNB)	* (Workspace: need MM + 4M, prefer MM + 3M + 2MNB)
*	*
CALL DGEBRD( M, M, WORK( IU ), LDWRKU, S,	CALL DGEBRD( M, M, WORK( IU ), LDWRKU, S,
$ WORK( IE ), WORK( ITAUQ ),	$ WORK( IE ), WORK( ITAUQ ),
Line 3156	Line 3258
$ LDU )	$ LDU )
*	*
* Generate right bidiagonalizing vectors in WORK(IU)	* Generate right bidiagonalizing vectors in WORK(IU)
* (Workspace: need MM+4M, prefer MM+3M+(M-1)*NB)	* (Workspace: need MM + 4M, prefer MM + 3M + (M-1)*NB)
*	*
CALL DORGBR( 'P', M, M, M, WORK( IU ), LDWRKU,	CALL DORGBR( 'P', M, M, M, WORK( IU ), LDWRKU,
$ WORK( ITAUP ), WORK( IWORK ),	$ WORK( ITAUP ), WORK( IWORK ),
$ LWORK-IWORK+1, IERR )	$ LWORK-IWORK+1, IERR )
*	*
* Generate left bidiagonalizing vectors in U	* Generate left bidiagonalizing vectors in U
* (Workspace: need MM+4M, prefer MM+3M+M*NB)	* (Workspace: need MM + 4M, prefer MM + 3M + M*NB)
*	*
CALL DORGBR( 'Q', M, M, M, U, LDU, WORK( ITAUQ ),	CALL DORGBR( 'Q', M, M, M, U, LDU, WORK( ITAUQ ),
$ WORK( IWORK ), LWORK-IWORK+1, IERR )	$ WORK( IWORK ), LWORK-IWORK+1, IERR )
Line 3172	Line 3274
* Perform bidiagonal QR iteration, computing left	* Perform bidiagonal QR iteration, computing left
* singular vectors of L in U and computing right	* singular vectors of L in U and computing right
* singular vectors of L in WORK(IU)	* singular vectors of L in WORK(IU)
* (Workspace: need M*M+BDSPAC)	* (Workspace: need M*M + BDSPAC)
*	*
CALL DBDSQR( 'U', M, M, M, 0, S, WORK( IE ),	CALL DBDSQR( 'U', M, M, M, 0, S, WORK( IE ),
$ WORK( IU ), LDWRKU, U, LDU, DUM, 1,	$ WORK( IU ), LDWRKU, U, LDU, DUM, 1,
Line 3197	Line 3299
IWORK = ITAU + M	IWORK = ITAU + M
*	*
* Compute A=L*Q, copying result to VT	* Compute A=L*Q, copying result to VT
* (Workspace: need 2M, prefer M+MNB)	* (Workspace: need 2M, prefer M + MNB)
*	*
CALL DGELQF( M, N, A, LDA, WORK( ITAU ),	CALL DGELQF( M, N, A, LDA, WORK( ITAU ),
$ WORK( IWORK ), LWORK-IWORK+1, IERR )	$ WORK( IWORK ), LWORK-IWORK+1, IERR )
CALL DLACPY( 'U', M, N, A, LDA, VT, LDVT )	CALL DLACPY( 'U', M, N, A, LDA, VT, LDVT )
*	*
* Generate Q in VT	* Generate Q in VT
* (Workspace: need M+N, prefer M+N*NB)	* (Workspace: need M + N, prefer M + N*NB)
*	*
CALL DORGLQ( N, N, M, VT, LDVT, WORK( ITAU ),	CALL DORGLQ( N, N, M, VT, LDVT, WORK( ITAU ),
$ WORK( IWORK ), LWORK-IWORK+1, IERR )	$ WORK( IWORK ), LWORK-IWORK+1, IERR )
Line 3220	Line 3322
IWORK = ITAUP + M	IWORK = ITAUP + M
*	*
* Bidiagonalize L in U	* Bidiagonalize L in U
* (Workspace: need 4M, prefer 3M+2MNB)	* (Workspace: need 4M, prefer 3M + 2MNB)
*	*
CALL DGEBRD( M, M, U, LDU, S, WORK( IE ),	CALL DGEBRD( M, M, U, LDU, S, WORK( IE ),
$ WORK( ITAUQ ), WORK( ITAUP ),	$ WORK( ITAUQ ), WORK( ITAUP ),
Line 3228	Line 3330
*	*
* Multiply right bidiagonalizing vectors in U by Q	* Multiply right bidiagonalizing vectors in U by Q
* in VT	* in VT
* (Workspace: need 3M+N, prefer 3M+N*NB)	* (Workspace: need 3M + N, prefer 3M + N*NB)
*	*
CALL DORMBR( 'P', 'L', 'T', M, N, M, U, LDU,	CALL DORMBR( 'P', 'L', 'T', M, N, M, U, LDU,
$ WORK( ITAUP ), VT, LDVT,	$ WORK( ITAUP ), VT, LDVT,
$ WORK( IWORK ), LWORK-IWORK+1, IERR )	$ WORK( IWORK ), LWORK-IWORK+1, IERR )
*	*
* Generate left bidiagonalizing vectors in U	* Generate left bidiagonalizing vectors in U
* (Workspace: need 4M, prefer 3M+M*NB)	* (Workspace: need 4M, prefer 3M + M*NB)
*	*
CALL DORGBR( 'Q', M, M, M, U, LDU, WORK( ITAUQ ),	CALL DORGBR( 'Q', M, M, M, U, LDU, WORK( ITAUQ ),
$ WORK( IWORK ), LWORK-IWORK+1, IERR )	$ WORK( IWORK ), LWORK-IWORK+1, IERR )
Line 3269	Line 3371
IWORK = ITAUP + M	IWORK = ITAUP + M
*	*
* Bidiagonalize A	* Bidiagonalize A
* (Workspace: need 3M+N, prefer 3M+(M+N)*NB)	* (Workspace: need 3M + N, prefer 3M + (M + N)*NB)
*	*
CALL DGEBRD( M, N, A, LDA, S, WORK( IE ), WORK( ITAUQ ),	CALL DGEBRD( M, N, A, LDA, S, WORK( IE ), WORK( ITAUQ ),
$ WORK( ITAUP ), WORK( IWORK ), LWORK-IWORK+1,	$ WORK( ITAUP ), WORK( IWORK ), LWORK-IWORK+1,
Line 3278	Line 3380
*	*
* If left singular vectors desired in U, copy result to U	* If left singular vectors desired in U, copy result to U
* and generate left bidiagonalizing vectors in U	* and generate left bidiagonalizing vectors in U
* (Workspace: need 4M-1, prefer 3M+(M-1)*NB)	* (Workspace: need 4M-1, prefer 3M + (M-1)*NB)
*	*
CALL DLACPY( 'L', M, M, A, LDA, U, LDU )	CALL DLACPY( 'L', M, M, A, LDA, U, LDU )
CALL DORGBR( 'Q', M, M, N, U, LDU, WORK( ITAUQ ),	CALL DORGBR( 'Q', M, M, N, U, LDU, WORK( ITAUQ ),
Line 3288	Line 3390
*	*
* If right singular vectors desired in VT, copy result to	* If right singular vectors desired in VT, copy result to
* VT and generate right bidiagonalizing vectors in VT	* VT and generate right bidiagonalizing vectors in VT
* (Workspace: need 3M+NRVT, prefer 3M+NRVT*NB)	* (Workspace: need 3M + NRVT, prefer 3M + NRVT*NB)
*	*
CALL DLACPY( 'U', M, N, A, LDA, VT, LDVT )	CALL DLACPY( 'U', M, N, A, LDA, VT, LDVT )
IF( WNTVA )	IF( WNTVA )
Line 3302	Line 3404
*	*
* If left singular vectors desired in A, generate left	* If left singular vectors desired in A, generate left
* bidiagonalizing vectors in A	* bidiagonalizing vectors in A
* (Workspace: need 4M-1, prefer 3M+(M-1)*NB)	* (Workspace: need 4M-1, prefer 3M + (M-1)*NB)
*	*
CALL DORGBR( 'Q', M, M, N, A, LDA, WORK( ITAUQ ),	CALL DORGBR( 'Q', M, M, N, A, LDA, WORK( ITAUQ ),
$ WORK( IWORK ), LWORK-IWORK+1, IERR )	$ WORK( IWORK ), LWORK-IWORK+1, IERR )
Line 3311	Line 3413
*	*
* If right singular vectors desired in A, generate right	* If right singular vectors desired in A, generate right
* bidiagonalizing vectors in A	* bidiagonalizing vectors in A
* (Workspace: need 4M, prefer 3M+M*NB)	* (Workspace: need 4M, prefer 3M + M*NB)
*	*
CALL DORGBR( 'P', M, N, M, A, LDA, WORK( ITAUP ),	CALL DORGBR( 'P', M, N, M, A, LDA, WORK( ITAUP ),
$ WORK( IWORK ), LWORK-IWORK+1, IERR )	$ WORK( IWORK ), LWORK-IWORK+1, IERR )

CVSweb interface <joel.bertrand@systella.fr>

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