--- rpl/lapack/lapack/dgesvd.f 2011/07/22 07:38:05 1.8
+++ rpl/lapack/lapack/dgesvd.f 2023/08/07 08:38:50 1.20
@@ -1,10 +1,217 @@
- SUBROUTINE DGESVD( JOBU, JOBVT, M, N, A, LDA, S, U, LDU, VT, LDVT,
- $ WORK, LWORK, INFO )
+*> \brief DGESVD computes the singular value decomposition (SVD) for GE matrices
*
-* -- LAPACK driver routine (version 3.3.1) --
+* =========== DOCUMENTATION ===========
+*
+* Online html documentation available at
+* http://www.netlib.org/lapack/explore-html/
+*
+*> \htmlonly
+*> Download DGESVD + dependencies
+*>
+*> [TGZ]
+*>
+*> [ZIP]
+*>
+*> [TXT]
+*> \endhtmlonly
+*
+* Definition:
+* ===========
+*
+* SUBROUTINE DGESVD( JOBU, JOBVT, M, N, A, LDA, S, U, LDU, VT, LDVT,
+* WORK, LWORK, INFO )
+*
+* .. Scalar Arguments ..
+* CHARACTER JOBU, JOBVT
+* INTEGER INFO, LDA, LDU, LDVT, LWORK, M, N
+* ..
+* .. Array Arguments ..
+* DOUBLE PRECISION A( LDA, * ), S( * ), U( LDU, * ),
+* $ VT( LDVT, * ), WORK( * )
+* ..
+*
+*
+*> \par Purpose:
+* =============
+*>
+*> \verbatim
+*>
+*> DGESVD computes the singular value decomposition (SVD) of a real
+*> M-by-N matrix A, optionally computing the left and/or right singular
+*> vectors. The SVD is written
+*>
+*> A = U * SIGMA * transpose(V)
+*>
+*> where SIGMA is an M-by-N matrix which is zero except for its
+*> min(m,n) diagonal elements, U is an M-by-M orthogonal matrix, and
+*> V is an N-by-N orthogonal matrix. The diagonal elements of SIGMA
+*> are the singular values of A; they are real and non-negative, and
+*> are returned in descending order. The first min(m,n) columns of
+*> U and V are the left and right singular vectors of A.
+*>
+*> Note that the routine returns V**T, not V.
+*> \endverbatim
+*
+* Arguments:
+* ==========
+*
+*> \param[in] JOBU
+*> \verbatim
+*> JOBU is CHARACTER*1
+*> Specifies options for computing all or part of the matrix U:
+*> = 'A': all M columns of U are returned in array U:
+*> = 'S': the first min(m,n) columns of U (the left singular
+*> vectors) are returned in the array U;
+*> = 'O': the first min(m,n) columns of U (the left singular
+*> vectors) are overwritten on the array A;
+*> = 'N': no columns of U (no left singular vectors) are
+*> computed.
+*> \endverbatim
+*>
+*> \param[in] JOBVT
+*> \verbatim
+*> JOBVT is CHARACTER*1
+*> Specifies options for computing all or part of the matrix
+*> V**T:
+*> = 'A': all N rows of V**T are returned in the array VT;
+*> = 'S': the first min(m,n) rows of V**T (the right singular
+*> vectors) are returned in the array VT;
+*> = 'O': the first min(m,n) rows of V**T (the right singular
+*> vectors) are overwritten on the array A;
+*> = 'N': no rows of V**T (no right singular vectors) are
+*> computed.
+*>
+*> JOBVT and JOBU cannot both be 'O'.
+*> \endverbatim
+*>
+*> \param[in] M
+*> \verbatim
+*> M is INTEGER
+*> The number of rows of the input matrix A. M >= 0.
+*> \endverbatim
+*>
+*> \param[in] N
+*> \verbatim
+*> N is INTEGER
+*> The number of columns of the input matrix A. N >= 0.
+*> \endverbatim
+*>
+*> \param[in,out] A
+*> \verbatim
+*> A is DOUBLE PRECISION array, dimension (LDA,N)
+*> On entry, the M-by-N matrix A.
+*> On exit,
+*> if JOBU = 'O', A is overwritten with the first min(m,n)
+*> columns of U (the left singular vectors,
+*> stored columnwise);
+*> if JOBVT = 'O', A is overwritten with the first min(m,n)
+*> rows of V**T (the right singular vectors,
+*> stored rowwise);
+*> if JOBU .ne. 'O' and JOBVT .ne. 'O', the contents of A
+*> are destroyed.
+*> \endverbatim
+*>
+*> \param[in] LDA
+*> \verbatim
+*> LDA is INTEGER
+*> The leading dimension of the array A. LDA >= max(1,M).
+*> \endverbatim
+*>
+*> \param[out] S
+*> \verbatim
+*> S is DOUBLE PRECISION array, dimension (min(M,N))
+*> The singular values of A, sorted so that S(i) >= S(i+1).
+*> \endverbatim
+*>
+*> \param[out] U
+*> \verbatim
+*> U is DOUBLE PRECISION array, dimension (LDU,UCOL)
+*> (LDU,M) if JOBU = 'A' or (LDU,min(M,N)) if JOBU = 'S'.
+*> If JOBU = 'A', U contains the M-by-M orthogonal matrix U;
+*> if JOBU = 'S', U contains the first min(m,n) columns of U
+*> (the left singular vectors, stored columnwise);
+*> if JOBU = 'N' or 'O', U is not referenced.
+*> \endverbatim
+*>
+*> \param[in] LDU
+*> \verbatim
+*> LDU is INTEGER
+*> The leading dimension of the array U. LDU >= 1; if
+*> JOBU = 'S' or 'A', LDU >= M.
+*> \endverbatim
+*>
+*> \param[out] VT
+*> \verbatim
+*> VT is DOUBLE PRECISION array, dimension (LDVT,N)
+*> If JOBVT = 'A', VT contains the N-by-N orthogonal matrix
+*> V**T;
+*> if JOBVT = 'S', VT contains the first min(m,n) rows of
+*> V**T (the right singular vectors, stored rowwise);
+*> if JOBVT = 'N' or 'O', VT is not referenced.
+*> \endverbatim
+*>
+*> \param[in] LDVT
+*> \verbatim
+*> LDVT is INTEGER
+*> The leading dimension of the array VT. LDVT >= 1; if
+*> JOBVT = 'A', LDVT >= N; if JOBVT = 'S', LDVT >= min(M,N).
+*> \endverbatim
+*>
+*> \param[out] WORK
+*> \verbatim
+*> WORK is DOUBLE PRECISION array, dimension (MAX(1,LWORK))
+*> On exit, if INFO = 0, WORK(1) returns the optimal LWORK;
+*> if INFO > 0, WORK(2:MIN(M,N)) contains the unconverged
+*> superdiagonal elements of an upper bidiagonal matrix B
+*> whose diagonal is in S (not necessarily sorted). B
+*> satisfies A = U * B * VT, so it has the same singular values
+*> as A, and singular vectors related by U and VT.
+*> \endverbatim
+*>
+*> \param[in] LWORK
+*> \verbatim
+*> LWORK is INTEGER
+*> The dimension of the array WORK.
+*> LWORK >= MAX(1,5*MIN(M,N)) for the paths (see comments inside code):
+*> - PATH 1 (M much larger than N, JOBU='N')
+*> - PATH 1t (N much larger than M, JOBVT='N')
+*> LWORK >= MAX(1,3*MIN(M,N) + MAX(M,N),5*MIN(M,N)) for the other paths
+*> For good performance, LWORK should generally be larger.
+*>
+*> If LWORK = -1, then a workspace query is assumed; the routine
+*> only calculates the optimal size of the WORK array, returns
+*> this value as the first entry of the WORK array, and no error
+*> message related to LWORK is issued by XERBLA.
+*> \endverbatim
+*>
+*> \param[out] INFO
+*> \verbatim
+*> INFO is INTEGER
+*> = 0: successful exit.
+*> < 0: if INFO = -i, the i-th argument had an illegal value.
+*> > 0: if DBDSQR did not converge, INFO specifies how many
+*> superdiagonals of an intermediate bidiagonal form B
+*> did not converge to zero. See the description of WORK
+*> above for details.
+*> \endverbatim
+*
+* Authors:
+* ========
+*
+*> \author Univ. of Tennessee
+*> \author Univ. of California Berkeley
+*> \author Univ. of Colorado Denver
+*> \author NAG Ltd.
+*
+*> \ingroup doubleGEsing
+*
+* =====================================================================
+ SUBROUTINE DGESVD( JOBU, JOBVT, M, N, A, LDA, S, U, LDU,
+ $ VT, LDVT, WORK, LWORK, INFO )
+*
+* -- LAPACK driver routine --
* -- LAPACK is a software package provided by Univ. of Tennessee, --
* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
-* -- April 2011 --
*
* .. Scalar Arguments ..
CHARACTER JOBU, JOBVT
@@ -15,125 +222,6 @@
$ VT( LDVT, * ), WORK( * )
* ..
*
-* Purpose
-* =======
-*
-* DGESVD computes the singular value decomposition (SVD) of a real
-* M-by-N matrix A, optionally computing the left and/or right singular
-* vectors. The SVD is written
-*
-* A = U * SIGMA * transpose(V)
-*
-* where SIGMA is an M-by-N matrix which is zero except for its
-* min(m,n) diagonal elements, U is an M-by-M orthogonal matrix, and
-* V is an N-by-N orthogonal matrix. The diagonal elements of SIGMA
-* are the singular values of A; they are real and non-negative, and
-* are returned in descending order. The first min(m,n) columns of
-* U and V are the left and right singular vectors of A.
-*
-* Note that the routine returns V**T, not V.
-*
-* Arguments
-* =========
-*
-* JOBU (input) CHARACTER*1
-* Specifies options for computing all or part of the matrix U:
-* = 'A': all M columns of U are returned in array U:
-* = 'S': the first min(m,n) columns of U (the left singular
-* vectors) are returned in the array U;
-* = 'O': the first min(m,n) columns of U (the left singular
-* vectors) are overwritten on the array A;
-* = 'N': no columns of U (no left singular vectors) are
-* computed.
-*
-* JOBVT (input) CHARACTER*1
-* Specifies options for computing all or part of the matrix
-* V**T:
-* = 'A': all N rows of V**T are returned in the array VT;
-* = 'S': the first min(m,n) rows of V**T (the right singular
-* vectors) are returned in the array VT;
-* = 'O': the first min(m,n) rows of V**T (the right singular
-* vectors) are overwritten on the array A;
-* = 'N': no rows of V**T (no right singular vectors) are
-* computed.
-*
-* JOBVT and JOBU cannot both be 'O'.
-*
-* M (input) INTEGER
-* The number of rows of the input matrix A. M >= 0.
-*
-* N (input) INTEGER
-* The number of columns of the input matrix A. N >= 0.
-*
-* A (input/output) DOUBLE PRECISION array, dimension (LDA,N)
-* On entry, the M-by-N matrix A.
-* On exit,
-* if JOBU = 'O', A is overwritten with the first min(m,n)
-* columns of U (the left singular vectors,
-* stored columnwise);
-* if JOBVT = 'O', A is overwritten with the first min(m,n)
-* rows of V**T (the right singular vectors,
-* stored rowwise);
-* if JOBU .ne. 'O' and JOBVT .ne. 'O', the contents of A
-* are destroyed.
-*
-* LDA (input) INTEGER
-* The leading dimension of the array A. LDA >= max(1,M).
-*
-* S (output) DOUBLE PRECISION array, dimension (min(M,N))
-* The singular values of A, sorted so that S(i) >= S(i+1).
-*
-* U (output) DOUBLE PRECISION array, dimension (LDU,UCOL)
-* (LDU,M) if JOBU = 'A' or (LDU,min(M,N)) if JOBU = 'S'.
-* If JOBU = 'A', U contains the M-by-M orthogonal matrix U;
-* if JOBU = 'S', U contains the first min(m,n) columns of U
-* (the left singular vectors, stored columnwise);
-* if JOBU = 'N' or 'O', U is not referenced.
-*
-* LDU (input) INTEGER
-* The leading dimension of the array U. LDU >= 1; if
-* JOBU = 'S' or 'A', LDU >= M.
-*
-* VT (output) DOUBLE PRECISION array, dimension (LDVT,N)
-* If JOBVT = 'A', VT contains the N-by-N orthogonal matrix
-* V**T;
-* if JOBVT = 'S', VT contains the first min(m,n) rows of
-* V**T (the right singular vectors, stored rowwise);
-* if JOBVT = 'N' or 'O', VT is not referenced.
-*
-* LDVT (input) INTEGER
-* The leading dimension of the array VT. LDVT >= 1; if
-* JOBVT = 'A', LDVT >= N; if JOBVT = 'S', LDVT >= min(M,N).
-*
-* WORK (workspace/output) DOUBLE PRECISION array, dimension (MAX(1,LWORK))
-* On exit, if INFO = 0, WORK(1) returns the optimal LWORK;
-* if INFO > 0, WORK(2:MIN(M,N)) contains the unconverged
-* superdiagonal elements of an upper bidiagonal matrix B
-* whose diagonal is in S (not necessarily sorted). B
-* satisfies A = U * B * VT, so it has the same singular values
-* as A, and singular vectors related by U and VT.
-*
-* LWORK (input) INTEGER
-* The dimension of the array WORK.
-* LWORK >= MAX(1,5*MIN(M,N)) for the paths (see comments inside code):
-* - PATH 1 (M much larger than N, JOBU='N')
-* - PATH 1t (N much larger than M, JOBVT='N')
-* LWORK >= MAX(1,3*MIN(M,N)+MAX(M,N),5*MIN(M,N)) for the other paths
-* For good performance, LWORK should generally be larger.
-*
-* If LWORK = -1, then a workspace query is assumed; the routine
-* only calculates the optimal size of the WORK array, returns
-* this value as the first entry of the WORK array, and no error
-* message related to LWORK is issued by XERBLA.
-*
-* INFO (output) INTEGER
-* = 0: successful exit.
-* < 0: if INFO = -i, the i-th argument had an illegal value.
-* > 0: if DBDSQR did not converge, INFO specifies how many
-* superdiagonals of an intermediate bidiagonal form B
-* did not converge to zero. See the description of WORK
-* above for details.
-*
* =====================================================================
*
* .. Parameters ..
@@ -147,6 +235,9 @@
$ ITAU, ITAUP, ITAUQ, IU, IWORK, LDWRKR, LDWRKU,
$ MAXWRK, MINMN, MINWRK, MNTHR, NCU, NCVT, NRU,
$ NRVT, WRKBL
+ INTEGER LWORK_DGEQRF, LWORK_DORGQR_N, LWORK_DORGQR_M,
+ $ LWORK_DGEBRD, LWORK_DORGBR_P, LWORK_DORGBR_Q,
+ $ LWORK_DGELQF, LWORK_DORGLQ_N, LWORK_DORGLQ_M
DOUBLE PRECISION ANRM, BIGNUM, EPS, SMLNUM
* ..
* .. Local Arrays ..
@@ -218,163 +309,160 @@
*
MNTHR = ILAENV( 6, 'DGESVD', JOBU // JOBVT, M, N, 0, 0 )
BDSPAC = 5*N
+* Compute space needed for DGEQRF
+ CALL DGEQRF( M, N, A, LDA, DUM(1), DUM(1), -1, IERR )
+ LWORK_DGEQRF = INT( DUM(1) )
+* Compute space needed for DORGQR
+ CALL DORGQR( M, N, N, A, LDA, DUM(1), DUM(1), -1, IERR )
+ LWORK_DORGQR_N = INT( DUM(1) )
+ CALL DORGQR( M, M, N, A, LDA, DUM(1), DUM(1), -1, IERR )
+ LWORK_DORGQR_M = INT( DUM(1) )
+* Compute space needed for DGEBRD
+ CALL DGEBRD( N, N, A, LDA, S, DUM(1), DUM(1),
+ $ DUM(1), DUM(1), -1, IERR )
+ LWORK_DGEBRD = INT( DUM(1) )
+* Compute space needed for DORGBR P
+ CALL DORGBR( 'P', N, N, N, A, LDA, DUM(1),
+ $ DUM(1), -1, IERR )
+ LWORK_DORGBR_P = INT( DUM(1) )
+* Compute space needed for DORGBR Q
+ CALL DORGBR( 'Q', N, N, N, A, LDA, DUM(1),
+ $ DUM(1), -1, IERR )
+ LWORK_DORGBR_Q = INT( DUM(1) )
+*
IF( M.GE.MNTHR ) THEN
IF( WNTUN ) THEN
*
* Path 1 (M much larger than N, JOBU='N')
*
- MAXWRK = N + N*ILAENV( 1, 'DGEQRF', ' ', M, N, -1,
- $ -1 )
- MAXWRK = MAX( MAXWRK, 3*N+2*N*
- $ ILAENV( 1, 'DGEBRD', ' ', N, N, -1, -1 ) )
+ MAXWRK = N + LWORK_DGEQRF
+ MAXWRK = MAX( MAXWRK, 3*N + LWORK_DGEBRD )
IF( WNTVO .OR. WNTVAS )
- $ MAXWRK = MAX( MAXWRK, 3*N+( N-1 )*
- $ ILAENV( 1, 'DORGBR', 'P', N, N, N, -1 ) )
+ $ MAXWRK = MAX( MAXWRK, 3*N + LWORK_DORGBR_P )
MAXWRK = MAX( MAXWRK, BDSPAC )
MINWRK = MAX( 4*N, BDSPAC )
ELSE IF( WNTUO .AND. WNTVN ) THEN
*
* Path 2 (M much larger than N, JOBU='O', JOBVT='N')
*
- WRKBL = N + N*ILAENV( 1, 'DGEQRF', ' ', M, N, -1, -1 )
- WRKBL = MAX( WRKBL, N+N*ILAENV( 1, 'DORGQR', ' ', M,
- $ N, N, -1 ) )
- WRKBL = MAX( WRKBL, 3*N+2*N*
- $ ILAENV( 1, 'DGEBRD', ' ', N, N, -1, -1 ) )
- WRKBL = MAX( WRKBL, 3*N+N*
- $ ILAENV( 1, 'DORGBR', 'Q', N, N, N, -1 ) )
+ WRKBL = N + LWORK_DGEQRF
+ WRKBL = MAX( WRKBL, N + LWORK_DORGQR_N )
+ WRKBL = MAX( WRKBL, 3*N + LWORK_DGEBRD )
+ WRKBL = MAX( WRKBL, 3*N + LWORK_DORGBR_Q )
WRKBL = MAX( WRKBL, BDSPAC )
- MAXWRK = MAX( N*N+WRKBL, N*N+M*N+N )
- MINWRK = MAX( 3*N+M, BDSPAC )
+ MAXWRK = MAX( N*N + WRKBL, N*N + M*N + N )
+ MINWRK = MAX( 3*N + M, BDSPAC )
ELSE IF( WNTUO .AND. WNTVAS ) THEN
*
* Path 3 (M much larger than N, JOBU='O', JOBVT='S' or
* 'A')
*
- WRKBL = N + N*ILAENV( 1, 'DGEQRF', ' ', M, N, -1, -1 )
- WRKBL = MAX( WRKBL, N+N*ILAENV( 1, 'DORGQR', ' ', M,
- $ N, N, -1 ) )
- WRKBL = MAX( WRKBL, 3*N+2*N*
- $ ILAENV( 1, 'DGEBRD', ' ', N, N, -1, -1 ) )
- WRKBL = MAX( WRKBL, 3*N+N*
- $ ILAENV( 1, 'DORGBR', 'Q', N, N, N, -1 ) )
- WRKBL = MAX( WRKBL, 3*N+( N-1 )*
- $ ILAENV( 1, 'DORGBR', 'P', N, N, N, -1 ) )
+ WRKBL = N + LWORK_DGEQRF
+ WRKBL = MAX( WRKBL, N + LWORK_DORGQR_N )
+ WRKBL = MAX( WRKBL, 3*N + LWORK_DGEBRD )
+ WRKBL = MAX( WRKBL, 3*N + LWORK_DORGBR_Q )
+ WRKBL = MAX( WRKBL, 3*N + LWORK_DORGBR_P )
WRKBL = MAX( WRKBL, BDSPAC )
- MAXWRK = MAX( N*N+WRKBL, N*N+M*N+N )
- MINWRK = MAX( 3*N+M, BDSPAC )
+ MAXWRK = MAX( N*N + WRKBL, N*N + M*N + N )
+ MINWRK = MAX( 3*N + M, BDSPAC )
ELSE IF( WNTUS .AND. WNTVN ) THEN
*
* Path 4 (M much larger than N, JOBU='S', JOBVT='N')
*
- WRKBL = N + N*ILAENV( 1, 'DGEQRF', ' ', M, N, -1, -1 )
- WRKBL = MAX( WRKBL, N+N*ILAENV( 1, 'DORGQR', ' ', M,
- $ N, N, -1 ) )
- WRKBL = MAX( WRKBL, 3*N+2*N*
- $ ILAENV( 1, 'DGEBRD', ' ', N, N, -1, -1 ) )
- WRKBL = MAX( WRKBL, 3*N+N*
- $ ILAENV( 1, 'DORGBR', 'Q', N, N, N, -1 ) )
+ WRKBL = N + LWORK_DGEQRF
+ WRKBL = MAX( WRKBL, N + LWORK_DORGQR_N )
+ WRKBL = MAX( WRKBL, 3*N + LWORK_DGEBRD )
+ WRKBL = MAX( WRKBL, 3*N + LWORK_DORGBR_Q )
WRKBL = MAX( WRKBL, BDSPAC )
MAXWRK = N*N + WRKBL
- MINWRK = MAX( 3*N+M, BDSPAC )
+ MINWRK = MAX( 3*N + M, BDSPAC )
ELSE IF( WNTUS .AND. WNTVO ) THEN
*
* Path 5 (M much larger than N, JOBU='S', JOBVT='O')
*
- WRKBL = N + N*ILAENV( 1, 'DGEQRF', ' ', M, N, -1, -1 )
- WRKBL = MAX( WRKBL, N+N*ILAENV( 1, 'DORGQR', ' ', M,
- $ N, N, -1 ) )
- WRKBL = MAX( WRKBL, 3*N+2*N*
- $ ILAENV( 1, 'DGEBRD', ' ', N, N, -1, -1 ) )
- WRKBL = MAX( WRKBL, 3*N+N*
- $ ILAENV( 1, 'DORGBR', 'Q', N, N, N, -1 ) )
- WRKBL = MAX( WRKBL, 3*N+( N-1 )*
- $ ILAENV( 1, 'DORGBR', 'P', N, N, N, -1 ) )
+ WRKBL = N + LWORK_DGEQRF
+ WRKBL = MAX( WRKBL, N + LWORK_DORGQR_N )
+ WRKBL = MAX( WRKBL, 3*N + LWORK_DGEBRD )
+ WRKBL = MAX( WRKBL, 3*N + LWORK_DORGBR_Q )
+ WRKBL = MAX( WRKBL, 3*N + LWORK_DORGBR_P )
WRKBL = MAX( WRKBL, BDSPAC )
MAXWRK = 2*N*N + WRKBL
- MINWRK = MAX( 3*N+M, BDSPAC )
+ MINWRK = MAX( 3*N + M, BDSPAC )
ELSE IF( WNTUS .AND. WNTVAS ) THEN
*
* Path 6 (M much larger than N, JOBU='S', JOBVT='S' or
* 'A')
*
- WRKBL = N + N*ILAENV( 1, 'DGEQRF', ' ', M, N, -1, -1 )
- WRKBL = MAX( WRKBL, N+N*ILAENV( 1, 'DORGQR', ' ', M,
- $ N, N, -1 ) )
- WRKBL = MAX( WRKBL, 3*N+2*N*
- $ ILAENV( 1, 'DGEBRD', ' ', N, N, -1, -1 ) )
- WRKBL = MAX( WRKBL, 3*N+N*
- $ ILAENV( 1, 'DORGBR', 'Q', N, N, N, -1 ) )
- WRKBL = MAX( WRKBL, 3*N+( N-1 )*
- $ ILAENV( 1, 'DORGBR', 'P', N, N, N, -1 ) )
+ WRKBL = N + LWORK_DGEQRF
+ WRKBL = MAX( WRKBL, N + LWORK_DORGQR_N )
+ WRKBL = MAX( WRKBL, 3*N + LWORK_DGEBRD )
+ WRKBL = MAX( WRKBL, 3*N + LWORK_DORGBR_Q )
+ WRKBL = MAX( WRKBL, 3*N + LWORK_DORGBR_P )
WRKBL = MAX( WRKBL, BDSPAC )
MAXWRK = N*N + WRKBL
- MINWRK = MAX( 3*N+M, BDSPAC )
+ MINWRK = MAX( 3*N + M, BDSPAC )
ELSE IF( WNTUA .AND. WNTVN ) THEN
*
* Path 7 (M much larger than N, JOBU='A', JOBVT='N')
*
- WRKBL = N + N*ILAENV( 1, 'DGEQRF', ' ', M, N, -1, -1 )
- WRKBL = MAX( WRKBL, N+M*ILAENV( 1, 'DORGQR', ' ', M,
- $ M, N, -1 ) )
- WRKBL = MAX( WRKBL, 3*N+2*N*
- $ ILAENV( 1, 'DGEBRD', ' ', N, N, -1, -1 ) )
- WRKBL = MAX( WRKBL, 3*N+N*
- $ ILAENV( 1, 'DORGBR', 'Q', N, N, N, -1 ) )
+ WRKBL = N + LWORK_DGEQRF
+ WRKBL = MAX( WRKBL, N + LWORK_DORGQR_M )
+ WRKBL = MAX( WRKBL, 3*N + LWORK_DGEBRD )
+ WRKBL = MAX( WRKBL, 3*N + LWORK_DORGBR_Q )
WRKBL = MAX( WRKBL, BDSPAC )
MAXWRK = N*N + WRKBL
- MINWRK = MAX( 3*N+M, BDSPAC )
+ MINWRK = MAX( 3*N + M, BDSPAC )
ELSE IF( WNTUA .AND. WNTVO ) THEN
*
* Path 8 (M much larger than N, JOBU='A', JOBVT='O')
*
- WRKBL = N + N*ILAENV( 1, 'DGEQRF', ' ', M, N, -1, -1 )
- WRKBL = MAX( WRKBL, N+M*ILAENV( 1, 'DORGQR', ' ', M,
- $ M, N, -1 ) )
- WRKBL = MAX( WRKBL, 3*N+2*N*
- $ ILAENV( 1, 'DGEBRD', ' ', N, N, -1, -1 ) )
- WRKBL = MAX( WRKBL, 3*N+N*
- $ ILAENV( 1, 'DORGBR', 'Q', N, N, N, -1 ) )
- WRKBL = MAX( WRKBL, 3*N+( N-1 )*
- $ ILAENV( 1, 'DORGBR', 'P', N, N, N, -1 ) )
+ WRKBL = N + LWORK_DGEQRF
+ WRKBL = MAX( WRKBL, N + LWORK_DORGQR_M )
+ WRKBL = MAX( WRKBL, 3*N + LWORK_DGEBRD )
+ WRKBL = MAX( WRKBL, 3*N + LWORK_DORGBR_Q )
+ WRKBL = MAX( WRKBL, 3*N + LWORK_DORGBR_P )
WRKBL = MAX( WRKBL, BDSPAC )
MAXWRK = 2*N*N + WRKBL
- MINWRK = MAX( 3*N+M, BDSPAC )
+ MINWRK = MAX( 3*N + M, BDSPAC )
ELSE IF( WNTUA .AND. WNTVAS ) THEN
*
* Path 9 (M much larger than N, JOBU='A', JOBVT='S' or
* 'A')
*
- WRKBL = N + N*ILAENV( 1, 'DGEQRF', ' ', M, N, -1, -1 )
- WRKBL = MAX( WRKBL, N+M*ILAENV( 1, 'DORGQR', ' ', M,
- $ M, N, -1 ) )
- WRKBL = MAX( WRKBL, 3*N+2*N*
- $ ILAENV( 1, 'DGEBRD', ' ', N, N, -1, -1 ) )
- WRKBL = MAX( WRKBL, 3*N+N*
- $ ILAENV( 1, 'DORGBR', 'Q', N, N, N, -1 ) )
- WRKBL = MAX( WRKBL, 3*N+( N-1 )*
- $ ILAENV( 1, 'DORGBR', 'P', N, N, N, -1 ) )
+ WRKBL = N + LWORK_DGEQRF
+ WRKBL = MAX( WRKBL, N + LWORK_DORGQR_M )
+ WRKBL = MAX( WRKBL, 3*N + LWORK_DGEBRD )
+ WRKBL = MAX( WRKBL, 3*N + LWORK_DORGBR_Q )
+ WRKBL = MAX( WRKBL, 3*N + LWORK_DORGBR_P )
WRKBL = MAX( WRKBL, BDSPAC )
MAXWRK = N*N + WRKBL
- MINWRK = MAX( 3*N+M, BDSPAC )
+ MINWRK = MAX( 3*N + M, BDSPAC )
END IF
ELSE
*
* Path 10 (M at least N, but not much larger)
*
- MAXWRK = 3*N + ( M+N )*ILAENV( 1, 'DGEBRD', ' ', M, N,
- $ -1, -1 )
- IF( WNTUS .OR. WNTUO )
- $ MAXWRK = MAX( MAXWRK, 3*N+N*
- $ ILAENV( 1, 'DORGBR', 'Q', M, N, N, -1 ) )
- IF( WNTUA )
- $ MAXWRK = MAX( MAXWRK, 3*N+M*
- $ ILAENV( 1, 'DORGBR', 'Q', M, M, N, -1 ) )
- IF( .NOT.WNTVN )
- $ MAXWRK = MAX( MAXWRK, 3*N+( N-1 )*
- $ ILAENV( 1, 'DORGBR', 'P', N, N, N, -1 ) )
+ CALL DGEBRD( M, N, A, LDA, S, DUM(1), DUM(1),
+ $ DUM(1), DUM(1), -1, IERR )
+ LWORK_DGEBRD = INT( DUM(1) )
+ MAXWRK = 3*N + LWORK_DGEBRD
+ IF( WNTUS .OR. WNTUO ) THEN
+ CALL DORGBR( 'Q', M, N, N, A, LDA, DUM(1),
+ $ DUM(1), -1, IERR )
+ LWORK_DORGBR_Q = INT( DUM(1) )
+ MAXWRK = MAX( MAXWRK, 3*N + LWORK_DORGBR_Q )
+ END IF
+ IF( WNTUA ) THEN
+ CALL DORGBR( 'Q', M, M, N, A, LDA, DUM(1),
+ $ DUM(1), -1, IERR )
+ LWORK_DORGBR_Q = INT( DUM(1) )
+ MAXWRK = MAX( MAXWRK, 3*N + LWORK_DORGBR_Q )
+ END IF
+ IF( .NOT.WNTVN ) THEN
+ MAXWRK = MAX( MAXWRK, 3*N + LWORK_DORGBR_P )
+ END IF
MAXWRK = MAX( MAXWRK, BDSPAC )
- MINWRK = MAX( 3*N+M, BDSPAC )
+ MINWRK = MAX( 3*N + M, BDSPAC )
END IF
ELSE IF( MINMN.GT.0 ) THEN
*
@@ -382,163 +470,160 @@
*
MNTHR = ILAENV( 6, 'DGESVD', JOBU // JOBVT, M, N, 0, 0 )
BDSPAC = 5*M
+* Compute space needed for DGELQF
+ CALL DGELQF( M, N, A, LDA, DUM(1), DUM(1), -1, IERR )
+ LWORK_DGELQF = INT( DUM(1) )
+* Compute space needed for DORGLQ
+ CALL DORGLQ( N, N, M, DUM(1), N, DUM(1), DUM(1), -1, IERR )
+ LWORK_DORGLQ_N = INT( DUM(1) )
+ CALL DORGLQ( M, N, M, A, LDA, DUM(1), DUM(1), -1, IERR )
+ LWORK_DORGLQ_M = INT( DUM(1) )
+* Compute space needed for DGEBRD
+ CALL DGEBRD( M, M, A, LDA, S, DUM(1), DUM(1),
+ $ DUM(1), DUM(1), -1, IERR )
+ LWORK_DGEBRD = INT( DUM(1) )
+* Compute space needed for DORGBR P
+ CALL DORGBR( 'P', M, M, M, A, N, DUM(1),
+ $ DUM(1), -1, IERR )
+ LWORK_DORGBR_P = INT( DUM(1) )
+* Compute space needed for DORGBR Q
+ CALL DORGBR( 'Q', M, M, M, A, N, DUM(1),
+ $ DUM(1), -1, IERR )
+ LWORK_DORGBR_Q = INT( DUM(1) )
IF( N.GE.MNTHR ) THEN
IF( WNTVN ) THEN
*
* Path 1t(N much larger than M, JOBVT='N')
*
- MAXWRK = M + M*ILAENV( 1, 'DGELQF', ' ', M, N, -1,
- $ -1 )
- MAXWRK = MAX( MAXWRK, 3*M+2*M*
- $ ILAENV( 1, 'DGEBRD', ' ', M, M, -1, -1 ) )
+ MAXWRK = M + LWORK_DGELQF
+ MAXWRK = MAX( MAXWRK, 3*M + LWORK_DGEBRD )
IF( WNTUO .OR. WNTUAS )
- $ MAXWRK = MAX( MAXWRK, 3*M+M*
- $ ILAENV( 1, 'DORGBR', 'Q', M, M, M, -1 ) )
+ $ MAXWRK = MAX( MAXWRK, 3*M + LWORK_DORGBR_Q )
MAXWRK = MAX( MAXWRK, BDSPAC )
MINWRK = MAX( 4*M, BDSPAC )
ELSE IF( WNTVO .AND. WNTUN ) THEN
*
* Path 2t(N much larger than M, JOBU='N', JOBVT='O')
*
- WRKBL = M + M*ILAENV( 1, 'DGELQF', ' ', M, N, -1, -1 )
- WRKBL = MAX( WRKBL, M+M*ILAENV( 1, 'DORGLQ', ' ', M,
- $ N, M, -1 ) )
- WRKBL = MAX( WRKBL, 3*M+2*M*
- $ ILAENV( 1, 'DGEBRD', ' ', M, M, -1, -1 ) )
- WRKBL = MAX( WRKBL, 3*M+( M-1 )*
- $ ILAENV( 1, 'DORGBR', 'P', M, M, M, -1 ) )
+ WRKBL = M + LWORK_DGELQF
+ WRKBL = MAX( WRKBL, M + LWORK_DORGLQ_M )
+ WRKBL = MAX( WRKBL, 3*M + LWORK_DGEBRD )
+ WRKBL = MAX( WRKBL, 3*M + LWORK_DORGBR_P )
WRKBL = MAX( WRKBL, BDSPAC )
- MAXWRK = MAX( M*M+WRKBL, M*M+M*N+M )
- MINWRK = MAX( 3*M+N, BDSPAC )
+ MAXWRK = MAX( M*M + WRKBL, M*M + M*N + M )
+ MINWRK = MAX( 3*M + N, BDSPAC )
ELSE IF( WNTVO .AND. WNTUAS ) THEN
*
* Path 3t(N much larger than M, JOBU='S' or 'A',
* JOBVT='O')
*
- WRKBL = M + M*ILAENV( 1, 'DGELQF', ' ', M, N, -1, -1 )
- WRKBL = MAX( WRKBL, M+M*ILAENV( 1, 'DORGLQ', ' ', M,
- $ N, M, -1 ) )
- WRKBL = MAX( WRKBL, 3*M+2*M*
- $ ILAENV( 1, 'DGEBRD', ' ', M, M, -1, -1 ) )
- WRKBL = MAX( WRKBL, 3*M+( M-1 )*
- $ ILAENV( 1, 'DORGBR', 'P', M, M, M, -1 ) )
- WRKBL = MAX( WRKBL, 3*M+M*
- $ ILAENV( 1, 'DORGBR', 'Q', M, M, M, -1 ) )
+ WRKBL = M + LWORK_DGELQF
+ WRKBL = MAX( WRKBL, M + LWORK_DORGLQ_M )
+ WRKBL = MAX( WRKBL, 3*M + LWORK_DGEBRD )
+ WRKBL = MAX( WRKBL, 3*M + LWORK_DORGBR_P )
+ WRKBL = MAX( WRKBL, 3*M + LWORK_DORGBR_Q )
WRKBL = MAX( WRKBL, BDSPAC )
- MAXWRK = MAX( M*M+WRKBL, M*M+M*N+M )
- MINWRK = MAX( 3*M+N, BDSPAC )
+ MAXWRK = MAX( M*M + WRKBL, M*M + M*N + M )
+ MINWRK = MAX( 3*M + N, BDSPAC )
ELSE IF( WNTVS .AND. WNTUN ) THEN
*
* Path 4t(N much larger than M, JOBU='N', JOBVT='S')
*
- WRKBL = M + M*ILAENV( 1, 'DGELQF', ' ', M, N, -1, -1 )
- WRKBL = MAX( WRKBL, M+M*ILAENV( 1, 'DORGLQ', ' ', M,
- $ N, M, -1 ) )
- WRKBL = MAX( WRKBL, 3*M+2*M*
- $ ILAENV( 1, 'DGEBRD', ' ', M, M, -1, -1 ) )
- WRKBL = MAX( WRKBL, 3*M+( M-1 )*
- $ ILAENV( 1, 'DORGBR', 'P', M, M, M, -1 ) )
+ WRKBL = M + LWORK_DGELQF
+ WRKBL = MAX( WRKBL, M + LWORK_DORGLQ_M )
+ WRKBL = MAX( WRKBL, 3*M + LWORK_DGEBRD )
+ WRKBL = MAX( WRKBL, 3*M + LWORK_DORGBR_P )
WRKBL = MAX( WRKBL, BDSPAC )
MAXWRK = M*M + WRKBL
- MINWRK = MAX( 3*M+N, BDSPAC )
+ MINWRK = MAX( 3*M + N, BDSPAC )
ELSE IF( WNTVS .AND. WNTUO ) THEN
*
* Path 5t(N much larger than M, JOBU='O', JOBVT='S')
*
- WRKBL = M + M*ILAENV( 1, 'DGELQF', ' ', M, N, -1, -1 )
- WRKBL = MAX( WRKBL, M+M*ILAENV( 1, 'DORGLQ', ' ', M,
- $ N, M, -1 ) )
- WRKBL = MAX( WRKBL, 3*M+2*M*
- $ ILAENV( 1, 'DGEBRD', ' ', M, M, -1, -1 ) )
- WRKBL = MAX( WRKBL, 3*M+( M-1 )*
- $ ILAENV( 1, 'DORGBR', 'P', M, M, M, -1 ) )
- WRKBL = MAX( WRKBL, 3*M+M*
- $ ILAENV( 1, 'DORGBR', 'Q', M, M, M, -1 ) )
+ WRKBL = M + LWORK_DGELQF
+ WRKBL = MAX( WRKBL, M + LWORK_DORGLQ_M )
+ WRKBL = MAX( WRKBL, 3*M + LWORK_DGEBRD )
+ WRKBL = MAX( WRKBL, 3*M + LWORK_DORGBR_P )
+ WRKBL = MAX( WRKBL, 3*M + LWORK_DORGBR_Q )
WRKBL = MAX( WRKBL, BDSPAC )
MAXWRK = 2*M*M + WRKBL
- MINWRK = MAX( 3*M+N, BDSPAC )
+ MINWRK = MAX( 3*M + N, BDSPAC )
ELSE IF( WNTVS .AND. WNTUAS ) THEN
*
* Path 6t(N much larger than M, JOBU='S' or 'A',
* JOBVT='S')
*
- WRKBL = M + M*ILAENV( 1, 'DGELQF', ' ', M, N, -1, -1 )
- WRKBL = MAX( WRKBL, M+M*ILAENV( 1, 'DORGLQ', ' ', M,
- $ N, M, -1 ) )
- WRKBL = MAX( WRKBL, 3*M+2*M*
- $ ILAENV( 1, 'DGEBRD', ' ', M, M, -1, -1 ) )
- WRKBL = MAX( WRKBL, 3*M+( M-1 )*
- $ ILAENV( 1, 'DORGBR', 'P', M, M, M, -1 ) )
- WRKBL = MAX( WRKBL, 3*M+M*
- $ ILAENV( 1, 'DORGBR', 'Q', M, M, M, -1 ) )
+ WRKBL = M + LWORK_DGELQF
+ WRKBL = MAX( WRKBL, M + LWORK_DORGLQ_M )
+ WRKBL = MAX( WRKBL, 3*M + LWORK_DGEBRD )
+ WRKBL = MAX( WRKBL, 3*M + LWORK_DORGBR_P )
+ WRKBL = MAX( WRKBL, 3*M + LWORK_DORGBR_Q )
WRKBL = MAX( WRKBL, BDSPAC )
MAXWRK = M*M + WRKBL
- MINWRK = MAX( 3*M+N, BDSPAC )
+ MINWRK = MAX( 3*M + N, BDSPAC )
ELSE IF( WNTVA .AND. WNTUN ) THEN
*
* Path 7t(N much larger than M, JOBU='N', JOBVT='A')
*
- WRKBL = M + M*ILAENV( 1, 'DGELQF', ' ', M, N, -1, -1 )
- WRKBL = MAX( WRKBL, M+N*ILAENV( 1, 'DORGLQ', ' ', N,
- $ N, M, -1 ) )
- WRKBL = MAX( WRKBL, 3*M+2*M*
- $ ILAENV( 1, 'DGEBRD', ' ', M, M, -1, -1 ) )
- WRKBL = MAX( WRKBL, 3*M+( M-1 )*
- $ ILAENV( 1, 'DORGBR', 'P', M, M, M, -1 ) )
+ WRKBL = M + LWORK_DGELQF
+ WRKBL = MAX( WRKBL, M + LWORK_DORGLQ_N )
+ WRKBL = MAX( WRKBL, 3*M + LWORK_DGEBRD )
+ WRKBL = MAX( WRKBL, 3*M + LWORK_DORGBR_P )
WRKBL = MAX( WRKBL, BDSPAC )
MAXWRK = M*M + WRKBL
- MINWRK = MAX( 3*M+N, BDSPAC )
+ MINWRK = MAX( 3*M + N, BDSPAC )
ELSE IF( WNTVA .AND. WNTUO ) THEN
*
* Path 8t(N much larger than M, JOBU='O', JOBVT='A')
*
- WRKBL = M + M*ILAENV( 1, 'DGELQF', ' ', M, N, -1, -1 )
- WRKBL = MAX( WRKBL, M+N*ILAENV( 1, 'DORGLQ', ' ', N,
- $ N, M, -1 ) )
- WRKBL = MAX( WRKBL, 3*M+2*M*
- $ ILAENV( 1, 'DGEBRD', ' ', M, M, -1, -1 ) )
- WRKBL = MAX( WRKBL, 3*M+( M-1 )*
- $ ILAENV( 1, 'DORGBR', 'P', M, M, M, -1 ) )
- WRKBL = MAX( WRKBL, 3*M+M*
- $ ILAENV( 1, 'DORGBR', 'Q', M, M, M, -1 ) )
+ WRKBL = M + LWORK_DGELQF
+ WRKBL = MAX( WRKBL, M + LWORK_DORGLQ_N )
+ WRKBL = MAX( WRKBL, 3*M + LWORK_DGEBRD )
+ WRKBL = MAX( WRKBL, 3*M + LWORK_DORGBR_P )
+ WRKBL = MAX( WRKBL, 3*M + LWORK_DORGBR_Q )
WRKBL = MAX( WRKBL, BDSPAC )
MAXWRK = 2*M*M + WRKBL
- MINWRK = MAX( 3*M+N, BDSPAC )
+ MINWRK = MAX( 3*M + N, BDSPAC )
ELSE IF( WNTVA .AND. WNTUAS ) THEN
*
* Path 9t(N much larger than M, JOBU='S' or 'A',
* JOBVT='A')
*
- WRKBL = M + M*ILAENV( 1, 'DGELQF', ' ', M, N, -1, -1 )
- WRKBL = MAX( WRKBL, M+N*ILAENV( 1, 'DORGLQ', ' ', N,
- $ N, M, -1 ) )
- WRKBL = MAX( WRKBL, 3*M+2*M*
- $ ILAENV( 1, 'DGEBRD', ' ', M, M, -1, -1 ) )
- WRKBL = MAX( WRKBL, 3*M+( M-1 )*
- $ ILAENV( 1, 'DORGBR', 'P', M, M, M, -1 ) )
- WRKBL = MAX( WRKBL, 3*M+M*
- $ ILAENV( 1, 'DORGBR', 'Q', M, M, M, -1 ) )
+ WRKBL = M + LWORK_DGELQF
+ WRKBL = MAX( WRKBL, M + LWORK_DORGLQ_N )
+ WRKBL = MAX( WRKBL, 3*M + LWORK_DGEBRD )
+ WRKBL = MAX( WRKBL, 3*M + LWORK_DORGBR_P )
+ WRKBL = MAX( WRKBL, 3*M + LWORK_DORGBR_Q )
WRKBL = MAX( WRKBL, BDSPAC )
MAXWRK = M*M + WRKBL
- MINWRK = MAX( 3*M+N, BDSPAC )
+ MINWRK = MAX( 3*M + N, BDSPAC )
END IF
ELSE
*
* Path 10t(N greater than M, but not much larger)
*
- MAXWRK = 3*M + ( M+N )*ILAENV( 1, 'DGEBRD', ' ', M, N,
- $ -1, -1 )
- IF( WNTVS .OR. WNTVO )
- $ MAXWRK = MAX( MAXWRK, 3*M+M*
- $ ILAENV( 1, 'DORGBR', 'P', M, N, M, -1 ) )
- IF( WNTVA )
- $ MAXWRK = MAX( MAXWRK, 3*M+N*
- $ ILAENV( 1, 'DORGBR', 'P', N, N, M, -1 ) )
- IF( .NOT.WNTUN )
- $ MAXWRK = MAX( MAXWRK, 3*M+( M-1 )*
- $ ILAENV( 1, 'DORGBR', 'Q', M, M, M, -1 ) )
+ CALL DGEBRD( M, N, A, LDA, S, DUM(1), DUM(1),
+ $ DUM(1), DUM(1), -1, IERR )
+ LWORK_DGEBRD = INT( DUM(1) )
+ MAXWRK = 3*M + LWORK_DGEBRD
+ IF( WNTVS .OR. WNTVO ) THEN
+* Compute space needed for DORGBR P
+ CALL DORGBR( 'P', M, N, M, A, N, DUM(1),
+ $ DUM(1), -1, IERR )
+ LWORK_DORGBR_P = INT( DUM(1) )
+ MAXWRK = MAX( MAXWRK, 3*M + LWORK_DORGBR_P )
+ END IF
+ IF( WNTVA ) THEN
+ CALL DORGBR( 'P', N, N, M, A, N, DUM(1),
+ $ DUM(1), -1, IERR )
+ LWORK_DORGBR_P = INT( DUM(1) )
+ MAXWRK = MAX( MAXWRK, 3*M + LWORK_DORGBR_P )
+ END IF
+ IF( .NOT.WNTUN ) THEN
+ MAXWRK = MAX( MAXWRK, 3*M + LWORK_DORGBR_Q )
+ END IF
MAXWRK = MAX( MAXWRK, BDSPAC )
- MINWRK = MAX( 3*M+N, BDSPAC )
+ MINWRK = MAX( 3*M + N, BDSPAC )
END IF
END IF
MAXWRK = MAX( MAXWRK, MINWRK )
@@ -597,21 +682,24 @@
IWORK = ITAU + N
*
* Compute A=Q*R
-* (Workspace: need 2*N, prefer N+N*NB)
+* (Workspace: need 2*N, prefer N + N*NB)
*
CALL DGEQRF( M, N, A, LDA, WORK( ITAU ), WORK( IWORK ),
$ LWORK-IWORK+1, IERR )
*
* 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
ITAUQ = IE + N
ITAUP = ITAUQ + N
IWORK = ITAUP + N
*
* Bidiagonalize R in A
-* (Workspace: need 4*N, prefer 3*N+2*N*NB)
+* (Workspace: need 4*N, prefer 3*N + 2*N*NB)
*
CALL DGEBRD( N, N, A, LDA, S, WORK( IE ), WORK( ITAUQ ),
$ WORK( ITAUP ), WORK( IWORK ), LWORK-IWORK+1,
@@ -620,7 +708,7 @@
IF( WNTVO .OR. WNTVAS ) THEN
*
* If right singular vectors desired, generate P'.
-* (Workspace: need 4*N-1, prefer 3*N+(N-1)*NB)
+* (Workspace: need 4*N-1, prefer 3*N + (N-1)*NB)
*
CALL DORGBR( 'P', N, N, N, A, LDA, WORK( ITAUP ),
$ WORK( IWORK ), LWORK-IWORK+1, IERR )
@@ -651,13 +739,13 @@
* Sufficient workspace for a fast algorithm
*
IR = 1
- IF( LWORK.GE.MAX( WRKBL, LDA*N+N )+LDA*N ) THEN
+ IF( LWORK.GE.MAX( WRKBL, LDA*N + N ) + LDA*N ) THEN
*
* WORK(IU) is LDA by N, WORK(IR) is LDA by N
*
LDWRKU = LDA
LDWRKR = LDA
- ELSE IF( LWORK.GE.MAX( WRKBL, LDA*N+N )+N*N ) THEN
+ ELSE IF( LWORK.GE.MAX( WRKBL, LDA*N + N ) + N*N ) THEN
*
* WORK(IU) is LDA by N, WORK(IR) is N by N
*
@@ -674,7 +762,7 @@
IWORK = ITAU + N
*
* Compute A=Q*R
-* (Workspace: need N*N+2*N, prefer N*N+N+N*NB)
+* (Workspace: need N*N + 2*N, prefer N*N + N + N*NB)
*
CALL DGEQRF( M, N, A, LDA, WORK( ITAU ),
$ WORK( IWORK ), LWORK-IWORK+1, IERR )
@@ -686,7 +774,7 @@
$ LDWRKR )
*
* Generate Q in A
-* (Workspace: need N*N+2*N, prefer N*N+N+N*NB)
+* (Workspace: need N*N + 2*N, prefer N*N + N + N*NB)
*
CALL DORGQR( M, N, N, A, LDA, WORK( ITAU ),
$ WORK( IWORK ), LWORK-IWORK+1, IERR )
@@ -696,14 +784,14 @@
IWORK = ITAUP + N
*
* Bidiagonalize R in WORK(IR)
-* (Workspace: need N*N+4*N, prefer N*N+3*N+2*N*NB)
+* (Workspace: need N*N + 4*N, prefer N*N + 3*N + 2*N*NB)
*
CALL DGEBRD( N, N, WORK( IR ), LDWRKR, S, WORK( IE ),
$ WORK( ITAUQ ), WORK( ITAUP ),
$ WORK( IWORK ), LWORK-IWORK+1, IERR )
*
* Generate left vectors bidiagonalizing R
-* (Workspace: need N*N+4*N, prefer N*N+3*N+N*NB)
+* (Workspace: need N*N + 4*N, prefer N*N + 3*N + N*NB)
*
CALL DORGBR( 'Q', N, N, N, WORK( IR ), LDWRKR,
$ WORK( ITAUQ ), WORK( IWORK ),
@@ -712,7 +800,7 @@
*
* Perform bidiagonal QR iteration, computing left
* 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,
$ WORK( IR ), LDWRKR, DUM, 1,
@@ -721,7 +809,7 @@
*
* Multiply Q in A by left singular vectors of R in
* WORK(IR), storing result in WORK(IU) and copying to A
-* (Workspace: need N*N+2*N, prefer N*N+M*N+N)
+* (Workspace: need N*N + 2*N, prefer N*N + M*N + N)
*
DO 10 I = 1, M, LDWRKU
CHUNK = MIN( M-I+1, LDWRKU )
@@ -742,14 +830,14 @@
IWORK = ITAUP + N
*
* Bidiagonalize A
-* (Workspace: need 3*N+M, prefer 3*N+(M+N)*NB)
+* (Workspace: need 3*N + M, prefer 3*N + (M + N)*NB)
*
CALL DGEBRD( M, N, A, LDA, S, WORK( IE ),
$ WORK( ITAUQ ), WORK( ITAUP ),
$ WORK( IWORK ), LWORK-IWORK+1, IERR )
*
* Generate left vectors bidiagonalizing A
-* (Workspace: need 4*N, prefer 3*N+N*NB)
+* (Workspace: need 4*N, prefer 3*N + N*NB)
*
CALL DORGBR( 'Q', M, N, N, A, LDA, WORK( ITAUQ ),
$ WORK( IWORK ), LWORK-IWORK+1, IERR )
@@ -775,13 +863,13 @@
* Sufficient workspace for a fast algorithm
*
IR = 1
- IF( LWORK.GE.MAX( WRKBL, LDA*N+N )+LDA*N ) THEN
+ IF( LWORK.GE.MAX( WRKBL, LDA*N + N ) + LDA*N ) THEN
*
* WORK(IU) is LDA by N and WORK(IR) is LDA by N
*
LDWRKU = LDA
LDWRKR = LDA
- ELSE IF( LWORK.GE.MAX( WRKBL, LDA*N+N )+N*N ) THEN
+ ELSE IF( LWORK.GE.MAX( WRKBL, LDA*N + N ) + N*N ) THEN
*
* WORK(IU) is LDA by N and WORK(IR) is N by N
*
@@ -798,7 +886,7 @@
IWORK = ITAU + N
*
* Compute A=Q*R
-* (Workspace: need N*N+2*N, prefer N*N+N+N*NB)
+* (Workspace: need N*N + 2*N, prefer N*N + N + N*NB)
*
CALL DGEQRF( M, N, A, LDA, WORK( ITAU ),
$ WORK( IWORK ), LWORK-IWORK+1, IERR )
@@ -811,7 +899,7 @@
$ VT( 2, 1 ), LDVT )
*
* Generate Q in A
-* (Workspace: need N*N+2*N, prefer N*N+N+N*NB)
+* (Workspace: need N*N + 2*N, prefer N*N + N + N*NB)
*
CALL DORGQR( M, N, N, A, LDA, WORK( ITAU ),
$ WORK( IWORK ), LWORK-IWORK+1, IERR )
@@ -821,7 +909,7 @@
IWORK = ITAUP + N
*
* Bidiagonalize R in VT, copying result to WORK(IR)
-* (Workspace: need N*N+4*N, prefer N*N+3*N+2*N*NB)
+* (Workspace: need N*N + 4*N, prefer N*N + 3*N + 2*N*NB)
*
CALL DGEBRD( N, N, VT, LDVT, S, WORK( IE ),
$ WORK( ITAUQ ), WORK( ITAUP ),
@@ -829,14 +917,14 @@
CALL DLACPY( 'L', N, N, VT, LDVT, WORK( IR ), LDWRKR )
*
* Generate left vectors bidiagonalizing R in WORK(IR)
-* (Workspace: need N*N+4*N, prefer N*N+3*N+N*NB)
+* (Workspace: need N*N + 4*N, prefer N*N + 3*N + N*NB)
*
CALL DORGBR( 'Q', N, N, N, WORK( IR ), LDWRKR,
$ WORK( ITAUQ ), WORK( IWORK ),
$ LWORK-IWORK+1, IERR )
*
* Generate right vectors bidiagonalizing R in VT
-* (Workspace: need N*N+4*N-1, prefer N*N+3*N+(N-1)*NB)
+* (Workspace: need N*N + 4*N-1, prefer N*N + 3*N + (N-1)*NB)
*
CALL DORGBR( 'P', N, N, N, VT, LDVT, WORK( ITAUP ),
$ WORK( IWORK ), LWORK-IWORK+1, IERR )
@@ -845,7 +933,7 @@
* Perform bidiagonal QR iteration, computing left
* singular vectors of R in WORK(IR) and computing 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, LDVT,
$ WORK( IR ), LDWRKR, DUM, 1,
@@ -854,7 +942,7 @@
*
* Multiply Q in A by left singular vectors of R in
* WORK(IR), storing result in WORK(IU) and copying to A
-* (Workspace: need N*N+2*N, prefer N*N+M*N+N)
+* (Workspace: need N*N + 2*N, prefer N*N + M*N + N)
*
DO 20 I = 1, M, LDWRKU
CHUNK = MIN( M-I+1, LDWRKU )
@@ -873,7 +961,7 @@
IWORK = ITAU + N
*
* Compute A=Q*R
-* (Workspace: need 2*N, prefer N+N*NB)
+* (Workspace: need 2*N, prefer N + N*NB)
*
CALL DGEQRF( M, N, A, LDA, WORK( ITAU ),
$ WORK( IWORK ), LWORK-IWORK+1, IERR )
@@ -886,7 +974,7 @@
$ VT( 2, 1 ), LDVT )
*
* Generate Q in A
-* (Workspace: need 2*N, prefer N+N*NB)
+* (Workspace: need 2*N, prefer N + N*NB)
*
CALL DORGQR( M, N, N, A, LDA, WORK( ITAU ),
$ WORK( IWORK ), LWORK-IWORK+1, IERR )
@@ -896,21 +984,21 @@
IWORK = ITAUP + N
*
* Bidiagonalize R in VT
-* (Workspace: need 4*N, prefer 3*N+2*N*NB)
+* (Workspace: need 4*N, prefer 3*N + 2*N*NB)
*
CALL DGEBRD( N, N, VT, LDVT, S, WORK( IE ),
$ WORK( ITAUQ ), WORK( ITAUP ),
$ WORK( IWORK ), LWORK-IWORK+1, IERR )
*
* Multiply Q in A by left vectors bidiagonalizing R
-* (Workspace: need 3*N+M, prefer 3*N+M*NB)
+* (Workspace: need 3*N + M, prefer 3*N + M*NB)
*
CALL DORMBR( 'Q', 'R', 'N', M, N, N, VT, LDVT,
$ WORK( ITAUQ ), A, LDA, WORK( IWORK ),
$ LWORK-IWORK+1, IERR )
*
* Generate right vectors bidiagonalizing R in VT
-* (Workspace: need 4*N-1, prefer 3*N+(N-1)*NB)
+* (Workspace: need 4*N-1, prefer 3*N + (N-1)*NB)
*
CALL DORGBR( 'P', N, N, N, VT, LDVT, WORK( ITAUP ),
$ WORK( IWORK ), LWORK-IWORK+1, IERR )
@@ -954,7 +1042,7 @@
IWORK = ITAU + N
*
* Compute A=Q*R
-* (Workspace: need N*N+2*N, prefer N*N+N+N*NB)
+* (Workspace: need N*N + 2*N, prefer N*N + N + N*NB)
*
CALL DGEQRF( M, N, A, LDA, WORK( ITAU ),
$ WORK( IWORK ), LWORK-IWORK+1, IERR )
@@ -967,7 +1055,7 @@
$ WORK( IR+1 ), LDWRKR )
*
* Generate Q in A
-* (Workspace: need N*N+2*N, prefer N*N+N+N*NB)
+* (Workspace: need N*N + 2*N, prefer N*N + N + N*NB)
*
CALL DORGQR( M, N, N, A, LDA, WORK( ITAU ),
$ WORK( IWORK ), LWORK-IWORK+1, IERR )
@@ -977,7 +1065,7 @@
IWORK = ITAUP + N
*
* Bidiagonalize R in WORK(IR)
-* (Workspace: need N*N+4*N, prefer N*N+3*N+2*N*NB)
+* (Workspace: need N*N + 4*N, prefer N*N + 3*N + 2*N*NB)
*
CALL DGEBRD( N, N, WORK( IR ), LDWRKR, S,
$ WORK( IE ), WORK( ITAUQ ),
@@ -985,7 +1073,7 @@
$ LWORK-IWORK+1, IERR )
*
* Generate left vectors bidiagonalizing R in WORK(IR)
-* (Workspace: need N*N+4*N, prefer N*N+3*N+N*NB)
+* (Workspace: need N*N + 4*N, prefer N*N + 3*N + N*NB)
*
CALL DORGBR( 'Q', N, N, N, WORK( IR ), LDWRKR,
$ WORK( ITAUQ ), WORK( IWORK ),
@@ -994,7 +1082,7 @@
*
* Perform bidiagonal QR iteration, computing left
* 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, WORK( IR ), LDWRKR, DUM, 1,
@@ -1015,14 +1103,14 @@
IWORK = ITAU + N
*
* Compute A=Q*R, copying result to U
-* (Workspace: need 2*N, prefer N+N*NB)
+* (Workspace: need 2*N, prefer N + N*NB)
*
CALL DGEQRF( M, N, A, LDA, WORK( ITAU ),
$ WORK( IWORK ), LWORK-IWORK+1, IERR )
CALL DLACPY( 'L', M, N, A, LDA, U, LDU )
*
* Generate Q in U
-* (Workspace: need 2*N, prefer N+N*NB)
+* (Workspace: need 2*N, prefer N + N*NB)
*
CALL DORGQR( M, N, N, U, LDU, WORK( ITAU ),
$ WORK( IWORK ), LWORK-IWORK+1, IERR )
@@ -1033,18 +1121,20 @@
*
* Zero out below R in A
*
- 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
*
* Bidiagonalize R in A
-* (Workspace: need 4*N, prefer 3*N+2*N*NB)
+* (Workspace: need 4*N, prefer 3*N + 2*N*NB)
*
CALL DGEBRD( N, N, A, LDA, S, WORK( IE ),
$ WORK( ITAUQ ), WORK( ITAUP ),
$ WORK( IWORK ), LWORK-IWORK+1, IERR )
*
* Multiply Q in U by left vectors bidiagonalizing R
-* (Workspace: need 3*N+M, prefer 3*N+M*NB)
+* (Workspace: need 3*N + M, prefer 3*N + M*NB)
*
CALL DORMBR( 'Q', 'R', 'N', M, N, N, A, LDA,
$ WORK( ITAUQ ), U, LDU, WORK( IWORK ),
@@ -1079,7 +1169,7 @@
LDWRKU = LDA
IR = IU + LDWRKU*N
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
*
@@ -1098,7 +1188,7 @@
IWORK = ITAU + N
*
* Compute A=Q*R
-* (Workspace: need 2*N*N+2*N, prefer 2*N*N+N+N*NB)
+* (Workspace: need 2*N*N + 2*N, prefer 2*N*N + N + N*NB)
*
CALL DGEQRF( M, N, A, LDA, WORK( ITAU ),
$ WORK( IWORK ), LWORK-IWORK+1, IERR )
@@ -1111,7 +1201,7 @@
$ WORK( IU+1 ), LDWRKU )
*
* Generate Q in A
-* (Workspace: need 2*N*N+2*N, prefer 2*N*N+N+N*NB)
+* (Workspace: need 2*N*N + 2*N, prefer 2*N*N + N + N*NB)
*
CALL DORGQR( M, N, N, A, LDA, WORK( ITAU ),
$ WORK( IWORK ), LWORK-IWORK+1, IERR )
@@ -1122,7 +1212,7 @@
*
* Bidiagonalize R in WORK(IU), copying result to
* WORK(IR)
-* (Workspace: need 2*N*N+4*N,
+* (Workspace: need 2*N*N + 4*N,
* prefer 2*N*N+3*N+2*N*NB)
*
CALL DGEBRD( N, N, WORK( IU ), LDWRKU, S,
@@ -1133,14 +1223,14 @@
$ WORK( IR ), LDWRKR )
*
* Generate left bidiagonalizing vectors in WORK(IU)
-* (Workspace: need 2*N*N+4*N, prefer 2*N*N+3*N+N*NB)
+* (Workspace: need 2*N*N + 4*N, prefer 2*N*N + 3*N + N*NB)
*
CALL DORGBR( 'Q', N, N, N, WORK( IU ), LDWRKU,
$ WORK( ITAUQ ), WORK( IWORK ),
$ LWORK-IWORK+1, IERR )
*
* Generate right bidiagonalizing vectors in WORK(IR)
-* (Workspace: need 2*N*N+4*N-1,
+* (Workspace: need 2*N*N + 4*N-1,
* prefer 2*N*N+3*N+(N-1)*NB)
*
CALL DORGBR( 'P', N, N, N, WORK( IR ), LDWRKR,
@@ -1151,7 +1241,7 @@
* Perform bidiagonal QR iteration, computing left
* singular vectors of R in WORK(IU) and computing
* right singular vectors of R in WORK(IR)
-* (Workspace: need 2*N*N+BDSPAC)
+* (Workspace: need 2*N*N + BDSPAC)
*
CALL DBDSQR( 'U', N, N, N, 0, S, WORK( IE ),
$ WORK( IR ), LDWRKR, WORK( IU ),
@@ -1178,14 +1268,14 @@
IWORK = ITAU + N
*
* Compute A=Q*R, copying result to U
-* (Workspace: need 2*N, prefer N+N*NB)
+* (Workspace: need 2*N, prefer N + N*NB)
*
CALL DGEQRF( M, N, A, LDA, WORK( ITAU ),
$ WORK( IWORK ), LWORK-IWORK+1, IERR )
CALL DLACPY( 'L', M, N, A, LDA, U, LDU )
*
* Generate Q in U
-* (Workspace: need 2*N, prefer N+N*NB)
+* (Workspace: need 2*N, prefer N + N*NB)
*
CALL DORGQR( M, N, N, U, LDU, WORK( ITAU ),
$ WORK( IWORK ), LWORK-IWORK+1, IERR )
@@ -1196,25 +1286,27 @@
*
* Zero out below R in A
*
- 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
*
* Bidiagonalize R in A
-* (Workspace: need 4*N, prefer 3*N+2*N*NB)
+* (Workspace: need 4*N, prefer 3*N + 2*N*NB)
*
CALL DGEBRD( N, N, A, LDA, S, WORK( IE ),
$ WORK( ITAUQ ), WORK( ITAUP ),
$ WORK( IWORK ), LWORK-IWORK+1, IERR )
*
* Multiply Q in U by left vectors bidiagonalizing R
-* (Workspace: need 3*N+M, prefer 3*N+M*NB)
+* (Workspace: need 3*N + M, prefer 3*N + M*NB)
*
CALL DORMBR( 'Q', 'R', 'N', M, N, N, A, LDA,
$ WORK( ITAUQ ), U, LDU, WORK( IWORK ),
$ LWORK-IWORK+1, IERR )
*
* Generate right vectors bidiagonalizing R in A
-* (Workspace: need 4*N-1, prefer 3*N+(N-1)*NB)
+* (Workspace: need 4*N-1, prefer 3*N + (N-1)*NB)
*
CALL DORGBR( 'P', N, N, N, A, LDA, WORK( ITAUP ),
$ WORK( IWORK ), LWORK-IWORK+1, IERR )
@@ -1258,7 +1350,7 @@
IWORK = ITAU + N
*
* Compute A=Q*R
-* (Workspace: need N*N+2*N, prefer N*N+N+N*NB)
+* (Workspace: need N*N + 2*N, prefer N*N + N + N*NB)
*
CALL DGEQRF( M, N, A, LDA, WORK( ITAU ),
$ WORK( IWORK ), LWORK-IWORK+1, IERR )
@@ -1271,7 +1363,7 @@
$ WORK( IU+1 ), LDWRKU )
*
* Generate Q in A
-* (Workspace: need N*N+2*N, prefer N*N+N+N*NB)
+* (Workspace: need N*N + 2*N, prefer N*N + N + N*NB)
*
CALL DORGQR( M, N, N, A, LDA, WORK( ITAU ),
$ WORK( IWORK ), LWORK-IWORK+1, IERR )
@@ -1281,7 +1373,7 @@
IWORK = ITAUP + N
*
* Bidiagonalize R in WORK(IU), copying result to VT
-* (Workspace: need N*N+4*N, prefer N*N+3*N+2*N*NB)
+* (Workspace: need N*N + 4*N, prefer N*N + 3*N + 2*N*NB)
*
CALL DGEBRD( N, N, WORK( IU ), LDWRKU, S,
$ WORK( IE ), WORK( ITAUQ ),
@@ -1291,14 +1383,14 @@
$ LDVT )
*
* Generate left bidiagonalizing vectors in WORK(IU)
-* (Workspace: need N*N+4*N, prefer N*N+3*N+N*NB)
+* (Workspace: need N*N + 4*N, prefer N*N + 3*N + N*NB)
*
CALL DORGBR( 'Q', N, N, N, WORK( IU ), LDWRKU,
$ WORK( ITAUQ ), WORK( IWORK ),
$ LWORK-IWORK+1, IERR )
*
* Generate right bidiagonalizing vectors in VT
-* (Workspace: need N*N+4*N-1,
+* (Workspace: need N*N + 4*N-1,
* prefer N*N+3*N+(N-1)*NB)
*
CALL DORGBR( 'P', N, N, N, VT, LDVT, WORK( ITAUP ),
@@ -1308,7 +1400,7 @@
* Perform bidiagonal QR iteration, computing left
* singular vectors of R in WORK(IU) and computing
* 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,
$ LDVT, WORK( IU ), LDWRKU, DUM, 1,
@@ -1329,14 +1421,14 @@
IWORK = ITAU + N
*
* Compute A=Q*R, copying result to U
-* (Workspace: need 2*N, prefer N+N*NB)
+* (Workspace: need 2*N, prefer N + N*NB)
*
CALL DGEQRF( M, N, A, LDA, WORK( ITAU ),
$ WORK( IWORK ), LWORK-IWORK+1, IERR )
CALL DLACPY( 'L', M, N, A, LDA, U, LDU )
*
* Generate Q in U
-* (Workspace: need 2*N, prefer N+N*NB)
+* (Workspace: need 2*N, prefer N + N*NB)
*
CALL DORGQR( M, N, N, U, LDU, WORK( ITAU ),
$ WORK( IWORK ), LWORK-IWORK+1, IERR )
@@ -1353,7 +1445,7 @@
IWORK = ITAUP + N
*
* Bidiagonalize R in VT
-* (Workspace: need 4*N, prefer 3*N+2*N*NB)
+* (Workspace: need 4*N, prefer 3*N + 2*N*NB)
*
CALL DGEBRD( N, N, VT, LDVT, S, WORK( IE ),
$ WORK( ITAUQ ), WORK( ITAUP ),
@@ -1361,14 +1453,14 @@
*
* Multiply Q in U by left bidiagonalizing vectors
* in VT
-* (Workspace: need 3*N+M, prefer 3*N+M*NB)
+* (Workspace: need 3*N + M, prefer 3*N + M*NB)
*
CALL DORMBR( 'Q', 'R', 'N', M, N, N, VT, LDVT,
$ WORK( ITAUQ ), U, LDU, WORK( IWORK ),
$ LWORK-IWORK+1, IERR )
*
* Generate right bidiagonalizing vectors in VT
-* (Workspace: need 4*N-1, prefer 3*N+(N-1)*NB)
+* (Workspace: need 4*N-1, prefer 3*N + (N-1)*NB)
*
CALL DORGBR( 'P', N, N, N, VT, LDVT, WORK( ITAUP ),
$ WORK( IWORK ), LWORK-IWORK+1, IERR )
@@ -1415,7 +1507,7 @@
IWORK = ITAU + N
*
* Compute A=Q*R, copying result to U
-* (Workspace: need N*N+2*N, prefer N*N+N+N*NB)
+* (Workspace: need N*N + 2*N, prefer N*N + N + N*NB)
*
CALL DGEQRF( M, N, A, LDA, WORK( ITAU ),
$ WORK( IWORK ), LWORK-IWORK+1, IERR )
@@ -1429,7 +1521,7 @@
$ WORK( IR+1 ), LDWRKR )
*
* Generate Q in U
-* (Workspace: need N*N+N+M, prefer N*N+N+M*NB)
+* (Workspace: need N*N + N + M, prefer N*N + N + M*NB)
*
CALL DORGQR( M, M, N, U, LDU, WORK( ITAU ),
$ WORK( IWORK ), LWORK-IWORK+1, IERR )
@@ -1439,7 +1531,7 @@
IWORK = ITAUP + N
*
* Bidiagonalize R in WORK(IR)
-* (Workspace: need N*N+4*N, prefer N*N+3*N+2*N*NB)
+* (Workspace: need N*N + 4*N, prefer N*N + 3*N + 2*N*NB)
*
CALL DGEBRD( N, N, WORK( IR ), LDWRKR, S,
$ WORK( IE ), WORK( ITAUQ ),
@@ -1447,7 +1539,7 @@
$ LWORK-IWORK+1, IERR )
*
* Generate left bidiagonalizing vectors in WORK(IR)
-* (Workspace: need N*N+4*N, prefer N*N+3*N+N*NB)
+* (Workspace: need N*N + 4*N, prefer N*N + 3*N + N*NB)
*
CALL DORGBR( 'Q', N, N, N, WORK( IR ), LDWRKR,
$ WORK( ITAUQ ), WORK( IWORK ),
@@ -1456,7 +1548,7 @@
*
* Perform bidiagonal QR iteration, computing left
* 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, WORK( IR ), LDWRKR, DUM, 1,
@@ -1481,14 +1573,14 @@
IWORK = ITAU + N
*
* Compute A=Q*R, copying result to U
-* (Workspace: need 2*N, prefer N+N*NB)
+* (Workspace: need 2*N, prefer N + N*NB)
*
CALL DGEQRF( M, N, A, LDA, WORK( ITAU ),
$ WORK( IWORK ), LWORK-IWORK+1, IERR )
CALL DLACPY( 'L', M, N, A, LDA, U, LDU )
*
* 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 ),
$ WORK( IWORK ), LWORK-IWORK+1, IERR )
@@ -1499,11 +1591,13 @@
*
* Zero out below R in A
*
- 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
*
* Bidiagonalize R in A
-* (Workspace: need 4*N, prefer 3*N+2*N*NB)
+* (Workspace: need 4*N, prefer 3*N + 2*N*NB)
*
CALL DGEBRD( N, N, A, LDA, S, WORK( IE ),
$ WORK( ITAUQ ), WORK( ITAUP ),
@@ -1511,7 +1605,7 @@
*
* Multiply Q in U by left bidiagonalizing vectors
* in A
-* (Workspace: need 3*N+M, prefer 3*N+M*NB)
+* (Workspace: need 3*N + M, prefer 3*N + M*NB)
*
CALL DORMBR( 'Q', 'R', 'N', M, N, N, A, LDA,
$ WORK( ITAUQ ), U, LDU, WORK( IWORK ),
@@ -1546,7 +1640,7 @@
LDWRKU = LDA
IR = IU + LDWRKU*N
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
*
@@ -1565,14 +1659,14 @@
IWORK = ITAU + N
*
* Compute A=Q*R, copying result to U
-* (Workspace: need 2*N*N+2*N, prefer 2*N*N+N+N*NB)
+* (Workspace: need 2*N*N + 2*N, prefer 2*N*N + N + N*NB)
*
CALL DGEQRF( M, N, A, LDA, WORK( ITAU ),
$ WORK( IWORK ), LWORK-IWORK+1, IERR )
CALL DLACPY( 'L', M, N, A, LDA, U, LDU )
*
* Generate Q in U
-* (Workspace: need 2*N*N+N+M, prefer 2*N*N+N+M*NB)
+* (Workspace: need 2*N*N + N + M, prefer 2*N*N + N + M*NB)
*
CALL DORGQR( M, M, N, U, LDU, WORK( ITAU ),
$ WORK( IWORK ), LWORK-IWORK+1, IERR )
@@ -1590,7 +1684,7 @@
*
* Bidiagonalize R in WORK(IU), copying result to
* WORK(IR)
-* (Workspace: need 2*N*N+4*N,
+* (Workspace: need 2*N*N + 4*N,
* prefer 2*N*N+3*N+2*N*NB)
*
CALL DGEBRD( N, N, WORK( IU ), LDWRKU, S,
@@ -1601,14 +1695,14 @@
$ WORK( IR ), LDWRKR )
*
* Generate left bidiagonalizing vectors in WORK(IU)
-* (Workspace: need 2*N*N+4*N, prefer 2*N*N+3*N+N*NB)
+* (Workspace: need 2*N*N + 4*N, prefer 2*N*N + 3*N + N*NB)
*
CALL DORGBR( 'Q', N, N, N, WORK( IU ), LDWRKU,
$ WORK( ITAUQ ), WORK( IWORK ),
$ LWORK-IWORK+1, IERR )
*
* Generate right bidiagonalizing vectors in WORK(IR)
-* (Workspace: need 2*N*N+4*N-1,
+* (Workspace: need 2*N*N + 4*N-1,
* prefer 2*N*N+3*N+(N-1)*NB)
*
CALL DORGBR( 'P', N, N, N, WORK( IR ), LDWRKR,
@@ -1619,7 +1713,7 @@
* Perform bidiagonal QR iteration, computing left
* singular vectors of R in WORK(IU) and computing
* right singular vectors of R in WORK(IR)
-* (Workspace: need 2*N*N+BDSPAC)
+* (Workspace: need 2*N*N + BDSPAC)
*
CALL DBDSQR( 'U', N, N, N, 0, S, WORK( IE ),
$ WORK( IR ), LDWRKR, WORK( IU ),
@@ -1649,14 +1743,14 @@
IWORK = ITAU + N
*
* Compute A=Q*R, copying result to U
-* (Workspace: need 2*N, prefer N+N*NB)
+* (Workspace: need 2*N, prefer N + N*NB)
*
CALL DGEQRF( M, N, A, LDA, WORK( ITAU ),
$ WORK( IWORK ), LWORK-IWORK+1, IERR )
CALL DLACPY( 'L', M, N, A, LDA, U, LDU )
*
* 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 ),
$ WORK( IWORK ), LWORK-IWORK+1, IERR )
@@ -1667,11 +1761,13 @@
*
* Zero out below R in A
*
- 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
*
* Bidiagonalize R in A
-* (Workspace: need 4*N, prefer 3*N+2*N*NB)
+* (Workspace: need 4*N, prefer 3*N + 2*N*NB)
*
CALL DGEBRD( N, N, A, LDA, S, WORK( IE ),
$ WORK( ITAUQ ), WORK( ITAUP ),
@@ -1679,14 +1775,14 @@
*
* Multiply Q in U by left bidiagonalizing vectors
* in A
-* (Workspace: need 3*N+M, prefer 3*N+M*NB)
+* (Workspace: need 3*N + M, prefer 3*N + M*NB)
*
CALL DORMBR( 'Q', 'R', 'N', M, N, N, A, LDA,
$ WORK( ITAUQ ), U, LDU, WORK( IWORK ),
$ LWORK-IWORK+1, IERR )
*
* Generate right bidiagonalizing vectors in A
-* (Workspace: need 4*N-1, prefer 3*N+(N-1)*NB)
+* (Workspace: need 4*N-1, prefer 3*N + (N-1)*NB)
*
CALL DORGBR( 'P', N, N, N, A, LDA, WORK( ITAUP ),
$ WORK( IWORK ), LWORK-IWORK+1, IERR )
@@ -1730,14 +1826,14 @@
IWORK = ITAU + N
*
* Compute A=Q*R, copying result to U
-* (Workspace: need N*N+2*N, prefer N*N+N+N*NB)
+* (Workspace: need N*N + 2*N, prefer N*N + N + N*NB)
*
CALL DGEQRF( M, N, A, LDA, WORK( ITAU ),
$ WORK( IWORK ), LWORK-IWORK+1, IERR )
CALL DLACPY( 'L', M, N, A, LDA, U, LDU )
*
* Generate Q in U
-* (Workspace: need N*N+N+M, prefer N*N+N+M*NB)
+* (Workspace: need N*N + N + M, prefer N*N + N + M*NB)
*
CALL DORGQR( M, M, N, U, LDU, WORK( ITAU ),
$ WORK( IWORK ), LWORK-IWORK+1, IERR )
@@ -1754,7 +1850,7 @@
IWORK = ITAUP + N
*
* Bidiagonalize R in WORK(IU), copying result to VT
-* (Workspace: need N*N+4*N, prefer N*N+3*N+2*N*NB)
+* (Workspace: need N*N + 4*N, prefer N*N + 3*N + 2*N*NB)
*
CALL DGEBRD( N, N, WORK( IU ), LDWRKU, S,
$ WORK( IE ), WORK( ITAUQ ),
@@ -1764,14 +1860,14 @@
$ LDVT )
*
* Generate left bidiagonalizing vectors in WORK(IU)
-* (Workspace: need N*N+4*N, prefer N*N+3*N+N*NB)
+* (Workspace: need N*N + 4*N, prefer N*N + 3*N + N*NB)
*
CALL DORGBR( 'Q', N, N, N, WORK( IU ), LDWRKU,
$ WORK( ITAUQ ), WORK( IWORK ),
$ LWORK-IWORK+1, IERR )
*
* Generate right bidiagonalizing vectors in VT
-* (Workspace: need N*N+4*N-1,
+* (Workspace: need N*N + 4*N-1,
* prefer N*N+3*N+(N-1)*NB)
*
CALL DORGBR( 'P', N, N, N, VT, LDVT, WORK( ITAUP ),
@@ -1781,7 +1877,7 @@
* Perform bidiagonal QR iteration, computing left
* singular vectors of R in WORK(IU) and computing
* 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,
$ LDVT, WORK( IU ), LDWRKU, DUM, 1,
@@ -1806,14 +1902,14 @@
IWORK = ITAU + N
*
* Compute A=Q*R, copying result to U
-* (Workspace: need 2*N, prefer N+N*NB)
+* (Workspace: need 2*N, prefer N + N*NB)
*
CALL DGEQRF( M, N, A, LDA, WORK( ITAU ),
$ WORK( IWORK ), LWORK-IWORK+1, IERR )
CALL DLACPY( 'L', M, N, A, LDA, U, LDU )
*
* 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 ),
$ WORK( IWORK ), LWORK-IWORK+1, IERR )
@@ -1830,7 +1926,7 @@
IWORK = ITAUP + N
*
* Bidiagonalize R in VT
-* (Workspace: need 4*N, prefer 3*N+2*N*NB)
+* (Workspace: need 4*N, prefer 3*N + 2*N*NB)
*
CALL DGEBRD( N, N, VT, LDVT, S, WORK( IE ),
$ WORK( ITAUQ ), WORK( ITAUP ),
@@ -1838,14 +1934,14 @@
*
* Multiply Q in U by left bidiagonalizing vectors
* in VT
-* (Workspace: need 3*N+M, prefer 3*N+M*NB)
+* (Workspace: need 3*N + M, prefer 3*N + M*NB)
*
CALL DORMBR( 'Q', 'R', 'N', M, N, N, VT, LDVT,
$ WORK( ITAUQ ), U, LDU, WORK( IWORK ),
$ LWORK-IWORK+1, IERR )
*
* Generate right bidiagonalizing vectors in VT
-* (Workspace: need 4*N-1, prefer 3*N+(N-1)*NB)
+* (Workspace: need 4*N-1, prefer 3*N + (N-1)*NB)
*
CALL DORGBR( 'P', N, N, N, VT, LDVT, WORK( ITAUP ),
$ WORK( IWORK ), LWORK-IWORK+1, IERR )
@@ -1879,7 +1975,7 @@
IWORK = ITAUP + N
*
* Bidiagonalize A
-* (Workspace: need 3*N+M, prefer 3*N+(M+N)*NB)
+* (Workspace: need 3*N + M, prefer 3*N + (M + N)*NB)
*
CALL DGEBRD( M, N, A, LDA, S, WORK( IE ), WORK( ITAUQ ),
$ WORK( ITAUP ), WORK( IWORK ), LWORK-IWORK+1,
@@ -1888,7 +1984,7 @@
*
* If left singular vectors desired in U, copy result to U
* and generate left bidiagonalizing vectors in U
-* (Workspace: need 3*N+NCU, prefer 3*N+NCU*NB)
+* (Workspace: need 3*N + NCU, prefer 3*N + NCU*NB)
*
CALL DLACPY( 'L', M, N, A, LDA, U, LDU )
IF( WNTUS )
@@ -1902,7 +1998,7 @@
*
* If right singular vectors desired in VT, copy result to
* VT and generate right bidiagonalizing vectors in VT
-* (Workspace: need 4*N-1, prefer 3*N+(N-1)*NB)
+* (Workspace: need 4*N-1, prefer 3*N + (N-1)*NB)
*
CALL DLACPY( 'U', N, N, A, LDA, VT, LDVT )
CALL DORGBR( 'P', N, N, N, VT, LDVT, WORK( ITAUP ),
@@ -1912,7 +2008,7 @@
*
* If left singular vectors desired in A, generate left
* bidiagonalizing vectors in A
-* (Workspace: need 4*N, prefer 3*N+N*NB)
+* (Workspace: need 4*N, prefer 3*N + N*NB)
*
CALL DORGBR( 'Q', M, N, N, A, LDA, WORK( ITAUQ ),
$ WORK( IWORK ), LWORK-IWORK+1, IERR )
@@ -1921,7 +2017,7 @@
*
* If right singular vectors desired in A, generate right
* bidiagonalizing vectors in A
-* (Workspace: need 4*N-1, prefer 3*N+(N-1)*NB)
+* (Workspace: need 4*N-1, prefer 3*N + (N-1)*NB)
*
CALL DORGBR( 'P', N, N, N, A, LDA, WORK( ITAUP ),
$ WORK( IWORK ), LWORK-IWORK+1, IERR )
@@ -1983,7 +2079,7 @@
IWORK = ITAU + M
*
* Compute A=L*Q
-* (Workspace: need 2*M, prefer M+M*NB)
+* (Workspace: need 2*M, prefer M + M*NB)
*
CALL DGELQF( M, N, A, LDA, WORK( ITAU ), WORK( IWORK ),
$ LWORK-IWORK+1, IERR )
@@ -1997,7 +2093,7 @@
IWORK = ITAUP + M
*
* Bidiagonalize L in A
-* (Workspace: need 4*M, prefer 3*M+2*M*NB)
+* (Workspace: need 4*M, prefer 3*M + 2*M*NB)
*
CALL DGEBRD( M, M, A, LDA, S, WORK( IE ), WORK( ITAUQ ),
$ WORK( ITAUP ), WORK( IWORK ), LWORK-IWORK+1,
@@ -2005,7 +2101,7 @@
IF( WNTUO .OR. WNTUAS ) THEN
*
* If left singular vectors desired, generate Q
-* (Workspace: need 4*M, prefer 3*M+M*NB)
+* (Workspace: need 4*M, prefer 3*M + M*NB)
*
CALL DORGBR( 'Q', M, M, M, A, LDA, WORK( ITAUQ ),
$ WORK( IWORK ), LWORK-IWORK+1, IERR )
@@ -2038,14 +2134,14 @@
* Sufficient workspace for a fast algorithm
*
IR = 1
- IF( LWORK.GE.MAX( WRKBL, LDA*N+M )+LDA*M ) THEN
+ IF( LWORK.GE.MAX( WRKBL, LDA*N + M ) + LDA*M ) THEN
*
* WORK(IU) is LDA by N and WORK(IR) is LDA by M
*
LDWRKU = LDA
CHUNK = N
LDWRKR = LDA
- ELSE IF( LWORK.GE.MAX( WRKBL, LDA*N+M )+M*M ) THEN
+ ELSE IF( LWORK.GE.MAX( WRKBL, LDA*N + M ) + M*M ) THEN
*
* WORK(IU) is LDA by N and WORK(IR) is M by M
*
@@ -2064,7 +2160,7 @@
IWORK = ITAU + M
*
* Compute A=L*Q
-* (Workspace: need M*M+2*M, prefer M*M+M+M*NB)
+* (Workspace: need M*M + 2*M, prefer M*M + M + M*NB)
*
CALL DGELQF( M, N, A, LDA, WORK( ITAU ),
$ WORK( IWORK ), LWORK-IWORK+1, IERR )
@@ -2076,7 +2172,7 @@
$ WORK( IR+LDWRKR ), LDWRKR )
*
* Generate Q in A
-* (Workspace: need M*M+2*M, prefer M*M+M+M*NB)
+* (Workspace: need M*M + 2*M, prefer M*M + M + M*NB)
*
CALL DORGLQ( M, N, M, A, LDA, WORK( ITAU ),
$ WORK( IWORK ), LWORK-IWORK+1, IERR )
@@ -2086,14 +2182,14 @@
IWORK = ITAUP + M
*
* Bidiagonalize L in WORK(IR)
-* (Workspace: need M*M+4*M, prefer M*M+3*M+2*M*NB)
+* (Workspace: need M*M + 4*M, prefer M*M + 3*M + 2*M*NB)
*
CALL DGEBRD( M, M, WORK( IR ), LDWRKR, S, WORK( IE ),
$ WORK( ITAUQ ), WORK( ITAUP ),
$ WORK( IWORK ), LWORK-IWORK+1, IERR )
*
* Generate right vectors bidiagonalizing L
-* (Workspace: need M*M+4*M-1, prefer M*M+3*M+(M-1)*NB)
+* (Workspace: need M*M + 4*M-1, prefer M*M + 3*M + (M-1)*NB)
*
CALL DORGBR( 'P', M, M, M, WORK( IR ), LDWRKR,
$ WORK( ITAUP ), WORK( IWORK ),
@@ -2102,7 +2198,7 @@
*
* Perform bidiagonal QR iteration, computing right
* 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 ),
$ WORK( IR ), LDWRKR, DUM, 1, DUM, 1,
@@ -2111,7 +2207,7 @@
*
* Multiply right singular vectors of L in WORK(IR) by Q
* in A, storing result in WORK(IU) and copying to A
-* (Workspace: need M*M+2*M, prefer M*M+M*N+M)
+* (Workspace: need M*M + 2*M, prefer M*M + M*N + M)
*
DO 30 I = 1, N, CHUNK
BLK = MIN( N-I+1, CHUNK )
@@ -2132,14 +2228,14 @@
IWORK = ITAUP + M
*
* Bidiagonalize A
-* (Workspace: need 3*M+N, prefer 3*M+(M+N)*NB)
+* (Workspace: need 3*M + N, prefer 3*M + (M + N)*NB)
*
CALL DGEBRD( M, N, A, LDA, S, WORK( IE ),
$ WORK( ITAUQ ), WORK( ITAUP ),
$ WORK( IWORK ), LWORK-IWORK+1, IERR )
*
* Generate right vectors bidiagonalizing A
-* (Workspace: need 4*M, prefer 3*M+M*NB)
+* (Workspace: need 4*M, prefer 3*M + M*NB)
*
CALL DORGBR( 'P', M, N, M, A, LDA, WORK( ITAUP ),
$ WORK( IWORK ), LWORK-IWORK+1, IERR )
@@ -2165,14 +2261,14 @@
* Sufficient workspace for a fast algorithm
*
IR = 1
- IF( LWORK.GE.MAX( WRKBL, LDA*N+M )+LDA*M ) THEN
+ IF( LWORK.GE.MAX( WRKBL, LDA*N + M ) + LDA*M ) THEN
*
* WORK(IU) is LDA by N and WORK(IR) is LDA by M
*
LDWRKU = LDA
CHUNK = N
LDWRKR = LDA
- ELSE IF( LWORK.GE.MAX( WRKBL, LDA*N+M )+M*M ) THEN
+ ELSE IF( LWORK.GE.MAX( WRKBL, LDA*N + M ) + M*M ) THEN
*
* WORK(IU) is LDA by N and WORK(IR) is M by M
*
@@ -2191,7 +2287,7 @@
IWORK = ITAU + M
*
* Compute A=L*Q
-* (Workspace: need M*M+2*M, prefer M*M+M+M*NB)
+* (Workspace: need M*M + 2*M, prefer M*M + M + M*NB)
*
CALL DGELQF( M, N, A, LDA, WORK( ITAU ),
$ WORK( IWORK ), LWORK-IWORK+1, IERR )
@@ -2203,7 +2299,7 @@
$ LDU )
*
* Generate Q in A
-* (Workspace: need M*M+2*M, prefer M*M+M+M*NB)
+* (Workspace: need M*M + 2*M, prefer M*M + M + M*NB)
*
CALL DORGLQ( M, N, M, A, LDA, WORK( ITAU ),
$ WORK( IWORK ), LWORK-IWORK+1, IERR )
@@ -2213,7 +2309,7 @@
IWORK = ITAUP + M
*
* Bidiagonalize L in U, copying result to WORK(IR)
-* (Workspace: need M*M+4*M, prefer M*M+3*M+2*M*NB)
+* (Workspace: need M*M + 4*M, prefer M*M + 3*M + 2*M*NB)
*
CALL DGEBRD( M, M, U, LDU, S, WORK( IE ),
$ WORK( ITAUQ ), WORK( ITAUP ),
@@ -2221,14 +2317,14 @@
CALL DLACPY( 'U', M, M, U, LDU, WORK( IR ), LDWRKR )
*
* Generate right vectors bidiagonalizing L in WORK(IR)
-* (Workspace: need M*M+4*M-1, prefer M*M+3*M+(M-1)*NB)
+* (Workspace: need M*M + 4*M-1, prefer M*M + 3*M + (M-1)*NB)
*
CALL DORGBR( 'P', M, M, M, WORK( IR ), LDWRKR,
$ WORK( ITAUP ), WORK( IWORK ),
$ LWORK-IWORK+1, IERR )
*
* Generate left vectors bidiagonalizing L in U
-* (Workspace: need M*M+4*M, prefer M*M+3*M+M*NB)
+* (Workspace: need M*M + 4*M, prefer M*M + 3*M + M*NB)
*
CALL DORGBR( 'Q', M, M, M, U, LDU, WORK( ITAUQ ),
$ WORK( IWORK ), LWORK-IWORK+1, IERR )
@@ -2237,7 +2333,7 @@
* Perform bidiagonal QR iteration, computing left
* singular vectors of L in U, and computing right
* 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 ),
$ WORK( IR ), LDWRKR, U, LDU, DUM, 1,
@@ -2246,7 +2342,7 @@
*
* Multiply right singular vectors of L in WORK(IR) by Q
* in A, storing result in WORK(IU) and copying to A
-* (Workspace: need M*M+2*M, prefer M*M+M*N+M))
+* (Workspace: need M*M + 2*M, prefer M*M + M*N + M))
*
DO 40 I = 1, N, CHUNK
BLK = MIN( N-I+1, CHUNK )
@@ -2265,7 +2361,7 @@
IWORK = ITAU + M
*
* Compute A=L*Q
-* (Workspace: need 2*M, prefer M+M*NB)
+* (Workspace: need 2*M, prefer M + M*NB)
*
CALL DGELQF( M, N, A, LDA, WORK( ITAU ),
$ WORK( IWORK ), LWORK-IWORK+1, IERR )
@@ -2277,7 +2373,7 @@
$ LDU )
*
* Generate Q in A
-* (Workspace: need 2*M, prefer M+M*NB)
+* (Workspace: need 2*M, prefer M + M*NB)
*
CALL DORGLQ( M, N, M, A, LDA, WORK( ITAU ),
$ WORK( IWORK ), LWORK-IWORK+1, IERR )
@@ -2287,21 +2383,21 @@
IWORK = ITAUP + M
*
* Bidiagonalize L in U
-* (Workspace: need 4*M, prefer 3*M+2*M*NB)
+* (Workspace: need 4*M, prefer 3*M + 2*M*NB)
*
CALL DGEBRD( M, M, U, LDU, S, WORK( IE ),
$ WORK( ITAUQ ), WORK( ITAUP ),
$ WORK( IWORK ), LWORK-IWORK+1, IERR )
*
* Multiply right vectors bidiagonalizing L by Q in A
-* (Workspace: need 3*M+N, prefer 3*M+N*NB)
+* (Workspace: need 3*M + N, prefer 3*M + N*NB)
*
CALL DORMBR( 'P', 'L', 'T', M, N, M, U, LDU,
$ WORK( ITAUP ), A, LDA, WORK( IWORK ),
$ LWORK-IWORK+1, IERR )
*
* Generate left vectors bidiagonalizing L in U
-* (Workspace: need 4*M, prefer 3*M+M*NB)
+* (Workspace: need 4*M, prefer 3*M + M*NB)
*
CALL DORGBR( 'Q', M, M, M, U, LDU, WORK( ITAUQ ),
$ WORK( IWORK ), LWORK-IWORK+1, IERR )
@@ -2345,7 +2441,7 @@
IWORK = ITAU + M
*
* Compute A=L*Q
-* (Workspace: need M*M+2*M, prefer M*M+M+M*NB)
+* (Workspace: need M*M + 2*M, prefer M*M + M + M*NB)
*
CALL DGELQF( M, N, A, LDA, WORK( ITAU ),
$ WORK( IWORK ), LWORK-IWORK+1, IERR )
@@ -2358,7 +2454,7 @@
$ WORK( IR+LDWRKR ), LDWRKR )
*
* Generate Q in A
-* (Workspace: need M*M+2*M, prefer M*M+M+M*NB)
+* (Workspace: need M*M + 2*M, prefer M*M + M + M*NB)
*
CALL DORGLQ( M, N, M, A, LDA, WORK( ITAU ),
$ WORK( IWORK ), LWORK-IWORK+1, IERR )
@@ -2368,7 +2464,7 @@
IWORK = ITAUP + M
*
* Bidiagonalize L in WORK(IR)
-* (Workspace: need M*M+4*M, prefer M*M+3*M+2*M*NB)
+* (Workspace: need M*M + 4*M, prefer M*M + 3*M + 2*M*NB)
*
CALL DGEBRD( M, M, WORK( IR ), LDWRKR, S,
$ WORK( IE ), WORK( ITAUQ ),
@@ -2377,7 +2473,7 @@
*
* Generate right vectors bidiagonalizing L in
* WORK(IR)
-* (Workspace: need M*M+4*M, prefer M*M+3*M+(M-1)*NB)
+* (Workspace: need M*M + 4*M, prefer M*M + 3*M + (M-1)*NB)
*
CALL DORGBR( 'P', M, M, M, WORK( IR ), LDWRKR,
$ WORK( ITAUP ), WORK( IWORK ),
@@ -2386,7 +2482,7 @@
*
* Perform bidiagonal QR iteration, computing right
* 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 ),
$ WORK( IR ), LDWRKR, DUM, 1, DUM, 1,
@@ -2407,7 +2503,7 @@
IWORK = ITAU + M
*
* Compute A=L*Q
-* (Workspace: need 2*M, prefer M+M*NB)
+* (Workspace: need 2*M, prefer M + M*NB)
*
CALL DGELQF( M, N, A, LDA, WORK( ITAU ),
$ WORK( IWORK ), LWORK-IWORK+1, IERR )
@@ -2417,7 +2513,7 @@
CALL DLACPY( 'U', M, N, A, LDA, VT, LDVT )
*
* Generate Q in VT
-* (Workspace: need 2*M, prefer M+M*NB)
+* (Workspace: need 2*M, prefer M + M*NB)
*
CALL DORGLQ( M, N, M, VT, LDVT, WORK( ITAU ),
$ WORK( IWORK ), LWORK-IWORK+1, IERR )
@@ -2432,14 +2528,14 @@
$ LDA )
*
* Bidiagonalize L in A
-* (Workspace: need 4*M, prefer 3*M+2*M*NB)
+* (Workspace: need 4*M, prefer 3*M + 2*M*NB)
*
CALL DGEBRD( M, M, A, LDA, S, WORK( IE ),
$ WORK( ITAUQ ), WORK( ITAUP ),
$ WORK( IWORK ), LWORK-IWORK+1, IERR )
*
* Multiply right vectors bidiagonalizing L by Q in VT
-* (Workspace: need 3*M+N, prefer 3*M+N*NB)
+* (Workspace: need 3*M + N, prefer 3*M + N*NB)
*
CALL DORMBR( 'P', 'L', 'T', M, N, M, A, LDA,
$ WORK( ITAUP ), VT, LDVT,
@@ -2474,7 +2570,7 @@
LDWRKU = LDA
IR = IU + LDWRKU*M
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
*
@@ -2493,7 +2589,7 @@
IWORK = ITAU + M
*
* Compute A=L*Q
-* (Workspace: need 2*M*M+2*M, prefer 2*M*M+M+M*NB)
+* (Workspace: need 2*M*M + 2*M, prefer 2*M*M + M + M*NB)
*
CALL DGELQF( M, N, A, LDA, WORK( ITAU ),
$ WORK( IWORK ), LWORK-IWORK+1, IERR )
@@ -2506,7 +2602,7 @@
$ WORK( IU+LDWRKU ), LDWRKU )
*
* Generate Q in A
-* (Workspace: need 2*M*M+2*M, prefer 2*M*M+M+M*NB)
+* (Workspace: need 2*M*M + 2*M, prefer 2*M*M + M + M*NB)
*
CALL DORGLQ( M, N, M, A, LDA, WORK( ITAU ),
$ WORK( IWORK ), LWORK-IWORK+1, IERR )
@@ -2517,7 +2613,7 @@
*
* Bidiagonalize L in WORK(IU), copying result to
* WORK(IR)
-* (Workspace: need 2*M*M+4*M,
+* (Workspace: need 2*M*M + 4*M,
* prefer 2*M*M+3*M+2*M*NB)
*
CALL DGEBRD( M, M, WORK( IU ), LDWRKU, S,
@@ -2528,7 +2624,7 @@
$ WORK( IR ), LDWRKR )
*
* Generate right bidiagonalizing vectors in WORK(IU)
-* (Workspace: need 2*M*M+4*M-1,
+* (Workspace: need 2*M*M + 4*M-1,
* prefer 2*M*M+3*M+(M-1)*NB)
*
CALL DORGBR( 'P', M, M, M, WORK( IU ), LDWRKU,
@@ -2536,7 +2632,7 @@
$ LWORK-IWORK+1, IERR )
*
* Generate left bidiagonalizing vectors in WORK(IR)
-* (Workspace: need 2*M*M+4*M, prefer 2*M*M+3*M+M*NB)
+* (Workspace: need 2*M*M + 4*M, prefer 2*M*M + 3*M + M*NB)
*
CALL DORGBR( 'Q', M, M, M, WORK( IR ), LDWRKR,
$ WORK( ITAUQ ), WORK( IWORK ),
@@ -2546,7 +2642,7 @@
* Perform bidiagonal QR iteration, computing left
* singular vectors of L in WORK(IR) and computing
* right singular vectors of L in WORK(IU)
-* (Workspace: need 2*M*M+BDSPAC)
+* (Workspace: need 2*M*M + BDSPAC)
*
CALL DBDSQR( 'U', M, M, M, 0, S, WORK( IE ),
$ WORK( IU ), LDWRKU, WORK( IR ),
@@ -2573,14 +2669,14 @@
IWORK = ITAU + M
*
* Compute A=L*Q, copying result to VT
-* (Workspace: need 2*M, prefer M+M*NB)
+* (Workspace: need 2*M, prefer M + M*NB)
*
CALL DGELQF( M, N, A, LDA, WORK( ITAU ),
$ WORK( IWORK ), LWORK-IWORK+1, IERR )
CALL DLACPY( 'U', M, N, A, LDA, VT, LDVT )
*
* Generate Q in VT
-* (Workspace: need 2*M, prefer M+M*NB)
+* (Workspace: need 2*M, prefer M + M*NB)
*
CALL DORGLQ( M, N, M, VT, LDVT, WORK( ITAU ),
$ WORK( IWORK ), LWORK-IWORK+1, IERR )
@@ -2595,21 +2691,21 @@
$ LDA )
*
* Bidiagonalize L in A
-* (Workspace: need 4*M, prefer 3*M+2*M*NB)
+* (Workspace: need 4*M, prefer 3*M + 2*M*NB)
*
CALL DGEBRD( M, M, A, LDA, S, WORK( IE ),
$ WORK( ITAUQ ), WORK( ITAUP ),
$ WORK( IWORK ), LWORK-IWORK+1, IERR )
*
* Multiply right vectors bidiagonalizing L by Q in VT
-* (Workspace: need 3*M+N, prefer 3*M+N*NB)
+* (Workspace: need 3*M + N, prefer 3*M + N*NB)
*
CALL DORMBR( 'P', 'L', 'T', M, N, M, A, LDA,
$ WORK( ITAUP ), VT, LDVT,
$ WORK( IWORK ), LWORK-IWORK+1, IERR )
*
* Generate left bidiagonalizing vectors of L in A
-* (Workspace: need 4*M, prefer 3*M+M*NB)
+* (Workspace: need 4*M, prefer 3*M + M*NB)
*
CALL DORGBR( 'Q', M, M, M, A, LDA, WORK( ITAUQ ),
$ WORK( IWORK ), LWORK-IWORK+1, IERR )
@@ -2653,7 +2749,7 @@
IWORK = ITAU + M
*
* Compute A=L*Q
-* (Workspace: need M*M+2*M, prefer M*M+M+M*NB)
+* (Workspace: need M*M + 2*M, prefer M*M + M + M*NB)
*
CALL DGELQF( M, N, A, LDA, WORK( ITAU ),
$ WORK( IWORK ), LWORK-IWORK+1, IERR )
@@ -2666,7 +2762,7 @@
$ WORK( IU+LDWRKU ), LDWRKU )
*
* Generate Q in A
-* (Workspace: need M*M+2*M, prefer M*M+M+M*NB)
+* (Workspace: need M*M + 2*M, prefer M*M + M + M*NB)
*
CALL DORGLQ( M, N, M, A, LDA, WORK( ITAU ),
$ WORK( IWORK ), LWORK-IWORK+1, IERR )
@@ -2676,7 +2772,7 @@
IWORK = ITAUP + M
*
* Bidiagonalize L in WORK(IU), copying result to U
-* (Workspace: need M*M+4*M, prefer M*M+3*M+2*M*NB)
+* (Workspace: need M*M + 4*M, prefer M*M + 3*M + 2*M*NB)
*
CALL DGEBRD( M, M, WORK( IU ), LDWRKU, S,
$ WORK( IE ), WORK( ITAUQ ),
@@ -2686,7 +2782,7 @@
$ LDU )
*
* Generate right bidiagonalizing vectors in WORK(IU)
-* (Workspace: need M*M+4*M-1,
+* (Workspace: need M*M + 4*M-1,
* prefer M*M+3*M+(M-1)*NB)
*
CALL DORGBR( 'P', M, M, M, WORK( IU ), LDWRKU,
@@ -2694,7 +2790,7 @@
$ LWORK-IWORK+1, IERR )
*
* Generate left bidiagonalizing vectors in U
-* (Workspace: need M*M+4*M, prefer M*M+3*M+M*NB)
+* (Workspace: need M*M + 4*M, prefer M*M + 3*M + M*NB)
*
CALL DORGBR( 'Q', M, M, M, U, LDU, WORK( ITAUQ ),
$ WORK( IWORK ), LWORK-IWORK+1, IERR )
@@ -2703,7 +2799,7 @@
* Perform bidiagonal QR iteration, computing left
* singular vectors of L in U and computing right
* 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 ),
$ WORK( IU ), LDWRKU, U, LDU, DUM, 1,
@@ -2724,14 +2820,14 @@
IWORK = ITAU + M
*
* Compute A=L*Q, copying result to VT
-* (Workspace: need 2*M, prefer M+M*NB)
+* (Workspace: need 2*M, prefer M + M*NB)
*
CALL DGELQF( M, N, A, LDA, WORK( ITAU ),
$ WORK( IWORK ), LWORK-IWORK+1, IERR )
CALL DLACPY( 'U', M, N, A, LDA, VT, LDVT )
*
* Generate Q in VT
-* (Workspace: need 2*M, prefer M+M*NB)
+* (Workspace: need 2*M, prefer M + M*NB)
*
CALL DORGLQ( M, N, M, VT, LDVT, WORK( ITAU ),
$ WORK( IWORK ), LWORK-IWORK+1, IERR )
@@ -2747,7 +2843,7 @@
IWORK = ITAUP + M
*
* Bidiagonalize L in U
-* (Workspace: need 4*M, prefer 3*M+2*M*NB)
+* (Workspace: need 4*M, prefer 3*M + 2*M*NB)
*
CALL DGEBRD( M, M, U, LDU, S, WORK( IE ),
$ WORK( ITAUQ ), WORK( ITAUP ),
@@ -2755,14 +2851,14 @@
*
* Multiply right bidiagonalizing vectors in U by Q
* in VT
-* (Workspace: need 3*M+N, prefer 3*M+N*NB)
+* (Workspace: need 3*M + N, prefer 3*M + N*NB)
*
CALL DORMBR( 'P', 'L', 'T', M, N, M, U, LDU,
$ WORK( ITAUP ), VT, LDVT,
$ WORK( IWORK ), LWORK-IWORK+1, IERR )
*
* Generate left bidiagonalizing vectors in U
-* (Workspace: need 4*M, prefer 3*M+M*NB)
+* (Workspace: need 4*M, prefer 3*M + M*NB)
*
CALL DORGBR( 'Q', M, M, M, U, LDU, WORK( ITAUQ ),
$ WORK( IWORK ), LWORK-IWORK+1, IERR )
@@ -2789,7 +2885,7 @@
* N right singular vectors to be computed in VT and
* no left singular vectors to be computed
*
- IF( LWORK.GE.M*M+MAX( N+M, 4*M, BDSPAC ) ) THEN
+ IF( LWORK.GE.M*M+MAX( N + M, 4*M, BDSPAC ) ) THEN
*
* Sufficient workspace for a fast algorithm
*
@@ -2809,7 +2905,7 @@
IWORK = ITAU + M
*
* Compute A=L*Q, copying result to VT
-* (Workspace: need M*M+2*M, prefer M*M+M+M*NB)
+* (Workspace: need M*M + 2*M, prefer M*M + M + M*NB)
*
CALL DGELQF( M, N, A, LDA, WORK( ITAU ),
$ WORK( IWORK ), LWORK-IWORK+1, IERR )
@@ -2823,7 +2919,7 @@
$ WORK( IR+LDWRKR ), LDWRKR )
*
* Generate Q in VT
-* (Workspace: need M*M+M+N, prefer M*M+M+N*NB)
+* (Workspace: need M*M + M + N, prefer M*M + M + N*NB)
*
CALL DORGLQ( N, N, M, VT, LDVT, WORK( ITAU ),
$ WORK( IWORK ), LWORK-IWORK+1, IERR )
@@ -2833,7 +2929,7 @@
IWORK = ITAUP + M
*
* Bidiagonalize L in WORK(IR)
-* (Workspace: need M*M+4*M, prefer M*M+3*M+2*M*NB)
+* (Workspace: need M*M + 4*M, prefer M*M + 3*M + 2*M*NB)
*
CALL DGEBRD( M, M, WORK( IR ), LDWRKR, S,
$ WORK( IE ), WORK( ITAUQ ),
@@ -2841,7 +2937,7 @@
$ LWORK-IWORK+1, IERR )
*
* Generate right bidiagonalizing vectors in WORK(IR)
-* (Workspace: need M*M+4*M-1,
+* (Workspace: need M*M + 4*M-1,
* prefer M*M+3*M+(M-1)*NB)
*
CALL DORGBR( 'P', M, M, M, WORK( IR ), LDWRKR,
@@ -2851,7 +2947,7 @@
*
* Perform bidiagonal QR iteration, computing right
* 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 ),
$ WORK( IR ), LDWRKR, DUM, 1, DUM, 1,
@@ -2876,14 +2972,14 @@
IWORK = ITAU + M
*
* Compute A=L*Q, copying result to VT
-* (Workspace: need 2*M, prefer M+M*NB)
+* (Workspace: need 2*M, prefer M + M*NB)
*
CALL DGELQF( M, N, A, LDA, WORK( ITAU ),
$ WORK( IWORK ), LWORK-IWORK+1, IERR )
CALL DLACPY( 'U', M, N, A, LDA, VT, LDVT )
*
* 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 ),
$ WORK( IWORK ), LWORK-IWORK+1, IERR )
@@ -2898,7 +2994,7 @@
$ LDA )
*
* Bidiagonalize L in A
-* (Workspace: need 4*M, prefer 3*M+2*M*NB)
+* (Workspace: need 4*M, prefer 3*M + 2*M*NB)
*
CALL DGEBRD( M, M, A, LDA, S, WORK( IE ),
$ WORK( ITAUQ ), WORK( ITAUP ),
@@ -2906,7 +3002,7 @@
*
* Multiply right bidiagonalizing vectors in A by Q
* in VT
-* (Workspace: need 3*M+N, prefer 3*M+N*NB)
+* (Workspace: need 3*M + N, prefer 3*M + N*NB)
*
CALL DORMBR( 'P', 'L', 'T', M, N, M, A, LDA,
$ WORK( ITAUP ), VT, LDVT,
@@ -2929,7 +3025,7 @@
* N right singular vectors to be computed in VT and
* M left singular vectors to be overwritten on A
*
- IF( LWORK.GE.2*M*M+MAX( N+M, 4*M, BDSPAC ) ) THEN
+ IF( LWORK.GE.2*M*M+MAX( N + M, 4*M, BDSPAC ) ) THEN
*
* Sufficient workspace for a fast algorithm
*
@@ -2941,7 +3037,7 @@
LDWRKU = LDA
IR = IU + LDWRKU*M
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
*
@@ -2960,14 +3056,14 @@
IWORK = ITAU + M
*
* Compute A=L*Q, copying result to VT
-* (Workspace: need 2*M*M+2*M, prefer 2*M*M+M+M*NB)
+* (Workspace: need 2*M*M + 2*M, prefer 2*M*M + M + M*NB)
*
CALL DGELQF( M, N, A, LDA, WORK( ITAU ),
$ WORK( IWORK ), LWORK-IWORK+1, IERR )
CALL DLACPY( 'U', M, N, A, LDA, VT, LDVT )
*
* Generate Q in VT
-* (Workspace: need 2*M*M+M+N, prefer 2*M*M+M+N*NB)
+* (Workspace: need 2*M*M + M + N, prefer 2*M*M + M + N*NB)
*
CALL DORGLQ( N, N, M, VT, LDVT, WORK( ITAU ),
$ WORK( IWORK ), LWORK-IWORK+1, IERR )
@@ -2985,7 +3081,7 @@
*
* Bidiagonalize L in WORK(IU), copying result to
* WORK(IR)
-* (Workspace: need 2*M*M+4*M,
+* (Workspace: need 2*M*M + 4*M,
* prefer 2*M*M+3*M+2*M*NB)
*
CALL DGEBRD( M, M, WORK( IU ), LDWRKU, S,
@@ -2996,7 +3092,7 @@
$ WORK( IR ), LDWRKR )
*
* Generate right bidiagonalizing vectors in WORK(IU)
-* (Workspace: need 2*M*M+4*M-1,
+* (Workspace: need 2*M*M + 4*M-1,
* prefer 2*M*M+3*M+(M-1)*NB)
*
CALL DORGBR( 'P', M, M, M, WORK( IU ), LDWRKU,
@@ -3004,7 +3100,7 @@
$ LWORK-IWORK+1, IERR )
*
* Generate left bidiagonalizing vectors in WORK(IR)
-* (Workspace: need 2*M*M+4*M, prefer 2*M*M+3*M+M*NB)
+* (Workspace: need 2*M*M + 4*M, prefer 2*M*M + 3*M + M*NB)
*
CALL DORGBR( 'Q', M, M, M, WORK( IR ), LDWRKR,
$ WORK( ITAUQ ), WORK( IWORK ),
@@ -3014,7 +3110,7 @@
* Perform bidiagonal QR iteration, computing left
* singular vectors of L in WORK(IR) and computing
* right singular vectors of L in WORK(IU)
-* (Workspace: need 2*M*M+BDSPAC)
+* (Workspace: need 2*M*M + BDSPAC)
*
CALL DBDSQR( 'U', M, M, M, 0, S, WORK( IE ),
$ WORK( IU ), LDWRKU, WORK( IR ),
@@ -3044,14 +3140,14 @@
IWORK = ITAU + M
*
* Compute A=L*Q, copying result to VT
-* (Workspace: need 2*M, prefer M+M*NB)
+* (Workspace: need 2*M, prefer M + M*NB)
*
CALL DGELQF( M, N, A, LDA, WORK( ITAU ),
$ WORK( IWORK ), LWORK-IWORK+1, IERR )
CALL DLACPY( 'U', M, N, A, LDA, VT, LDVT )
*
* 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 ),
$ WORK( IWORK ), LWORK-IWORK+1, IERR )
@@ -3066,7 +3162,7 @@
$ LDA )
*
* Bidiagonalize L in A
-* (Workspace: need 4*M, prefer 3*M+2*M*NB)
+* (Workspace: need 4*M, prefer 3*M + 2*M*NB)
*
CALL DGEBRD( M, M, A, LDA, S, WORK( IE ),
$ WORK( ITAUQ ), WORK( ITAUP ),
@@ -3074,14 +3170,14 @@
*
* Multiply right bidiagonalizing vectors in A by Q
* in VT
-* (Workspace: need 3*M+N, prefer 3*M+N*NB)
+* (Workspace: need 3*M + N, prefer 3*M + N*NB)
*
CALL DORMBR( 'P', 'L', 'T', M, N, M, A, LDA,
$ WORK( ITAUP ), VT, LDVT,
$ WORK( IWORK ), LWORK-IWORK+1, IERR )
*
* Generate left bidiagonalizing vectors in A
-* (Workspace: need 4*M, prefer 3*M+M*NB)
+* (Workspace: need 4*M, prefer 3*M + M*NB)
*
CALL DORGBR( 'Q', M, M, M, A, LDA, WORK( ITAUQ ),
$ WORK( IWORK ), LWORK-IWORK+1, IERR )
@@ -3105,7 +3201,7 @@
* N right singular vectors to be computed in VT and
* M left singular vectors to be computed in U
*
- IF( LWORK.GE.M*M+MAX( N+M, 4*M, BDSPAC ) ) THEN
+ IF( LWORK.GE.M*M+MAX( N + M, 4*M, BDSPAC ) ) THEN
*
* Sufficient workspace for a fast algorithm
*
@@ -3125,14 +3221,14 @@
IWORK = ITAU + M
*
* Compute A=L*Q, copying result to VT
-* (Workspace: need M*M+2*M, prefer M*M+M+M*NB)
+* (Workspace: need M*M + 2*M, prefer M*M + M + M*NB)
*
CALL DGELQF( M, N, A, LDA, WORK( ITAU ),
$ WORK( IWORK ), LWORK-IWORK+1, IERR )
CALL DLACPY( 'U', M, N, A, LDA, VT, LDVT )
*
* Generate Q in VT
-* (Workspace: need M*M+M+N, prefer M*M+M+N*NB)
+* (Workspace: need M*M + M + N, prefer M*M + M + N*NB)
*
CALL DORGLQ( N, N, M, VT, LDVT, WORK( ITAU ),
$ WORK( IWORK ), LWORK-IWORK+1, IERR )
@@ -3149,7 +3245,7 @@
IWORK = ITAUP + M
*
* Bidiagonalize L in WORK(IU), copying result to U
-* (Workspace: need M*M+4*M, prefer M*M+3*M+2*M*NB)
+* (Workspace: need M*M + 4*M, prefer M*M + 3*M + 2*M*NB)
*
CALL DGEBRD( M, M, WORK( IU ), LDWRKU, S,
$ WORK( IE ), WORK( ITAUQ ),
@@ -3159,14 +3255,14 @@
$ LDU )
*
* Generate right bidiagonalizing vectors in WORK(IU)
-* (Workspace: need M*M+4*M, prefer M*M+3*M+(M-1)*NB)
+* (Workspace: need M*M + 4*M, prefer M*M + 3*M + (M-1)*NB)
*
CALL DORGBR( 'P', M, M, M, WORK( IU ), LDWRKU,
$ WORK( ITAUP ), WORK( IWORK ),
$ LWORK-IWORK+1, IERR )
*
* Generate left bidiagonalizing vectors in U
-* (Workspace: need M*M+4*M, prefer M*M+3*M+M*NB)
+* (Workspace: need M*M + 4*M, prefer M*M + 3*M + M*NB)
*
CALL DORGBR( 'Q', M, M, M, U, LDU, WORK( ITAUQ ),
$ WORK( IWORK ), LWORK-IWORK+1, IERR )
@@ -3175,7 +3271,7 @@
* Perform bidiagonal QR iteration, computing left
* singular vectors of L in U and computing right
* 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 ),
$ WORK( IU ), LDWRKU, U, LDU, DUM, 1,
@@ -3200,14 +3296,14 @@
IWORK = ITAU + M
*
* Compute A=L*Q, copying result to VT
-* (Workspace: need 2*M, prefer M+M*NB)
+* (Workspace: need 2*M, prefer M + M*NB)
*
CALL DGELQF( M, N, A, LDA, WORK( ITAU ),
$ WORK( IWORK ), LWORK-IWORK+1, IERR )
CALL DLACPY( 'U', M, N, A, LDA, VT, LDVT )
*
* 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 ),
$ WORK( IWORK ), LWORK-IWORK+1, IERR )
@@ -3223,7 +3319,7 @@
IWORK = ITAUP + M
*
* Bidiagonalize L in U
-* (Workspace: need 4*M, prefer 3*M+2*M*NB)
+* (Workspace: need 4*M, prefer 3*M + 2*M*NB)
*
CALL DGEBRD( M, M, U, LDU, S, WORK( IE ),
$ WORK( ITAUQ ), WORK( ITAUP ),
@@ -3231,14 +3327,14 @@
*
* Multiply right bidiagonalizing vectors in U by Q
* in VT
-* (Workspace: need 3*M+N, prefer 3*M+N*NB)
+* (Workspace: need 3*M + N, prefer 3*M + N*NB)
*
CALL DORMBR( 'P', 'L', 'T', M, N, M, U, LDU,
$ WORK( ITAUP ), VT, LDVT,
$ WORK( IWORK ), LWORK-IWORK+1, IERR )
*
* Generate left bidiagonalizing vectors in U
-* (Workspace: need 4*M, prefer 3*M+M*NB)
+* (Workspace: need 4*M, prefer 3*M + M*NB)
*
CALL DORGBR( 'Q', M, M, M, U, LDU, WORK( ITAUQ ),
$ WORK( IWORK ), LWORK-IWORK+1, IERR )
@@ -3272,7 +3368,7 @@
IWORK = ITAUP + M
*
* Bidiagonalize A
-* (Workspace: need 3*M+N, prefer 3*M+(M+N)*NB)
+* (Workspace: need 3*M + N, prefer 3*M + (M + N)*NB)
*
CALL DGEBRD( M, N, A, LDA, S, WORK( IE ), WORK( ITAUQ ),
$ WORK( ITAUP ), WORK( IWORK ), LWORK-IWORK+1,
@@ -3281,7 +3377,7 @@
*
* If left singular vectors desired in U, copy result to U
* and generate left bidiagonalizing vectors in U
-* (Workspace: need 4*M-1, prefer 3*M+(M-1)*NB)
+* (Workspace: need 4*M-1, prefer 3*M + (M-1)*NB)
*
CALL DLACPY( 'L', M, M, A, LDA, U, LDU )
CALL DORGBR( 'Q', M, M, N, U, LDU, WORK( ITAUQ ),
@@ -3291,7 +3387,7 @@
*
* If right singular vectors desired in VT, copy result to
* VT and generate right bidiagonalizing vectors in VT
-* (Workspace: need 3*M+NRVT, prefer 3*M+NRVT*NB)
+* (Workspace: need 3*M + NRVT, prefer 3*M + NRVT*NB)
*
CALL DLACPY( 'U', M, N, A, LDA, VT, LDVT )
IF( WNTVA )
@@ -3305,7 +3401,7 @@
*
* If left singular vectors desired in A, generate left
* bidiagonalizing vectors in A
-* (Workspace: need 4*M-1, prefer 3*M+(M-1)*NB)
+* (Workspace: need 4*M-1, prefer 3*M + (M-1)*NB)
*
CALL DORGBR( 'Q', M, M, N, A, LDA, WORK( ITAUQ ),
$ WORK( IWORK ), LWORK-IWORK+1, IERR )
@@ -3314,7 +3410,7 @@
*
* If right singular vectors desired in A, generate right
* bidiagonalizing vectors in A
-* (Workspace: need 4*M, prefer 3*M+M*NB)
+* (Workspace: need 4*M, prefer 3*M + M*NB)
*
CALL DORGBR( 'P', M, N, M, A, LDA, WORK( ITAUP ),
$ WORK( IWORK ), LWORK-IWORK+1, IERR )