version 1.1, 2015/11/26 11:44:15
|
version 1.8, 2023/08/07 08:38:50
|
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 DGESVDX + dependencies |
*> Download DGESVDX + dependencies |
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dgesvdx.f"> |
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dgesvdx.f"> |
*> [TGZ]</a> |
*> [TGZ]</a> |
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dgesvdx.f"> |
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dgesvdx.f"> |
*> [ZIP]</a> |
*> [ZIP]</a> |
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dgesvdx.f"> |
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dgesvdx.f"> |
*> [TXT]</a> |
*> [TXT]</a> |
*> \endhtmlonly |
*> \endhtmlonly |
* |
* |
* Definition: |
* Definition: |
* =========== |
* =========== |
* |
* |
* SUBROUTINE DGESVDX( JOBU, JOBVT, RANGE, M, N, A, LDA, VL, VU, |
* SUBROUTINE DGESVDX( 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, IWORK, INFO ) |
* $ LWORK, IWORK, INFO ) |
* |
* |
* |
* |
* .. Scalar Arguments .. |
* .. Scalar Arguments .. |
* CHARACTER JOBU, JOBVT, RANGE |
* CHARACTER JOBU, JOBVT, RANGE |
Line 33
|
Line 33
|
* DOUBLE PRECISION A( LDA, * ), S( * ), U( LDU, * ), |
* DOUBLE PRECISION A( LDA, * ), S( * ), U( LDU, * ), |
* $ VT( LDVT, * ), WORK( * ) |
* $ VT( LDVT, * ), WORK( * ) |
* .. |
* .. |
* |
* |
* |
* |
*> \par Purpose: |
*> \par Purpose: |
* ============= |
* ============= |
Line 43
|
Line 43
|
*> DGESVDX computes the singular value decomposition (SVD) of a real |
*> DGESVDX computes the singular value decomposition (SVD) of a real |
*> M-by-N matrix A, optionally computing the left and/or right singular |
*> M-by-N matrix A, optionally computing the left and/or right singular |
*> vectors. The SVD is written |
*> vectors. The SVD is written |
*> |
*> |
*> A = U * SIGMA * transpose(V) |
*> A = U * SIGMA * transpose(V) |
*> |
*> |
*> where SIGMA is an M-by-N matrix which is zero except for its |
*> 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 |
*> 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 |
*> 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 the singular values of A; they are real and non-negative, and |
*> are returned in descending order. The first min(m,n) columns of |
*> are returned in descending order. The first min(m,n) columns of |
*> U and V are the left and right singular vectors of A. |
*> U and V are the left and right singular vectors of A. |
*> |
*> |
*> DGESVDX uses an eigenvalue problem for obtaining the SVD, which |
*> DGESVDX uses an eigenvalue problem for obtaining the SVD, which |
*> allows for the computation of a subset of singular values and |
*> allows for the computation of a subset of singular values and |
*> vectors. See DBDSVDX for details. |
*> vectors. See DBDSVDX for details. |
*> |
*> |
*> Note that the routine returns V**T, not V. |
*> Note that the routine returns V**T, not V. |
*> \endverbatim |
*> \endverbatim |
* |
* |
* Arguments: |
* Arguments: |
* ========== |
* ========== |
* |
* |
Line 68
|
Line 68
|
*> 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 80
|
Line 80
|
*> 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 92
|
Line 92
|
*> = '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 123
|
Line 123
|
*> \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 137
|
Line 139
|
*> \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 151
|
Line 157
|
*> \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 165
|
Line 171
|
*> \param[out] U |
*> \param[out] U |
*> \verbatim |
*> \verbatim |
*> U is DOUBLE PRECISION array, dimension (LDU,UCOL) |
*> U is DOUBLE PRECISION 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 183
|
Line 189
|
*> \param[out] VT |
*> \param[out] VT |
*> \verbatim |
*> \verbatim |
*> VT is DOUBLE PRECISION array, dimension (LDVT,N) |
*> VT is DOUBLE PRECISION 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 208
|
Line 214
|
*> \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 224
|
Line 230
|
*> \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 243
|
Line 249
|
* 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 |
|
* |
* |
*> \ingroup doubleGEsing |
*> \ingroup doubleGEsing |
* |
* |
* ===================================================================== |
* ===================================================================== |
SUBROUTINE DGESVDX( JOBU, JOBVT, RANGE, M, N, A, LDA, VL, VU, |
SUBROUTINE DGESVDX( 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, IWORK, INFO ) |
$ LWORK, IWORK, INFO ) |
* |
* |
* -- LAPACK driver routine (version 3.6.0) -- |
* -- LAPACK driver routine -- |
* -- 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 |
|
* |
* |
* .. Scalar Arguments .. |
* .. Scalar Arguments .. |
CHARACTER JOBU, JOBVT, RANGE |
CHARACTER JOBU, JOBVT, RANGE |
Line 283
|
Line 286
|
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, ITGKZ, IUTGK, |
$ J, MAXWRK, MINMN, MINWRK, MNTHR |
$ J, MAXWRK, MINMN, MINWRK, MNTHR |
DOUBLE PRECISION ABSTOL, ANRM, BIGNUM, EPS, SMLNUM |
DOUBLE PRECISION ABSTOL, ANRM, BIGNUM, EPS, SMLNUM |
* .. |
* .. |
Line 293
|
Line 296
|
* .. External Subroutines .. |
* .. External Subroutines .. |
EXTERNAL DBDSVDX, DGEBRD, DGELQF, DGEQRF, DLACPY, |
EXTERNAL DBDSVDX, DGEBRD, DGELQF, DGEQRF, DLACPY, |
$ DLASCL, DLASET, DORMBR, DORMLQ, DORMQR, |
$ DLASCL, DLASET, DORMBR, DORMLQ, DORMQR, |
$ DSCAL, XERBLA |
$ DCOPY, XERBLA |
* .. |
* .. |
* .. External Functions .. |
* .. External Functions .. |
LOGICAL LSAME |
LOGICAL LSAME |
INTEGER ILAENV |
INTEGER ILAENV |
DOUBLE PRECISION DLAMCH, DLANGE, DNRM2 |
DOUBLE PRECISION DLAMCH, DLANGE |
EXTERNAL LSAME, ILAENV, DLAMCH, DLANGE, DNRM2 |
EXTERNAL LSAME, ILAENV, DLAMCH, DLANGE |
* .. |
* .. |
* .. Intrinsic Functions .. |
* .. Intrinsic Functions .. |
INTRINSIC MAX, MIN, SQRT |
INTRINSIC MAX, MIN, SQRT |
Line 357
|
Line 360
|
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 380
|
Line 389
|
* |
* |
* Path 1 (M much larger than N) |
* Path 1 (M much larger than N) |
* |
* |
MAXWRK = N*(N*2+16) + |
MAXWRK = N + |
$ N*ILAENV( 1, 'DGEQRF', ' ', M, N, -1, -1 ) |
$ 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 ) ) |
$ 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 |
ELSE |
* |
* |
* Path 2 (M at least N, but not much larger) |
* 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 ) |
$ 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 |
END IF |
ELSE |
ELSE |
MNTHR = ILAENV( 6, 'DGESVD', JOBU // JOBVT, M, N, 0, 0 ) |
MNTHR = ILAENV( 6, 'DGESVD', JOBU // JOBVT, M, N, 0, 0 ) |
Line 399
|
Line 424
|
* |
* |
* Path 1t (N much larger than M) |
* Path 1t (N much larger than M) |
* |
* |
MAXWRK = M*(M*2+16) + |
MAXWRK = M + |
$ M*ILAENV( 1, 'DGELQF', ' ', M, N, -1, -1 ) |
$ 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 ) ) |
$ 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 |
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 ) |
$ 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 |
END IF |
END IF |
END IF |
Line 445
|
Line 486
|
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 489
|
Line 530
|
ITEMP = ITAU + N |
ITEMP = ITAU + N |
CALL DGEQRF( M, N, A, LDA, WORK( ITAU ), WORK( ITEMP ), |
CALL DGEQRF( 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+5*N, prefer N*N+4*N+2*N*NB) |
* (Workspace: need N*N+5*N, prefer N*N+4*N+2*N*NB) |
* |
* |
Line 498
|
Line 539
|
IE = ID + N |
IE = ID + N |
ITAUQ = IE + N |
ITAUQ = IE + N |
ITAUP = ITAUQ + N |
ITAUP = ITAUQ + N |
ITEMP = ITAUP + N |
ITEMP = ITAUP + N |
CALL DLACPY( 'U', N, N, A, LDA, WORK( IQRF ), 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 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 ), |
$ WORK( ITAUQ ), WORK( ITAUP ), WORK( ITEMP ), |
$ LWORK-ITEMP+1, INFO ) |
$ LWORK-ITEMP+1, INFO ) |
* |
* |
* Solve eigenvalue problem TGK*Z=Z*S. |
* 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 |
ITGKZ = ITEMP |
ITEMP = ITGKZ + N*(N*2+1) |
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 ), |
$ VL, VU, ILTGK, IUTGK, NS, S, WORK( ITGKZ ), |
$ N*2, WORK( ITEMP ), IWORK, INFO) |
$ N*2, WORK( ITEMP ), IWORK, INFO) |
* |
* |
Line 522
|
Line 563
|
CALL DCOPY( N, WORK( J ), 1, U( 1,I ), 1 ) |
CALL DCOPY( N, WORK( J ), 1, U( 1,I ), 1 ) |
J = J + N*2 |
J = J + N*2 |
END DO |
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. |
* Call DORMBR to compute QB*UB. |
* (Workspace in WORK( ITEMP ): need N, prefer N*NB) |
* (Workspace in WORK( ITEMP ): need N, prefer N*NB) |
* |
* |
CALL DORMBR( 'Q', 'L', 'N', N, NS, N, WORK( IQRF ), N, |
CALL DORMBR( '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 DORMQR to compute Q*(QB*UB). |
* Call DORMQR to compute Q*(QB*UB). |
* (Workspace in WORK( ITEMP ): need N, prefer N*NB) |
* (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 ), |
$ 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 551
|
Line 592
|
* Call DORMBR to compute VB**T * PB**T |
* Call DORMBR 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 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 ), |
$ WORK( ITAUP ), VT, LDVT, WORK( ITEMP ), |
$ LWORK-ITEMP+1, INFO ) |
$ LWORK-ITEMP+1, INFO ) |
END IF |
END IF |
Line 569
|
Line 610
|
IE = ID + N |
IE = ID + N |
ITAUQ = IE + N |
ITAUQ = IE + N |
ITAUP = ITAUQ + N |
ITAUP = ITAUQ + N |
ITEMP = ITAUP + N |
ITEMP = ITAUP + N |
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 ), |
$ WORK( ITAUQ ), WORK( ITAUP ), WORK( ITEMP ), |
$ LWORK-ITEMP+1, INFO ) |
$ LWORK-ITEMP+1, INFO ) |
* |
* |
* Solve eigenvalue problem TGK*Z=Z*S. |
* 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 |
ITGKZ = ITEMP |
ITEMP = ITGKZ + N*(N*2+1) |
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 ), |
$ VL, VU, ILTGK, IUTGK, NS, S, WORK( ITGKZ ), |
$ N*2, WORK( ITEMP ), IWORK, INFO) |
$ N*2, WORK( ITEMP ), IWORK, INFO) |
* |
* |
Line 591
|
Line 632
|
CALL DCOPY( N, WORK( J ), 1, U( 1,I ), 1 ) |
CALL DCOPY( N, WORK( J ), 1, U( 1,I ), 1 ) |
J = J + N*2 |
J = J + N*2 |
END DO |
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. |
* Call DORMBR to compute QB*UB. |
* (Workspace in WORK( ITEMP ): need N, prefer N*NB) |
* (Workspace in WORK( ITEMP ): need N, prefer N*NB) |
* |
* |
CALL DORMBR( 'Q', 'L', 'N', M, NS, N, A, LDA, |
CALL DORMBR( '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 613
|
Line 654
|
* Call DORMBR to compute VB**T * PB**T |
* Call DORMBR 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 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 ), |
$ 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 626
|
Line 667
|
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 649
|
Line 690
|
ITEMP = ITAUP + M |
ITEMP = ITAUP + M |
CALL DLACPY( 'L', M, M, A, LDA, WORK( ILQF ), 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 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 ), |
$ WORK( ITAUQ ), WORK( ITAUP ), WORK( ITEMP ), |
$ LWORK-ITEMP+1, INFO ) |
$ LWORK-ITEMP+1, INFO ) |
* |
* |
* 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) |
* |
* |
ITGKZ = ITEMP |
ITGKZ = ITEMP |
ITEMP = ITGKZ + M*(M*2+1) |
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 ), |
$ VL, VU, ILTGK, IUTGK, NS, S, WORK( ITGKZ ), |
$ M*2, WORK( ITEMP ), IWORK, INFO) |
$ M*2, WORK( ITEMP ), IWORK, INFO) |
* |
* |
Line 674
|
Line 715
|
* Call DORMBR to compute QB*UB. |
* Call DORMBR to compute QB*UB. |
* (Workspace in WORK( ITEMP ): need M, prefer M*NB) |
* (Workspace in WORK( ITEMP ): need M, prefer M*NB) |
* |
* |
CALL DORMBR( 'Q', 'L', 'N', M, NS, M, WORK( ILQF ), M, |
CALL DORMBR( '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 687
|
Line 728
|
CALL DCOPY( M, WORK( J ), 1, VT( I,1 ), LDVT ) |
CALL DCOPY( M, WORK( J ), 1, VT( I,1 ), LDVT ) |
J = J + M*2 |
J = J + M*2 |
END DO |
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) |
* Call DORMBR 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 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 ), |
$ WORK( ITAUP ), VT, LDVT, WORK( ITEMP ), |
$ LWORK-ITEMP+1, INFO ) |
$ LWORK-ITEMP+1, INFO ) |
* |
* |
* Call DORMLQ to compute ((VB**T)*(PB**T))*Q. |
* Call DORMLQ 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 DORMLQ( 'R', 'N', NS, N, M, A, LDA, |
CALL DORMLQ( '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 4*M+N, prefer 4*M+(M+N)*NB) |
* (Workspace: need 4*M+N, prefer 4*M+(M+N)*NB) |
Line 718
|
Line 759
|
ITAUQ = IE + M |
ITAUQ = IE + M |
ITAUP = ITAUQ + M |
ITAUP = ITAUQ + M |
ITEMP = ITAUP + 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 ), |
$ WORK( ITAUQ ), WORK( ITAUP ), WORK( ITEMP ), |
$ LWORK-ITEMP+1, INFO ) |
$ LWORK-ITEMP+1, INFO ) |
* |
* |
* 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) |
* |
* |
ITGKZ = ITEMP |
ITGKZ = ITEMP |
ITEMP = ITGKZ + M*(M*2+1) |
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 ), |
$ VL, VU, ILTGK, IUTGK, NS, S, WORK( ITGKZ ), |
$ M*2, WORK( ITEMP ), IWORK, INFO) |
$ M*2, WORK( ITEMP ), IWORK, INFO) |
* |
* |
* If needed, compute left singular vectors. |
* If needed, compute left singular vectors. |
* |
* |
IF( WANTU ) THEN |
IF( WANTU ) THEN |
Line 743
|
Line 784
|
* Call DORMBR to compute QB*UB. |
* Call DORMBR to compute QB*UB. |
* (Workspace in WORK( ITEMP ): need M, prefer M*NB) |
* (Workspace in WORK( ITEMP ): need M, prefer M*NB) |
* |
* |
CALL DORMBR( 'Q', 'L', 'N', M, NS, N, A, LDA, |
CALL DORMBR( '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 756
|
Line 797
|
CALL DCOPY( M, WORK( J ), 1, VT( I,1 ), LDVT ) |
CALL DCOPY( M, WORK( J ), 1, VT( I,1 ), LDVT ) |
J = J + M*2 |
J = J + M*2 |
END DO |
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 |
* Call DORMBR 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 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 ), |
$ 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 |
* |
* |