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.14 and 1.15

version 1.14, 2014/01/27 09:28:17	version 1.15, 2016/08/27 15:27:08
Line 175	Line 175
> LWORK >= MAX(1,5MIN(M,N)) for the paths (see comments inside code):	> LWORK >= MAX(1,5MIN(M,N)) for the paths (see comments inside code):
*> - PATH 1 (M much larger than N, JOBU='N')	*> - PATH 1 (M much larger than N, JOBU='N')
*> - PATH 1t (N much larger than M, JOBVT='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	> 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.	*> For good performance, LWORK should generally be larger.
*>	*>
*> If LWORK = -1, then a workspace query is assumed; the routine	*> If LWORK = -1, then a workspace query is assumed; the routine
Line 211	Line 211
SUBROUTINE DGESVD( JOBU, JOBVT, M, N, A, LDA, S, U, LDU,	SUBROUTINE DGESVD( JOBU, JOBVT, M, N, A, LDA, S, U, LDU,
$ VT, LDVT, WORK, LWORK, INFO )	$ VT, LDVT, WORK, LWORK, INFO )
*	*
* -- LAPACK driver routine (version 3.4.1) --	* -- 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..--
* April 2012	* April 2012
Line 314	Line 314
BDSPAC = 5*N	BDSPAC = 5*N
* Compute space needed for DGEQRF	* Compute space needed for DGEQRF
CALL DGEQRF( M, N, A, LDA, DUM(1), DUM(1), -1, IERR )	CALL DGEQRF( M, N, A, LDA, DUM(1), DUM(1), -1, IERR )
LWORK_DGEQRF=DUM(1)	LWORK_DGEQRF = INT( DUM(1) )
* Compute space needed for DORGQR	* Compute space needed for DORGQR
CALL DORGQR( M, N, N, A, LDA, DUM(1), DUM(1), -1, IERR )	CALL DORGQR( M, N, N, A, LDA, DUM(1), DUM(1), -1, IERR )
LWORK_DORGQR_N=DUM(1)	LWORK_DORGQR_N = INT( DUM(1) )
CALL DORGQR( M, M, N, A, LDA, DUM(1), DUM(1), -1, IERR )	CALL DORGQR( M, M, N, A, LDA, DUM(1), DUM(1), -1, IERR )
LWORK_DORGQR_M=DUM(1)	LWORK_DORGQR_M = INT( DUM(1) )
* Compute space needed for DGEBRD	* Compute space needed for DGEBRD
CALL DGEBRD( N, N, A, LDA, S, DUM(1), DUM(1),	CALL DGEBRD( N, N, A, LDA, S, DUM(1), DUM(1),
$ DUM(1), DUM(1), -1, IERR )	$ DUM(1), DUM(1), -1, IERR )
LWORK_DGEBRD=DUM(1)	LWORK_DGEBRD = INT( DUM(1) )
* Compute space needed for DORGBR P	* Compute space needed for DORGBR P
CALL DORGBR( 'P', N, N, N, A, LDA, DUM(1),	CALL DORGBR( 'P', N, N, N, A, LDA, DUM(1),
$ DUM(1), -1, IERR )	$ DUM(1), -1, IERR )
LWORK_DORGBR_P=DUM(1)	LWORK_DORGBR_P = INT( DUM(1) )
* Compute space needed for DORGBR Q	* Compute space needed for DORGBR Q
CALL DORGBR( 'Q', N, N, N, A, LDA, DUM(1),	CALL DORGBR( 'Q', N, N, N, A, LDA, DUM(1),
$ DUM(1), -1, IERR )	$ DUM(1), -1, IERR )
LWORK_DORGBR_Q=DUM(1)	LWORK_DORGBR_Q = INT( DUM(1) )
*	*
IF( M.GE.MNTHR ) THEN	IF( M.GE.MNTHR ) THEN
IF( WNTUN ) THEN	IF( WNTUN ) THEN
Line 339	Line 339
* Path 1 (M much larger than N, JOBU='N')	* Path 1 (M much larger than N, JOBU='N')
*	*
MAXWRK = N + LWORK_DGEQRF	MAXWRK = N + LWORK_DGEQRF
MAXWRK = MAX( MAXWRK, 3*N+LWORK_DGEBRD )	MAXWRK = MAX( MAXWRK, 3*N + LWORK_DGEBRD )
IF( WNTVO .OR. WNTVAS )	IF( WNTVO .OR. WNTVAS )
$ MAXWRK = MAX( MAXWRK, 3*N+LWORK_DORGBR_P )	$ MAXWRK = MAX( MAXWRK, 3*N + LWORK_DORGBR_P )
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
Line 349	Line 349
* Path 2 (M much larger than N, JOBU='O', JOBVT='N')	* Path 2 (M much larger than N, JOBU='O', JOBVT='N')
*	*
WRKBL = N + LWORK_DGEQRF	WRKBL = N + LWORK_DGEQRF
WRKBL = MAX( WRKBL, N+LWORK_DORGQR_N )	WRKBL = MAX( WRKBL, N + LWORK_DORGQR_N )
WRKBL = MAX( WRKBL, 3*N+LWORK_DGEBRD )	WRKBL = MAX( WRKBL, 3*N + LWORK_DGEBRD )
WRKBL = MAX( WRKBL, 3*N+LWORK_DORGBR_Q )	WRKBL = MAX( WRKBL, 3*N + LWORK_DORGBR_Q )
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 + LWORK_DGEQRF	WRKBL = N + LWORK_DGEQRF
WRKBL = MAX( WRKBL, N+LWORK_DORGQR_N )	WRKBL = MAX( WRKBL, N + LWORK_DORGQR_N )
WRKBL = MAX( WRKBL, 3*N+LWORK_DGEBRD )	WRKBL = MAX( WRKBL, 3*N + LWORK_DGEBRD )
WRKBL = MAX( WRKBL, 3*N+LWORK_DORGBR_Q )	WRKBL = MAX( WRKBL, 3*N + LWORK_DORGBR_Q )
WRKBL = MAX( WRKBL, 3*N+LWORK_DORGBR_P )	WRKBL = MAX( WRKBL, 3*N + LWORK_DORGBR_P )
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 + LWORK_DGEQRF	WRKBL = N + LWORK_DGEQRF
WRKBL = MAX( WRKBL, N+LWORK_DORGQR_N )	WRKBL = MAX( WRKBL, N + LWORK_DORGQR_N )
WRKBL = MAX( WRKBL, 3*N+LWORK_DGEBRD )	WRKBL = MAX( WRKBL, 3*N + LWORK_DGEBRD )
WRKBL = MAX( WRKBL, 3*N+LWORK_DORGBR_Q )	WRKBL = MAX( WRKBL, 3*N + LWORK_DORGBR_Q )
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 + LWORK_DGEQRF	WRKBL = N + LWORK_DGEQRF
WRKBL = MAX( WRKBL, N+LWORK_DORGQR_N )	WRKBL = MAX( WRKBL, N + LWORK_DORGQR_N )
WRKBL = MAX( WRKBL, 3*N+LWORK_DGEBRD )	WRKBL = MAX( WRKBL, 3*N + LWORK_DGEBRD )
WRKBL = MAX( WRKBL, 3*N+LWORK_DORGBR_Q )	WRKBL = MAX( WRKBL, 3*N + LWORK_DORGBR_Q )
WRKBL = MAX( WRKBL, 3*N+LWORK_DORGBR_P )	WRKBL = MAX( WRKBL, 3*N + LWORK_DORGBR_P )
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 + LWORK_DGEQRF	WRKBL = N + LWORK_DGEQRF
WRKBL = MAX( WRKBL, N+LWORK_DORGQR_N )	WRKBL = MAX( WRKBL, N + LWORK_DORGQR_N )
WRKBL = MAX( WRKBL, 3*N+LWORK_DGEBRD )	WRKBL = MAX( WRKBL, 3*N + LWORK_DGEBRD )
WRKBL = MAX( WRKBL, 3*N+LWORK_DORGBR_Q )	WRKBL = MAX( WRKBL, 3*N + LWORK_DORGBR_Q )
WRKBL = MAX( WRKBL, 3*N+LWORK_DORGBR_P )	WRKBL = MAX( WRKBL, 3*N + LWORK_DORGBR_P )
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 + LWORK_DGEQRF	WRKBL = N + LWORK_DGEQRF
WRKBL = MAX( WRKBL, N+LWORK_DORGQR_M )	WRKBL = MAX( WRKBL, N + LWORK_DORGQR_M )
WRKBL = MAX( WRKBL, 3*N+LWORK_DGEBRD )	WRKBL = MAX( WRKBL, 3*N + LWORK_DGEBRD )
WRKBL = MAX( WRKBL, 3*N+LWORK_DORGBR_Q )	WRKBL = MAX( WRKBL, 3*N + LWORK_DORGBR_Q )
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 + LWORK_DGEQRF	WRKBL = N + LWORK_DGEQRF
WRKBL = MAX( WRKBL, N+LWORK_DORGQR_M )	WRKBL = MAX( WRKBL, N + LWORK_DORGQR_M )
WRKBL = MAX( WRKBL, 3*N+LWORK_DGEBRD )	WRKBL = MAX( WRKBL, 3*N + LWORK_DGEBRD )
WRKBL = MAX( WRKBL, 3*N+LWORK_DORGBR_Q )	WRKBL = MAX( WRKBL, 3*N + LWORK_DORGBR_Q )
WRKBL = MAX( WRKBL, 3*N+LWORK_DORGBR_P )	WRKBL = MAX( WRKBL, 3*N + LWORK_DORGBR_P )
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 + LWORK_DGEQRF	WRKBL = N + LWORK_DGEQRF
WRKBL = MAX( WRKBL, N+LWORK_DORGQR_M )	WRKBL = MAX( WRKBL, N + LWORK_DORGQR_M )
WRKBL = MAX( WRKBL, 3*N+LWORK_DGEBRD )	WRKBL = MAX( WRKBL, 3*N + LWORK_DGEBRD )
WRKBL = MAX( WRKBL, 3*N+LWORK_DORGBR_Q )	WRKBL = MAX( WRKBL, 3*N + LWORK_DORGBR_Q )
WRKBL = MAX( WRKBL, 3*N+LWORK_DORGBR_P )	WRKBL = MAX( WRKBL, 3*N + LWORK_DORGBR_P )
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
*	*
Line 447	Line 447
*	*
CALL DGEBRD( M, N, A, LDA, S, DUM(1), DUM(1),	CALL DGEBRD( M, N, A, LDA, S, DUM(1), DUM(1),
$ DUM(1), DUM(1), -1, IERR )	$ DUM(1), DUM(1), -1, IERR )
LWORK_DGEBRD=DUM(1)	LWORK_DGEBRD = INT( DUM(1) )
MAXWRK = 3*N + LWORK_DGEBRD	MAXWRK = 3*N + LWORK_DGEBRD
IF( WNTUS .OR. WNTUO ) THEN	IF( WNTUS .OR. WNTUO ) THEN
CALL DORGBR( 'Q', M, N, N, A, LDA, DUM(1),	CALL DORGBR( 'Q', M, N, N, A, LDA, DUM(1),
$ DUM(1), -1, IERR )	$ DUM(1), -1, IERR )
LWORK_DORGBR_Q=DUM(1)	LWORK_DORGBR_Q = INT( DUM(1) )
MAXWRK = MAX( MAXWRK, 3*N+LWORK_DORGBR_Q )	MAXWRK = MAX( MAXWRK, 3*N + LWORK_DORGBR_Q )
END IF	END IF
IF( WNTUA ) THEN	IF( WNTUA ) THEN
CALL DORGBR( 'Q', M, M, N, A, LDA, DUM(1),	CALL DORGBR( 'Q', M, M, N, A, LDA, DUM(1),
$ DUM(1), -1, IERR )	$ DUM(1), -1, IERR )
LWORK_DORGBR_Q=DUM(1)	LWORK_DORGBR_Q = INT( DUM(1) )
MAXWRK = MAX( MAXWRK, 3*N+LWORK_DORGBR_Q )	MAXWRK = MAX( MAXWRK, 3*N + LWORK_DORGBR_Q )
END IF	END IF
IF( .NOT.WNTVN ) THEN	IF( .NOT.WNTVN ) THEN
MAXWRK = MAX( MAXWRK, 3*N+LWORK_DORGBR_P )	MAXWRK = MAX( MAXWRK, 3*N + LWORK_DORGBR_P )
END IF	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 475	Line 475
BDSPAC = 5*M	BDSPAC = 5*M
* Compute space needed for DGELQF	* Compute space needed for DGELQF
CALL DGELQF( M, N, A, LDA, DUM(1), DUM(1), -1, IERR )	CALL DGELQF( M, N, A, LDA, DUM(1), DUM(1), -1, IERR )
LWORK_DGELQF=DUM(1)	LWORK_DGELQF = INT( DUM(1) )
* Compute space needed for DORGLQ	* Compute space needed for DORGLQ
CALL DORGLQ( N, N, M, DUM(1), N, DUM(1), DUM(1), -1, IERR )	CALL DORGLQ( N, N, M, DUM(1), N, DUM(1), DUM(1), -1, IERR )
LWORK_DORGLQ_N=DUM(1)	LWORK_DORGLQ_N = INT( DUM(1) )
CALL DORGLQ( M, N, M, A, LDA, DUM(1), DUM(1), -1, IERR )	CALL DORGLQ( M, N, M, A, LDA, DUM(1), DUM(1), -1, IERR )
LWORK_DORGLQ_M=DUM(1)	LWORK_DORGLQ_M = INT( DUM(1) )
* Compute space needed for DGEBRD	* Compute space needed for DGEBRD
CALL DGEBRD( M, M, A, LDA, S, DUM(1), DUM(1),	CALL DGEBRD( M, M, A, LDA, S, DUM(1), DUM(1),
$ DUM(1), DUM(1), -1, IERR )	$ DUM(1), DUM(1), -1, IERR )
LWORK_DGEBRD=DUM(1)	LWORK_DGEBRD = INT( DUM(1) )
* Compute space needed for DORGBR P	* Compute space needed for DORGBR P
CALL DORGBR( 'P', M, M, M, A, N, DUM(1),	CALL DORGBR( 'P', M, M, M, A, N, DUM(1),
$ DUM(1), -1, IERR )	$ DUM(1), -1, IERR )
LWORK_DORGBR_P=DUM(1)	LWORK_DORGBR_P = INT( DUM(1) )
* Compute space needed for DORGBR Q	* Compute space needed for DORGBR Q
CALL DORGBR( 'Q', M, M, M, A, N, DUM(1),	CALL DORGBR( 'Q', M, M, M, A, N, DUM(1),
$ DUM(1), -1, IERR )	$ DUM(1), -1, IERR )
LWORK_DORGBR_Q=DUM(1)	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 + LWORK_DGELQF	MAXWRK = M + LWORK_DGELQF
MAXWRK = MAX( MAXWRK, 3*M+LWORK_DGEBRD )	MAXWRK = MAX( MAXWRK, 3*M + LWORK_DGEBRD )
IF( WNTUO .OR. WNTUAS )	IF( WNTUO .OR. WNTUAS )
$ MAXWRK = MAX( MAXWRK, 3*M+LWORK_DORGBR_Q )	$ MAXWRK = MAX( MAXWRK, 3*M + LWORK_DORGBR_Q )
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
Line 509	Line 509
* Path 2t(N much larger than M, JOBU='N', JOBVT='O')	* Path 2t(N much larger than M, JOBU='N', JOBVT='O')
*	*
WRKBL = M + LWORK_DGELQF	WRKBL = M + LWORK_DGELQF
WRKBL = MAX( WRKBL, M+LWORK_DORGLQ_M )	WRKBL = MAX( WRKBL, M + LWORK_DORGLQ_M )
WRKBL = MAX( WRKBL, 3*M+LWORK_DGEBRD )	WRKBL = MAX( WRKBL, 3*M + LWORK_DGEBRD )
WRKBL = MAX( WRKBL, 3*M+LWORK_DORGBR_P )	WRKBL = MAX( WRKBL, 3*M + LWORK_DORGBR_P )
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 + LWORK_DGELQF	WRKBL = M + LWORK_DGELQF
WRKBL = MAX( WRKBL, M+LWORK_DORGLQ_M )	WRKBL = MAX( WRKBL, M + LWORK_DORGLQ_M )
WRKBL = MAX( WRKBL, 3*M+LWORK_DGEBRD )	WRKBL = MAX( WRKBL, 3*M + LWORK_DGEBRD )
WRKBL = MAX( WRKBL, 3*M+LWORK_DORGBR_P )	WRKBL = MAX( WRKBL, 3*M + LWORK_DORGBR_P )
WRKBL = MAX( WRKBL, 3*M+LWORK_DORGBR_Q )	WRKBL = MAX( WRKBL, 3*M + LWORK_DORGBR_Q )
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 + LWORK_DGELQF	WRKBL = M + LWORK_DGELQF
WRKBL = MAX( WRKBL, M+LWORK_DORGLQ_M )	WRKBL = MAX( WRKBL, M + LWORK_DORGLQ_M )
WRKBL = MAX( WRKBL, 3*M+LWORK_DGEBRD )	WRKBL = MAX( WRKBL, 3*M + LWORK_DGEBRD )
WRKBL = MAX( WRKBL, 3*M+LWORK_DORGBR_P )	WRKBL = MAX( WRKBL, 3*M + LWORK_DORGBR_P )
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 + LWORK_DGELQF	WRKBL = M + LWORK_DGELQF
WRKBL = MAX( WRKBL, M+LWORK_DORGLQ_M )	WRKBL = MAX( WRKBL, M + LWORK_DORGLQ_M )
WRKBL = MAX( WRKBL, 3*M+LWORK_DGEBRD )	WRKBL = MAX( WRKBL, 3*M + LWORK_DGEBRD )
WRKBL = MAX( WRKBL, 3*M+LWORK_DORGBR_P )	WRKBL = MAX( WRKBL, 3*M + LWORK_DORGBR_P )
WRKBL = MAX( WRKBL, 3*M+LWORK_DORGBR_Q )	WRKBL = MAX( WRKBL, 3*M + LWORK_DORGBR_Q )
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 + LWORK_DGELQF	WRKBL = M + LWORK_DGELQF
WRKBL = MAX( WRKBL, M+LWORK_DORGLQ_M )	WRKBL = MAX( WRKBL, M + LWORK_DORGLQ_M )
WRKBL = MAX( WRKBL, 3*M+LWORK_DGEBRD )	WRKBL = MAX( WRKBL, 3*M + LWORK_DGEBRD )
WRKBL = MAX( WRKBL, 3*M+LWORK_DORGBR_P )	WRKBL = MAX( WRKBL, 3*M + LWORK_DORGBR_P )
WRKBL = MAX( WRKBL, 3*M+LWORK_DORGBR_Q )	WRKBL = MAX( WRKBL, 3*M + LWORK_DORGBR_Q )
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 + LWORK_DGELQF	WRKBL = M + LWORK_DGELQF
WRKBL = MAX( WRKBL, M+LWORK_DORGLQ_N )	WRKBL = MAX( WRKBL, M + LWORK_DORGLQ_N )
WRKBL = MAX( WRKBL, 3*M+LWORK_DGEBRD )	WRKBL = MAX( WRKBL, 3*M + LWORK_DGEBRD )
WRKBL = MAX( WRKBL, 3*M+LWORK_DORGBR_P )	WRKBL = MAX( WRKBL, 3*M + LWORK_DORGBR_P )
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 + LWORK_DGELQF	WRKBL = M + LWORK_DGELQF
WRKBL = MAX( WRKBL, M+LWORK_DORGLQ_N )	WRKBL = MAX( WRKBL, M + LWORK_DORGLQ_N )
WRKBL = MAX( WRKBL, 3*M+LWORK_DGEBRD )	WRKBL = MAX( WRKBL, 3*M + LWORK_DGEBRD )
WRKBL = MAX( WRKBL, 3*M+LWORK_DORGBR_P )	WRKBL = MAX( WRKBL, 3*M + LWORK_DORGBR_P )
WRKBL = MAX( WRKBL, 3*M+LWORK_DORGBR_Q )	WRKBL = MAX( WRKBL, 3*M + LWORK_DORGBR_Q )
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 + LWORK_DGELQF	WRKBL = M + LWORK_DGELQF
WRKBL = MAX( WRKBL, M+LWORK_DORGLQ_N )	WRKBL = MAX( WRKBL, M + LWORK_DORGLQ_N )
WRKBL = MAX( WRKBL, 3*M+LWORK_DGEBRD )	WRKBL = MAX( WRKBL, 3*M + LWORK_DGEBRD )
WRKBL = MAX( WRKBL, 3*M+LWORK_DORGBR_P )	WRKBL = MAX( WRKBL, 3*M + LWORK_DORGBR_P )
WRKBL = MAX( WRKBL, 3*M+LWORK_DORGBR_Q )	WRKBL = MAX( WRKBL, 3*M + LWORK_DORGBR_Q )
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
*	*
Line 607	Line 607
*	*
CALL DGEBRD( M, N, A, LDA, S, DUM(1), DUM(1),	CALL DGEBRD( M, N, A, LDA, S, DUM(1), DUM(1),
$ DUM(1), DUM(1), -1, IERR )	$ DUM(1), DUM(1), -1, IERR )
LWORK_DGEBRD=DUM(1)	LWORK_DGEBRD = INT( DUM(1) )
MAXWRK = 3*M + LWORK_DGEBRD	MAXWRK = 3*M + LWORK_DGEBRD
IF( WNTVS .OR. WNTVO ) THEN	IF( WNTVS .OR. WNTVO ) THEN
* Compute space needed for DORGBR P	* Compute space needed for DORGBR P
CALL DORGBR( 'P', M, N, M, A, N, DUM(1),	CALL DORGBR( 'P', M, N, M, A, N, DUM(1),
$ DUM(1), -1, IERR )	$ DUM(1), -1, IERR )
LWORK_DORGBR_P=DUM(1)	LWORK_DORGBR_P = INT( DUM(1) )
MAXWRK = MAX( MAXWRK, 3*M+LWORK_DORGBR_P )	MAXWRK = MAX( MAXWRK, 3*M + LWORK_DORGBR_P )
END IF	END IF
IF( WNTVA ) THEN	IF( WNTVA ) THEN
CALL DORGBR( 'P', N, N, M, A, N, DUM(1),	CALL DORGBR( 'P', N, N, M, A, N, DUM(1),
$ DUM(1), -1, IERR )	$ DUM(1), -1, IERR )
LWORK_DORGBR_P=DUM(1)	LWORK_DORGBR_P = INT( DUM(1) )
MAXWRK = MAX( MAXWRK, 3*M+LWORK_DORGBR_P )	MAXWRK = MAX( MAXWRK, 3*M + LWORK_DORGBR_P )
END IF	END IF
IF( .NOT.WNTUN ) THEN	IF( .NOT.WNTUN ) THEN
MAXWRK = MAX( MAXWRK, 3*M+LWORK_DORGBR_Q )	MAXWRK = MAX( MAXWRK, 3*M + LWORK_DORGBR_Q )
END IF	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 685	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 708	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 739	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 762	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 774	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 784	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 800	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 809	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 830	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 863	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 886	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 899	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 909	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 917	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 933	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 942	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 961	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 974	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 984	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 1042	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 1055	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 1065	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 1073	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 1082	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 1103	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 1121	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 1167	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 1186	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 1199	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 1210	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 1221	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 1239	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 1266	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 1284	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 1346	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 1359	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 1369	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 1379	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 1396	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 1417	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 1441	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 1449	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 1503	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 1517	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 1527	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 1535	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 1544	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 1569	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 1587	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 1599	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 1634	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 1653	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 1678	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 1689	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 1707	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 1737	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 1755	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 1767	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 1818	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 1842	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 1852	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 1869	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 1894	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 1918	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 1926	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 1967	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 1976	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 1990	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 2000	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 2009	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 2071	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 2085	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 2093	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 2126	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 2152	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 2164	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 2174	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 2190	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 2199	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 2220	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 2253	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 2279	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 2291	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 2301	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 2309	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 2325	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 2334	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 2353	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 2365	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 2375	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 2433	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 2446	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 2456	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 2465	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 2474	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 2495	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 2505	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 2520	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 2562	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 2581	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 2594	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 2605	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 2616	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 2624	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 2634	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 2661	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 2683	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 2741	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 2754	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 2764	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 2774	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 2782	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 2791	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 2812	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 2835	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 2843	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 2877	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 2897	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 2911	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 2921	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 2929	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 2939	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 2964	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 2986	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 2994	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 3017	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 3029	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 3048	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 3073	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 3084	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 3092	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 3102	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 3132	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 3154	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 3162	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 3193	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 3213	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 3237	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 3247	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 3263	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 3288	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 3311	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 3319	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 3360	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 3369	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 3379	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 3393	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 3402	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.14
changed lines
	Added in v.1.15