--- rpl/lapack/lapack/dgesvdx.f 2015/11/26 11:44:15 1.1
+++ rpl/lapack/lapack/dgesvdx.f 2023/08/07 08:38:50 1.8
@@ -2,26 +2,26 @@
*
* =========== DOCUMENTATION ===========
*
-* Online html documentation available at
-* http://www.netlib.org/lapack/explore-html/
+* Online html documentation available at
+* http://www.netlib.org/lapack/explore-html/
*
*> \htmlonly
-*> Download DGESVDX + dependencies
-*>
-*> [TGZ]
-*>
-*> [ZIP]
-*>
+*> Download DGESVDX + dependencies
+*>
+*> [TGZ]
+*>
+*> [ZIP]
+*>
*> [TXT]
-*> \endhtmlonly
+*> \endhtmlonly
*
* Definition:
* ===========
*
-* SUBROUTINE DGESVDX( JOBU, JOBVT, RANGE, M, N, A, LDA, VL, VU,
-* $ IL, IU, NS, S, U, LDU, VT, LDVT, WORK,
+* SUBROUTINE DGESVDX( JOBU, JOBVT, RANGE, M, N, A, LDA, VL, VU,
+* $ IL, IU, NS, S, U, LDU, VT, LDVT, WORK,
* $ LWORK, IWORK, INFO )
-*
+*
*
* .. Scalar Arguments ..
* CHARACTER JOBU, JOBVT, RANGE
@@ -33,7 +33,7 @@
* DOUBLE PRECISION A( LDA, * ), S( * ), U( LDU, * ),
* $ VT( LDVT, * ), WORK( * )
* ..
-*
+*
*
*> \par Purpose:
* =============
@@ -43,23 +43,23 @@
*> DGESVDX 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.
-*>
-*> DGESVDX uses an eigenvalue problem for obtaining the SVD, which
-*> allows for the computation of a subset of singular values and
+*>
+*> DGESVDX uses an eigenvalue problem for obtaining the SVD, which
+*> allows for the computation of a subset of singular values and
*> vectors. See DBDSVDX for details.
-*>
+*>
*> Note that the routine returns V**T, not V.
*> \endverbatim
-*
+*
* Arguments:
* ==========
*
@@ -68,7 +68,7 @@
*> JOBU is CHARACTER*1
*> Specifies options for computing all or part of the matrix U:
*> = 'V': the first min(m,n) columns of U (the left singular
-*> vectors) or as specified by RANGE are returned in
+*> vectors) or as specified by RANGE are returned in
*> the array U;
*> = 'N': no columns of U (no left singular vectors) are
*> computed.
@@ -80,7 +80,7 @@
*> Specifies options for computing all or part of the matrix
*> V**T:
*> = 'V': the first min(m,n) rows of V**T (the right singular
-*> vectors) or as specified by RANGE are returned in
+*> vectors) or as specified by RANGE are returned in
*> the array VT;
*> = 'N': no rows of V**T (no right singular vectors) are
*> computed.
@@ -92,7 +92,7 @@
*> = 'A': all singular values will be found.
*> = 'V': all singular values in the half-open interval (VL,VU]
*> will be found.
-*> = 'I': the IL-th through IU-th singular values will be found.
+*> = 'I': the IL-th through IU-th singular values will be found.
*> \endverbatim
*>
*> \param[in] M
@@ -123,13 +123,15 @@
*> \param[in] VL
*> \verbatim
*> VL is DOUBLE PRECISION
-*> VL >=0.
+*> If RANGE='V', the lower bound of the interval to
+*> be searched for singular values. VU > VL.
+*> Not referenced if RANGE = 'A' or 'I'.
*> \endverbatim
*>
*> \param[in] VU
*> \verbatim
*> VU is DOUBLE PRECISION
-*> If RANGE='V', the lower and upper bounds of the interval to
+*> If RANGE='V', the upper bound of the interval to
*> be searched for singular values. VU > VL.
*> Not referenced if RANGE = 'A' or 'I'.
*> \endverbatim
@@ -137,13 +139,17 @@
*> \param[in] IL
*> \verbatim
*> IL is INTEGER
+*> If RANGE='I', the index of the
+*> smallest singular value to be returned.
+*> 1 <= IL <= IU <= min(M,N), if min(M,N) > 0.
+*> Not referenced if RANGE = 'A' or 'V'.
*> \endverbatim
*>
*> \param[in] IU
*> \verbatim
*> IU is INTEGER
-*> If RANGE='I', the indices (in ascending order) of the
-*> smallest and largest singular values to be returned.
+*> If RANGE='I', the index of the
+*> largest singular value to be returned.
*> 1 <= IL <= IU <= min(M,N), if min(M,N) > 0.
*> Not referenced if RANGE = 'A' or 'V'.
*> \endverbatim
@@ -151,7 +157,7 @@
*> \param[out] NS
*> \verbatim
*> NS is INTEGER
-*> The total number of singular values found,
+*> The total number of singular values found,
*> 0 <= NS <= min(M,N).
*> If RANGE = 'A', NS = min(M,N); if RANGE = 'I', NS = IU-IL+1.
*> \endverbatim
@@ -165,11 +171,11 @@
*> \param[out] U
*> \verbatim
*> U is DOUBLE PRECISION array, dimension (LDU,UCOL)
-*> If JOBU = 'V', U contains columns of U (the left singular
-*> vectors, stored columnwise) as specified by RANGE; if
+*> If JOBU = 'V', U contains columns of U (the left singular
+*> vectors, stored columnwise) as specified by RANGE; if
*> JOBU = 'N', U is not referenced.
-*> Note: The user must ensure that UCOL >= NS; if RANGE = 'V',
-*> the exact value of NS is not known ILQFin advance and an upper
+*> Note: The user must ensure that UCOL >= NS; if RANGE = 'V',
+*> the exact value of NS is not known in advance and an upper
*> bound must be used.
*> \endverbatim
*>
@@ -183,11 +189,11 @@
*> \param[out] VT
*> \verbatim
*> VT is DOUBLE PRECISION array, dimension (LDVT,N)
-*> If JOBVT = 'V', VT contains the rows of V**T (the right singular
-*> vectors, stored rowwise) as specified by RANGE; if JOBVT = 'N',
+*> If JOBVT = 'V', VT contains the rows of V**T (the right singular
+*> vectors, stored rowwise) as specified by RANGE; if JOBVT = 'N',
*> VT is not referenced.
-*> Note: The user must ensure that LDVT >= NS; if RANGE = 'V',
-*> the exact value of NS is not known in advance and an upper
+*> Note: The user must ensure that LDVT >= NS; if RANGE = 'V',
+*> the exact value of NS is not known in advance and an upper
*> bound must be used.
*> \endverbatim
*>
@@ -208,9 +214,9 @@
*> \verbatim
*> LWORK is INTEGER
*> The dimension of the array WORK.
-*> LWORK >= MAX(1,MIN(M,N)*(MIN(M,N)+4)) for the paths (see
+*> LWORK >= MAX(1,MIN(M,N)*(MIN(M,N)+4)) for the paths (see
*> comments inside the code):
-*> - PATH 1 (M much larger than N)
+*> - PATH 1 (M much larger than N)
*> - PATH 1t (N much larger than M)
*> LWORK >= MAX(1,MIN(M,N)*2+MAX(M,N)) for the other paths.
*> For good performance, LWORK should generally be larger.
@@ -224,8 +230,8 @@
*> \param[out] IWORK
*> \verbatim
*> IWORK is INTEGER array, dimension (12*MIN(M,N))
-*> If INFO = 0, the first NS elements of IWORK are zero. If INFO > 0,
-*> then IWORK contains the indices of the eigenvectors that failed
+*> If INFO = 0, the first NS elements of IWORK are zero. If INFO > 0,
+*> then IWORK contains the indices of the eigenvectors that failed
*> to converge in DBDSVDX/DSTEVX.
*> \endverbatim
*>
@@ -243,24 +249,21 @@
* Authors:
* ========
*
-*> \author Univ. of Tennessee
-*> \author Univ. of California Berkeley
-*> \author Univ. of Colorado Denver
-*> \author NAG Ltd.
-*
-*> \date November 2015
+*> \author Univ. of Tennessee
+*> \author Univ. of California Berkeley
+*> \author Univ. of Colorado Denver
+*> \author NAG Ltd.
*
*> \ingroup doubleGEsing
*
* =====================================================================
- SUBROUTINE DGESVDX( JOBU, JOBVT, RANGE, M, N, A, LDA, VL, VU,
- $ IL, IU, NS, S, U, LDU, VT, LDVT, WORK,
+ SUBROUTINE DGESVDX( JOBU, JOBVT, RANGE, M, N, A, LDA, VL, VU,
+ $ IL, IU, NS, S, U, LDU, VT, LDVT, WORK,
$ LWORK, IWORK, INFO )
*
-* -- LAPACK driver routine (version 3.6.0) --
+* -- LAPACK driver routine --
* -- LAPACK is a software package provided by Univ. of Tennessee, --
* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
-* November 2015
*
* .. Scalar Arguments ..
CHARACTER JOBU, JOBVT, RANGE
@@ -283,7 +286,7 @@
CHARACTER JOBZ, RNGTGK
LOGICAL ALLS, INDS, LQUERY, VALS, WANTU, WANTVT
INTEGER I, ID, IE, IERR, ILQF, ILTGK, IQRF, ISCL,
- $ ITAU, ITAUP, ITAUQ, ITEMP, ITGKZ, IUTGK,
+ $ ITAU, ITAUP, ITAUQ, ITEMP, ITGKZ, IUTGK,
$ J, MAXWRK, MINMN, MINWRK, MNTHR
DOUBLE PRECISION ABSTOL, ANRM, BIGNUM, EPS, SMLNUM
* ..
@@ -293,13 +296,13 @@
* .. External Subroutines ..
EXTERNAL DBDSVDX, DGEBRD, DGELQF, DGEQRF, DLACPY,
$ DLASCL, DLASET, DORMBR, DORMLQ, DORMQR,
- $ DSCAL, XERBLA
+ $ DCOPY, XERBLA
* ..
* .. External Functions ..
LOGICAL LSAME
INTEGER ILAENV
- DOUBLE PRECISION DLAMCH, DLANGE, DNRM2
- EXTERNAL LSAME, ILAENV, DLAMCH, DLANGE, DNRM2
+ DOUBLE PRECISION DLAMCH, DLANGE
+ EXTERNAL LSAME, ILAENV, DLAMCH, DLANGE
* ..
* .. Intrinsic Functions ..
INTRINSIC MAX, MIN, SQRT
@@ -357,8 +360,14 @@
IF( INFO.EQ.0 ) THEN
IF( WANTU .AND. LDU.LT.M ) THEN
INFO = -15
- ELSE IF( WANTVT .AND. LDVT.LT.MINMN ) THEN
- INFO = -16
+ ELSE IF( WANTVT ) THEN
+ IF( INDS ) THEN
+ IF( LDVT.LT.IU-IL+1 ) THEN
+ INFO = -17
+ END IF
+ ELSE IF( LDVT.LT.MINMN ) THEN
+ INFO = -17
+ END IF
END IF
END IF
END IF
@@ -380,18 +389,34 @@
*
* Path 1 (M much larger than N)
*
- MAXWRK = N*(N*2+16) +
+ MAXWRK = N +
$ N*ILAENV( 1, 'DGEQRF', ' ', M, N, -1, -1 )
- MAXWRK = MAX( MAXWRK, N*(N*2+20) + 2*N*
+ MAXWRK = MAX( MAXWRK, N*(N+5) + 2*N*
$ ILAENV( 1, 'DGEBRD', ' ', N, N, -1, -1 ) )
- MINWRK = N*(N*2+21)
+ IF (WANTU) THEN
+ MAXWRK = MAX(MAXWRK,N*(N*3+6)+N*
+ $ ILAENV( 1, 'DORMQR', ' ', N, N, -1, -1 ) )
+ END IF
+ IF (WANTVT) THEN
+ MAXWRK = MAX(MAXWRK,N*(N*3+6)+N*
+ $ ILAENV( 1, 'DORMLQ', ' ', N, N, -1, -1 ) )
+ END IF
+ MINWRK = N*(N*3+20)
ELSE
*
* Path 2 (M at least N, but not much larger)
*
- MAXWRK = N*(N*2+19) + ( M+N )*
+ MAXWRK = 4*N + ( M+N )*
$ ILAENV( 1, 'DGEBRD', ' ', M, N, -1, -1 )
- MINWRK = N*(N*2+20) + M
+ IF (WANTU) THEN
+ MAXWRK = MAX(MAXWRK,N*(N*2+5)+N*
+ $ ILAENV( 1, 'DORMQR', ' ', N, N, -1, -1 ) )
+ END IF
+ IF (WANTVT) THEN
+ MAXWRK = MAX(MAXWRK,N*(N*2+5)+N*
+ $ ILAENV( 1, 'DORMLQ', ' ', N, N, -1, -1 ) )
+ END IF
+ MINWRK = MAX(N*(N*2+19),4*N+M)
END IF
ELSE
MNTHR = ILAENV( 6, 'DGESVD', JOBU // JOBVT, M, N, 0, 0 )
@@ -399,18 +424,34 @@
*
* Path 1t (N much larger than M)
*
- MAXWRK = M*(M*2+16) +
+ MAXWRK = M +
$ M*ILAENV( 1, 'DGELQF', ' ', M, N, -1, -1 )
- MAXWRK = MAX( MAXWRK, M*(M*2+20) + 2*M*
+ MAXWRK = MAX( MAXWRK, M*(M+5) + 2*M*
$ ILAENV( 1, 'DGEBRD', ' ', M, M, -1, -1 ) )
- MINWRK = M*(M*2+21)
+ IF (WANTU) THEN
+ MAXWRK = MAX(MAXWRK,M*(M*3+6)+M*
+ $ ILAENV( 1, 'DORMQR', ' ', M, M, -1, -1 ) )
+ END IF
+ IF (WANTVT) THEN
+ MAXWRK = MAX(MAXWRK,M*(M*3+6)+M*
+ $ ILAENV( 1, 'DORMLQ', ' ', M, M, -1, -1 ) )
+ END IF
+ MINWRK = M*(M*3+20)
ELSE
*
-* Path 2t (N greater than M, but not much larger)
+* Path 2t (N at least M, but not much larger)
*
- MAXWRK = M*(M*2+19) + ( M+N )*
+ MAXWRK = 4*M + ( M+N )*
$ ILAENV( 1, 'DGEBRD', ' ', M, N, -1, -1 )
- MINWRK = M*(M*2+20) + N
+ IF (WANTU) THEN
+ MAXWRK = MAX(MAXWRK,M*(M*2+5)+M*
+ $ ILAENV( 1, 'DORMQR', ' ', M, M, -1, -1 ) )
+ END IF
+ IF (WANTVT) THEN
+ MAXWRK = MAX(MAXWRK,M*(M*2+5)+M*
+ $ ILAENV( 1, 'DORMLQ', ' ', M, M, -1, -1 ) )
+ END IF
+ MINWRK = MAX(M*(M*2+19),4*M+N)
END IF
END IF
END IF
@@ -445,7 +486,7 @@
RNGTGK = 'I'
ILTGK = IL
IUTGK = IU
- ELSE
+ ELSE
RNGTGK = 'V'
ILTGK = 0
IUTGK = 0
@@ -489,7 +530,7 @@
ITEMP = ITAU + N
CALL DGEQRF( M, N, A, LDA, WORK( ITAU ), WORK( ITEMP ),
$ LWORK-ITEMP+1, INFO )
-*
+*
* Copy R into WORK and bidiagonalize it:
* (Workspace: need N*N+5*N, prefer N*N+4*N+2*N*NB)
*
@@ -498,19 +539,19 @@
IE = ID + N
ITAUQ = IE + N
ITAUP = ITAUQ + N
- ITEMP = ITAUP + N
+ ITEMP = ITAUP + N
CALL DLACPY( 'U', N, N, A, LDA, WORK( IQRF ), N )
CALL DLASET( 'L', N-1, N-1, ZERO, ZERO, WORK( IQRF+1 ), N )
- CALL DGEBRD( N, N, WORK( IQRF ), N, WORK( ID ), WORK( IE ),
+ CALL DGEBRD( N, N, WORK( IQRF ), N, WORK( ID ), WORK( IE ),
$ WORK( ITAUQ ), WORK( ITAUP ), WORK( ITEMP ),
$ LWORK-ITEMP+1, INFO )
*
* Solve eigenvalue problem TGK*Z=Z*S.
-* (Workspace: need 14*N + 2*N*(N+1))
-*
+* (Workspace: need 14*N + 2*N*(N+1))
+*
ITGKZ = ITEMP
ITEMP = ITGKZ + N*(N*2+1)
- CALL DBDSVDX( 'U', JOBZ, RNGTGK, N, WORK( ID ), WORK( IE ),
+ CALL DBDSVDX( 'U', JOBZ, RNGTGK, N, WORK( ID ), WORK( IE ),
$ VL, VU, ILTGK, IUTGK, NS, S, WORK( ITGKZ ),
$ N*2, WORK( ITEMP ), IWORK, INFO)
*
@@ -522,23 +563,23 @@
CALL DCOPY( N, WORK( J ), 1, U( 1,I ), 1 )
J = J + N*2
END DO
- CALL DLASET( 'A', M-N, N, ZERO, ZERO, U( N+1,1 ), LDU )
+ CALL DLASET( 'A', M-N, NS, ZERO, ZERO, U( N+1,1 ), LDU )
*
* Call DORMBR to compute QB*UB.
* (Workspace in WORK( ITEMP ): need N, prefer N*NB)
*
- CALL DORMBR( 'Q', 'L', 'N', N, NS, N, WORK( IQRF ), N,
- $ WORK( ITAUQ ), U, LDU, WORK( ITEMP ),
+ CALL DORMBR( 'Q', 'L', 'N', N, NS, N, WORK( IQRF ), N,
+ $ WORK( ITAUQ ), U, LDU, WORK( ITEMP ),
$ LWORK-ITEMP+1, INFO )
*
* Call DORMQR to compute Q*(QB*UB).
* (Workspace in WORK( ITEMP ): need N, prefer N*NB)
*
- CALL DORMQR( 'L', 'N', M, NS, N, A, LDA,
+ CALL DORMQR( 'L', 'N', M, NS, N, A, LDA,
$ WORK( ITAU ), U, LDU, WORK( ITEMP ),
$ LWORK-ITEMP+1, INFO )
- END IF
-*
+ END IF
+*
* If needed, compute right singular vectors.
*
IF( WANTVT) THEN
@@ -551,7 +592,7 @@
* Call DORMBR to compute VB**T * PB**T
* (Workspace in WORK( ITEMP ): need N, prefer N*NB)
*
- CALL DORMBR( 'P', 'R', 'T', NS, N, N, WORK( IQRF ), N,
+ CALL DORMBR( 'P', 'R', 'T', NS, N, N, WORK( IQRF ), N,
$ WORK( ITAUP ), VT, LDVT, WORK( ITEMP ),
$ LWORK-ITEMP+1, INFO )
END IF
@@ -569,17 +610,17 @@
IE = ID + N
ITAUQ = IE + N
ITAUP = ITAUQ + N
- ITEMP = ITAUP + N
- CALL DGEBRD( M, N, A, LDA, WORK( ID ), WORK( IE ),
+ ITEMP = ITAUP + N
+ CALL DGEBRD( M, N, A, LDA, WORK( ID ), WORK( IE ),
$ WORK( ITAUQ ), WORK( ITAUP ), WORK( ITEMP ),
$ LWORK-ITEMP+1, INFO )
*
* Solve eigenvalue problem TGK*Z=Z*S.
-* (Workspace: need 14*N + 2*N*(N+1))
-*
+* (Workspace: need 14*N + 2*N*(N+1))
+*
ITGKZ = ITEMP
ITEMP = ITGKZ + N*(N*2+1)
- CALL DBDSVDX( 'U', JOBZ, RNGTGK, N, WORK( ID ), WORK( IE ),
+ CALL DBDSVDX( 'U', JOBZ, RNGTGK, N, WORK( ID ), WORK( IE ),
$ VL, VU, ILTGK, IUTGK, NS, S, WORK( ITGKZ ),
$ N*2, WORK( ITEMP ), IWORK, INFO)
*
@@ -591,16 +632,16 @@
CALL DCOPY( N, WORK( J ), 1, U( 1,I ), 1 )
J = J + N*2
END DO
- CALL DLASET( 'A', M-N, N, ZERO, ZERO, U( N+1,1 ), LDU )
+ CALL DLASET( 'A', M-N, NS, ZERO, ZERO, U( N+1,1 ), LDU )
*
* Call DORMBR to compute QB*UB.
* (Workspace in WORK( ITEMP ): need N, prefer N*NB)
-*
- CALL DORMBR( 'Q', 'L', 'N', M, NS, N, A, LDA,
- $ WORK( ITAUQ ), U, LDU, WORK( ITEMP ),
+*
+ CALL DORMBR( 'Q', 'L', 'N', M, NS, N, A, LDA,
+ $ WORK( ITAUQ ), U, LDU, WORK( ITEMP ),
$ LWORK-ITEMP+1, IERR )
- END IF
-*
+ END IF
+*
* If needed, compute right singular vectors.
*
IF( WANTVT) THEN
@@ -613,11 +654,11 @@
* Call DORMBR to compute VB**T * PB**T
* (Workspace in WORK( ITEMP ): need N, prefer N*NB)
*
- CALL DORMBR( 'P', 'R', 'T', NS, N, N, A, LDA,
+ CALL DORMBR( 'P', 'R', 'T', NS, N, N, A, LDA,
$ WORK( ITAUP ), VT, LDVT, WORK( ITEMP ),
$ LWORK-ITEMP+1, IERR )
END IF
- END IF
+ END IF
ELSE
*
* A has more columns than rows. If A has sufficiently more
@@ -626,7 +667,7 @@
IF( N.GE.MNTHR ) THEN
*
* Path 1t (N much larger than M):
-* A = L * Q = ( QB * B * PB**T ) * Q
+* A = L * Q = ( QB * B * PB**T ) * Q
* = ( QB * ( UB * S * VB**T ) * PB**T ) * Q
* U = QB * UB ; V**T = VB**T * PB**T * Q
*
@@ -649,16 +690,16 @@
ITEMP = ITAUP + M
CALL DLACPY( 'L', M, M, A, LDA, WORK( ILQF ), M )
CALL DLASET( 'U', M-1, M-1, ZERO, ZERO, WORK( ILQF+M ), M )
- CALL DGEBRD( M, M, WORK( ILQF ), M, WORK( ID ), WORK( IE ),
+ CALL DGEBRD( M, M, WORK( ILQF ), M, WORK( ID ), WORK( IE ),
$ WORK( ITAUQ ), WORK( ITAUP ), WORK( ITEMP ),
$ LWORK-ITEMP+1, INFO )
*
* Solve eigenvalue problem TGK*Z=Z*S.
-* (Workspace: need 2*M*M+14*M)
+* (Workspace: need 2*M*M+14*M)
*
ITGKZ = ITEMP
ITEMP = ITGKZ + M*(M*2+1)
- CALL DBDSVDX( 'U', JOBZ, RNGTGK, M, WORK( ID ), WORK( IE ),
+ CALL DBDSVDX( 'U', JOBZ, RNGTGK, M, WORK( ID ), WORK( IE ),
$ VL, VU, ILTGK, IUTGK, NS, S, WORK( ITGKZ ),
$ M*2, WORK( ITEMP ), IWORK, INFO)
*
@@ -674,11 +715,11 @@
* Call DORMBR to compute QB*UB.
* (Workspace in WORK( ITEMP ): need M, prefer M*NB)
*
- CALL DORMBR( 'Q', 'L', 'N', M, NS, M, WORK( ILQF ), M,
- $ WORK( ITAUQ ), U, LDU, WORK( ITEMP ),
+ CALL DORMBR( 'Q', 'L', 'N', M, NS, M, WORK( ILQF ), M,
+ $ WORK( ITAUQ ), U, LDU, WORK( ITEMP ),
$ LWORK-ITEMP+1, INFO )
- END IF
-*
+ END IF
+*
* If needed, compute right singular vectors.
*
IF( WANTVT) THEN
@@ -687,28 +728,28 @@
CALL DCOPY( M, WORK( J ), 1, VT( I,1 ), LDVT )
J = J + M*2
END DO
- CALL DLASET( 'A', M, N-M, ZERO, ZERO, VT( 1,M+1 ), LDVT )
+ CALL DLASET( 'A', NS, N-M, ZERO, ZERO, VT( 1,M+1 ), LDVT)
*
* Call DORMBR to compute (VB**T)*(PB**T)
* (Workspace in WORK( ITEMP ): need M, prefer M*NB)
*
- CALL DORMBR( 'P', 'R', 'T', NS, M, M, WORK( ILQF ), M,
+ CALL DORMBR( 'P', 'R', 'T', NS, M, M, WORK( ILQF ), M,
$ WORK( ITAUP ), VT, LDVT, WORK( ITEMP ),
$ LWORK-ITEMP+1, INFO )
*
* Call DORMLQ to compute ((VB**T)*(PB**T))*Q.
* (Workspace in WORK( ITEMP ): need M, prefer M*NB)
*
- CALL DORMLQ( 'R', 'N', NS, N, M, A, LDA,
+ CALL DORMLQ( 'R', 'N', NS, N, M, A, LDA,
$ WORK( ITAU ), VT, LDVT, WORK( ITEMP ),
$ LWORK-ITEMP+1, INFO )
- END IF
+ END IF
ELSE
*
* Path 2t (N greater than M, but not much larger)
* Reduce to bidiagonal form without LQ decomposition
* A = QB * B * PB**T = QB * ( UB * S * VB**T ) * PB**T
-* U = QB * UB; V**T = VB**T * PB**T
+* U = QB * UB; V**T = VB**T * PB**T
*
* Bidiagonalize A
* (Workspace: need 4*M+N, prefer 4*M+(M+N)*NB)
@@ -718,19 +759,19 @@
ITAUQ = IE + M
ITAUP = ITAUQ + M
ITEMP = ITAUP + M
- CALL DGEBRD( M, N, A, LDA, WORK( ID ), WORK( IE ),
+ CALL DGEBRD( M, N, A, LDA, WORK( ID ), WORK( IE ),
$ WORK( ITAUQ ), WORK( ITAUP ), WORK( ITEMP ),
$ LWORK-ITEMP+1, INFO )
*
* Solve eigenvalue problem TGK*Z=Z*S.
-* (Workspace: need 2*M*M+14*M)
+* (Workspace: need 2*M*M+14*M)
*
ITGKZ = ITEMP
ITEMP = ITGKZ + M*(M*2+1)
- CALL DBDSVDX( 'L', JOBZ, RNGTGK, M, WORK( ID ), WORK( IE ),
+ CALL DBDSVDX( 'L', JOBZ, RNGTGK, M, WORK( ID ), WORK( IE ),
$ VL, VU, ILTGK, IUTGK, NS, S, WORK( ITGKZ ),
$ M*2, WORK( ITEMP ), IWORK, INFO)
-*
+*
* If needed, compute left singular vectors.
*
IF( WANTU ) THEN
@@ -743,11 +784,11 @@
* Call DORMBR to compute QB*UB.
* (Workspace in WORK( ITEMP ): need M, prefer M*NB)
*
- CALL DORMBR( 'Q', 'L', 'N', M, NS, N, A, LDA,
- $ WORK( ITAUQ ), U, LDU, WORK( ITEMP ),
+ CALL DORMBR( 'Q', 'L', 'N', M, NS, N, A, LDA,
+ $ WORK( ITAUQ ), U, LDU, WORK( ITEMP ),
$ LWORK-ITEMP+1, INFO )
- END IF
-*
+ END IF
+*
* If needed, compute right singular vectors.
*
IF( WANTVT) THEN
@@ -756,15 +797,15 @@
CALL DCOPY( M, WORK( J ), 1, VT( I,1 ), LDVT )
J = J + M*2
END DO
- CALL DLASET( 'A', M, N-M, ZERO, ZERO, VT( 1,M+1 ), LDVT )
+ CALL DLASET( 'A', NS, N-M, ZERO, ZERO, VT( 1,M+1 ), LDVT)
*
* Call DORMBR to compute VB**T * PB**T
* (Workspace in WORK( ITEMP ): need M, prefer M*NB)
*
- CALL DORMBR( 'P', 'R', 'T', NS, N, M, A, LDA,
+ CALL DORMBR( 'P', 'R', 'T', NS, N, M, A, LDA,
$ WORK( ITAUP ), VT, LDVT, WORK( ITEMP ),
$ LWORK-ITEMP+1, INFO )
- END IF
+ END IF
END IF
END IF
*