--- rpl/lapack/lapack/zstemr.f 2011/11/21 20:43:21 1.8
+++ rpl/lapack/lapack/zstemr.f 2017/06/17 10:54:28 1.18
@@ -2,18 +2,18 @@
*
* =========== DOCUMENTATION ===========
*
-* Online html documentation available at
-* http://www.netlib.org/lapack/explore-html/
+* Online html documentation available at
+* http://www.netlib.org/lapack/explore-html/
*
*> \htmlonly
-*> Download ZSTEMR + dependencies
-*>
-*> [TGZ]
-*>
-*> [ZIP]
-*>
+*> Download ZSTEMR + dependencies
+*>
+*> [TGZ]
+*>
+*> [ZIP]
+*>
*> [TXT]
-*> \endhtmlonly
+*> \endhtmlonly
*
* Definition:
* ===========
@@ -21,7 +21,7 @@
* SUBROUTINE ZSTEMR( JOBZ, RANGE, N, D, E, VL, VU, IL, IU,
* M, W, Z, LDZ, NZC, ISUPPZ, TRYRAC, WORK, LWORK,
* IWORK, LIWORK, INFO )
-*
+*
* .. Scalar Arguments ..
* CHARACTER JOBZ, RANGE
* LOGICAL TRYRAC
@@ -33,7 +33,7 @@
* DOUBLE PRECISION D( * ), E( * ), W( * ), WORK( * )
* COMPLEX*16 Z( LDZ, * )
* ..
-*
+*
*
*> \par Purpose:
* =============
@@ -153,13 +153,17 @@
*> \param[in] VL
*> \verbatim
*> 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
*>
*> \param[in] VU
*> \verbatim
*> 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.
*> Not referenced if RANGE = 'A' or 'I'.
*> \endverbatim
@@ -167,14 +171,19 @@
*> \param[in] IL
*> \verbatim
*> 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
*>
*> \param[in] IU
*> \verbatim
*> IU is INTEGER
*>
-*> If RANGE='I', the indices (in ascending order) of the
-*> smallest and largest eigenvalues to be returned.
+*> If RANGE='I', the index of the
+*> largest eigenvalue to be returned.
*> 1 <= IL <= IU <= N, if N > 0.
*> Not referenced if RANGE = 'A' or 'V'.
*> \endverbatim
@@ -306,12 +315,12 @@
* Authors:
* ========
*
-*> \author Univ. of Tennessee
-*> \author Univ. of California Berkeley
-*> \author Univ. of Colorado Denver
-*> \author NAG Ltd.
+*> \author Univ. of Tennessee
+*> \author Univ. of California Berkeley
+*> \author Univ. of Colorado Denver
+*> \author NAG Ltd.
*
-*> \date November 2011
+*> \date June 2016
*
*> \ingroup complex16OTHERcomputational
*
@@ -329,10 +338,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.7.0) --
* -- LAPACK is a software package provided by Univ. of Tennessee, --
* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
-* November 2011
+* June 2016
*
* .. Scalar Arguments ..
CHARACTER JOBZ, RANGE
@@ -408,6 +417,7 @@
WU = ZERO
IIL = 0
IIU = 0
+ NSPLIT = 0
IF( VALEIG ) THEN
* We do not reference VL, VU in the cases RANGE = 'I','A'
@@ -525,10 +535,10 @@
IF (SN.NE.ZERO) THEN
IF (CS.NE.ZERO) THEN
ISUPPZ(2*M-1) = 1
- ISUPPZ(2*M-1) = 2
+ ISUPPZ(2*M) = 2
ELSE
ISUPPZ(2*M-1) = 1
- ISUPPZ(2*M-1) = 1
+ ISUPPZ(2*M) = 1
END IF
ELSE
ISUPPZ(2*M-1) = 2
@@ -549,10 +559,10 @@
IF (SN.NE.ZERO) THEN
IF (CS.NE.ZERO) THEN
ISUPPZ(2*M-1) = 1
- ISUPPZ(2*M-1) = 2
+ ISUPPZ(2*M) = 2
ELSE
ISUPPZ(2*M-1) = 1
- ISUPPZ(2*M-1) = 1
+ ISUPPZ(2*M) = 1
END IF
ELSE
ISUPPZ(2*M-1) = 2
@@ -560,184 +570,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