rpl/lapack/lapack/dgesdd.f - diff

Return to dgesdd.f CVS log

Up to [local] / rpl / lapack / lapack

Diff for /rpl/lapack/lapack/dgesdd.f between versions 1.15 and 1.16

version 1.15, 2015/11/26 11:44:15	version 1.16, 2016/08/27 15:27:08
Line 18	Line 18
* Definition:	* Definition:
* ===========	* ===========
*	*
* SUBROUTINE DGESDD( JOBZ, M, N, A, LDA, S, U, LDU, VT, LDVT, WORK,	* SUBROUTINE DGESDD( JOBZ, M, N, A, LDA, S, U, LDU, VT, LDVT,
* LWORK, IWORK, INFO )	* WORK, LWORK, IWORK, INFO )
*	*
* .. Scalar Arguments ..	* .. Scalar Arguments ..
* CHARACTER JOBZ	* CHARACTER JOBZ
Line 154	Line 154
*> \param[in] LDVT	*> \param[in] LDVT
*> \verbatim	*> \verbatim
*> LDVT is INTEGER	*> LDVT is INTEGER
*> The leading dimension of the array VT. LDVT >= 1; if	*> The leading dimension of the array VT. LDVT >= 1;
*> JOBZ = 'A' or JOBZ = 'O' and M >= N, LDVT >= N;	*> if JOBZ = 'A' or JOBZ = 'O' and M >= N, LDVT >= N;
*> if JOBZ = 'S', LDVT >= min(M,N).	*> if JOBZ = 'S', LDVT >= min(M,N).
*> \endverbatim	*> \endverbatim
*>	*>
Line 169	Line 169
*> \verbatim	*> \verbatim
*> LWORK is INTEGER	*> LWORK is INTEGER
*> The dimension of the array WORK. LWORK >= 1.	*> The dimension of the array WORK. LWORK >= 1.
*> If JOBZ = 'N',	*> If LWORK = -1, a workspace query is assumed. The optimal
> LWORK >= 3min(M,N) + max(max(M,N),7*min(M,N)).	*> size for the WORK array is calculated and stored in WORK(1),
*> If JOBZ = 'O',	*> and no other work except argument checking is performed.
> LWORK >= 3min(M,N) +	*>
> max(max(M,N),5min(M,N)min(M,N)+4min(M,N)).	*> Let mx = max(M,N) and mn = min(M,N).
*> If JOBZ = 'S' or 'A'	> If JOBZ = 'N', LWORK >= 3mn + max( mx, 7*mn ).
> LWORK >= min(M,N)(7+4*min(M,N))	> If JOBZ = 'O', LWORK >= 3mn + max( mx, 5mnmn + 4*mn ).
*> For good performance, LWORK should generally be larger.	> If JOBZ = 'S', LWORK >= 4mnmn + 7mn.
*> If LWORK = -1 but other input arguments are legal, WORK(1)	> If JOBZ = 'A', LWORK >= 4mnmn + 6mn + mx.
*> returns the optimal LWORK.	*> These are not tight minimums in all cases; see comments inside code.
	*> For good performance, LWORK should generally be larger;
	*> a query is recommended.
*> \endverbatim	*> \endverbatim
*>	*>
*> \param[out] IWORK	*> \param[out] IWORK
Line 202	Line 204
*> \author Univ. of Colorado Denver	*> \author Univ. of Colorado Denver
*> \author NAG Ltd.	*> \author NAG Ltd.
*	*
*> \date November 2015	*> \date June 2016
*	*
*> \ingroup doubleGEsing	*> \ingroup doubleGEsing
*	*
Line 212	Line 214
*> Ming Gu and Huan Ren, Computer Science Division, University of	*> Ming Gu and Huan Ren, Computer Science Division, University of
*> California at Berkeley, USA	*> California at Berkeley, USA
*>	*>
	*> @precisions fortran d -> s
* =====================================================================	* =====================================================================
SUBROUTINE DGESDD( JOBZ, M, N, A, LDA, S, U, LDU, VT, LDVT, WORK,	SUBROUTINE DGESDD( JOBZ, M, N, A, LDA, S, U, LDU, VT, LDVT,
$ LWORK, IWORK, INFO )	$ WORK, LWORK, IWORK, INFO )
	implicit none
*	*
* -- LAPACK driver routine (version 3.6.0) --	* -- LAPACK driver routine (version 3.6.1) --
* -- 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 2015	* June 2016
*	*
* .. Scalar Arguments ..	* .. Scalar Arguments ..
CHARACTER JOBZ	CHARACTER JOBZ
Line 243	Line 247
$ IR, ISCL, ITAU, ITAUP, ITAUQ, IU, IVT, LDWKVT,	$ IR, ISCL, ITAU, ITAUP, ITAUQ, IU, IVT, LDWKVT,
$ LDWRKL, LDWRKR, LDWRKU, MAXWRK, MINMN, MINWRK,	$ LDWRKL, LDWRKR, LDWRKU, MAXWRK, MINMN, MINWRK,
$ MNTHR, NWORK, WRKBL	$ MNTHR, NWORK, WRKBL
	INTEGER LWORK_DGEBRD_MN, LWORK_DGEBRD_MM,
	$ LWORK_DGEBRD_NN, LWORK_DGELQF_MN,
	$ LWORK_DGEQRF_MN,
	$ LWORK_DORGBR_P_MM, LWORK_DORGBR_Q_NN,
	$ LWORK_DORGLQ_MN, LWORK_DORGLQ_NN,
	$ LWORK_DORGQR_MM, LWORK_DORGQR_MN,
	$ LWORK_DORMBR_PRT_MM, LWORK_DORMBR_QLN_MM,
	$ LWORK_DORMBR_PRT_MN, LWORK_DORMBR_QLN_MN,
	$ LWORK_DORMBR_PRT_NN, LWORK_DORMBR_QLN_NN
DOUBLE PRECISION ANRM, BIGNUM, EPS, SMLNUM	DOUBLE PRECISION ANRM, BIGNUM, EPS, SMLNUM
* ..	* ..
* .. Local Arrays ..	* .. Local Arrays ..
Line 256	Line 269
* ..	* ..
* .. External Functions ..	* .. External Functions ..
LOGICAL LSAME	LOGICAL LSAME
INTEGER ILAENV
DOUBLE PRECISION DLAMCH, DLANGE	DOUBLE PRECISION DLAMCH, DLANGE
EXTERNAL DLAMCH, DLANGE, ILAENV, LSAME	EXTERNAL DLAMCH, DLANGE, LSAME
* ..	* ..
* .. Intrinsic Functions ..	* .. Intrinsic Functions ..
INTRINSIC INT, MAX, MIN, SQRT	INTRINSIC INT, MAX, MIN, SQRT
Line 267	Line 279
*	*
* Test the input arguments	* Test the input arguments
*	*
INFO = 0	INFO = 0
MINMN = MIN( M, N )	MINMN = MIN( M, N )
WNTQA = LSAME( JOBZ, 'A' )	WNTQA = LSAME( JOBZ, 'A' )
WNTQS = LSAME( JOBZ, 'S' )	WNTQS = LSAME( JOBZ, 'S' )
WNTQAS = WNTQA .OR. WNTQS	WNTQAS = WNTQA .OR. WNTQS
WNTQO = LSAME( JOBZ, 'O' )	WNTQO = LSAME( JOBZ, 'O' )
WNTQN = LSAME( JOBZ, 'N' )	WNTQN = LSAME( JOBZ, 'N' )
LQUERY = ( LWORK.EQ.-1 )	LQUERY = ( LWORK.EQ.-1 )
*	*
IF( .NOT.( WNTQA .OR. WNTQS .OR. WNTQO .OR. WNTQN ) ) THEN	IF( .NOT.( WNTQA .OR. WNTQS .OR. WNTQO .OR. WNTQN ) ) THEN
Line 294	Line 306
END IF	END IF
*	*
* Compute workspace	* Compute workspace
* (Note: Comments in the code beginning "Workspace:" describe the	* Note: Comments in the code beginning "Workspace:" describe the
* minimal amount of workspace needed at that point in the code,	* minimal amount of workspace allocated at that point in the code,
* as well as the preferred amount for good performance.	* as well as the preferred amount for good performance.
* NB refers to the optimal block size for the immediately	* NB refers to the optimal block size for the immediately
* following subroutine, as returned by ILAENV.)	* following subroutine, as returned by ILAENV.
*	*
IF( INFO.EQ.0 ) THEN	IF( INFO.EQ.0 ) THEN
MINWRK = 1	MINWRK = 1
MAXWRK = 1	MAXWRK = 1
	BDSPAC = 0
	MNTHR = INT( MINMN*11.0D0 / 6.0D0 )
IF( M.GE.N .AND. MINMN.GT.0 ) THEN	IF( M.GE.N .AND. MINMN.GT.0 ) THEN
*	*
* Compute space needed for DBDSDC	* Compute space needed for DBDSDC
*	*
MNTHR = INT( MINMN*11.0D0 / 6.0D0 )
IF( WNTQN ) THEN	IF( WNTQN ) THEN
	* dbdsdc needs only 4N (or 6N for uplo=L for LAPACK <= 3.6)
	* keep 7*N for backwards compatability.
BDSPAC = 7*N	BDSPAC = 7*N
ELSE	ELSE
BDSPAC = 3NN + 4*N	BDSPAC = 3NN + 4*N
END IF	END IF
	*
	* Compute space preferred for each routine
	CALL DGEBRD( M, N, DUM(1), M, DUM(1), DUM(1), DUM(1),
	$ DUM(1), DUM(1), -1, IERR )
	LWORK_DGEBRD_MN = INT( DUM(1) )
	*
	CALL DGEBRD( N, N, DUM(1), N, DUM(1), DUM(1), DUM(1),
	$ DUM(1), DUM(1), -1, IERR )
	LWORK_DGEBRD_NN = INT( DUM(1) )
	*
	CALL DGEQRF( M, N, DUM(1), M, DUM(1), DUM(1), -1, IERR )
	LWORK_DGEQRF_MN = INT( DUM(1) )
	*
	CALL DORGBR( 'Q', N, N, N, DUM(1), N, DUM(1), DUM(1), -1,
	$ IERR )
	LWORK_DORGBR_Q_NN = INT( DUM(1) )
	*
	CALL DORGQR( M, M, N, DUM(1), M, DUM(1), DUM(1), -1, IERR )
	LWORK_DORGQR_MM = INT( DUM(1) )
	*
	CALL DORGQR( M, N, N, DUM(1), M, DUM(1), DUM(1), -1, IERR )
	LWORK_DORGQR_MN = INT( DUM(1) )
	*
	CALL DORMBR( 'P', 'R', 'T', N, N, N, DUM(1), N,
	$ DUM(1), DUM(1), N, DUM(1), -1, IERR )
	LWORK_DORMBR_PRT_NN = INT( DUM(1) )
	*
	CALL DORMBR( 'Q', 'L', 'N', N, N, N, DUM(1), N,
	$ DUM(1), DUM(1), N, DUM(1), -1, IERR )
	LWORK_DORMBR_QLN_NN = INT( DUM(1) )
	*
	CALL DORMBR( 'Q', 'L', 'N', M, N, N, DUM(1), M,
	$ DUM(1), DUM(1), M, DUM(1), -1, IERR )
	LWORK_DORMBR_QLN_MN = INT( DUM(1) )
	*
	CALL DORMBR( 'Q', 'L', 'N', M, M, N, DUM(1), M,
	$ DUM(1), DUM(1), M, DUM(1), -1, IERR )
	LWORK_DORMBR_QLN_MM = INT( DUM(1) )
	*
IF( M.GE.MNTHR ) THEN	IF( M.GE.MNTHR ) THEN
IF( WNTQN ) THEN	IF( WNTQN ) THEN
*	*
* Path 1 (M much larger than N, JOBZ='N')	* Path 1 (M >> N, JOBZ='N')
*	*
WRKBL = N + N*ILAENV( 1, 'DGEQRF', ' ', M, N, -1,	WRKBL = N + LWORK_DGEQRF_MN
$ -1 )	WRKBL = MAX( WRKBL, 3*N + LWORK_DGEBRD_NN )
WRKBL = MAX( WRKBL, 3N+2N*	MAXWRK = MAX( WRKBL, BDSPAC + N )
$ ILAENV( 1, 'DGEBRD', ' ', N, N, -1, -1 ) )
MAXWRK = MAX( WRKBL, BDSPAC+N )
MINWRK = BDSPAC + N	MINWRK = BDSPAC + N
ELSE IF( WNTQO ) THEN	ELSE IF( WNTQO ) THEN
*	*
* Path 2 (M much larger than N, JOBZ='O')	* Path 2 (M >> N, JOBZ='O')
*	*
WRKBL = N + N*ILAENV( 1, 'DGEQRF', ' ', M, N, -1, -1 )	WRKBL = N + LWORK_DGEQRF_MN
WRKBL = MAX( WRKBL, N+N*ILAENV( 1, 'DORGQR', ' ', M,	WRKBL = MAX( WRKBL, N + LWORK_DORGQR_MN )
$ N, N, -1 ) )	WRKBL = MAX( WRKBL, 3*N + LWORK_DGEBRD_NN )
WRKBL = MAX( WRKBL, 3N+2N*	WRKBL = MAX( WRKBL, 3*N + LWORK_DORMBR_QLN_NN )
$ ILAENV( 1, 'DGEBRD', ' ', N, N, -1, -1 ) )	WRKBL = MAX( WRKBL, 3*N + LWORK_DORMBR_PRT_NN )
WRKBL = MAX( WRKBL, 3N+N	WRKBL = MAX( WRKBL, 3*N + BDSPAC )
$ ILAENV( 1, 'DORMBR', 'QLN', N, N, N, -1 ) )
WRKBL = MAX( WRKBL, 3N+N
$ ILAENV( 1, 'DORMBR', 'PRT', N, N, N, -1 ) )
WRKBL = MAX( WRKBL, BDSPAC+3*N )
MAXWRK = WRKBL + 2NN	MAXWRK = WRKBL + 2NN
MINWRK = BDSPAC + 2NN + 3*N	MINWRK = BDSPAC + 2NN + 3*N
ELSE IF( WNTQS ) THEN	ELSE IF( WNTQS ) THEN
*	*
* Path 3 (M much larger than N, JOBZ='S')	* Path 3 (M >> N, JOBZ='S')
*	*
WRKBL = N + N*ILAENV( 1, 'DGEQRF', ' ', M, N, -1, -1 )	WRKBL = N + LWORK_DGEQRF_MN
WRKBL = MAX( WRKBL, N+N*ILAENV( 1, 'DORGQR', ' ', M,	WRKBL = MAX( WRKBL, N + LWORK_DORGQR_MN )
$ N, N, -1 ) )	WRKBL = MAX( WRKBL, 3*N + LWORK_DGEBRD_NN )
WRKBL = MAX( WRKBL, 3N+2N*	WRKBL = MAX( WRKBL, 3*N + LWORK_DORMBR_QLN_NN )
$ ILAENV( 1, 'DGEBRD', ' ', N, N, -1, -1 ) )	WRKBL = MAX( WRKBL, 3*N + LWORK_DORMBR_PRT_NN )
WRKBL = MAX( WRKBL, 3N+N	WRKBL = MAX( WRKBL, 3*N + BDSPAC )
$ ILAENV( 1, 'DORMBR', 'QLN', N, N, N, -1 ) )
WRKBL = MAX( WRKBL, 3N+N
$ ILAENV( 1, 'DORMBR', 'PRT', N, N, N, -1 ) )
WRKBL = MAX( WRKBL, BDSPAC+3*N )
MAXWRK = WRKBL + N*N	MAXWRK = WRKBL + N*N
MINWRK = BDSPAC + NN + 3N	MINWRK = BDSPAC + NN + 3N
ELSE IF( WNTQA ) THEN	ELSE IF( WNTQA ) THEN
*	*
* Path 4 (M much larger than N, JOBZ='A')	* Path 4 (M >> N, JOBZ='A')
*	*
WRKBL = N + N*ILAENV( 1, 'DGEQRF', ' ', M, N, -1, -1 )	WRKBL = N + LWORK_DGEQRF_MN
WRKBL = MAX( WRKBL, N+M*ILAENV( 1, 'DORGQR', ' ', M,	WRKBL = MAX( WRKBL, N + LWORK_DORGQR_MM )
$ M, N, -1 ) )	WRKBL = MAX( WRKBL, 3*N + LWORK_DGEBRD_NN )
WRKBL = MAX( WRKBL, 3N+2N*	WRKBL = MAX( WRKBL, 3*N + LWORK_DORMBR_QLN_NN )
$ ILAENV( 1, 'DGEBRD', ' ', N, N, -1, -1 ) )	WRKBL = MAX( WRKBL, 3*N + LWORK_DORMBR_PRT_NN )
WRKBL = MAX( WRKBL, 3N+N	WRKBL = MAX( WRKBL, 3*N + BDSPAC )
$ ILAENV( 1, 'DORMBR', 'QLN', N, N, N, -1 ) )
WRKBL = MAX( WRKBL, 3N+N
$ ILAENV( 1, 'DORMBR', 'PRT', N, N, N, -1 ) )
WRKBL = MAX( WRKBL, BDSPAC+3*N )
MAXWRK = WRKBL + N*N	MAXWRK = WRKBL + N*N
MINWRK = BDSPAC + NN + 2N + M	MINWRK = NN + MAX( 3N + BDSPAC, N + M )
END IF	END IF
ELSE	ELSE
*	*
* Path 5 (M at least N, but not much larger)	* Path 5 (M >= N, but not much larger)
*	*
WRKBL = 3N + ( M+N )ILAENV( 1, 'DGEBRD', ' ', M, N, -1,	WRKBL = 3*N + LWORK_DGEBRD_MN
$ -1 )
IF( WNTQN ) THEN	IF( WNTQN ) THEN
MAXWRK = MAX( WRKBL, BDSPAC+3*N )	* Path 5n (M >= N, jobz='N')
	MAXWRK = MAX( WRKBL, 3*N + BDSPAC )
MINWRK = 3*N + MAX( M, BDSPAC )	MINWRK = 3*N + MAX( M, BDSPAC )
ELSE IF( WNTQO ) THEN	ELSE IF( WNTQO ) THEN
WRKBL = MAX( WRKBL, 3N+N	* Path 5o (M >= N, jobz='O')
$ ILAENV( 1, 'DORMBR', 'QLN', M, N, N, -1 ) )	WRKBL = MAX( WRKBL, 3*N + LWORK_DORMBR_PRT_NN )
WRKBL = MAX( WRKBL, 3N+N	WRKBL = MAX( WRKBL, 3*N + LWORK_DORMBR_QLN_MN )
$ ILAENV( 1, 'DORMBR', 'PRT', N, N, N, -1 ) )	WRKBL = MAX( WRKBL, 3*N + BDSPAC )
WRKBL = MAX( WRKBL, BDSPAC+3*N )
MAXWRK = WRKBL + M*N	MAXWRK = WRKBL + M*N
MINWRK = 3N + MAX( M, NN+BDSPAC )	MINWRK = 3N + MAX( M, NN + BDSPAC )
ELSE IF( WNTQS ) THEN	ELSE IF( WNTQS ) THEN
WRKBL = MAX( WRKBL, 3N+N	* Path 5s (M >= N, jobz='S')
$ ILAENV( 1, 'DORMBR', 'QLN', M, N, N, -1 ) )	WRKBL = MAX( WRKBL, 3*N + LWORK_DORMBR_QLN_MN )
WRKBL = MAX( WRKBL, 3N+N	WRKBL = MAX( WRKBL, 3*N + LWORK_DORMBR_PRT_NN )
$ ILAENV( 1, 'DORMBR', 'PRT', N, N, N, -1 ) )	MAXWRK = MAX( WRKBL, 3*N + BDSPAC )
MAXWRK = MAX( WRKBL, BDSPAC+3*N )
MINWRK = 3*N + MAX( M, BDSPAC )	MINWRK = 3*N + MAX( M, BDSPAC )
ELSE IF( WNTQA ) THEN	ELSE IF( WNTQA ) THEN
WRKBL = MAX( WRKBL, 3N+M	* Path 5a (M >= N, jobz='A')
$ ILAENV( 1, 'DORMBR', 'QLN', M, M, N, -1 ) )	WRKBL = MAX( WRKBL, 3*N + LWORK_DORMBR_QLN_MM )
WRKBL = MAX( WRKBL, 3N+N	WRKBL = MAX( WRKBL, 3*N + LWORK_DORMBR_PRT_NN )
$ ILAENV( 1, 'DORMBR', 'PRT', N, N, N, -1 ) )	MAXWRK = MAX( WRKBL, 3*N + BDSPAC )
MAXWRK = MAX( MAXWRK, BDSPAC+3*N )
MINWRK = 3*N + MAX( M, BDSPAC )	MINWRK = 3*N + MAX( M, BDSPAC )
END IF	END IF
END IF	END IF
Line 410	Line 447
*	*
* Compute space needed for DBDSDC	* Compute space needed for DBDSDC
*	*
MNTHR = INT( MINMN*11.0D0 / 6.0D0 )
IF( WNTQN ) THEN	IF( WNTQN ) THEN
	* dbdsdc needs only 4N (or 6N for uplo=L for LAPACK <= 3.6)
	* keep 7*N for backwards compatability.
BDSPAC = 7*M	BDSPAC = 7*M
ELSE	ELSE
BDSPAC = 3MM + 4*M	BDSPAC = 3MM + 4*M
END IF	END IF
	*
	* Compute space preferred for each routine
	CALL DGEBRD( M, N, DUM(1), M, DUM(1), DUM(1), DUM(1),
	$ DUM(1), DUM(1), -1, IERR )
	LWORK_DGEBRD_MN = INT( DUM(1) )
	*
	CALL DGEBRD( M, M, A, M, S, DUM(1), DUM(1),
	$ DUM(1), DUM(1), -1, IERR )
	LWORK_DGEBRD_MM = INT( DUM(1) )
	*
	CALL DGELQF( M, N, A, M, DUM(1), DUM(1), -1, IERR )
	LWORK_DGELQF_MN = INT( DUM(1) )
	*
	CALL DORGLQ( N, N, M, DUM(1), N, DUM(1), DUM(1), -1, IERR )
	LWORK_DORGLQ_NN = INT( DUM(1) )
	*
	CALL DORGLQ( M, N, M, A, M, DUM(1), DUM(1), -1, IERR )
	LWORK_DORGLQ_MN = INT( DUM(1) )
	*
	CALL DORGBR( 'P', M, M, M, A, N, DUM(1), DUM(1), -1, IERR )
	LWORK_DORGBR_P_MM = INT( DUM(1) )
	*
	CALL DORMBR( 'P', 'R', 'T', M, M, M, DUM(1), M,
	$ DUM(1), DUM(1), M, DUM(1), -1, IERR )
	LWORK_DORMBR_PRT_MM = INT( DUM(1) )
	*
	CALL DORMBR( 'P', 'R', 'T', M, N, M, DUM(1), M,
	$ DUM(1), DUM(1), M, DUM(1), -1, IERR )
	LWORK_DORMBR_PRT_MN = INT( DUM(1) )
	*
	CALL DORMBR( 'P', 'R', 'T', N, N, M, DUM(1), N,
	$ DUM(1), DUM(1), N, DUM(1), -1, IERR )
	LWORK_DORMBR_PRT_NN = INT( DUM(1) )
	*
	CALL DORMBR( 'Q', 'L', 'N', M, M, M, DUM(1), M,
	$ DUM(1), DUM(1), M, DUM(1), -1, IERR )
	LWORK_DORMBR_QLN_MM = INT( DUM(1) )
	*
IF( N.GE.MNTHR ) THEN	IF( N.GE.MNTHR ) THEN
IF( WNTQN ) THEN	IF( WNTQN ) THEN
*	*
* Path 1t (N much larger than M, JOBZ='N')	* Path 1t (N >> M, JOBZ='N')
*	*
WRKBL = M + M*ILAENV( 1, 'DGELQF', ' ', M, N, -1,	WRKBL = M + LWORK_DGELQF_MN
$ -1 )	WRKBL = MAX( WRKBL, 3*M + LWORK_DGEBRD_MM )
WRKBL = MAX( WRKBL, 3M+2M*	MAXWRK = MAX( WRKBL, BDSPAC + M )
$ ILAENV( 1, 'DGEBRD', ' ', M, M, -1, -1 ) )
MAXWRK = MAX( WRKBL, BDSPAC+M )
MINWRK = BDSPAC + M	MINWRK = BDSPAC + M
ELSE IF( WNTQO ) THEN	ELSE IF( WNTQO ) THEN
*	*
* Path 2t (N much larger than M, JOBZ='O')	* Path 2t (N >> M, JOBZ='O')
*	*
WRKBL = M + M*ILAENV( 1, 'DGELQF', ' ', M, N, -1, -1 )	WRKBL = M + LWORK_DGELQF_MN
WRKBL = MAX( WRKBL, M+M*ILAENV( 1, 'DORGLQ', ' ', M,	WRKBL = MAX( WRKBL, M + LWORK_DORGLQ_MN )
$ N, M, -1 ) )	WRKBL = MAX( WRKBL, 3*M + LWORK_DGEBRD_MM )
WRKBL = MAX( WRKBL, 3M+2M*	WRKBL = MAX( WRKBL, 3*M + LWORK_DORMBR_QLN_MM )
$ ILAENV( 1, 'DGEBRD', ' ', M, M, -1, -1 ) )	WRKBL = MAX( WRKBL, 3*M + LWORK_DORMBR_PRT_MM )
WRKBL = MAX( WRKBL, 3M+M	WRKBL = MAX( WRKBL, 3*M + BDSPAC )
$ ILAENV( 1, 'DORMBR', 'QLN', M, M, M, -1 ) )
WRKBL = MAX( WRKBL, 3M+M
$ ILAENV( 1, 'DORMBR', 'PRT', M, M, M, -1 ) )
WRKBL = MAX( WRKBL, BDSPAC+3*M )
MAXWRK = WRKBL + 2MM	MAXWRK = WRKBL + 2MM
MINWRK = BDSPAC + 2MM + 3*M	MINWRK = BDSPAC + 2MM + 3*M
ELSE IF( WNTQS ) THEN	ELSE IF( WNTQS ) THEN
*	*
* Path 3t (N much larger than M, JOBZ='S')	* Path 3t (N >> M, JOBZ='S')
*	*
WRKBL = M + M*ILAENV( 1, 'DGELQF', ' ', M, N, -1, -1 )	WRKBL = M + LWORK_DGELQF_MN
WRKBL = MAX( WRKBL, M+M*ILAENV( 1, 'DORGLQ', ' ', M,	WRKBL = MAX( WRKBL, M + LWORK_DORGLQ_MN )
$ N, M, -1 ) )	WRKBL = MAX( WRKBL, 3*M + LWORK_DGEBRD_MM )
WRKBL = MAX( WRKBL, 3M+2M*	WRKBL = MAX( WRKBL, 3*M + LWORK_DORMBR_QLN_MM )
$ ILAENV( 1, 'DGEBRD', ' ', M, M, -1, -1 ) )	WRKBL = MAX( WRKBL, 3*M + LWORK_DORMBR_PRT_MM )
WRKBL = MAX( WRKBL, 3M+M	WRKBL = MAX( WRKBL, 3*M + BDSPAC )
$ ILAENV( 1, 'DORMBR', 'QLN', M, M, M, -1 ) )
WRKBL = MAX( WRKBL, 3M+M
$ ILAENV( 1, 'DORMBR', 'PRT', M, M, M, -1 ) )
WRKBL = MAX( WRKBL, BDSPAC+3*M )
MAXWRK = WRKBL + M*M	MAXWRK = WRKBL + M*M
MINWRK = BDSPAC + MM + 3M	MINWRK = BDSPAC + MM + 3M
ELSE IF( WNTQA ) THEN	ELSE IF( WNTQA ) THEN
*	*
* Path 4t (N much larger than M, JOBZ='A')	* Path 4t (N >> M, JOBZ='A')
*	*
WRKBL = M + M*ILAENV( 1, 'DGELQF', ' ', M, N, -1, -1 )	WRKBL = M + LWORK_DGELQF_MN
WRKBL = MAX( WRKBL, M+N*ILAENV( 1, 'DORGLQ', ' ', N,	WRKBL = MAX( WRKBL, M + LWORK_DORGLQ_NN )
$ N, M, -1 ) )	WRKBL = MAX( WRKBL, 3*M + LWORK_DGEBRD_MM )
WRKBL = MAX( WRKBL, 3M+2M*	WRKBL = MAX( WRKBL, 3*M + LWORK_DORMBR_QLN_MM )
$ ILAENV( 1, 'DGEBRD', ' ', M, M, -1, -1 ) )	WRKBL = MAX( WRKBL, 3*M + LWORK_DORMBR_PRT_MM )
WRKBL = MAX( WRKBL, 3M+M	WRKBL = MAX( WRKBL, 3*M + BDSPAC )
$ ILAENV( 1, 'DORMBR', 'QLN', M, M, M, -1 ) )
WRKBL = MAX( WRKBL, 3M+M
$ ILAENV( 1, 'DORMBR', 'PRT', M, M, M, -1 ) )
WRKBL = MAX( WRKBL, BDSPAC+3*M )
MAXWRK = WRKBL + M*M	MAXWRK = WRKBL + M*M
MINWRK = BDSPAC + MM + 3M	MINWRK = MM + MAX( 3M + BDSPAC, M + N )
END IF	END IF
ELSE	ELSE
*	*
* Path 5t (N greater than M, but not much larger)	* Path 5t (N > M, but not much larger)
*	*
WRKBL = 3M + ( M+N )ILAENV( 1, 'DGEBRD', ' ', M, N, -1,	WRKBL = 3*M + LWORK_DGEBRD_MN
$ -1 )
IF( WNTQN ) THEN	IF( WNTQN ) THEN
MAXWRK = MAX( WRKBL, BDSPAC+3*M )	* Path 5tn (N > M, jobz='N')
	MAXWRK = MAX( WRKBL, 3*M + BDSPAC )
MINWRK = 3*M + MAX( N, BDSPAC )	MINWRK = 3*M + MAX( N, BDSPAC )
ELSE IF( WNTQO ) THEN	ELSE IF( WNTQO ) THEN
WRKBL = MAX( WRKBL, 3M+M	* Path 5to (N > M, jobz='O')
$ ILAENV( 1, 'DORMBR', 'QLN', M, M, N, -1 ) )	WRKBL = MAX( WRKBL, 3*M + LWORK_DORMBR_QLN_MM )
WRKBL = MAX( WRKBL, 3M+M	WRKBL = MAX( WRKBL, 3*M + LWORK_DORMBR_PRT_MN )
$ ILAENV( 1, 'DORMBR', 'PRT', M, N, M, -1 ) )	WRKBL = MAX( WRKBL, 3*M + BDSPAC )
WRKBL = MAX( WRKBL, BDSPAC+3*M )
MAXWRK = WRKBL + M*N	MAXWRK = WRKBL + M*N
MINWRK = 3M + MAX( N, MM+BDSPAC )	MINWRK = 3M + MAX( N, MM + BDSPAC )
ELSE IF( WNTQS ) THEN	ELSE IF( WNTQS ) THEN
WRKBL = MAX( WRKBL, 3M+M	* Path 5ts (N > M, jobz='S')
$ ILAENV( 1, 'DORMBR', 'QLN', M, M, N, -1 ) )	WRKBL = MAX( WRKBL, 3*M + LWORK_DORMBR_QLN_MM )
WRKBL = MAX( WRKBL, 3M+M	WRKBL = MAX( WRKBL, 3*M + LWORK_DORMBR_PRT_MN )
$ ILAENV( 1, 'DORMBR', 'PRT', M, N, M, -1 ) )	MAXWRK = MAX( WRKBL, 3*M + BDSPAC )
MAXWRK = MAX( WRKBL, BDSPAC+3*M )
MINWRK = 3*M + MAX( N, BDSPAC )	MINWRK = 3*M + MAX( N, BDSPAC )
ELSE IF( WNTQA ) THEN	ELSE IF( WNTQA ) THEN
WRKBL = MAX( WRKBL, 3M+M	* Path 5ta (N > M, jobz='A')
$ ILAENV( 1, 'DORMBR', 'QLN', M, M, N, -1 ) )	WRKBL = MAX( WRKBL, 3*M + LWORK_DORMBR_QLN_MM )
WRKBL = MAX( WRKBL, 3M+M	WRKBL = MAX( WRKBL, 3*M + LWORK_DORMBR_PRT_NN )
$ ILAENV( 1, 'DORMBR', 'PRT', N, N, M, -1 ) )	MAXWRK = MAX( WRKBL, 3*M + BDSPAC )
MAXWRK = MAX( WRKBL, BDSPAC+3*M )
MINWRK = 3*M + MAX( N, BDSPAC )	MINWRK = 3*M + MAX( N, BDSPAC )
END IF	END IF
END IF	END IF
END IF	END IF

MAXWRK = MAX( MAXWRK, MINWRK )	MAXWRK = MAX( MAXWRK, MINWRK )
WORK( 1 ) = MAXWRK	WORK( 1 ) = MAXWRK
*	*
Line 559	Line 619
*	*
IF( WNTQN ) THEN	IF( WNTQN ) THEN
*	*
* Path 1 (M much larger than N, JOBZ='N')	* Path 1 (M >> N, JOBZ='N')
* No singular vectors to be computed	* No singular vectors to be computed
*	*
ITAU = 1	ITAU = 1
NWORK = ITAU + N	NWORK = ITAU + N
*	*
* Compute A=Q*R	* Compute A=Q*R
* (Workspace: need 2N, prefer N+NNB)	* Workspace: need N [tau] + N [work]
	* Workspace: prefer N [tau] + N*NB [work]
*	*
CALL DGEQRF( M, N, A, LDA, WORK( ITAU ), WORK( NWORK ),	CALL DGEQRF( M, N, A, LDA, WORK( ITAU ), WORK( NWORK ),
$ LWORK-NWORK+1, IERR )	$ LWORK - NWORK + 1, IERR )
*	*
* Zero out below R	* Zero out below R
*	*
Line 580	Line 641
NWORK = ITAUP + N	NWORK = ITAUP + N
*	*
* Bidiagonalize R in A	* Bidiagonalize R in A
* (Workspace: need 4N, prefer 3N+2NNB)	* Workspace: need 3*N [e, tauq, taup] + N [work]
	* Workspace: prefer 3N [e, tauq, taup] + 2N*NB [work]
*	*
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( NWORK ), LWORK-NWORK+1,	$ WORK( ITAUP ), WORK( NWORK ), LWORK-NWORK+1,
Line 588	Line 650
NWORK = IE + N	NWORK = IE + N
*	*
* Perform bidiagonal SVD, computing singular values only	* Perform bidiagonal SVD, computing singular values only
* (Workspace: need N+BDSPAC)	* Workspace: need N [e] + BDSPAC
*	*
CALL DBDSDC( 'U', 'N', N, S, WORK( IE ), DUM, 1, DUM, 1,	CALL DBDSDC( 'U', 'N', N, S, WORK( IE ), DUM, 1, DUM, 1,
$ DUM, IDUM, WORK( NWORK ), IWORK, INFO )	$ DUM, IDUM, WORK( NWORK ), IWORK, INFO )
*	*
ELSE IF( WNTQO ) THEN	ELSE IF( WNTQO ) THEN
*	*
* Path 2 (M much larger than N, JOBZ = 'O')	* Path 2 (M >> N, JOBZ = 'O')
* N left singular vectors to be overwritten on A and	* N left singular vectors to be overwritten on A and
* N right singular vectors to be computed in VT	* N right singular vectors to be computed in VT
*	*
Line 603	Line 665
*	*
* WORK(IR) is LDWRKR by N	* WORK(IR) is LDWRKR by N
*	*
IF( LWORK.GE.LDAN+NN+3*N+BDSPAC ) THEN	IF( LWORK .GE. LDAN + NN + 3*N + BDSPAC ) THEN
LDWRKR = LDA	LDWRKR = LDA
ELSE	ELSE
LDWRKR = ( LWORK-NN-3N-BDSPAC ) / N	LDWRKR = ( LWORK - NN - 3N - BDSPAC ) / N
END IF	END IF
ITAU = IR + LDWRKR*N	ITAU = IR + LDWRKR*N
NWORK = ITAU + N	NWORK = ITAU + N
*	*
* Compute A=Q*R	* Compute A=Q*R
* (Workspace: need NN+2N, prefer NN+N+NNB)	* Workspace: need N*N [R] + N [tau] + N [work]
	* Workspace: prefer NN [R] + N [tau] + NNB [work]
*	*
CALL DGEQRF( M, N, A, LDA, WORK( ITAU ), WORK( NWORK ),	CALL DGEQRF( M, N, A, LDA, WORK( ITAU ), WORK( NWORK ),
$ LWORK-NWORK+1, IERR )	$ LWORK - NWORK + 1, IERR )
*	*
* Copy R to WORK(IR), zeroing out below it	* Copy R to WORK(IR), zeroing out below it
*	*
CALL DLACPY( 'U', N, N, A, LDA, WORK( IR ), LDWRKR )	CALL DLACPY( 'U', N, N, A, LDA, WORK( IR ), LDWRKR )
CALL DLASET( 'L', N-1, N-1, ZERO, ZERO, WORK( IR+1 ),	CALL DLASET( 'L', N - 1, N - 1, ZERO, ZERO, WORK(IR+1),
$ LDWRKR )	$ LDWRKR )
*	*
* Generate Q in A	* Generate Q in A
* (Workspace: need NN+2N, prefer NN+N+NNB)	* Workspace: need N*N [R] + N [tau] + N [work]
	* Workspace: prefer NN [R] + N [tau] + NNB [work]
*	*
CALL DORGQR( M, N, N, A, LDA, WORK( ITAU ),	CALL DORGQR( M, N, N, A, LDA, WORK( ITAU ),
$ WORK( NWORK ), LWORK-NWORK+1, IERR )	$ WORK( NWORK ), LWORK - NWORK + 1, IERR )
IE = ITAU	IE = ITAU
ITAUQ = IE + N	ITAUQ = IE + N
ITAUP = ITAUQ + N	ITAUP = ITAUQ + N
NWORK = ITAUP + N	NWORK = ITAUP + N
*	*
* Bidiagonalize R in VT, copying result to WORK(IR)	* Bidiagonalize R in WORK(IR)
* (Workspace: need NN+4N, prefer NN+3N+2NNB)	* Workspace: need NN [R] + 3N [e, tauq, taup] + N [work]
	* Workspace: prefer NN [R] + 3N [e, tauq, taup] + 2NNB [work]
*	*
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( NWORK ),	$ WORK( ITAUQ ), WORK( ITAUP ), WORK( NWORK ),
$ LWORK-NWORK+1, IERR )	$ LWORK - NWORK + 1, IERR )
*	*
* WORK(IU) is N by N	* WORK(IU) is N by N
*	*
Line 648	Line 713
* Perform bidiagonal SVD, computing left singular vectors	* Perform bidiagonal SVD, computing left singular vectors
* of bidiagonal matrix in WORK(IU) and computing right	* of bidiagonal matrix in WORK(IU) and computing right
* singular vectors of bidiagonal matrix in VT	* singular vectors of bidiagonal matrix in VT
* (Workspace: need N+N*N+BDSPAC)	* Workspace: need NN [R] + 3N [e, tauq, taup] + N*N [U] + BDSPAC
*	*
CALL DBDSDC( 'U', 'I', N, S, WORK( IE ), WORK( IU ), N,	CALL DBDSDC( 'U', 'I', N, S, WORK( IE ), WORK( IU ), N,
$ VT, LDVT, DUM, IDUM, WORK( NWORK ), IWORK,	$ VT, LDVT, DUM, IDUM, WORK( NWORK ), IWORK,
Line 656	Line 721
*	*
* Overwrite WORK(IU) by left singular vectors of R	* Overwrite WORK(IU) by left singular vectors of R
* and VT by right singular vectors of R	* and VT by right singular vectors of R
* (Workspace: need 2NN+3N, prefer 2NN+2N+N*NB)	* Workspace: need NN [R] + 3N [e, tauq, taup] + N*N [U] + N [work]
	* Workspace: prefer NN [R] + 3N [e, tauq, taup] + NN [U] + NNB [work]
*	*
CALL DORMBR( 'Q', 'L', 'N', N, N, N, WORK( IR ), LDWRKR,	CALL DORMBR( 'Q', 'L', 'N', N, N, N, WORK( IR ), LDWRKR,
$ WORK( ITAUQ ), WORK( IU ), N, WORK( NWORK ),	$ WORK( ITAUQ ), WORK( IU ), N, WORK( NWORK ),
$ LWORK-NWORK+1, IERR )	$ LWORK - NWORK + 1, IERR )
CALL DORMBR( 'P', 'R', 'T', N, N, N, WORK( IR ), LDWRKR,	CALL DORMBR( 'P', 'R', 'T', N, N, N, WORK( IR ), LDWRKR,
$ WORK( ITAUP ), VT, LDVT, WORK( NWORK ),	$ WORK( ITAUP ), VT, LDVT, WORK( NWORK ),
$ LWORK-NWORK+1, IERR )	$ LWORK - NWORK + 1, IERR )
*	*
* Multiply Q in A by left singular vectors of R in	* Multiply Q in A by left singular vectors of R in
* WORK(IU), storing result in WORK(IR) and copying to A	* WORK(IU), storing result in WORK(IR) and copying to A
* (Workspace: need 2NN, prefer NN+MN)	* Workspace: need NN [R] + 3N [e, tauq, taup] + N*N [U]
	* Workspace: prefer MN [R] + 3N [e, tauq, taup] + N*N [U]
*	*
DO 10 I = 1, M, LDWRKR	DO 10 I = 1, M, LDWRKR
CHUNK = MIN( M-I+1, LDWRKR )	CHUNK = MIN( M - I + 1, LDWRKR )
CALL DGEMM( 'N', 'N', CHUNK, N, N, ONE, A( I, 1 ),	CALL DGEMM( 'N', 'N', CHUNK, N, N, ONE, A( I, 1 ),
$ LDA, WORK( IU ), N, ZERO, WORK( IR ),	$ LDA, WORK( IU ), N, ZERO, WORK( IR ),
$ LDWRKR )	$ LDWRKR )
Line 680	Line 747
*	*
ELSE IF( WNTQS ) THEN	ELSE IF( WNTQS ) THEN
*	*
* Path 3 (M much larger than N, JOBZ='S')	* Path 3 (M >> N, JOBZ='S')
* N left singular vectors to be computed in U and	* N left singular vectors to be computed in U and
* N right singular vectors to be computed in VT	* N right singular vectors to be computed in VT
*	*
Line 693	Line 760
NWORK = ITAU + N	NWORK = ITAU + N
*	*
* Compute A=Q*R	* Compute A=Q*R
* (Workspace: need NN+2N, prefer NN+N+NNB)	* Workspace: need N*N [R] + N [tau] + N [work]
	* Workspace: prefer NN [R] + N [tau] + NNB [work]
*	*
CALL DGEQRF( M, N, A, LDA, WORK( ITAU ), WORK( NWORK ),	CALL DGEQRF( M, N, A, LDA, WORK( ITAU ), WORK( NWORK ),
$ LWORK-NWORK+1, IERR )	$ LWORK - NWORK + 1, IERR )
*	*
* Copy R to WORK(IR), zeroing out below it	* Copy R to WORK(IR), zeroing out below it
*	*
CALL DLACPY( 'U', N, N, A, LDA, WORK( IR ), LDWRKR )	CALL DLACPY( 'U', N, N, A, LDA, WORK( IR ), LDWRKR )
CALL DLASET( 'L', N-1, N-1, ZERO, ZERO, WORK( IR+1 ),	CALL DLASET( 'L', N - 1, N - 1, ZERO, ZERO, WORK(IR+1),
$ LDWRKR )	$ LDWRKR )
*	*
* Generate Q in A	* Generate Q in A
* (Workspace: need NN+2N, prefer NN+N+NNB)	* Workspace: need N*N [R] + N [tau] + N [work]
	* Workspace: prefer NN [R] + N [tau] + NNB [work]
*	*
CALL DORGQR( M, N, N, A, LDA, WORK( ITAU ),	CALL DORGQR( M, N, N, A, LDA, WORK( ITAU ),
$ WORK( NWORK ), LWORK-NWORK+1, IERR )	$ WORK( NWORK ), LWORK - NWORK + 1, IERR )
IE = ITAU	IE = ITAU
ITAUQ = IE + N	ITAUQ = IE + N
ITAUP = ITAUQ + N	ITAUP = ITAUQ + N
NWORK = ITAUP + N	NWORK = ITAUP + N
*	*
* Bidiagonalize R in WORK(IR)	* Bidiagonalize R in WORK(IR)
* (Workspace: need NN+4N, prefer NN+3N+2NNB)	* Workspace: need NN [R] + 3N [e, tauq, taup] + N [work]
	* Workspace: prefer NN [R] + 3N [e, tauq, taup] + 2NNB [work]
*	*
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( NWORK ),	$ WORK( ITAUQ ), WORK( ITAUP ), WORK( NWORK ),
$ LWORK-NWORK+1, IERR )	$ LWORK - NWORK + 1, IERR )
*	*
* Perform bidiagonal SVD, computing left singular vectors	* Perform bidiagonal SVD, computing left singular vectors
* of bidiagoal matrix in U and computing right singular	* of bidiagoal matrix in U and computing right singular
* vectors of bidiagonal matrix in VT	* vectors of bidiagonal matrix in VT
* (Workspace: need N+BDSPAC)	* Workspace: need NN [R] + 3N [e, tauq, taup] + BDSPAC
*	*
CALL DBDSDC( 'U', 'I', N, S, WORK( IE ), U, LDU, VT,	CALL DBDSDC( 'U', 'I', N, S, WORK( IE ), U, LDU, VT,
$ LDVT, DUM, IDUM, WORK( NWORK ), IWORK,	$ LDVT, DUM, IDUM, WORK( NWORK ), IWORK,
Line 732	Line 802
*	*
* Overwrite U by left singular vectors of R and VT	* Overwrite U by left singular vectors of R and VT
* by right singular vectors of R	* by right singular vectors of R
* (Workspace: need NN+3N, prefer NN+2N+N*NB)	* Workspace: need NN [R] + 3N [e, tauq, taup] + N [work]
	* Workspace: prefer NN [R] + 3N [e, tauq, taup] + N*NB [work]
*	*
CALL DORMBR( 'Q', 'L', 'N', N, N, N, WORK( IR ), LDWRKR,	CALL DORMBR( 'Q', 'L', 'N', N, N, N, WORK( IR ), LDWRKR,
$ WORK( ITAUQ ), U, LDU, WORK( NWORK ),	$ WORK( ITAUQ ), U, LDU, WORK( NWORK ),
$ LWORK-NWORK+1, IERR )	$ LWORK - NWORK + 1, IERR )
*	*
CALL DORMBR( 'P', 'R', 'T', N, N, N, WORK( IR ), LDWRKR,	CALL DORMBR( 'P', 'R', 'T', N, N, N, WORK( IR ), LDWRKR,
$ WORK( ITAUP ), VT, LDVT, WORK( NWORK ),	$ WORK( ITAUP ), VT, LDVT, WORK( NWORK ),
$ LWORK-NWORK+1, IERR )	$ LWORK - NWORK + 1, IERR )
*	*
* Multiply Q in A by left singular vectors of R in	* Multiply Q in A by left singular vectors of R in
* WORK(IR), storing result in U	* WORK(IR), storing result in U
* (Workspace: need N*N)	* Workspace: need N*N [R]
*	*
CALL DLACPY( 'F', N, N, U, LDU, WORK( IR ), LDWRKR )	CALL DLACPY( 'F', N, N, U, LDU, WORK( IR ), LDWRKR )
CALL DGEMM( 'N', 'N', M, N, N, ONE, A, LDA, WORK( IR ),	CALL DGEMM( 'N', 'N', M, N, N, ONE, A, LDA, WORK( IR ),
Line 752	Line 823
*	*
ELSE IF( WNTQA ) THEN	ELSE IF( WNTQA ) THEN
*	*
* Path 4 (M much larger than N, JOBZ='A')	* Path 4 (M >> N, JOBZ='A')
* M left singular vectors to be computed in U and	* M left singular vectors to be computed in U and
* N right singular vectors to be computed in VT	* N right singular vectors to be computed in VT
*	*
Line 765	Line 836
NWORK = ITAU + N	NWORK = ITAU + N
*	*
* Compute A=Q*R, copying result to U	* Compute A=Q*R, copying result to U
* (Workspace: need NN+N+M, prefer NN+N+M*NB)	* Workspace: need N*N [U] + N [tau] + N [work]
	* Workspace: prefer NN [U] + N [tau] + NNB [work]
*	*
CALL DGEQRF( M, N, A, LDA, WORK( ITAU ), WORK( NWORK ),	CALL DGEQRF( M, N, A, LDA, WORK( ITAU ), WORK( NWORK ),
$ LWORK-NWORK+1, IERR )	$ LWORK - NWORK + 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 N*N [U] + N [tau] + M [work]
	* Workspace: prefer NN [U] + N [tau] + MNB [work]
CALL DORGQR( M, M, N, U, LDU, WORK( ITAU ),	CALL DORGQR( M, M, N, U, LDU, WORK( ITAU ),
$ WORK( NWORK ), LWORK-NWORK+1, IERR )	$ WORK( NWORK ), LWORK - NWORK + 1, IERR )
*	*
* Produce R in A, zeroing out other entries	* Produce R in A, zeroing out other entries
*	*
Line 785	Line 858
NWORK = ITAUP + N	NWORK = ITAUP + N
*	*
* Bidiagonalize R in A	* Bidiagonalize R in A
* (Workspace: need NN+4N, prefer NN+3N+2NNB)	* Workspace: need NN [U] + 3N [e, tauq, taup] + N [work]
	* Workspace: prefer NN [U] + 3N [e, tauq, taup] + 2NNB [work]
*	*
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( NWORK ), LWORK-NWORK+1,	$ WORK( ITAUP ), WORK( NWORK ), LWORK-NWORK+1,
Line 794	Line 868
* Perform bidiagonal SVD, computing left singular vectors	* Perform bidiagonal SVD, computing left singular vectors
* of bidiagonal matrix in WORK(IU) and computing right	* of bidiagonal matrix in WORK(IU) and computing right
* singular vectors of bidiagonal matrix in VT	* singular vectors of bidiagonal matrix in VT
* (Workspace: need N+N*N+BDSPAC)	* Workspace: need NN [U] + 3N [e, tauq, taup] + BDSPAC
*	*
CALL DBDSDC( 'U', 'I', N, S, WORK( IE ), WORK( IU ), N,	CALL DBDSDC( 'U', 'I', N, S, WORK( IE ), WORK( IU ), N,
$ VT, LDVT, DUM, IDUM, WORK( NWORK ), IWORK,	$ VT, LDVT, DUM, IDUM, WORK( NWORK ), IWORK,
Line 802	Line 876
*	*
* Overwrite WORK(IU) by left singular vectors of R and VT	* Overwrite WORK(IU) by left singular vectors of R and VT
* by right singular vectors of R	* by right singular vectors of R
* (Workspace: need NN+3N, prefer NN+2N+N*NB)	* Workspace: need NN [U] + 3N [e, tauq, taup] + N [work]
	* Workspace: prefer NN [U] + 3N [e, tauq, taup] + N*NB [work]
*	*
CALL DORMBR( 'Q', 'L', 'N', N, N, N, A, LDA,	CALL DORMBR( 'Q', 'L', 'N', N, N, N, A, LDA,
$ WORK( ITAUQ ), WORK( IU ), LDWRKU,	$ WORK( ITAUQ ), WORK( IU ), LDWRKU,
$ WORK( NWORK ), LWORK-NWORK+1, IERR )	$ WORK( NWORK ), LWORK - NWORK + 1, IERR )
CALL DORMBR( 'P', 'R', 'T', N, N, N, A, LDA,	CALL DORMBR( 'P', 'R', 'T', N, N, N, A, LDA,
$ WORK( ITAUP ), VT, LDVT, WORK( NWORK ),	$ WORK( ITAUP ), VT, LDVT, WORK( NWORK ),
$ LWORK-NWORK+1, IERR )	$ LWORK - NWORK + 1, IERR )
*	*
* Multiply Q in U by left singular vectors of R in	* Multiply Q in U by left singular vectors of R in
* WORK(IU), storing result in A	* WORK(IU), storing result in A
* (Workspace: need N*N)	* Workspace: need N*N [U]
*	*
CALL DGEMM( 'N', 'N', M, N, N, ONE, U, LDU, WORK( IU ),	CALL DGEMM( 'N', 'N', M, N, N, ONE, U, LDU, WORK( IU ),
$ LDWRKU, ZERO, A, LDA )	$ LDWRKU, ZERO, A, LDA )
Line 828	Line 903
*	*
* M .LT. MNTHR	* M .LT. MNTHR
*	*
* Path 5 (M at least N, but not much larger)	* Path 5 (M >= N, but not much larger)
* Reduce to bidiagonal form without QR decomposition	* Reduce to bidiagonal form without QR decomposition
*	*
IE = 1	IE = 1
Line 837	Line 912
NWORK = ITAUP + N	NWORK = ITAUP + N
*	*
* Bidiagonalize A	* Bidiagonalize A
* (Workspace: need 3N+M, prefer 3N+(M+N)*NB)	* Workspace: need 3*N [e, tauq, taup] + M [work]
	* Workspace: prefer 3N [e, tauq, taup] + (M+N)NB [work]
*	*
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( NWORK ), LWORK-NWORK+1,	$ WORK( ITAUP ), WORK( NWORK ), LWORK-NWORK+1,
$ IERR )	$ IERR )
IF( WNTQN ) THEN	IF( WNTQN ) THEN
*	*
	* Path 5n (M >= N, JOBZ='N')
* Perform bidiagonal SVD, only computing singular values	* Perform bidiagonal SVD, only computing singular values
* (Workspace: need N+BDSPAC)	* Workspace: need 3*N [e, tauq, taup] + BDSPAC
*	*
CALL DBDSDC( 'U', 'N', N, S, WORK( IE ), DUM, 1, DUM, 1,	CALL DBDSDC( 'U', 'N', N, S, WORK( IE ), DUM, 1, DUM, 1,
$ DUM, IDUM, WORK( NWORK ), IWORK, INFO )	$ DUM, IDUM, WORK( NWORK ), IWORK, INFO )
ELSE IF( WNTQO ) THEN	ELSE IF( WNTQO ) THEN
	* Path 5o (M >= N, JOBZ='O')
IU = NWORK	IU = NWORK
IF( LWORK.GE.MN+3N+BDSPAC ) THEN	IF( LWORK .GE. MN + 3N + BDSPAC ) THEN
*	*
* WORK( IU ) is M by N	* WORK( IU ) is M by N
*	*
Line 859	Line 937
NWORK = IU + LDWRKU*N	NWORK = IU + LDWRKU*N
CALL DLASET( 'F', M, N, ZERO, ZERO, WORK( IU ),	CALL DLASET( 'F', M, N, ZERO, ZERO, WORK( IU ),
$ LDWRKU )	$ LDWRKU )
	* IR is unused; silence compile warnings
	IR = -1
ELSE	ELSE
*	*
* WORK( IU ) is N by N	* WORK( IU ) is N by N
Line 869	Line 949
* WORK(IR) is LDWRKR by N	* WORK(IR) is LDWRKR by N
*	*
IR = NWORK	IR = NWORK
LDWRKR = ( LWORK-NN-3N ) / N	LDWRKR = ( LWORK - NN - 3N ) / N
END IF	END IF
NWORK = IU + LDWRKU*N	NWORK = IU + LDWRKU*N
*	*
* Perform bidiagonal SVD, computing left singular vectors	* Perform bidiagonal SVD, computing left singular vectors
* of bidiagonal matrix in WORK(IU) and computing right	* of bidiagonal matrix in WORK(IU) and computing right
* singular vectors of bidiagonal matrix in VT	* singular vectors of bidiagonal matrix in VT
* (Workspace: need N+N*N+BDSPAC)	* Workspace: need 3N [e, tauq, taup] + NN [U] + BDSPAC
*	*
CALL DBDSDC( 'U', 'I', N, S, WORK( IE ), WORK( IU ),	CALL DBDSDC( 'U', 'I', N, S, WORK( IE ), WORK( IU ),
$ LDWRKU, VT, LDVT, DUM, IDUM, WORK( NWORK ),	$ LDWRKU, VT, LDVT, DUM, IDUM, WORK( NWORK ),
$ IWORK, INFO )	$ IWORK, INFO )
*	*
* Overwrite VT by right singular vectors of A	* Overwrite VT by right singular vectors of A
* (Workspace: need NN+2N, prefer NN+N+NNB)	* Workspace: need 3N [e, tauq, taup] + NN [U] + N [work]
	* Workspace: prefer 3N [e, tauq, taup] + NN [U] + N*NB [work]
*	*
CALL DORMBR( 'P', 'R', 'T', N, N, N, A, LDA,	CALL DORMBR( 'P', 'R', 'T', N, N, N, A, LDA,
$ WORK( ITAUP ), VT, LDVT, WORK( NWORK ),	$ WORK( ITAUP ), VT, LDVT, WORK( NWORK ),
$ LWORK-NWORK+1, IERR )	$ LWORK - NWORK + 1, IERR )
*	*
IF( LWORK.GE.MN+3N+BDSPAC ) THEN	IF( LWORK .GE. MN + 3N + BDSPAC ) THEN
*	*
	* Path 5o-fast
* Overwrite WORK(IU) by left singular vectors of A	* Overwrite WORK(IU) by left singular vectors of A
* (Workspace: need NN+2N, prefer NN+N+NNB)	* Workspace: need 3N [e, tauq, taup] + MN [U] + N [work]
	* Workspace: prefer 3N [e, tauq, taup] + MN [U] + N*NB [work]
*	*
CALL DORMBR( 'Q', 'L', 'N', M, N, N, A, LDA,	CALL DORMBR( 'Q', 'L', 'N', M, N, N, A, LDA,
$ WORK( ITAUQ ), WORK( IU ), LDWRKU,	$ WORK( ITAUQ ), WORK( IU ), LDWRKU,
$ WORK( NWORK ), LWORK-NWORK+1, IERR )	$ WORK( NWORK ), LWORK - NWORK + 1, IERR )
*	*
* Copy left singular vectors of A from WORK(IU) to A	* Copy left singular vectors of A from WORK(IU) to A
*	*
CALL DLACPY( 'F', M, N, WORK( IU ), LDWRKU, A, LDA )	CALL DLACPY( 'F', M, N, WORK( IU ), LDWRKU, A, LDA )
ELSE	ELSE
*	*
	* Path 5o-slow
* Generate Q in A	* Generate Q in A
* (Workspace: need NN+2N, prefer NN+N+NNB)	* Workspace: need 3N [e, tauq, taup] + NN [U] + N [work]
	* Workspace: prefer 3N [e, tauq, taup] + NN [U] + N*NB [work]
*	*
CALL DORGBR( 'Q', M, N, N, A, LDA, WORK( ITAUQ ),	CALL DORGBR( 'Q', M, N, N, A, LDA, WORK( ITAUQ ),
$ WORK( NWORK ), LWORK-NWORK+1, IERR )	$ WORK( NWORK ), LWORK - NWORK + 1, IERR )
*	*
* Multiply Q in A by left singular vectors of	* Multiply Q in A by left singular vectors of
* bidiagonal matrix in WORK(IU), storing result in	* bidiagonal matrix in WORK(IU), storing result in
* WORK(IR) and copying to A	* WORK(IR) and copying to A
* (Workspace: need 2NN, prefer NN+MN)	* Workspace: need 3N [e, tauq, taup] + NN [U] + NB*N [R]
	* Workspace: prefer 3N [e, tauq, taup] + NN [U] + M*N [R]
*	*
DO 20 I = 1, M, LDWRKR	DO 20 I = 1, M, LDWRKR
CHUNK = MIN( M-I+1, LDWRKR )	CHUNK = MIN( M - I + 1, LDWRKR )
CALL DGEMM( 'N', 'N', CHUNK, N, N, ONE, A( I, 1 ),	CALL DGEMM( 'N', 'N', CHUNK, N, N, ONE, A( I, 1 ),
$ LDA, WORK( IU ), LDWRKU, ZERO,	$ LDA, WORK( IU ), LDWRKU, ZERO,
$ WORK( IR ), LDWRKR )	$ WORK( IR ), LDWRKR )
Line 926	Line 1012
*	*
ELSE IF( WNTQS ) THEN	ELSE IF( WNTQS ) THEN
*	*
	* Path 5s (M >= N, JOBZ='S')
* Perform bidiagonal SVD, computing left singular vectors	* Perform bidiagonal SVD, computing left singular vectors
* of bidiagonal matrix in U and computing right singular	* of bidiagonal matrix in U and computing right singular
* vectors of bidiagonal matrix in VT	* vectors of bidiagonal matrix in VT
* (Workspace: need N+BDSPAC)	* Workspace: need 3*N [e, tauq, taup] + BDSPAC
*	*
CALL DLASET( 'F', M, N, ZERO, ZERO, U, LDU )	CALL DLASET( 'F', M, N, ZERO, ZERO, U, LDU )
CALL DBDSDC( 'U', 'I', N, S, WORK( IE ), U, LDU, VT,	CALL DBDSDC( 'U', 'I', N, S, WORK( IE ), U, LDU, VT,
Line 938	Line 1025
*	*
* Overwrite U by left singular vectors of A and VT	* Overwrite U by left singular vectors of A and VT
* by right singular vectors of A	* by right singular vectors of A
* (Workspace: need 3N, prefer 2N+N*NB)	* Workspace: need 3*N [e, tauq, taup] + N [work]
	* Workspace: prefer 3N [e, tauq, taup] + NNB [work]
*	*
CALL DORMBR( 'Q', 'L', 'N', M, N, N, A, LDA,	CALL DORMBR( 'Q', 'L', 'N', M, N, N, A, LDA,
$ WORK( ITAUQ ), U, LDU, WORK( NWORK ),	$ WORK( ITAUQ ), U, LDU, WORK( NWORK ),
$ LWORK-NWORK+1, IERR )	$ LWORK - NWORK + 1, IERR )
CALL DORMBR( 'P', 'R', 'T', N, N, N, A, LDA,	CALL DORMBR( 'P', 'R', 'T', N, N, N, A, LDA,
$ WORK( ITAUP ), VT, LDVT, WORK( NWORK ),	$ WORK( ITAUP ), VT, LDVT, WORK( NWORK ),
$ LWORK-NWORK+1, IERR )	$ LWORK - NWORK + 1, IERR )
ELSE IF( WNTQA ) THEN	ELSE IF( WNTQA ) THEN
*	*
	* Path 5a (M >= N, JOBZ='A')
* Perform bidiagonal SVD, computing left singular vectors	* Perform bidiagonal SVD, computing left singular vectors
* of bidiagonal matrix in U and computing right singular	* of bidiagonal matrix in U and computing right singular
* vectors of bidiagonal matrix in VT	* vectors of bidiagonal matrix in VT
* (Workspace: need N+BDSPAC)	* Workspace: need 3*N [e, tauq, taup] + BDSPAC
*	*
CALL DLASET( 'F', M, M, ZERO, ZERO, U, LDU )	CALL DLASET( 'F', M, M, ZERO, ZERO, U, LDU )
CALL DBDSDC( 'U', 'I', N, S, WORK( IE ), U, LDU, VT,	CALL DBDSDC( 'U', 'I', N, S, WORK( IE ), U, LDU, VT,
Line 961	Line 1050
* Set the right corner of U to identity matrix	* Set the right corner of U to identity matrix
*	*
IF( M.GT.N ) THEN	IF( M.GT.N ) THEN
CALL DLASET( 'F', M-N, M-N, ZERO, ONE, U( N+1, N+1 ),	CALL DLASET( 'F', M - N, M - N, ZERO, ONE, U(N+1,N+1),
$ LDU )	$ LDU )
END IF	END IF
*	*
* Overwrite U by left singular vectors of A and VT	* Overwrite U by left singular vectors of A and VT
* by right singular vectors of A	* by right singular vectors of A
* (Workspace: need NN+2N+M, prefer NN+2N+M*NB)	* Workspace: need 3*N [e, tauq, taup] + M [work]
	* Workspace: prefer 3N [e, tauq, taup] + MNB [work]
*	*
CALL DORMBR( 'Q', 'L', 'N', M, M, N, A, LDA,	CALL DORMBR( 'Q', 'L', 'N', M, M, N, A, LDA,
$ WORK( ITAUQ ), U, LDU, WORK( NWORK ),	$ WORK( ITAUQ ), U, LDU, WORK( NWORK ),
$ LWORK-NWORK+1, IERR )	$ LWORK - NWORK + 1, IERR )
CALL DORMBR( 'P', 'R', 'T', N, N, M, A, LDA,	CALL DORMBR( 'P', 'R', 'T', N, N, M, A, LDA,
$ WORK( ITAUP ), VT, LDVT, WORK( NWORK ),	$ WORK( ITAUP ), VT, LDVT, WORK( NWORK ),
$ LWORK-NWORK+1, IERR )	$ LWORK - NWORK + 1, IERR )
END IF	END IF
*	*
END IF	END IF
Line 989	Line 1079
*	*
IF( WNTQN ) THEN	IF( WNTQN ) THEN
*	*
* Path 1t (N much larger than M, JOBZ='N')	* Path 1t (N >> M, JOBZ='N')
* No singular vectors to be computed	* No singular vectors to be computed
*	*
ITAU = 1	ITAU = 1
NWORK = ITAU + M	NWORK = ITAU + M
*	*
* Compute A=L*Q	* Compute A=L*Q
* (Workspace: need 2M, prefer M+MNB)	* Workspace: need M [tau] + M [work]
	* Workspace: prefer M [tau] + M*NB [work]
*	*
CALL DGELQF( M, N, A, LDA, WORK( ITAU ), WORK( NWORK ),	CALL DGELQF( M, N, A, LDA, WORK( ITAU ), WORK( NWORK ),
$ LWORK-NWORK+1, IERR )	$ LWORK - NWORK + 1, IERR )
*	*
* Zero out above L	* Zero out above L
*	*
Line 1010	Line 1101
NWORK = ITAUP + M	NWORK = ITAUP + M
*	*
* Bidiagonalize L in A	* Bidiagonalize L in A
* (Workspace: need 4M, prefer 3M+2MNB)	* Workspace: need 3*M [e, tauq, taup] + M [work]
	* Workspace: prefer 3M [e, tauq, taup] + 2M*NB [work]
*	*
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( NWORK ), LWORK-NWORK+1,	$ WORK( ITAUP ), WORK( NWORK ), LWORK-NWORK+1,
Line 1018	Line 1110
NWORK = IE + M	NWORK = IE + M
*	*
* Perform bidiagonal SVD, computing singular values only	* Perform bidiagonal SVD, computing singular values only
* (Workspace: need M+BDSPAC)	* Workspace: need M [e] + BDSPAC
*	*
CALL DBDSDC( 'U', 'N', M, S, WORK( IE ), DUM, 1, DUM, 1,	CALL DBDSDC( 'U', 'N', M, S, WORK( IE ), DUM, 1, DUM, 1,
$ DUM, IDUM, WORK( NWORK ), IWORK, INFO )	$ DUM, IDUM, WORK( NWORK ), IWORK, INFO )
*	*
ELSE IF( WNTQO ) THEN	ELSE IF( WNTQO ) THEN
*	*
* Path 2t (N much larger than M, JOBZ='O')	* Path 2t (N >> M, JOBZ='O')
* M right singular vectors to be overwritten on A and	* M right singular vectors to be overwritten on A and
* M left singular vectors to be computed in U	* M left singular vectors to be computed in U
*	*
IVT = 1	IVT = 1
*	*
* IVT is M by M	* WORK(IVT) is M by M
	* WORK(IL) is M by M; it is later resized to M by chunk for gemm
*	*
IL = IVT + M*M	IL = IVT + M*M
IF( LWORK.GE.MN+MM+3*M+BDSPAC ) THEN	IF( LWORK .GE. MN + MM + 3*M + BDSPAC ) THEN
*
* WORK(IL) is M by N
*
LDWRKL = M	LDWRKL = M
CHUNK = N	CHUNK = N
ELSE	ELSE
LDWRKL = M	LDWRKL = M
CHUNK = ( LWORK-M*M ) / M	CHUNK = ( LWORK - M*M ) / M
END IF	END IF
ITAU = IL + LDWRKL*M	ITAU = IL + LDWRKL*M
NWORK = ITAU + M	NWORK = ITAU + M
*	*
* Compute A=L*Q	* Compute A=L*Q
* (Workspace: need MM+2M, prefer MM+M+MNB)	* Workspace: need MM [VT] + MM [L] + M [tau] + M [work]
	* Workspace: prefer MM [VT] + MM [L] + M [tau] + M*NB [work]
*	*
CALL DGELQF( M, N, A, LDA, WORK( ITAU ), WORK( NWORK ),	CALL DGELQF( M, N, A, LDA, WORK( ITAU ), WORK( NWORK ),
$ LWORK-NWORK+1, IERR )	$ LWORK - NWORK + 1, IERR )
*	*
* Copy L to WORK(IL), zeroing about above it	* Copy L to WORK(IL), zeroing about above it
*	*
CALL DLACPY( 'L', M, M, A, LDA, WORK( IL ), LDWRKL )	CALL DLACPY( 'L', M, M, A, LDA, WORK( IL ), LDWRKL )
CALL DLASET( 'U', M-1, M-1, ZERO, ZERO,	CALL DLASET( 'U', M - 1, M - 1, ZERO, ZERO,
$ WORK( IL+LDWRKL ), LDWRKL )	$ WORK( IL + LDWRKL ), LDWRKL )
*	*
* Generate Q in A	* Generate Q in A
* (Workspace: need MM+2M, prefer MM+M+MNB)	* Workspace: need MM [VT] + MM [L] + M [tau] + M [work]
	* Workspace: prefer MM [VT] + MM [L] + M [tau] + M*NB [work]
*	*
CALL DORGLQ( M, N, M, A, LDA, WORK( ITAU ),	CALL DORGLQ( M, N, M, A, LDA, WORK( ITAU ),
$ WORK( NWORK ), LWORK-NWORK+1, IERR )	$ WORK( NWORK ), LWORK - NWORK + 1, IERR )
IE = ITAU	IE = ITAU
ITAUQ = IE + M	ITAUQ = IE + M
ITAUP = ITAUQ + M	ITAUP = ITAUQ + M
NWORK = ITAUP + M	NWORK = ITAUP + M
*	*
* Bidiagonalize L in WORK(IL)	* Bidiagonalize L in WORK(IL)
* (Workspace: need MM+4M, prefer MM+3M+2MNB)	* Workspace: need MM [VT] + MM [L] + 3*M [e, tauq, taup] + M [work]
	* Workspace: prefer MM [VT] + MM [L] + 3M [e, tauq, taup] + 2M*NB [work]
*	*
CALL DGEBRD( M, M, WORK( IL ), LDWRKL, S, WORK( IE ),	CALL DGEBRD( M, M, WORK( IL ), LDWRKL, S, WORK( IE ),
$ WORK( ITAUQ ), WORK( ITAUP ), WORK( NWORK ),	$ WORK( ITAUQ ), WORK( ITAUP ), WORK( NWORK ),
$ LWORK-NWORK+1, IERR )	$ LWORK - NWORK + 1, IERR )
*	*
* Perform bidiagonal SVD, computing left singular vectors	* Perform bidiagonal SVD, computing left singular vectors
* of bidiagonal matrix in U, and computing right singular	* of bidiagonal matrix in U, and computing right singular
* vectors of bidiagonal matrix in WORK(IVT)	* vectors of bidiagonal matrix in WORK(IVT)
* (Workspace: need M+M*M+BDSPAC)	* Workspace: need MM [VT] + MM [L] + 3*M [e, tauq, taup] + BDSPAC
*	*
CALL DBDSDC( 'U', 'I', M, S, WORK( IE ), U, LDU,	CALL DBDSDC( 'U', 'I', M, S, WORK( IE ), U, LDU,
$ WORK( IVT ), M, DUM, IDUM, WORK( NWORK ),	$ WORK( IVT ), M, DUM, IDUM, WORK( NWORK ),
Line 1087	Line 1180
*	*
* Overwrite U by left singular vectors of L and WORK(IVT)	* Overwrite U by left singular vectors of L and WORK(IVT)
* by right singular vectors of L	* by right singular vectors of L
* (Workspace: need 2MM+3M, prefer 2MM+2M+M*NB)	* Workspace: need MM [VT] + MM [L] + 3*M [e, tauq, taup] + M [work]
	* Workspace: prefer MM [VT] + MM [L] + 3M [e, tauq, taup] + MNB [work]
*	*
CALL DORMBR( 'Q', 'L', 'N', M, M, M, WORK( IL ), LDWRKL,	CALL DORMBR( 'Q', 'L', 'N', M, M, M, WORK( IL ), LDWRKL,
$ WORK( ITAUQ ), U, LDU, WORK( NWORK ),	$ WORK( ITAUQ ), U, LDU, WORK( NWORK ),
$ LWORK-NWORK+1, IERR )	$ LWORK - NWORK + 1, IERR )
CALL DORMBR( 'P', 'R', 'T', M, M, M, WORK( IL ), LDWRKL,	CALL DORMBR( 'P', 'R', 'T', M, M, M, WORK( IL ), LDWRKL,
$ WORK( ITAUP ), WORK( IVT ), M,	$ WORK( ITAUP ), WORK( IVT ), M,
$ WORK( NWORK ), LWORK-NWORK+1, IERR )	$ WORK( NWORK ), LWORK - NWORK + 1, IERR )
*	*
* Multiply right singular vectors of L in WORK(IVT) by Q	* Multiply right singular vectors of L in WORK(IVT) by Q
* in A, storing result in WORK(IL) and copying to A	* in A, storing result in WORK(IL) and copying to A
* (Workspace: need 2MM, prefer MM+MN)	* Workspace: need MM [VT] + MM [L]
	* Workspace: prefer MM [VT] + MN [L]
	* At this point, L is resized as M by chunk.
*	*
DO 30 I = 1, N, CHUNK	DO 30 I = 1, N, CHUNK
BLK = MIN( N-I+1, CHUNK )	BLK = MIN( N - I + 1, CHUNK )
CALL DGEMM( 'N', 'N', M, BLK, M, ONE, WORK( IVT ), M,	CALL DGEMM( 'N', 'N', M, BLK, M, ONE, WORK( IVT ), M,
$ A( 1, I ), LDA, ZERO, WORK( IL ), LDWRKL )	$ A( 1, I ), LDA, ZERO, WORK( IL ), LDWRKL )
CALL DLACPY( 'F', M, BLK, WORK( IL ), LDWRKL,	CALL DLACPY( 'F', M, BLK, WORK( IL ), LDWRKL,
Line 1110	Line 1206
*	*
ELSE IF( WNTQS ) THEN	ELSE IF( WNTQS ) THEN
*	*
* Path 3t (N much larger than M, JOBZ='S')	* Path 3t (N >> M, JOBZ='S')
* M right singular vectors to be computed in VT and	* M right singular vectors to be computed in VT and
* M left singular vectors to be computed in U	* M left singular vectors to be computed in U
*	*
Line 1123	Line 1219
NWORK = ITAU + M	NWORK = ITAU + M
*	*
* Compute A=L*Q	* Compute A=L*Q
* (Workspace: need MM+2M, prefer MM+M+MNB)	* Workspace: need M*M [L] + M [tau] + M [work]
	* Workspace: prefer MM [L] + M [tau] + MNB [work]
*	*
CALL DGELQF( M, N, A, LDA, WORK( ITAU ), WORK( NWORK ),	CALL DGELQF( M, N, A, LDA, WORK( ITAU ), WORK( NWORK ),
$ LWORK-NWORK+1, IERR )	$ LWORK - NWORK + 1, IERR )
*	*
* Copy L to WORK(IL), zeroing out above it	* Copy L to WORK(IL), zeroing out above it
*	*
CALL DLACPY( 'L', M, M, A, LDA, WORK( IL ), LDWRKL )	CALL DLACPY( 'L', M, M, A, LDA, WORK( IL ), LDWRKL )
CALL DLASET( 'U', M-1, M-1, ZERO, ZERO,	CALL DLASET( 'U', M - 1, M - 1, ZERO, ZERO,
$ WORK( IL+LDWRKL ), LDWRKL )	$ WORK( IL + LDWRKL ), LDWRKL )
*	*
* Generate Q in A	* Generate Q in A
* (Workspace: need MM+2M, prefer MM+M+MNB)	* Workspace: need M*M [L] + M [tau] + M [work]
	* Workspace: prefer MM [L] + M [tau] + MNB [work]
*	*
CALL DORGLQ( M, N, M, A, LDA, WORK( ITAU ),	CALL DORGLQ( M, N, M, A, LDA, WORK( ITAU ),
$ WORK( NWORK ), LWORK-NWORK+1, IERR )	$ WORK( NWORK ), LWORK - NWORK + 1, IERR )
IE = ITAU	IE = ITAU
ITAUQ = IE + M	ITAUQ = IE + M
ITAUP = ITAUQ + M	ITAUP = ITAUQ + M
NWORK = ITAUP + M	NWORK = ITAUP + M
*	*
* Bidiagonalize L in WORK(IU), copying result to U	* Bidiagonalize L in WORK(IU).
* (Workspace: need MM+4M, prefer MM+3M+2MNB)	* Workspace: need MM [L] + 3M [e, tauq, taup] + M [work]
	* Workspace: prefer MM [L] + 3M [e, tauq, taup] + 2MNB [work]
*	*
CALL DGEBRD( M, M, WORK( IL ), LDWRKL, S, WORK( IE ),	CALL DGEBRD( M, M, WORK( IL ), LDWRKL, S, WORK( IE ),
$ WORK( ITAUQ ), WORK( ITAUP ), WORK( NWORK ),	$ WORK( ITAUQ ), WORK( ITAUP ), WORK( NWORK ),
$ LWORK-NWORK+1, IERR )	$ LWORK - NWORK + 1, IERR )
*	*
* Perform bidiagonal SVD, computing left singular vectors	* Perform bidiagonal SVD, computing left singular vectors
* of bidiagonal matrix in U and computing right singular	* of bidiagonal matrix in U and computing right singular
* vectors of bidiagonal matrix in VT	* vectors of bidiagonal matrix in VT
* (Workspace: need M+BDSPAC)	* Workspace: need MM [L] + 3M [e, tauq, taup] + BDSPAC
*	*
CALL DBDSDC( 'U', 'I', M, S, WORK( IE ), U, LDU, VT,	CALL DBDSDC( 'U', 'I', M, S, WORK( IE ), U, LDU, VT,
$ LDVT, DUM, IDUM, WORK( NWORK ), IWORK,	$ LDVT, DUM, IDUM, WORK( NWORK ), IWORK,
Line 1162	Line 1261
*	*
* Overwrite U by left singular vectors of L and VT	* Overwrite U by left singular vectors of L and VT
* by right singular vectors of L	* by right singular vectors of L
* (Workspace: need MM+3M, prefer MM+2M+M*NB)	* Workspace: need MM [L] + 3M [e, tauq, taup] + M [work]
	* Workspace: prefer MM [L] + 3M [e, tauq, taup] + M*NB [work]
*	*
CALL DORMBR( 'Q', 'L', 'N', M, M, M, WORK( IL ), LDWRKL,	CALL DORMBR( 'Q', 'L', 'N', M, M, M, WORK( IL ), LDWRKL,
$ WORK( ITAUQ ), U, LDU, WORK( NWORK ),	$ WORK( ITAUQ ), U, LDU, WORK( NWORK ),
$ LWORK-NWORK+1, IERR )	$ LWORK - NWORK + 1, IERR )
CALL DORMBR( 'P', 'R', 'T', M, M, M, WORK( IL ), LDWRKL,	CALL DORMBR( 'P', 'R', 'T', M, M, M, WORK( IL ), LDWRKL,
$ WORK( ITAUP ), VT, LDVT, WORK( NWORK ),	$ WORK( ITAUP ), VT, LDVT, WORK( NWORK ),
$ LWORK-NWORK+1, IERR )	$ LWORK - NWORK + 1, IERR )
*	*
* Multiply right singular vectors of L in WORK(IL) by	* Multiply right singular vectors of L in WORK(IL) by
* Q in A, storing result in VT	* Q in A, storing result in VT
* (Workspace: need M*M)	* Workspace: need M*M [L]
*	*
CALL DLACPY( 'F', M, M, VT, LDVT, WORK( IL ), LDWRKL )	CALL DLACPY( 'F', M, M, VT, LDVT, WORK( IL ), LDWRKL )
CALL DGEMM( 'N', 'N', M, N, M, ONE, WORK( IL ), LDWRKL,	CALL DGEMM( 'N', 'N', M, N, M, ONE, WORK( IL ), LDWRKL,
Line 1181	Line 1281
*	*
ELSE IF( WNTQA ) THEN	ELSE IF( WNTQA ) THEN
*	*
* Path 4t (N much larger than M, JOBZ='A')	* Path 4t (N >> M, JOBZ='A')
* N right singular vectors to be computed in VT and	* N right singular vectors to be computed in VT and
* M left singular vectors to be computed in U	* M left singular vectors to be computed in U
*	*
Line 1194	Line 1294
NWORK = ITAU + M	NWORK = 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 M*M [VT] + M [tau] + M [work]
	* Workspace: prefer MM [VT] + M [tau] + MNB [work]
*	*
CALL DGELQF( M, N, A, LDA, WORK( ITAU ), WORK( NWORK ),	CALL DGELQF( M, N, A, LDA, WORK( ITAU ), WORK( NWORK ),
$ LWORK-NWORK+1, IERR )	$ LWORK - NWORK + 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+2M, prefer MM+M+MNB)	* Workspace: need M*M [VT] + M [tau] + N [work]
	* Workspace: prefer MM [VT] + M [tau] + NNB [work]
*	*
CALL DORGLQ( N, N, M, VT, LDVT, WORK( ITAU ),	CALL DORGLQ( N, N, M, VT, LDVT, WORK( ITAU ),
$ WORK( NWORK ), LWORK-NWORK+1, IERR )	$ WORK( NWORK ), LWORK - NWORK + 1, IERR )
*	*
* Produce L in A, zeroing out other entries	* Produce L in A, zeroing out other entries
*	*
Line 1215	Line 1317
NWORK = ITAUP + M	NWORK = ITAUP + M
*	*
* Bidiagonalize L in A	* Bidiagonalize L in A
* (Workspace: need MM+4M, prefer MM+3M+2MNB)	* Workspace: need MM [VT] + 3M [e, tauq, taup] + M [work]
	* Workspace: prefer MM [VT] + 3M [e, tauq, taup] + 2MNB [work]
*	*
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( NWORK ), LWORK-NWORK+1,	$ WORK( ITAUP ), WORK( NWORK ), LWORK-NWORK+1,
Line 1224	Line 1327
* Perform bidiagonal SVD, computing left singular vectors	* Perform bidiagonal SVD, computing left singular vectors
* of bidiagonal matrix in U and computing right singular	* of bidiagonal matrix in U and computing right singular
* vectors of bidiagonal matrix in WORK(IVT)	* vectors of bidiagonal matrix in WORK(IVT)
* (Workspace: need M+M*M+BDSPAC)	* Workspace: need MM [VT] + 3M [e, tauq, taup] + BDSPAC
*	*
CALL DBDSDC( 'U', 'I', M, S, WORK( IE ), U, LDU,	CALL DBDSDC( 'U', 'I', M, S, WORK( IE ), U, LDU,
$ WORK( IVT ), LDWKVT, DUM, IDUM,	$ WORK( IVT ), LDWKVT, DUM, IDUM,
Line 1232	Line 1335
*	*
* Overwrite U by left singular vectors of L and WORK(IVT)	* Overwrite U by left singular vectors of L and WORK(IVT)
* by right singular vectors of L	* by right singular vectors of L
* (Workspace: need MM+3M, prefer MM+2M+M*NB)	* Workspace: need MM [VT] + 3M [e, tauq, taup]+ M [work]
	* Workspace: prefer MM [VT] + 3M [e, tauq, taup]+ M*NB [work]
*	*
CALL DORMBR( 'Q', 'L', 'N', M, M, M, A, LDA,	CALL DORMBR( 'Q', 'L', 'N', M, M, M, A, LDA,
$ WORK( ITAUQ ), U, LDU, WORK( NWORK ),	$ WORK( ITAUQ ), U, LDU, WORK( NWORK ),
$ LWORK-NWORK+1, IERR )	$ LWORK - NWORK + 1, IERR )
CALL DORMBR( 'P', 'R', 'T', M, M, M, A, LDA,	CALL DORMBR( 'P', 'R', 'T', M, M, M, A, LDA,
$ WORK( ITAUP ), WORK( IVT ), LDWKVT,	$ WORK( ITAUP ), WORK( IVT ), LDWKVT,
$ WORK( NWORK ), LWORK-NWORK+1, IERR )	$ WORK( NWORK ), LWORK - NWORK + 1, IERR )
*	*
* Multiply right singular vectors of L in WORK(IVT) by	* Multiply right singular vectors of L in WORK(IVT) by
* Q in VT, storing result in A	* Q in VT, storing result in A
* (Workspace: need M*M)	* Workspace: need M*M [VT]
*	*
CALL DGEMM( 'N', 'N', M, N, M, ONE, WORK( IVT ), LDWKVT,	CALL DGEMM( 'N', 'N', M, N, M, ONE, WORK( IVT ), LDWKVT,
$ VT, LDVT, ZERO, A, LDA )	$ VT, LDVT, ZERO, A, LDA )
Line 1258	Line 1362
*	*
* N .LT. MNTHR	* N .LT. MNTHR
*	*
* Path 5t (N greater than M, but not much larger)	* Path 5t (N > M, but not much larger)
* Reduce to bidiagonal form without LQ decomposition	* Reduce to bidiagonal form without LQ decomposition
*	*
IE = 1	IE = 1
Line 1267	Line 1371
NWORK = ITAUP + M	NWORK = ITAUP + M
*	*
* Bidiagonalize A	* Bidiagonalize A
* (Workspace: need 3M+N, prefer 3M+(M+N)*NB)	* Workspace: need 3*M [e, tauq, taup] + N [work]
	* Workspace: prefer 3M [e, tauq, taup] + (M+N)NB [work]
*	*
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( NWORK ), LWORK-NWORK+1,	$ WORK( ITAUP ), WORK( NWORK ), LWORK-NWORK+1,
$ IERR )	$ IERR )
IF( WNTQN ) THEN	IF( WNTQN ) THEN
*	*
	* Path 5tn (N > M, JOBZ='N')
* Perform bidiagonal SVD, only computing singular values	* Perform bidiagonal SVD, only computing singular values
* (Workspace: need M+BDSPAC)	* Workspace: need 3*M [e, tauq, taup] + BDSPAC
*	*
CALL DBDSDC( 'L', 'N', M, S, WORK( IE ), DUM, 1, DUM, 1,	CALL DBDSDC( 'L', 'N', M, S, WORK( IE ), DUM, 1, DUM, 1,
$ DUM, IDUM, WORK( NWORK ), IWORK, INFO )	$ DUM, IDUM, WORK( NWORK ), IWORK, INFO )
ELSE IF( WNTQO ) THEN	ELSE IF( WNTQO ) THEN
	* Path 5to (N > M, JOBZ='O')
LDWKVT = M	LDWKVT = M
IVT = NWORK	IVT = NWORK
IF( LWORK.GE.MN+3M+BDSPAC ) THEN	IF( LWORK .GE. MN + 3M + BDSPAC ) THEN
*	*
* WORK( IVT ) is M by N	* WORK( IVT ) is M by N
*	*
CALL DLASET( 'F', M, N, ZERO, ZERO, WORK( IVT ),	CALL DLASET( 'F', M, N, ZERO, ZERO, WORK( IVT ),
$ LDWKVT )	$ LDWKVT )
NWORK = IVT + LDWKVT*N	NWORK = IVT + LDWKVT*N
	* IL is unused; silence compile warnings
	IL = -1
ELSE	ELSE
*	*
* WORK( IVT ) is M by M	* WORK( IVT ) is M by M
Line 1298	Line 1407
*	*
* WORK(IL) is M by CHUNK	* WORK(IL) is M by CHUNK
*	*
CHUNK = ( LWORK-MM-3M ) / M	CHUNK = ( LWORK - MM - 3M ) / M
END IF	END IF
*	*
* Perform bidiagonal SVD, computing left singular vectors	* Perform bidiagonal SVD, computing left singular vectors
* of bidiagonal matrix in U and computing right singular	* of bidiagonal matrix in U and computing right singular
* vectors of bidiagonal matrix in WORK(IVT)	* vectors of bidiagonal matrix in WORK(IVT)
* (Workspace: need M*M+BDSPAC)	* Workspace: need 3M [e, tauq, taup] + MM [VT] + BDSPAC
*	*
CALL DBDSDC( 'L', 'I', M, S, WORK( IE ), U, LDU,	CALL DBDSDC( 'L', 'I', M, S, WORK( IE ), U, LDU,
$ WORK( IVT ), LDWKVT, DUM, IDUM,	$ WORK( IVT ), LDWKVT, DUM, IDUM,
$ WORK( NWORK ), IWORK, INFO )	$ WORK( NWORK ), IWORK, INFO )
*	*
* Overwrite U by left singular vectors of A	* Overwrite U by left singular vectors of A
* (Workspace: need MM+2M, prefer MM+M+MNB)	* Workspace: need 3M [e, tauq, taup] + MM [VT] + M [work]
	* Workspace: prefer 3M [e, tauq, taup] + MM [VT] + M*NB [work]
*	*
CALL DORMBR( 'Q', 'L', 'N', M, M, N, A, LDA,	CALL DORMBR( 'Q', 'L', 'N', M, M, N, A, LDA,
$ WORK( ITAUQ ), U, LDU, WORK( NWORK ),	$ WORK( ITAUQ ), U, LDU, WORK( NWORK ),
$ LWORK-NWORK+1, IERR )	$ LWORK - NWORK + 1, IERR )
*	*
IF( LWORK.GE.MN+3M+BDSPAC ) THEN	IF( LWORK .GE. MN + 3M + BDSPAC ) THEN
*	*
	* Path 5to-fast
* Overwrite WORK(IVT) by left singular vectors of A	* Overwrite WORK(IVT) by left singular vectors of A
* (Workspace: need MM+2M, prefer MM+M+MNB)	* Workspace: need 3M [e, tauq, taup] + MN [VT] + M [work]
	* Workspace: prefer 3M [e, tauq, taup] + MN [VT] + M*NB [work]
*	*
CALL DORMBR( 'P', 'R', 'T', M, N, M, A, LDA,	CALL DORMBR( 'P', 'R', 'T', M, N, M, A, LDA,
$ WORK( ITAUP ), WORK( IVT ), LDWKVT,	$ WORK( ITAUP ), WORK( IVT ), LDWKVT,
$ WORK( NWORK ), LWORK-NWORK+1, IERR )	$ WORK( NWORK ), LWORK - NWORK + 1, IERR )
*	*
* Copy right singular vectors of A from WORK(IVT) to A	* Copy right singular vectors of A from WORK(IVT) to A
*	*
CALL DLACPY( 'F', M, N, WORK( IVT ), LDWKVT, A, LDA )	CALL DLACPY( 'F', M, N, WORK( IVT ), LDWKVT, A, LDA )
ELSE	ELSE
*	*
	* Path 5to-slow
* Generate P**T in A	* Generate P**T in A
* (Workspace: need MM+2M, prefer MM+M+MNB)	* Workspace: need 3M [e, tauq, taup] + MM [VT] + M [work]
	* Workspace: prefer 3M [e, tauq, taup] + MM [VT] + M*NB [work]
*	*
CALL DORGBR( 'P', M, N, M, A, LDA, WORK( ITAUP ),	CALL DORGBR( 'P', M, N, M, A, LDA, WORK( ITAUP ),
$ WORK( NWORK ), LWORK-NWORK+1, IERR )	$ WORK( NWORK ), LWORK - NWORK + 1, IERR )
*	*
* Multiply Q in A by right singular vectors of	* Multiply Q in A by right singular vectors of
* bidiagonal matrix in WORK(IVT), storing result in	* bidiagonal matrix in WORK(IVT), storing result in
* WORK(IL) and copying to A	* WORK(IL) and copying to A
* (Workspace: need 2MM, prefer MM+MN)	* Workspace: need 3M [e, tauq, taup] + MM [VT] + M*NB [L]
	* Workspace: prefer 3M [e, tauq, taup] + MM [VT] + M*N [L]
*	*
DO 40 I = 1, N, CHUNK	DO 40 I = 1, N, CHUNK
BLK = MIN( N-I+1, CHUNK )	BLK = MIN( N - I + 1, CHUNK )
CALL DGEMM( 'N', 'N', M, BLK, M, ONE, WORK( IVT ),	CALL DGEMM( 'N', 'N', M, BLK, M, ONE, WORK( IVT ),
$ LDWKVT, A( 1, I ), LDA, ZERO,	$ LDWKVT, A( 1, I ), LDA, ZERO,
$ WORK( IL ), M )	$ WORK( IL ), M )
Line 1353	Line 1468
END IF	END IF
ELSE IF( WNTQS ) THEN	ELSE IF( WNTQS ) THEN
*	*
	* Path 5ts (N > M, JOBZ='S')
* Perform bidiagonal SVD, computing left singular vectors	* Perform bidiagonal SVD, computing left singular vectors
* of bidiagonal matrix in U and computing right singular	* of bidiagonal matrix in U and computing right singular
* vectors of bidiagonal matrix in VT	* vectors of bidiagonal matrix in VT
* (Workspace: need M+BDSPAC)	* Workspace: need 3*M [e, tauq, taup] + BDSPAC
*	*
CALL DLASET( 'F', M, N, ZERO, ZERO, VT, LDVT )	CALL DLASET( 'F', M, N, ZERO, ZERO, VT, LDVT )
CALL DBDSDC( 'L', 'I', M, S, WORK( IE ), U, LDU, VT,	CALL DBDSDC( 'L', 'I', M, S, WORK( IE ), U, LDU, VT,
Line 1365	Line 1481
*	*
* Overwrite U by left singular vectors of A and VT	* Overwrite U by left singular vectors of A and VT
* by right singular vectors of A	* by right singular vectors of A
* (Workspace: need 3M, prefer 2M+M*NB)	* Workspace: need 3*M [e, tauq, taup] + M [work]
	* Workspace: prefer 3M [e, tauq, taup] + MNB [work]
*	*
CALL DORMBR( 'Q', 'L', 'N', M, M, N, A, LDA,	CALL DORMBR( 'Q', 'L', 'N', M, M, N, A, LDA,
$ WORK( ITAUQ ), U, LDU, WORK( NWORK ),	$ WORK( ITAUQ ), U, LDU, WORK( NWORK ),
$ LWORK-NWORK+1, IERR )	$ LWORK - NWORK + 1, IERR )
CALL DORMBR( 'P', 'R', 'T', M, N, M, A, LDA,	CALL DORMBR( 'P', 'R', 'T', M, N, M, A, LDA,
$ WORK( ITAUP ), VT, LDVT, WORK( NWORK ),	$ WORK( ITAUP ), VT, LDVT, WORK( NWORK ),
$ LWORK-NWORK+1, IERR )	$ LWORK - NWORK + 1, IERR )
ELSE IF( WNTQA ) THEN	ELSE IF( WNTQA ) THEN
*	*
	* Path 5ta (N > M, JOBZ='A')
* Perform bidiagonal SVD, computing left singular vectors	* Perform bidiagonal SVD, computing left singular vectors
* of bidiagonal matrix in U and computing right singular	* of bidiagonal matrix in U and computing right singular
* vectors of bidiagonal matrix in VT	* vectors of bidiagonal matrix in VT
* (Workspace: need M+BDSPAC)	* Workspace: need 3*M [e, tauq, taup] + BDSPAC
*	*
CALL DLASET( 'F', N, N, ZERO, ZERO, VT, LDVT )	CALL DLASET( 'F', N, N, ZERO, ZERO, VT, LDVT )
CALL DBDSDC( 'L', 'I', M, S, WORK( IE ), U, LDU, VT,	CALL DBDSDC( 'L', 'I', M, S, WORK( IE ), U, LDU, VT,
Line 1388	Line 1506
* Set the right corner of VT to identity matrix	* Set the right corner of VT to identity matrix
*	*
IF( N.GT.M ) THEN	IF( N.GT.M ) THEN
CALL DLASET( 'F', N-M, N-M, ZERO, ONE, VT( M+1, M+1 ),	CALL DLASET( 'F', N-M, N-M, ZERO, ONE, VT(M+1,M+1),
$ LDVT )	$ LDVT )
END IF	END IF
*	*
* Overwrite U by left singular vectors of A and VT	* Overwrite U by left singular vectors of A and VT
* by right singular vectors of A	* by right singular vectors of A
* (Workspace: need 2M+N, prefer 2M+N*NB)	* Workspace: need 3*M [e, tauq, taup] + N [work]
	* Workspace: prefer 3M [e, tauq, taup] + NNB [work]
*	*
CALL DORMBR( 'Q', 'L', 'N', M, M, N, A, LDA,	CALL DORMBR( 'Q', 'L', 'N', M, M, N, A, LDA,
$ WORK( ITAUQ ), U, LDU, WORK( NWORK ),	$ WORK( ITAUQ ), U, LDU, WORK( NWORK ),
$ LWORK-NWORK+1, IERR )	$ LWORK - NWORK + 1, IERR )
CALL DORMBR( 'P', 'R', 'T', N, N, M, A, LDA,	CALL DORMBR( 'P', 'R', 'T', N, N, M, A, LDA,
$ WORK( ITAUP ), VT, LDVT, WORK( NWORK ),	$ WORK( ITAUP ), VT, LDVT, WORK( NWORK ),
$ LWORK-NWORK+1, IERR )	$ LWORK - NWORK + 1, IERR )
END IF	END IF
*	*
END IF	END IF

CVSweb interface <joel.bertrand@systella.fr>

Removed from v.1.15
changed lines
	Added in v.1.16