version 1.1, 2015/11/26 11:44:21
|
version 1.5, 2017/06/17 11:06:44
|
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( * ) |
* .. |
* .. |
* |
* |
* |
* |
* Purpose |
*> \par Purpose: |
* ======= |
* ============= |
|
*> |
|
*> \verbatim |
|
*> |
|
*> ZGESVDX computes the singular value decomposition (SVD) of a complex |
|
*> 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 unitary matrix, and |
|
*> V is an N-by-N unitary 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. |
|
*> |
|
*> ZGESVDX 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 |
* |
* |
* ZGESVDX computes the singular value decomposition (SVD) of a complex |
|
* 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 unitary matrix, and |
|
* V is an N-by-N unitary 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. |
|
* |
|
* ZGESVDX 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. |
|
* |
|
* Arguments: |
* Arguments: |
* ========== |
* ========== |
* |
* |
Line 66
|
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 78
|
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 90
|
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 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 149
|
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 163
|
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 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 181
|
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 206
|
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 228
|
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 247
|
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 November 2015 |
*> \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.0) -- |
* -- LAPACK driver routine (version 3.7.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..-- |
* November 2015 |
* June 2016 |
* |
* |
* .. Scalar Arguments .. |
* .. Scalar Arguments .. |
CHARACTER JOBU, JOBVT, RANGE |
CHARACTER JOBU, JOBVT, RANGE |
Line 274
|
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 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 454
|
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 498
|
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 ) |
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) |
* |
* |
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 568
|
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 584
|
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 ) |
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 606
|
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 |
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) |
* |
* |
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 637
|
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 650
|
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 665
|
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 673
|
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 ) |
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 703
|
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 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) |
* (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 ) |
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. |
* |
* |
IF( WANTU ) THEN |
IF( WANTU ) THEN |
Line 780
|
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 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 |
* (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 |
* |
* |