--- rpl/lapack/lapack/zstemr.f 2012/08/22 09:48:40 1.10 +++ rpl/lapack/lapack/zstemr.f 2012/12/14 12:30:34 1.11 @@ -311,7 +311,7 @@ *> \author Univ. of Colorado Denver *> \author NAG Ltd. * -*> \date November 2011 +*> \date September 2012 * *> \ingroup complex16OTHERcomputational * @@ -329,10 +329,10 @@ $ M, W, Z, LDZ, NZC, ISUPPZ, TRYRAC, WORK, LWORK, $ 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, -- * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- -* November 2011 +* September 2012 * * .. Scalar Arguments .. CHARACTER JOBZ, RANGE @@ -560,184 +560,184 @@ END IF ENDIF ENDIF - RETURN - END IF + ELSE -* Continue with general N +* Continue with general N - INDGRS = 1 - INDERR = 2*N + 1 - INDGP = 3*N + 1 - INDD = 4*N + 1 - INDE2 = 5*N + 1 - INDWRK = 6*N + 1 -* - IINSPL = 1 - IINDBL = N + 1 - IINDW = 2*N + 1 - IINDWK = 3*N + 1 -* -* Scale matrix to allowable range, if necessary. -* The allowable range is related to the PIVMIN parameter; see the -* comments in DLARRD. The preference for scaling small values -* up is heuristic; we expect users' matrices not to be close to the -* RMAX threshold. -* - SCALE = ONE - TNRM = DLANST( 'M', N, D, E ) - IF( TNRM.GT.ZERO .AND. TNRM.LT.RMIN ) THEN - SCALE = RMIN / TNRM - ELSE IF( TNRM.GT.RMAX ) THEN - SCALE = RMAX / TNRM - END IF - IF( SCALE.NE.ONE ) THEN - CALL DSCAL( N, SCALE, D, 1 ) - CALL DSCAL( N-1, SCALE, E, 1 ) - TNRM = TNRM*SCALE - IF( VALEIG ) THEN -* If eigenvalues in interval have to be found, -* scale (WL, WU] accordingly - WL = WL*SCALE - WU = WU*SCALE - ENDIF - END IF + INDGRS = 1 + INDERR = 2*N + 1 + INDGP = 3*N + 1 + INDD = 4*N + 1 + INDE2 = 5*N + 1 + INDWRK = 6*N + 1 +* + IINSPL = 1 + IINDBL = N + 1 + IINDW = 2*N + 1 + IINDWK = 3*N + 1 +* +* Scale matrix to allowable range, if necessary. +* The allowable range is related to the PIVMIN parameter; see the +* comments in DLARRD. The preference for scaling small values +* up is heuristic; we expect users' matrices not to be close to the +* RMAX threshold. +* + SCALE = ONE + TNRM = DLANST( 'M', N, D, E ) + IF( TNRM.GT.ZERO .AND. TNRM.LT.RMIN ) THEN + SCALE = RMIN / TNRM + ELSE IF( TNRM.GT.RMAX ) THEN + SCALE = RMAX / TNRM + END IF + IF( SCALE.NE.ONE ) THEN + CALL DSCAL( N, SCALE, D, 1 ) + CALL DSCAL( N-1, SCALE, E, 1 ) + TNRM = TNRM*SCALE + IF( VALEIG ) THEN +* If eigenvalues in interval have to be found, +* scale (WL, WU] accordingly + WL = WL*SCALE + WU = WU*SCALE + ENDIF + END IF * -* Compute the desired eigenvalues of the tridiagonal after splitting -* into smaller subblocks if the corresponding off-diagonal elements -* are small -* THRESH is the splitting parameter for DLARRE -* A negative THRESH forces the old splitting criterion based on the -* size of the off-diagonal. A positive THRESH switches to splitting -* which preserves relative accuracy. -* - IF( TRYRAC ) THEN -* Test whether the matrix warrants the more expensive relative approach. - CALL DLARRR( N, D, E, IINFO ) - ELSE -* The user does not care about relative accurately eigenvalues - IINFO = -1 - ENDIF -* Set the splitting criterion - IF (IINFO.EQ.0) THEN - THRESH = EPS - ELSE - THRESH = -EPS -* relative accuracy is desired but T does not guarantee it - TRYRAC = .FALSE. - ENDIF +* Compute the desired eigenvalues of the tridiagonal after splitting +* into smaller subblocks if the corresponding off-diagonal elements +* are small +* THRESH is the splitting parameter for DLARRE +* A negative THRESH forces the old splitting criterion based on the +* size of the off-diagonal. A positive THRESH switches to splitting +* which preserves relative accuracy. +* + IF( TRYRAC ) THEN +* Test whether the matrix warrants the more expensive relative approach. + CALL DLARRR( N, D, E, IINFO ) + ELSE +* The user does not care about relative accurately eigenvalues + IINFO = -1 + ENDIF +* Set the splitting criterion + IF (IINFO.EQ.0) THEN + THRESH = EPS + ELSE + THRESH = -EPS +* relative accuracy is desired but T does not guarantee it + TRYRAC = .FALSE. + ENDIF * - IF( TRYRAC ) THEN -* Copy original diagonal, needed to guarantee relative accuracy - CALL DCOPY(N,D,1,WORK(INDD),1) - ENDIF -* Store the squares of the offdiagonal values of T - DO 5 J = 1, N-1 - WORK( INDE2+J-1 ) = E(J)**2 + IF( TRYRAC ) THEN +* Copy original diagonal, needed to guarantee relative accuracy + CALL DCOPY(N,D,1,WORK(INDD),1) + ENDIF +* Store the squares of the offdiagonal values of T + DO 5 J = 1, N-1 + WORK( INDE2+J-1 ) = E(J)**2 5 CONTINUE -* Set the tolerance parameters for bisection - IF( .NOT.WANTZ ) THEN -* DLARRE computes the eigenvalues to full precision. - RTOL1 = FOUR * EPS - RTOL2 = FOUR * EPS - ELSE -* DLARRE computes the eigenvalues to less than full precision. -* ZLARRV will refine the eigenvalue approximations, and we only -* need less accurate initial bisection in DLARRE. -* Note: these settings do only affect the subset case and DLARRE - RTOL1 = SQRT(EPS) - RTOL2 = MAX( SQRT(EPS)*5.0D-3, FOUR * EPS ) - ENDIF - CALL DLARRE( RANGE, N, WL, WU, IIL, IIU, D, E, +* Set the tolerance parameters for bisection + IF( .NOT.WANTZ ) THEN +* DLARRE computes the eigenvalues to full precision. + RTOL1 = FOUR * EPS + RTOL2 = FOUR * EPS + ELSE +* DLARRE computes the eigenvalues to less than full precision. +* ZLARRV will refine the eigenvalue approximations, and we only +* need less accurate initial bisection in DLARRE. +* Note: these settings do only affect the subset case and DLARRE + RTOL1 = SQRT(EPS) + RTOL2 = MAX( SQRT(EPS)*5.0D-3, FOUR * EPS ) + ENDIF + CALL DLARRE( RANGE, N, WL, WU, IIL, IIU, D, E, $ WORK(INDE2), RTOL1, RTOL2, THRESH, NSPLIT, $ IWORK( IINSPL ), M, W, WORK( INDERR ), $ WORK( INDGP ), IWORK( IINDBL ), $ IWORK( IINDW ), WORK( INDGRS ), PIVMIN, $ WORK( INDWRK ), IWORK( IINDWK ), IINFO ) - IF( IINFO.NE.0 ) THEN - INFO = 10 + ABS( IINFO ) - RETURN - END IF -* Note that if RANGE .NE. 'V', DLARRE computes bounds on the desired -* part of the spectrum. All desired eigenvalues are contained in -* (WL,WU] + IF( IINFO.NE.0 ) THEN + INFO = 10 + ABS( IINFO ) + RETURN + END IF +* Note that if RANGE .NE. 'V', DLARRE computes bounds on the desired +* part of the spectrum. All desired eigenvalues are contained in +* (WL,WU] - IF( WANTZ ) THEN + IF( WANTZ ) THEN * -* Compute the desired eigenvectors corresponding to the computed -* eigenvalues +* Compute the desired eigenvectors corresponding to the computed +* eigenvalues * - CALL ZLARRV( N, WL, WU, D, E, + CALL ZLARRV( N, WL, WU, D, E, $ PIVMIN, IWORK( IINSPL ), M, $ 1, M, MINRGP, RTOL1, RTOL2, $ W, WORK( INDERR ), WORK( INDGP ), IWORK( IINDBL ), $ IWORK( IINDW ), WORK( INDGRS ), Z, LDZ, $ ISUPPZ, WORK( INDWRK ), IWORK( IINDWK ), IINFO ) - IF( IINFO.NE.0 ) THEN - INFO = 20 + ABS( IINFO ) - RETURN - END IF - ELSE -* DLARRE computes eigenvalues of the (shifted) root representation -* ZLARRV returns the eigenvalues of the unshifted matrix. -* However, if the eigenvectors are not desired by the user, we need -* to apply the corresponding shifts from DLARRE to obtain the -* eigenvalues of the original matrix. - DO 20 J = 1, M - ITMP = IWORK( IINDBL+J-1 ) - W( J ) = W( J ) + E( IWORK( IINSPL+ITMP-1 ) ) + IF( IINFO.NE.0 ) THEN + INFO = 20 + ABS( IINFO ) + RETURN + END IF + ELSE +* DLARRE computes eigenvalues of the (shifted) root representation +* ZLARRV returns the eigenvalues of the unshifted matrix. +* However, if the eigenvectors are not desired by the user, we need +* to apply the corresponding shifts from DLARRE to obtain the +* eigenvalues of the original matrix. + DO 20 J = 1, M + ITMP = IWORK( IINDBL+J-1 ) + W( J ) = W( J ) + E( IWORK( IINSPL+ITMP-1 ) ) 20 CONTINUE - END IF + END IF * - IF ( TRYRAC ) THEN -* Refine computed eigenvalues so that they are relatively accurate -* with respect to the original matrix T. - IBEGIN = 1 - WBEGIN = 1 - DO 39 JBLK = 1, IWORK( IINDBL+M-1 ) - IEND = IWORK( IINSPL+JBLK-1 ) - IN = IEND - IBEGIN + 1 - WEND = WBEGIN - 1 -* check if any eigenvalues have to be refined in this block + IF ( TRYRAC ) THEN +* Refine computed eigenvalues so that they are relatively accurate +* with respect to the original matrix T. + IBEGIN = 1 + WBEGIN = 1 + DO 39 JBLK = 1, IWORK( IINDBL+M-1 ) + IEND = IWORK( IINSPL+JBLK-1 ) + IN = IEND - IBEGIN + 1 + WEND = WBEGIN - 1 +* check if any eigenvalues have to be refined in this block 36 CONTINUE - IF( WEND.LT.M ) THEN - IF( IWORK( IINDBL+WEND ).EQ.JBLK ) THEN - WEND = WEND + 1 - GO TO 36 + IF( WEND.LT.M ) THEN + IF( IWORK( IINDBL+WEND ).EQ.JBLK ) THEN + WEND = WEND + 1 + GO TO 36 + END IF + END IF + IF( WEND.LT.WBEGIN ) THEN + IBEGIN = IEND + 1 + GO TO 39 END IF - END IF - IF( WEND.LT.WBEGIN ) THEN - IBEGIN = IEND + 1 - GO TO 39 - END IF - OFFSET = IWORK(IINDW+WBEGIN-1)-1 - IFIRST = IWORK(IINDW+WBEGIN-1) - ILAST = IWORK(IINDW+WEND-1) - RTOL2 = FOUR * EPS - CALL DLARRJ( IN, + OFFSET = IWORK(IINDW+WBEGIN-1)-1 + IFIRST = IWORK(IINDW+WBEGIN-1) + ILAST = IWORK(IINDW+WEND-1) + RTOL2 = FOUR * EPS + CALL DLARRJ( IN, $ WORK(INDD+IBEGIN-1), WORK(INDE2+IBEGIN-1), $ IFIRST, ILAST, RTOL2, OFFSET, W(WBEGIN), $ WORK( INDERR+WBEGIN-1 ), $ WORK( INDWRK ), IWORK( IINDWK ), PIVMIN, $ TNRM, IINFO ) - IBEGIN = IEND + 1 - WBEGIN = WEND + 1 + IBEGIN = IEND + 1 + WBEGIN = WEND + 1 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 - CALL DSCAL( M, ONE / SCALE, W, 1 ) + IF( SCALE.NE.ONE ) THEN + CALL DSCAL( M, ONE / SCALE, W, 1 ) + END IF END IF * * If eigenvalues are not in increasing order, then sort them, * possibly along with eigenvectors. * - IF( NSPLIT.GT.1 ) THEN + IF( NSPLIT.GT.1 .OR. N.EQ.2 ) THEN IF( .NOT. WANTZ ) THEN CALL DLASRT( 'I', M, W, IINFO ) IF( IINFO.NE.0 ) THEN