version 1.1, 2015/11/26 11:44:15
|
version 1.7, 2018/05/29 07:17:49
|
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 DBDSVDX + dependencies |
*> Download DBDSVDX + dependencies |
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dbdsvdx.f"> |
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dbdsvdx.f"> |
*> [TGZ]</a> |
*> [TGZ]</a> |
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dbdsvdx.f"> |
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dbdsvdx.f"> |
*> [ZIP]</a> |
*> [ZIP]</a> |
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dbdsvdx.f"> |
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dbdsvdx.f"> |
*> [TXT]</a> |
*> [TXT]</a> |
*> \endhtmlonly |
*> \endhtmlonly |
* |
* |
* Definition: |
* Definition: |
* =========== |
* =========== |
* |
* |
* SUBROUTINE DBDSVDX( UPLO, JOBZ, RANGE, N, D, E, VL, VU, IL, IU, |
* SUBROUTINE DBDSVDX( UPLO, JOBZ, RANGE, N, D, E, VL, VU, IL, IU, |
* $ NS, S, Z, LDZ, WORK, IWORK, INFO ) |
* $ NS, S, Z, LDZ, WORK, IWORK, INFO ) |
* |
* |
* .. Scalar Arguments .. |
* .. Scalar Arguments .. |
Line 28
|
Line 28
|
* .. |
* .. |
* .. Array Arguments .. |
* .. Array Arguments .. |
* INTEGER IWORK( * ) |
* INTEGER IWORK( * ) |
* DOUBLE PRECISION D( * ), E( * ), S( * ), WORK( * ), |
* DOUBLE PRECISION D( * ), E( * ), S( * ), WORK( * ), |
* Z( LDZ, * ) |
* Z( LDZ, * ) |
* .. |
* .. |
* |
* |
*> \par Purpose: |
*> \par Purpose: |
* ============= |
* ============= |
*> |
*> |
*> \verbatim |
*> \verbatim |
*> |
*> |
*> DBDSVDX computes the singular value decomposition (SVD) of a real |
*> DBDSVDX computes the singular value decomposition (SVD) of a real |
*> N-by-N (upper or lower) bidiagonal matrix B, B = U * S * VT, |
*> N-by-N (upper or lower) bidiagonal matrix B, B = U * S * VT, |
*> where S is a diagonal matrix with non-negative diagonal elements |
*> where S is a diagonal matrix with non-negative diagonal elements |
*> (the singular values of B), and U and VT are orthogonal matrices |
*> (the singular values of B), and U and VT are orthogonal matrices |
*> of left and right singular vectors, respectively. |
*> of left and right singular vectors, respectively. |
*> |
*> |
*> Given an upper bidiagonal B with diagonal D = [ d_1 d_2 ... d_N ] |
*> Given an upper bidiagonal B with diagonal D = [ d_1 d_2 ... d_N ] |
*> and superdiagonal E = [ e_1 e_2 ... e_N-1 ], DBDSVDX computes the |
*> and superdiagonal E = [ e_1 e_2 ... e_N-1 ], DBDSVDX computes the |
*> singular value decompositon of B through the eigenvalues and |
*> singular value decompositon of B through the eigenvalues and |
*> eigenvectors of the N*2-by-N*2 tridiagonal matrix |
*> eigenvectors of the N*2-by-N*2 tridiagonal matrix |
*> |
*> |
*> | 0 d_1 | |
*> | 0 d_1 | |
*> | d_1 0 e_1 | |
*> | d_1 0 e_1 | |
*> TGK = | e_1 0 d_2 | |
*> TGK = | e_1 0 d_2 | |
*> | d_2 . . | |
*> | d_2 . . | |
*> | . . . | |
*> | . . . | |
*> |
*> |
*> If (s,u,v) is a singular triplet of B with ||u|| = ||v|| = 1, then |
*> If (s,u,v) is a singular triplet of B with ||u|| = ||v|| = 1, then |
*> (+/-s,q), ||q|| = 1, are eigenpairs of TGK, with q = P * ( u' +/-v' ) / |
*> (+/-s,q), ||q|| = 1, are eigenpairs of TGK, with q = P * ( u' +/-v' ) / |
*> sqrt(2) = ( v_1 u_1 v_2 u_2 ... v_n u_n ) / sqrt(2), and |
*> sqrt(2) = ( v_1 u_1 v_2 u_2 ... v_n u_n ) / sqrt(2), and |
*> P = [ e_{n+1} e_{1} e_{n+2} e_{2} ... ]. |
*> P = [ e_{n+1} e_{1} e_{n+2} e_{2} ... ]. |
*> |
*> |
*> Given a TGK matrix, one can either a) compute -s,-v and change signs |
*> Given a TGK matrix, one can either a) compute -s,-v and change signs |
*> so that the singular values (and corresponding vectors) are already in |
*> so that the singular values (and corresponding vectors) are already in |
*> descending order (as in DGESVD/DGESDD) or b) compute s,v and reorder |
*> descending order (as in DGESVD/DGESDD) or b) compute s,v and reorder |
*> the values (and corresponding vectors). DBDSVDX implements a) by |
*> the values (and corresponding vectors). DBDSVDX implements a) by |
*> calling DSTEVX (bisection plus inverse iteration, to be replaced |
*> calling DSTEVX (bisection plus inverse iteration, to be replaced |
*> with a version of the Multiple Relative Robust Representation |
*> with a version of the Multiple Relative Robust Representation |
*> algorithm. (See P. Willems and B. Lang, A framework for the MR^3 |
*> algorithm. (See P. Willems and B. Lang, A framework for the MR^3 |
*> algorithm: theory and implementation, SIAM J. Sci. Comput., |
*> algorithm: theory and implementation, SIAM J. Sci. Comput., |
*> 35:740-766, 2013.) |
*> 35:740-766, 2013.) |
*> \endverbatim |
*> \endverbatim |
* |
* |
Line 80
|
Line 80
|
*> = 'L': B is lower bidiagonal. |
*> = 'L': B is lower bidiagonal. |
*> \endverbatim |
*> \endverbatim |
*> |
*> |
*> \param[in] JOBXZ |
*> \param[in] JOBZ |
*> \verbatim |
*> \verbatim |
*> JOBZ is CHARACTER*1 |
*> JOBZ is CHARACTER*1 |
*> = 'N': Compute singular values only; |
*> = 'N': Compute singular values only; |
Line 101
|
Line 101
|
*> N is INTEGER |
*> N is INTEGER |
*> The order of the bidiagonal matrix. N >= 0. |
*> The order of the bidiagonal matrix. N >= 0. |
*> \endverbatim |
*> \endverbatim |
*> |
*> |
*> \param[in] D |
*> \param[in] D |
*> \verbatim |
*> \verbatim |
*> D is DOUBLE PRECISION array, dimension (N) |
*> D is DOUBLE PRECISION array, dimension (N) |
*> The n diagonal elements of the bidiagonal matrix B. |
*> The n diagonal elements of the bidiagonal matrix B. |
*> \endverbatim |
*> \endverbatim |
*> |
*> |
*> \param[in] E |
*> \param[in] E |
*> \verbatim |
*> \verbatim |
*> E is DOUBLE PRECISION array, dimension (max(1,N-1)) |
*> E is DOUBLE PRECISION array, dimension (max(1,N-1)) |
Line 117
|
Line 117
|
*> |
*> |
*> \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 132
|
Line 134
|
*> \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 161
|
Line 167
|
*> \verbatim |
*> \verbatim |
*> Z is DOUBLE PRECISION array, dimension (2*N,K) ) |
*> Z is DOUBLE PRECISION array, dimension (2*N,K) ) |
*> If JOBZ = 'V', then if INFO = 0 the first NS columns of Z |
*> If JOBZ = 'V', then if INFO = 0 the first NS columns of Z |
*> contain the singular vectors of the matrix B corresponding to |
*> contain the singular vectors of the matrix B corresponding to |
*> the selected singular values, with U in rows 1 to N and V |
*> the selected singular values, with U in rows 1 to N and V |
*> in rows N+1 to N*2, i.e. |
*> in rows N+1 to N*2, i.e. |
*> Z = [ U ] |
*> Z = [ U ] |
*> [ V ] |
*> [ V ] |
*> If JOBZ = 'N', then Z is not referenced. |
*> If JOBZ = 'N', then Z is not referenced. |
*> Note: The user must ensure that at least K = NS+1 columns are |
*> Note: The user must ensure that at least K = NS+1 columns are |
*> supplied in the array Z; if RANGE = 'V', the exact value of |
*> supplied in the array Z; if RANGE = 'V', the exact value of |
*> NS is not known in advance and an upper bound must be used. |
*> NS is not known in advance and an upper bound must be used. |
*> \endverbatim |
*> \endverbatim |
*> |
*> |
Line 188
|
Line 194
|
*> \verbatim |
*> \verbatim |
*> IWORK is INTEGER array, dimension (12*N) |
*> IWORK is INTEGER array, dimension (12*N) |
*> If JOBZ = 'V', then if INFO = 0, the first NS elements of |
*> If JOBZ = 'V', then if INFO = 0, the first NS elements of |
*> IWORK are zero. If INFO > 0, then IWORK contains the indices |
*> IWORK are zero. If INFO > 0, then IWORK contains the indices |
*> of the eigenvectors that failed to converge in DSTEVX. |
*> of the eigenvectors that failed to converge in DSTEVX. |
|
*> \endverbatim |
*> |
*> |
|
*> \param[out] INFO |
|
*> \verbatim |
*> INFO is INTEGER |
*> INFO is INTEGER |
*> = 0: successful exit |
*> = 0: successful exit |
*> < 0: if INFO = -i, the i-th argument had an illegal value |
*> < 0: if INFO = -i, the i-th argument had an illegal value |
Line 204
|
Line 213
|
* 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 2011 |
*> \date June 2016 |
* |
* |
*> \ingroup doubleOTHEReigen |
*> \ingroup doubleOTHEReigen |
* |
* |
* ===================================================================== |
* ===================================================================== |
SUBROUTINE DBDSVDX( UPLO, JOBZ, RANGE, N, D, E, VL, VU, IL, IU, |
SUBROUTINE DBDSVDX( UPLO, JOBZ, RANGE, N, D, E, VL, VU, IL, IU, |
$ NS, S, Z, LDZ, WORK, IWORK, INFO) |
$ NS, S, Z, LDZ, WORK, IWORK, INFO) |
* |
* |
* -- LAPACK driver routine (version 3.6.0) -- |
* -- 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..-- |
* November 2016 |
* November 2017 |
* |
* |
* .. Scalar Arguments .. |
* .. Scalar Arguments .. |
CHARACTER JOBZ, RANGE, UPLO |
CHARACTER JOBZ, RANGE, UPLO |
INTEGER IL, INFO, IU, LDZ, N, NS |
INTEGER IL, INFO, IU, LDZ, N, NS |
Line 229
|
Line 238
|
* .. |
* .. |
* .. Array Arguments .. |
* .. Array Arguments .. |
INTEGER IWORK( * ) |
INTEGER IWORK( * ) |
DOUBLE PRECISION D( * ), E( * ), S( * ), WORK( * ), |
DOUBLE PRECISION D( * ), E( * ), S( * ), WORK( * ), |
$ Z( LDZ, * ) |
$ Z( LDZ, * ) |
* .. |
* .. |
* |
* |
* ===================================================================== |
* ===================================================================== |
* |
* |
* .. Parameters .. |
* .. Parameters .. |
DOUBLE PRECISION ZERO, ONE, TEN, HNDRD, MEIGTH |
DOUBLE PRECISION ZERO, ONE, TEN, HNDRD, MEIGTH |
PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0, TEN = 10.0D0, |
PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0, TEN = 10.0D0, |
$ HNDRD = 100.0D0, MEIGTH = -0.1250D0 ) |
$ HNDRD = 100.0D0, MEIGTH = -0.1250D0 ) |
DOUBLE PRECISION FUDGE |
DOUBLE PRECISION FUDGE |
PARAMETER ( FUDGE = 2.0D0 ) |
PARAMETER ( FUDGE = 2.0D0 ) |
* .. |
* .. |
* .. Local Scalars .. |
* .. Local Scalars .. |
CHARACTER RNGVX |
CHARACTER RNGVX |
LOGICAL ALLSV, INDSV, LOWER, SPLIT, SVEQ0, VALSV, WANTZ |
LOGICAL ALLSV, INDSV, LOWER, SPLIT, SVEQ0, VALSV, WANTZ |
INTEGER I, ICOLZ, IDBEG, IDEND, IDTGK, IDPTR, IEPTR, |
INTEGER I, ICOLZ, IDBEG, IDEND, IDTGK, IDPTR, IEPTR, |
$ IETGK, IIFAIL, IIWORK, ILTGK, IROWU, IROWV, |
$ IETGK, IIFAIL, IIWORK, ILTGK, IROWU, IROWV, |
$ IROWZ, ISBEG, ISPLT, ITEMP, IUTGK, J, K, |
$ IROWZ, ISBEG, ISPLT, ITEMP, IUTGK, J, K, |
$ NTGK, NRU, NRV, NSL |
$ NTGK, NRU, NRV, NSL |
DOUBLE PRECISION ABSTOL, EPS, EMIN, MU, NRMU, NRMV, ORTOL, SMAX, |
DOUBLE PRECISION ABSTOL, EPS, EMIN, MU, NRMU, NRMV, ORTOL, SMAX, |
$ SMIN, SQRT2, THRESH, TOL, ULP, |
$ SMIN, SQRT2, THRESH, TOL, ULP, |
$ VLTGK, VUTGK, ZJTJI |
$ VLTGK, VUTGK, ZJTJI |
* .. |
* .. |
* .. External Functions .. |
* .. External Functions .. |
Line 260
|
Line 269
|
EXTERNAL IDAMAX, LSAME, DAXPY, DDOT, DLAMCH, DNRM2 |
EXTERNAL IDAMAX, LSAME, DAXPY, DDOT, DLAMCH, DNRM2 |
* .. |
* .. |
* .. External Subroutines .. |
* .. External Subroutines .. |
EXTERNAL DCOPY, DLASET, DSCAL, DSWAP |
EXTERNAL DSTEVX, DCOPY, DLASET, DSCAL, DSWAP, XERBLA |
* .. |
* .. |
* .. Intrinsic Functions .. |
* .. Intrinsic Functions .. |
INTRINSIC ABS, DBLE, SIGN, SQRT |
INTRINSIC ABS, DBLE, SIGN, SQRT |
* .. |
* .. |
* .. Executable Statements .. |
* .. Executable Statements .. |
* |
* |
* Test the input parameters. |
* Test the input parameters. |
* |
* |
Line 312
|
Line 321
|
* |
* |
NS = 0 |
NS = 0 |
IF( N.EQ.0 ) RETURN |
IF( N.EQ.0 ) RETURN |
* |
* |
IF( N.EQ.1 ) THEN |
IF( N.EQ.1 ) THEN |
IF( ALLSV .OR. INDSV ) THEN |
IF( ALLSV .OR. INDSV ) THEN |
NS = 1 |
NS = 1 |
Line 330
|
Line 339
|
RETURN |
RETURN |
END IF |
END IF |
* |
* |
ABSTOL = 2*DLAMCH( 'Safe Minimum' ) |
ABSTOL = 2*DLAMCH( 'Safe Minimum' ) |
ULP = DLAMCH( 'Precision' ) |
ULP = DLAMCH( 'Precision' ) |
EPS = DLAMCH( 'Epsilon' ) |
EPS = DLAMCH( 'Epsilon' ) |
SQRT2 = SQRT( 2.0D0 ) |
SQRT2 = SQRT( 2.0D0 ) |
ORTOL = SQRT( ULP ) |
ORTOL = SQRT( ULP ) |
* |
* |
* Criterion for splitting is taken from DBDSQR when singular |
* Criterion for splitting is taken from DBDSQR when singular |
* values are computed to relative accuracy TOL. (See J. Demmel and |
* values are computed to relative accuracy TOL. (See J. Demmel and |
* W. Kahan, Accurate singular values of bidiagonal matrices, SIAM |
* W. Kahan, Accurate singular values of bidiagonal matrices, SIAM |
* J. Sci. and Stat. Comput., 11:873–912, 1990.) |
* J. Sci. and Stat. Comput., 11:873–912, 1990.) |
* |
* |
TOL = MAX( TEN, MIN( HNDRD, EPS**MEIGTH ) )*EPS |
TOL = MAX( TEN, MIN( HNDRD, EPS**MEIGTH ) )*EPS |
* |
* |
* Compute approximate maximum, minimum singular values. |
* Compute approximate maximum, minimum singular values. |
Line 371
|
Line 380
|
IF( ABS( E( I ) ).LE.THRESH ) E( I ) = ZERO |
IF( ABS( E( I ) ).LE.THRESH ) E( I ) = ZERO |
END DO |
END DO |
IF( ABS( D( N ) ).LE.THRESH ) D( N ) = ZERO |
IF( ABS( D( N ) ).LE.THRESH ) D( N ) = ZERO |
E( N ) = ZERO |
|
* |
* |
* Pointers for arrays used by DSTEVX. |
* Pointers for arrays used by DSTEVX. |
* |
* |
Line 382
|
Line 390
|
IIWORK = IIFAIL + N*2 |
IIWORK = IIFAIL + N*2 |
* |
* |
* Set RNGVX, which corresponds to RANGE for DSTEVX in TGK mode. |
* Set RNGVX, which corresponds to RANGE for DSTEVX in TGK mode. |
* VL,VU or IL,IU are redefined to conform to implementation a) |
* VL,VU or IL,IU are redefined to conform to implementation a) |
* described in the leading comments. |
* described in the leading comments. |
* |
* |
ILTGK = 0 |
ILTGK = 0 |
IUTGK = 0 |
IUTGK = 0 |
VLTGK = ZERO |
VLTGK = ZERO |
VUTGK = ZERO |
VUTGK = ZERO |
* |
* |
IF( ALLSV ) THEN |
IF( ALLSV ) THEN |
* |
* |
* All singular values will be found. We aim at -s (see |
* All singular values will be found. We aim at -s (see |
* leading comments) with RNGVX = 'I'. IL and IU are set |
* leading comments) with RNGVX = 'I'. IL and IU are set |
* later (as ILTGK and IUTGK) according to the dimension |
* later (as ILTGK and IUTGK) according to the dimension |
* of the active submatrix. |
* of the active submatrix. |
* |
* |
RNGVX = 'I' |
RNGVX = 'I' |
CALL DLASET( 'F', N*2, N+1, ZERO, ZERO, Z, LDZ ) |
IF( WANTZ ) CALL DLASET( 'F', N*2, N+1, ZERO, ZERO, Z, LDZ ) |
ELSE IF( VALSV ) THEN |
ELSE IF( VALSV ) THEN |
* |
* |
* Find singular values in a half-open interval. We aim |
* Find singular values in a half-open interval. We aim |
Line 411
|
Line 419
|
WORK( IDTGK:IDTGK+2*N-1 ) = ZERO |
WORK( IDTGK:IDTGK+2*N-1 ) = ZERO |
CALL DCOPY( N, D, 1, WORK( IETGK ), 2 ) |
CALL DCOPY( N, D, 1, WORK( IETGK ), 2 ) |
CALL DCOPY( N-1, E, 1, WORK( IETGK+1 ), 2 ) |
CALL DCOPY( N-1, E, 1, WORK( IETGK+1 ), 2 ) |
CALL DSTEVX( 'N', 'V', N*2, WORK( IDTGK ), WORK( IETGK ), |
CALL DSTEVX( 'N', 'V', N*2, WORK( IDTGK ), WORK( IETGK ), |
$ VLTGK, VUTGK, ILTGK, ILTGK, ABSTOL, NS, S, |
$ VLTGK, VUTGK, ILTGK, ILTGK, ABSTOL, NS, S, |
$ Z, LDZ, WORK( ITEMP ), IWORK( IIWORK ), |
$ Z, LDZ, WORK( ITEMP ), IWORK( IIWORK ), |
$ IWORK( IIFAIL ), INFO ) |
$ IWORK( IIFAIL ), INFO ) |
IF( NS.EQ.0 ) THEN |
IF( NS.EQ.0 ) THEN |
RETURN |
RETURN |
ELSE |
ELSE |
CALL DLASET( 'F', N*2, NS, ZERO, ZERO, Z, LDZ ) |
IF( WANTZ ) CALL DLASET( 'F', N*2, NS, ZERO, ZERO, Z, LDZ ) |
END IF |
END IF |
ELSE IF( INDSV ) THEN |
ELSE IF( INDSV ) THEN |
* |
* |
* Find the IL-th through the IU-th singular values. We aim |
* Find the IL-th through the IU-th singular values. We aim |
* at -s (see leading comments) and indices are mapped into |
* at -s (see leading comments) and indices are mapped into |
* values, therefore mimicking DSTEBZ, where |
* values, therefore mimicking DSTEBZ, where |
* |
* |
* GL = GL - FUDGE*TNORM*ULP*N - FUDGE*TWO*PIVMIN |
* GL = GL - FUDGE*TNORM*ULP*N - FUDGE*TWO*PIVMIN |
* GU = GU + FUDGE*TNORM*ULP*N + FUDGE*PIVMIN |
* GU = GU + FUDGE*TNORM*ULP*N + FUDGE*PIVMIN |
* |
* |
ILTGK = IL |
ILTGK = IL |
IUTGK = IU |
IUTGK = IU |
RNGVX = 'V' |
RNGVX = 'V' |
WORK( IDTGK:IDTGK+2*N-1 ) = ZERO |
WORK( IDTGK:IDTGK+2*N-1 ) = ZERO |
CALL DCOPY( N, D, 1, WORK( IETGK ), 2 ) |
CALL DCOPY( N, D, 1, WORK( IETGK ), 2 ) |
CALL DCOPY( N-1, E, 1, WORK( IETGK+1 ), 2 ) |
CALL DCOPY( N-1, E, 1, WORK( IETGK+1 ), 2 ) |
CALL DSTEVX( 'N', 'I', N*2, WORK( IDTGK ), WORK( IETGK ), |
CALL DSTEVX( 'N', 'I', N*2, WORK( IDTGK ), WORK( IETGK ), |
$ VLTGK, VLTGK, ILTGK, ILTGK, ABSTOL, NS, S, |
$ VLTGK, VLTGK, ILTGK, ILTGK, ABSTOL, NS, S, |
$ Z, LDZ, WORK( ITEMP ), IWORK( IIWORK ), |
$ Z, LDZ, WORK( ITEMP ), IWORK( IIWORK ), |
$ IWORK( IIFAIL ), INFO ) |
$ IWORK( IIFAIL ), INFO ) |
Line 443
|
Line 451
|
WORK( IDTGK:IDTGK+2*N-1 ) = ZERO |
WORK( IDTGK:IDTGK+2*N-1 ) = ZERO |
CALL DCOPY( N, D, 1, WORK( IETGK ), 2 ) |
CALL DCOPY( N, D, 1, WORK( IETGK ), 2 ) |
CALL DCOPY( N-1, E, 1, WORK( IETGK+1 ), 2 ) |
CALL DCOPY( N-1, E, 1, WORK( IETGK+1 ), 2 ) |
CALL DSTEVX( 'N', 'I', N*2, WORK( IDTGK ), WORK( IETGK ), |
CALL DSTEVX( 'N', 'I', N*2, WORK( IDTGK ), WORK( IETGK ), |
$ VUTGK, VUTGK, IUTGK, IUTGK, ABSTOL, NS, S, |
$ VUTGK, VUTGK, IUTGK, IUTGK, ABSTOL, NS, S, |
$ Z, LDZ, WORK( ITEMP ), IWORK( IIWORK ), |
$ Z, LDZ, WORK( ITEMP ), IWORK( IIWORK ), |
$ IWORK( IIFAIL ), INFO ) |
$ IWORK( IIFAIL ), INFO ) |
Line 451
|
Line 459
|
VUTGK = MIN( VUTGK, ZERO ) |
VUTGK = MIN( VUTGK, ZERO ) |
* |
* |
* If VLTGK=VUTGK, DSTEVX returns an error message, |
* If VLTGK=VUTGK, DSTEVX returns an error message, |
* so if needed we change VUTGK slightly. |
* so if needed we change VUTGK slightly. |
* |
* |
IF( VLTGK.EQ.VUTGK ) VLTGK = VLTGK - TOL |
IF( VLTGK.EQ.VUTGK ) VLTGK = VLTGK - TOL |
* |
* |
CALL DLASET( 'F', N*2, IU-IL+1, ZERO, ZERO, Z, LDZ ) |
IF( WANTZ ) CALL DLASET( 'F', N*2, IU-IL+1, ZERO, ZERO, Z, LDZ) |
END IF |
END IF |
* |
* |
* Initialize variables and pointers for S, Z, and WORK. |
* Initialize variables and pointers for S, Z, and WORK. |
* |
* |
Line 475
|
Line 483
|
IROWU = 2 |
IROWU = 2 |
IROWV = 1 |
IROWV = 1 |
SPLIT = .FALSE. |
SPLIT = .FALSE. |
SVEQ0 = .FALSE. |
SVEQ0 = .FALSE. |
* |
* |
* Form the tridiagonal TGK matrix. |
* Form the tridiagonal TGK matrix. |
* |
* |
Line 486
|
Line 494
|
CALL DCOPY( N-1, E, 1, WORK( IETGK+1 ), 2 ) |
CALL DCOPY( N-1, E, 1, WORK( IETGK+1 ), 2 ) |
* |
* |
* |
* |
* Check for splits in two levels, outer level |
* Check for splits in two levels, outer level |
* in E and inner level in D. |
* in E and inner level in D. |
* |
* |
DO IEPTR = 2, N*2, 2 |
DO IEPTR = 2, N*2, 2 |
IF( WORK( IETGK+IEPTR-1 ).EQ.ZERO ) THEN |
IF( WORK( IETGK+IEPTR-1 ).EQ.ZERO ) THEN |
* |
* |
* Split in E (this piece of B is square) or bottom |
* Split in E (this piece of B is square) or bottom |
* of the (input bidiagonal) matrix. |
* of the (input bidiagonal) matrix. |
* |
* |
ISPLT = IDBEG |
ISPLT = IDBEG |
IDEND = IEPTR - 1 |
IDEND = IEPTR - 1 |
DO IDPTR = IDBEG, IDEND, 2 |
DO IDPTR = IDBEG, IDEND, 2 |
Line 511
|
Line 519
|
IF( IDBEG.EQ.IDEND) THEN |
IF( IDBEG.EQ.IDEND) THEN |
NRU = 1 |
NRU = 1 |
NRV = 1 |
NRV = 1 |
END IF |
END IF |
ELSE IF( IDPTR.EQ.IDEND ) THEN |
ELSE IF( IDPTR.EQ.IDEND ) THEN |
* |
* |
* D=0 at the bottom. |
* D=0 at the bottom. |
* |
* |
SVEQ0 = .TRUE. |
SVEQ0 = .TRUE. |
NRU = (IDEND-ISPLT)/2 + 1 |
NRU = (IDEND-ISPLT)/2 + 1 |
NRV = NRU |
NRV = NRU |
IF( ISPLT.NE.IDBEG ) THEN |
IF( ISPLT.NE.IDBEG ) THEN |
NRU = NRU + 1 |
NRU = NRU + 1 |
END IF |
END IF |
ELSE |
ELSE |
IF( ISPLT.EQ.IDBEG ) THEN |
IF( ISPLT.EQ.IDBEG ) THEN |
* |
* |
* Split: top rectangular submatrix. |
* Split: top rectangular submatrix. |
* |
* |
NRU = (IDPTR-IDBEG)/2 |
NRU = (IDPTR-IDBEG)/2 |
NRV = NRU + 1 |
NRV = NRU + 1 |
ELSE |
ELSE |
Line 534
|
Line 542
|
* Split: middle square submatrix. |
* Split: middle square submatrix. |
* |
* |
NRU = (IDPTR-ISPLT)/2 + 1 |
NRU = (IDPTR-ISPLT)/2 + 1 |
NRV = NRU |
NRV = NRU |
END IF |
END IF |
END IF |
END IF |
ELSE IF( IDPTR.EQ.IDEND ) THEN |
ELSE IF( IDPTR.EQ.IDEND ) THEN |
Line 552
|
Line 560
|
* Split: bottom rectangular submatrix. |
* Split: bottom rectangular submatrix. |
* |
* |
NRV = (IDEND-ISPLT)/2 + 1 |
NRV = (IDEND-ISPLT)/2 + 1 |
NRU = NRV + 1 |
NRU = NRV + 1 |
END IF |
END IF |
END IF |
END IF |
* |
* |
Line 560
|
Line 568
|
* |
* |
IF( NTGK.GT.0 ) THEN |
IF( NTGK.GT.0 ) THEN |
* |
* |
* Compute eigenvalues/vectors of the active |
* Compute eigenvalues/vectors of the active |
* submatrix according to RANGE: |
* submatrix according to RANGE: |
* if RANGE='A' (ALLSV) then RNGVX = 'I' |
* if RANGE='A' (ALLSV) then RNGVX = 'I' |
* if RANGE='V' (VALSV) then RNGVX = 'V' |
* if RANGE='V' (VALSV) then RNGVX = 'V' |
* if RANGE='I' (INDSV) then RNGVX = 'V' |
* if RANGE='I' (INDSV) then RNGVX = 'V' |
* |
* |
ILTGK = 1 |
ILTGK = 1 |
IUTGK = NTGK / 2 |
IUTGK = NTGK / 2 |
IF( ALLSV .OR. VUTGK.EQ.ZERO ) THEN |
IF( ALLSV .OR. VUTGK.EQ.ZERO ) THEN |
IF( SVEQ0 .OR. |
IF( SVEQ0 .OR. |
$ SMIN.LT.EPS .OR. |
$ SMIN.LT.EPS .OR. |
$ MOD(NTGK,2).GT.0 ) THEN |
$ MOD(NTGK,2).GT.0 ) THEN |
* Special case: eigenvalue equal to zero or very |
* Special case: eigenvalue equal to zero or very |
* small, additional eigenvector is needed. |
* small, additional eigenvector is needed. |
IUTGK = IUTGK + 1 |
IUTGK = IUTGK + 1 |
END IF |
END IF |
END IF |
END IF |
* |
* |
* Workspace needed by DSTEVX: |
* Workspace needed by DSTEVX: |
* WORK( ITEMP: ): 2*5*NTGK |
* WORK( ITEMP: ): 2*5*NTGK |
* IWORK( 1: ): 2*6*NTGK |
* IWORK( 1: ): 2*6*NTGK |
* |
* |
CALL DSTEVX( JOBZ, RNGVX, NTGK, WORK( IDTGK+ISPLT-1 ), |
CALL DSTEVX( JOBZ, RNGVX, NTGK, WORK( IDTGK+ISPLT-1 ), |
$ WORK( IETGK+ISPLT-1 ), VLTGK, VUTGK, |
$ WORK( IETGK+ISPLT-1 ), VLTGK, VUTGK, |
$ ILTGK, IUTGK, ABSTOL, NSL, S( ISBEG ), |
$ ILTGK, IUTGK, ABSTOL, NSL, S( ISBEG ), |
$ Z( IROWZ,ICOLZ ), LDZ, WORK( ITEMP ), |
$ Z( IROWZ,ICOLZ ), LDZ, WORK( ITEMP ), |
$ IWORK( IIWORK ), IWORK( IIFAIL ), |
$ IWORK( IIWORK ), IWORK( IIFAIL ), |
$ INFO ) |
$ INFO ) |
IF( INFO.NE.0 ) THEN |
IF( INFO.NE.0 ) THEN |
Line 593
|
Line 601
|
RETURN |
RETURN |
END IF |
END IF |
EMIN = ABS( MAXVAL( S( ISBEG:ISBEG+NSL-1 ) ) ) |
EMIN = ABS( MAXVAL( S( ISBEG:ISBEG+NSL-1 ) ) ) |
* |
* |
IF( NSL.GT.0 .AND. WANTZ ) THEN |
IF( NSL.GT.0 .AND. WANTZ ) THEN |
* |
* |
* Normalize u=Z([2,4,...],:) and v=Z([1,3,...],:), |
* Normalize u=Z([2,4,...],:) and v=Z([1,3,...],:), |
Line 607
|
Line 615
|
IF( NSL.GT.1 .AND. |
IF( NSL.GT.1 .AND. |
$ VUTGK.EQ.ZERO .AND. |
$ VUTGK.EQ.ZERO .AND. |
$ MOD(NTGK,2).EQ.0 .AND. |
$ MOD(NTGK,2).EQ.0 .AND. |
$ EMIN.EQ.0 .AND. .NOT.SPLIT ) THEN |
$ EMIN.EQ.0 .AND. .NOT.SPLIT ) THEN |
* |
* |
* D=0 at the top or bottom of the active submatrix: |
* D=0 at the top or bottom of the active submatrix: |
* one eigenvalue is equal to zero; concatenate the |
* one eigenvalue is equal to zero; concatenate the |
* eigenvectors corresponding to the two smallest |
* eigenvectors corresponding to the two smallest |
* eigenvalues. |
* eigenvalues. |
* |
* |
Z( IROWZ:IROWZ+NTGK-1,ICOLZ+NSL-2 ) = |
Z( IROWZ:IROWZ+NTGK-1,ICOLZ+NSL-2 ) = |
$ Z( IROWZ:IROWZ+NTGK-1,ICOLZ+NSL-2 ) + |
$ Z( IROWZ:IROWZ+NTGK-1,ICOLZ+NSL-2 ) + |
$ Z( IROWZ:IROWZ+NTGK-1,ICOLZ+NSL-1 ) |
$ Z( IROWZ:IROWZ+NTGK-1,ICOLZ+NSL-1 ) |
Z( IROWZ:IROWZ+NTGK-1,ICOLZ+NSL-1 ) = |
Z( IROWZ:IROWZ+NTGK-1,ICOLZ+NSL-1 ) = |
$ ZERO |
$ ZERO |
* IF( IUTGK*2.GT.NTGK ) THEN |
* IF( IUTGK*2.GT.NTGK ) THEN |
* Eigenvalue equal to zero or very small. |
* Eigenvalue equal to zero or very small. |
* NSL = NSL - 1 |
* NSL = NSL - 1 |
* END IF |
* END IF |
END IF |
END IF |
* |
* |
DO I = 0, MIN( NSL-1, NRU-1 ) |
DO I = 0, MIN( NSL-1, NRU-1 ) |
Line 631
|
Line 639
|
INFO = N*2 + 1 |
INFO = N*2 + 1 |
RETURN |
RETURN |
END IF |
END IF |
CALL DSCAL( NRU, ONE/NRMU, |
CALL DSCAL( NRU, ONE/NRMU, |
$ Z( IROWU,ICOLZ+I ), 2 ) |
$ Z( IROWU,ICOLZ+I ), 2 ) |
IF( NRMU.NE.ONE .AND. |
IF( NRMU.NE.ONE .AND. |
$ ABS( NRMU-ORTOL )*SQRT2.GT.ONE ) |
$ ABS( NRMU-ORTOL )*SQRT2.GT.ONE ) |
$ THEN |
$ THEN |
DO J = 0, I-1 |
DO J = 0, I-1 |
ZJTJI = -DDOT( NRU, Z( IROWU, ICOLZ+J ), |
ZJTJI = -DDOT( NRU, Z( IROWU, ICOLZ+J ), |
$ 2, Z( IROWU, ICOLZ+I ), 2 ) |
$ 2, Z( IROWU, ICOLZ+I ), 2 ) |
CALL DAXPY( NRU, ZJTJI, |
CALL DAXPY( NRU, ZJTJI, |
$ Z( IROWU, ICOLZ+J ), 2, |
$ Z( IROWU, ICOLZ+J ), 2, |
$ Z( IROWU, ICOLZ+I ), 2 ) |
$ Z( IROWU, ICOLZ+I ), 2 ) |
END DO |
END DO |
NRMU = DNRM2( NRU, Z( IROWU, ICOLZ+I ), 2 ) |
NRMU = DNRM2( NRU, Z( IROWU, ICOLZ+I ), 2 ) |
CALL DSCAL( NRU, ONE/NRMU, |
CALL DSCAL( NRU, ONE/NRMU, |
$ Z( IROWU,ICOLZ+I ), 2 ) |
$ Z( IROWU,ICOLZ+I ), 2 ) |
END IF |
END IF |
END DO |
END DO |
Line 654
|
Line 662
|
INFO = N*2 + 1 |
INFO = N*2 + 1 |
RETURN |
RETURN |
END IF |
END IF |
CALL DSCAL( NRV, -ONE/NRMV, |
CALL DSCAL( NRV, -ONE/NRMV, |
$ Z( IROWV,ICOLZ+I ), 2 ) |
$ Z( IROWV,ICOLZ+I ), 2 ) |
IF( NRMV.NE.ONE .AND. |
IF( NRMV.NE.ONE .AND. |
$ ABS( NRMV-ORTOL )*SQRT2.GT.ONE ) |
$ ABS( NRMV-ORTOL )*SQRT2.GT.ONE ) |
Line 662
|
Line 670
|
DO J = 0, I-1 |
DO J = 0, I-1 |
ZJTJI = -DDOT( NRV, Z( IROWV, ICOLZ+J ), |
ZJTJI = -DDOT( NRV, Z( IROWV, ICOLZ+J ), |
$ 2, Z( IROWV, ICOLZ+I ), 2 ) |
$ 2, Z( IROWV, ICOLZ+I ), 2 ) |
CALL DAXPY( NRU, ZJTJI, |
CALL DAXPY( NRU, ZJTJI, |
$ Z( IROWV, ICOLZ+J ), 2, |
$ Z( IROWV, ICOLZ+J ), 2, |
$ Z( IROWV, ICOLZ+I ), 2 ) |
$ Z( IROWV, ICOLZ+I ), 2 ) |
END DO |
END DO |
NRMV = DNRM2( NRV, Z( IROWV, ICOLZ+I ), 2 ) |
NRMV = DNRM2( NRV, Z( IROWV, ICOLZ+I ), 2 ) |
CALL DSCAL( NRV, ONE/NRMV, |
CALL DSCAL( NRV, ONE/NRMV, |
$ Z( IROWV,ICOLZ+I ), 2 ) |
$ Z( IROWV,ICOLZ+I ), 2 ) |
END IF |
END IF |
END DO |
END DO |
Line 676
|
Line 684
|
$ MOD(NTGK,2).GT.0 ) THEN |
$ MOD(NTGK,2).GT.0 ) THEN |
* |
* |
* D=0 in the middle of the active submatrix (one |
* D=0 in the middle of the active submatrix (one |
* eigenvalue is equal to zero): save the corresponding |
* eigenvalue is equal to zero): save the corresponding |
* eigenvector for later use (when bottom of the |
* eigenvector for later use (when bottom of the |
* active submatrix is reached). |
* active submatrix is reached). |
* |
* |
SPLIT = .TRUE. |
SPLIT = .TRUE. |
Z( IROWZ:IROWZ+NTGK-1,N+1 ) = |
Z( IROWZ:IROWZ+NTGK-1,N+1 ) = |
$ Z( IROWZ:IROWZ+NTGK-1,NS+NSL ) |
$ Z( IROWZ:IROWZ+NTGK-1,NS+NSL ) |
Z( IROWZ:IROWZ+NTGK-1,NS+NSL ) = |
Z( IROWZ:IROWZ+NTGK-1,NS+NSL ) = |
$ ZERO |
$ ZERO |
END IF |
END IF |
END IF !** WANTZ **! |
END IF !** WANTZ **! |
* |
* |
NSL = MIN( NSL, NRU ) |
NSL = MIN( NSL, NRU ) |
SVEQ0 = .FALSE. |
SVEQ0 = .FALSE. |
* |
* |
Line 698
|
Line 706
|
END DO |
END DO |
* |
* |
* Update pointers for TGK, S and Z. |
* Update pointers for TGK, S and Z. |
* |
* |
ISBEG = ISBEG + NSL |
ISBEG = ISBEG + NSL |
IROWZ = IROWZ + NTGK |
IROWZ = IROWZ + NTGK |
ICOLZ = ICOLZ + NSL |
ICOLZ = ICOLZ + NSL |
IROWU = IROWZ |
IROWU = IROWZ |
IROWV = IROWZ + 1 |
IROWV = IROWZ + 1 |
ISPLT = IDPTR + 1 |
ISPLT = IDPTR + 1 |
NS = NS + NSL |
NS = NS + NSL |
NRU = 0 |
NRU = 0 |
NRV = 0 |
NRV = 0 |
END IF !** NTGK.GT.0 **! |
END IF !** NTGK.GT.0 **! |
IF( IROWZ.LT.N*2 ) Z( 1:IROWZ-1, ICOLZ ) = ZERO |
IF( IROWZ.LT.N*2 .AND. WANTZ ) THEN |
|
Z( 1:IROWZ-1, ICOLZ ) = ZERO |
|
END IF |
END DO !** IDPTR loop **! |
END DO !** IDPTR loop **! |
IF( SPLIT ) THEN |
IF( SPLIT .AND. WANTZ ) THEN |
* |
* |
* Bring back eigenvector corresponding |
* Bring back eigenvector corresponding |
* to eigenvalue equal to zero. |
* to eigenvalue equal to zero. |
* |
* |
Z( IDBEG:IDEND-NTGK+1,ISBEG-1 ) = |
Z( IDBEG:IDEND-NTGK+1,ISBEG-1 ) = |
$ Z( IDBEG:IDEND-NTGK+1,ISBEG-1 ) + |
$ Z( IDBEG:IDEND-NTGK+1,ISBEG-1 ) + |
$ Z( IDBEG:IDEND-NTGK+1,N+1 ) |
$ Z( IDBEG:IDEND-NTGK+1,N+1 ) |
Z( IDBEG:IDEND-NTGK+1,N+1 ) = 0 |
Z( IDBEG:IDEND-NTGK+1,N+1 ) = 0 |
Line 725
|
Line 735
|
IROWU = IROWU + 1 |
IROWU = IROWU + 1 |
IDBEG = IEPTR + 1 |
IDBEG = IEPTR + 1 |
SVEQ0 = .FALSE. |
SVEQ0 = .FALSE. |
SPLIT = .FALSE. |
SPLIT = .FALSE. |
END IF !** Check for split in E **! |
END IF !** Check for split in E **! |
END DO !** IEPTR loop **! |
END DO !** IEPTR loop **! |
* |
* |
Line 744
|
Line 754
|
IF( K.NE.NS+1-I ) THEN |
IF( K.NE.NS+1-I ) THEN |
S( K ) = S( NS+1-I ) |
S( K ) = S( NS+1-I ) |
S( NS+1-I ) = SMIN |
S( NS+1-I ) = SMIN |
CALL DSWAP( N*2, Z( 1,K ), 1, Z( 1,NS+1-I ), 1 ) |
IF( WANTZ ) CALL DSWAP( N*2, Z( 1,K ), 1, Z( 1,NS+1-I ), 1 ) |
END IF |
END IF |
END DO |
END DO |
* |
* |
* If RANGE=I, check for singular values/vectors to be discarded. |
* If RANGE=I, check for singular values/vectors to be discarded. |
* |
* |
IF( INDSV ) THEN |
IF( INDSV ) THEN |
K = IU - IL + 1 |
K = IU - IL + 1 |
IF( K.LT.NS ) THEN |
IF( K.LT.NS ) THEN |
S( K+1:NS ) = ZERO |
S( K+1:NS ) = ZERO |
Z( 1:N*2,K+1:NS ) = ZERO |
IF( WANTZ ) Z( 1:N*2,K+1:NS ) = ZERO |
NS = K |
NS = K |
END IF |
END IF |
END IF |
END IF |
* |
* |
* Reorder Z: U = Z( 1:N,1:NS ), V = Z( N+1:N*2,1:NS ). |
* Reorder Z: U = Z( 1:N,1:NS ), V = Z( N+1:N*2,1:NS ). |
* If B is a lower diagonal, swap U and V. |
* If B is a lower diagonal, swap U and V. |
* |
* |
|
IF( WANTZ ) THEN |
DO I = 1, NS |
DO I = 1, NS |
CALL DCOPY( N*2, Z( 1,I ), 1, WORK, 1 ) |
CALL DCOPY( N*2, Z( 1,I ), 1, WORK, 1 ) |
IF( LOWER ) THEN |
IF( LOWER ) THEN |
Line 772
|
Line 783
|
CALL DCOPY( N, WORK( 1 ), 2, Z( N+1,I ), 1 ) |
CALL DCOPY( N, WORK( 1 ), 2, Z( N+1,I ), 1 ) |
END IF |
END IF |
END DO |
END DO |
|
END IF |
* |
* |
RETURN |
RETURN |
* |
* |