version 1.1, 2015/11/26 11:44:21
|
version 1.3, 2016/08/27 15:34:47
|
Line 36
|
Line 36
|
* .. |
* .. |
* |
* |
* |
* |
* Purpose |
*> \par Purpose: |
* ======= |
* ============= |
* |
*> |
* ZGESVDX computes the singular value decomposition (SVD) of a complex |
*> \verbatim |
* M-by-N matrix A, optionally computing the left and/or right singular |
*> |
* vectors. The SVD is written |
*> ZGESVDX computes the singular value decomposition (SVD) of a complex |
* |
*> M-by-N matrix A, optionally computing the left and/or right singular |
* A = U * SIGMA * transpose(V) |
*> vectors. The SVD is written |
* |
*> |
* where SIGMA is an M-by-N matrix which is zero except for its |
*> A = U * SIGMA * transpose(V) |
* min(m,n) diagonal elements, U is an M-by-M unitary matrix, and |
*> |
* V is an N-by-N unitary matrix. The diagonal elements of SIGMA |
*> where SIGMA is an M-by-N matrix which is zero except for its |
* are the singular values of A; they are real and non-negative, and |
*> min(m,n) diagonal elements, U is an M-by-M unitary matrix, and |
* are returned in descending order. The first min(m,n) columns of |
*> V is an N-by-N unitary matrix. The diagonal elements of SIGMA |
* U and V are the left and right singular vectors of A. |
*> 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 |
* ZGESVDX uses an eigenvalue problem for obtaining the SVD, which |
*> U and V are the left and right singular vectors of A. |
* allows for the computation of a subset of singular values and |
*> |
* vectors. See DBDSVDX for details. |
*> ZGESVDX uses an eigenvalue problem for obtaining the SVD, which |
* |
*> allows for the computation of a subset of singular values and |
* Note that the routine returns V**T, not V. |
*> vectors. See DBDSVDX for details. |
|
*> |
|
*> Note that the routine returns V**T, not V. |
|
*> \endverbatim |
* |
* |
* Arguments: |
* Arguments: |
* ========== |
* ========== |
Line 107
|
Line 110
|
*> |
*> |
*> \param[in,out] A |
*> \param[in,out] A |
*> \verbatim |
*> \verbatim |
*> A is COMPLEX array, dimension (LDA,N) |
*> A is COMPLEX*16 array, dimension (LDA,N) |
*> On entry, the M-by-N matrix A. |
*> On entry, the M-by-N matrix A. |
*> On exit, the contents of A are destroyed. |
*> On exit, the contents of A are destroyed. |
*> \endverbatim |
*> \endverbatim |
Line 121
|
Line 124
|
*> \param[in] VL |
*> \param[in] VL |
*> \verbatim |
*> \verbatim |
*> VL is DOUBLE PRECISION |
*> 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 |
*> \endverbatim |
*> |
*> |
*> \param[in] VU |
*> \param[in] VU |
*> \verbatim |
*> \verbatim |
*> VU is DOUBLE PRECISION |
*> 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. |
*> be searched for singular values. VU > VL. |
*> Not referenced if RANGE = 'A' or 'I'. |
*> Not referenced if RANGE = 'A' or 'I'. |
*> \endverbatim |
*> \endverbatim |
Line 135
|
Line 140
|
*> \param[in] IL |
*> \param[in] IL |
*> \verbatim |
*> \verbatim |
*> IL is INTEGER |
*> 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 |
*> \endverbatim |
*> |
*> |
*> \param[in] IU |
*> \param[in] IU |
*> \verbatim |
*> \verbatim |
*> IU is INTEGER |
*> IU is INTEGER |
*> If RANGE='I', the indices (in ascending order) of the |
*> If RANGE='I', the index of the |
*> smallest and largest singular values to be returned. |
*> largest singular value to be returned. |
*> 1 <= IL <= IU <= min(M,N), if min(M,N) > 0. |
*> 1 <= IL <= IU <= min(M,N), if min(M,N) > 0. |
*> Not referenced if RANGE = 'A' or 'V'. |
*> Not referenced if RANGE = 'A' or 'V'. |
*> \endverbatim |
*> \endverbatim |
Line 167
|
Line 176
|
*> vectors, stored columnwise) as specified by RANGE; if |
*> vectors, stored columnwise) as specified by RANGE; if |
*> JOBU = 'N', U is not referenced. |
*> JOBU = 'N', U is not referenced. |
*> Note: The user must ensure that UCOL >= NS; if RANGE = 'V', |
*> Note: The user must ensure that UCOL >= NS; if RANGE = 'V', |
*> the exact value of NS is not known ILQFin advance and an upper |
*> the exact value of NS is not known in advance and an upper |
*> bound must be used. |
*> bound must be used. |
*> \endverbatim |
*> \endverbatim |
*> |
*> |
Line 252
|
Line 261
|
*> \author Univ. of Colorado Denver |
*> \author Univ. of Colorado Denver |
*> \author NAG Ltd. |
*> \author NAG Ltd. |
* |
* |
*> \date November 2015 |
*> \date June 2016 |
* |
* |
*> \ingroup complex16GEsing |
*> \ingroup complex16GEsing |
* |
* |
Line 261
|
Line 270
|
$ IL, IU, NS, S, U, LDU, VT, LDVT, WORK, |
$ IL, IU, NS, S, U, LDU, VT, LDVT, WORK, |
$ LWORK, RWORK, IWORK, INFO ) |
$ LWORK, RWORK, IWORK, INFO ) |
* |
* |
* -- LAPACK driver routine (version 3.6.0) -- |
* -- LAPACK driver routine (version 3.6.1) -- |
* -- LAPACK is a software package provided by Univ. of Tennessee, -- |
* -- LAPACK is a software package provided by Univ. of Tennessee, -- |
* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- |
* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- |
* November 2015 |
* June 2016 |
* |
* |
* .. Scalar Arguments .. |
* .. Scalar Arguments .. |
CHARACTER JOBU, JOBVT, RANGE |
CHARACTER JOBU, JOBVT, RANGE |
Line 291
|
Line 300
|
CHARACTER JOBZ, RNGTGK |
CHARACTER JOBZ, RNGTGK |
LOGICAL ALLS, INDS, LQUERY, VALS, WANTU, WANTVT |
LOGICAL ALLS, INDS, LQUERY, VALS, WANTU, WANTVT |
INTEGER I, ID, IE, IERR, ILQF, ILTGK, IQRF, ISCL, |
INTEGER I, ID, IE, IERR, ILQF, ILTGK, IQRF, ISCL, |
$ ITAU, ITAUP, ITAUQ, ITEMP, ITGKZ, IUTGK, |
$ ITAU, ITAUP, ITAUQ, ITEMP, ITEMPR, ITGKZ, |
$ J, K, MAXWRK, MINMN, MINWRK, MNTHR |
$ IUTGK, J, K, MAXWRK, MINMN, MINWRK, MNTHR |
DOUBLE PRECISION ABSTOL, ANRM, BIGNUM, EPS, SMLNUM |
DOUBLE PRECISION ABSTOL, ANRM, BIGNUM, EPS, SMLNUM |
* .. |
* .. |
* .. Local Arrays .. |
* .. Local Arrays .. |
Line 364
|
Line 373
|
IF( INFO.EQ.0 ) THEN |
IF( INFO.EQ.0 ) THEN |
IF( WANTU .AND. LDU.LT.M ) THEN |
IF( WANTU .AND. LDU.LT.M ) THEN |
INFO = -15 |
INFO = -15 |
ELSE IF( WANTVT .AND. LDVT.LT.MINMN ) THEN |
ELSE IF( WANTVT ) THEN |
INFO = -16 |
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 |
END IF |
END IF |
END IF |
Line 387
|
Line 402
|
* |
* |
* Path 1 (M much larger than N) |
* Path 1 (M much larger than N) |
* |
* |
MAXWRK = N + N* |
MINWRK = N*(N+5) |
$ ILAENV( 1, 'DGEQRF', ' ', M, N, -1, -1 ) |
MAXWRK = N + N*ILAENV(1,'ZGEQRF',' ',M,N,-1,-1) |
MAXWRK = MAX( MAXWRK, N*N + N + 2*N* |
MAXWRK = MAX(MAXWRK, |
$ ILAENV( 1, 'DGEBRD', ' ', N, N, -1, -1 ) ) |
$ N*N+2*N+2*N*ILAENV(1,'ZGEBRD',' ',N,N,-1,-1)) |
MINWRK = N*(N+4) |
IF (WANTU .OR. WANTVT) THEN |
|
MAXWRK = MAX(MAXWRK, |
|
$ N*N+2*N+N*ILAENV(1,'ZUNMQR','LN',N,N,N,-1)) |
|
END IF |
ELSE |
ELSE |
* |
* |
* Path 2 (M at least N, but not much larger) |
* Path 2 (M at least N, but not much larger) |
* |
* |
MAXWRK = 2*N + ( M+N )* |
MINWRK = 3*N + M |
$ ILAENV( 1, 'ZGEBRD', ' ', M, N, -1, -1 ) |
MAXWRK = 2*N + (M+N)*ILAENV(1,'ZGEBRD',' ',M,N,-1,-1) |
MINWRK = 2*N + M |
IF (WANTU .OR. WANTVT) THEN |
|
MAXWRK = MAX(MAXWRK, |
|
$ 2*N+N*ILAENV(1,'ZUNMQR','LN',N,N,N,-1)) |
|
END IF |
END IF |
END IF |
ELSE |
ELSE |
MNTHR = ILAENV( 6, 'ZGESVD', JOBU // JOBVT, M, N, 0, 0 ) |
MNTHR = ILAENV( 6, 'ZGESVD', JOBU // JOBVT, M, N, 0, 0 ) |
Line 406
|
Line 427
|
* |
* |
* Path 1t (N much larger than M) |
* Path 1t (N much larger than M) |
* |
* |
MAXWRK = M + M* |
MINWRK = M*(M+5) |
$ ILAENV( 1, 'ZGELQF', ' ', M, N, -1, -1 ) |
MAXWRK = M + M*ILAENV(1,'ZGELQF',' ',M,N,-1,-1) |
MAXWRK = MAX( MAXWRK, M*M + M + 2*M* |
MAXWRK = MAX(MAXWRK, |
$ ILAENV( 1, 'ZGEBRD', ' ', M, M, -1, -1 ) ) |
$ M*M+2*M+2*M*ILAENV(1,'ZGEBRD',' ',M,M,-1,-1)) |
MINWRK = M*(M+4) |
IF (WANTU .OR. WANTVT) THEN |
|
MAXWRK = MAX(MAXWRK, |
|
$ M*M+2*M+M*ILAENV(1,'ZUNMQR','LN',M,M,M,-1)) |
|
END IF |
ELSE |
ELSE |
* |
* |
* Path 2t (N greater than M, but not much larger) |
* Path 2t (N greater than M, but not much larger) |
* |
* |
MAXWRK = M*(M*2+19) + ( M+N )* |
* |
$ ILAENV( 1, 'ZGEBRD', ' ', M, N, -1, -1 ) |
MINWRK = 3*M + N |
MINWRK = 2*M + N |
MAXWRK = 2*M + (M+N)*ILAENV(1,'ZGEBRD',' ',M,N,-1,-1) |
|
IF (WANTU .OR. WANTVT) THEN |
|
MAXWRK = MAX(MAXWRK, |
|
$ 2*M+M*ILAENV(1,'ZUNMQR','LN',M,M,M,-1)) |
|
END IF |
END IF |
END IF |
END IF |
END IF |
END IF |
END IF |
Line 444
|
Line 472
|
* |
* |
* Set singular values indices accord to RANGE='A'. |
* Set singular values indices accord to RANGE='A'. |
* |
* |
ALLS = LSAME( RANGE, 'A' ) |
|
INDS = LSAME( RANGE, 'I' ) |
|
IF( ALLS ) THEN |
IF( ALLS ) THEN |
RNGTGK = 'I' |
RNGTGK = 'I' |
ILTGK = 1 |
ILTGK = 1 |
Line 515
|
Line 541
|
CALL ZGEBRD( N, N, WORK( IQRF ), N, RWORK( ID ), |
CALL ZGEBRD( N, N, WORK( IQRF ), N, RWORK( ID ), |
$ RWORK( IE ), WORK( ITAUQ ), WORK( ITAUP ), |
$ RWORK( IE ), WORK( ITAUQ ), WORK( ITAUP ), |
$ WORK( ITEMP ), LWORK-ITEMP+1, INFO ) |
$ WORK( ITEMP ), LWORK-ITEMP+1, INFO ) |
ITEMP = ITGKZ + N*(N*2+1) |
ITEMPR = ITGKZ + N*(N*2+1) |
* |
* |
* Solve eigenvalue problem TGK*Z=Z*S. |
* Solve eigenvalue problem TGK*Z=Z*S. |
* (Workspace: need 2*N*N+14*N) |
* (Workspace: need 2*N*N+14*N) |
* |
* |
CALL DBDSVDX( 'U', JOBZ, RNGTGK, N, RWORK( ID ), |
CALL DBDSVDX( 'U', JOBZ, RNGTGK, N, RWORK( ID ), |
$ RWORK( IE ), VL, VU, ILTGK, IUTGK, NS, S, |
$ RWORK( IE ), VL, VU, ILTGK, IUTGK, NS, S, |
$ RWORK( ITGKZ ), N*2, RWORK( ITEMP ), |
$ RWORK( ITGKZ ), N*2, RWORK( ITEMPR ), |
$ IWORK, INFO) |
$ IWORK, INFO) |
* |
* |
* If needed, compute left singular vectors. |
* If needed, compute left singular vectors. |
Line 536
|
Line 562
|
END DO |
END DO |
K = K + N |
K = K + N |
END DO |
END DO |
CALL ZLASET( 'A', M-N, N, CZERO, CZERO, U( N+1,1 ), LDU ) |
CALL ZLASET( 'A', M-N, NS, CZERO, CZERO, U( N+1,1 ), LDU) |
* |
* |
* Call ZUNMBR to compute QB*UB. |
* Call ZUNMBR to compute QB*UB. |
* (Workspace in WORK( ITEMP ): need N, prefer N*NB) |
* (Workspace in WORK( ITEMP ): need N, prefer N*NB) |
Line 591
|
Line 617
|
CALL ZGEBRD( M, N, A, LDA, RWORK( ID ), RWORK( IE ), |
CALL ZGEBRD( M, N, A, LDA, RWORK( ID ), RWORK( IE ), |
$ WORK( ITAUQ ), WORK( ITAUP ), WORK( ITEMP ), |
$ WORK( ITAUQ ), WORK( ITAUP ), WORK( ITEMP ), |
$ LWORK-ITEMP+1, INFO ) |
$ LWORK-ITEMP+1, INFO ) |
ITEMP = ITGKZ + N*(N*2+1) |
ITEMPR = ITGKZ + N*(N*2+1) |
* |
* |
* Solve eigenvalue problem TGK*Z=Z*S. |
* Solve eigenvalue problem TGK*Z=Z*S. |
* (Workspace: need 2*N*N+14*N) |
* (Workspace: need 2*N*N+14*N) |
* |
* |
CALL DBDSVDX( 'U', JOBZ, RNGTGK, N, RWORK( ID ), |
CALL DBDSVDX( 'U', JOBZ, RNGTGK, N, RWORK( ID ), |
$ RWORK( IE ), VL, VU, ILTGK, IUTGK, NS, S, |
$ RWORK( IE ), VL, VU, ILTGK, IUTGK, NS, S, |
$ RWORK( ITGKZ ), N*2, RWORK( ITEMP ), |
$ RWORK( ITGKZ ), N*2, RWORK( ITEMPR ), |
$ IWORK, INFO) |
$ IWORK, INFO) |
* |
* |
* If needed, compute left singular vectors. |
* If needed, compute left singular vectors. |
Line 612
|
Line 638
|
END DO |
END DO |
K = K + N |
K = K + N |
END DO |
END DO |
CALL ZLASET( 'A', M-N, N, CZERO, CZERO, U( N+1,1 ), LDU ) |
CALL ZLASET( 'A', M-N, NS, CZERO, CZERO, U( N+1,1 ), LDU) |
* |
* |
* Call ZUNMBR to compute QB*UB. |
* Call ZUNMBR to compute QB*UB. |
* (Workspace in WORK( ITEMP ): need N, prefer N*NB) |
* (Workspace in WORK( ITEMP ): need N, prefer N*NB) |
Line 678
|
Line 704
|
CALL ZGEBRD( M, M, WORK( ILQF ), M, RWORK( ID ), |
CALL ZGEBRD( M, M, WORK( ILQF ), M, RWORK( ID ), |
$ RWORK( IE ), WORK( ITAUQ ), WORK( ITAUP ), |
$ RWORK( IE ), WORK( ITAUQ ), WORK( ITAUP ), |
$ WORK( ITEMP ), LWORK-ITEMP+1, INFO ) |
$ WORK( ITEMP ), LWORK-ITEMP+1, INFO ) |
ITEMP = ITGKZ + M*(M*2+1) |
ITEMPR = ITGKZ + M*(M*2+1) |
* |
* |
* Solve eigenvalue problem TGK*Z=Z*S. |
* Solve eigenvalue problem TGK*Z=Z*S. |
* (Workspace: need 2*M*M+14*M) |
* (Workspace: need 2*M*M+14*M) |
* |
* |
CALL DBDSVDX( 'U', JOBZ, RNGTGK, M, RWORK( ID ), |
CALL DBDSVDX( 'U', JOBZ, RNGTGK, M, RWORK( ID ), |
$ RWORK( IE ), VL, VU, ILTGK, IUTGK, NS, S, |
$ RWORK( IE ), VL, VU, ILTGK, IUTGK, NS, S, |
$ RWORK( ITGKZ ), M*2, RWORK( ITEMP ), |
$ RWORK( ITGKZ ), M*2, RWORK( ITEMPR ), |
$ IWORK, INFO) |
$ IWORK, INFO) |
* |
* |
* If needed, compute left singular vectors. |
* If needed, compute left singular vectors. |
Line 719
|
Line 745
|
END DO |
END DO |
K = K + M |
K = K + M |
END DO |
END DO |
CALL ZLASET( 'A', M, N-M, CZERO, CZERO, |
CALL ZLASET( 'A', NS, N-M, CZERO, CZERO, |
$ VT( 1,M+1 ), LDVT ) |
$ VT( 1,M+1 ), LDVT ) |
* |
* |
* Call ZUNMBR to compute (VB**T)*(PB**T) |
* Call ZUNMBR to compute (VB**T)*(PB**T) |
Line 755
|
Line 781
|
CALL ZGEBRD( M, N, A, LDA, RWORK( ID ), RWORK( IE ), |
CALL ZGEBRD( M, N, A, LDA, RWORK( ID ), RWORK( IE ), |
$ WORK( ITAUQ ), WORK( ITAUP ), WORK( ITEMP ), |
$ WORK( ITAUQ ), WORK( ITAUP ), WORK( ITEMP ), |
$ LWORK-ITEMP+1, INFO ) |
$ LWORK-ITEMP+1, INFO ) |
ITEMP = ITGKZ + M*(M*2+1) |
ITEMPR = ITGKZ + M*(M*2+1) |
* |
* |
* Solve eigenvalue problem TGK*Z=Z*S. |
* Solve eigenvalue problem TGK*Z=Z*S. |
* (Workspace: need 2*M*M+14*M) |
* (Workspace: need 2*M*M+14*M) |
* |
* |
CALL DBDSVDX( 'L', JOBZ, RNGTGK, M, RWORK( ID ), |
CALL DBDSVDX( 'L', JOBZ, RNGTGK, M, RWORK( ID ), |
$ RWORK( IE ), VL, VU, ILTGK, IUTGK, NS, S, |
$ RWORK( IE ), VL, VU, ILTGK, IUTGK, NS, S, |
$ RWORK( ITGKZ ), M*2, RWORK( ITEMP ), |
$ RWORK( ITGKZ ), M*2, RWORK( ITEMPR ), |
$ IWORK, INFO) |
$ IWORK, INFO) |
* |
* |
* If needed, compute left singular vectors. |
* If needed, compute left singular vectors. |
Line 796
|
Line 822
|
END DO |
END DO |
K = K + M |
K = K + M |
END DO |
END DO |
CALL ZLASET( 'A', M, N-M, CZERO, CZERO, |
CALL ZLASET( 'A', NS, N-M, CZERO, CZERO, |
$ VT( 1,M+1 ), LDVT ) |
$ VT( 1,M+1 ), LDVT ) |
* |
* |
* Call ZUNMBR to compute VB**T * PB**T |
* Call ZUNMBR to compute VB**T * PB**T |