version 1.9, 2011/11/21 22:19:57
|
version 1.21, 2018/05/29 07:18:35
|
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 ZSTEMR + dependencies |
*> Download ZSTEMR + dependencies |
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zstemr.f"> |
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zstemr.f"> |
*> [TGZ]</a> |
*> [TGZ]</a> |
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zstemr.f"> |
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zstemr.f"> |
*> [ZIP]</a> |
*> [ZIP]</a> |
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zstemr.f"> |
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zstemr.f"> |
*> [TXT]</a> |
*> [TXT]</a> |
*> \endhtmlonly |
*> \endhtmlonly |
* |
* |
* Definition: |
* Definition: |
* =========== |
* =========== |
Line 21
|
Line 21
|
* SUBROUTINE ZSTEMR( JOBZ, RANGE, N, D, E, VL, VU, IL, IU, |
* SUBROUTINE ZSTEMR( JOBZ, RANGE, N, D, E, VL, VU, IL, IU, |
* M, W, Z, LDZ, NZC, ISUPPZ, TRYRAC, WORK, LWORK, |
* M, W, Z, LDZ, NZC, ISUPPZ, TRYRAC, WORK, LWORK, |
* IWORK, LIWORK, INFO ) |
* IWORK, LIWORK, INFO ) |
* |
* |
* .. Scalar Arguments .. |
* .. Scalar Arguments .. |
* CHARACTER JOBZ, RANGE |
* CHARACTER JOBZ, RANGE |
* LOGICAL TRYRAC |
* LOGICAL TRYRAC |
Line 33
|
Line 33
|
* DOUBLE PRECISION D( * ), E( * ), W( * ), WORK( * ) |
* DOUBLE PRECISION D( * ), E( * ), W( * ), WORK( * ) |
* COMPLEX*16 Z( LDZ, * ) |
* COMPLEX*16 Z( LDZ, * ) |
* .. |
* .. |
* |
* |
* |
* |
*> \par Purpose: |
*> \par Purpose: |
* ============= |
* ============= |
Line 153
|
Line 153
|
*> \param[in] VL |
*> \param[in] VL |
*> \verbatim |
*> \verbatim |
*> VL is DOUBLE PRECISION |
*> VL is DOUBLE PRECISION |
|
*> |
|
*> If RANGE='V', the lower bound of the interval to |
|
*> be searched for eigenvalues. VL < VU. |
|
*> 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 eigenvalues. VL < VU. |
*> be searched for eigenvalues. VL < VU. |
*> Not referenced if RANGE = 'A' or 'I'. |
*> Not referenced if RANGE = 'A' or 'I'. |
*> \endverbatim |
*> \endverbatim |
Line 167
|
Line 171
|
*> \param[in] IL |
*> \param[in] IL |
*> \verbatim |
*> \verbatim |
*> IL is INTEGER |
*> IL is INTEGER |
|
*> |
|
*> If RANGE='I', the index of the |
|
*> smallest eigenvalue to be returned. |
|
*> 1 <= IL <= IU <= N, if 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 eigenvalues to be returned. |
*> largest eigenvalue to be returned. |
*> 1 <= IL <= IU <= N, if N > 0. |
*> 1 <= IL <= IU <= N, if N > 0. |
*> Not referenced if RANGE = 'A' or 'V'. |
*> Not referenced if RANGE = 'A' or 'V'. |
*> \endverbatim |
*> \endverbatim |
Line 230
|
Line 239
|
*> |
*> |
*> \param[out] ISUPPZ |
*> \param[out] ISUPPZ |
*> \verbatim |
*> \verbatim |
*> ISUPPZ is INTEGER ARRAY, dimension ( 2*max(1,M) ) |
*> ISUPPZ is INTEGER array, dimension ( 2*max(1,M) ) |
*> The support of the eigenvectors in Z, i.e., the indices |
*> The support of the eigenvectors in Z, i.e., the indices |
*> indicating the nonzero elements in Z. The i-th computed eigenvector |
*> indicating the nonzero elements in Z. The i-th computed eigenvector |
*> is nonzero only in elements ISUPPZ( 2*i-1 ) through |
*> is nonzero only in elements ISUPPZ( 2*i-1 ) through |
Line 306
|
Line 315
|
* 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 complex16OTHERcomputational |
*> \ingroup complex16OTHERcomputational |
* |
* |
Line 329
|
Line 338
|
$ M, W, Z, LDZ, NZC, ISUPPZ, TRYRAC, WORK, LWORK, |
$ M, W, Z, LDZ, NZC, ISUPPZ, TRYRAC, WORK, LWORK, |
$ IWORK, LIWORK, INFO ) |
$ IWORK, LIWORK, INFO ) |
* |
* |
* -- LAPACK computational routine (version 3.4.0) -- |
* -- LAPACK computational routine (version 3.7.1) -- |
* -- LAPACK is a software package provided by Univ. of Tennessee, -- |
* -- LAPACK is a software package provided by Univ. of Tennessee, -- |
* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- |
* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- |
* November 2011 |
* June 2016 |
* |
* |
* .. Scalar Arguments .. |
* .. Scalar Arguments .. |
CHARACTER JOBZ, RANGE |
CHARACTER JOBZ, RANGE |
Line 408
|
Line 417
|
WU = ZERO |
WU = ZERO |
IIL = 0 |
IIL = 0 |
IIU = 0 |
IIU = 0 |
|
NSPLIT = 0 |
|
|
IF( VALEIG ) THEN |
IF( VALEIG ) THEN |
* We do not reference VL, VU in the cases RANGE = 'I','A' |
* We do not reference VL, VU in the cases RANGE = 'I','A' |
Line 525
|
Line 535
|
IF (SN.NE.ZERO) THEN |
IF (SN.NE.ZERO) THEN |
IF (CS.NE.ZERO) THEN |
IF (CS.NE.ZERO) THEN |
ISUPPZ(2*M-1) = 1 |
ISUPPZ(2*M-1) = 1 |
ISUPPZ(2*M-1) = 2 |
ISUPPZ(2*M) = 2 |
ELSE |
ELSE |
ISUPPZ(2*M-1) = 1 |
ISUPPZ(2*M-1) = 1 |
ISUPPZ(2*M-1) = 1 |
ISUPPZ(2*M) = 1 |
END IF |
END IF |
ELSE |
ELSE |
ISUPPZ(2*M-1) = 2 |
ISUPPZ(2*M-1) = 2 |
Line 549
|
Line 559
|
IF (SN.NE.ZERO) THEN |
IF (SN.NE.ZERO) THEN |
IF (CS.NE.ZERO) THEN |
IF (CS.NE.ZERO) THEN |
ISUPPZ(2*M-1) = 1 |
ISUPPZ(2*M-1) = 1 |
ISUPPZ(2*M-1) = 2 |
ISUPPZ(2*M) = 2 |
ELSE |
ELSE |
ISUPPZ(2*M-1) = 1 |
ISUPPZ(2*M-1) = 1 |
ISUPPZ(2*M-1) = 1 |
ISUPPZ(2*M) = 1 |
END IF |
END IF |
ELSE |
ELSE |
ISUPPZ(2*M-1) = 2 |
ISUPPZ(2*M-1) = 2 |
Line 560
|
Line 570
|
END IF |
END IF |
ENDIF |
ENDIF |
ENDIF |
ENDIF |
RETURN |
ELSE |
END IF |
|
|
|
* Continue with general N |
* Continue with general N |
|
|
INDGRS = 1 |
INDGRS = 1 |
INDERR = 2*N + 1 |
INDERR = 2*N + 1 |
INDGP = 3*N + 1 |
INDGP = 3*N + 1 |
INDD = 4*N + 1 |
INDD = 4*N + 1 |
INDE2 = 5*N + 1 |
INDE2 = 5*N + 1 |
INDWRK = 6*N + 1 |
INDWRK = 6*N + 1 |
* |
* |
IINSPL = 1 |
IINSPL = 1 |
IINDBL = N + 1 |
IINDBL = N + 1 |
IINDW = 2*N + 1 |
IINDW = 2*N + 1 |
IINDWK = 3*N + 1 |
IINDWK = 3*N + 1 |
* |
* |
* Scale matrix to allowable range, if necessary. |
* Scale matrix to allowable range, if necessary. |
* The allowable range is related to the PIVMIN parameter; see the |
* The allowable range is related to the PIVMIN parameter; see the |
* comments in DLARRD. The preference for scaling small values |
* comments in DLARRD. The preference for scaling small values |
* up is heuristic; we expect users' matrices not to be close to the |
* up is heuristic; we expect users' matrices not to be close to the |
* RMAX threshold. |
* RMAX threshold. |
* |
* |
SCALE = ONE |
SCALE = ONE |
TNRM = DLANST( 'M', N, D, E ) |
TNRM = DLANST( 'M', N, D, E ) |
IF( TNRM.GT.ZERO .AND. TNRM.LT.RMIN ) THEN |
IF( TNRM.GT.ZERO .AND. TNRM.LT.RMIN ) THEN |
SCALE = RMIN / TNRM |
SCALE = RMIN / TNRM |
ELSE IF( TNRM.GT.RMAX ) THEN |
ELSE IF( TNRM.GT.RMAX ) THEN |
SCALE = RMAX / TNRM |
SCALE = RMAX / TNRM |
END IF |
END IF |
IF( SCALE.NE.ONE ) THEN |
IF( SCALE.NE.ONE ) THEN |
CALL DSCAL( N, SCALE, D, 1 ) |
CALL DSCAL( N, SCALE, D, 1 ) |
CALL DSCAL( N-1, SCALE, E, 1 ) |
CALL DSCAL( N-1, SCALE, E, 1 ) |
TNRM = TNRM*SCALE |
TNRM = TNRM*SCALE |
IF( VALEIG ) THEN |
IF( VALEIG ) THEN |
* If eigenvalues in interval have to be found, |
* If eigenvalues in interval have to be found, |
* scale (WL, WU] accordingly |
* scale (WL, WU] accordingly |
WL = WL*SCALE |
WL = WL*SCALE |
WU = WU*SCALE |
WU = WU*SCALE |
ENDIF |
ENDIF |
END IF |
END IF |
* |
* |
* Compute the desired eigenvalues of the tridiagonal after splitting |
* Compute the desired eigenvalues of the tridiagonal after splitting |
* into smaller subblocks if the corresponding off-diagonal elements |
* into smaller subblocks if the corresponding off-diagonal elements |
* are small |
* are small |
* THRESH is the splitting parameter for DLARRE |
* THRESH is the splitting parameter for DLARRE |
* A negative THRESH forces the old splitting criterion based on the |
* A negative THRESH forces the old splitting criterion based on the |
* size of the off-diagonal. A positive THRESH switches to splitting |
* size of the off-diagonal. A positive THRESH switches to splitting |
* which preserves relative accuracy. |
* which preserves relative accuracy. |
* |
* |
IF( TRYRAC ) THEN |
IF( TRYRAC ) THEN |
* Test whether the matrix warrants the more expensive relative approach. |
* Test whether the matrix warrants the more expensive relative approach. |
CALL DLARRR( N, D, E, IINFO ) |
CALL DLARRR( N, D, E, IINFO ) |
ELSE |
ELSE |
* The user does not care about relative accurately eigenvalues |
* The user does not care about relative accurately eigenvalues |
IINFO = -1 |
IINFO = -1 |
ENDIF |
ENDIF |
* Set the splitting criterion |
* Set the splitting criterion |
IF (IINFO.EQ.0) THEN |
IF (IINFO.EQ.0) THEN |
THRESH = EPS |
THRESH = EPS |
ELSE |
ELSE |
THRESH = -EPS |
THRESH = -EPS |
* relative accuracy is desired but T does not guarantee it |
* relative accuracy is desired but T does not guarantee it |
TRYRAC = .FALSE. |
TRYRAC = .FALSE. |
ENDIF |
ENDIF |
* |
* |
IF( TRYRAC ) THEN |
IF( TRYRAC ) THEN |
* Copy original diagonal, needed to guarantee relative accuracy |
* Copy original diagonal, needed to guarantee relative accuracy |
CALL DCOPY(N,D,1,WORK(INDD),1) |
CALL DCOPY(N,D,1,WORK(INDD),1) |
ENDIF |
ENDIF |
* Store the squares of the offdiagonal values of T |
* Store the squares of the offdiagonal values of T |
DO 5 J = 1, N-1 |
DO 5 J = 1, N-1 |
WORK( INDE2+J-1 ) = E(J)**2 |
WORK( INDE2+J-1 ) = E(J)**2 |
5 CONTINUE |
5 CONTINUE |
|
|
* Set the tolerance parameters for bisection |
* Set the tolerance parameters for bisection |
IF( .NOT.WANTZ ) THEN |
IF( .NOT.WANTZ ) THEN |
* DLARRE computes the eigenvalues to full precision. |
* DLARRE computes the eigenvalues to full precision. |
RTOL1 = FOUR * EPS |
RTOL1 = FOUR * EPS |
RTOL2 = FOUR * EPS |
RTOL2 = FOUR * EPS |
ELSE |
ELSE |
* DLARRE computes the eigenvalues to less than full precision. |
* DLARRE computes the eigenvalues to less than full precision. |
* ZLARRV will refine the eigenvalue approximations, and we only |
* ZLARRV will refine the eigenvalue approximations, and we only |
* need less accurate initial bisection in DLARRE. |
* need less accurate initial bisection in DLARRE. |
* Note: these settings do only affect the subset case and DLARRE |
* Note: these settings do only affect the subset case and DLARRE |
RTOL1 = SQRT(EPS) |
RTOL1 = SQRT(EPS) |
RTOL2 = MAX( SQRT(EPS)*5.0D-3, FOUR * EPS ) |
RTOL2 = MAX( SQRT(EPS)*5.0D-3, FOUR * EPS ) |
ENDIF |
ENDIF |
CALL DLARRE( RANGE, N, WL, WU, IIL, IIU, D, E, |
CALL DLARRE( RANGE, N, WL, WU, IIL, IIU, D, E, |
$ WORK(INDE2), RTOL1, RTOL2, THRESH, NSPLIT, |
$ WORK(INDE2), RTOL1, RTOL2, THRESH, NSPLIT, |
$ IWORK( IINSPL ), M, W, WORK( INDERR ), |
$ IWORK( IINSPL ), M, W, WORK( INDERR ), |
$ WORK( INDGP ), IWORK( IINDBL ), |
$ WORK( INDGP ), IWORK( IINDBL ), |
$ IWORK( IINDW ), WORK( INDGRS ), PIVMIN, |
$ IWORK( IINDW ), WORK( INDGRS ), PIVMIN, |
$ WORK( INDWRK ), IWORK( IINDWK ), IINFO ) |
$ WORK( INDWRK ), IWORK( IINDWK ), IINFO ) |
IF( IINFO.NE.0 ) THEN |
IF( IINFO.NE.0 ) THEN |
INFO = 10 + ABS( IINFO ) |
INFO = 10 + ABS( IINFO ) |
RETURN |
RETURN |
END IF |
END IF |
* Note that if RANGE .NE. 'V', DLARRE computes bounds on the desired |
* Note that if RANGE .NE. 'V', DLARRE computes bounds on the desired |
* part of the spectrum. All desired eigenvalues are contained in |
* part of the spectrum. All desired eigenvalues are contained in |
* (WL,WU] |
* (WL,WU] |
|
|
|
|
IF( WANTZ ) THEN |
IF( WANTZ ) THEN |
* |
* |
* Compute the desired eigenvectors corresponding to the computed |
* Compute the desired eigenvectors corresponding to the computed |
* eigenvalues |
* eigenvalues |
* |
* |
CALL ZLARRV( N, WL, WU, D, E, |
CALL ZLARRV( N, WL, WU, D, E, |
$ PIVMIN, IWORK( IINSPL ), M, |
$ PIVMIN, IWORK( IINSPL ), M, |
$ 1, M, MINRGP, RTOL1, RTOL2, |
$ 1, M, MINRGP, RTOL1, RTOL2, |
$ W, WORK( INDERR ), WORK( INDGP ), IWORK( IINDBL ), |
$ W, WORK( INDERR ), WORK( INDGP ), IWORK( IINDBL ), |
$ IWORK( IINDW ), WORK( INDGRS ), Z, LDZ, |
$ IWORK( IINDW ), WORK( INDGRS ), Z, LDZ, |
$ ISUPPZ, WORK( INDWRK ), IWORK( IINDWK ), IINFO ) |
$ ISUPPZ, WORK( INDWRK ), IWORK( IINDWK ), IINFO ) |
IF( IINFO.NE.0 ) THEN |
IF( IINFO.NE.0 ) THEN |
INFO = 20 + ABS( IINFO ) |
INFO = 20 + ABS( IINFO ) |
RETURN |
RETURN |
END IF |
END IF |
ELSE |
ELSE |
* DLARRE computes eigenvalues of the (shifted) root representation |
* DLARRE computes eigenvalues of the (shifted) root representation |
* ZLARRV returns the eigenvalues of the unshifted matrix. |
* ZLARRV returns the eigenvalues of the unshifted matrix. |
* However, if the eigenvectors are not desired by the user, we need |
* However, if the eigenvectors are not desired by the user, we need |
* to apply the corresponding shifts from DLARRE to obtain the |
* to apply the corresponding shifts from DLARRE to obtain the |
* eigenvalues of the original matrix. |
* eigenvalues of the original matrix. |
DO 20 J = 1, M |
DO 20 J = 1, M |
ITMP = IWORK( IINDBL+J-1 ) |
ITMP = IWORK( IINDBL+J-1 ) |
W( J ) = W( J ) + E( IWORK( IINSPL+ITMP-1 ) ) |
W( J ) = W( J ) + E( IWORK( IINSPL+ITMP-1 ) ) |
20 CONTINUE |
20 CONTINUE |
END IF |
END IF |
* |
* |
|
|
IF ( TRYRAC ) THEN |
IF ( TRYRAC ) THEN |
* Refine computed eigenvalues so that they are relatively accurate |
* Refine computed eigenvalues so that they are relatively accurate |
* with respect to the original matrix T. |
* with respect to the original matrix T. |
IBEGIN = 1 |
IBEGIN = 1 |
WBEGIN = 1 |
WBEGIN = 1 |
DO 39 JBLK = 1, IWORK( IINDBL+M-1 ) |
DO 39 JBLK = 1, IWORK( IINDBL+M-1 ) |
IEND = IWORK( IINSPL+JBLK-1 ) |
IEND = IWORK( IINSPL+JBLK-1 ) |
IN = IEND - IBEGIN + 1 |
IN = IEND - IBEGIN + 1 |
WEND = WBEGIN - 1 |
WEND = WBEGIN - 1 |
* check if any eigenvalues have to be refined in this block |
* check if any eigenvalues have to be refined in this block |
36 CONTINUE |
36 CONTINUE |
IF( WEND.LT.M ) THEN |
IF( WEND.LT.M ) THEN |
IF( IWORK( IINDBL+WEND ).EQ.JBLK ) THEN |
IF( IWORK( IINDBL+WEND ).EQ.JBLK ) THEN |
WEND = WEND + 1 |
WEND = WEND + 1 |
GO TO 36 |
GO TO 36 |
|
END IF |
|
END IF |
|
IF( WEND.LT.WBEGIN ) THEN |
|
IBEGIN = IEND + 1 |
|
GO TO 39 |
END IF |
END IF |
END IF |
|
IF( WEND.LT.WBEGIN ) THEN |
|
IBEGIN = IEND + 1 |
|
GO TO 39 |
|
END IF |
|
|
|
OFFSET = IWORK(IINDW+WBEGIN-1)-1 |
OFFSET = IWORK(IINDW+WBEGIN-1)-1 |
IFIRST = IWORK(IINDW+WBEGIN-1) |
IFIRST = IWORK(IINDW+WBEGIN-1) |
ILAST = IWORK(IINDW+WEND-1) |
ILAST = IWORK(IINDW+WEND-1) |
RTOL2 = FOUR * EPS |
RTOL2 = FOUR * EPS |
CALL DLARRJ( IN, |
CALL DLARRJ( IN, |
$ WORK(INDD+IBEGIN-1), WORK(INDE2+IBEGIN-1), |
$ WORK(INDD+IBEGIN-1), WORK(INDE2+IBEGIN-1), |
$ IFIRST, ILAST, RTOL2, OFFSET, W(WBEGIN), |
$ IFIRST, ILAST, RTOL2, OFFSET, W(WBEGIN), |
$ WORK( INDERR+WBEGIN-1 ), |
$ WORK( INDERR+WBEGIN-1 ), |
$ WORK( INDWRK ), IWORK( IINDWK ), PIVMIN, |
$ WORK( INDWRK ), IWORK( IINDWK ), PIVMIN, |
$ TNRM, IINFO ) |
$ TNRM, IINFO ) |
IBEGIN = IEND + 1 |
IBEGIN = IEND + 1 |
WBEGIN = WEND + 1 |
WBEGIN = WEND + 1 |
39 CONTINUE |
39 CONTINUE |
ENDIF |
ENDIF |
* |
* |
* If matrix was scaled, then rescale eigenvalues appropriately. |
* If matrix was scaled, then rescale eigenvalues appropriately. |
* |
* |
IF( SCALE.NE.ONE ) THEN |
IF( SCALE.NE.ONE ) THEN |
CALL DSCAL( M, ONE / SCALE, W, 1 ) |
CALL DSCAL( M, ONE / SCALE, W, 1 ) |
|
END IF |
END IF |
END IF |
* |
* |
* If eigenvalues are not in increasing order, then sort them, |
* If eigenvalues are not in increasing order, then sort them, |
* possibly along with eigenvectors. |
* possibly along with eigenvectors. |
* |
* |
IF( NSPLIT.GT.1 ) THEN |
IF( NSPLIT.GT.1 .OR. N.EQ.2 ) THEN |
IF( .NOT. WANTZ ) THEN |
IF( .NOT. WANTZ ) THEN |
CALL DLASRT( 'I', M, W, IINFO ) |
CALL DLASRT( 'I', M, W, IINFO ) |
IF( IINFO.NE.0 ) THEN |
IF( IINFO.NE.0 ) THEN |