version 1.3, 2016/08/27 15:34:47
|
version 1.6, 2018/05/29 06:55:22
|
Line 2
|
Line 2
|
* |
* |
* =========== DOCUMENTATION =========== |
* =========== DOCUMENTATION =========== |
* |
* |
* Online html documentation available at |
* Online html documentation available at |
* http://www.netlib.org/lapack/explore-html/ |
* http://www.netlib.org/lapack/explore-html/ |
* |
* |
*> \htmlonly |
*> \htmlonly |
*> Download ZGESVDX + dependencies |
*> Download ZGESVDX + dependencies |
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zgesvdx.f"> |
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zgesvdx.f"> |
*> [TGZ]</a> |
*> [TGZ]</a> |
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zgesvdx.f"> |
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zgesvdx.f"> |
*> [ZIP]</a> |
*> [ZIP]</a> |
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zgesvdx.f"> |
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zgesvdx.f"> |
*> [TXT]</a> |
*> [TXT]</a> |
*> \endhtmlonly |
*> \endhtmlonly |
* |
* |
* Definition: |
* Definition: |
* =========== |
* =========== |
* |
* |
* SUBROUTINE CGESVDX( JOBU, JOBVT, RANGE, M, N, A, LDA, VL, VU, |
* SUBROUTINE CGESVDX( JOBU, JOBVT, RANGE, M, N, A, LDA, VL, VU, |
* $ 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 ) |
* |
* |
* |
* |
* .. Scalar Arguments .. |
* .. Scalar Arguments .. |
* CHARACTER JOBU, JOBVT, RANGE |
* CHARACTER JOBU, JOBVT, RANGE |
Line 31
|
Line 31
|
* .. Array Arguments .. |
* .. Array Arguments .. |
* INTEGER IWORK( * ) |
* INTEGER IWORK( * ) |
* DOUBLE PRECISION S( * ), RWORK( * ) |
* DOUBLE PRECISION S( * ), RWORK( * ) |
* COMPLEX*16 A( LDA, * ), U( LDU, * ), VT( LDVT, * ), |
* COMPLEX*16 A( LDA, * ), U( LDU, * ), VT( LDVT, * ), |
* $ WORK( * ) |
* $ WORK( * ) |
* .. |
* .. |
* |
* |
Line 60
|
Line 60
|
*> |
*> |
*> Note that the routine returns V**T, not V. |
*> Note that the routine returns V**T, not V. |
*> \endverbatim |
*> \endverbatim |
* |
* |
* Arguments: |
* Arguments: |
* ========== |
* ========== |
* |
* |
Line 69
|
Line 69
|
*> JOBU is CHARACTER*1 |
*> JOBU is CHARACTER*1 |
*> Specifies options for computing all or part of the matrix U: |
*> Specifies options for computing all or part of the matrix U: |
*> = 'V': the first min(m,n) columns of U (the left singular |
*> = '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; |
*> the array U; |
*> = 'N': no columns of U (no left singular vectors) are |
*> = 'N': no columns of U (no left singular vectors) are |
*> computed. |
*> computed. |
Line 81
|
Line 81
|
*> Specifies options for computing all or part of the matrix |
*> Specifies options for computing all or part of the matrix |
*> V**T: |
*> V**T: |
*> = 'V': the first min(m,n) rows of V**T (the right singular |
*> = '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; |
*> the array VT; |
*> = 'N': no rows of V**T (no right singular vectors) are |
*> = 'N': no rows of V**T (no right singular vectors) are |
*> computed. |
*> computed. |
Line 93
|
Line 93
|
*> = 'A': all singular values will be found. |
*> = 'A': all singular values will be found. |
*> = 'V': all singular values in the half-open interval (VL,VU] |
*> = 'V': all singular values in the half-open interval (VL,VU] |
*> will be found. |
*> 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 |
*> \endverbatim |
*> |
*> |
*> \param[in] M |
*> \param[in] M |
Line 158
|
Line 158
|
*> \param[out] NS |
*> \param[out] NS |
*> \verbatim |
*> \verbatim |
*> NS is INTEGER |
*> NS is INTEGER |
*> The total number of singular values found, |
*> The total number of singular values found, |
*> 0 <= NS <= min(M,N). |
*> 0 <= NS <= min(M,N). |
*> If RANGE = 'A', NS = min(M,N); if RANGE = 'I', NS = IU-IL+1. |
*> If RANGE = 'A', NS = min(M,N); if RANGE = 'I', NS = IU-IL+1. |
*> \endverbatim |
*> \endverbatim |
Line 172
|
Line 172
|
*> \param[out] U |
*> \param[out] U |
*> \verbatim |
*> \verbatim |
*> U is COMPLEX*16 array, dimension (LDU,UCOL) |
*> U is COMPLEX*16 array, dimension (LDU,UCOL) |
*> If JOBU = 'V', U contains columns of U (the left singular |
*> If JOBU = 'V', U contains columns of U (the left singular |
*> 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 in 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 190
|
Line 190
|
*> \param[out] VT |
*> \param[out] VT |
*> \verbatim |
*> \verbatim |
*> VT is COMPLEX*16 array, dimension (LDVT,N) |
*> VT is COMPLEX*16 array, dimension (LDVT,N) |
*> If JOBVT = 'V', VT contains the rows of V**T (the right singular |
*> If JOBVT = 'V', VT contains the rows of V**T (the right singular |
*> vectors, stored rowwise) as specified by RANGE; if JOBVT = 'N', |
*> vectors, stored rowwise) as specified by RANGE; if JOBVT = 'N', |
*> VT is not referenced. |
*> VT is not referenced. |
*> Note: The user must ensure that LDVT >= NS; if RANGE = 'V', |
*> Note: The user must ensure that LDVT >= NS; if RANGE = 'V', |
*> the exact value of NS is not known in 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 215
|
Line 215
|
*> \verbatim |
*> \verbatim |
*> LWORK is INTEGER |
*> LWORK is INTEGER |
*> The dimension of the array WORK. |
*> 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): |
*> 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) |
*> - PATH 1t (N much larger than M) |
*> LWORK >= MAX(1,MIN(M,N)*2+MAX(M,N)) for the other paths. |
*> LWORK >= MAX(1,MIN(M,N)*2+MAX(M,N)) for the other paths. |
*> For good performance, LWORK should generally be larger. |
*> For good performance, LWORK should generally be larger. |
Line 237
|
Line 237
|
*> \param[out] IWORK |
*> \param[out] IWORK |
*> \verbatim |
*> \verbatim |
*> IWORK is INTEGER array, dimension (12*MIN(M,N)) |
*> IWORK is INTEGER array, dimension (12*MIN(M,N)) |
*> If INFO = 0, the first NS elements of IWORK are zero. If INFO > 0, |
*> If INFO = 0, the first NS elements of IWORK are zero. If INFO > 0, |
*> then IWORK contains the indices of the eigenvectors that failed |
*> then IWORK contains the indices of the eigenvectors that failed |
*> to converge in DBDSVDX/DSTEVX. |
*> to converge in DBDSVDX/DSTEVX. |
*> \endverbatim |
*> \endverbatim |
*> |
*> |
Line 256
|
Line 256
|
* Authors: |
* Authors: |
* ======== |
* ======== |
* |
* |
*> \author Univ. of Tennessee |
*> \author Univ. of Tennessee |
*> \author Univ. of California Berkeley |
*> \author Univ. of California Berkeley |
*> \author Univ. of Colorado Denver |
*> \author Univ. of Colorado Denver |
*> \author NAG Ltd. |
*> \author NAG Ltd. |
* |
* |
*> \date June 2016 |
*> \date June 2016 |
* |
* |
*> \ingroup complex16GEsing |
*> \ingroup complex16GEsing |
* |
* |
* ===================================================================== |
* ===================================================================== |
SUBROUTINE ZGESVDX( JOBU, JOBVT, RANGE, M, N, A, LDA, VL, VU, |
SUBROUTINE ZGESVDX( JOBU, JOBVT, RANGE, M, N, A, LDA, VL, VU, |
$ 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.1) -- |
* -- LAPACK driver routine (version 3.8.0) -- |
* -- 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..-- |
* June 2016 |
* June 2016 |
Line 283
|
Line 283
|
* .. Array Arguments .. |
* .. Array Arguments .. |
INTEGER IWORK( * ) |
INTEGER IWORK( * ) |
DOUBLE PRECISION S( * ), RWORK( * ) |
DOUBLE PRECISION S( * ), RWORK( * ) |
COMPLEX*16 A( LDA, * ), U( LDU, * ), VT( LDVT, * ), |
COMPLEX*16 A( LDA, * ), U( LDU, * ), VT( LDVT, * ), |
$ WORK( * ) |
$ WORK( * ) |
* .. |
* .. |
* |
* |
Line 300
|
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, ITEMPR, ITGKZ, |
$ ITAU, ITAUP, ITAUQ, ITEMP, ITEMPR, ITGKZ, |
$ IUTGK, 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 |
* .. |
* .. |
Line 308
|
Line 308
|
DOUBLE PRECISION DUM( 1 ) |
DOUBLE PRECISION DUM( 1 ) |
* .. |
* .. |
* .. External Subroutines .. |
* .. External Subroutines .. |
EXTERNAL ZGEBRD, ZGELQF, ZGEQRF, ZLASCL, ZLASET, |
EXTERNAL ZGEBRD, ZGELQF, ZGEQRF, ZLASCL, ZLASET, ZLACPY, |
$ DLASCL, XERBLA |
$ ZUNMLQ, ZUNMBR, ZUNMQR, DBDSVDX, DLASCL, XERBLA |
* .. |
* .. |
* .. External Functions .. |
* .. External Functions .. |
LOGICAL LSAME |
LOGICAL LSAME |
Line 480
|
Line 480
|
RNGTGK = 'I' |
RNGTGK = 'I' |
ILTGK = IL |
ILTGK = IL |
IUTGK = IU |
IUTGK = IU |
ELSE |
ELSE |
RNGTGK = 'V' |
RNGTGK = 'V' |
ILTGK = 0 |
ILTGK = 0 |
IUTGK = 0 |
IUTGK = 0 |
Line 524
|
Line 524
|
ITEMP = ITAU + N |
ITEMP = ITAU + N |
CALL ZGEQRF( M, N, A, LDA, WORK( ITAU ), WORK( ITEMP ), |
CALL ZGEQRF( M, N, A, LDA, WORK( ITAU ), WORK( ITEMP ), |
$ LWORK-ITEMP+1, INFO ) |
$ LWORK-ITEMP+1, INFO ) |
* |
* |
* Copy R into WORK and bidiagonalize it: |
* Copy R into WORK and bidiagonalize it: |
* (Workspace: need N*N+3*N, prefer N*N+N+2*N*NB) |
* (Workspace: need N*N+3*N, prefer N*N+N+2*N*NB) |
* |
* |
IQRF = ITEMP |
IQRF = ITEMP |
ITAUQ = ITEMP + N*N |
ITAUQ = ITEMP + N*N |
ITAUP = ITAUQ + N |
ITAUP = ITAUQ + N |
ITEMP = ITAUP + N |
ITEMP = ITAUP + N |
ID = 1 |
ID = 1 |
IE = ID + N |
IE = ID + N |
ITGKZ = IE + N |
ITGKZ = IE + N |
CALL ZLACPY( 'U', N, N, A, LDA, WORK( IQRF ), N ) |
CALL ZLACPY( 'U', N, N, A, LDA, WORK( IQRF ), N ) |
CALL ZLASET( 'L', N-1, N-1, CZERO, CZERO, |
CALL ZLASET( 'L', N-1, N-1, CZERO, CZERO, |
$ WORK( IQRF+1 ), N ) |
$ WORK( IQRF+1 ), N ) |
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 ) |
ITEMPR = 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( ITEMPR ), |
$ RWORK( ITGKZ ), N*2, RWORK( ITEMPR ), |
Line 567
|
Line 567
|
* 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) |
* |
* |
CALL ZUNMBR( 'Q', 'L', 'N', N, NS, N, WORK( IQRF ), N, |
CALL ZUNMBR( 'Q', 'L', 'N', N, NS, N, WORK( IQRF ), N, |
$ WORK( ITAUQ ), U, LDU, WORK( ITEMP ), |
$ WORK( ITAUQ ), U, LDU, WORK( ITEMP ), |
$ LWORK-ITEMP+1, INFO ) |
$ LWORK-ITEMP+1, INFO ) |
* |
* |
* Call ZUNMQR to compute Q*(QB*UB). |
* Call ZUNMQR to compute Q*(QB*UB). |
* (Workspace in WORK( ITEMP ): need N, prefer N*NB) |
* (Workspace in WORK( ITEMP ): need N, prefer N*NB) |
* |
* |
CALL ZUNMQR( 'L', 'N', M, NS, N, A, LDA, |
CALL ZUNMQR( 'L', 'N', M, NS, N, A, LDA, |
$ WORK( ITAU ), U, LDU, WORK( ITEMP ), |
$ WORK( ITAU ), U, LDU, WORK( ITEMP ), |
$ LWORK-ITEMP+1, INFO ) |
$ LWORK-ITEMP+1, INFO ) |
END IF |
END IF |
* |
* |
* If needed, compute right singular vectors. |
* If needed, compute right singular vectors. |
* |
* |
IF( WANTVT) THEN |
IF( WANTVT) THEN |
Line 594
|
Line 594
|
* Call ZUNMBR to compute VB**T * PB**T |
* Call ZUNMBR to compute VB**T * PB**T |
* (Workspace in WORK( ITEMP ): need N, prefer N*NB) |
* (Workspace in WORK( ITEMP ): need N, prefer N*NB) |
* |
* |
CALL ZUNMBR( 'P', 'R', 'C', NS, N, N, WORK( IQRF ), N, |
CALL ZUNMBR( 'P', 'R', 'C', NS, N, N, WORK( IQRF ), N, |
$ WORK( ITAUP ), VT, LDVT, WORK( ITEMP ), |
$ WORK( ITAUP ), VT, LDVT, WORK( ITEMP ), |
$ LWORK-ITEMP+1, INFO ) |
$ LWORK-ITEMP+1, INFO ) |
END IF |
END IF |
Line 610
|
Line 610
|
* |
* |
ITAUQ = 1 |
ITAUQ = 1 |
ITAUP = ITAUQ + N |
ITAUP = ITAUQ + N |
ITEMP = ITAUP + N |
ITEMP = ITAUP + N |
ID = 1 |
ID = 1 |
IE = ID + N |
IE = ID + N |
ITGKZ = IE + N |
ITGKZ = IE + N |
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 ) |
ITEMPR = 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( ITEMPR ), |
$ RWORK( ITGKZ ), N*2, RWORK( ITEMPR ), |
$ IWORK, INFO) |
$ IWORK, INFO) |
* |
* |
* If needed, compute left singular vectors. |
* If needed, compute left singular vectors. |
Line 632
|
Line 632
|
IF( WANTU ) THEN |
IF( WANTU ) THEN |
K = ITGKZ |
K = ITGKZ |
DO I = 1, NS |
DO I = 1, NS |
DO J = 1, N |
DO J = 1, N |
U( J, I ) = DCMPLX( RWORK( K ), ZERO ) |
U( J, I ) = DCMPLX( RWORK( K ), ZERO ) |
K = K + 1 |
K = K + 1 |
END DO |
END DO |
Line 642
|
Line 642
|
* |
* |
* 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) |
* |
* |
CALL ZUNMBR( 'Q', 'L', 'N', M, NS, N, A, LDA, |
CALL ZUNMBR( 'Q', 'L', 'N', M, NS, N, A, LDA, |
$ WORK( ITAUQ ), U, LDU, WORK( ITEMP ), |
$ WORK( ITAUQ ), U, LDU, WORK( ITEMP ), |
$ LWORK-ITEMP+1, IERR ) |
$ LWORK-ITEMP+1, IERR ) |
END IF |
END IF |
* |
* |
* If needed, compute right singular vectors. |
* If needed, compute right singular vectors. |
* |
* |
IF( WANTVT) THEN |
IF( WANTVT) THEN |
Line 663
|
Line 663
|
* Call ZUNMBR to compute VB**T * PB**T |
* Call ZUNMBR to compute VB**T * PB**T |
* (Workspace in WORK( ITEMP ): need N, prefer N*NB) |
* (Workspace in WORK( ITEMP ): need N, prefer N*NB) |
* |
* |
CALL ZUNMBR( 'P', 'R', 'C', NS, N, N, A, LDA, |
CALL ZUNMBR( 'P', 'R', 'C', NS, N, N, A, LDA, |
$ WORK( ITAUP ), VT, LDVT, WORK( ITEMP ), |
$ WORK( ITAUP ), VT, LDVT, WORK( ITEMP ), |
$ LWORK-ITEMP+1, IERR ) |
$ LWORK-ITEMP+1, IERR ) |
END IF |
END IF |
END IF |
END IF |
ELSE |
ELSE |
* |
* |
* A has more columns than rows. If A has sufficiently more |
* A has more columns than rows. If A has sufficiently more |
Line 676
|
Line 676
|
IF( N.GE.MNTHR ) THEN |
IF( N.GE.MNTHR ) THEN |
* |
* |
* Path 1t (N much larger than M): |
* 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 |
* = ( QB * ( UB * S * VB**T ) * PB**T ) * Q |
* U = QB * UB ; V**T = VB**T * PB**T * Q |
* U = QB * UB ; V**T = VB**T * PB**T * Q |
* |
* |
Line 691
|
Line 691
|
* Copy L into WORK and bidiagonalize it: |
* Copy L into WORK and bidiagonalize it: |
* (Workspace in WORK( ITEMP ): need M*M+3*M, prefer M*M+M+2*M*NB) |
* (Workspace in WORK( ITEMP ): need M*M+3*M, prefer M*M+M+2*M*NB) |
* |
* |
ILQF = ITEMP |
ILQF = ITEMP |
ITAUQ = ILQF + M*M |
ITAUQ = ILQF + M*M |
ITAUP = ITAUQ + M |
ITAUP = ITAUQ + M |
ITEMP = ITAUP + M |
ITEMP = ITAUP + M |
Line 699
|
Line 699
|
IE = ID + M |
IE = ID + M |
ITGKZ = IE + M |
ITGKZ = IE + M |
CALL ZLACPY( 'L', M, M, A, LDA, WORK( ILQF ), M ) |
CALL ZLACPY( 'L', M, M, A, LDA, WORK( ILQF ), M ) |
CALL ZLASET( 'U', M-1, M-1, CZERO, CZERO, |
CALL ZLASET( 'U', M-1, M-1, CZERO, CZERO, |
$ WORK( ILQF+M ), M ) |
$ WORK( ILQF+M ), M ) |
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 ) |
ITEMPR = 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( ITEMPR ), |
$ RWORK( ITGKZ ), M*2, RWORK( ITEMPR ), |
$ IWORK, INFO) |
$ IWORK, INFO) |
* |
* |
* If needed, compute left singular vectors. |
* If needed, compute left singular vectors. |
Line 729
|
Line 729
|
* Call ZUNMBR to compute QB*UB. |
* Call ZUNMBR to compute QB*UB. |
* (Workspace in WORK( ITEMP ): need M, prefer M*NB) |
* (Workspace in WORK( ITEMP ): need M, prefer M*NB) |
* |
* |
CALL ZUNMBR( 'Q', 'L', 'N', M, NS, M, WORK( ILQF ), M, |
CALL ZUNMBR( 'Q', 'L', 'N', M, NS, M, WORK( ILQF ), M, |
$ WORK( ITAUQ ), U, LDU, WORK( ITEMP ), |
$ WORK( ITAUQ ), U, LDU, WORK( ITEMP ), |
$ LWORK-ITEMP+1, INFO ) |
$ LWORK-ITEMP+1, INFO ) |
END IF |
END IF |
* |
* |
* If needed, compute right singular vectors. |
* If needed, compute right singular vectors. |
* |
* |
IF( WANTVT) THEN |
IF( WANTVT) THEN |
Line 751
|
Line 751
|
* Call ZUNMBR to compute (VB**T)*(PB**T) |
* Call ZUNMBR to compute (VB**T)*(PB**T) |
* (Workspace in WORK( ITEMP ): need M, prefer M*NB) |
* (Workspace in WORK( ITEMP ): need M, prefer M*NB) |
* |
* |
CALL ZUNMBR( 'P', 'R', 'C', NS, M, M, WORK( ILQF ), M, |
CALL ZUNMBR( 'P', 'R', 'C', NS, M, M, WORK( ILQF ), M, |
$ WORK( ITAUP ), VT, LDVT, WORK( ITEMP ), |
$ WORK( ITAUP ), VT, LDVT, WORK( ITEMP ), |
$ LWORK-ITEMP+1, INFO ) |
$ LWORK-ITEMP+1, INFO ) |
* |
* |
* Call ZUNMLQ to compute ((VB**T)*(PB**T))*Q. |
* Call ZUNMLQ to compute ((VB**T)*(PB**T))*Q. |
* (Workspace in WORK( ITEMP ): need M, prefer M*NB) |
* (Workspace in WORK( ITEMP ): need M, prefer M*NB) |
* |
* |
CALL ZUNMLQ( 'R', 'N', NS, N, M, A, LDA, |
CALL ZUNMLQ( 'R', 'N', NS, N, M, A, LDA, |
$ WORK( ITAU ), VT, LDVT, WORK( ITEMP ), |
$ WORK( ITAU ), VT, LDVT, WORK( ITEMP ), |
$ LWORK-ITEMP+1, INFO ) |
$ LWORK-ITEMP+1, INFO ) |
END IF |
END IF |
ELSE |
ELSE |
* |
* |
* Path 2t (N greater than M, but not much larger) |
* Path 2t (N greater than M, but not much larger) |
* Reduce to bidiagonal form without LQ decomposition |
* Reduce to bidiagonal form without LQ decomposition |
* A = QB * B * PB**T = QB * ( UB * S * VB**T ) * PB**T |
* 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 |
* Bidiagonalize A |
* (Workspace: need 2*M+N, prefer 2*M+(M+N)*NB) |
* (Workspace: need 2*M+N, prefer 2*M+(M+N)*NB) |
* |
* |
ITAUQ = 1 |
ITAUQ = 1 |
ITAUP = ITAUQ + M |
ITAUP = ITAUQ + M |
ITEMP = ITAUP + M |
ITEMP = ITAUP + M |
ID = 1 |
ID = 1 |
IE = ID + M |
IE = ID + M |
ITGKZ = IE + M |
ITGKZ = IE + M |
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 ) |
ITEMPR = 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( ITEMPR ), |
$ RWORK( ITGKZ ), M*2, RWORK( ITEMPR ), |
$ IWORK, INFO) |
$ IWORK, INFO) |
* |
* |
* If needed, compute left singular vectors. |
* If needed, compute left singular vectors. |
* |
* |
IF( WANTU ) THEN |
IF( WANTU ) THEN |
Line 806
|
Line 806
|
* Call ZUNMBR to compute QB*UB. |
* Call ZUNMBR to compute QB*UB. |
* (Workspace in WORK( ITEMP ): need M, prefer M*NB) |
* (Workspace in WORK( ITEMP ): need M, prefer M*NB) |
* |
* |
CALL ZUNMBR( 'Q', 'L', 'N', M, NS, N, A, LDA, |
CALL ZUNMBR( 'Q', 'L', 'N', M, NS, N, A, LDA, |
$ WORK( ITAUQ ), U, LDU, WORK( ITEMP ), |
$ WORK( ITAUQ ), U, LDU, WORK( ITEMP ), |
$ LWORK-ITEMP+1, INFO ) |
$ LWORK-ITEMP+1, INFO ) |
END IF |
END IF |
* |
* |
* If needed, compute right singular vectors. |
* If needed, compute right singular vectors. |
* |
* |
IF( WANTVT) THEN |
IF( WANTVT) THEN |
Line 828
|
Line 828
|
* Call ZUNMBR to compute VB**T * PB**T |
* Call ZUNMBR to compute VB**T * PB**T |
* (Workspace in WORK( ITEMP ): need M, prefer M*NB) |
* (Workspace in WORK( ITEMP ): need M, prefer M*NB) |
* |
* |
CALL ZUNMBR( 'P', 'R', 'C', NS, N, M, A, LDA, |
CALL ZUNMBR( 'P', 'R', 'C', NS, N, M, A, LDA, |
$ WORK( ITAUP ), VT, LDVT, WORK( ITEMP ), |
$ WORK( ITAUP ), VT, LDVT, WORK( ITEMP ), |
$ LWORK-ITEMP+1, INFO ) |
$ LWORK-ITEMP+1, INFO ) |
END IF |
END IF |
END IF |
END IF |
END IF |
END IF |
* |
* |