version 1.9, 2011/11/21 22:19:57
|
version 1.12, 2012/12/14 14:22:54
|
Line 311
|
Line 311
|
*> \author Univ. of Colorado Denver |
*> \author Univ. of Colorado Denver |
*> \author NAG Ltd. |
*> \author NAG Ltd. |
* |
* |
*> \date November 2011 |
*> \date September 2012 |
* |
* |
*> \ingroup complex16OTHERcomputational |
*> \ingroup complex16OTHERcomputational |
* |
* |
Line 329
|
Line 329
|
$ 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.4.2) -- |
* -- 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 |
* September 2012 |
* |
* |
* .. Scalar Arguments .. |
* .. Scalar Arguments .. |
CHARACTER JOBZ, RANGE |
CHARACTER JOBZ, RANGE |
Line 560
|
Line 560
|
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 |