1 *> \brief <b> DSTEVR computes the eigenvalues and, optionally, the left and/or right eigenvectors for OTHER matrices</b>
3 * =========== DOCUMENTATION ===========
5 * Online html documentation available at
6 * http://www.netlib.org/lapack/explore-html/
9 *> Download DSTEVR + dependencies
10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dstevr.f">
12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dstevr.f">
14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dstevr.f">
21 * SUBROUTINE DSTEVR( JOBZ, RANGE, N, D, E, VL, VU, IL, IU, ABSTOL,
22 * M, W, Z, LDZ, ISUPPZ, WORK, LWORK, IWORK,
25 * .. Scalar Arguments ..
26 * CHARACTER JOBZ, RANGE
27 * INTEGER IL, INFO, IU, LDZ, LIWORK, LWORK, M, N
28 * DOUBLE PRECISION ABSTOL, VL, VU
30 * .. Array Arguments ..
31 * INTEGER ISUPPZ( * ), IWORK( * )
32 * DOUBLE PRECISION D( * ), E( * ), W( * ), WORK( * ), Z( LDZ, * )
41 *> DSTEVR computes selected eigenvalues and, optionally, eigenvectors
42 *> of a real symmetric tridiagonal matrix T. Eigenvalues and
43 *> eigenvectors can be selected by specifying either a range of values
44 *> or a range of indices for the desired eigenvalues.
46 *> Whenever possible, DSTEVR calls DSTEMR to compute the
47 *> eigenspectrum using Relatively Robust Representations. DSTEMR
48 *> computes eigenvalues by the dqds algorithm, while orthogonal
49 *> eigenvectors are computed from various "good" L D L^T representations
50 *> (also known as Relatively Robust Representations). Gram-Schmidt
51 *> orthogonalization is avoided as far as possible. More specifically,
52 *> the various steps of the algorithm are as follows. For the i-th
53 *> unreduced block of T,
54 *> (a) Compute T - sigma_i = L_i D_i L_i^T, such that L_i D_i L_i^T
55 *> is a relatively robust representation,
56 *> (b) Compute the eigenvalues, lambda_j, of L_i D_i L_i^T to high
57 *> relative accuracy by the dqds algorithm,
58 *> (c) If there is a cluster of close eigenvalues, "choose" sigma_i
59 *> close to the cluster, and go to step (a),
60 *> (d) Given the approximate eigenvalue lambda_j of L_i D_i L_i^T,
61 *> compute the corresponding eigenvector by forming a
62 *> rank-revealing twisted factorization.
63 *> The desired accuracy of the output can be specified by the input
66 *> For more details, see "A new O(n^2) algorithm for the symmetric
67 *> tridiagonal eigenvalue/eigenvector problem", by Inderjit Dhillon,
68 *> Computer Science Division Technical Report No. UCB//CSD-97-971,
69 *> UC Berkeley, May 1997.
72 *> Note 1 : DSTEVR calls DSTEMR when the full spectrum is requested
73 *> on machines which conform to the ieee-754 floating point standard.
74 *> DSTEVR calls DSTEBZ and DSTEIN on non-ieee machines and
75 *> when partial spectrum requests are made.
77 *> Normal execution of DSTEMR may create NaNs and infinities and
78 *> hence may abort due to a floating point exception in environments
79 *> which do not handle NaNs and infinities in the ieee standard default
88 *> JOBZ is CHARACTER*1
89 *> = 'N': Compute eigenvalues only;
90 *> = 'V': Compute eigenvalues and eigenvectors.
95 *> RANGE is CHARACTER*1
96 *> = 'A': all eigenvalues will be found.
97 *> = 'V': all eigenvalues in the half-open interval (VL,VU]
99 *> = 'I': the IL-th through IU-th eigenvalues will be found.
100 *> For RANGE = 'V' or 'I' and IU - IL < N - 1, DSTEBZ and
107 *> The order of the matrix. N >= 0.
112 *> D is DOUBLE PRECISION array, dimension (N)
113 *> On entry, the n diagonal elements of the tridiagonal matrix
115 *> On exit, D may be multiplied by a constant factor chosen
116 *> to avoid over/underflow in computing the eigenvalues.
121 *> E is DOUBLE PRECISION array, dimension (max(1,N-1))
122 *> On entry, the (n-1) subdiagonal elements of the tridiagonal
123 *> matrix A in elements 1 to N-1 of E.
124 *> On exit, E may be multiplied by a constant factor chosen
125 *> to avoid over/underflow in computing the eigenvalues.
130 *> VL is DOUBLE PRECISION
131 *> If RANGE='V', the lower bound of the interval to
132 *> be searched for eigenvalues. VL < VU.
133 *> Not referenced if RANGE = 'A' or 'I'.
138 *> VU is DOUBLE PRECISION
139 *> If RANGE='V', the upper bound of the interval to
140 *> be searched for eigenvalues. VL < VU.
141 *> Not referenced if RANGE = 'A' or 'I'.
147 *> If RANGE='I', the index of the
148 *> smallest eigenvalue to be returned.
149 *> 1 <= IL <= IU <= N, if N > 0; IL = 1 and IU = 0 if N = 0.
150 *> Not referenced if RANGE = 'A' or 'V'.
156 *> If RANGE='I', the index of the
157 *> largest eigenvalue to be returned.
158 *> 1 <= IL <= IU <= N, if N > 0; IL = 1 and IU = 0 if N = 0.
159 *> Not referenced if RANGE = 'A' or 'V'.
164 *> ABSTOL is DOUBLE PRECISION
165 *> The absolute error tolerance for the eigenvalues.
166 *> An approximate eigenvalue is accepted as converged
167 *> when it is determined to lie in an interval [a,b]
168 *> of width less than or equal to
170 *> ABSTOL + EPS * max( |a|,|b| ) ,
172 *> where EPS is the machine precision. If ABSTOL is less than
173 *> or equal to zero, then EPS*|T| will be used in its place,
174 *> where |T| is the 1-norm of the tridiagonal matrix obtained
175 *> by reducing A to tridiagonal form.
177 *> See "Computing Small Singular Values of Bidiagonal Matrices
178 *> with Guaranteed High Relative Accuracy," by Demmel and
179 *> Kahan, LAPACK Working Note #3.
181 *> If high relative accuracy is important, set ABSTOL to
182 *> DLAMCH( 'Safe minimum' ). Doing so will guarantee that
183 *> eigenvalues are computed to high relative accuracy when
184 *> possible in future releases. The current code does not
185 *> make any guarantees about high relative accuracy, but
186 *> future releases will. See J. Barlow and J. Demmel,
187 *> "Computing Accurate Eigensystems of Scaled Diagonally
188 *> Dominant Matrices", LAPACK Working Note #7, for a discussion
189 *> of which matrices define their eigenvalues to high relative
196 *> The total number of eigenvalues found. 0 <= M <= N.
197 *> If RANGE = 'A', M = N, and if RANGE = 'I', M = IU-IL+1.
202 *> W is DOUBLE PRECISION array, dimension (N)
203 *> The first M elements contain the selected eigenvalues in
209 *> Z is DOUBLE PRECISION array, dimension (LDZ, max(1,M) )
210 *> If JOBZ = 'V', then if INFO = 0, the first M columns of Z
211 *> contain the orthonormal eigenvectors of the matrix A
212 *> corresponding to the selected eigenvalues, with the i-th
213 *> column of Z holding the eigenvector associated with W(i).
214 *> Note: the user must ensure that at least max(1,M) columns are
215 *> supplied in the array Z; if RANGE = 'V', the exact value of M
216 *> is not known in advance and an upper bound must be used.
222 *> The leading dimension of the array Z. LDZ >= 1, and if
223 *> JOBZ = 'V', LDZ >= max(1,N).
226 *> \param[out] ISUPPZ
228 *> ISUPPZ is INTEGER array, dimension ( 2*max(1,M) )
229 *> The support of the eigenvectors in Z, i.e., the indices
230 *> indicating the nonzero elements in Z. The i-th eigenvector
231 *> is nonzero only in elements ISUPPZ( 2*i-1 ) through
233 *> Implemented only for RANGE = 'A' or 'I' and IU - IL = N - 1
238 *> WORK is DOUBLE PRECISION array, dimension (MAX(1,LWORK))
239 *> On exit, if INFO = 0, WORK(1) returns the optimal (and
246 *> The dimension of the array WORK. LWORK >= max(1,20*N).
248 *> If LWORK = -1, then a workspace query is assumed; the routine
249 *> only calculates the optimal sizes of the WORK and IWORK
250 *> arrays, returns these values as the first entries of the WORK
251 *> and IWORK arrays, and no error message related to LWORK or
252 *> LIWORK is issued by XERBLA.
257 *> IWORK is INTEGER array, dimension (MAX(1,LIWORK))
258 *> On exit, if INFO = 0, IWORK(1) returns the optimal (and
265 *> The dimension of the array IWORK. LIWORK >= max(1,10*N).
267 *> If LIWORK = -1, then a workspace query is assumed; the
268 *> routine only calculates the optimal sizes of the WORK and
269 *> IWORK arrays, returns these values as the first entries of
270 *> the WORK and IWORK arrays, and no error message related to
271 *> LWORK or LIWORK is issued by XERBLA.
277 *> = 0: successful exit
278 *> < 0: if INFO = -i, the i-th argument had an illegal value
279 *> > 0: Internal error
285 *> \author Univ. of Tennessee
286 *> \author Univ. of California Berkeley
287 *> \author Univ. of Colorado Denver
292 *> \ingroup doubleOTHEReigen
294 *> \par Contributors:
297 *> Inderjit Dhillon, IBM Almaden, USA \n
298 *> Osni Marques, LBNL/NERSC, USA \n
299 *> Ken Stanley, Computer Science Division, University of
300 *> California at Berkeley, USA \n
302 * =====================================================================
303 SUBROUTINE DSTEVR( JOBZ, RANGE, N, D, E, VL, VU, IL, IU, ABSTOL,
304 $ M, W, Z, LDZ, ISUPPZ, WORK, LWORK, IWORK,
307 * -- LAPACK driver routine (version 3.6.1) --
308 * -- LAPACK is a software package provided by Univ. of Tennessee, --
309 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
312 * .. Scalar Arguments ..
313 CHARACTER JOBZ, RANGE
314 INTEGER IL, INFO, IU, LDZ, LIWORK, LWORK, M, N
315 DOUBLE PRECISION ABSTOL, VL, VU
317 * .. Array Arguments ..
318 INTEGER ISUPPZ( * ), IWORK( * )
319 DOUBLE PRECISION D( * ), E( * ), W( * ), WORK( * ), Z( LDZ, * )
322 * =====================================================================
325 DOUBLE PRECISION ZERO, ONE, TWO
326 PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0, TWO = 2.0D+0 )
328 * .. Local Scalars ..
329 LOGICAL ALLEIG, INDEIG, TEST, LQUERY, VALEIG, WANTZ,
332 INTEGER I, IEEEOK, IMAX, INDIBL, INDIFL, INDISP,
333 $ INDIWO, ISCALE, ITMP1, J, JJ, LIWMIN, LWMIN,
335 DOUBLE PRECISION BIGNUM, EPS, RMAX, RMIN, SAFMIN, SIGMA, SMLNUM,
336 $ TMP1, TNRM, VLL, VUU
338 * .. External Functions ..
341 DOUBLE PRECISION DLAMCH, DLANST
342 EXTERNAL LSAME, ILAENV, DLAMCH, DLANST
344 * .. External Subroutines ..
345 EXTERNAL DCOPY, DSCAL, DSTEBZ, DSTEMR, DSTEIN, DSTERF,
348 * .. Intrinsic Functions ..
349 INTRINSIC MAX, MIN, SQRT
351 * .. Executable Statements ..
354 * Test the input parameters.
356 IEEEOK = ILAENV( 10, 'DSTEVR', 'N', 1, 2, 3, 4 )
358 WANTZ = LSAME( JOBZ, 'V' )
359 ALLEIG = LSAME( RANGE, 'A' )
360 VALEIG = LSAME( RANGE, 'V' )
361 INDEIG = LSAME( RANGE, 'I' )
363 LQUERY = ( ( LWORK.EQ.-1 ) .OR. ( LIWORK.EQ.-1 ) )
364 LWMIN = MAX( 1, 20*N )
365 LIWMIN = MAX( 1, 10*N )
369 IF( .NOT.( WANTZ .OR. LSAME( JOBZ, 'N' ) ) ) THEN
371 ELSE IF( .NOT.( ALLEIG .OR. VALEIG .OR. INDEIG ) ) THEN
373 ELSE IF( N.LT.0 ) THEN
377 IF( N.GT.0 .AND. VU.LE.VL )
379 ELSE IF( INDEIG ) THEN
380 IF( IL.LT.1 .OR. IL.GT.MAX( 1, N ) ) THEN
382 ELSE IF( IU.LT.MIN( N, IL ) .OR. IU.GT.N ) THEN
388 IF( LDZ.LT.1 .OR. ( WANTZ .AND. LDZ.LT.N ) ) THEN
397 IF( LWORK.LT.LWMIN .AND. .NOT.LQUERY ) THEN
399 ELSE IF( LIWORK.LT.LIWMIN .AND. .NOT.LQUERY ) THEN
405 CALL XERBLA( 'DSTEVR', -INFO )
407 ELSE IF( LQUERY ) THEN
411 * Quick return if possible
418 IF( ALLEIG .OR. INDEIG ) THEN
422 IF( VL.LT.D( 1 ) .AND. VU.GE.D( 1 ) ) THEN
432 * Get machine constants.
434 SAFMIN = DLAMCH( 'Safe minimum' )
435 EPS = DLAMCH( 'Precision' )
436 SMLNUM = SAFMIN / EPS
437 BIGNUM = ONE / SMLNUM
438 RMIN = SQRT( SMLNUM )
439 RMAX = MIN( SQRT( BIGNUM ), ONE / SQRT( SQRT( SAFMIN ) ) )
442 * Scale matrix to allowable range, if necessary.
450 TNRM = DLANST( 'M', N, D, E )
451 IF( TNRM.GT.ZERO .AND. TNRM.LT.RMIN ) THEN
454 ELSE IF( TNRM.GT.RMAX ) THEN
458 IF( ISCALE.EQ.1 ) THEN
459 CALL DSCAL( N, SIGMA, D, 1 )
460 CALL DSCAL( N-1, SIGMA, E( 1 ), 1 )
467 * Initialize indices into workspaces. Note: These indices are used only
468 * if DSTERF or DSTEMR fail.
470 * IWORK(INDIBL:INDIBL+M-1) corresponds to IBLOCK in DSTEBZ and
471 * stores the block indices of each of the M<=N eigenvalues.
473 * IWORK(INDISP:INDISP+NSPLIT-1) corresponds to ISPLIT in DSTEBZ and
474 * stores the starting and finishing indices of each block.
476 * IWORK(INDIFL:INDIFL+N-1) stores the indices of eigenvectors
477 * that corresponding to eigenvectors that fail to converge in
478 * DSTEIN. This information is discarded; if any fail, the driver
481 * INDIWO is the offset of the remaining integer workspace.
484 * If all eigenvalues are desired, then
485 * call DSTERF or DSTEMR. If this fails for some eigenvalue, then
491 IF( IL.EQ.1 .AND. IU.EQ.N ) THEN
495 IF( ( ALLEIG .OR. TEST ) .AND. IEEEOK.EQ.1 ) THEN
496 CALL DCOPY( N-1, E( 1 ), 1, WORK( 1 ), 1 )
497 IF( .NOT.WANTZ ) THEN
498 CALL DCOPY( N, D, 1, W, 1 )
499 CALL DSTERF( N, W, WORK, INFO )
501 CALL DCOPY( N, D, 1, WORK( N+1 ), 1 )
502 IF (ABSTOL .LE. TWO*N*EPS) THEN
507 CALL DSTEMR( JOBZ, 'A', N, WORK( N+1 ), WORK, VL, VU, IL,
508 $ IU, M, W, Z, LDZ, N, ISUPPZ, TRYRAC,
509 $ WORK( 2*N+1 ), LWORK-2*N, IWORK, LIWORK, INFO )
519 * Otherwise, call DSTEBZ and, if eigenvectors are desired, DSTEIN.
527 CALL DSTEBZ( RANGE, ORDER, N, VLL, VUU, IL, IU, ABSTOL, D, E, M,
528 $ NSPLIT, W, IWORK( INDIBL ), IWORK( INDISP ), WORK,
529 $ IWORK( INDIWO ), INFO )
532 CALL DSTEIN( N, D, E, M, W, IWORK( INDIBL ), IWORK( INDISP ),
533 $ Z, LDZ, WORK, IWORK( INDIWO ), IWORK( INDIFL ),
537 * If matrix was scaled, then rescale eigenvalues appropriately.
540 IF( ISCALE.EQ.1 ) THEN
546 CALL DSCAL( IMAX, ONE / SIGMA, W, 1 )
549 * If eigenvalues are not in order, then sort them, along with
557 IF( W( JJ ).LT.TMP1 ) THEN
566 IWORK( I ) = IWORK( J )
569 CALL DSWAP( N, Z( 1, I ), 1, Z( 1, J ), 1 )
574 * Causes problems with tests 19 & 20:
575 * IF (wantz .and. INDEIG ) Z( 1,1) = Z(1,1) / 1.002 + .002